A cell-based framework for modeling cardiac mechanics

Cardiomyocytes are the functional building blocks of the heart—yet most models developed to simulate cardiac mechanics do not represent the individual cells and their surrounding matrix. Instead, they work on a homogenized tissue level, assuming that cellular and subcellular structures and processes scale uniformly. Here we present a mathematical and numerical framework for exploring tissue-level cardiac mechanics on a microscale given an explicit three-dimensional geometrical representation of cells embedded in a matrix. We defined a mathematical model over such a geometry and parametrized our model using publicly available data from tissue stretching and shearing experiments. We then used the model to explore mechanical differences between the extracellular and the intracellular space. Through sensitivity analysis, we found the stiffness in the extracellular matrix to be most important for the intracellular stress values under contraction. Strain and stress values were observed to follow a normal-tangential pattern concentrated along the membrane, with substantial spatial variations both under contraction and stretching. We also examined how it scales to larger size simulations, considering multicellular domains. Our work extends existing continuum models, providing a new geometrical-based framework for exploring complex cell–cell and cell–matrix interactions. Supplementary Information The online version contains supplementary material available at 10.1007/s10237-022-01660-8.


Introduction
Myocardium contraction is known to be affected both by subcellular (Borbély et al. 2005;Azeloglu and Costa 2010) and extracellular (Qin et al. 2007;Deckx et al. 2019) mechanisms. Cellular geometrical configurations have been demonstrated to be important for key mechanical features (Stein et al. 2011;Humphries et al. 2017). Existing models of cardiac tissue mechanics are usually based on extensive homogenization-useful for many purposes, yet unsuitable for capturing smaller-scale effects and interactions. These are widely used on tissue and organ level, see e.g., Guccione et al. (1991); Holzapfel and Ogden (2009), and have been used extensively for interpreting clinical data (Xi et al. 2012;Sack et al. 2016;Finsberg et al. 2018). One of the limitations of this approach, however, is that intracellular and extracellular processes are assumed to take place everywhere instead of being organized into discrete structures.
The cells and their extracellular material each have unique biochemical constituents and structure (Fomovsky et al. 2010;Avazmohammadi et al. 2019) and can be expected to have different mechanical properties. Furthermore, the force that drives cardiac contraction is only generated within the cells and not in the surrounding matrix. There are models that describe how this force is generated from electromechanical subprocesses on a sarcomere level (Land et al. 2017;Rice et al. 2008). When coupled to spatially resolved models, these are usually assumed to scale homogeneously through the tissue without taking into account the cellular geometries. Higher-resolution imaging techniques have been developed, see e.g., Pinali and Kitmitto (2014); Bensley et al. (2016), making it possible to extract exact geometrical representations of how cells are embedded in the extracellular matrix. To make use of this information in cardiac modeling, there is a need for developing spatially resolved modeling frameworks that explicitly capture both the cells and their surroundings.
Some mechanical continuum-based models that do take into account these explicit geometries have been developed, considering single isolated cells. These are modeled using a hyperelastic and isotropic (Tracqui et al. 2008;Okada et al. 2005;Garcia-Canadilla et al. 2017;Lenarda et al. 2018), or anisotropic material model (Tracqui and Ohayon 2009;Ruiz-Baier et al. 2014;Gizzi et al. 2015). The impact of the cell membrane and the extracellular matrix can be explicitly implemented using a combination of Dirichlet and Neumann boundary conditions (Ruiz-Baier et al. 2014;Gizzi et al. 2015;Lenarda et al. 2018). In Lenarda et al. (2018), cell-cell interactions were investigated by considering two cells connected through a surface representing a gap junction. The study presented in Tracqui et al. (2008) considered a single cell embedded in a substrate, investigating the impact of substrate stiffness on intracellular dynamics. None of these, however, include the extracellular matrix as a part of the domain.
Going somewhat broader, there are models for contractile skeletal cells considering cell-matrix interactions. An analytical model was presented in Sharafi and Blemker (2011), backed up with numerical simulations considering up to nine cells connected in a bundle, including an endomysium layer separating the cells from each other. Circular, elliptical and spherical geometries were used in multiple subsequent models, see e.g., Abhilash et al. (2014); Liang et al. (2016); Sopher et al. (2018); Mann et al. (2019), considering single cells and pairs of cells embedded in different kinds of polymer matrices. These works report that considerable amounts of force transmission occur through shear stresses and prominent force chains. Although all of these models have been developed for general contractile cells, the main findings most likely hold for cardiomyocytes as well.
Like mechanics, cardiac electrophysiology has often been modeled in a homogenized manner. Some of these models have been further refined into models that take into account the cells explicitly represented in the domain, surrounded by the extracellular matrix (Hogues et al. 1992;Stinstra et al. 2010;Tveito et al. 2017). Following the terminology in Tveito et al. (2017), such a model can be referred to as an EMI model-a model that separates the domain into an extracellular subdomain, a membrane and an intracellular subdomain.
In this work, we present the extension and characterization of the model presented in Telle et al. (2021), in which the extracellular and the intracellular subdomains are explicitly represented in the geometry. The main purpose behind this model is to capture interactions arising due to differences in structures and properties on a microscale, considering smaller tissue samples built up from individual cells, which potentially could differ in their geometries, material properties, or contraction dynamics. Here, we use experimental data to parametrize our model, considering stretching and shearing experiments. We were able to capture the full orthogonality of cardiac tissue through the geometry rather than imposing it in the strain energy function. Using this parametrized model as a baseline, we explored the parameter space subject to fiber direction stretch and contraction. Utilizing high-performance computing (HPC), we are able to move from the single cells to multicellular domains, representing small cubical tissue samples.

