A computational study of right ventricular mechanics in a rat model of pulmonary arterial hypertension

Pulmonary arterial hypertension (PAH) presents a significant challenge to right ventricular (RV) function due to progressive pressure overload, necessitating adaptive remodeling in the form of increased wall thickness, enhanced myocardial contractility and stiffness to maintain cardiac performance. However, the impact of these remodeling mechanisms on RV mechanics in not clearly understood. In addition, there is a lack of quantitative understanding of how each mechanism individually influences RV mechanics. Utilizing experimental data from a rat model of PAH at three distinct time points, we developed biventricular finite element models to investigate how RV stress and strain evolved with PAH progression. The finite element models were fitted to hemodynamic and morphological data to represent different disease stages and used to analyze the impact of RV remodeling as well as the altered RV pressure. Furthermore, we performed a number of theoretical simulation studies with different combinations of morphological and physiological remodeling, to assess and quantify their individual impact on overall RV load and function. Our findings revealed a substantial 4-fold increase in RV stiffness and a transient 2-fold rise in contractility, which returned to baseline by week 12. These changes in RV material properties in addition to the 2-fold increase in wall thickness significantly mitigated the increase in wall stress and strain caused by the progressive increase in RV afterload. Despite the PAH-induced cases showing increased wall stress and strain at end-diastole and end-systole compared to the control, our simulations suggest that without the observed remodeling mechanisms, the increase in stress and strain would have been much more pronounced. Our model analysis also indicated that while changes in the RV's material properties-particularly increased RV stiffness - have a notable effect on its mechanics, the primary compensatory factor limiting the stress and strain increase in the early stages of PAH was the significant increase in wall thickness. These findings underscore the importance of RV remodeling in managing the mechanical burden on the right ventricle due to pressure overload.

Stretch Harmonizes Sarcomere Strain Across the Cardiomyocyte

BACKGROUND

Increasing cardiomyocyte contraction during myocardial stretch serves as the basis for the Frank-Starling mechanism in the heart. However, it remains unclear how this phenomenon occurs regionally within cardiomyocytes, at the level of individual sarcomeres. We investigated sarcomere contractile synchrony and how intersarcomere dynamics contribute to increasing contractility during cell lengthening.

METHODS

Sarcomere strain and Ca2+ were simultaneously recorded in isolated left ventricular cardiomyocytes during 1 Hz field stimulation at 37 °C, at resting length and following stepwise stretch.

RESULTS

We observed that in unstretched rat cardiomyocytes, differential sarcomere deformation occurred during each beat. Specifically, while most sarcomeres shortened during the stimulus, ≈10% to 20% of sarcomeres were stretched or remained stationary. This nonuniform strain was not traced to regional Ca2+ disparities but rather shorter resting lengths and lower force production in systolically stretched sarcomeres. Lengthening of the cell recruited additional shortening sarcomeres, which increased contractile efficiency as less negative, wasted work was performed by stretched sarcomeres. Given the known role of titin in setting sarcomere dimensions, we next hypothesized that modulating titin expression would alter intersarcomere dynamics. Indeed, in cardiomyocytes from mice with titin haploinsufficiency, we observed greater variability in resting sarcomere length, lower recruitment of shortening sarcomeres, and impaired work performance during cell lengthening.

CONCLUSIONS

Graded sarcomere recruitment directs cardiomyocyte work performance, and harmonization of sarcomere strain increases contractility during cell stretch. By setting sarcomere dimensions, titin controls sarcomere recruitment, and its lowered expression in haploinsufficiency mutations impairs cardiomyocyte contractility.

A verified and validated moving domain computational fluid dynamics solver with applications to cardiovascular flows

Computational fluid dynamics (CFD) in combination with patient-specific medical images has been used to correlate flow phenotypes with disease initiation, progression and outcome, in search of a prospective clinical tool. A large number of CFD software packages are available, but are typically based on rigid domains and low-order finite volume methods, and are often implemented in massive low-level C++ libraries. Furthermore, only a handful of solvers have been appropriately verified and validated for their intended use. Our goal was to develop, verify and validate an open-source CFD solver for moving domains, with applications to cardiovascular flows. The solver is an extension of the CFD solver Oasis, which is based on the finite element method and implemented using the FEniCS open source framework. The new solver, named OasisMove, extends Oasis by expressing the Navier-Stokes equations in the arbitrary Lagrangian-Eulerian formulation, which is suitable for handling moving domains. For code verification we used the method of manufactured solutions for a moving 2D vortex problem, and for validation we compared our results against existing high-resolution simulations and laboratory experiments for two moving domain problems of varying complexity. Verification results showed that the L 2 error followed the theoretical convergence rates. The temporal accuracy was second-order, while the spatial accuracy was second- and third-order using ℙ 1 / ℙ 1 and ℙ 2 / ℙ 1 finite elements, respectively. Validation results showed good agreement with existing benchmark results, by reproducing lift and drag coefficients with less than 1% error, and demonstrating the solver's ability to capture vortex patterns in transitional and turbulent-like flow regimes. In conclusion, we have shown that OasisMove is an open-source, accurate and reliable solver for cardiovascular flows in moving domains.

Regional Left Ventricular Fiber Stress Analysis for Cardiac Resynchronization Therapy Response

Cardiac resynchronization therapy (CRT) is an effective treatment for a subgroup of heart failure (HF) patients, but more than 30% of those selected do not improve after CRT implantation. Imperfect pre-procedural criteria for patient selection and optimization are the main causes of the high non-response rate. In this study, we evaluated a novel measure for assessing CRT response. We used a computational modeling framework to calculate the regional stress of the left ventricular wall of seven CRT patients and seven healthy controls. The standard deviation of regional wall stress at the time of mitral valve closure (SD_MVC) was used to quantify dyssynchrony and compared between patients and controls and among the patients. The results show that SD_MVC is significantly lower in controls than patients and correlates with long-term response in patients, based on end-diastolic volume reduction. In contrast to our initial hypothesis, patients with lower SD_MVC respond better to therapy. The patient with the highest SD_MVC was the only non-responder in the patient cohort. The distribution of fiber stress at the beginning of the isovolumetric phase seems to correlate with the degree of response and the use of this measurement could potentially improve selection criteria for CRT implantation. Further studies with a larger cohort of patients are needed to validate these results.

Computational cardiac physiology for new modelers: Origins, foundations, and future

