# A 3D-shell model of left atrial electromechanics

Oscar Ruz, Carlos Brito-Pacheco, Marina Vidrascu, Dominique Chapelle, Miguel Angel Fernández

► To cite this version:

Oscar Ruz, Carlos Brito-Pacheco, Marina Vidrascu, Dominique Chapelle, Miguel Angel Fernández. A 3D-shell model of left atrial electromechanics. 2026. hal-05483567v2

HAL Id: hal-05483567

https://inria.hal.science/hal-05483567v2

Preprint submitted on 30 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A 3D-Shell Model of Left Atrial Electromechanics

Oscar Ruz<sup>1,2</sup>, Carlos Brito-Pacheco<sup>2</sup>, Marina Vidrascu<sup>2</sup>,
Dominique Chapelle<sup>3</sup>, Miguel A. Fernández<sup>2*</sup>

<sup>1</sup>Universidad de Santiago de Chile, Santiago, Chile.

<sup>2</sup>Sorbonne Université, CNRS, LJLL, Inria, Paris, France.

<sup>3</sup>Inria & CMAP, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France.

*Corresponding author(s). E-mail(s): miguel.fernandez@inria.fr;

Contributing authors: oscar.ruz@usach.cl; carlos.brito-pacheco@inria.fr;

marina.vidrascu@inria.fr; dominique.chapelle@inria.fr;

## Abstract

The thin-walled nature of the atrial myocardium can lead to artificial stiffening when full 3D electromechanical models are discretized using standard finite elements. In this work, we propose an electromechanical model of the left atrium based on a 3D-shell formulation that overcomes these limitations. The model incorporates both passive and active components of atrial tissue mechanics, while atrioventricular interaction is described by the coupling with a 0D electromechanical model of the left ventricle. The proposed approach is assessed under physiological and pathological conditions and systematically compared with the standard full 3D formulation. The results demonstrate the superior robustness and computational efficiency of the proposed 3D-shell electromechanical model.

**Keywords:** Cardiac mechanics, Left atrium, 3D-shell model, Lumped-parameter model, Finite element method.

## 1 Introduction

Computational modeling of cardiac electromechanics has become an essential tool in the analysis of heart function. Such simulations provide quantitative insights into internal mechanical variables that cannot be accessed non-invasively and enable in-silico evaluation of therapeutic strategies and device performance (see, e.g., [1–5]).
The underlying mathematical models must capture the detailed multiscale and multiphysics processes that govern cardiac function, such as electrical activation, active and passive tissue mechanics and coupling with blood flow. The atria contribute to ventricular filling in three distinct ways: by acting as a reservoir during ventricular systole, by allowing passive blood flow into the ventricles during early diastole, and by augmenting ventricular filling during late diastole thanks to atrial contraction. Atrial anatomy is highly complex, with markedly heterogeneous wall thickness that can be as low as 0.5 mm (see, e.g., [6, 7]). Anatomically detailed geometric models with spatially varying wall thickness can be obtained from CT imaging data (see, e.g., [6]). In contrast to the ventricles, atrial myocardial fiber architecture is considerably more complex and does not appear to follow a simple pattern. Moreover, the thin-walled nature of the atrial myocardium complicates the acquisition of patient-specific fiber data (see, e.g., [8]). Rule-based methods are a widely used approach for prescribing realistic fiber orientations in computational atrial geometries (see, e.g., [9–16]).

While ventricular electromechanical modeling has been extensively studied over the past two decades (see, e.g., [17–28]), the literature on atrial mechanics remains comparatively scarce. The first mechanical models of the atria treated the myocardial tissue as a purely passive material, neglecting its contractile function (see, e.g., [29–32]). Only recently has atrial contraction begun to receive significant attention (see [8, 33–46]). Several of these studies have demonstrated good agreement with key physiological indicators, such as the characteristic figure-eight shape of the pressure-volume loop. However, the vast majority of these works disregard the numerical challenges associated with the thin-walled nature of the atrial myocardium. Indeed, with the exception of the studies reported in [29, 33], which are based on a Reissner-Mindlin shell formulation, all the aforementioned approaches model atrial tissue using a standard displacement-based 3D formulation, as commonly done for ventricular myocardium.

Standard finite element approximations of displacement-based 3D solid mechanics models are known to suffer from severe robustness issues, namely artificial stiffening, when the thickness is small compared to the other characteristic dimensions. This so-called numerical locking phenomenon often necessitates highly refined meshes, thereby leading to prohibitive high computational costs. A well-established strategy in solid mechanics to circumvent these issues is to adopt reduced-order modeling approaches based on shell theory, combined with appropriate techniques to mitigate numerical locking (see, e.g., [47, 48] and the references therein). In particular, 3D-shell elements with quadratic shape functions through the thickness are particularly appealing, as they allow for an accurate representation of transverse strains, enable the use of general constitutive laws and facilitate coupling with surrounding media (see [48]).

In this paper, we introduce an electromechanical model of the left atrium based on the 3D-shell formulation introduced in [48]. Additionally, a Mixed Interpolation of Tensorial Components (MITC) reinterpolation technique is considered to mitigate potential locking effects, which are common in thin-walled solid mechanics models (see, e.g., [49]). In order to account for atrioventricular interaction, the atrial model is coupled with a reduced 0D model of the left ventricle electromechanics [50]. The capabilities of the proposed approach are investigated through a series of numerical
experiments under physiological and pathological conditions, and comparisons are also made with a standard displacement-based 3D solid model.

A preliminary version of this work was presented in the FIMH conference paper [51]. The present study significantly extends that contribution in several aspects. First, atrial mechanics are coupled with a 0D reduced-order electromechanical model of the left ventricle, resulting in a more physiologically consistent description of the inter-chamber pressure-volume interaction. Second, mitral valve motion is incorporated through viscoelastic Robin-type boundary conditions. Finally, the proposed 3D-shell formulation is systematically assessed through qualitative and quantitative comparisons with standard full 3D finite element formulations. This includes a comprehensive set of physiological biomarkers, such as atrial stroke work, ejection fraction and pressure-volume loops.

# 2 Materials and Methods

In this section, we provide the detailed formulation of the left atrial electromechanical model. We introduce the geometric setting, kinematic variables, and the Lagrangian framework, followed by the passive and active constitutive laws and the equations governing the contractile dynamics. Finally, the numerical discretization of the resulting coupled nonlinear problem is presented.

## 2.1 Electromechanical model of the left atrium

We consider a 3D anatomical geometry of the left atrium represented by a domain $\Omega \subset \mathbb{R}^3$ in its reference configuration (see Figure 1). Its boundary $\partial\Omega$ is assumed to be partitioned into four disjoint parts, $\partial\Omega = \Gamma_{vs} \cup \Gamma_{mv} \cup \Gamma_{en} \cup \Gamma_{ep}$, where $\Gamma_{vs}$ corresponds to the pulmonary-vein inlets, $\Gamma_{mv}$ to the mitral-valve annulus interfacing with the left ventricle, and $\Gamma_{en}$ and $\Gamma_{ep}$ denote the endocardial and epicardial surfaces, respectively. The exterior unit normal vector to $\partial\Omega$ is denoted by $\mathbf{n}$. The mechanical response of the atrial cardiac tissue will be described using a total Lagrangian framework. To this purpose, we let $\varphi: \Omega \times \mathbb{R}^+ \to \mathbb{R}^3$ denote the motion map, assigning at each material point $\mathbf{x} \in \Omega$ its position $\varphi(\mathbf{x}, t)$ in the current configuration at time $t \ge 0$. The associated displacement field is defined by $\mathbf{y}(\mathbf{x}, t) \stackrel{\text{def}}{=} \varphi(\mathbf{x}, t) - \mathbf{x}$. The associated deformation gradient, characterizing the local stretch and rotation of material line elements, is given by $\mathbf{F} \stackrel{\text{def}}{=} \nabla\varphi$, while the Jacobian $J \stackrel{\text{def}}{=} \det \mathbf{F}$ represents the local volume ratio between the deformed and reference configurations. We shall also make use of the right Cauchy-Green strain tensor given by $\mathbf{C} \stackrel{\text{def}}{=} \mathbf{F}^\top \mathbf{F}$.

To model the constitutive material behavior, we adopt the multiscale electromechanical framework introduced in [23], which couples the passive and active mechanisms of the cardiac tissue, capturing their interaction through the nonlinear Hill-Maxwell extension. This results in the following additive decomposition of the second Piola-Kirchhoff stress tensor $\Sigma$:

$$ \Sigma \stackrel{\text{def}}{=} \Sigma_p + \Sigma_a \quad \text{in } \Omega, $$
<figure><img src="image_2.png" /><figcaption>Fig. 1: Geometrical description of the left atrium.</figcaption></figure>

where $\Sigma_p$ accounts for the passive stress behavior of the myocardium and $\Sigma_a$ denotes the active stress.

The passive response is modeled as a combination of hyperelastic and viscous contributions, formulated as

$$ \Sigma_p \stackrel{\text{def}}{=} \frac{\partial W_e}{\partial \dot{E}} + \frac{\partial W_v}{\partial \dot{E}}, $$

where $E \stackrel{\text{def}}{=} \frac{1}{2}(C - I)$ stands for the Green–Lagrange strain tensor (with $I$ denoting the identity matrix in $\mathbb{R}^{3 \times 3}$), $\dot{E} \stackrel{\text{def}}{=} \partial_t E$ and $W_e$ and $W_v$ are the hyperelastic potential and the viscous pseudo-potential, respectively. The myocardium is assumed to be a nearly incompressible, transversely isotropic material with a preferred direction aligned with the local fiber orientation, represented by unit tangent field $f$. Following previous studies on atrial mechanics [34, 37, 52], we consider the Demiray-Gasser strain-energy density function [53, 54]

$$ W_e \stackrel{\text{def}}{=} \frac{a_p}{2b_p} (e^{b_p(J_1-3)} - 1) + \frac{a_f}{2b_f} (e^{b_f(J_4-1)^2} - 1) + \kappa \log^2 J, \quad (1) $$

where $a_p, b_p, a_f, b_f, \kappa > 0$ are given constitutive parameters. The strain energy depends on the reduced invariants $J_1 \stackrel{\text{def}}{=} J^{-2/3} \text{tr}C$ and $J_4 \stackrel{\text{def}}{=} J^{-2/3} f^\top C f$. The viscous behavior is described through the pseudo-potential

$$ W_v \stackrel{\text{def}}{=} \frac{\nu}{2} \dot{E} : \dot{E}, $$

where $\nu > 0$ denotes the passive viscosity coefficient.

The active behavior is modeled by the multiscale approach outlined in [23, 55]. The contractile response arises from chemical processes in myocytes that induce sarcomere contraction and relaxation, producing the subsequent deformation of the myocardial
tissue along the fibers. This microscopic deformation induces a change in the fiber extension $e_c$, which drives both the active tension $\tau_c$ and the active stiffness $k_c$ of the sarcomeres. The dynamics of these contractile quantities are described by the following extension of the model originally proposed in [56]:

$$ \begin{cases} \dot{k}_c = -(|u_e| + \alpha_d|\dot{e}_c|)k_c + n_0(e_c)k_0|u_e|_+ & \text{in } \Omega, \\ \dot{\tau}_c = -(|u_e| + \alpha_d|\dot{e}_c|)\tau_c + k_c \dot{e}_c + n_0(e_c)\sigma_0 |u_e|_+ & \text{in } \Omega, \end{cases} $$

where $k_0$ and $\sigma_0$ respectively denote the maximum stiffness and contractility. The term $\alpha_d|\dot{e}_c|$ accounts for the bridges destruction upon rapid length changes. The parameter $0 \le n_0(e_c) \le 1$ is a reduction factor that accounts for the Frank-Starling effect, which relates the force-generating capacity of the cardiac tissue based on its degree of stretch. Finally, $|x|_+ \stackrel{\text{def}}{=} \max\{x, 0\}$ and $u_e$ represents the chemical activation. This source term is modeled phenomenologically as an affine function of the tissue transmembrane potential $V_m$:

$$ u_e = \alpha V_m + \beta \quad \text{in } \Omega, \qquad (2) $$

where $\alpha$ and $\beta$ are given constants (see, e.g., [20]). The transmembrane potential $V_m$ is obtained from the bidomain model of cardiac electrophysiology (see Appendix A). The contractile and elastic stresses within the sarcomeres are respectively given by

$$ \sigma_c = \tau_c + \mu_c \dot{e}_c, \qquad (3) $$

and

$$ \sigma_s = E_s e_s, \qquad (4) $$

where $\mu_c > 0$ is the viscous dissipation parameter and $e_s$ for the elastic extension of the sarcomeres, while $E_s$ denotes the elastic modulus of the passive sarcomere response. The total stress along the fibers, denoted by $\sigma_a$, can be expressed in terms of contractile and elastic stresses via the relation (see [23])

$$ \sigma_a = \frac{\sigma_c}{1 + 2e_s} = \frac{\sigma_s}{1 + 2e_c}. \qquad (5) $$

Furthermore, the total fiber extension $e_f \stackrel{\text{def}}{=} f^T E f$ satisfies

$$ 1 + 2e_f = (1 + 2e_s)(1 + 2e_c). \qquad (6) $$

From (5) and (6), the macroscopic active stress can be written as

$$ \sigma_a = E_s \frac{e_f - e_c}{(1 + 2e_c)^2}. \qquad (7) $$

By combining (3) and (5) with (7), the evolution of the contractile extension $e_c$ is given by the relation

$$ \mu_c \dot{e}_c + \tau_c = E_s \frac{(e_f - e_c)(1 + 2e_f)}{(1 + 2e_c)^3}. \qquad (8) $$
Finally, the active part of the second Piola–Kirchhoff stress tensor is expressed as

$$ \Sigma_a = \sigma_a \boldsymbol{f} \otimes \boldsymbol{f} = \frac{E_s (\boldsymbol{f}^\mathrm{T} \boldsymbol{E} \boldsymbol{f} - e_c)}{(1 + 2e_c)^2} \boldsymbol{f} \otimes \boldsymbol{f}. \qquad (9) $$

Considering the standard 3D space of admissible solid displacements

$$ V \stackrel{\mathrm{def}}{=} \{ \boldsymbol{v} \in [H^1(\Omega)]^3 : \boldsymbol{v}|_{\Gamma_{\mathrm{vs}}} = \mathbf{0} \}, \qquad (10) $$

the weak form of the resulting full 3D electromechanical model of the left atrium reads: For $t > 0$, find the displacements $\mathbf{y}(t) \in V$, with $\dot{\mathbf{y}} = \partial_t \mathbf{y}$, the extension along the fibers $e_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}$, the active stiffness $k_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}^+$ and the active tension $\tau_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{l} \rho \int_{\Omega} \partial_t \dot{\mathbf{y}} \cdot \boldsymbol{v} + \int_{\Omega} \boldsymbol{\Sigma} : D_y \boldsymbol{E} \cdot \boldsymbol{v} + \int_{\Gamma_{\mathrm{mv}}} (a\mathbf{y} + b\dot{\mathbf{y}}) \cdot \boldsymbol{v} = - \int_{\Gamma_{\mathrm{en}}} p_{\mathrm{la}} J \boldsymbol{F}^{-\mathrm{T}} \boldsymbol{n} \cdot \boldsymbol{v}, \\ \mu_c \dot{e}_c + \tau_c = E_s \frac{(\boldsymbol{f}^\mathrm{T} \boldsymbol{E} \boldsymbol{f} - e_c)(1 + 2\boldsymbol{f}^\mathrm{T} \boldsymbol{E} \boldsymbol{f})}{(1 + 2e_c)^3} \quad \text{in } \Omega, \\ \dot{k}_c = -(|u_e| + \alpha_d |\dot{e}_c|) k_c + n_0(e_c) k_0 |u_e|_+ \quad \text{in } \Omega, \\ \dot{\tau}_c = -(|u_e| + \alpha_d |\dot{e}_c|) \tau_c + \dot{e}_c k_c + n_0(e_c) \sigma_0 |u_e|_+ \quad \text{in } \Omega \end{array} \right. \qquad (11) $$

for all $\boldsymbol{v} \in V$ and where $\rho > 0$ denotes the volumic mass of the atrial wall. Note that the annuli of the pulmonary veins inlets $\Gamma_{\mathrm{vs}}$ are assumed to be fixed and a viscoelastic boundary condition is enforced on the the mitral valve annulus $\Gamma_{\mathrm{mv}}$. The dynamics of the unknown internal pressure $p_{\mathrm{la}}$ in the atrial cavity will be described in Section 2.1.2 below using a lumped parameter modeling.

### 2.1.1 3D-shell modeling of the left atrium