The mathematical framework
Following the geometrical framework presented in Tveito et al. (2017) and employed in Telle et al. (2021), we consider a three-dimensional domain consisting of two volumes-the cells and the surrounding matrix. The intracellular space ( Ω i ) is surrounded by the extracellular space ( Ω e ), separated by a surface ( Γ ) representing the cell membrane; see Fig. 1. Geometrically Ω i and Ω e are both closed volumes, while Γ represents the intersection Ω i ∩ Ω e , where again Ω i and Ω e represent the boundaries of Ω i and Ω e . We let Ω = Ω i ∪ Ω e denote the whole domain and Ω the outer boundary.
Using a fully Lagrangian formulation, let u denote the displacement of a given point in the domain Ω , defined by u ∶= x − X , where x is a point in the current configuration and X a point in the reference configuration. Let F ∶= ∇u + I denote the deformation gradient, where I is the identity tensor. The stress-strain relationship in the material is given by a strain energy function (F) , and the first Piola-Kirchhoff stress tensor by P ∶= (F) F . An equilibrium solution, where all forces are balanced, can be found by solving Fig. 1 Subdomains Ω i and Ω e and their boundaries. Schematic drawing of the subdomains; the intracellular subdomain Ω i is surrounded by a surface Γ representing the membrane, separating it from the extracellular subdomain Ω e . The whole domain is surrounded by an outer boundary Ω . and denote normal vectors of surfaces Γ and Ω . In 3D, Ω i and Ω e are volumes, separated by a surface Γ over the whole domain Ω , subject to imposed boundary conditions.
We assumed continuity of displacement and stresses across the membrane, mathematically expressed as and on Γ where and denote normal vectors of Ω i and Ω e respectively.
To incorporate contraction of the cells, we used an active strain approach, as described, e.g., in Ambrosi and Pezzuto (2012); Rossi et al. (2012). In this approach, the deformation gradient F is decomposed in a multiplicative manner, F = . The active component, , quantifies the strain caused by the cell contraction over time, while = F −1 gives the passive elastic deformation. Following Ruiz-Baier et al. (2014) and Rossi et al. (2012), we let be transversely isotropic, given by where is a scalar time-dependent function.
In our model, active contraction was only prescribed within Ω i , i.e., within the cells. We imposed the discretization by letting vary over time within Ω i while being set to zero in Ω e : For simplicity, we took i to be a scalar function dependent on time only.
For the strain energy function , we defined different strain energy functions i and e for the intracellular and extracellular spaces. We then combined these into one common strain energy function, For both domains, we used a hyperelastic invariant-based strain energy function, based on the one proposed in Holzapfel and Ogden (2009). These were based on the invariants I 1 and I 4f , given by Here tr denotes the trace operator, C ∶= J −2∕3 T the modified isochoric Cauchy-Green deformation tensor, and the fiber direction, i.e., the longitudinal direction of the cells.
For the intracellular subdomain, we let Here ‖ ⋅ ‖ + denotes a conditional term, given by ‖I‖ + = max(I, 0) . The parameters a i , b i , a if and b if determine the cellular stiffness. The second component gives significant increasing stiffness in the fiber direction subject to stretching, and none under compression. For the extracellular domain, we used an isotropic formulation, given by Again, the parameters a e and b e determine the stiffness of the matrix surrounding the cells. Other variants for this strain energy functions were also explored. In particular, we tried including shear terms and terms increasing the stiffness in transverse (sheet and normal) directions. These additional terms were, however, found to be redundant.
Cardiac tissue is known to be fully orthotropic, with the sheet direction determined by distinct alternating layers of cells and perimysium (Holzapfel and Ogden 2009;Costa et al. 1999). Rather than imposing this in the strain energy function, we let the full orthotropy of myocardium be imposed by including layers in the geometry used in the simulations-see detailed description in Sect. 2.3.
In this work, we assumed both subdomains to be incompressible, i.e., we required for all X in Ω during deformation. This restriction was imposed using a Lagrange multiplier, see full derivation in Holzapfel (2000). This gives us a modified strain energy function, where p can be interpreted as the hydrostatic pressure.

Weak form and implementation details
To solve the above systems of equations numerically, we used the finite element method. Here, we solved (1) in the weak sense by solving the following problem: Find displacement u and the hydrostatic pressure p such that for all test functions v and q from suitable test spaces. A full derivation can be found in Appendix A.1. The numerical experiments were implemented using FEniCS (Alnaes et al. 2015) (version 2019.1), using a Taylor-Hood discretization ( P 2 -P 1 ) to represent the displacement and the pressure, respectively.
Equation (13) is solved as a stationary problem for each time step. Before we solve the problem for a given step, either boundary displacement (for stretch/shear experiments) or active tension (5) (for contraction experiments) values are updated. The results from the previous step are used as an initial guess for the next step, with zero values being used for the very first step. Each such step can hence be considered a continuation step, for which Newton's method is used to find an equilibrium solution.

Geometries and meshes
For our experiments, we considered meshes representing a single cardiac cell embedded in an extracellular matrix, as well as tiled meshes consisting of multiple cells; see Fig. 2.
We take the fiber direction f 0 to be aligned with the cells in the longitudinal direction. Following the convention for (13) ∫ Ω P ∶ ∇v + q(J − 1)dX = 0 several tissue-level models, we define sheet and normal directions s 0 and n 0 , perpendicular to the fiber direction and each other. We follow the usual convention of defining the microstructural unit vectors in the reference configuration using a subscript 0, and the corresponding vector in the current configuration without the subscript.
The cell itself has a cylindrical shape, having a length of 102 μm (100 μm + 1 μm for each of the connections) and a diameter of 18 μ m, with rounded edges at each end; see Fig. 2 (c, f). The extracellular material was added around to form a box, such that the matrix surrounds the cell, with the thinnest layer being 1 μm in the sheet direction and 3 μm in the normal direction. The thicker layers in the normal direction were added to resemble layers of perimysium. Together with the cell orientations, these give rise to the characteristic local three-dimensional structure of myocardium. These proportions give us a total volume of 24.42 ⋅ 10 3 μm 3 for Ω i , and 24.54 ⋅ 10 3 μm 3 for Ω e for a single-cell geometry.
The regular cubical shape of the domain was explicitly developed for the purpose of extending the framework to multicellular domains representing varying numbers of cardiac cells. Here, copies of the single-cell geometry were simply tiled next to each other in width, length, and height. The cells were connected in the fiber direction, sharing a common surface in the mesh at the connections.
Properties of meshes of various resolutions for a single-cell geometry are reported in Table 1 The single-cell mesh was used as a base for the tiled mesh, combined into a common geometry by being tiled in a grid fashion. Note that the padding in the normal direction n is larger than the padding in the sheet direction s , meant to resemble layers of perimysium. The dimensions are indicated in (c) and (f) 1 3 used for convergence studies, and for most other experiments, the middle (5.0 μm ) mesh was used. Tiled meshes, generated using the 5.0 μm as a base, were used for scaling experiments. Properties of these meshes are listed in Table 2. The single-cell mesh (5.0 μm ) can be seen in Fig. 2a-c, and a corresponding tiled mesh is displayed in Fig. 2d-f. All meshes are tetrahedral and were generated using Gmsh (Geuzaine and Remacle 2009).

Deformation modes
Virtual stretching and shearing experiments were used first to parametrize our EMI model. This parameterized model was then used to explore convergence, sensitivity, spatial variation, and scalability. For these simulations, we considered nine distinct deformation modes-stretching in the fiber, sheet and normal directions (FF, SS and NN), and shearing in all combinations of fiber/sheet/normal directions (FS, FN, SF, SN, NF, NS). See Fig. 3c for a schematic overview. For each deformation mode, displacement on pairwise opposite surfaces (see Fig. 3a) was enforced using Dirichlet boundary conditions while the other surfaces were allowed to move freely. A complete mathematical description of all boundary conditions can be found in Appendices A.2 and A.3. Cellular contraction was modeled using a precomputed time-dependent active strain transient, as displayed in Fig. 3d. For these experiments, the traction along every surface was set to zero. Rigid motion was avoided by enforcing orthogonality of the solution with respect to the kernel using Lagrangian multipliers, as explained, e.g., in Kuchta et al. (2016). See Appendix A.4 for more details.