Mathematical models of the cardiovascular system have come a long way since they were first introduced in the early 19th century. Driven by a rapid development of experimental techniques, numerical methods, and computer hardware, detailed models that describe physical scales from the molecular level up to organs and organ systems have been derived and used for physiological research. Mathematical and computational models can be seen as condensed and quantitative formulations of extensive physiological knowledge and are used for formulating and testing hypotheses, interpreting and directing experimental research, and have contributed substantially to our understanding of cardiovascular physiology. However, in spite of the strengths of mathematics to precisely describe complex relationships and the obvious need for the mathematical and computational models to be informed by experimental data, there still exist considerable barriers between experimental and computational physiological research. In this review, we present a historical overview of the development of mathematical and computational models in cardiovascular physiology, including the current state of the art. We further argue why a tighter integration is needed between experimental and computational scientists in physiology, and point out important obstacles and challenges that must be overcome in order to fully realize the synergy of experimental and computational physiological research.

Computational models of ventricular mechanics and adaptation in response to right-ventricular pressure overload

Pulmonary arterial hypertension (PAH) is associated with substantial remodeling of the right ventricle (RV), which may at first be compensatory but at a later stage becomes detrimental to RV function and patient survival. Unlike the left ventricle (LV), the RV remains understudied, and with its thin-walled crescent shape, it is often modeled simply as an appendage of the LV. Furthermore, PAH diagnosis is challenging because it often leaves the LV and systemic circulation largely unaffected. Several treatment strategies such as atrial septostomy, right ventricular assist devices (RVADs) or RV resynchronization therapy have been shown to improve RV function and the quality of life in patients with PAH. However, evidence of their long-term efficacy is limited and lung transplantation is still the most effective and curative treatment option. As such, the clinical need for improved diagnosis and treatment of PAH drives a strong need for increased understanding of drivers and mechanisms of RV growth and remodeling (G&R), and more generally for targeted research into RV mechanics pathology. Computational models stand out as a valuable supplement to experimental research, offering detailed analysis of the drivers and consequences of G&R, as well as a virtual test bench for exploring and refining hypotheses of growth mechanisms. In this review we summarize the current efforts towards understanding RV G&R processes using computational approaches such as reduced-order models, three dimensional (3D) finite element (FE) models, and G&R models. In addition to an overview of the relevant literature of RV computational models, we discuss how the models have contributed to increased scientific understanding and to potential clinical treatment of PAH patients.

Nationwide rollout reveals efficacy of epidemic control through digital contact tracing

Fuelled by epidemiological studies of SARS-CoV-2, contact tracing by mobile phones has been put to use in many countries. Over a year into the pandemic, we lack conclusive evidence on its effectiveness. To address this gap, we used a unique real world contact data set, collected during the rollout of the first Norwegian contact tracing app in the Spring of 2020. Our dataset involves millions of contacts between 12.5% of the adult population, which enabled us to measure the real-world app performance. The technological tracing efficacy was measured at 80%, and we estimated that at least 11.0% of the discovered close contacts could not have been identified by manual contact tracing. Our results also indicated that digital contact tracing can flag individuals with excessive contacts, which can help contain superspreading related outbreaks. The overall effectiveness of digital tracing depends strongly on app uptake, but significant impact can be achieved for moderate uptake numbers. Used as a supplement to manual tracing and other measures, digital tracing can be instrumental in controlling the pandemic. Our findings can thus help informing public health policies in the coming months.

Frequent blood flow restricted training not to failure and to failure induces similar gains in myonuclei and muscle mass

The purpose of the present study was to compare the effects of short-term high-frequency failure vs non-failure blood flow-restricted resistance exercise (BFRRE) on changes in satellite cells (SCs), myonuclei, muscle size, and strength. Seventeen untrained men performed four sets of BFRRE to failure (Failure) with one leg and not to failure (Non-failure; 30-15-15-15 repetitions) with the other leg using knee-extensions at 20% of one repetition maximum (1RM). Fourteen sessions were distributed over two 5-day blocks, separated by a 10-day rest period. Muscle samples obtained before, at mid-training, and 10-day post-intervention (Post10) were analyzed for muscle fiber area (MFA), myonuclei, and SC. Muscle size and echo intensity of m.rectus femoris (RF) and m.vastus lateralis (VL) were measured by ultrasonography, and knee extension strength with 1RM and maximal isometric contraction (MVC) up until Post24. Both protocols increased myonuclear numbers in type-1 (12%-17%) and type-2 fibers (20%-23%), and SC in type-1 (92%-134%) and type-2 fibers (23%-48%) at Post10 (p < 0.05). RF and VL size increased by 5%-10% in both legs at Post10 to Post24, whereas the MFA of type-1 fibers in Failure was decreased at Post10 (-10 ± 16%; p = 0.02). Echo intensity increased by ~20% in both legs during Block1 (p < 0.001) and was ~8 to 11% below baseline at Post24 (p = 0.001-0.002). MVC and 1RM decreased by 5%-10% after Block1, but increased in both legs by 6%-11% at Post24 (p < 0.05). In conclusion, both short-term high-frequency failure and non-failure BFRRE induced increases in SCs, in myonuclei content, muscle size, and strength, concomitant with decreased echo intensity. Intriguingly, the responses were delayed and peaked 10-24 days after the training intervention. Our findings may shed light on the mechanisms involved in resistance exercise-induced overreaching and supercompensation.

Computational Modeling Studies of the Roles of Left Ventricular Geometry, Afterload, and Muscle Contractility on Myocardial Strains in Heart Failure with Preserved Ejection Fraction

Global longitudinal strain and circumferential strain are found to be reduced in HFpEF, which some have interpreted that the global left ventricular (LV) contractility is impaired. This finding is, however, contradicted by a preserved ejection fraction (EF) and confounded by changes in LV geometry and afterload resistance that may also affect the global strains. To reconcile these issues, we used a validated computational framework consisting of a finite element LV model to isolate the effects of HFpEF features in affecting systolic function metrics. Simulations were performed to quantify the effects on myocardial strains due to changes in LV geometry, active tension developed by the tissue, and afterload. We found that only a reduction in myocardial contractility and an increase in afterload can simultaneously reproduce the blood pressures, EF and strains measured in HFpEF patients. This finding suggests that it is likely that the myocardial contractility is reduced in HFpEF patients. Graphical abstract.

A Computational Study of the Effects of Tachycardia-Induced Remodeling on Calcium Wave Propagation in Rabbit Atrial Myocytes