A fundamental drawback of (11)$_1$ is that due to the intrinsic thin-walled nature of $\Omega$, standard 3D finite element discretizations are prone to become highly sensitive to wall thickness, leading to well-known robustness issues. In particular, the resulting solutions exhibit artificial stiffening effects (see, e.g., [47]). Accurate results typically require substantial mesh refinement, which leads to a significantly more expensive computational problem. These numerical difficulties can be mitigated by adopting reduced kinematic formulations based on shell theory. When combined with appropriate numerical techniques to alleviate locking at the discrete level (see, e.g., [47, 48]), shell formulations considerably reduce the sensitivity of the method to wall thickness and provide a robust alternative to full 3D discretizations. In this work, we consider a thin-wall modeling of the left atrium by building on the reduced 3D-shell model proposed in [48], which consists of approximating the 3D displacement field by a quadratic expansion across the thickness, without requiring a plane-stress assumption (in contrast to standard shell formulations). As demonstrated in [48], this approach yields a well-posed shell model and is asymptotically consistent with classical shell theories in the vanishing-thickness limit.
In what follows, the mid-surface $\mathcal{S}$ of the atrial wall $\Omega$ is assumed to be parametrized by a smooth mapping $\phi: \omega \subset \mathbb{R}^2 \to \mathbb{R}^3$, such that $\mathcal{S} = \phi(\omega)$, with $(\xi^1, \xi^2) \in \omega$ denoting the curvilinear coordinates in the reference 2D domain $\omega$. The local geometry of $\mathcal{S}$ is characterized by its covariant basis vectors

$$a_i \stackrel{\text{def}}{=} \partial_{\xi^i} \phi \quad i = 1, 2,$$

which span the tangent plane to the mid-surface $\mathcal{S}$. Consequently, the unit normal vector field $a_3$ is defined as

$$a_3 \stackrel{\text{def}}{=} \frac{a_1 \times a_2}{\|a_1 \times a_2\|}.$$

The reference configuration of the left atrium $\Omega$ is finally assumed to be parametrized as

$$\Omega = \Phi(\omega \times (-1, 1)),$$

with the mapping $\Phi: \omega \times (-1, 1) \to \mathbb{R}^3$ given by the relation

$$\Phi(\xi^1, \xi^2, \xi^3) \stackrel{\text{def}}{=} \phi(\xi^1, \xi^2) + \xi^3 \frac{\epsilon(\xi^1, \xi^2)}{2} a_3(\xi^1, \xi^2), \quad (12)$$

where $\xi^3$ denotes the curvilinear coordinate oriented along the unit normal to $\mathcal{S}$, and $\epsilon: \omega \to \mathbb{R}^+$ represents the spatially varying wall thickness. A schematic representation of the mapping $\Phi$ is provided in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2: The geometrical mapping $\Phi$.</figcaption></figure>

The covariant basis vectors describing the full atrial wall geometry are given by $g_i \stackrel{\text{def}}{=} \partial_{\xi^i} \Phi$ with $i = 1, 2, 3$. Finally, the components of the 3D metric tensor in the reference configuration are given by $G_{ij} \stackrel{\text{def}}{=} g_i \cdot g_j$. The corresponding contravariant basis vectors $g^i$ are defined through the orthogonality condition $g^i \cdot g_j = \delta^i_j$, yielding the contravariant metric components $G^{ij} \stackrel{\text{def}}{=} g^i \cdot g^j$.
The fundamental idea for the 3D-shell modeling of the left atrium consists in performing a Galerkin projection of the variational formulation (11)$_1$ by considering the following reduced space (see [48]):

$$V^{\text{sh}} \stackrel{\text{def}}{=} \{ \mathbf{v}_s \in \mathbf{v} : \mathbf{v} = \mathbf{v}_0 + \mathbf{v}_1 \xi^3 + \mathbf{v}_2 (\xi^3)^2, \ \mathbf{v}_i \in [H^1(\mathcal{S})]^3, \ i = 0, 1, 2 \}, \quad (13)$$

where $\xi^3$ denotes the transverse coordinate introduced in (12) through the geometric mapping $\Phi$. Note that the vector fields $\mathbf{v}_i$, $i = 0, 1, 2$, are defined solely on the mid-surface $\mathcal{S}$ and therefore depend only on the tangential coordinates $\xi^1$ and $\xi^2$. Therefore, after appropriate composition with the mappings $\Phi$ and $\phi$, the virtual displacement identity in (13) should be read as

$$\mathbf{v}(\xi^1, \xi^2, \xi^3) = \mathbf{v}_0(\xi^1, \xi^2) + \mathbf{v}_1(\xi^1, \xi^2)\xi^3 + \mathbf{v}_2(\xi^1, \xi^2)(\xi^3)^2$$

for all $(\xi^1, \xi^2) \in \omega$ and $\xi^3 \in (-1, 1)$.

The proposed 3D-shell model of the left atrium reads therefore as follows:

$$\left\{ \begin{array}{l} \text{For } t > 0, \text{ find } \mathbf{y}(t) \in \mathbf{V}^{\text{sh}}, \text{ with } \dot{\mathbf{y}} = \partial_t \mathbf{y}, e_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}, k_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}^+, \\ \text{and } \tau_c : \Omega \times \mathbb{R}^+ \to \mathbb{R} \text{ such that (11) holds for all } \mathbf{v} \in \mathbf{V}^{\text{sh}}. \end{array} \right. \qquad (14)$$

It is worth noting that the only difference between this reduced problem and the full 3D formulation of the left atrium presented above lies in the kinematic assumption introduced through the reduced space (13). Additionally, since no plane-stress assumption is imposed in (14), the incorporation of fully 3D, non-quadratic constitutive law (such as the Demiray-Gasser law (1)) is straightforward, in contrast to what occurs with standard shell formulations.

In this framework, the usual mechanical quantities such as the deformation gradient admit a mixed-component representation

$$F_j^i = \delta_j^i + g^i \cdot \mathbf{y}_{,j},$$

where $\mathbf{y}_{,j}$ denotes the partial derivative of the displacement field with respect to the curvilinear coordinate $\xi^j$, and the identity term arises from the duality relation $g^i \cdot g_j = \delta_j^i$. Similarly, the covariant components of the Green-Lagrange and right Cauchy-Green strain tensor are respectively given by

$$E_{ij} = \frac{1}{2}(\mathbf{y}_{,i} \cdot \mathbf{g}_j + \mathbf{y}_{,j} \cdot \mathbf{g}_i + \mathbf{y}_{,i} \cdot \mathbf{y}_{,j}), \quad C_{ij} = G_{ij} + \mathbf{g}_i \cdot \mathbf{y}_{,j} + \mathbf{g}_j \cdot \mathbf{y}_{,i} + \mathbf{y}_{,i} \cdot \mathbf{y}_{,j}.$$

### 2.1.2 0D blood flow modeling in the left cardiac cavities

In order to model the interactions between cavities at a reduced computational cost, while retaining the essential physiological features, we adopt a 0D lumped-parameter description of the circulatory compartments connected to the left atrium. In this framework, each chamber or vascular segment has an associated pressure, while the
flows between compartments are determined through algebraic relations that mimic the effect of valves and vascular resistances.

## Blood flow in the left atrium

The exchange of blood between the pulmonary veins, the left atrium and the left ventricle is controlled by the mitral valve, which acts as a unidirectional flow regulator driven by the pressure gradient between the cavities. The valve is considered open when the downstream pressure is lower than the upstream pressure, and closed otherwise. Accordingly, the net flux exiting the left atrium is given by the following expression:

$$Q_{la} \stackrel{\text{def}}{=} Q_{vs}(p_{la}, p_{vs}) + Q_{mv}(p_{la}, p_{lv}), \quad (15)$$

with $Q_{vs}(p_{la}, p_{vs}) \stackrel{\text{def}}{=} k_{vs}(p_{la} - p_{vs})$ and

$$Q_{mv}(p_{lv}, p_{la}) \stackrel{\text{def}}{=} \begin{cases} 0 & \text{if } p_{lv} \ge p_{la}, \\ k_{mv}(p_{la} - p_{lv}) & \text{if } p_{lv} < p_{la}. \end{cases}$$

Here, $p_{vs}(t)$ and $p_{lv}(t)$ respectively denote the spatially homogeneous pressures in the pulmonary veins and in the left ventricle, whereas $k_{vs}$, $k_{mv} > 0$ stand for the pulmonary and mitral conductance coefficients. Mass flow in the left atrium is preserved by enforcing the zero net flux condition

$$- \int_{\Gamma_{en} \cup \Gamma_{mc}} J \dot{y} \cdot \mathbf{F}^{-\top} \mathbf{n} + Q_{la}(p_{vs}, p_{lv}, p_{la}) = 0,$$

where $\Gamma_{mc}$ stands for an ad hoc surface closure of the mitral valve. Note that, this introduces a scalar constraint in the weak formulation (14)$_1$, with the atrial pressure $p_{la}(t)$ acting as Lagrange multiplier. In this work, the pulmonary venous pressure $p_v(t)$ is assumed to be given, while the dynamics of $p_{lv}$ are described in the next paragraph using a 0D electromechanical model of the left ventricle.

## Reduced electromechanical model of the left ventricle

To capture the hemodynamic coupling between the left atrium and the left ventricle, the ventricular pressure $p_{lv}$ is described using the 0D reduced-order model proposed in [50]. This formulation provides a consistent reduction of the full 3D electromechanical problem (11) proposed in [23] and is adapted to an idealized geometry. The reduced model assumes spherical symmetry for both the geometry and the constitutive behavior. Under this hypothesis, the ventricle is represented as a spherical shell defined by a radius $R_0$ and wall thickness $d_0$ in the stress-free reference configuration $\Omega_{lv}$. Therefore, the associated displacement field of the left ventricle $\mathbf{y}^{\text{lv}}(t) \stackrel{\text{def}}{=} \mathbf{y}^{\text{lv}}\mathbf{i}_r = (R - R_0)\mathbf{i}_r$, is restricted to the radial direction $\mathbf{i}_r$, where $R$ denotes the current radius.

The myocardium is modeled as an incompressible material, enforcing the isochoric constraint $J = 1$. The stress on the idealized sphere can be represented by an effective spherical stress $\Sigma_{\text{sph}}$. This term arises from projecting the full stress tensor onto the
reduced radial kinematics (see [50] for derivation details) and incorporates active, passive, and viscous contributions such that

$$ \Sigma_{\text{sph}} = \sigma_a^{\text{lv}} + 4(1 - C^{-3}) \left( \frac{\partial W_e^{\text{lv}}}{\partial J_1} + C \frac{\partial W_e^{\text{lv}}}{\partial J_2} \right) + 2 \frac{\partial W_e^{\text{lv}}}{\partial J_4} + \eta \dot{C} (1 + 2C^{-6}). $$

Here, $C \stackrel{\text{def}}{=} (1+y/R_0)^2$ is the squared circumferential stretch, and $J_1, J_2$, and $J_4$ denote the standard reduced invariants of the right Cauchy-Green deformation tensor, with $J_2 \stackrel{\text{def}}{=} \frac{1}{2}((\text{tr}\mathbf{C}^{\text{lv}})^2 - \text{tr}(\mathbf{C}^{\text{lv}})^2)J^{-4/3}$, where the right Cauchy-Green deformation tensor of the reduced sphere is defined by

$$ \mathbf{C}^{\text{lv}} \stackrel{\text{def}}{=} \begin{pmatrix} C_{rr} & 0 & 0 \\ 0 & C & 0 \\ 0 & 0 & C \end{pmatrix}, $$

with $C_{rr} \stackrel{\text{def}}{=} C^{-2}$. The active contribution of the reduced sphere $\sigma_a^{\text{lv}}$ is defined similarly to the atrial tissue in (7) but assuming a circumferential fiber orientation aligned with the orthoradial direction $i_{\phi_1}$ of the sphere, such that $f^{\text{lv}} = i_{\phi_1}$. This results in the following strain along the fibers $e_f^{\text{lv}} \stackrel{\text{def}}{=} (f^{\text{lv}})^{\top} \mathbf{E}^{\text{lv}} f^{\text{lv}}$, where $\mathbf{E}^{\text{lv}} \stackrel{\text{def}}{=} 1/2(\mathbf{C}^{\text{lv}} - \mathbf{I})$ represents the Green-Lagrange strain tensor of the reduced sphere.

Although the general expression for $\Sigma_{\text{sph}}$ admits a dependence on the second invariant $J_2$, the specific strain-energy density function $W_e$ adopted in this work for the left ventricle relies exclusively on the isotropic invariant $J_1$ and the anisotropic invariant $J_4$. Therefore, following the formulation proposed in [50] and inspired from [57], the elastic potential is defined as

$$ W_e^{\text{lv}} = C_0 \exp C_1(J_1 - 3)^2 + C_2 \exp C_3(J_4 - 1)^2. $$

The parameters $C_0, C_1, C_2$ and $C_3$ are selected to reproduce a physiological response of the left ventricle.