Reported quantities
For many experiments, we report values derived from the Green-Lagrange strain tensor, given by and values derived from the Cauchy stress tensor, given by In particular, we considered the normal components along the fiber, sheet, and normal directions of each tensor, i.e., E ff , E ss and E nn for the strain, and ff , ss and nn for stress. Except from when displayed as spatial plots, these were taken as averaged quantities over the separate subdomains. These were calculated as volume integrals, given by and Here is equal to , or , in the reference configuration.
Subdomain Ω j is either the (combined) intracellular space Ω i , the intracellular space in a single cell ( Ω i,k for k = 1...N , considering N cells) or the extracellular space Ω e . The integrated quantities were computed using fourthorder Gaussian quadrature, which coincided with second and third order for the strain values. For visualization of the spatial distributions, these were projected to a discontinuous Lagrangian function space ( DG 2 ), the discontinuity being helpful for capturing strain and stress components close to the membrane.

Parameter estimation
To parametrize our model, we used the same experimental data as explored in Kakaletsis et al. (2021), made publicly available by the authors (Kakaletsis et al. 2020). In their study, cardiac tissue blocks from sheep hearts were extracted, then shearing and stretching experiments were performed on each sample. We restricted the scope to the 11 samples taken from the left ventricle. In order to make the optimization tractable with our given computation resources, we only considered the positive values, i.e., from reported displacement zero and upwards.
For each sample, there were data on forces, in both normal and shear directions, and resulting displacement, as well as length, width and height measurements for each of the samples. We derived load and stretch values for each sample and used these for direct comparison to the model. We performed corresponding virtual stretch and shear experiments by deforming the domain using the same range of stretch and shear magnitudes used in the experimental data. For each stretch value and each deformation mode M, the total load on the surface, imposed by the applied Dirichlet boundary conditions, was calculated by Here e is either the mesh normal direction or the shear direction, depending on the tracked component. S i is either S 2 , S 4 or S 6 , depending on the deformation mode (see again Fig. 3), and N denotes the normal vector to this surface.
The parameter estimation was then formulated as an optimization problem, expressing the difference between the experimental load values and the virtual load values as an L 2 norm. We performed the optimization using SciPy's minimize function, allowing material parameters a i , b i , a e , b e , a if and b if to vary within the interval [0.01, 40].

Parameter sensitivity
We next explored the sensitivity of each parameter by computing Sobol indices. This analysis was done for all deformation modes used in the stretching and shear experiments, tracking load values of interest, as well as for averaged stresses across the entire domain under fiber direction stretch experiment FF and contraction. The sensitivity analysis was performed using the Python library SALib (Herman and Usher 2017;Iwanaga et al. 2022).
For each of these cases, sensitivity analysis was performed by sampling the parameter space with N = 512 values, allowing each material parameter to vary in the interval [0.1, 30]. We considered D = 6 material parameters, which resulted in N(D + 2) = 512(6 + 2) = 4096 different parameter combinations per deformation mode. For each parameter combination, we simulated virtual stretch up to 10% and virtual shear up to 40%, dictated by the range of strain values of the experimental dataset used for the parameterization. Contraction was simulated to a peak of approximately 20% shortening. These experiments were all performed using a single-cell geometry. Finally, based on the resulting load and stress values, we computed the first order and total sensitivity indices.

Convergence
To test the convergence properties of our problem, we generated a set of meshes with decreasing maximum mesh . element size (circumradius), as listed in Table 1 (Max c).
Here h max and h min are defined as the largest and smallest cell diameter in the mesh. The cell diameter is defined as twice the circumradius. We performed convergence experiments by tracking normal stress components for each subdomain, as given in (17), considering fiber direction stretch and contraction. These were also performed using a single-cell geometry.

Spatial plots
To explore spatial distributions we performed fiber direction stretch and contraction experiments, processed in Paraview. The fiber direction stretch was simulated up to 10% stretch. For the contraction experiment, we simulated half a cardiac cycle-and report spatial values for the maximal contracted state (138 ms). For the contraction experiments, we also compared individual components of stress and strain values averaged over each subdomain. We considered both single cell meshes and tiled meshes, consisting of 3 × 3 × 3 cells.

Scalability
Finally, we investigated the performance and scalability of our current solver implementation with the aim of tackling larger, tissue-scale problems using meshes of tiled cells. For these experiments, we stretched the domain by 15 % in the fiber direction using ten continuation steps while varying the problem size and number of CPUs used to perform the calculations. We considered one weak scaling problem, where we used tiled meshes such that there is one cardiac cell per CPU, and one strong scaling problem, where a fixed mesh representing 4 × 4 × 4 cells was used. Properties of meshes used in these experiments are listed in Table 2.
The current solver is based on using Newton's method combined with a distributed-memory parallel direct solver provided by SuperLU_dist (Li and Demmel 2003). A known limitation of this approach is the considerable memory usage that is associated with the LU factorization stage of the direct solver (see, e.g., Dongarra et al. (1998)). To increase the amount of available memory, in order to make the problem feasible to solve, we employed six compute nodes for the strong scaling experiments. Each compute node consisted of two Intel Xeon Gold 6138 CPUs with 40 CPU cores and 192 GiB of memory. Consequently, we limited our current experiments to at most 240 CPU cores and a total of 1 152 GB of memory. We configured SuperLU_dist to use the serial MC64 algorithm to compute a row permutation. An alternative parallel algorithm for based on Approximate Weight Perfect Matching (AWPM) (Azad et al. 2020) was tested, but found to be significantly slower for our case.

Parameter estimation
Optimized parameter values for (9) and (10) found for each experimental sample, as well as the average across all of them, are listed in Table 3. The experimental data are plotted together with the model fit, using average parameters, in Fig. 4. We observe that most of the load values fall within the range given by the experimental data-in particular, the load values calculated for the FF, FS, and FN experiments all are close to the middle of the range defined by the experimental data. On the other hand, the load values for the sheet and normal direction stretch (SS, NN) are both somewhat too high.