In atrial cardiomyocytes without a well-developed T-tubule system, calcium diffuses from the periphery toward the center creating a centripetal wave pattern. During atrial fibrillation, rapid activation of atrial myocytes induces complex remodeling in diffusion properties that result in failure of calcium to propagate in a fully regenerative manner toward the center; a phenomenon termed “calcium silencing.” This has been observed in rabbit atrial myocytes after exposure to prolonged rapid pacing. Although experimental studies have pointed to possible mechanisms underlying calcium silencing, their individual effects and relative importance remain largely unknown. In this study we used computational modeling of the rabbit atrial cardiomyocyte to query the individual and combined effects of the proposed mechanisms leading to calcium silencing and abnormal calcium wave propagation. We employed a population of models obtained from a newly developed model of the rabbit atrial myocyte with spatial representation of intracellular calcium handling. We selected parameters in the model that represent experimentally observed cellular remodeling which have been implicated in calcium silencing, and scaled their values in the population to match experimental observations. In particular, we changed the maximum conductances of ICaL, INCX, and INaK, RyR open probability, RyR density, Serca2a density, and calcium buffering strength. We incorporated remodeling in a population of 16 models by independently varying parameters that reproduce experimentally observed cellular remodeling, and quantified the resulting alterations in calcium dynamics and wave propagation patterns. The results show a strong effect of ICaL in driving calcium silencing, with INCX, INaK, and RyR density also resulting in calcium silencing in some models. Calcium alternans was observed in some models where INCX and Serca2a density had been changed. Simultaneously incorporating changes in all remodeled parameters resulted in calcium silencing in all models, indicating the predominant role of decreasing ICaL in the population phenotype.

A Novel Computational Model of the Rabbit Atrial Cardiomyocyte With Spatial Calcium Dynamics

Models of cardiac electrophysiology are widely used to supplement experimental results and to provide insight into mechanisms of cardiac function and pathology. The rabbit has been a particularly important animal model for studying mechanisms of atrial pathophysiology and atrial fibrillation, which has motivated the development of models for the rabbit atrial cardiomyocyte electrophysiology. Previously developed models include detailed representations of membrane currents and intracellular ionic concentrations, but these so-called “common-pool” models lack a spatially distributed description of the calcium handling system, which reflects the detailed ultrastructure likely found in cells in vivo. Because of the less well-developed T-tubular system in atrial compared to ventricular cardiomyocytes, spatial gradients in intracellular calcium concentrations may play a more significant role in atrial cardiomyocyte pathophysiology, rendering common-pool models less suitable for investigating underlying electrophysiological mechanisms. In this study, we developed a novel computational model of the rabbit atrial cardiomyocyte incorporating detailed compartmentalization of intracellular calcium dynamics, in addition to a description of membrane currents and intracellular processes. The spatial representation of calcium was based on dividing the intracellular space into eighteen different compartments in the transversal direction, each with separate systems for internal calcium storage and release, and tracking ionic fluxes between compartments in addition to the dynamics driven by membrane currents and calcium release. The model was parameterized employing a population-of-models approach using experimental data from different sources. The parameterization of this novel model resulted in a reduced population of models with inherent variability in calcium dynamics and electrophysiological properties, all of which fall within the range of observed experimental values. As such, the population of models may represent natural variability in cardiomyocyte electrophysiology or inherent uncertainty in the underlying experimental data. The ionic model population was also able to reproduce the U-shaped waveform observed in line-scans of triggered calcium waves in atrial cardiomyocytes, characteristic of the absence of T-tubules, resulting in a centripetal calcium wave due to subcellular calcium diffusion. This novel spatial model of the rabbit atrial cardiomyocyte can be used to integrate experimental findings, offering the potential to enhance our understanding of the pathophysiological role of calcium-handling abnormalities under diseased conditions, such as atrial fibrillation.

Uncertainty quantification and sensitivity analysis of left ventricular function during the full cardiac cycle

Patient-specific computer simulations can be a powerful tool in clinical applications, helping in diagnostics and the development of new treatments. However, its practical use depends on the reliability of the models. The construction of cardiac simulations involves several steps with inherent uncertainties, including model parameters, the generation of personalized geometry and fibre orientation assignment, which are semi-manual processes subject to errors. Thus, it is important to quantify how these uncertainties impact model predictions. The present work performs uncertainty quantification and sensitivity analyses to assess the variability in important quantities of interest (QoI). Clinical quantities are analysed in terms of overall variability and to identify which parameters are the major contributors. The analyses are performed for simulations of the left ventricle function during the entire cardiac cycle. Uncertainties are incorporated in several model parameters, including regional wall thickness, fibre orientation, passive material parameters, active stress and the circulatory model. The results show that the QoI are very sensitive to active stress, wall thickness and fibre direction, where ejection fraction and ventricular torsion are the most impacted outputs. Thus, to improve the precision of models of cardiac mechanics, new methods should be considered to decrease uncertainties associated with geometrical reconstruction, estimation of active stress and of fibre orientation. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Multisite pacing and myocardial scars: a computational study

Cardiac resynchronization therapy (CRT) is a frequently effective treatment modality for dyssynchronous heart failure, however, 30% of patients do not respond, usually due to suboptimal activation of the left ventricle (LV). Multisite pacing (MSP) may increase the response rate, but its effect in the presence of myocardial scars is not fully understood. We use a computational model to study the outcome of MSP in an LV with scars in two different locations and of two different sizes. The LV was stimulated from anterior, posterior and lateral locations individually and in pairs, while a septal stimulation site represented right ventricular (RV) pacing. Intraventricular pressures were measured, and outcomes evaluated in terms of maximum LV pressure gradient (dP/dt_max)- change compared to isolated RV pacing. The best result obtained using various LV pacing locations included a combination of sites remote from scars and the septum. The highest dP/dt_max increase was achieved, regardless of scar size, using MSP with one pacing site located on the LV free wall opposite to the scar and one site opposite to the septum. These in silico modelling results suggest that making placement of pacing electrodes dependent on location of scarring, may alter acute haemodynamics and that such modelling may contribute to future CRT optimization.

Arrhythmogenic Current Generation by Myofilament-Triggered Ca2+ Release and Sarcomere Heterogeneity