Given the geometrical and kinematical assumptions, the terms of the dynamic balance equation (analogous to (11)$_1$) integrated over the reference domain are space-wise. After some algebraic manipulations, the complete 0D electromechanical formulation for the left ventricle read as follows: Find the radial displacement $y^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, the extension along the fibers $e_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, the active stiffness $k_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}^+$, the active tension $\tau_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$ and the internal pressure $p_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, such that
$$
\left\{
\begin{aligned}
& \rho d_0 \ddot{y}^{\text{lv}} + \frac{d_0}{R_0} \left(1 + \frac{y^{\text{lv}}}{R_0}\right) \Sigma_{\text{sph}} = p_{\text{lv}} \left(1 + \frac{y^{\text{lv}}}{R_0}\right)^2, \\
& \Sigma_{\text{sph}} = \sigma_a^{\text{lv}} + 4(1 - C^{-3}) \left(\frac{\partial W_e^{\text{lv}}}{\partial J_1} + C \frac{\partial W_e^{\text{lv}}}{\partial J_2}\right) + 2 \frac{\partial W_e^{\text{lv}}}{\partial J_4} + 2\eta \dot{C}(1 - 2C^{-6}), \\
& \sigma_a^{\text{lv}} = E_s^{\text{lv}} \frac{e_f^{\text{lv}} - e_c^{\text{lv}}}{(1 + 2e_c^{\text{lv}})}, \\
& \tau_c^{\text{lv}} + \mu^{\text{lv}} \dot{e}_c^{\text{lv}} = E_s \frac{(e_f^{\text{lv}} - e_c^{\text{lv}})(1 + 2e_f^{\text{lv}})}{(1 + 2e_c^{\text{lv}})^3}, \\
& \dot{k}_c^{\text{lv}} = -(|u_e^{\text{lv}}| + \alpha_d^{\text{lv}} |\dot{e}_c^{\text{lv}}|) k_c^{\text{lv}} + n_0 (e_c^{\text{lv}}) k_0^{\text{lv}} |u_e^{\text{lv}}|_+, \\
& \dot{\tau}_c^{\text{lv}} = -(|u_e^{\text{lv}}| + \alpha_d^{\text{lv}} |\dot{e}_c^{\text{lv}}|) \tau_c^{\text{lv}} + n_0 (e_c^{\text{lv}}) \sigma_0^{\text{lv}} |u_e^{\text{lv}}|_+ + k_c^{\text{lv}} \dot{e}_c^{\text{lv}}.
\end{aligned}
\right.
\qquad (16)
$$

In this formulation, all variables marked with the superscript $^{\text{lv}}$ denote quantities specific to the reduced-order electromechanical model of the left ventricle.

## Blood flow in the left ventricle

Blood flow in the left ventricle is regulated by the coordinated action of the mitral and aortic valves. Depending on the cardiac phase, namely, filling, isovolumetric contraction, ejection, or isovolumetric relaxation, the status of these valves (mitral-aortic) transitions sequentially: open-closed, closed-closed, closed-open, and closed-closed, respectively. Consequently, the total flux on the left ventricle, $Q_{\text{lv}}$ can be modeled as

$$ Q_{\text{lv}}(p_{\text{lv}}, p_{\text{ar}}, p_{\text{la}}) \stackrel{\text{def}}{=} \begin{cases} -Q_{\text{mv}}(p_{\text{la}}, p_{\text{lv}}) & \text{if } p_{\text{lv}} \le p_{\text{la}}, \\ 0 & \text{if } p_{\text{la}} \le p_{\text{lv}} \le p_{\text{ar}}, \\ k_{\text{ar}}(p_{\text{lv}} - p_{\text{ar}}) & \text{if } p_{\text{lv}} \ge p_{\text{ar}}, \end{cases} \qquad (17) $$

where $k_{\text{ar}} > 0$ is the arterialconductance, while $p_{\text{ar}}$ denotes the pressure in the aorta.

Mass conservation requires that the net flux entering or leaving the left ventricle balances the rate of volume change of the cavity. Expressed in terms of spherical shell kinematics, this dynamic continuity equation is defined as

$$ \dot{V} = 4\pi R_0^2 \left(1 + \frac{y^{\text{lv}}}{R_0}\right)^2 \dot{y}^{\text{lv}} = -Q_{\text{lv}}(p_{\text{lv}}, p_{\text{ar}}, p_{\text{la}}). $$

The aortic pressure $p_{\text{ar}}$ is determined using a lumped-parameter model that describes the systemic circulation downstream of the left ventricle. Specifically, we employ the two-element Windkessel model (see, e.g., [58]), which characterizes the external
circulatory dynamics as

$$ egin{cases} C_p \dot{p}_{\text{ar}} + \frac{p_{\text{ar}} - p_{\text{lv}}}{R_p} = |Q_{\text{lv}}(p_{\text{lv}}, p_{\text{ar}}, p_{\text{la}})|_+, \\ C_d \dot{p}_{\text{d}} + \frac{p_{\text{d}} - p_{\text{ar}}}{R_p} = \frac{p_{\text{v}} - p_{\text{d}}}{R_d}. \end{cases} \qquad (18) $$

where $C_p$ and $C_d$ denote the proximal and distal compliances, respectively, while $R_p$ and $R_d$ correspond to the proximal and distal resistances. The symbols $p_v$ and $p_d$ represent the venous and distal pressures, respectively.

### 2.1.3 Complete mathematical model

The complete formulation of the proposed 3D-shell electromechanical model of the left atrium read as follows: For $t > 0$, find $\mathbf{y}(t) \in \mathbf{V}^{\text{sh}}$, with $\dot{\mathbf{y}} = \partial_t \mathbf{y}$, $e_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}$, $k_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}^+$, $\tau_c : \Omega \times \mathbb{R}^+ \to \mathbb{R}$, $p_{\text{at}} : \mathbb{R}^+ \to \mathbb{R}$, $y^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, $e_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, $k_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}^+$, $\tau_c^{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$ and $p_{\text{lv}} : \mathbb{R}^+ \to \mathbb{R}$, such that

$$ egin{cases} 
\rho \int_{\Omega} \partial_t \dot{\mathbf{y}} \cdot \mathbf{v} + \int_{\Omega} \boldsymbol{\Sigma} : D_{\mathbf{y}} \mathbf{E} \cdot \mathbf{v} + \int_{\Gamma_{\text{mv}}} (a\mathbf{y} + b\dot{\mathbf{y}}) \cdot \mathbf{v} = - \int_{\Gamma_{\text{en}}} p_{\text{la}} J \mathbf{F}^{-\top} \mathbf{n} \cdot \mathbf{v} \quad \forall \mathbf{v} \in \mathbf{V}^{\text{sh}}, \\
\mu_c \dot{e}_c + \tau_c = E_s \frac{(\mathbf{f}^{\top} \mathbf{E} \mathbf{f} - e_c)(1 + 2\mathbf{f}^{\top} \mathbf{E} \mathbf{f})}{(1 + 2e_c)^3} \quad \text{in } \Omega, \\
\dot{k}_c = - (|u_e| + \alpha_d |\dot{e}_c|) k_c + n_0(e_c) k_0 |u_e|_+ \quad \text{in } \Omega, \\
\dot{\tau}_c = - (|u_e| + \alpha_d |\dot{e}_c|) \tau_c + \dot{e}_c k_c + n_0(e_c) \sigma_0 |u_e|_+ \quad \text{in } \Omega, \\
- \int_{\Gamma_{\text{en}} \cup \Gamma_{\text{mc}}} J \dot{\mathbf{y}} \cdot \mathbf{F}^{-\top} \mathbf{n} + Q(p_{\text{vs}}, p_{\text{lv}}, p_{\text{la}}) = 0, \\
\rho d_0 \ddot{y}^{\text{lv}} + \frac{d_0}{R_0} \left(1 + \frac{y^{\text{lv}}}{R_0}\right) \Sigma_{\text{sph}} = p_{\text{lv}} \left(1 + \frac{y^{\text{lv}}}{R_0}\right)^2, \\
\sigma_a^{\text{lv}} = E_s^{\text{lv}} \frac{e_f^{\text{lv}} - e_c^{\text{lv}}}{(1 + 2e_c^{\text{lv}})}, \\
\tau_c^{\text{lv}} + \mu^{\text{lv}} \dot{e}_c^{\text{lv}} = E_s^{\text{lv}} \frac{(e_f^{\text{lv}} - e_c^{\text{lv}})(1 + 2e_f^{\text{lv}})}{(1 + 2e_c^{\text{lv}})^3}, \\
\dot{k}_c^{\text{lv}} = - (|u_e^{\text{lv}}| + \alpha_d^{\text{lv}} |\dot{e}_c^{\text{lv}}|) k_c^{\text{lv}} + n_0(e_c^{\text{lv}}) k_0^{\text{lv}} |u_e^{\text{lv}}|_+, \\
\dot{\tau}_c^{\text{lv}} = - (|u_e^{\text{lv}}| + \alpha_d^{\text{lv}} |\dot{e}_c^{\text{lv}}|) \tau_c^{\text{lv}} + n_0(e_c^{\text{lv}}) \sigma_0^{\text{lv}} |u_e^{\text{lv}}|_+ + k_c^{\text{lv}} \dot{e}_c^{\text{lv}}, \\
4\pi R_0^2 \left(1 + \frac{y^{\text{lv}}}{R_0}\right)^2 \dot{y}^{\text{lv}} + Q_{\text{lv}}(p_{\text{lv}}, p_{\text{ar}}, p_{\text{la}}) = 0, \\
C_p \dot{p}_{\text{ar}} + \frac{p_{\text{ar}} - p_{\text{lv}}}{R_p} = |Q_{\text{lv}}(p_{\text{lv}}, p_{\text{ar}}, p_{\text{la}})|_+, \\
C_d \dot{p}_{\text{d}} + \frac{p_{\text{d}} - p_{\text{ar}}}{R_p} = \frac{p_v - p_d}{R_d}. 
\end{cases} \qquad (19) $$
This geometrically heterogeneous coupled non-linear system is schematically illustrated in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Illustration of the different components involved in the 3D-shell model of the left atrium (19).</figcaption></figure>

## 2.2 Numerical approximation

In this section, we detail the numerical approximation of the coupled electromechanical model (19). We first present the time semi-discrete formulation, using a modified midpoint scheme for the mechanical and hemodynamic equations, together with a backward Euler method for the active variables. We then introduce the spatial discretization of the 3D-shell atrial problem using $P_1 \otimes P_2$ prismatic finite elements.

### 2.2.1 Time discretization

The time interval $(0, T_f)$ is partitioned into $N$ sub-intervals $(t_n, t_{n+1})$, where $t_n \stackrel{\text{def}}{=} \tau n$, $T_f = t_N$, and $\tau > 0$ denotes the time-step length. In what follows, the following standard notation will be employed:

$$ \partial_{\tau} f^{n+1} \stackrel{\text{def}}{=} \frac{1}{\tau}(f^{n+1} - f^n), \quad f^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \frac{f^n + f^{n+1}}{2}, $$

for the first backward difference and the midpoint approximation, respectively.
The time discretization of the coupled electromechanical system (19) relies on a mixed mid-point/back-Euler scheme. Specifically, in the spirit of [23], the 3D-shell electromechanical model (14) and the reduced valve relations (15) are discretized in time using a mid-point approximation, while the internal variables governing active contraction are updated using a first-order Backward Euler scheme written in terms of the discrete unknowns

$$e_c^{n+1}, \quad \varrho^{n+1} \stackrel{\text{def}}{=} \tau_c^{n+1} / \sqrt{k_c^{n+1}}, \quad \gamma^{n+1} \stackrel{\text{def}}{=} \sqrt{k_c^{n+1}}. \qquad (20)$$

This time discretization yields an energy-preserving scheme (see [23]).

Similarly, for the reduced ventricular model (16), the momentum equation, the valve relations (17), and the circulatory subsystem (18) are discretized using the mid-point approximation. The corresponding active variables (16)$_{4,5,6}$ are directly discretized in time using a first-order Backward Euler approximation. In both systems, the contributions coming from the reduction factor $n_0$ are treated in an explicit fashion.

Finally, the coupling between the two chambers is discretized in time using staggered explicit scheme. The left-ventricular pressure involved in the volume constraint of 3D-shell atrial formulation is evaluated explicitly as $p_{\text{lv}}^{n-\frac{1}{2}}$, viz.,

$$\int_{\Gamma_{\text{en}} \cup \Gamma_{\text{mc}}} J^{n+\frac{1}{2}} \dot{\boldsymbol{y}}^{n+\frac{1}{2}} \cdot (\boldsymbol{F}^{n+\frac{1}{2}})^{-\top} \boldsymbol{n} = Q_{\text{la}} (p_{\text{vs}}(t_{n+\frac{1}{2}}), p_{\text{lv}}^{n-\frac{1}{2}}, p_{\text{la}}^{n+\frac{1}{2}}),$$

whereas the new atrial pressure $p_{\text{la}}^{n+1}$ obtained from the 3D-shell model is used as input for the 0D ventricular model, namely,

$$\left\{
\begin{array}{l}
4\pi R_0^2 \left( 1 + \frac{y^{\text{lv}, n+\frac{1}{2}}}{R_0} \right)^2 \dot{y}^{\text{lv}, n+\frac{1}{2}} + Q_{\text{lv}}(p_{\text{lv}}^{n+\frac{1}{2}}, p_{\text{ar}}^{n+\frac{1}{2}}, p_{\text{la}}^{n+\frac{1}{2}}) = 0, \\
C_p \partial_{\tau} p_{\text{ar}}^{n+1} + \frac{p_{\text{ar}}^{n+\frac{1}{2}} - p_{\text{lv}}^{n+\frac{1}{2}}}{R_p} = |Q_{\text{lv}}(p_{\text{lv}}^{n+\frac{1}{2}}, p_{\text{ar}}^{n+\frac{1}{2}}, p_{\text{la}}^{n+\frac{1}{2}})|_{+}.
\end{array}
\right.$$

This staggered solution procedure avoids the computational cost of a fully monolithic approach and has been found numerically stable for all the regimes considered (see Section 3). The resulting time semi-discrete approximation of (19) reads therefore as follows: For $n \ge 0$,
1. 3D-shell atrial electromechanics: Find $(\mathbf{y}^{n+1}, \dot{\mathbf{y}}^{n+1}, e_c^{n+1}, \varrho^{n+1}, \gamma^{n+1}, p_{la}^{n+1}) \in \mathbf{V}^{\text{sh}} \times \mathbf{V}^{\text{sh}} \times \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{+} \times \mathbb{R}$, with $\dot{\mathbf{y}}^{n+\frac{1}{2}} = \partial_{\tau} \mathbf{y}^{n+1}$, such that

$$
\left\{
\begin{aligned}
& \rho_s \int_{\Omega} \partial_{\tau} \dot{\mathbf{y}}^{n+1} \cdot \mathbf{v} + \int_{\Omega} \mathbf{\Sigma}^{n+\frac{1}{2}} : D_y \mathbf{E}^{n+\frac{1}{2}} \cdot \mathbf{v} + \int_{\Gamma_{mv}} (a \mathbf{y}^{n+\frac{1}{2}} + b \dot{\mathbf{y}}^{n+\frac{1}{2}}) \cdot \mathbf{v} \\
& = - \int_{\Gamma_{en}} J^{n+\frac{1}{2}} p_{la}^{n+\frac{1}{2}} (\mathbf{F}^{n+\frac{1}{2}})^{-\top} \mathbf{n} \cdot \mathbf{v} \quad \forall \mathbf{v} \in \mathbf{V}^{\text{sh}}, \\
& \mathbf{\Sigma}^{n+\frac{1}{2}} = \mathbf{\Sigma}_{\text{p}}^{n+\frac{1}{2}} + \frac{E_s (\mathbf{f}^{\top} \mathbf{E}^{n+\frac{1}{2}} \mathbf{f} - e_c^{n+\frac{1}{2}})}{(1 + 2e_c^{n+\frac{1}{2}})^2} \mathbf{f} \otimes \mathbf{f}, \\
& \gamma^{n+1} \varrho^{n+1} + \mu_c \partial_{\tau} e_c^{n+1} = E_s \frac{(\mathbf{f}^{\top} \mathbf{E}^{n+\frac{1}{2}} \mathbf{f} - e_c^{n+\frac{1}{2}})(1 + 2\mathbf{f}^{\top} \mathbf{E}^{n+\frac{1}{2}} \mathbf{f})}{(1 + 2e_c^{n+\frac{1}{2}})^3}, \quad (21) \\
& \partial_{\tau}(\gamma^{n+1})^2 = - (|u_e^{n+1}| + \alpha_d |\partial_{\tau} e_c^{n+1}|)(\gamma^{n+1})^2 + n_0(e_c^n)k_0 |u_e^{n+1}|_+, \\
& \partial_{\tau} \varrho^{n+1} = \gamma^{n+1} \partial_{\tau} e_c^{n+1} + \frac{n_0(e_c^n)\sigma_0 |u_{e,l}^{n+1}|_+}{\gamma^{n+1}} \left(1 - \frac{k_0}{2\sigma_0\gamma^{n+1}} \varrho^{n+1}\right) \\
& \qquad - \frac{1}{2} (|u_e^{n+\frac{1}{2}}| + |\partial_{\tau} e_c^{n+1}|) \varrho^{n+1}, \\
& \int_{\Gamma_{en} \cup \Gamma_{mc}} J^{n+\frac{1}{2}} \dot{\mathbf{y}}^{n+\frac{1}{2}} \cdot (\mathbf{F}^{n+\frac{1}{2}})^{-\top} \mathbf{n} = Q_{la} (p_{vs}(t_{n+\frac{1}{2}}), p_{lv}^{n-\frac{1}{2}}, p_{la}^{n+\frac{1}{2}}),
\end{aligned}
\right.
$$
2. 0D ventricular electromechanics: Find $(y^{\text{lv},n+1}, \dot{y}^{\text{lv},n+1}, e_c^{\text{lv},n+1}, k_c^{\text{lv},n+1}, \tau_c^{\text{lv},n+1}, p_{\text{lv}}^{n+1}) \in \mathbb{R}^6$ with $\dot{y}^{\text{lv},n+1/2} = \partial_{\tau} y^{\text{lv},n+1}$, such that

$$
\left\{
\begin{aligned}
& \rho d_0 \partial_{\tau} \dot{y}^{\text{lv}, n+1} + \frac{d_0}{R_0} \left( 1 + \frac{y^{\text{lv}, n+1/2}}{R_0} \right) \Sigma_{\text{sph}}^{n+1/2} = p_{\text{lv}}^{n+1/2} \left( 1 + \frac{y^{\text{lv}, n+1/2}}{R_0} \right)^2, \\
& \Sigma_{\text{sph}}^{n+1/2} = \sigma_{\text{a}}^{\text{lv}, n+1/2} + 4(1 - (C^{n+1/2})^{-3}) \left( \frac{\partial W_{\text{e}}^{\text{lv}, n+1/2}}{\partial J_1} + C^{n+1/2} \frac{\partial W_{\text{e}}^{\text{lv}, n+1/2}}{\partial J_2} \right) \\
& \qquad + 2 \frac{\partial W_{\text{e}}^{\text{lv}, n+1/2}}{\partial J_4} + 2\eta \dot{C}^{n+1/2} (1 - 2(C^{n+1/2})^{-6}), \\
& \sigma_{\text{a}}^{\text{lv}, n+1/2} = E_{\text{s}}^{\text{lv}} \frac{e_{\text{f}}^{\text{lv}, n+1/2} - e_{\text{c}}^{\text{lv}, n+1/2}}{(1 + 2e_{\text{c}}^{\text{lv}, n+1/2})}, \\
& \tau_{\text{c}}^{\text{lv}, n+1} + \mu^{\text{lv}} \partial_{\tau} e_{\text{c}}^{\text{lv}, n+1} = E_{\text{s}}^{\text{lv}} \frac{(e_{\text{f}}^{\text{lv}, n+1} - e_{\text{c}}^{\text{lv}, n+1})(1 + 2e_{\text{f}}^{\text{lv}, n+1})}{(1 + 2e_{\text{c}}^{\text{lv}, n+1})^3}, \\
& \partial_{\tau} k_{\text{c}}^{\text{lv}, n+1} = -(|u_{\text{e}}^{\text{lv}}(t_{n+1})| + \alpha_{\text{d}}^{\text{lv}} |\partial_{\tau} e_{\text{c}}^{\text{lv}, n+1}|) k_{\text{c}}^{\text{lv}, n+1} + n_0 (e_{\text{c}}^{\text{lv}, n}) k_0^{\text{lv}} |u_{\text{e}}^{\text{lv}}(t_{n+1})|_+, \\
& \partial_{\tau} \tau_{\text{c}}^{\text{lv}, n+1} = -(|u_{\text{e}}^{\text{lv}}(t_{n+1})| + \alpha_{\text{d}}^{\text{lv}} |\partial_{\tau} e_{\text{c}}^{\text{lv}, n+1}|) \tau_{\text{c}}^{\text{lv}, n+1} + n_0 (e_{\text{c}}^{\text{lv}, n}) \sigma_0^{\text{lv}} |u_{\text{e}}^{\text{lv}}(t_{n+1})|_+ \\
& \qquad + k_{\text{c}}^{\text{lv}, n+1} \dot{e}_{\text{c}}^{\text{lv}, n+1}, \\
& 4\pi R_0^2 \left( 1 + \frac{y^{\text{lv}, n+1/2}}{R_0} \right)^2 \dot{y}^{\text{lv}, n+1/2} + Q_{\text{lv}}(p_{\text{lv}}^{n+1/2}, p_{\text{ar}}^{n+1/2}, p_{\text{la}}^{n+1/2}) = 0, \\
& C_{\text{p}} \partial_{\tau} p_{\text{ar}}^{n+1} + \frac{p_{\text{ar}}^{n+1/2} - p_{\text{lv}}^{n+1/2}}{R_{\text{p}}} = |Q_{\text{lv}}(p_{\text{lv}}^{n+1/2}, p_{\text{ar}}^{n+1/2}, p_{\text{la}}^{n+1/2})|_+, \\
& C_{\text{d}} \partial_{\tau} p_{\text{d}}^{n+1} + \frac{p_{\text{d}}^{n+1/2} - p_{\text{ar}}^{n+1/2}}{R_{\text{p}}} = \frac{p_{\text{v}}(t_{n+1/2}) - p_{\text{d}}^{n+1/2}}{R_{\text{d}}}. \end{aligned}
\right.
\tag{22}
$$

### 2.2.2 Fully discrete method

Having specified the time discretization and the explicit coupling strategy in (21)-(22), we now introduce the spatial discretization used for the 3D-shell atrial sub-problem (21). The reference domain $\Omega$ is discretized using a computational mesh $\mathcal{T}_h$ composed of triangular prismatic elements (see Fig. 4). The construction of this mesh relies on the following geometrical steps:

* An approximation of the atrial mid-surface $\mathcal{S}$, represented by a triangulation $\mathcal{T}_h^S$ made of linear triangles;

* An approximation of the unit normal to $\mathcal{S}$ at each vertex of $\mathcal{T}_h^S$;

* An approximation of the local atrial wall thickness associated with each vertex of $\mathcal{T}_h^S$.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Left: Prismatic mesh of the left-atrial wall Ω. Right: The map Φh,K between the reference and current element of the mesh Th.</figcaption></figure>

Let $\hat{K} \stackrel{\text{def}}{=} \hat{T} \times (-1, 1)$ be the reference unit prism, where $\hat{T}$ stands for the reference unit triangle in $\mathbb{R}^2$. Each element $K$ of the computational prismatic mesh $\mathcal{T}_h$ is obtained as the image of the local mapping $\Phi_{h,K}: \hat{K} \to K \subset \mathbb{R}^3$ defined by

$$ \Phi_{h,K}(r, s, z) \stackrel{\text{def}}{=} \sum_{i=1}^{3} \lambda_i(r, s) \left( \mathbf{x}_T^{(i)} + z \frac{\epsilon_K^{(i)}}{2} \mathbf{a}_{3,T}^{(i)} \right) \quad \forall (r, s, z) \in \hat{K}, \quad (23) $$

where $\{\lambda_i\}_{i=1}^3$ stand for the 2D affine shape functions of the reference element $\hat{T}$, $\{\mathbf{x}_T^{(i)}\}_{i=1}^3$ denote the positions of the nodes of the mid-surface triangle $T \in \mathcal{T}_h^S$ associated to the prismatic element $K \in \mathcal{T}_h$ and $\{\mathbf{a}_{3,T}^{(i)}\}_{i=1}^3$, $\{\epsilon_T^{(i)}\}_{i=1}^3$ respectively represent the normal vectors to the mid-surface $S$ and the thickness associated to each node of the triangle $T$. All these geometrical ingredients are illustrated in Fig. 4.

It is worth noting that the mapping $\Phi_{h,K}$ defined in (23) provides a piecewise parametrization of the computational domain, where $(r, s)$ represent the local coordinates tangential to the mid-surface and $z$ denotes the coordinate in the direction normal to it. This construction can be regarded as the discrete analogue of the continuous mapping (12). However, it offers a practical advantage, instead of requiring an explicit parametrization $\phi$ of the mid-surface $S$, it only relies on its discrete approximation given by the triangulation $\mathcal{T}_h^S$ together with the nodal normal vectors $\mathbf{a}_{3,T}^{(i)}$ and the nodal thickness values $\epsilon_T^{(i)}$.

The 3D-shell electromechanical formulation of the left atrium (21) is spatially discretized using the following finite element approximation of the space $V^{\text{sh}}$:

$$ V_h^{\text{sh}} \stackrel{\text{def}}{=} \{ v_h \in V : v_h \circ \Phi_{h,K} = \sum_{i=1}^{3} \lambda_i(r, s) \left( v_0^{(i)} + v_1^{(i)}z + v_2^{(i)}(z)^2 \right) \quad \forall K \in \mathcal{T}_h, $$
$$ v_j^{(i)} \in \mathbb{R}^3, \quad i = 1, 2, 3, \quad j = 0, 1, 2 \}, \quad (24) $$

which corresponds to the tensorization of a $\mathbb{P}_1$ approximation in the in-plane coordinates $(r, s)$ with a $\mathbb{P}_2$ approximation along the transverse coordinate $z$. Therefore, the computer implementation of this 3D-shell element is analogous to that of a standard 3D finite element. The set $\{v_j^{(i)}\}_{j=0,1,2}^{i=1,2,3}$ in (24) represents the 27 degrees of freedom of the displacement in each prismatic element $K \in \mathcal{T}_h$. Furthermore, as outlined in [48], for each $v_h \in \mathbf{V}_h^{\text{sh}}$ the following equivalent representation holds

$$ v_h \circ \Phi_{h,K} = \sum_{i=1}^{3} \lambda_i(r, s) \left[ \frac{(z-1)z}{2} v_{\text{bot}}^{(i)} + (1-(z)^2) v_{\text{mid}}^{(i)} + \frac{(z+1)z}{2} v_{\text{upp}}^{(i)} \right], $$

where the subscripts `bot`, `mid`, `upp` denote the displacement degrees of freedom associated with the bottom, middle, and upper triangles of the prismatic element (see Fig. 4).

A Mixed Interpolation of Tensorial Components (MITC) technique is employed to treat the transverse strain components, as required in shell formulations to avoid numerical locking induced by asymptotic constraints on the discrete solution (see [47, 48, 59, 60]). Specifically, the MITC3 formulation is applied within each triangular layer of the prismatic element [60–62]. In this approach, the in-layer strain components

$$ \tilde{E}_{13}, \quad \tilde{E}_{23}, \quad \tilde{E}_{q3} \stackrel{\text{def}}{=} \tilde{E}_{23} - \tilde{E}_{13}, $$

are evaluated at the quadrature points by reinterpolation from the associated tying points (denoted by the tilde $\tilde{\cdot}$), while the in-plane strain components are left unmodified. This treatment is standard in 3D-shell discretizations and is necessary to obtain a mechanically consistent and thickness-robust response.

The contractile unknowns $e_c^{n+1}$, $\varrho^{n+1}$ and $\gamma^{n+1}$ in (20) are discretized in space using a collocation scheme at each quadrature point $q_l$, $l = 1, \dots, N_{q,h}$, where $N_{q,h} \stackrel{\text{def}}{=} |\mathcal{T}_h|n_q$ and $n_q$ denotes the number of quadrature points in each element $K$ of the mesh $\mathcal{T}_h$. This leads to $N_{q,h}$ degrees of freedom for each active variable, collected in the arrays:

$$ e_{c,h}^{n+1} \stackrel{\text{def}}{=} \{e_{c,l}^{n+1}\}_{l=1}^{N_{q,h}}, \quad \varrho_h^{n+1} \stackrel{\text{def}}{=} \{\varrho_l^{n+1}\}_{l=1}^{N_{q,h}}, \quad \gamma_h^{n+1} = \{\gamma_l^{n+1}\}_{l=1}^{N_{q,h}}. $$

The proposed fully discrete approximation of the electromechanical problem (19) is summarized in Algorithm 1. The subscript $l$ in the quantities $f_l$, $E_l^{n+\frac{1}{2}}$, $\tilde{\Sigma}_{p,l}^{n+\frac{1}{2}}$, $\tilde{\Sigma}_l^{n+\frac{1}{2}}$, $\Sigma_a^{n+\frac{1}{2}}$ indicates their evaluation at the quadrature points $q_l$, $l = 1, \dots, N_{q,h}$. It is also noted that the passive stress $\tilde{\Sigma}_{p,l}^{n+\frac{1}{2}}$ depends on the re-interpolated Green-Lagrange strain tensor $\tilde{\Sigma}_l^{n+\frac{1}{2}}$. Furthermore, the discrete active variables $e_{c,h}^{n+1}$, $\varrho_h^{n+1}$ and $\gamma_h^{n+1}$ can be eliminated at the element level through a Schur complement procedure (see [63]).
**Algorithm 1** Fully discrete approximation of the 3D-shell electromechanical model (19).

For $n \ge 0$:

1. Find $(\mathbf{y}_h^{n+1}, \dot{\mathbf{y}}_h^{n+1}, e_{c,h}^{n+1}, \varrho_h^{n+1}, \gamma_h^{n+1}, p_{at}^{n+1}) \in \mathbf{V}_h^{\text{sh}} \times \mathbf{V}_h^{\text{sh}} \times \mathbb{R}^{N_{q,h}} \times \mathbb{R}^{N_{q,h}} \times (\mathbb{R}^+)^{N_{q,h}} \times \mathbb{R}$ with $\dot{\mathbf{y}}_h^{n+\frac{1}{2}} = \partial_{\tau} \mathbf{y}_h^{n+1}$ and such that

$$
\left\{
\begin{aligned}
& \rho_s \int_{\Omega} \partial_{\tau} \mathbf{u}_h^{n+1} \cdot \mathbf{v}_h + \int_{\Omega} \tilde{\mathbf{\Sigma}}_l^{n+\frac{1}{2}} : D_y \tilde{\mathbf{E}}^{n+\frac{1}{2}} \cdot \mathbf{v}_h + \int_{\Omega} \mathbf{\Sigma}_a^{n+\frac{1}{2}} : D_y \mathbf{E}^{n+\frac{1}{2}} \cdot \mathbf{v}_h \\
& + \int_{\Gamma_{mv}} (a \mathbf{y}_h^{n+\frac{1}{2}} + b \dot{\mathbf{y}}_h^{n+\frac{1}{2}}) \cdot \mathbf{v}_h = - \int_{\Gamma_{en}} J^{n+\frac{1}{2}} p_{at}^{n+\frac{1}{2}} (\mathbf{F}^{n+\frac{1}{2}})^{-\top} \mathbf{n} \cdot \mathbf{v}_h, \\
& \tilde{\mathbf{\Sigma}}_l^{n+\frac{1}{2}} = \tilde{\mathbf{\Sigma}}_{p,l}^{n+\frac{1}{2}} + \mathbf{\Sigma}_{a,l}^{n+\frac{1}{2}}, \quad \mathbf{\Sigma}_{a,l}^{n+\frac{1}{2}} = \frac{E_s (\mathbf{f}_l^{\top} \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{f}_l - e_{c,l}^{n+\frac{1}{2}})}{(1 + 2e_{c,l}^{n+\frac{1}{2}})^2} \mathbf{f}_l \otimes \mathbf{f}_l, \\
& \gamma_l^{n+1} \varrho_l^{n+1} + \mu_c \partial_{\tau} e_{c,l}^{n+1} = E_s \frac{(\mathbf{f}_l^{\top} \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{f}_l - e_{c,l}^{n+\frac{1}{2}})(1 + 2\mathbf{f}_l^{\top} \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{f}_l)}{(1 + 2e_{c,l}^{n+\frac{1}{2}})^3}, \quad (25) \\
& \partial_{\tau} (\gamma_l^{n+1})^2 = - (|u_{e,l}^{n+1}| + \alpha_d |\partial_{\tau} e_{c,l}^{n+1}|) (\gamma_l^{n+1})^2 + n_0 (e_c^n) k_0 |u_{e,l}^{n+1}|_+, \\
& \partial_{\tau} \varrho_l^{n+1} = \gamma_l^{n+1} \partial_{\tau} e_{c,l}^{n+1} + \frac{n_0 (e_c^n) \sigma_0 |u_{e,l}^{n+1}|_+}{\gamma_l^{n+1}} \left(1 - \frac{k_0}{2\sigma_0 \gamma_l^{n+1}} \varrho_l^{n+1}\right) \\
& \qquad - \frac{1}{2} (|u_{e,l}^{n+\frac{1}{2}}| + |\partial_{\tau} e_{c,l}^{n+1}|) \varrho_l^{n+1} = 0, \\
& \int_{\Gamma_{en} \cup \Gamma_{mc}} J^{n+\frac{1}{2}} \dot{\mathbf{y}}_h^{n+\frac{1}{2}} \cdot (\mathbf{F}^{n+\frac{1}{2}})^{-\top} \mathbf{n} = Q_{la} (p_{vs}(t_{n+\frac{1}{2}}), p_{lv}^{n-\frac{1}{2}}, p_{at,l}^{n+\frac{1}{2}}),
\end{aligned}
\right.
$$

for all $\mathbf{v}_h \in \mathbf{V}_h^{\text{sh}}$ and $l = 1, \dots, N_{q,h}$.

2. Find $(\mathbf{y}^{\text{lv},n+1}, \dot{\mathbf{y}}^{\text{lv},n+1}, e_c^{\text{lv},n+1}, k_c^{\text{lv},n+1}, \tau_c^{\text{lv},n+1}, p_{\text{lv}}^{n+1}) \in \mathbb{R}^6$ with $\dot{\mathbf{y}}^{\text{lv},n+\frac{1}{2}} = \partial_{\tau} \mathbf{y}^{\text{lv},n+1}$, such that (22) holds.

## 3 Results

This section presents the numerical results of the electromechanical simulations of the left atrium in physiological (Section 3.1) and pathological (Section 3.2) conditions. A series of numerical experiments has been carried out with the purpose of assessing the physiological relevance and computational benefits of the 3D-shell electromechanical model (19). To this purpose, the results obtained with the numerical method reported in Algorithm 1 are compared with those obtained with a standard $\mathbb{P}_1$ finite element approximation of the full 3D formulation (11). To ensure a rigorous comparison, constitutive electromechanical parameters and boundary conditions remain identical across all the 3D-shell and full 3D simulations.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Geometry and fiber architecture of the left-atrium.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: The three types of spatial grids used in the numerical experiments. Left: $P_1 \otimes P_2$ (nine-nodes) prismatic elements (quadrilateral elements on the later faces). Middle: linear tetrahedral mesh with two layer of elements. Right: linear tetrahedral mesh with one layer of elements.</figcaption></figure>

<table><caption>Table 1: Mesh details, execution times, memory and entity count per test case.</caption><thead><tr><th></th><th>3D-shell</th><th>3D (2-layer)</th><th>3D (1-layer)</th></tr></thead><tbody><tr><td>Number of nodes</td><td>59,574</td><td>59,574</td><td>39,716</td></tr><tr><td>Number of cells</td><td>39,131</td><td>234,786</td><td>117,393</td></tr><tr><td>Hours per cycle</td><td>2.25</td><td>3.7</td><td>2.11</td></tr><tr><td>Total compute hours</td><td>15.8</td><td>26.06</td><td>14.77</td></tr><tr><td>Memory per MPI process</td><td>300 MB</td><td>600 MB</td><td>350 MB</td></tr><tr><td>Avg. Newton iterations</td><td>4</td><td>6</td><td>5</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Imposed data curves for the pulmonary veins pressure, the electrical activation for the left-ventricle model, and the profile for $n_0$.</figcaption></figure>

The considered left atrial geometry was extracted from the Zygote heart model<fnref n="1" />, from which the internal and external surfaces where smoothed and refined (see Fig. 5a). Specifically, localized smoothing was applied to the endocardial and epicardial surfaces of the left atrial appendage (LAA) and the junction between the LAA and the pulmonary veins, to prevent the generation of sharp elements. These geometric treatments were performed using the 3-Matic<fnref n="2" /> software. The final computational mesh consists of a single layer of triangular prismatic elements, generated by extruding the epicardial surface along its local inward normal direction. The spatial heterogeneity of the atrial wall was captured by defining the local thickness as the point-wise Euclidean distance between the reference epicardial and endocardial surfaces. This distance is computed using the closest point projection algorithm (see, e.g., [64]) and is explicitly imposed during the mesh generation process to ensure a geometrically accurate reconstruction. Since anatomical geometries such as the Zygote model typically correspond to a loaded state, a pre-stressing procedure is employed to recover an approximate stress-free reference configuration of the left atrium (see, e.g., [65, 66]), ensuring that the application of physiological loads reproduces the target geometry with a consistent

<footnote n="1">¹https://www.zygote.com</footnote>

<footnote n="2">²https://www materialise.com</footnote>
<table><caption>**Table 2:** Parameters of the lumped-parameter left-ventricle model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>R<sub>0</sub></td><td>Reference LV radius</td><td>2.6 × 10<sup>-2</sup></td><td>m</td></tr><tr><td>d<sub>0</sub></td><td>Reference wall thickness</td><td>1.7 × 10<sup>-2</sup></td><td>m</td></tr><tr><td>ρ</td><td>Myocardial density</td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>η</td><td>Viscous coefficient in Σ<sub>sph</sub></td><td>70</td><td>–</td></tr><tr><td>C<sub>0</sub></td><td>Passive stiffness parameter (isotropic)</td><td>1.9 × 10<sup>3</sup></td><td>Pa</td></tr><tr><td>C<sub>1</sub></td><td>Passive exponent parameter (isotropic)</td><td>1.1 × 10<sup>-1</sup></td><td>–</td></tr><tr><td>C<sub>2</sub></td><td>Passive stiffness parameter (anisotropic)</td><td>1.9 × 10<sup>3</sup></td><td>Pa</td></tr><tr><td>C<sub>3</sub></td><td>Passive exponent parameter (anisotropic)</td><td>1.1 × 10<sup>-1</sup></td><td>–</td></tr><tr><td>E<sup>lv</sup><sub>s</sub></td><td>Active stiffness modulus</td><td>3 × 10<sup>5</sup></td><td>Pa</td></tr><tr><td>k<sub>ar</sub></td><td>Aortic valve conductance</td><td>1.3 × 10<sup>-5</sup></td><td>m<sup>3</sup>/Pa·s</td></tr><tr><td>R<sub>p</sub></td><td>Proximal resistance</td><td>8 × 10<sup>6</sup></td><td>Pa·s/m<sup>3</sup></td></tr><tr><td>R<sub>d</sub></td><td>Distal resistance</td><td>1 × 10<sup>8</sup></td><td>Pa·s/m<sup>3</sup></td></tr><tr><td>C<sub>p</sub></td><td>Proximal compliance</td><td>5 × 10<sup>-9</sup></td><td>m<sup>3</sup>/Pa</td></tr><tr><td>C<sub>d</sub></td><td>Distal compliance</td><td>1 × 10<sup>-8</sup></td><td>m<sup>3</sup>/Pa</td></tr><tr><td>t<sub>d</sub></td><td>Activation delay</td><td>130</td><td>ms</td></tr><tr><td>t<sub>dep</sub></td><td>Depolarization duration</td><td>11</td><td>ms</td></tr><tr><td>t<sub>plat</sub></td><td>Plateau duration</td><td>140</td><td>ms</td></tr><tr><td>t<sub>rep</sub></td><td>Repolarization duration</td><td>80</td><td>ms</td></tr><tr><td>u<sub>max</sub></td><td>Maximum activation potential</td><td>35</td><td>mV</td></tr><tr><td>u<sub>min</sub></td><td>Minimum activation potential</td><td>-12</td><td>mV</td></tr><tr><td>c<sub>w</sub></td><td>Conduction wave velocity</td><td>0.6</td><td>m/s</td></tr></tbody></table>

stress state, resulting in the final mesh reported in Fig. 5a. Finally, the planar surface closure $\Gamma_{mc}$ is triangulated without node insertion, using only the inner edges of $\Gamma_{mv}$.

The fiber architecture (see Fig. 5b) is subsequently generated on the computational mesh by using an adaptation of the Laplace-Dirichlet Rule-Based (LDRB) framework proposed in [16]. Unlike the standard volumetric approach, the algorithm is applied separately to the epicardial and endocardial surfaces to capture the significant transmural variations in fiber orientation [67, 68], with the wall-transverse distribution interpolated according to [69] (see Appendix B).

In this work, three spatial discretizations of the atrial wall are evaluated and compared (see Fig. 6):

* the prismatic $P_1 \otimes P_2$ approximation provided by the 3D-shell formulation of Algorithm 1;

* a standard $P_1$ approximation of the full 3D formulation (11) using single-layer tetrahedral mesh;

* a standard $P_1$ approximation of the full 3D formulation (11) using a two-layer tetrahedral mesh.

To guarantee a consistent comparison between the shell and solid models, both tetrahedral meshes are derived from the same base 9-node prismatic mesh. The two-layer configuration is generated by subdividing each prism through its thickness (see, e.g.,
<table><caption>**Table 3:** Constitutive parameters of the 3D-shell electromechanical model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>a<sub>p</sub></td><td>Passive matrix stiffness parameter</td><td>6 &middot; 10<sup>3</sup></td><td>Pa</td></tr><tr><td>b<sub>p</sub></td><td>Passive matrix exponential coefficient</td><td>5.5</td><td>–</td></tr><tr><td>a<sub>f</sub></td><td>Fiber stiffness parameter</td><td>12 &middot; 10<sup>3</sup></td><td>Pa</td></tr><tr><td>b<sub>f</sub></td><td>Fiber exponential coefficient</td><td>11.75</td><td>–</td></tr><tr><td>&kappa;</td><td>Bulk modulus penalization parameter</td><td>4 &middot; 10<sup>5</sup></td><td>Pa</td></tr><tr><td>&nu;</td><td>Passive viscosity coefficient</td><td>1000</td><td>N &middot; s &middot; m<sup>-2</sup></td></tr><tr><td>E<sub>s</sub></td><td>Elastic stiffness of sarcomere element</td><td>5 &middot; 10<sup>4</sup></td><td>Pa</td></tr><tr><td>k<sub>0</sub></td><td>Maximum active stiffness</td><td>1 &middot; 10<sup>5</sup></td><td>Pa</td></tr><tr><td>&sigma;<sub>0</sub></td><td>Intrinsic active stress (contractility constant)</td><td>2 &middot; 10<sup>6</sup></td><td>Pa</td></tr><tr><td>&mu;<sub>c</sub></td><td>Active viscous coefficient (contractile branch)</td><td>100</td><td>N &middot; s &middot; m<sup>-1</sup></td></tr><tr><td>k<sub>mv</sub></td><td>Conductance of mitral valve</td><td>8 &middot; 10<sup>-7</sup></td><td>m<sup>3</sup> &middot; s<sup>-1</sup> &middot; Pa<sup>-1</sup></td></tr><tr><td>k<sub>vs</sub></td><td>Conductance of pulmonary veins</td><td>5 &middot; 10<sup>-7</sup></td><td>m<sup>3</sup> &middot; s<sup>-1</sup> &middot; Pa<sup>-1</sup></td></tr><tr><td>&alpha;<sub>d</sub></td><td>Cross-bridge destruction coefficient</td><td>3.0</td><td>–</td></tr><tr><td>&alpha;</td><td>Activation scaling for V<sub>m</sub></td><td>0</td><td>–</td></tr><tr><td>&beta;</td><td>Activation offset for V<sub>m</sub></td><td>1</td><td>–</td></tr><tr><td>a</td><td>Viscoelastic displacement parameter</td><td>1 &middot; 10<sup>5</sup></td><td>–</td></tr><tr><td>b</td><td>Viscoelastic velocity parameter</td><td>1 &middot; 10<sup>7</sup></td><td>–</td></tr></tbody></table>

[70]), while the single-layer mesh is constructed by first reducing the prism to a 6-node wedge prior to subdivision. This strategy ensures that all discretizations share an identical mid-surface and boundary triangulation, providing a coherent framework for analyzing their mechanical responses. The resulting element and node counts are detailed in Table 1.

All simulations employ the same passive and active material parameters and the same mechanical boundary conditions and have been run at a time resolution of $\tau = 10^{-3}$ s. The constitutive parameters of the Demiray-Gasser law, the viscous regularization, and the active contraction model are identical in all test cases (see Table 3). The imposed analytic data curves for the pulmonary veins pressure of the electromechanical model (14), the imposed electrical activation of the left ventricle model (16), and the scaling factor profile $n_0$ relative to the fiber elongation are reported in Fig. 7.

Each simulation is run over seven cardiac cycles to eliminate initial transients and to reach a periodic solution. The electromechanical coupling strategy, the time-integration scheme, as well as the nonlinear solvers and tolerances are kept fixed throughout, thereby providing a controlled setting for analyzing the influence of the spatial discretization on the predicted deformation and computational cost.

All simulations were performed using the open-source finite element library FELiScE<fnref n="3" /> on a 2023 Mac Studio equipped with 128 GB of RAM and an Apple M2 Ultra

<footnote n="3">³https://gitlab.inria.fr/felisce/felisce</footnote>
processor, using 12 cores. The resulting execution times, memory usage, discretiza-
tion sizes and average number of Newton iterations per timestep for each mesh are
reported in Table 1.

**3.1 Normal condition**

In this section, we present the numerical results of the electromechanical simulations
of the left atrium under physiological conditions.

Electrical propagation

The transmembrane potential $V_m$ involved in the electrical input (2) is obtained by using the bidomain model (see Appendix A for further details). The atrial tissue is assumed to be initially repolarized, with $V_m(\cdot, 0) = -81.2$ mV in $\Omega$. Electrical activation is initiated at Bachmann's bundle. Fig. 8 reports some snapshots of the

<figure><img src="image_9.png" /><figcaption>**Fig. 8:** Left atrium transmembrane potential $V_m$ used as electrical activation in all physiological simulations.</figcaption></figure>

transmembrane potential at different time instants. The activation wave originates in the anterior wall near Bachmann's bundle and propagates throughout the chamber. At $t = 50$ ms, the wavefront has advanced significantly across the anterior region, with accelerated propagation along Bachmann's bundle. By this time, the wave has also reached the atrial roof and starts its progression towards the middle section of the posterior part of the left atrium. At $t = 100$ ms, depolarization covers nearly the entire domain, with only the apex of the left atrial appendage remaining inactivated. Finally, repolarization of the atrial wall starts at around $t = 150$ ms near the Bachmann's bundle, and the tissue is fully repolarized by $t = 250$ ms. These activation patterns are in good agreement with physiological observations reported in the literature (see, e.g., [16, 71, 72]).
### 3D-shell model

The temporal evolution of left atrial pressure, volume and pressure-volume loop obtained with the 3D-shell model (Algorithm 1) are reported in Fig. 9. These results show reasonable agreement with experimental and numerical data reported in the literature (see, e.g., [46, 73, 74]). In particular, they accurately capture the characteristic three-phase function of the atrium. Initially, the booster pump phase is observed,

<figure><img src="image_10.png" /><figcaption>**Fig. 9:** Left atrial pressure, volume and pressure-volume loop obtained the 3D-shell model (Algorithm 1) under healthy conditions.</figcaption></figure>

characterized by a reduction in volume driven by active atrial contraction (a-wave). This is accompanied by a sustained increase in atrial pressure. As the atrium starts to repolarize, the pressure begins to decrease. Subsequently, the closure of the mitral valve marks the beginning of the reservoir phase. This phase is characterized by a temporary pressure spike (c-wave), followed by a pressure decrease as the atrium continues to relax. Then, both volume and pressure begin to rise due to inflow from the pulmonary veins (v-wave). As the atrium approaches its maximum capacity, the rate of filling slows down while the pressure increases more rapidly. Finally, the mitral valve reopens to start the conduit phase. This results in a sharp drop in both volume and pressure (y-wave). As the left ventricle fills, the emptying rate slows, eventually reaching a point where both variables remain nearly constant. The pressure-volume loop
(Fig. 9c) exhibits the characteristic eight-shape of the left atrium, with both the a-loop and v-loop displaying shapes consistent with the literature, though with a bigger a-loop than v-loop which is usually expected to be the opposite.

<figure><img src="image_11.png" /><figcaption>(a) Fluxes through pulmonary veins (PV) and (b) Pressure curves for left atrium (LA) and mitral valve (MV). (c) Pressure curves for left ventricle (LV).</figcaption></figure>

**Fig. 10:** Fluxes and pressures obtained with the 3D-shell model (Algorithm 1) under healthy conditions.

Fig. 10a shows the flux evolution across the pulmonary veins and the mitral valve. During atrial contraction (a-wave), the flow through the pulmonary veins is reduced. This is followed by the systolic forward flow (c-wave), characterized by a continuous decrease in pulmonary inflow immediately after the mitral valve closes. This flux reaches a minimum at the late systolic phase of the left ventricle (see Fig. 10b), marking the transition to diastolic forward flow (v-wave). This phase features a gradual increase in flow as the left atrial pressure converges with the pulmonary venous pressure, which produces a sustained decrease in the flux eventually reaching a plateau phase. Following this stabilization, the opening of the mitral valve (y-wave) marks a sudden decrease in the flux across the pulmonary veins followed by a sustained increase as the atrial pressure drops. Finally, the flux in the pulmonary veins reaches a plateau phase just before the next contraction begins. On the other hand, flow through the mitral valve becomes maximal during atrial contraction. This flow decreases as the left atrium begins to relax, eventually reaching zero when the mitral valve closes (at the beginning of the c-wave). The closure of the valve, induced by the pressure increase within the left ventricle, marks the start of the left atrial filling phase, characterized by zero flux through the mitral valve under physiological conditions. Finally, the mitral valve opens when the atrial pressure exceeds the ventricular pressure. This gradient drives a sustained increase in flow across the valve, followed by a plateau phase as the atrial and ventricular pressures converge (see Fig. 10a). The pressure evolution within the left ventricle (see Fig. 10b) exhibits physiological behavior, where the opening and closing of the mitral valve are dictated by the pressure gradient between the atrium and the ventricle, showing the hemodynamic coupling between these two chambers.
The left atrial pressure evolution remains within physiological ranges, reaching a minimum of 390.26 Pa at the end of the relaxation phase and a maximum of 1277.07 Pa at the end of the filling phase (see, e.g., [67, 73, 75, 76]). However, the volume variation is slightly lower than values reported in the literature (see, e.g., [77]). The atrial ejection fraction (EF) is defined as

$$EF \stackrel{\text{def}}{=} \frac{V_{\text{max}} - V_{\text{min}}}{V_{\text{max}}}$$

where $V_{\text{max}}$ and $V_{\text{min}}$ represent the maximum and minimum internal volumes of the left atrium, respectively. The computed ejection fraction is $EF = 32.08\%$, with minimum and maximum volumes of $V_{\text{min}} = 69.48$ mL and $V_{\text{max}} = 102.29$ mL. This value falls slightly below the physiological range of [45%–60%] typically reported in the experimental literature (see, e.g., [78–80]).

Another relevant parameter is the stroke volume (SV), which quantifies the total volume of blood ejected by the left atrium during a single cardiac cycle, and is formally defined as $SV \stackrel{\text{def}}{=} V_{\text{max}} - V_{\text{min}}$. Furthermore, this metric can be decomposed into passive and active contributions, corresponding to the atrial conduit (early diastolic) and contractile (booster pump) phases, respectively. The stroke volume for the atrial conduit phase is $SV_p = 12.48$ mL, while the contribution from the atrial pump phase is $SV_a = 20.41$ mL, yielding a total stroke volume of $SV = 32.89$ mL (see, e.g., [80]).

Finally, a key parameter to quantify the mechanical energy generated by the left atrium during the conduit and booster pump phases is the stroke work (SW). This metric corresponds to the area enclosed by the pressure-volume loop and is calculated as $SW \stackrel{\text{def}}{=} \int p_{\text{at}}(V)dV$. Similar to the stroke volume, the total stroke work is partitioned into the specific work performed during the active atrial pump phase and the passive conduit phase. The computed values are $SW_p = 16.06$ mmHg mL and $SW_a = 60.85$ mmHg mL, respectively, producing a total stroke work of $SW = 76.91$ mmHg mL.

### 3D-shell vs. full 3D formulation

The numerical results show that the spatial discretization has a significant impact on the predicted electromechanical response of the left atrium under physiological conditions. In particular, the single-layer tetrahedral discretization fails to reproduce a physically plausible evolution of the atrial volume, pressure, and flow (see Fig. 11).

This phenomenon is evidenced throughout the left atrial cardiac cycle, which is characterized by active contraction followed by relaxation and passive filling. While this sequence is accurately captured by the 3D-shell formulation, the single-layer tetrahedral mesh exhibits an artificially stiff response that delays atrial filling. This results in a non-physiological phase characterized by negligible volume changes—essentially a spurious plateau—which fails to represent the compliant nature of the atrial wall during diastole. Consequently, this lack of compliance prevents the model from accurately mimicking the pressure-volume relationship required for realistic hemodynamic simulations. On the other hand, introducing a second layer of tetrahedral elements substantially improves the mechanical response. The two-layer tetrahedral mesh recovers realistic pressure, volume, and flow evolutions that are in close agreement with the
<figure><img src="image_12.png" /><figcaption>**Fig. 11:** Comparison of pressures, volumes, and fluxes for the three discretizations (3D-shell, 3D two-layer tetrahedral, 3D one-layer tetrahedral) under healthy conditions.</figcaption></figure>

shell results. However, this improvement comes at the cost of increasing the number of elements from 39,131 prismatic elements to 234,786 tetrahedra, resulting in a 1.6-fold increase in execution time (see Table 1). Similar observations, can be inferred from Figs. 12 and 13, which reports some snapshots of the displacement in the deformed configuration obtained with the 3D-shell and the full 3D formulations.

Finally, Table 4 presents several relevant physiological biomarkers comparing the 3D-shell model with the one- and two-layer 3D tetrahedral models. As observed, the
<figure><img src="image_13.png" /><figcaption>**Fig. 12:** Snapshots of the atrial displacement in deformed configuration. Left: 3D-shell model; Middle: 3D (2-layers); Right: 3D (1-layer).</figcaption></figure>

3D-shell and the two-layer tetrahedral formulations yield lower minimum and higher maximum volumes of similar magnitude, contrasting with the restricted range of the single-layer case. This trend is also preserved in the pressure evolution, where the 3D-shell and two-layer models achieve lower minimum and higher maximum pressures. Notably, the 3D-shell model retrieves the highest stroke volume and ejection fraction, exceeding the single-layer 3D model by around 30.2%. Finally, the comparison among the stroke works reveals that the single-layer model provides a notably inferior stroke
<figure><img src="image_14.png" /><figcaption>**Fig. 13:** Snapshots of the atrial displacement in deformed configuration. Left: 3D-shell model; Middle: 3D (2-layers); Right: 3D (1-layer).</figcaption></figure>

work during the active pumping phase, representing only approximately 20% of the work performed by the 3D-shell and 3D two-layer models.

## 3.2 Pathological condition

This section presents the numerical results of left atrial electromechanics under the pathological condition of atrial fibrillation, a cardiac pathology that disrupts normal electrical propagation, leading to irregular and chaotic activity. Specifically, this condition alters the standard activation wavefronts, generating rotors or spiral waves throughout the left atrium (see, e.g., [81]).
<table><caption>**Table 4:** Main biomarkers for different spatial approximations (healthy condition).</caption><thead><tr><th>Variable</th><th>3D-Shell</th><th>3D (2-layer)</th><th>3D (1-layer)</th><th>Unit</th></tr></thead><tbody><tr><td>V<sub>min</sub></td><td>69.48</td><td>68.32</td><td>71.53</td><td>mL</td></tr><tr><td>V<sub>max</sub></td><td>102.30</td><td>100.75</td><td>96.62</td><td>mL</td></tr><tr><td>p<sub>min</sub></td><td>390.26</td><td>400.53</td><td>457.16</td><td>Pa</td></tr><tr><td>p<sub>max</sub></td><td>1277.07</td><td>1279.14</td><td>1167.03</td><td>Pa</td></tr><tr><td>SV<sub>p</sub></td><td>12.48</td><td>12.56</td><td>10.74</td><td>mL</td></tr><tr><td>SV<sub>a</sub></td><td>20.41</td><td>19.93</td><td>14.52</td><td>mL</td></tr><tr><td>SV</td><td>32.89</td><td>32.50</td><td>25.26</td><td>mL</td></tr><tr><td>EF</td><td>32.08</td><td>32.18</td><td>26.07</td><td>-</td></tr><tr><td>SW<sub>p</sub></td><td>16.06</td><td>16.23</td><td>24.67</td><td>mmHg mL</td></tr><tr><td>SW<sub>a</sub></td><td>60.85</td><td>57.97</td><td>14.22</td><td>mmHg mL</td></tr></tbody></table>

## Electrical activation

The disruption of electrical propagation results in a weak and disorganized wavefront that, from a biomechanical perspective, is insufficient to trigger a significant mechanical contraction, leading to the absence of the atrial kick. Accordingly, this state is modeled using (19) as for the physiological case, but with the left atrial active tension set to zero, viz., $\Sigma_a = 0$. This approach is consistent with numerous electromechanical studies on atrial fibrillation (see, e.g., [39, 40, 82–85]). Furthermore, the impact of atrial fibrillation extends to the pulmonary veins. Accordingly, the pulmonary venous pressure profile imposed in these pathological example lacks of a- and c-waves (see dotted line in Fig. 7a).

## 3D-shell model

Fig. 14 compares the evolution of pressure, volume, the pressure-volume loop, the left ventricular pressure, and the transmitral and veins fluxes obtained with the 3D-shell model (Algorithm 1) for the healthy and AF conditions. As observed, during atrial fibrillation (AF), the atrium loses its ability to contract effectively, leading to the absence of the atrial kick and a subsequent alteration of all hemodynamic quantities. Specifically, the pressure evolution in the AF case lacks the characteristic pressure increase associated with the atrial kick (a-wave).

Furthermore, following mitral valve closure, the pressure in the AF model exhibits a continuous rise, contrasting with the physiological case where a relaxation phase typically lowers the pressure before the reservoir phase begins. Conversely, upon mitral valve opening, the y-wave remains similar in both scenarios, indicating that AF primarily impacts the active contractile phase rather than passive filling. The evolution of the volume shows that during the atrial kick phase, the atrium under AF conditions fails to contracts, and thus does not reduce its volume. Consequently, immediately following mitral valve closure, the volume continues to increase during the v-wave, reaching a slightly larger maximum volume compared to the physiological case. When the mitral valve subsequently opens, the passive volume decrease (early filling) remains similar to the physiological baseline. These alterations are distinctly reflected in the
<figure><img src="image_15.png" /><figcaption>**Fig. 14:** Comparison between healthy (continuous lines) and atrial-fibrillation (dashed lines) conditions obtained with the 3d-shell model (Algorithm 1).</figcaption></figure>

pressure-volume loop, which notably lacks the active a-loop component during the pathological case.

Furthermore, the simulation of AF demonstrates a decrease in left ventricular preload, as indicated by the absence of the atrial kick. This reduction in filling subsequently lowers the maximum pressure developed within the left ventricle during the pathological cycle. Finally, the flux through the mitral valve shows that there is a decrease in the mitral flux during the atrial kick for the case of AF, which is expected since the atrium is solely passively filled during this pathology. Then, upon the opening
<figure><img src="image_16.png" /><figcaption>(a) Pressure-volume loop.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Left atrial pressure.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Pressure curve for the left ventricle.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) Volume curve for the left atrium.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) Flux across the mitral valve.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) Flux across the pulmonary veins.</figcaption></figure>

<caption>**Fig. 15:** Comparison of pressures, volumes, and fluxes for the three discretizations (3D-shell, 3D two-layer tetrahedral, 3D one-layer tetrahedral) under atrial-fibrillation.</caption>

of the mitral valve, the transmitral flux is nearly identical in both cases. In contrast, the flux through the pulmonary veins remains nearly constant during the late diastolic phase in the AF model. This stability results from the absence of the atrial kick, which impulse significant pressure variations within the pulmonary veins or the left atrium during this interval. Subsequently, upon mitral valve closure, the venous flux suffers a sudden drop as the left atrial pressure rises and approaches to the pressure in the pulmonary veins. Finally, once the mitral valve opens, the flow profile resembles the physiological case, initially slightly decreasing due to the predominance of the high transmitral flux, before increasing again as the pressure in the left atrium decreases.
### 3D-shell model vs. full 3D

Under atrial fibrillation, the absence of coordinated atrial contraction leads to a mechanical response dominated by low-amplitude, pressure-driven deformations. In this regime, all discretizations yield nearly identical pressure, volume, and flow evolutions, as shown in Fig. 15, up to minor shifts attributable to compliance differences between discretizations.

## 4 Discussion

The results reported in Section 3 show that the proposed 3D-shell formulation reproduces the essential electromechanical behavior of the left atrium at a substantially lower computational cost than standard 3D solid formulations. This is supported by the pressure-volume dynamics, deformation patterns, and flow profiles shown in Figs. 9 and 11, together with the execution times and memory usage reported in Table 1. The proposed 3D-shell model of atrial electromechanics integrates physiological electrical activation sequences, adratic through-thickness shell kinematics with MITC3 re-interpolation, transmural fiber architecture, atrium-specific active constitutive law, and the coupling with a dedicated 0D electromechanical model of the left ventricle.

As mentioned in Section 2.1.1, the thin-walled nature of the left atrium introduces a structural limitation in standard displacement-based 3D formulations such as (11)$_1$. When the wall thickness is small relative to the in-plane dimensions, low-order kinematic approximations lead to an artificially stiff response unless several elements are introduced across the thickness. This limitation is clearly observed in the single-layer tetrahedral discretization, which fails to reproduce a realistic evolution of atrial volume, pressure, and flow under physiological conditions (see Figs. 11 and 12). As shown in [48], this behavior is not a numerical locking artifact *per se*, but rather a modeling inconsistency induced by the fact that linear through-thickness kinematics are not asymptotically consistent with classical shell models as the thickness tends to zero. The artificial stiffening observed in the results of Section 3 therefore provides numerical evidence of this model inconsistency.

The 3D-shell formulation resolves this issue by adopting a quadratic behavior of the displacement field through the thickness. This kinematic assumption restores asymptotic consistency while retaining a fully 3D constitutive framework. In practice, this leads to a robust mechanical response, evidenced by larger displacement amplitudes, a realistic relaxation phase, and physiologically consistent pressure-volume loops (see Figs. 9 and 11), without requiring through-thickness mesh refinement.

Under physiological conditions, where active contraction induces large deformation gradients, the mechanical response is therefore critically dependent on the spatial discretization. Single-layer tetrahedral meshes systematically underestimate deformation and pressure amplitudes, leading to delayed relaxation, distorted pressure-volume dynamics, and altered flow patterns (Fig. 11). Introducing a second tetrahedral layer recovers pressure, volume, and flow evolutions comparable to those of the 3D-shell formulation, but at the cost of a nearly six-fold increase in element count and a 1.6-fold increase in execution time (Table 1). The 3D-shell formulation achieves comparable accuracy at a reduced computational cost.
In contrast, under atrial fibrillation, the absence of coordinated atrial contraction results in a mechanical response dominated by low-amplitude, pressure-driven deformations. In this regime, all discretizations yield nearly identical pressure, volume, and flow evolutions (Fig. 15), up to minor shifts attributable to compliance differences. This behavior follows directly from the kinematic regime involved: when deformation amplitudes remain small, the asymptotic inconsistency of low-order through-thickness kinematics identified in [48] is not activated. These results are in agreement with clinical observations and with previously reported numerical studies of AF (see, e.g., the previous work [51]).

Accurate estimation of left ventricular preload relies on a compliant mechanical response of the left atrium [86–89]. Inaccurate atrial pressure predictions therefore propagate directly to the ventricular response. The single-layer tetrahedral discretization systematically underestimates preload relative to both the two-layer and 3D-shell models (see Fig. 11c), resulting in reduced ventricular filling. Consistently, the hemodynamic analysis shows that the 3D-shell and two-layer discretizations produce nearly identical mitral and pulmonary venous flow profiles, whereas the single-layer mesh underestimates mitral inflow and delays the pulmonary venous v-wave peak, reflecting an artificial attenuation of atrial compliance.

While the proposed 3D-shell formulation incorporates the MITC3 technique to prevent shear locking, numerical evidence (not reported here) indicates that similar mechanical responses are obtained with and without the MITC3 reinterpolation. This suggests that shear locking is not a dominant numerical issue in the present application. Indeed, since the atrial tissue undergoes significant inflation and stretching, the deformation is likely membrane-dominated, such that shear locking is not a fundamental here. Nevertheless, owing to the complex and non-uniform kinematics of the left atrium under physiological and pathological conditions, the occurrence of other locking modes cannot be entirely excluded. A truly robust model of the left atrium would therefore benefit from advanced finite element formulations capable of accounting for all potential locking phenomena (see, e.g., [48, 49, 90]).

From a clinical perspective, 3D-shell modeling of the left atrium constitutes a meaningful step toward numerically reliable patient-specific simulations, as required by emerging digital-twin approaches (see, e.g., [91–93]). Applications such as atrial fibrillation analysis or catheter ablation planning require mechanically accurate predictions at a computational cost compatible with routine use, making spatial discretization a decisive modeling choice rather than a secondary numerical detail. The present results show that single-layer tetrahedral discretizations, although inexpensive, fail to deliver physiologically consistent atrial mechanics and systematically underestimate ventricular filling, stroke volume, and cardiac output. While two-layer tetrahedral meshes recover accuracy, their computational and memory costs limit their practical applicability [94, 95]. By contrast, the 3D-shell formulation provides a compliant and mechanically consistent atrial response without through-thickness refinement, offering a more favorable balance between accuracy and efficiency for clinically oriented electromechanical simulations.
## 5 Limitations

Despite the advantages of the proposed 3D-shell electromechanical model of the left atrium, several limitations of the present work warrant discussion. First, the model does not include a fully coupled 3D-0D closed-loop model of circulation. In the current formulation, venous inflow is prescribed through an imposed pressure profile (Fig. 7a), which lacks the dynamic feedback mechanisms inherent to a closed-loop system. Consequently, venous pressure does not adapt to changes in atrial or ventricular hemodynamics. A natural extension of this work would be the incorporation of a complete closed-loop circulation model, similar to the approach reported in [44]. Nevertheless, such an extension lies beyond the primary scope of this study, whose central objective is to investigate the benefits of the 3D-shell modeling of the left atrium. Second, mechanical contact with the pericardial sac is neglected. The pericardium provides a physical constraint that can significantly limit atrial motion, particularly during rapid filling phases or under pathological conditions such as constrictive pericarditis (see, e.g., [30, 36, 96]). Third, the left atrium and left ventricle exhibit strong mechanical interdependence. From this perspective, coupling the present atrial 3D-shell model with a 3D electromechanical model of the left ventricle would be desirable to improve both the predicted atrial and ventricular function. Such a multi-chamber framework would also enable the representation of longitudinal mitral annular motion, which plays a key role during atrial filling. A major strength of the proposed formulation lies in its suitability as a calibration platform for more advanced modeling frameworks, including fluid-structure interaction (FSI) models of the left atrium. The incorporation of 3D fluid dynamics is essential for a comprehensive description of both chamber hemodynamics and tissue mechanics. While several studies have investigated atrial hemodynamics and their influence on left ventricular flow patterns (see [97–99]), only a limited number have addressed the fully coupled FSI problem, owing to its substantial complexity and computational cost (see [52, 100–102]). In this context, the proposed 3D-shell electromechanical model represents a natural step toward a robust and computationally efficient FSI framework (see [103]).

## 6 Conclusion

This work introduces a 3D-shell electromechanical model of the left atrium that combines physiological fidelity with substantially reduced computational cost. By adopting quadratic through-thickness kinematics within a fully 3D constitutive framework, the proposed approach overcomes the artificial stiffening issues commonly associated with standard low-order 3D solid discretizations in thin-walled atrial geometries.

Numerical results under both physiological conditions and atrial fibrillation demonstrate that the 3D-shell model reproduces atrial deformation and pressure-volume loops with accuracy comparable to that of multi-layer 3D tetrahedral meshes, while requiring significantly lower computational cost. In particular, the numerical comparisons revealed that a single-layer tetrahedral mesh captures only about 20% of the stroke work performed during the active pumping phase, leading to a limited atrial
response and an underestimation of ventricular filling. Moreover, the 3D-shell formulation captures these compliant atrial mechanics without resorting to through-thickness mesh refinement, which is mandatory for standard 3D formulations.

These findings highlight spatial discretization as a critical numerical choice in atrial electromechanics and support the use of shell-based formulations for efficient and mechanically consistent simulations. The proposed 3D-shell model therefore provides a solid foundation for patient-specific applications and for future extensions toward multi-chamber and fluid-structure interaction frameworks.

## Declarations

**Competing interests.** The authors declare that they have no competing interests.

**Data availability.** No datasets were generated or analyzed during the current study.

**Funding.** This work was partially supported by Bpifrance through the MEDITWIN project sustained by France 2030 program of the French Government. OR was partially supported by the Vicerrectoria de Investigación, Desarrollo e Innovación (VRIDEI) of the Universidad de Santiago de Chile, and the National Agency for Research and Development (ANID) Doctorado Becas Chile/2020 - 72210101 of the Government of Chile.

**Authors' contributions.** Conceptualization: CB, DC, MF, MV, OR. Formal analysis and investigation: CB, DC, MF, MV, OR. Funding acquisition: MF, OR. Methodology: CB, DC, MF, MV, OR. Software: CB, MF, MV, OR. Supervision: DC, MF, MV. Writing - original draft preparation: CB, OR. Writing - review and editing: CB, DC, MF, MV, OR.

**Acknowledgments.** The authors wish to thank Philippe Moireau (École Polytechnique) for the valuable discussions on cardiac modeling.

## A Electrophysiological model

In this work, the electrical activity of the left atrium is described using the bidomain model (see, e.g., [104, 105]). Cardiac electrophysiology is hence formulated as a coupled PDE-ODE system involving the cellular ionic dynamics, a nonlinear reaction-diffusion equation for the transmembrane potential and an elliptic equation for the extracellular potential. The electrophysiological problem of the left atrium reads: Find the trans-membrane potential $V_m: \Omega \times \mathbb{R}^+ \to \mathbb{R}$, the extracellular potential $\phi_e: \Omega \times \mathbb{R}^+ \to \mathbb{R}$ and the ionic variable $c: \Omega \times \mathbb{R}^+ \to \mathbb{R}^{\text{ion}}$ such that

$$ 
\left\{
\begin{array}{ll}
A_m(C_m \partial_t V_m + I_{\text{ion}}(V_m, c)) - \text{div}(\boldsymbol{\sigma}_i \nabla V_m) - \text{div}(\boldsymbol{\sigma}_i \nabla \phi_e) = A_m I_{\text{app}} & \text{in } \Omega, \\
\text{div}((\boldsymbol{\sigma}_i + \boldsymbol{\sigma}_e) \nabla \phi_e) + \text{div}(\boldsymbol{\sigma}_i \nabla V_m) = 0 & \text{in } \Omega, \\
\partial_t c + g(V_m, c) = 0 & \text{in } \Omega, \\
\boldsymbol{\sigma}_i \nabla V_m \cdot \mathbf{n} + \boldsymbol{\sigma}_i \nabla \phi_e \cdot \mathbf{n} = 0 & \text{on } \partial\Omega, \\
(\boldsymbol{\sigma}_i + \boldsymbol{\sigma}_e) \nabla \phi_e \cdot \mathbf{n} + \boldsymbol{\sigma}_i \nabla V_m \cdot \mathbf{n} = 0 & \text{on } \partial\Omega.
\end{array}
\right.
\tag{26}
$$
Here, $A_m$ denotes the membrane surface/volume ratio, $C_m$ is the membrane capacitance per unit area, $I_{ion}$ represents the nonlinear ionic current describing the exchange of ions across the cell membrane, and $I_{app}$ denotes a prescribed external stimulus current.

The ionic current $I_{ion}$ is modeled based on the physiological ionic model proposed in [106], a widely adopted framework for simulating atrial electrical activity (see, e.g., [40, 72, 107–109]). The symbols $\sigma_i$ and $\sigma_e$ denote the intra- and extracellular conductivity tensors, respectively. The electrical conductivity properties of cardiac tissue are inherently anisotropic, exhibiting higher values along the fiber orientation compared to the transverse plane. This anisotropic behavior is incorporated into the conductivity tensors through the following relation:

$$ \sigma_i \stackrel{\text{def}}{=} \sigma_i^t \mathbb{I} + (\sigma_i^l - \sigma_i^t) \mathbf{f} \otimes \mathbf{f}, \quad \sigma_e \stackrel{\text{def}}{=} \sigma_e^t \mathbb{I} + (\sigma_e^l - \sigma_e^t) \mathbf{f} \otimes \mathbf{f}, $$

where $\mathbf{f}$, as in Section 2.1, stands for the unit vector aligned with the fiber direction, $\mathbb{I}$ is the identity matrix and $\sigma_k^l$ and $\sigma_k^t$, with $k \in \{i, e\}$, are the conductivity coefficients along and transverse to the fiber, respectively. The Neumann boundary conditions

<table><caption>**Table 5:** Conductivity values of the electrophysiological simulations. BB: Banchmann bundle. FO: Fossa ovalis.</caption><thead><tr><th>Parameter</th><th>Conductivity</th><th>Tissue</th><th>BB</th><th>FO</th><th>Unit</th></tr></thead><tbody><tr><td>$\sigma_i^t$</td><td>Interior transverse direction</td><td>$2.2 \cdot 10^{-4}$</td><td>$1.045 \cdot 10^{-3}$</td><td>$2.2 \cdot 10^{-4}$</td><td>S cm<sup>-1</sup></td></tr><tr><td>$\sigma_i^l$</td><td>Interiopr longitudinal direction</td><td>$2.2 \cdot 10^{-3}$</td><td>$1.639 \cdot 10^{-3}$</td><td>$2.2 \cdot 10^{-3}$</td><td>S cm<sup>-1</sup></td></tr><tr><td>$\sigma_e^t$</td><td>Exterior transverse direction</td><td>$9.0 \cdot 10^{-4}$</td><td>$4.275 \cdot 10^{-3}$</td><td>$9.0 \cdot 10^{-4}$</td><td>S cm<sup>-1</sup></td></tr><tr><td>$\sigma_e^l$</td><td>Exterior longitudinal direction</td><td>$2.2 \cdot 10^{-3}$</td><td>$1.435 \cdot 10^{-3}$</td><td>$2.07 \cdot 10^{-4}$</td><td>S cm<sup>-1</sup></td></tr></tbody></table>

(26)$_{4,5}$ enforce that the left atrium is electrically isolated, that is, the intracellular and extracellular currents vanish on its boundary $\Omega$. As the bidomain formulation (26) defines $\phi_e$ only up to an additive constant, a unique solution is obtained by imposing, for instance, the zero-mean constraint $\int_{\Omega} \phi_e = 0$.

## B Atrial fiber generation

The fiber architecture of the left atrium is modeled using a surface-based adaptation of the Laplace-Dirichlet Rule-Based (LDRB) framework proposed in [16]. To capture the significant transmural heterogeneity observed between the inner and outer layers [68, 110], the algorithm is applied separately to the endocardial $\Gamma_{en}$ and epicardial $\Gamma_{ep}$ surfaces, rather than solving for a single volumetric field.

First, each surface is partitioned into anatomically relevant boundaries, including the pulmonary vein rings, the mitral annulus, and the appendage apex. These boundaries serve as Dirichlet constraints for a set of harmonic potentials $\Psi_j$ (where $j \in \{\text{en, ep}\}$), governed by the surface Laplace-Beltrami operator. The gradients of
these scalar fields define the local fiber direction for specific anatomical bundles. Distinct regions, such as the circular fibers surrounding the pulmonary veins or the longitudinal structures of the septal wall, are identified through threshold-based selection rules adapted from [16]. Notably, the epicardial definition includes a specific correction for the Bachmann's Bundle to ensure fibers align parallel to the mitral valve ring, contrasting with the orientation of the underlying endocardium.

Once the fiber orientations are determined on both surfaces, the full 3D fiber architecture is constructed via transmural interpolation. Following the approach in [69], the fiber vector $f$ at any point across the wall thickness is defined as

$$f(\xi^1, \xi^2, \xi^3) \stackrel{\text{def}}{=} f_0(\xi^1, \xi^2) \cos \left( \frac{2\theta(\xi^1, \xi^2)\xi^3}{\epsilon(\xi^1, \xi^2)} \right) + f_0^{\perp}(\xi^1, \xi^2) \sin \left( \frac{2\theta(\xi^1, \xi^2)\xi^3}{\epsilon(\xi^1, \xi^2)} \right).$$

where $f_0$ denotes the fiber orientation at the mid-surface of the left atrium, which lies in the tangential plane, and $f_0^{\perp}$ is the fiber orthogonal to $f_0$, also within the tangential plane. Additionally, $\epsilon$ represents the thickness of the left atrium at a given point, and $\theta$ is half the angle between the fibers in the epicardium and endocardium, as illustrated in Fig. 16. The resulting fiber distribution is illustrated in Fig. 5b.

<figure><img src="image_22.png" /><figcaption>Fig. 16: Angle variation across the thickness of the left atrium.</figcaption></figure>

## References

[1] Trayanova, N.A., Winslow, R.: Whole-heart modeling. Circulation Research 108(1), 113–128 (2011)

[2] Chabiniok, R., Wang, V.Y., Hadjicharalambous, M., Asner, L., Lee, J., Sermesant, M., Kuhl, E., Young, A.A., Moireau, P., Nash, M.P., Chapelle, D., Nordsletten, D.A.: Multiphysics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics. Interface Focus 6(2), 20150083 (2016)

[3] Niederer, S.A., Lumens, J., Trayanova, N.A.: Computational models in cardiology. Nature reviews cardiology 16(2), 100–111 (2019)
[4] Corral-Acero, J., Margara, F., Marciniak, M., Rodero, C., Loncaric, F., Feng, Y., Gilbert, A., Fernandes, J.F., Bukhari, H.A., Wajdan, A., Martinez, M.V., Santos, M.S., Shamohammdi, M., Luo, H., Westphal, P., Leeson, P., DiAchille, P., Gurev, V., Mayr, M., Geris, L., Pathmanathan, P., Morrison, T., Cornelussen, R., Prinzen, F., Delhaas, T., Doltra, A., Sitges, M., Vigmond, E.J., Zacur, E., Grau, V., Rodriguez, B., Remme, E.W., Niederer, S., Mortier, P., McLeod, K., Potse, M., Pueyo, E., Bueno-Orovio, A., Lamata, P.: The ‘digital twin’ to enable the vision of precision cardiology. European Heart Journal 41(48), 4556–4564 (2020)

[5] Peirlinck, M., Costabal, F.S., Yao, J., Guccione, J., Tripathy, S., Wang, Y., Ozturk, D., Segars, P., Morrison, T., Levine, S., et al.: Precision medicine in human heart modeling: Perspectives, challenges, and opportunities. Biomechanics and modeling in mechanobiology 20, 803–831 (2021)

[6] Bishop, M., Rajani, R., Plank, G., Gaddum, N., Carr-White, G., Wright, M., O'Neill, M., Niederer, S.: Three-dimensional atrial wall thickness maps to inform catheter ablation procedures for atrial fibrillation. EP Europace 18(3), 376–383 (2015)

[7] Whitaker, J., Rajani, R., Chubb, H., Gabrawi, M., Varela, M., Wright, M., Niederer, S., O'Neill, M.D.: The role of myocardial wall thickness in atrial arrhythmogenesis. EP Europace 18(12), 1758–1772 (2016)

[8] Hoermann, J.M., Pfaller, M.R., Avena, L., Bertoglio, C., Wall, W.A.: Automatic mapping of atrial fiber orientations for patient-specific modeling of cardiac electromechanics using image registration. International Journal for Numerical Methods in Biomedical Engineering 35(6), 3190 (2019)

[9] Seemann, G., Höper, C., Sachse, F.B., Dössel, O., Holden, A.V., Zhang, H.: Heterogeneous three-dimensional anatomical and electrophysiological model of human atria. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 364(1843), 1465–1481 (2006)

[10] Plank, G., Prassl, A., Hofer, E., Trayanova, N.A.: Evaluating intramural virtual electrodes in the myocardial wedge preparation: Simulations of experimental conditions. Biophysical Journal 94(5), 1904–1915 (2008)

[11] Krueger, M.W., Schmidt, V., Tobón, C., Weber, F.M., Lorenz, C., Keller, D.U.J., Barschdorf, H., Burdumy, M., Neher, P., Plank, G., Rhode, K., Seemann, G., Sanchez-Quintana, D., Saiz, J., Razavi, R., Dössel, O.: Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach. In: Metaxas, D.N., Axel, L. (eds.) Functional Imaging and Modeling of the Heart, pp. 223–232. Springer, Berlin, Heidelberg (2011)

[12] Labarthe, S., Coudiere, Y., Henry, J., Cochet, H.: A semi-automatic method to construct atrial fibre structures: A tool for atrial simulations. In: 2012
[13] Tobón, C., Ruiz-Villa, C.A., Heidenreich, E., Romero, L., Hornero, F., Saiz, J.: A three-dimensional human atrial model with fiber orientation. electrograms and arrhythmic activation patterns relationship. PLOS ONE **8**(2), 1–13 (2013)

[14] Ferrer, A., Sebastián, R., Sánchez-Quintana, D., Rodríguez, J.F., Godoy, E.J., Martínez, L., Saiz, J.: Detailed anatomical and electrophysiological models of human atria and torso for the simulation of atrial activation. PLOS ONE **10**(11), 1–29 (2015)

[15] Rocher, S., López, A., Ferrer, A., Martínez, L., Sánchez, D., Saiz, J.: A highly-detailed 3d model of the human atria. In: Lhotska, L., Sukupova, L., Lacković, I., Ibbott, G.S. (eds.) World Congress on Medical Physics and Biomedical Engineering 2018, pp. 649–653. Springer, Singapore (2019)

[16] Piersanti, R., Africa, P.C., Fedele, M., Vergara, C., Dedè, L., Corno, A.F., Quarteroni, A.: Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations. Computer Methods in Applied Mechanics and Engineering **373**, 113468 (2021)

[17] Usyk, T.P., LeGrice, I.J., McCulloch, A.D.: Computational model of three-dimensional cardiac electromechanics. Computing and Visualization in Science **4**(4), 249–257 (2002)

[18] Smith, N.P., Nickerson, D.P., Crampin, E.J., Hunter, P.J.: Multiscale computational modelling of the heart. Acta Numerica **13**, 371–431 (2004)

[19] Crampin, E.J., Halstead, M., Hunter, P., Nielsen, P., Noble, D., Smith, N., Tawhai, M.: Computational physiology and the physiome project. Experimental Physiology **89**(1), 1–26 (2004)

[20] Chapelle, D., Fernández, M.A., Gerbeau, J.-F., Moireau, P., Sainte-Marie, J., Zemzemi, N.: Numerical simulation of the electromechanical activity of the heart. In: International Conference on Functional Imaging and Modeling of the Heart, pp. 357–365 (2009). Springer

[21] Göktepe, S., Kuhl, E.: Electromechanics of the heart: a unified approach to the strongly coupled excitation–contraction problem. Computational Mechanics **45**(2), 227–243 (2010)

[22] Nordsletten, D.A., Niederer, S.A., Nash, M.P., Hunter, P.J., Smith, N.P.: Coupling multi-physics models to cardiac mechanics. Progress in Biophysics and Molecular Biology **104**(1), 77–88 (2011)

[23] Chapelle, D., Le Tallec, P., Moireau, P., Sorine, M.: Energy-preserving muscle tissue model: formulation and compatible discretizations. International Journal
for Multiscale Computational Engineering 10(2) (2012)

[24] Sermesant, M., Chabiniok, R., Chinchapatnam, P., Mansi, T., Billet, F., Moireau, P., Peyrat, J.M., Wong, K., Relan, J., Rhode, K., Ginks, M., Lambiase, P., Delingette, H., Sorine, M., Rinaldi, C.A., Chapelle, D., Razavi, R., Ayache, N.: Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation. Medical Image Analysis 16(1), 201–215 (2012)

[25] Hirschvogel, M., Bassilious, M., Jagschies, L., Wildhirt, S.M., Gee, M.W.: A monolithic 3d-0d coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics. International Journal for Numerical Methods in Biomedical Engineering 33(8), 2842 (2017)

[26] Garcia-Blanco, E., Ortigosa, R., Gil, A.J., Bonet, J.: Towards an efficient computational strategy for electro-activation in cardiac mechanics. Computer Methods in Applied Mechanics and Engineering 356, 220–260 (2019)

[27] Augustin, C.M., Gsell, M.A.F., Karabelas, E., Willemen, E., Prinzen, F.W., Lumens, J., Vigmond, E.J., Plank, G.: A computationally efficient physiologically comprehensive 3d-0d closed-loop model of the heart and circulation. Computer Methods in Applied Mechanics and Engineering 386, 114092 (2021)

[28] Piersanti, R., Regazzoni, R., Salvador, S., Corno, A.F., Dede', L., Vergara, C., Quarteroni, A.: 3D-0D closed-loop model for the simulation of cardiac biventricular electromechanics. Computer Methods in Applied Mechanics and Engineering 391, 114607 (2022)

[29] Di Martino, E.S., Bellini, C., Schwartzman, D.S.: In vivo porcine left atrial wall stress: Computational model. Journal of Biomechanics 44(15), 2589–2594 (2011)

[30] Fritz, T., Wieners, C., Seemann, G., Steen, H., Dössel, O.: Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology 13(3), 627–641 (2014)

[31] Augustin, C.M., Neic, A., Liebmann, M., Prassl, A.A., Niederer, S.A., Haase, G., Plank, G.: Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation. Journal of Computational Physics 305, 622–646 (2016)

[32] Strocchi, M., Gsell, M.A.F., Augustin, C.M., Razeghi, O., Roney, C.H., Prassl, A.J., Vigmond, E.J., Behar, J.M., Gould, J.S., Rinaldi, C.A., Bishop, M.J., Plank, G., Niederer, S.A.: Simulating ventricular systolic motion in a four-chamber heart model with spatially varying robin boundary conditions to model
the effect of the pericardium. Journal of Biomechanics **101**, 109645 (2020)

[33] Satriano, A., Vigmond, E.J., Schwartzman, D.S., Di Martino, E.S.: Mechanoelectric finite element model of the left atrium. Computers in Biology and Medicine **96**, 24–31 (2018)

[34] Feng, L., Gao, H., Griffith, B., Niederer, S., Luo, X.: Analysis of a coupled fluid-structure interaction model of the left atrium and mitral valve. International Journal for Numerical Methods in Biomedical Engineering **35**(11), 3254 (2019)

[35] Land, S., Niederer, S.A.: Influence of atrial contraction dynamics on cardiac function. International journal for numerical methods in biomedical engineering **34**(3), 2931 (2018)

[36] Pfaller, M.R., Hörmann, J.M., Weigl, M., Nagler, A., Chabiniok, R., Bertoglio, C., Wall, W.A.: The importance of the pericardium for cardiac biomechanics: from physiology to computational modeling. Biomechanics and modeling in mechanobiology **18**, 503–529 (2019)

[37] Augustin, C.M., Fastl, T.E., Neic, A., Bellini, C., Whitaker, J., Rajani, R., O'Neill, M.D., Bishop, M.J., Plank, G., Niederer, S.A.: The impact of wall thickness and curvature on wall stress in patient-specific electromechanical models of the left atrium. Biomechanics and modeling in mechanobiology **19**(3), 1015–1034 (2020)

[38] Baillargeon, B., Rebelo, N., Fox, D.D., Taylor, R.L., Kuhl, E.: The living heart project: A robust and integrative simulator for human heart function. European Journal of Mechanics - A/Solids **48**, 38–47 (2014). Frontiers in Finite-Deformation Electromechanics

[39] Gerach, T., Schuler, S., Fröhlich, J., Lindner, L., Kovacheva, E., Moss, R., Wülfers, E.M., Seemann, G., Wieners, C., Loewe, A.: Electro-mechanical whole-heart digital twins: a fully coupled multi-physics approach. Mathematics **9**(11), 1247 (2021)

[40] Fedele, M., Piersanti, R., Regazzoni, F., Salvador, M., Africa, P.C., Bucelli, M., Zingaro, A., Quarteroni, A., et al.: A comprehensive and biophysically detailed computational model of the whole human heart electromechanics. Computer Methods in Applied Mechanics and Engineering **410**, 115983 (2023)

[41] Feng, L., Gao, H., Luo, X.: Whole-heart modelling with valves in a fluid-structure interaction framework. Computer Methods in Applied Mechanics and Engineering **420**, 116724 (2024)

[42] Gonzalo, A., Augustin, C.M., Bifulco, S.F., Telle, A., Chahine, Y., Kassar, A., Guerrero-Hurtado, M., Durán, E., Martínez-Legazpi, P., Flores, O., Bermejo, J., Plank, G., Akoum, N., Boyle, P.M., Alamo, J.C.: Multiphysics simulations
reveal haemodynamic impacts of patient-derived fibrosis-related changes in left atrial tissue mechanics. *The Journal of Physiology* **602**(24), 6789–6812 (2024)

[43] Viola, F., Del Corso, G., Verzicco, R.: High-fidelity model of the human heart: An immersed boundary implementation. *Phys. Rev. Fluids* **8**, 100502 (2023)

[44] Telle, A., Kassar, A., Chamoun, N., Haykal, R., Gonzalo, A., Hensley, T., Chahine, Y., Flores, O., Alamo, J.C., Akoum, N., Augustin, C.M., Boyle, P.M.: Systematic computational assessment of atrial function impairment due to fibrotic remodeling in electromechanical properties. *PLOS Computational Biology* **21**(12), 1–33 (2025)

[45] Arjoune, T., Bilas, C., Meierhofer, C., Stern, H., Ewert, P., Gee, M.W.: Inverse analysis of patient-specific parameters of a 3d-0d closed-loop cardiovascular model with an exemplary application to an adult tetralogy of fallot case. *Biomechanics and Modeling in Mechanobiology* **24**(6), 2039–2068 (2025)

[46] Shi, L., Gan, B., Chen, I.Y., Vedula, V.: Personalized multiscale modeling of left atrial mechanics and blood flow. *Computer Methods in Applied Mechanics and Engineering* **448**, 118412 (2026)

[47] Chapelle, D., Bathe, K.-J.: The Finite Element Analysis of Shells - Fundamentals, 2nd edn. Computational Fluid and Solid Mechanics, p. 410. Springer, Berlin, Heidelberg (2011). https://doi.org/10.1007/978-3-642-16408-8

[48] Chapelle, D., Ferent, A., Bathe, K.J.: 3D-shell elements and their underlying mathematical model. *Mathematical Models and Methods in Applied Sciences* **14**(01), 105–142 (2004)

[49] Bathe, K.-J., Iosilevich, A., Chapelle, D.: An evaluation of the MITC shell elements. *Computers & Structures* **75**(1), 1–30 (2000)

[50] Caruel, M., Chabiniok, R., Moireau, P., Lecarpentier, Y., Chapelle, D.: Dimensional reductions of a cardiac model for effective validation and calibration. *Biomechanics and modeling in mechanobiology* **13**(4), 897–914 (2014)

[51] Ruz, O., Vidrascu, M., Chapelle, D., Fernández, M.A.: 3d-shell electromechanical modeling of the left atrium. In: International Conference on Functional Imaging and Modeling of the Heart, pp. 48–59 (2025). Springer

[52] Feng, L., Gao, H., Qi, N., Danton, M., Hill, N.A., Luo, X.: Fluid-structure interaction in a fully coupled three-dimensional mitral-atrium-pulmonary model. *Biomechanics and modeling in mechanobiology* **20**(4), 1267–1295 (2021)

[53] Demiray, H.: A note on the elasticity of soft biological tissues. *Journal of biomechanics* **5**(3), 309–311 (1972)
[54] Gasser, T.C., Ogden, R.W., Holzapfel, G.A.: Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of the royal society interface 3(6), 15–35 (2006)

[55] Kimmig, F., Chapelle, D., Moireau, P.: Thermodynamic properties of muscle contraction models and associated discrete-time principles. Advanced Modeling and Simulation in Engineering Sciences 6(1), 1–36 (2019)

[56] Bestel, J., Clément, F., Sorine, M.: A biomechanical model of muscle contraction. In: International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 1159–1161 (2001). Springer

[57] Holzapfel, G.A., Ogden, R.W.: Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367(1902), 3445–3475 (2009)

[58] Hart, J.: Normal resting pulse rate ranges. Journal of Nursing Education and Practice 5(8), 95–98 (2015)

[59] Sussman, T., Bathe, K.-J.: 3d-shell elements for structures in large strains. Computers & structures 122, 2–12 (2013)

[60] Chapelle, D., Paris, I.: Detailed reliability assessment of triangular MITC elements for thin shells. Computers & Structures 86(23), 2192–2202 (2008)

[61] Jeon, H.-M., Lee, P.-S., Bathe, K.-J.: The mitc3 shell finite element enriched by interpolation covers. Computers & Structures 134, 128–142 (2014)

[62] Jeon, H.-M., Lee, Y., Lee, P.-S., Bathe, K.-J.: The mitc3+ shell element in geometric nonlinear analysis. Computers & Structures 146, 91–104 (2015)

[63] Sainte-Marie, J., Chapelle, D., Cimrman, R., Sorine, M.: Modeling and estimation of the cardiac electromechanical activity. Computers & Structures 84(28), 1743–1759 (2006)

[64] Besl, P.J., McKay, N.D.: Method for registration of 3-d shapes. In: Sensor Fusion IV: Control Paradigms and Data Structures, vol. 1611, pp. 586–606 (1992). Spie

[65] Gee, M.W., Förster, C., Wall, W.: A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. International Journal for Numerical Methods in Biomedical Engineering 26(1), 52–72 (2010)

[66] Moireau, P., Xiao, N., Astorino, M., Figueroa, C.A., Chapelle, D., Taylor, C.A., Gerbeau, J.-F.: External tissue support and fluid-structure simulation in blood flows. Biomechanics and modeling in mechanobiology 11, 1–18 (2012)

[67] Ho, S.Y., Cabrera, J.A., Sanchez-Quintana, D.: Left atrial anatomy revisited.
Circulation: Arrhythmia and Electrophysiology 5(1), 220-228 (2012)

[68] Pashakhanloo, F., Herzka, D.A., Ashikaga, H., Mori, S., Gai, N., Bluemke, D.A., Trayanova, N.A., McVeigh, E.R.: Myofiber architecture of the human atria as revealed by submillimeter diffusion tensor imaging. Circulation: arrhythmia and electrophysiology 9(4), 004133 (2016)

[69] Chapelle, D., Collin, A., Gerbeau, J.-F.: A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling. Math. Models Methods Appl. Sci. 23(14), 2749-2776 (2013)

[70] Porumbescu, S.D., Budge, B., Feng, L., Joy, K.I.: Shell maps. ACM Transactions on Graphics (TOG) 24(3), 626-633 (2005)

[71] Zhao, J., Stephenson, R.S., Sands, G.B., LeGrice, I.J., Zhang, H., Jarvis, J.C., Smaill, B.H.: Atrial fibrosis and atrial fibrillation: a computer simulation in the posterior left atrium. In: Functional Imaging and Modeling of the Heart: 7th International Conference, FIMH 2013, London, UK, June 20-22, 2013. Proceedings 7, pp. 400-408 (2013). Springer

[72] Schenone, E., Collin, A., Gerbeau, J.-F.: Numerical simulation of electrocardiograms for full cardiac cycles in healthy and pathological conditions. International journal for numerical methods in biomedical engineering 32(5), 02744 (2016)

[73] Pagel, P.S., Kehl, F., Gare, M., Hettrick, D.A., Kersten, J.R., Warltier, D.C.: Mechanical function of the left atrium: new insights based on analysis of pressure-volume relations and doppler echocardiography. Anesthesiology 98(4), 975-994 (2003)

[74] Stöcker, Y., Guerrero-Hurtado, M., Durán, E., Gonzalo, A., Ristić, Z., Telle, Å., Kassar, A., Haykal, R., Akoum, N., Boyle, P.M., et al.: Effect of wall motion sampling on cfd-derived left atrial flow metrics. bioRxiv, 2025-10 (2025)

[75] Kistler, P.M., Sanders, P., Fynn, S.P., Stevenson, I.H., Spence, S.J., Vohra, J.K., Sparks, P.B., Kalman, J.M.: Electrophysiologic and electroanatomic changes in the human atrium associated with age. Journal of the American College of Cardiology 44(1), 109-116 (2004)

[76] Roşca, M., Lancellotti, P., Popescu, B.A., Piérard, L.A.: Left atrial function: pathophysiology, echocardiographic assessment, and clinical applications. Heart 97(23), 1982-1989 (2011)

[77] Abhayaratna, W.P., Seward, J.B., Appleton, C.P., Douglas, P.S., Oh, J.K., Tajik, A.J., Tsang, T.S.: Left atrial size: physiologic determinants and clinical applications. Journal of the American College of Cardiology 47(12), 2357-2363 (2006)
[78] Iaizzo, P.A. (ed.): Handbook of Cardiac Anatomy, Physiology, and Devices, 4th edn., p. 1009. Springer, Cham, Switzerland (2024). https://doi.org/10.1007/978-3-031-72581-4

[79] Sánchez-Quintana, D., Doblado-Calatrava, M., Cabrera, J.A., Macías, Y., Saremi, F.: Anatomical basis for the cardiac interventional electrophysiologist. BioMed Research International **2015**(1), 547364 (2015)

[80] Sugimoto, T., Robinet, S., Dulgheru, R., Bernard, A., Ilardi, F., Contu, L., Addetia, K., Caballero, L., Kacharava, G., Athanassopoulos, G.D., et al.: Echocardiographic reference ranges for normal left atrial function parameters: results from the eacvi norre study. European Heart Journal-Cardiovascular Imaging **19**(6), 630–638 (2018)

[81] Jalife, J.: Rotors and spiral waves in atrial fibrillation. Journal of cardiovascular electrophysiology **14**(7), 776–780 (2003)

[82] Adeniran, I., MacIver, D.H., Garratt, C.J., Ye, J., Hancox, J.C., Zhang, H.: Effects of persistent atrial fibrillation-induced electrical remodeling on atrial electro-mechanics–insights from a 3d model of the human atria. PloS one **10**(11), 0142397 (2015)

[83] Alonso, S., Bär, M., Echebarria, B.: Nonlinear physics of electrical wave propagation in the heart: a review. Reports on Progress in Physics **79**(9), 096601 (2016)

[84] Pagani, S., Dede', L., Frontera, A., Salvador, M., Limite, L., Manzoni, A., Lipartiti, F., Tsitsinakis, G., Hadjis, A., Della Bella, P., et al.: A computational study of the electrophysiological substrate in patients suffering from atrial fibrillation. Frontiers in Physiology **12**, 673612 (2021)

[85] Corti, M., Zingaro, A., Quarteroni, A.M., et al.: Impact of atrial fibrillation on left atrium haemodynamics: A computational fluid dynamics study. Computers in Biology and Medicine **150**, 106143 (2022)

[86] Suga, H.: Importance of atrial compliance in cardiac performance. Circulation research **35**(1), 39–43 (1974)

[87] Hoit, B.D., Shao, Y., Tsai, L.-M., Patel, R., Gabel, M., Walsh, R.A.: Altered left atrial compliance after atrial appendectomy. influence on left atrial and ventricular filling. Circulation research **72**(1), 167–175 (1993)

[88] Nagano, T., Arakawa, M., Tanaka, T., Yamaguchi, M., Takaya, T., Noda, T., Miwa, H., Kagawa, K., Hirakawa, S.: Diastolic compliance of the left atrium in man: a determinant of preload of the left ventricle. Heart and vessels **5**(1), 25–32 (1989)
[89] Inciardi, R.M., Rossi, A.: Left atrium: a forgotten biomarker and a potential target in cardiovascular medicine. *Journal of Cardiovascular Medicine* **20**(12), 797-808 (2019)

[90] Aguirre, A., Codina, R., Baiges, J., Castañar, I.: Stress-displacement stabilized finite element analysis of thin structures using solid-shell elements, part ii: Finite strain hyperelasticity. *Finite Elements in Analysis and Design* **236**, 104179 (2024)

[91] Trayanova, N.A., Prakosa, A.: Up digital and personal: How heart digital twins can transform heart patient care. *Heart Rhythm* **21**(1), 89-99 (2024)

[92] Karakasis, P., Antoniadis, A.P., Theofilis, P., Vlachakis, P.K., Milaras, N., Patoulias, D., Karamitsos, T., Fragakis, N.: Digital twin models in atrial fibrillation: Charting the future of precision therapy? *Journal of Personalized Medicine* **15**(6), 256 (2025)

[93] Saiz-Vivó, M., Mill, J., Iriart, X., Cochet, H., Piella, G., Sermesant, M., Camara, O.: Digital twin integrating clinical, morphological and hemodynamic data to identify stroke risk factors. *npj Digital Medicine* **8**(1), 369 (2025)

[94] Hoit, B.D.: Left atrial size and function: role in prognosis. *Journal of the American College of Cardiology* **63**(6), 493-505 (2014)

[95] Sel, K., Osman, D., Zare, F., Masoumi Shahrbabak, S., Brattain, L., Hahn, J.-O., Inan, O.T., Mukkamala, R., Palmer, J., Paydarfar, D., et al.: Building digital twins for cardiovascular health: From principles to clinical impact. *Journal of the American Heart Association* **13**(19), e0131981 (2024)

[96] Nishimura, R.: Constrictive pericarditis in the modern era: a diagnostic dilemma. *Heart* **86**(6), 619-623 (2001)

[97] Vedula, V., George, R., Younes, L., Mittal, R.: Hemodynamics in the left atrium and its effect on ventricular flow patterns. *Journal of biomechanical engineering* **137**(11), 111003 (2015)

[98] Koizumi, R., Funamoto, K., Hayase, T., Kanke, Y., Shibata, M., Shiraishi, Y., Yambe, T.: Numerical analysis of hemodynamic changes in the left atrium due to atrial fibrillation. *Journal of biomechanics* **48**(3), 472-478 (2015)

[99] Masci, A., Barone, L., Dedè, L., Fedele, M., Tomasi, C., Quarteroni, A., Corsi, C.: The impact of left atrium appendage morphology on stroke risk assessment in atrial fibrillation: a computational fluid dynamics study. *Frontiers in physiology* **9**, 1938 (2019)

[100] Feng, L., Gao, H., Griffith, B., Niederer, S., Luo, X.: Analysis of a coupled fluid-structure interaction model of the left atrium and mitral valve. *International
journal for numerical methods in biomedical engineering **35**(11), 3254 (2019)

[101] Bucelli, M., Zingaro, A., Africa, P.C., Fumagalli, I., Dede', L., Quarteroni, A.: A mathematical model that integrates cardiac electrophysiology, mechanics, and fluid dynamics: Application to the human left heart. International journal for numerical methods in biomedical engineering **39**(3), 3678 (2023)

[102] Bucelli, M., Dede, L.: Coupling models of resistive valves to muscle mechanics in cardiac fluid-structure interaction simulations. International Journal for Numerical Methods in Biomedical Engineering **41**(12), 70119 (2025)

[103] Ruz Núñez, O.: Mathematical modeling and numerical simulation of left heart hemodynamics with fluid-structure interaction. PhD thesis, Sorbonne Université (2025). https://theses.hal.science/tel-04973025

[104] Sundnes, J., Lines, G.T., Cai, X., Nielsen, B.F., Mardal, K.-A., Tveito, A.: Computing the Electrical Activity in the Heart. Monographs in Computational Science and Engineering, vol. 1. Springer, Berlin, Heidelberg (2006)

[105] Neu, J.C., Krassowska, W.: Homogenization of syncytial tissues. Critical reviews in biomedical engineering **21**(2), 137–199 (1993)

[106] Courtemanche, M., Ramirez, R.J., Nattel, S.: Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. American Journal of Physiology-Heart and Circulatory Physiology **275**(1), 301–321 (1998)

[107] Dössel, O., Krueger, M.W., Weber, F.M., Wilhelms, M., Seemann, G.: Computational modeling of the human atrial anatomy and electrophysiology. Medical & biological engineering & computing **50**, 773–799 (2012)

[108] Collin, A., Gerbeau, J.-F., Hocini, M., Haïssaguerre, M., Chapelle, D.: Surface-based electrophysiology modeling and assessment of physiological simulations in atria. In: Functional Imaging and Modeling of the Heart: 7th International Conference, FIMH 2013, London, UK, June 20-22, 2013. Proceedings 7, pp. 352–359 (2013). Springer

[109] Trayanova, N.A., Lyon, A., Shade, J., Heijman, J.: Computational modeling of cardiac electrophysiology and arrhythmogenesis: toward clinical translation. Physiological Reviews **104**(3), 1265–1333 (2024)

[110] Ho, S., Sánchez-Quintana, D.: The importance of atrial structure and fibers. Clinical Anatomy: The Official Journal of the American Association of Clinical Anatomists and the British Association of Clinical Anatomists **22**(1), 52–63 (2009)