Parameter sensitivity
Sensitivity analysis was performed using variance-based sensitivity analysis, with the resulting Sobol indices displayed in Fig. 5. We here report the first order ( S 1 ) and total Table 3 Optimized material parameters for the strain energy function (6) to experimental data (Kakaletsis et al. 2020), for samples 1-11 The average values are highlighted in bold  Table 3, is dis-played in red. The schematic drawings (top right corners) display displacement direction (white), normal (pink/white) and shear (white) components; see more detailed explanation in Fig. 3 ( S T ) sensitivity, and for the sake of brevity we do not show the confidence intervals. Overall, the confidence intervals span a range up to 0.23 for the load experiments, up to 0.09 for the fiber direction stress experiments, and up to 0.12 for the contraction stress experiments. Across all experiments, we see that the general patterns of the first-order sensitivity are repeated in the total sensitivity, and more so for the averaged stress values than for the load values. For the load experiments, we see that most of the modes are most sensitive to changes in the exponential b e parameter, followed by the exponential b i parameter. Prominently, the exceptions are the three stretch experiments, which, in particular, have low sensitivity for these two parameters. Instead, the a if parameter appears to matter most for the fiber direction stretch (FF), the a i and a e parameters most for the sheet direction stretch (SS), and a i followed by a e for the normal direction stretch (NN).
Considering the fiber direction stress experiment, we can see that the intracellular fiber direction stress ( ff across Ω i ) is most sensitive to the parameter a if , as for the corresponding load experiment. The extracellular sheet direction stress ( ss across Ω e ) is by far most sensitive to the a i parameter. Across all tracked stresses, the b i parameter matters the least.
Finally, for the contraction stress experiment, across all tracked stresses, the b e parameter matters the most, seconded by a e . Across the other parameters, the sensitivity is, in comparison, marginal. This holds for both the first order and the total sensitivity, although we can see a small enhancement in all values for the total sensitivity, i.e., having somewhat higher magnitudes.

Convergence
Results from the convergence studies are displayed in Fig. 6, reporting subdomain stress values for different mesh resolutions. Many of the tracked components appear converged at a mesh resolution with a maximum element size of 5.0 μm . However, under contraction, component ff appears to still vary somewhat going from the mesh with maximum element size of 2.5 μm to the mesh with an element size of 1.25 μm.

Single cell versus multicellular domains
We compared the load and stress results for a single cell, which was used for the parametrization, and 3 × 3 × 3 cells. The results are displayed in Fig. 7. We note that across all stretching and shearing experiments, considering load values, the results are similar for both geometries. For stresses in the stretching experiment (subplots b, e), most quantities are comparable (in absolute values; relative values vary somewhat more), while for contraction, there is a notable difference. Here, the averaged fiber direction stress given by ff across Ω i , almost doubles from 7.87 kPa at peak for the single-cell domain to 12.75 kPa for the multicellular domain. Similarly, but reversed in magnitude, the fiber direction stresses given by ff across Ω e , decreases from −8.17 kPa at peak for the single-cell domain to −12.85 kPa for the multicellular domain. The individual contributions of the fiber direction stresses of each cell are also plotted in Fig. 8. Here, we can see that the fiber direction stresses are higher in all cells in the tiled mesh compared to the single-cell mesh. We also note that this increase is highest in the middle cells, i.e., the cells surrounded by cells on two sides in the fiber direction. Along the sheet and normal directions there is much less variation between the different cells.

Spatial distributions of strain and stress
Spatial distributions of the pressure, strain and stress values for fiber direction stretching and contraction experiments are presented in Figs. 9, 10, 11, and 13. All spatial plots are deformed according to the respective displacement fields. For the contraction experiment for the single cell, we also include plots of strain and stress values averaged over each subdomain in Fig. 12 For the stretching experiments, the values are plotted at maximum stretch, while for contraction, this is the state in which the cell reaches the highest contraction, corresponding to the peak of the active strain in Fig. 3d. We report the normal components of E (14) and (15). Note that, following the continuity assumption (3) the normal stresses are continuous across the membrane. For stresses in the other directions, however, we predominantly see a clear discontinuity here.
We note that for both stretching and contraction experiments, the pressure field, as displayed in Fig. 9, contributes to about half of the stress in the fiber direction, consistent in sign and magnitudes. The pressure contribution to the stress values is the same along the diagonal entries (per definition, see Eq. (15)), while being zero for all shear entries. However, we can observe that for most of the domain the total stress distributions in the sheer and normal directions are mostly zero-valued, implying that the stress arising from (6) and the stress from the pressure have opposite signs and mostly cancel each other out. a b Fig. 6 Convergence experiments-normal stress components. The plots display the effect of mesh refinement (see Table 1), measured by the change in integrated stress values subject to fiber direction stretch (a) and contraction (b). The lightest line corresponds to the coarsest mesh (20.0 μm), while the darkest line corresponds to the finest mesh (1.25 μm). For both deformation modes we display the normal stress components of (15) averaged over intracellular subdomain Ω i and the extracellular subdomain Ω e , respectively, as given by (17). The curves marked with stars display the resolution used for most of the experiments in this paper

Stretching experiments
As displayed in Fig. 10a and 11a, we can observe fairly uniform strain values everywhere in the domain for the stretching experiments. There is, however, a subtle discontinuity along the membrane for components E ss and E nn (visible inside the quarter removed, to the right and on the top, respectively). The stress values, on the other hand, vary considerably. Considering the fiber direction stretch, even though the E ff values along and inside the cell are almost similar, ff takes significantly higher values in the intracellular subdomain. For the tiled meshes, the patterns observed mostly generalize what we see for the single-cell experiments. However, for E ss , and to some degree for E nn , we see that the area between the cells is being stretched considerably in their respective directions-without any corresponding strong response in stress values. In the movies (Movie 1 and Movie 2), we can see the evolution of strain and stress patterns over different stretching levels. In particular, we see strain patterns develop quite early, while, on the scales plotted, the stress patterns develop later. We also see that the ff is highest along the middle of the cells. This remains true for all steps, but is outside the scale depicted for the final steps in Fig. 10 and 11.

Contraction experiments
In contraction, we observed large variations in all scalar fields ranging from − 0.3 to 0.3 in strain, and from −80 kPa to 40 kPa in stress values. We see this to some degree in the average traces presented in Fig. 12, with highest variation for stresses represented by ff . Here we also observe higher orthotropy for strain values than for stress values-the difference between E ss and E nn is more prominent than the difference between ss and nn .
As expected, and consistent with the average values plotted in Fig. 12, we observe mostly negative E ff and mostly positive E ss and E nn values for both domains. The strain values are fairly similar along and inside the cells. In contrast, ff has a clear discrete transition across the membrane, taking positive values in Ω i and negative values in the extracellular subdomain ( Ω e ). For the ss and nn components, we have low stress in most of the domain. However, there is a small region along the cell with positive values (barely visible as a red line) and a more prominent area at the end of the cell, which has negative values.
Tensor glyphs displaying strain and stress values for a single cell under contraction are plotted in Fig. 13. For intracellular strain values, the glyphs are elongated close to the membrane tangentially. This implies a more prominent expansion close to the membrane than in the middle, and is also visible for the E ss and the E nn components in Fig. 10somewhat more prominent for E nn than for E ss . For extracellular strain values, all glyphs are spherical and somewhat We observe far more prominent stress values in the extracellular subdomain than in the intracellular subdomain. These are mainly concentrated close to the membrane, also in a tangential pattern. We can observe that these are flat glyphs, indicating two prominent principal components. One of these components coincides with the fiber direction, while the other follows the cell shape circularly. These glyphs, highly prominent in this plot, correspond exactly to the subtle red lines barely visible for components ss and ff in Fig. 10. We observe flat glyphs at each end of the cell for the intracellular subdomain, also orientated perpendicular to the fiber direction. In the middle section of the cell, however, we have elongated glyphs oriented in the fiber direction, indicating that the main contribution to all the stresses is in this direction.
We observe general higher strain and stress values for the tiled meshes compared to the single-cell experiment. In particular, we have increased magnitude for ff values, both for the intracellular and extracellular space. We also see that strain values between the cells changes-here, they obtain negative ( E ff ) or positive ( E ss , E nn ) values, following the cells surrounding them.
Both color maps and tensor glyphs are also included in movies for contraction, see Movie 3 and Movie 4. Many of the same observations as outlined above can be done here, but we also see that we start with a prominent strain concentration close to the membrane, while the stress concentration in the extracellular domain emerges later on. We also note that the prominent intracellular ff component is emerging first and is taking highest values close to the connection between the cells.