Heterogeneous mechanical dyskinesis has been implicated in many arrhythmogenic phenotypes. Strain-dependent perturbations to cardiomyocyte electrophysiology may contribute to this arrhythmogenesis through processes referred to as mechanoelectric feedback. Although the role of stretch-activated ion currents has been investigated using computational models, experimental studies suggest that mechanical strain may also promote arrhythmia by facilitating calcium wave propagation. To investigate whether strain-dependent changes in calcium affinity to the myofilament may promote arrhythmogenic intracellular calcium waves, we modified a mathematical model of rabbit excitation-contraction coupling coupled to a model of myofilament activation and force development. In a one-dimensional compartmental analysis, we bidirectionally coupled 50 sarcomere models in series to model calcium diffusion and stress transfer between adjacent sarcomeres. These considerations enabled the model to capture 1) the effects of mechanical feedback on calcium homeostasis at the sarcomeric level and 2) the combined effects of mechanical and calcium heterogeneities at the cellular level. The results suggest that in conditions of calcium overload, the vulnerable window of stretch-release to trigger suprathreshold delayed afterdepolarizations can be affected by heterogeneity in sarcomere length. Furthermore, stretch and sarcomere heterogeneity may modulate the susceptibility threshold for delayed afterdepolarizations and the aftercontraction wave propagation velocity.

Computational quantification of patient-specific changes in ventricular dynamics associated with pulmonary hypertension

Pulmonary arterial hypertension (PAH) causes an increase in the mechanical loading imposed on the right ventricle (RV) that results in progressive changes to its mechanics and function. Here, we quantify the mechanical changes associated with PAH by assimilating clinical data consisting of reconstructed three-dimensional geometry, pressure, and volume waveforms, as well as regional strains measured in patients with PAH (n = 12) and controls (n = 6) within a computational modeling framework of the ventricles. Modeling parameters reflecting regional passive stiffness and load-independent contractility as indexed by the tissue active tension were optimized so that simulation results matched the measurements. The optimized parameters were compared with clinical metrics to find usable indicators associated with the underlying mechanical changes. Peak contractility of the RV free wall (RVFW) γ_RVFW,max was found to be strongly correlated and had an inverse relationship with the RV and left ventricle (LV) end-diastolic volume ratio (i.e., RVEDV/LVEDV) (RVEDV/LVEDV)+ 0.44, R² = 0.77). Correlation with RV ejection fraction (R² = 0.50) and end-diastolic volume index (R² = 0.40) were comparatively weaker. Patients with with RVEDV/LVEDV > 1.5 had 25% lower γ_RVFW,max (P < 0.05) than that of the control. On average, RVFW passive stiffness progressively increased with the degree of remodeling as indexed by RVEDV/LVEDV. These results suggest a mechanical basis of using RVEDV/LVEDV as a clinical index for delineating disease severity and estimating RVFW contractility in patients with PAH.NEW & NOTEWORTHY This article presents patient-specific data assimilation of a patient cohort and physical description of clinical observations.

Uncertainty in cardiac myofiber orientation and stiffnesses dominate the variability of left ventricle deformation response

Computational cardiac modelling is a mature area of biomedical computing and is currently evolving from a pure research tool to aiding in clinical decision making. Assessing the reliability of computational model predictions is a key factor for clinical use, and uncertainty quantification (UQ) and sensitivity analysis are important parts of such an assessment. In this study, we apply UQ in computational heart mechanics to study uncertainty both in material parameters characterizing global myocardial stiffness and in the local muscle fiber orientation that governs tissue anisotropy. The uncertainty analysis is performed using the polynomial chaos expansion (PCE) method, which is a nonintrusive meta-modeling technique that surrogates the original computational model with a series of orthonormal polynomials over the random input parameter space. In addition, in order to study variability in the muscle fiber architecture, we model the uncertainty in orientation of the fiber field as an approximated random field using a truncated Karhunen-Loéve expansion. The results from the UQ and sensitivity analysis identify clear differences in the impact of various material parameters on global output quantities. Furthermore, our analysis of random field variations in the fiber architecture demonstrate a substantial impact of fiber angle variations on the selected outputs, highlighting the need for accurate assignment of fiber orientation in computational heart mechanics models.

Electromechanical Model to Predict Cardiac Resynchronization Therapy

Cardiac resynchronization therapy (CRT) can substantially improve dyssynchronous heart failure and reduce mortality. However, one-third of the CRT patients derive no measurable benefit from CRT, due to suboptimal placement of the left ventricular (LV) lead. We introduce a pipeline for improved CRT-therapy by creating an electromechanical model using patient-specific geometric parameters allowing individualization of therapy. The model successfully mimics expected changes when variables for tension, stiffness, and conduction are entered. Changing LV pacing site had a notable effect on maximum pressure gradient (dP/dt_max) in the presence of cardiac scarring, causing non-uniform excitation propagation through the LV. Tailoring CRT to the individual requires simulations with patient-specific biventricular meshes including cardiac geometry and conductivity properties.

Computational Modeling of Electrophysiology and Pharmacotherapy of Atrial Fibrillation: Recent Advances and Future Challenges

The pathophysiology of atrial fibrillation (AF) is broad, with components related to the unique and diverse cellular electrophysiology of atrial myocytes, structural complexity, and heterogeneity of atrial tissue, and pronounced disease-associated remodeling of both cells and tissue. A major challenge for rational design of AF therapy, particularly pharmacotherapy, is integrating these multiscale characteristics to identify approaches that are both efficacious and independent of ventricular contraindications. Computational modeling has long been touted as a basis for achieving such integration in a rapid, economical, and scalable manner. However, computational pipelines for AF-specific drug screening are in their infancy, and while the field is progressing quite rapidly, major challenges remain before computational approaches can fill the role of workhorse in rational design of AF pharmacotherapies. In this review, we briefly detail the unique aspects of AF pathophysiology that determine requirements for compounds targeting AF rhythm control, with emphasis on delimiting mechanisms that promote AF triggers from those providing substrate or supporting reentry. We then describe modeling approaches that have been used to assess the outcomes of drugs acting on established AF targets, as well as on novel promising targets including the ultra-rapidly activating delayed rectifier potassium current, the acetylcholine-activated potassium current and the small conductance calcium-activated potassium channel. Finally, we describe how heterogeneity and variability are being incorporated into AF-specific models, and how these approaches are yielding novel insights into the basic physiology of disease, as well as aiding identification of the important molecular players in the complex AF etiology.

Efficient estimation of personalized biventricular mechanical function employing gradient-based optimization

Individually personalized computational models of heart mechanics can be used to estimate important physiological and clinically-relevant quantities that are difficult, if not impossible, to directly measure in the beating heart. Here, we present a novel and efficient framework for creating patient-specific biventricular models using a gradient-based data assimilation method for evaluating regional myocardial contractility and estimating myofiber stress. These simulations can be performed on a regular laptop in less than 2 h and produce excellent fit between measured and simulated volume and strain data through the entire cardiac cycle. By applying the framework using data obtained from 3 healthy human biventricles, we extracted clinically important quantities as well as explored the role of fiber angles on heart function. Our results show that steep fiber angles at the endocardium and epicardium are required to produce simulated motion compatible with measured strain and volume data. We also find that the contraction and subsequent systolic stresses in the right ventricle are significantly lower than that in the left ventricle. Variability of the estimated quantities with respect to both patient data and modeling choices are also found to be low. Because of its high efficiency, this framework may be applicable to modeling of patient specific cardiac mechanics for diagnostic purposes.

Organ-level validation of a cross-bridge cycling descriptor in a left ventricular finite element model: effects of ventricular loading on myocardial strains

Although detailed cell‐based descriptors of cross‐bridge cycling have been applied in finite element (FE) heart models to describe ventricular mechanics, these multiscale models have never been tested rigorously to determine if these descriptors, when scaled up to the organ‐level, are able to reproduce well‐established organ‐level physiological behaviors. To address this void, we here validate a left ventricular (LV) FE model that is driven by a cell‐based cross‐bridge cycling descriptor against key organ‐level heart physiology. The LV FE model was coupled to a closed‐loop lumped parameter circulatory model to simulate different ventricular loading conditions (preload and afterload) and contractilities. We show that our model is able to reproduce a linear end‐systolic pressure volume relationship, a curvilinear end‐diastolic pressure volume relationship and a linear relationship between myocardial oxygen consumption and pressure–volume area. We also show that the validated model can predict realistic LV strain‐time profiles in the longitudinal, circumferential, and radial directions. The predicted strain‐time profiles display key features that are consistent with those measured in humans, such as having similar peak strains, time‐to‐peak‐strain, and a rapid change in strain during atrial contraction at late‐diastole. Our model shows that the myocardial strains are sensitive to not only LV contractility, but also to the LV loading conditions, especially to a change in afterload. This result suggests that caution must be exercised when associating changes in myocardial strain with changes in LV contractility. The methodically validated multiscale model will be used in future studies to understand human heart diseases.

In vivo estimation of elastic heterogeneity in an infarcted human heart

In myocardial infarction, muscle tissue of the heart is damaged as a result of ceased or severely impaired blood flow. Survivors have an increased risk of further complications, possibly leading to heart failure. Material properties play an important role in determining post-infarction outcome. Due to spatial variation in scarring, material properties can be expected to vary throughout the tissue of a heart after an infarction. In this study we propose a data assimilation technique that can efficiently estimate heterogeneous elastic material properties in a personalized model of cardiac mechanics. The proposed data assimilation is tested on a clinical dataset consisting of regional left ventricular strains and in vivo pressures during atrial systole from a human with a myocardial infarction. Good matches to regional strains are obtained, and simulated equi-biaxial tests are carried out to demonstrate regional heterogeneities in stress–strain relationships. A synthetic data test shows a good match of estimated versus ground truth material parameter fields in the presence of no to low levels of noise. This study is the first to apply adjoint-based data assimilation to the important problem of estimating cardiac elastic heterogeneities in 3-D from medical images.

Preconditioned augmented Lagrangian formulation for nearly incompressible cardiac mechanics

Computational modeling of the heart is a subject of substantial medical and scientific interest, which may contribute to increase the understanding of several phenomena associated with cardiac physiological and pathological states. Modeling the mechanics of the heart have led to considerable insights, but it still represents a complex and a demanding computational problem, especially in a strongly coupled electromechanical setting. Passive cardiac tissue is commonly modeled as hyperelastic and is characterized by quasi-incompressible, orthotropic, and nonlinear material behavior. These factors are known to be very challenging for the numerical solution of the model. The near-incompressibility is known to cause numerical issues such as the well-known locking phenomenon and ill-conditioning of the stiffness matrix. In this work, the augmented Lagrangian method is used to handle the nearly incompressible condition. This approach can potentially improve computational performance by reducing the condition number of the stiffness matrix and thereby improving the convergence of iterative solvers. We also improve the performance of iterative solvers by the use of an algebraic multigrid preconditioner. Numerical results of the augmented Lagrangian method combined with a preconditioned iterative solver for a cardiac mechanics benchmark suite are presented to show its improved performance.

High-resolution data assimilation of cardiac mechanics applied to a dyssynchronous ventricle

Computational models of cardiac mechanics, personalized to a patient, offer access to mechanical information above and beyond direct medical imaging. Additionally, such models can be used to optimize and plan therapies in-silico, thereby reducing risks and improving patient outcome. Model personalization has traditionally been achieved by data assimilation, which is the tuning or optimization of model parameters to match patient observations. Current data assimilation procedures for cardiac mechanics are limited in their ability to efficiently handle high-dimensional parameters. This restricts parameter spatial resolution, and thereby the ability of a personalized model to account for heterogeneities that are often present in a diseased or injured heart. In this paper, we address this limitation by proposing an adjoint gradient-based data assimilation method that can efficiently handle high-dimensional parameters. We test this procedure on a synthetic data set and provide a clinical example with a dyssynchronous left ventricle with highly irregular motion. Our results show that the method efficiently handles a high-dimensional optimization parameter and produces an excellent agreement for personalized models to both synthetic and clinical data.

A Centerline-Based Model Morphing Algorithm for Patient-Specific Finite Element Modeling of the Left Ventricle

GOAL

Hexahedral automatic model generation is a recurrent problem in computer vision and computational biomechanics. It may even become a challenging problem when one wants to develop a patient-specific finite element (FE) model of the left ventricle (LV), particularly when only low resolution images are available. In the present study, a fast and efficient algorithm is presented and tested to address such a situation.

METHODS

A template FE hexahedral model was created for an LV geometry using a general electric ultrasound (US) system. A system of centerline was considered for this LV mesh. Then, the nodes located over the endocardial and epicardial surfaces are, respectively, projected from this centerline onto the actual endocardial and epicardial surfaces reconstructed from a patient's US data. Finally, the position of the internal nodes is derived by finding the deformations with minimal elastic energy. This approach was applied to eight patients suffering from congestive heart disease. An FE analysis was performed to derive the stress induced in the LV tissue by diastolic blood pressure on each of them.