Scalability
The performance and memory consumption for the weak and strong scaling experiments are displayed in Tables 4 and 5. The solver required 43-46 Newton iterations in total for every case.
The parallel efficiency and memory efficiency for the scaling experiments are shown in Fig. 14. For the weak scaling experiments, we report the baseline time over the measured time for each measurement and the baseline memory divided by the measured memory. For the strong scaling experiments, we report the baseline time, times the number of nodes (6), divided by the measured time and the number of processes, and baseline memory divided by the measured memory.   (14), averaged over the intracellular and extracellular subdomains Ω i and Ω e -see Equations (16) and (17) over the first 500 ms of a cardiac cycle For the weak scaling experiments, one would ideally hope for the CPU time to remain constant as the problem size and number of CPUs increase proportionally. Unfortunately, this is not the case, and both the CPU time and memory increase significantly as the problem size (number of cardiac cells) increases, despite the work per processor remaining constant. These results indicate that some parts of the direct solver do not scale well.
In the case of the strong scaling, there is clearly some speedup from increasing the number of CPUs, where in particular we observe fairly good scaling of the LU factorization stage. Most of the time is initially spent in this stage, which goes from about 420-480 seconds per Newton iteration when   Table 4 and Table 5 using 1 CPU per node to about 25 seconds per Newton iteration when using 32 CPUs per node. However, the solver must also compute row and column permutations. The row permutation, computed using a serial algorithm, takes about 14 seconds per Newton iteration regardless of how many CPUs are used. The calculation of the column permutation requires more time as the number of CPUs per node increases, taking up to 26 seconds per Newton iteration. Ultimately the computation of row and column permutations limits the scalability of the solver.

Discussion
We have presented a cardiac model for single cells and small collections of cells, where these cells are embedded in a matrix with differing material properties. It takes into account an explicit geometrical cell configuration that allows for refinement of cardiac tissue mechanics in a more physiologically relevant manner than one can achieve by simply using finer meshes and a homogenized tissue model. This work extends the mechanical aspect of cellbased modeling of single cells (Ruiz-Baier et al. 2014;Gizzi et al. 2015;Garcia-Canadilla et al. 2017;Lenarda et al. 2018;Tracqui et al. 2008) by including the extracellular material. Our results indicate that this inclusion is important for the quantification of intracellular stresses.

Model parameterization
We used experimental data to parameterize our model, matching six parameters to data from nine different deformation modes, finding a reasonable fit. There is a fairly large variance in the optimal parametrization of each individual sample, which may reflect on variations in the underlying experimental data. Some general patterns observed here is that typically either the a e or the a i parameter is high, and the other low. We also see that the fiber direction stiffness parameters a if and b if vary across almost the whole range of allowed parameters, which could arise from varying fiber dispersion in the samples. We also noted that starting the optimization from different starting values did not change the main patterns in the final results. There were certainly some differences indicating local minima, but when averaged out this resulted in no significant differences. As displayed in Fig. 7, the differences between the single cell geometry and the multicellular geometry in computed load during stretching and shearing simulations are minor. This indicates that performing parameter optimizations for passive behavior on a single-cell geometry can, within reason, be extrapolated to multicellular domains. We were able to capture differences in the stress-strain ratio in different directions transverse to the fiber direction by incorporating geometrical differences, i.e., by making the matrix thicker in the normal direction than in the sheet direction, rather than through additional terms in our strain energy function. The sensitivity analysis of separate intra-and extracellular stresses displays interesting dynamics between the two subdomains. The extracellular material parameters, a e and b e , proved to be most important for intracellular stresses. Considering the load values, based on the relative sensitivity, we see that both stretching and shearing experiments are important to perform and use in the parametrization to capture these two accurately. The sensitivity analysis also displays that across the stress values, the first order and total Sobol indices are fairly similar, indicating little interaction between the different parameters. For the load values we observed a larger difference between the first order and total indices, indicating a higher degree of interaction.
The formulation used for the intracellular strain energy function i (9) is very similar to the one used in other works for isolated cardiomyocytes (Tracqui and Ohayon 2009;Ruiz-Baier et al. 2014;Gizzi et al. 2015), and can physiologically be motivated by the mechanical contributions of the sarcomere structure found within the cells. Other models for cardiomyocytes (Tracqui et al. 2008;Okada et al. 2005;Garcia-Canadilla et al. 2017;Lenarda et al. 2018) have assumed isotropy in the intracellular domain, not assuming that there is any difference between stress-strain values in different directions. This could, in particular, be a good assumption for cells and cell collections developed in vitro, which as immature cells have a less developed sarcomere structure. For the extracellular strain energy function e (10) it could in particular be relevant to compare with Sharafi and Blemker (2011) and Zhang and Gao (2012). In Sharafi and Blemker (2011), a term enhancing the shear stiffness is included, while in Zhang and Gao (2012) the authors utilize an isotropic formulation, modeling both subdomains as Mooney-Rivlin materials. Inspired by the first one, we explored several formulations including additional shear components, as well as formulations with additional transverse stiffness. In all cases, however, these were found to be redundant-our optimization script found the corresponding parameters to be zero across most of the samples. The difference in norm found by the optimization script remained fairly similar with or without these extra terms, also for the non-zero cases. As such, we believe the dynamics captured with these potential extra terms is already captured in the current formulation. Further studies are, however, needed to better determine the most representative constitutive relationships in this model, for each subdomain.
Compared to the study in Gizzi et al. (2015), our parameters a i , b i , a if and b if are all stiffer. They considered a single cell, employing a similar model to the one used for the intracellular space. Data are obtained by curve fitting based on the polynomial strain energy function used in Tracqui and Ohayon (2009), which again is based on data published in other studies. Compared to the original study in Kakaletsis et al. (2021), from which we have the experimental data (for tissue samples), we observe an increase in the anisotropy for the intracellular space. We found parameters a if and b if to take values 19.83 kPa and 24.72, respectively. These are comparable to the parameters a f and b f as listed in the original paper, in their parameterization of a corresponding homogenized tissue-level model. Across their four strategies, the a f parameters ranged between 3.85 and 6.60 (kPa), and b f between 4.34 and 10.79. We see that the intracellular space has a significant increase in degree of anisotropy compared to the homogenized tissue-level model. This can be explained by the isotropic assumption for the extracellular space-to match tissue-level experiments with decreasing anisotropy in one domain, one needs to increase the anisotropy in the other.
The main limitation of the parameterization process might be that, although the model itself is designed for taking into account explicit cellular geometries, we are using an idealized one. In particular, we have a much larger volume ratio between the intracellular and extracellular space than realistic. Our intracellular to extracellular ratio is almost 1 to 1, while values in the literature report the volume ratio to be around 4 to 1 (Olivetti et al. 1991). This is likely to lead to an overestimation of a if and b if , as the stiffness is defined over a volume which is too small. In addition, we have added extra padding between every sheet of single cells. In reality there should be a few cells per sheet, surrounded by perimysium layers on both sides, which will have an impact differences between the sheet and normal directions. We also only connect the cells in the fiber direction, in an artificially regular manner. As displayed in Fig. 8, there is a huge degree of interaction between the cells in this direction and only marginal interaction in the other directions, tracking stress values under contraction. Connecting them in the sheet direction as well, which is, e.g., done in the electrophysiological model in Tveito et al. (2017), would most likely bring in more interaction in this direction as well. Capturing all of these closer-to-realistic aspects would be of major interest, bringing us closer to capturing real physiological features, but would also require us to use far more advanced meshing techniques.
Histological information is included in the data published (Kakaletsis et al. 2020). One could potentially make geometries based on these images, which would most likely lead to a significant improvement of the fit and new exciting results. The main problem here, however, is the size of the tissue samples. All of the samples were cut to have approximate dimensions 10 mm × 10 mm × 10 mm , for which one can estimate a couple of millions of cells, far more than the 128 cells we are, at max, considering in our study. An additional limitation of our optimization procedure is that we only consider the positive values, while the original data also included negative ones, i.e., compressive stress and shear in two directions. Including both sides could potentially lead to a better parameterization-in particular, because we have a conditional term in the fiber direction, dependent on whether we have stretching or compression. We also assumed continuity of stresses across the membrane (3), implying that the membrane has no stiffness. Incorporating some stiffness for the membrane would probably be more realistic and likely to affect the cell-matrix dynamics observed. It could potentially also give a better fit to the experimental data. For simplicity, however, we chose not to include it.
As described above, a study using a cell-based model, paired with stretch and image-based geometries and corresponding stretch/strain experiments on a smaller tissue sample would probably be the most rigorous and accurate way for proper separate characterizations of material properties of the cardiac cells and the matrix. It could also make sense to explore geometries generated artificially, capturing more realistic features, e.g., by allowing for a non-regular tiling pattern or matching the extracellular/intracellular volume ratio. Finally, it would be interesting to compare how the parametrization of our model would change when parametrized based on data from other experimental studies.

Toward a more physiologically relevant model
The results presented in this paper are highly dependent on the underlying geometrical cell configuration. The strain and stress plots display spatial patterns, following the cell geometries, with clear transitions along the membrane. Under contraction, this is, in particular, visible as a clear tangential-normal pattern. Here stress values in the extracellular domain are concentrated in directions tangential to the membrane, perpendicular to prominent intracellular strain concentrations, oriented normal to the cell membrane. Through differences in the geometry in the sheet and normal direction, we were able to capture a fully orthogonal behavior of the tissue. This applies both for stretching and shearing experiments, as displayed in Fig. 4, as well as for stress and especially strain values under contraction, as displayed in Fig. 12.
Real geometries are much more complicated-for a single cell, the domain has a more complex shape, and for multiple cells, the cells self-organize in more complex, non-regular patterns. Images of cells can be used to construct more realistic geometries, as done in, e.g., Ruiz-Baier et al. (2014); Gizzi et al. (2015); Garcia-Canadilla et al. (2017);Tracqui et al. (2008) considering single cells, and to some degree in in Sharafi and Blemker (2011), considering a cross-section through nine cells. For the latter study, force values reported are within a comparable range to stresses observed in our study; however their setup is quite different from ours, making it hard to compare these directly.
Most of the models mentioned above include calcium dynamics, with experimental spatial-temporal calcium measurements driving the model. A limitation of our model is that it does not couple strain dependencies back to the active tension, which would be more physiologically correct. In particular it would be of interest to include this for comparison to the work performed by Tracqui et al. (2008). This was designed to match the experimental setup explored in Qin et al. (2007), and explains key sub-processes in cell-material dynamics subject to increasing stiffness of the surrounding material. A fully coupled model would, however, only matter for capturing active contraction dynamics, while the stretching and shearing experiments determining the material parameters would remain unaffected.
The proposed model has the ability to represent microscale heterogeneities in intracellular and extracellular matrix properties, e.g., by allowing the individual cells to have different properties. Given an image-based geometry of actual cardiac cells, far more realistic simulations could potentially be achieved. Extending our purely mechanical model to take into account calcium and sarcomere length dynamics would also bring in intracellular heterogeneities for the active tension, with resulting new strain and stress patterns. Measurements of subcellular mechanical structures and calcium measurements, as explored in, e.g., Reichardt et al. (2020); Garcia-Canadilla et al. (2017); Blatter et al. (2003); Deckx et al. (2019), could be used within such a framework for exploration of heterogeneities within a single cell or a small collection of a few cells. A key focus here could, for example, be to differentiate between heterogeneities in material properties versus heterogeneities in the activation parameters. An alternative approach, focusing on cell-cell and cell-matrix interactions, would be to keep the mechanical properties and calcium levels homogeneous within each cell, and instead utilize cell-based calcium measurements, as reported, e.g., in Jones et al. (2018). The numerical framework, as developed for computational resources presently available to us, is currently far from capable of capturing the whole heart. It might, however, be possible to use it either to capture the mechanics of in vitro cell collections, such as cardiac microtissues (Zhang et al. 2015), or to examine in detail small sections of the myocardium, to understand how local interactions may propagate up to tissue and organ levels.
A particular area of interest for the model is in disease modeling, targeting diseases affecting mechanics on this scale. Our model could, for example, be used to study the impact of hypertrophy, characterized by changes in the cell geometries (Göktepe et al. 2010). By explicitly modeling the intracellular space, one could look at changes to the myocyte during eccentric or concentric remodeling by altering the length and diameter of the cells. This could be useful in helping to delineate how mechanical triggers such as stress and strain drive these remodeling processes, by more accurately determining the stress and strain the myocyte experiences. Here it would be interesting to understand how remodeling-based changes to the cell geometries can normalize altered load or contribute to continued remodeling stimuli. Another application could be in modeling scarring, in which the cell structures in the damaged part change upon healing (Rog-Zielinska et al. 2016). In the scarred regions, fibroblast cells differentiate to myofibroblast, which have contractile properties (Baum and Duffy 2011). In principle, there is no reason our cell-based model should only work for cardiomyocytes; it would be fairly straightforward to extend the geometrical framework to including different cell types with different mechanical and contractile properties.
Hypertrophy and scarring are both common causes of cardiac fibrosis (Maulik and Mishra 2015;Hinderer and Schenke-Layland 2019), which also can be represented more explicitly in a cell-based framework than in traditional tissue-level models. This includes interstitial fibrosis, in which the matrix stiffness increases, or replacement fibrosis, in which cells are replaced by a collagen network (Hinderer and Schenke-Layland 2019). In our model, we could capture this by, respectively, increasing the extracellular material parameters and by replacing certain cells with more matrix in the geometry. For the first one, in particular, we have demonstrated that the cardiomyocyte stresses are highly dependent on the matrix stiffness, so even a small increase in matrix stiffness could be expected to have a large impact. Modeling this explicitly on a cell-based level could provide physiological insights and understanding of how these diseases work on the relevant scale. In a long-term perspective, following the development of more efficient simulations, our modeling framework may be well-suited to represent larger tissue samples of clinically relevant sizes. Such an explicit representation could, for example, be used as an alternative to the more common statistical representations of fiber dispersion.