RESULTS

Our model morphing algorithm was applied successfully and the obtained meshes showed only marginal mismatches when compared to the corresponding US geometries. The diastolic FE analyses were successfully performed in seven patients to derive the distribution of principal stresses.

CONCLUSION

The original model morphing algorithm is fast and robust with low computational cost.

SIGNIFICANCE

This low-cost model morphing algorithm may be highly beneficial for future patient-specific reduced-order modeling of the LV with potential application to other crucial organs.

A computational framework for testing arrhythmia marker sensitivities to model parameters in functionally calibrated populations of atrial cells

Models of cardiac cell electrophysiology are complex non-linear systems which can be used to gain insight into mechanisms of cardiac dynamics in both healthy and pathological conditions. However, the complexity of cardiac models can make mechanistic insight difficult. Moreover, these are typically fitted to averaged experimental data which do not incorporate the variability in observations. Recently, building populations of models to incorporate inter- and intra-subject variability in simulations has been combined with sensitivity analysis (SA) to uncover novel ionic mechanisms and potentially clarify arrhythmogenic behaviors. We used the Koivumäki human atrial cell model to create two populations, representing normal Sinus Rhythm (nSR) and chronic Atrial Fibrillation (cAF), by varying 22 key model parameters. In each population, 14 biomarkers related to the action potential and dynamic restitution were extracted. Populations were calibrated based on distributions of biomarkers to obtain reasonable physiological behavior, and subjected to SA to quantify correlations between model parameters and pro-arrhythmia markers. The two populations showed distinct behaviors under steady state and dynamic pacing. The nSR population revealed greater variability, and more unstable dynamic restitution, as compared to the cAF population, suggesting that simulated cAF remodeling rendered cells more stable to parameter variation and rate adaptation. SA revealed that the biomarkers depended mainly on five ionic currents, with noted differences in sensitivities to these between nSR and cAF. Also, parameters could be selected to produce a model variant with no alternans and unaltered action potential morphology, highlighting that unstable dynamical behavior may be driven by specific cell parameter settings. These results ultimately suggest that arrhythmia maintenance in cAF may not be due to instability in cell membrane excitability, but rather due to tissue-level effects which promote initiation and maintenance of reentrant arrhythmia.

An integrative appraisal of mechano-electric feedback mechanisms in the heart

Mechanically-induced alterations in cardiac electrophysiology are referred to as mechano-electric feedback (MEF), and play an important role in electrical regulation of cardiac performance. The influence of mechanical stress and strain on electrophysiology has been investigated at all levels, however the role of MEF in arrhythmia remains poorly understood. During the normal contraction of the heart, mechano-sensitive processes are an implicit component of cardiac activity. Under abnormal mechanical events, stretch-activated mechanisms may contribute to local or global changes in electrophysiology (EP). While such mechanisms have been hypothesised to be involved in mechanically-initiated arrhythmias, the details of these mechanisms and their importance remain elusive. We assess the theoretical role of stretch mechanisms using coupled models of cellular electrophysiology and sarcomere contraction dynamics. Using models of single ventricular myocytes, we first investigated the potential MEF contributions of stretch-activated currents (SAC), and stretch-induced myofilament calcium release, to test how strain and fibrosis may alter cellular electrophysiology. For all models investigated, SACs were alone not sufficient to create a pro-arrhythmic perturbation of the action potential with stretch. However, when combined with stretch-induced myofilament calcium release, the action potential could be shortened depending on the timing of the strain. This effect was highly model dependent, with a canine epicardial EP model being the most sensitive. These model results suggest that known mechanisms of mechano-electric coupling in cardiac myocyte may be sufficient to be pro-arrhythmic, but only in combination and under specific strain patterns.

An Evaluation of the Accuracy of Classical Models for Computing the Membrane Potential and Extracellular Potential for Neurons

Two mathematical models are part of the foundation of Computational neurophysiology; (a) the Cable equation is used to compute the membrane potential of neurons, and, (b) volume-conductor theory describes the extracellular potential around neurons. In the standard procedure for computing extracellular potentials, the transmembrane currents are computed by means of (a) and the extracellular potentials are computed using an explicit sum over analytical point-current source solutions as prescribed by volume conductor theory. Both models are extremely useful as they allow huge simplifications of the computational efforts involved in computing extracellular potentials. However, there are more accurate, though computationally very expensive, models available where the potentials inside and outside the neurons are computed simultaneously in a self-consistent scheme. In the present work we explore the accuracy of the classical models (a) and (b) by comparing them to these more accurate schemes. The main assumption of (a) is that the ephaptic current can be ignored in the derivation of the Cable equation. We find, however, for our examples with stylized neurons, that the ephaptic current is comparable in magnitude to other currents involved in the computations, suggesting that it may be significant-at least in parts of the simulation. The magnitude of the error introduced in the membrane potential is several millivolts, and this error also translates into errors in the predicted extracellular potentials. While the error becomes negligible if we assume the extracellular conductivity to be very large, this assumption is, unfortunately, not easy to justify a priori for all situations of interest.

Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology

In cardiac electrophysiology, the propagation of the action potential may be described by a set of reaction-diffusion equations known as the bidomain model. The shape of the solution is determined by a balance of a strong reaction and a relatively weak diffusion, which leads to steep variations in space and time. From a numerical point of view, the sharp spatial gradients may be seen as particularly problematic, because computational grid resolution on the order of 0.1 mm or less is required, yielding considerable computational efforts on human geometries. In this paper, we discuss a number of well-known numerical schemes for the bidomain equation and show how the quality of the solution is affected by the spatial discretization. In particular, we study in detail the effect of discretization on the conduction velocity (CV), which is an important quantity from a physiological point of view. We show that commonly applied finite element techniques tend to overestimate the CV on coarse grids, while it tends to be underestimated by finite difference schemes. Furthermore, the choice of interpolation and discretization scheme for the nonlinear reaction term has a strong impact on the CV. Finally, we exploit the results of the error analysis to propose improved numerical methods, including a stabilized scheme that tends to correct the CV on coarse grids but converges to the correct solution as the grid is refined. Copyright © 2016 John Wiley & Sons, Ltd.

Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour

Models of cardiac mechanics are increasingly used to investigate cardiac physiology. These models are characterized by a high level of complexity, including the particular anisotropic material properties of biological tissue and the actively contracting material. A large number of independent simulation codes have been developed, but a consistent way of verifying the accuracy and replicability of simulations is lacking. To aid in the verification of current and future cardiac mechanics solvers, this study provides three benchmark problems for cardiac mechanics. These benchmark problems test the ability to accurately simulate pressure-type forces that depend on the deformed objects geometry, anisotropic and spatially varying material properties similar to those seen in the left ventricle and active contractile forces. The benchmark was solved by 11 different groups to generate consensus solutions, with typical differences in higher-resolution solutions at approximately 0.5%, and consistent results between linear, quadratic and cubic finite elements as well as different approaches to simulating incompressible materials. Online tools and solutions are made available to allow these tests to be effectively used in verification of future cardiac mechanics software.

An integrated electromechanical-growth heart model for simulating cardiac therapies

An emerging class of models has been developed in recent years to predict cardiac growth and remodeling (G&R). We recently developed a cardiac G&R constitutive model that predicts remodeling in response to elevated hemodynamics loading, and a subsequent reversal of the remodeling process when the loading is reduced. Here, we describe the integration of this G&R model to an existing strongly coupled electromechanical model of the heart. A separation of timescale between growth deformation and elastic deformation was invoked in this integrated electromechanical-growth heart model. To test our model, we applied the G&R scheme to simulate the effects of myocardial infarction in a realistic left ventricular (LV) geometry using the finite element method. We also simulate the effects of a novel therapy that is based on alteration of the infarct mechanical properties. We show that our proposed model is able to predict key features that are consistent with experiments. Specifically, we show that the presence of a non-contractile infarct leads to a dilation of the left ventricle that results in a rightward shift of the pressure volume loop. Our model also predicts that G&R is attenuated by a reduction in LV dilation when the infarct stiffness is increased.

Increased cell membrane capacitance is the dominant mechanism of stretch-dependent conduction slowing in the rabbit heart: a computational study

Volume loading of the cardiac ventricles is known to slow electrical conduction in the rabbit heart, but the mechanisms remain unclear. Previous experimental and modeling studies have investigated some of these mechanisms, including stretch-activated membrane currents, reduced gap junctional conductance, and altered cell membrane capacitance. In order to quantify the relative contributions of these mechanisms, we combined a monomain model of rabbit ventricular electrophysiology with a hyperelastic model of passive ventricular mechanics. First, a simplified geometric model with prescribed homogeneous deformation was used to fit model parameters and characterize individual MEF mechanisms, and showed good qualitative agreement with experimentally measured strain-CV relations. A 3D model of the rabbit left and right ventricles was then compared with experimental measurements from optical electrical mapping studies in the isolated rabbit heart. The model was inflated to an end-diastolic pressure of 30 mmHg, resulting in epicardial strains comparable to those measured in the anterior left ventricular free wall. While the effects of stretch activated channels did alter epicardial conduction velocity, an increase in cellular capacitance was required to explain previously reported experimental results. The new results suggest that for large strains, various mechanisms can combine and produce a biphasic relationship between strain and conduction velocity. However, at the moderate strains generated by high end-diastolic pressure, a stretch-induced increase in myocyte membrane capacitance is the dominant driver of conduction slowing during ventricular volume loading.

Effects of deformation on transmural dispersion of repolarization using in silico models of human left ventricular wedge

Mechanical deformation affects the electrical activity of the heart through multiple feedback loops. The purpose of this work is to study the effect of deformation on transmural dispersion of repolarization and on surface electrograms using an in silico human ventricular wedge. To achieve this purpose, we developed a strongly coupled electromechanical cell model by coupling a human left ventricle electrophysiology model and an active contraction model reparameterized for human cells. This model was then embedded in tissue simulations on the basis of bidomain equations and nonlinear solid mechanics. The coupled model was used to evaluate effects of mechanical deformation on important features of repolarization and electrograms. Our results indicate an increase in the T-wave amplitude of the surface electrograms in simulations that account for the effects of cardiac deformation. This increased T-wave amplitude can be explained by changes to the coupling between neighboring myocytes, also known as electrotonic effect. The thickening of the ventricular wall during repolarization contributes to the decoupling of cells in the transmural direction, enhancing action potential heterogeneity and increasing both transmural repolarization dispersion and T-wave amplitude of surface electrograms. The simulations suggest that a considerable percentage of the T-wave amplitude (15%) may be related to cardiac deformation.

Improved discretisation and linearisation of active tension in strongly coupled cardiac electro-mechanics simulations

Mathematical models of cardiac electro-mechanics typically consist of three tightly coupled parts: systems of ordinary differential equations describing electro-chemical reactions and cross-bridge dynamics in the muscle cells, a system of partial differential equations modelling the propagation of the electrical activation through the tissue and a nonlinear elasticity problem describing the mechanical deformations of the heart muscle. The complexity of the mathematical model motivates numerical methods based on operator splitting, but simple explicit splitting schemes have been shown to give severe stability problems for realistic models of cardiac electro-mechanical coupling. The stability may be improved by adopting semi-implicit schemes, but these give rise to challenges in updating and linearising the active tension. In this paper we present an operator splitting framework for strongly coupled electro-mechanical simulations and discuss alternative strategies for updating and linearising the active stress component. Numerical experiments demonstrate considerable performance increases from an update method based on a generalised Rush-Larsen scheme and a consistent linearisation of active stress based on the first elasticity tensor.

The development of a new computational model for the electromechanics of the human ventricular myocyte

In this work we present a new electromechanical cardiac myocyte model tailored to reproduce the electrical and force generating activities of human ventricular myocytes. The model was created by coupling two existing models: the ten Tusscher electrophysiology model and the Rice myofilament mechanics model. The parameters of the new model were adjusted in order to replicate the available experimental data for human myocytes. The main challenges in this work were the strong feedbacks between the models, the high non-linearity of the models and mainly the lack of human data to make the adjustments.

A computationally efficient method for determining the size and location of myocardial ischemia

The purpose of this paper is to introduce a new method for solving the inverse problem of locating ischemic regions in the heart. The electrical activity in the human heart is modeled by the bidomain equations, which can be expanded to compute the potentials on the body surface. The associated inverse problem is to use ECG recordings to gain information about ischemias. We propose an algorithm for doing this, combining the level set method with a simpler minimization problem. Instead of trying to determine the shape, as in the level set method, we simply make the approximation that the ischemia is spherical. The effects of ischemia on the electrical attributes of heart tissue are examined. The new method is tested with computer simulations on synthetic body surface potential maps (BSPMs) in a realistic geometry, using realistic values for the parameters. It is shown to be, in some respects, superior to the level set approach and may be a step toward a practical algorithm useful in medical diagnostics.