Implications of continuity and discontinuity assumptions at the membrane
In our implementation of the model, we solved for displacement u and pressure p in a Taylor-Hood discretization space, assuming continuity of both in the whole domain. Furthermore, we took (6) to be discontinuous across the membrane, implying that the first Piola-Kirchhoff stress tensor can be discontinuous as well-except from normal to the membrane, where we assumed continuity (3). We also assumed incompressibility in both subdomains. Our approach is very similar to the one considered in Tracqui et al. (2008), in which three subdomains are considered (the substrate, the cell, and the nucleus), and the whole system is solved simultaneously. As in our case, they assume continuity of displacement and stresses across the membranes-but in all directions. In Farsad and Vernerey (2012), on the other hand, continuity is neither assumed for displacement nor stresses. They also use the extended finite element method (XFEM) rather than FEM. Here, the discontinuities at the membrane are represented by Heaviside functions instead of explicitly representing the intersection in the mesh. They observe stress concentrations close to the membrane, consistent with our results. Yet another alternative could be to employ a discontinuous Galerkin (DG) method, as outlined, e.g., in Ten Eyck et al. (2008). We also note that in the electrophysiology cell-based model in Tveito et al. (2017), they employ a splitting scheme in which the equations for the intracellular space, the membrane, and the extracellular space are solved separately. We note that an incompressible formulation is used in Ruiz-Baier et al. (2014); Lenarda et al. (2018) for isolated myocytes, and a nearly incompressible formulation is used in Zhang and Gao (2012) for skeletal cells, surrounded by an endomysium layer, for both subdomains.
In our case, we found that the base mesh combined with the Taylor-Hood discretization remains a reasonable choice. Using a continuous function space for the pressure is, in particular, a limitation of our work, as there is no physical reason that neither the pressure nor the pressure-dependent stresses should be continuous. As displayed in Fig. 9, the pressure fields are somewhat blurred out across the membrane-but this appears to be a fairly minor artifact. XFEM and DG methods could be used to capture the discontinuity, however, they include a more complex mathematical framework and are more expensive discretizations, which do not seem necessary for our work. A splitting scheme could potentially also be developed for the mechanical modelhowever, again, this would make the methods more complex without necessarily being more precise at this stage. Cardiac tissue is known to be compressible (Yin et al. 1996;Nolan and McGarry 2016), in which the cells, which primarily consist of water, are close to incompressible, while the matrix' volume changes significantly under pressure. The matrix is estimated to be 100-1000 times more compressible than the cells it surrounds (Dolega et al. 2021). Alternatively one could use a nearly incompressible formulation, as done in Telle et al. (2021), but with a much higher penalty parameter for the intracellular subdomain. High penalty parameters are, however, associated with locking (Hadjicharalambous et al. 2014;Karabelas et al. 2022), which both can lead to numerical instabilities and underestimation of variables of interest. For these reasons we chose an incompressible formulation. This challenge could, however, be overcome using other kinds of elements, as widely explored in (Karabelas et al. 2022). A more rigorous comparison between different numerical schemes would probably lead to the development of more efficient and accurate methods. In particular, if one wants to couple the mechanics with the underlying electrophysiology, the splitting scheme could be more appropriate. Alternative formulations for either incompressible or nearly incompressible formulations would be prudent to explore in the future, and it would in particular be interesting to see the impact of defining these differently in each subdomain.

HPC considerations
Through our scaling experiments, we explored limitations of the presented numerical solution approach both with respect to time and memory consumption. We observe reasonable speedup for some parts of the direct solver, such as for the LU factorization stage. However, other parts become increasingly costly, such as the computation of row and column permutations, which are needed for robustness and to reduce fill-in of the computed factorization. As the problem size increases, these parts seem to dominate the execution time.
As considered in, e.g., Whiteley (2017); Brune et al. (2015), iterative methods combined with suitable preconditioners might give better performance for nonlinear elasticity problems. In terms of problem size, our largest meshes lie somewhere in between the second largest and largest problem size considered in Whiteley (2017), indicating that our problem is comparable to their experiments in terms of degrees of freedom. Operator splitting schemes have proven to give significant speedup for the electrophysiological EMI model, as considered in Jaeger et al. (2021), and similar approaches would be interesting to investigate from a mechanical perspective as well. Hexahedral meshes have been demonstrated to be more accurate than tetrahedral meshes, which we have used, as reported in, e.g., Karabelas et al. (2022); Oliveira and Sundnes (2016) for related problems. From an efficiency perspective, this implies one could get better solutions with coarser meshes, which could be worth considering for future work.
Limited scalability obviously affects the generalization of our simulations-using our framework, we can only consider small collections of cells. Implementing preconditioners and operator splitting schemes are, however, not always straightforward and can be considered a separate extensive research question.
In the future, it would be prudent to investigate alternative solvers. It would also be interesting to see whether changing, e.g., the mesh structure would lead to different numerical properties, both with respect to convergence and with respect to efficiency. Physiologically relevant extensions of the model presented in this paper-e.g., imagebased geometries, coupling to electrophysiological models, utilizing experimentally measured calcium cell-wise or high-resolution single-cell experiments-are dependent on high-resolution meshes. Anyone investigating such questions would presumably meet similar scalability limitations. Pushing these limits would enable investigation of larger problems, opening up for investigating a wide range of new research questions.