Towards a computational method for imaging the extracellular potassium concentration during regional ischemia

We investigate the possibility of using body surface potential maps to image the extracellular potassium concentration during regional ischemia. The problem is formulated as an inverse problem based on a linear approximation of the bidomain model, where we minimize the difference between the results of the model and observations of body surface potentials. The minimization problem is solved by a one-shot technique, where the original PDE system, an adjoint problem, and the relation describing the minimum, are solved simultaneously. This formulation of the problem requires the solution of a 5 x 5 system of linear partial differential equations. The performance of the model is investigated by performing tests based on synthetic data. We find that the model will in many cases detect the correct position and approximate size of the ischemic regions, while some cases are more difficult to locate. It is observed that a simple post-processing of the results produces images that are qualitatively very similar to the true solution.

Numerical solution of the bidomain equations

Knowledge of cardiac electrophysiology is efficiently formulated in terms of mathematical models. However, most of these models are very complex and thus defeat direct mathematical reasoning founded on classical and analytical considerations. This is particularly so for the celebrated bidomain model that was developed almost 40 years ago for the concurrent analysis of extra- and intracellular electrical activity. Numerical simulations based on this model represent an indispensable tool for studying electrophysiology. However, complex mathematical models, steep gradients in the solutions and complicated geometries lead to extremely challenging computational problems. The greatest achievement in scientific computing over the past 50 years has been to enable the solving of linear systems of algebraic equations that arise from discretizations of partial differential equations in an optimal manner, i.e. such that the central processing unit (CPU) effort increases linearly with the number of computational nodes. Over the past decade, such optimal methods have been introduced in the simulation of electrophysiology. This development, together with the development of affordable parallel computers, has enabled the solution of the bidomain model combined with accurate cellular models, on geometries resembling a human heart. However, in spite of recent progress, the full potential of modern computational methods has yet to be exploited for the solution of the bidomain model. This paper reviews the development of numerical methods for solving the bidomain model. However, the field is huge and we thus restrict our focus to developments that have been made since the year 2000.

A comparison of non-standard solvers for ODEs describing cellular reactions in the heart

Mathematical models for the electrical activity in cardiac cells are normally formulated as systems of ordinary differential equations (ODEs). The equations are nonlinear and describe processes occurring on a wide range of time scales. Under normal accuracy requirements, this makes the systems stiff and therefore challenging to solve numerically. As standard implicit solvers are difficult to implement, explicit solvers such as the forward Euler method are commonly used, despite their poor efficiency. Non-standard formulations of the forward Euler method, derived from the analytical solution of linear ODEs, can give significantly improved performance while maintaining simplicity of implementation. In this paper we study the performance of three non-standard methods on two different cell models with comparable complexity but very different stiffness characteristics.

An unconditionally stable numerical method for the Luo–Rudy 1 model used in simulations of defibrillation

Numerical simulations of defibrillation using the Bidomain model coupled to a model of membrane kinetics represent a serious numerical challenge. This is because very high voltages close to defibrillation electrodes demand that extreme time step restrictions be placed on standard numerical schemes, e.g. the forward Euler scheme. A common solution to this problem is to modify the cell model by simple if-tests applied to several equations and rate functions. These changes are motivated by numerical problems rather than physiology, and should therefore be avoided whenever possible. The purpose of this paper is to present a numerical scheme that handles the original model without modifications and which is unconditionally stable for the Luo-Rudy phase 1 model. This also shows that the cell model is mathematically well-behaved, even in the presence of very high voltages. Our theoretical results are illustrated by numerical computations.

A linear system of partial differential equations modeling the resting potential of a heart with regional ischemia

Ischemic ST-segment shift has been modeled using scalar stationary approximations of the bidomain model. Here, we propose an alternative simplification of the bidomain equations: a linear system modeling the resting potential, to be used in determining ischemic TP shift. Results of 2D and 3D simulations show that the linear system model is much more accurate than the scalar model. This improved accuracy is important if the model is to be used for solving the inverse problem of determining the size and location of an ischemic region. Furthermore, the model can provide insight into how the resting potential depends on the variations in the extracellular potassium concentration that characterize ischemic regions.

On the Computational Complexity of the Bidomain and the Monodomain Models of Electrophysiology

The bidomain model, coupled with accurate models of cell membrane kinetics, is generally believed to provide a reasonable basis for numerical simulations of cardiac electrophysiology. Because of changes occurring in very short time intervals and over small spatial domains, discretized versions of these models must be solved on fine computational grids, and small time-steps must be applied. This leads to huge computational challenges that have been addressed by several authors. One popular way of reducing the CPU demands is to approximate the bidomain model by the monodomain model, and thus reducing a two by two set of partial differential equations to one scalar partial differential equation; both of which are coupled to a set of ordinary differential equations modeling the cell membrane kinetics. A reduction in CPU time of two orders of magnitude has been reported. It is the purpose of the present paper to provide arguments that such a reduction is not present when order-optimal numerical methods are applied. Theoretical considerations and numerical experiments indicate that the reduction factor of the CPU requirements from bidomain to monodomain computations, using order-optimal methods, typically is about 10 for simple cell models and less than two for more complex cell models.

A mixed finite element formulation for a non-linear, transversely isotropic material model for the cardiac tissue

In this paper we present a mixed finite element method for modeling the passive properties of the myocardium. The passive properties are described by a non-linear, transversely isotropic, hyperelastic material model, and the myocardium is assumed to be almost incompressible. Single-field, pure displacement-based formulations are known to cause numerical difficulties when applied to incompressible or slightly compressible material cases. This paper presents an alternative approach in the form of a mixed formulation, where a separately interpolated pressure field is introduced as a primary unknown in addition to the displacement field. Moreover, a constraint term is included in the formulation to enforce (almost) incompressibility. Numerical results presented in the paper demonstrate the difficulties related to employing a pure displacement-based method, applying a set of physically relevant material parameter values for the cardiac tissue. The same problems are not experienced for the proposed mixed method. We show that the mixed formulation provides reasonable numerical results for compressible as well as nearly incompressible cases, also in situations of large fiber stretches. There is good agreement between the numerical results and the underlying analytical models.