Conclusions
We have introduced a cell-based framework for modeling cardiac tissue, parametrized based on experimental data. This was done by pairing stretching and shearing experiments performed on cardiac tissue samples with corresponding virtual experiments. We observe that the adaption to a geometrical framework, based on relatively simple extensions of existing ideas used for tissue level models, give rise to striking spatial strain and stress values patterns and opens up for numerous new interesting research questions. Utilizing our model, we have been able to differentiate between deformation and stresses in the cells from the matrix surrounding them. Moving to multicellular and finer meshes have, however, proven computationally expensive, as demonstrated by our scaling and convergence experiments.
Our work demonstrates that it is feasible to work with a discretized model which explicitly represents the cells, ranging from a single cell to small collections of cells. Our geometrical approach can further be extended by using more realistic geometries or be coupled with calcium dynamics, which could be used for new studies leading to an improved understanding of cardiac mechanics.

Supporting information
We include movies for the spatial distributions, considering fiber direction stretch and contraction. For all movies we display normal components of and tensor glyphs for the Green-Lagrange strain tensor E (14) (top row) and the Cauchy stress tensor (15) (bottom row).
Movie 1 Fiber direction stretch (FF), single cell. Strain and stress components for a single cell subject to stretching in the fiber direction. Maximum stretching state values, at 10% stretch, are also displayed in Fig. 10.
Movie 2 Fiber direction stretch (FF), tiled cells. Strain and stress components for a 3 × 3 × 3 cells subject to stretching in the fiber direction. Maximum stretching state values, at 10% stretch, are also displayed in Fig. 11.

Movie 3
Contracting cell, single cell. Strain and stress components for a single contracting cell under contraction, for half a cardiac cycle. Maximum contracted state values are also displayed in Fig. 10 and Fig. 13.
Movie 4 Contracting cells, tiled cells. Strain and stress components for 3 × 3 × 3 cardiac cells under contraction, for half a cardiac cycle. Maximum contracted state values are also displayed in Fig. 11.

Appendix: Weak form and boundary conditions
In this appendix, we provide a more rigorous mathematical framework for how the weak form was derived, and how the different boundary conditions for the different deformation modes were imposed. All deformation modes are sketched and explained in Fig. 3.

Derivation of the weak form
Equilibrium of stresses is, in this paper, expressed using the first Piola-Kirchhoff stress tensor, as given in Eq (1). In this work, we employed an active strain formulation and assumed incompressibility imposed by a Lagrangian multiplier. Following the derivations in Ambrosi and Pezzuto (2012) and Holzapfel (2000), we can incorporate these changes in the first Piola-Kirchhoff stress tensor, which then reads where an equilibrium solution still is given by (1). P is dependent on the displacement u and the hydrostatic pressure p, which represent the unknowns.
One can define appropriate function spaces V and Q, defined over Ω , for the displacement and hydrostatic pressure respectively. Starting with the strong form of our problem (1), one can consider test functions v ∈ V and take the integral over Ω i and Ω e separately: and then perform integration by parts in order to obtain We assume the continuity of stresses across Γ as given in (3), which implies that the membrane has no stiffness. If one wanted to incorporate some stiffness here, which physically probably would be meaningful, this could be an equation to change. For simplicity, though, we will keep the formulation as it is, implying that By adding (22) and (23), the integral terms on Γ cancel, and we can combine the integrals to a integral over all of Ω, On the outer boundary Ω , we apply either Dirichlet boundary conditions or zero Neumann conditions, which makes the surface integral over Ω vanish. For the incompressibility condition, given by Eq (5), we can multiply the equation with a test function q ∈ Q to obtain where no decomposition of the subdomain is needed.
By adding these two equations, our weak form becomes and our problem can be formulated as follows: Find a solution (u, p) such that the above equation holds for all v ∈ V and all q ∈ Q. The above equation is solved as a stationary problem, in which time dependent variables are kept constant. These variables either define the boundary displacement, for stretching/shearing experiments, or the active tension (5), for contraction experiments.

Boundary conditions for stretching experiments
We imposed Dirichlet boundary conditions on two of the surfaces, step-wise increasing one of the components, giving a proportional stretch of of the whole domain. If, e.g., if is increased up to 0.15 for the fiber direction stretch (FF), the deformed domain becomes 15% longer in this direction.
(26) ∫ Ω (J − 1)qdX = 0 (27) ∫ Ω P(u, p) ∶ ∇v + q(J(u) − 1)dX = 0 this implies a shortening in the direction perpendicular to the stretching direction. Mathematically the boundary conditions described above can be expressed as for the fiber direction stretch. Similarly, for the sheet direction stretch, we have and for the normal direction stretch, we have For all the surfaces not mentioned above, for each experiment, we impose vanishing Neumann boundary conditions, i.e.,

Boundary conditions for shear experiments
Similarly to the stretching experiments, shear experiments were performed by imposing Dirichlet boundary conditions on two of the surfaces. One of the components was then increased, giving a proportional shear deformation of . If, e.g., was increased up to 0.15, the surface on which we apply non-zero Dirichlet boundary conditions is moved a distance corresponding to 15% of the distance between the two surfaces being fixed.
Mathematically, we have, for the FS experiment for the FN experiment for the SF experiment  For all the surfaces not mentioned above, for each experiment, we impose vanishing Neumann boundary conditions, i.e.,

Boundary conditions for active contraction
For the active contraction experiments, the traction along every surface was set to 0: Due to the Neumann boundary conditions, this formulation does not possess a unique solution. In fact, the formulation has a six-dimensional kernel consisting of rigid deformations, i.e., translation and rotation in a three-dimensional space. We remove these modes by enforcing orthogonality of the solution with respect to the kernel by means of Lagrangian multipliers, as done in, e.g., Kuchta et al. (2016): Let R be the space of rigid motions. Using the finite element formulation, we then seek to find displacement u , pressure p and rigid motion r ∈ R such that and for all v ∈ V , all q ∈ Q and all s ∈ R.
Here r can be interpreted as a Lagrange multiplier enforcing the constraint that the solution u is free of (49) ∫ Ω (P ∶ ∇v + q(J − 1) + r ⋅ v)dx = 0 (50) ∫ Ω u ⋅ sdx = 0 translation or rotation with respect to the original configuration. Equation (50) implies that u is orthogonal to R. In the practical implementation we fix a chosen set of basis vectors, 6 i=1 , of the six-dimensional rigid motion space R and obtain r through its expansion coefficients in the basis, i.e., r = ∑ i c i . Here the requirement (50) is equivalent to imposing ∫ Ω u ⋅ dX = 0 for all 1 ≤ i ≤ 6.