Skip to main content
Advertisement
  • Loading metrics

Cell to whole organ global sensitivity analysis on a four-chamber heart electromechanics model using Gaussian processes emulators

  • Marina Strocchi ,

    Roles Conceptualization, Data curation, Formal analysis, Methodology, Software, Visualization, Writing – original draft

    marina.strocchi@kcl.ac.uk

    Affiliations School of Biomedical Engineering and Imaging Sciences, King’s College London, London, United Kingdom, National Heart and Lung Institute, Imperial College London, London, United Kingdom

  • Stefano Longobardi,

    Roles Software, Writing – review & editing

    Affiliation School of Biomedical Engineering and Imaging Sciences, King’s College London, London, United Kingdom

  • Christoph M. Augustin,

    Roles Software, Writing – review & editing

    Affiliations Medical University of Graz, Graz, Austria, BioTechMed-Graz, Graz, Austria

  • Matthias A. F. Gsell,

    Roles Software, Writing – review & editing

    Affiliation Medical University of Graz, Graz, Austria

  • Argyrios Petras,

    Roles Software, Writing – review & editing

    Affiliation Johann Radon Institute for Computational and Applied Mathematics (RICAM), Linz, Austria

  • Christopher A. Rinaldi,

    Roles Data curation, Writing – review & editing

    Affiliations School of Biomedical Engineering and Imaging Sciences, King’s College London, London, United Kingdom, Guy’s and St Thomas’ NHS Foundation Trust, London, United Kingdom

  • Edward J. Vigmond,

    Roles Software, Writing – review & editing

    Affiliations University of Bordeaux, CNRS, Bordeaux, Talence, France, IHU Liryc, Bordeaux, Talence, France

  • Gernot Plank,

    Roles Software, Writing – review & editing

    Affiliations Medical University of Graz, Graz, Austria, BioTechMed-Graz, Graz, Austria

  • Chris J. Oates,

    Roles Conceptualization, Methodology, Supervision, Writing – review & editing

    Affiliation Newcastle University, Newcastle upon Tyne, United Kingdom

  • Richard D. Wilkinson,

    Roles Conceptualization, Methodology, Supervision, Writing – review & editing

    Affiliation University of Nottingham, Nottingham, United Kingdom

  • Steven A. Niederer

    Roles Conceptualization, Funding acquisition, Methodology, Resources, Supervision, Writing – review & editing

    Affiliations School of Biomedical Engineering and Imaging Sciences, King’s College London, London, United Kingdom, National Heart and Lung Institute, Imperial College London, London, United Kingdom, Alan Turing Institute, London, United Kingdom

Abstract

Cardiac pump function arises from a series of highly orchestrated events across multiple scales. Computational electromechanics can encode these events in physics-constrained models. However, the large number of parameters in these models has made the systematic study of the link between cellular, tissue, and organ scale parameters to whole heart physiology challenging. A patient-specific anatomical heart model, or digital twin, was created. Cellular ionic dynamics and contraction were simulated with the Courtemanche-Land and the ToR-ORd-Land models for the atria and the ventricles, respectively. Whole heart contraction was coupled with the circulatory system, simulated with CircAdapt, while accounting for the effect of the pericardium on cardiac motion. The four-chamber electromechanics framework resulted in 117 parameters of interest. The model was broken into five hierarchical sub-models: tissue electrophysiology, ToR-ORd-Land model, Courtemanche-Land model, passive mechanics and CircAdapt. For each sub-model, we trained Gaussian processes emulators (GPEs) that were then used to perform a global sensitivity analysis (GSA) to retain parameters explaining 90% of the total sensitivity for subsequent analysis. We identified 45 out of 117 parameters that were important for whole heart function. We performed a GSA over these 45 parameters and identified the systemic and pulmonary peripheral resistance as being critical parameters for a wide range of volumetric and hemodynamic cardiac indexes across all four chambers. We have shown that GPEs provide a robust method for mapping between cellular properties and clinical measurements. This could be applied to identify parameters that can be calibrated in patient-specific models or digital twins, and to link cellular function to clinical indexes.

Author summary

Cardiac function relies on complex links between the single cell and the whole organ. Digital twins or patient-specific models, e.g. computer models that replicate a patient’s heart, can help understanding these links in healthy or diseased states, and improving cardiac patient care. To build a patient-specific model, first we need to quantify which model parameters affect model outputs, to discard those that have little effect and to understand input-output interactions. This normally requires a lot of expensive model evaluations, making this type of analysis very challenging. We used Gaussian processes emulators (GPEs) to reduce the computational costs of our model. The heart simulator we approximated had 117 initial parameters, and was able to simulate whole heart electrical excitation and contraction, cellular dynamics, as well as the circulatory system and the interaction of the heart with the pericardium. Thanks to the GPEs, we were able to identify the most important 45 parameters at a feasible cost, and to study their effect on a wide range of clinically-measured biomarkers for cardiac function. Our analysis provides a comprehensive assay of how cellular function can impact the whole heart, and can be used to investigate a wide range of cardiac pathologies and treatment.

Introduction

Heart function is a multi-scale process, going from sub-cellular mechanisms initiating single cell excitation and contraction, up to the whole organ and body. The heartbeat starts at the right atrium (RA) at the sino-atrial node, where the cells are able to excite spontaneously. The activation wave then travels to the left atrium (LA) and ventricles through the atrioventricular pathways and the His–Purkinje system, causing every myocyte to undergo the action potential, a rapid sequence of changes in the transmembrane potential. During the action potential, calcium ions are released into the cytosol, where they bind to the Troponin C located on the thin filament. This induces a conformational change in the thin filament which exposes the actin binding sites to the myosyn heads on the thick filament. The myosin heads then bind to the actin and pull the thin filament towards the centre of the sarcomere, leading to sarcomere shortening and ultimately myocyte contraction. Thanks to the structured myocyte orientation within the cardiac muscle, this results into a coordinated contraction of the whole heart, which then pumps oxygenated blood across the aortic valve and into the whole circulation.

Cardiac electromechanics models are increasingly used to study heart function in healthy and diseased states [1]. Initially, heart models focused on the left ventricle (LV) due to otherwise prohibitive computational costs. Recent developments in code and high performance computers have made it possible to simulate all four chambers, offering a more accurate representation of the interaction between the atria and the ventricles and a more physiological systolic motion [29]. Nevertheless, running large numbers of simulations remains a challenge due to numerical instabilities of highly non-linear mechanics, and due to high computational costs. This, combined with the large numbers of parameters of the model, makes performing global sensitivity analysis (GSA) and parameter inference very challenging. In this context, tools from statistics and machine learning can be used to approximate complex and expensive three-dimensional models [10, 11]. Fast model evaluations offered by these tools allow, firstly, to run a GSA to increase model credibility [12] and to identify important model parameters and, secondly, to infer the values of these parameters to fit clinical data at a fraction of the computational cost. This can finally enable the construction of patient-specific four-chamber models, or digital twins, that replicate how a specific subject progresses and responds to therapy. The ability of multi-scale simulations to link cellular processes to whole organ function has the potential to facilitate therapy decision-making, improve patient stratification, and provide novel mechanistic insight into treatment of a wide range of cardiac pathologies.

In order to perform parameter inference and build a digital twin of a human heart, we need to quantify the effect of the input parameters on the model outputs of interest. This not only allows us to exclude the parameters that have little to no effect on model outputs and can therefore be excluded from the analysis, but it also provides information about the extent of the model input-output interactions. In this paper, we used Gaussian processes emulators (GPE) to perform a GSA on a four-chamber electromechanics model. Our simulator accounts for the atrial and ventricular action potential, calcium transient, cross-bridge kinetics, whole heart electrical excitation and contraction, the effect of the pericardium and the coupling with the circulatory system. By applying GPEs and GSA on different components of the electromechanics simulation framework, we were able to reduce the number of parameters from 117 to 45. We used GPEs and history matching (HM) on the tissue electrophysiology and electromechanics cell simulators to identify the areas of the parameter space where the total ventricular and atrial activation times, and cellular calcium and active tension transients for atria and ventricles were physiological. Finally, we performed a GSA on the fully coupled framework with the remaining 45 parameters to quantify the effect of the model parameters on atrial and the ventricular function.

Methods

Below, we summarize the clinical data and the four-chamber electromechanics framework used in this study. Tables 1, 2 and 3 list all simulator parameters we initially considered. We measured 17 of these parameters directly from clinical data (for example the heart rate and the valve areas). We then performed a GSA on different simulator sub-components to identify unimportant parameters and exclude these from further analysis.

thumbnail
Table 1. Simulator parameters.

The first, second and third columns show the parameter name, range or value, and its meaning. The last column provides the original paper the symbol refers to. The blue rows indicate parameters that were estimated from clinical data, while the gray rows represent parameters the model outputs were insensitive to. HM in the second column indicates the parameters that were constrained by a history matching procedure on a sub-model (see S2, S3 and S4 Files). Abbreviations/symbols: CV = conduction velocity, FEC = fast endocardial conduction, Na+ = sodium, K+ = potassium, Ca2+ = calcium, Cl- = chloride, SR = sarcoplasmic reticulum, DS = diadic space, HM = history matching.

https://doi.org/10.1371/journal.pcbi.1011257.t001

thumbnail
Table 2. Simulator parameters (continued).

Abbreviations/symbols: Ca2+ = calcium, CaTRPN = calcium/troponin complex, LV = left ventricle, RV = right ventricle, fraction of unbound (U), and weakly (W) or strongly (S) bound binding sites.

https://doi.org/10.1371/journal.pcbi.1011257.t002

Clinical data

The long-term goal of this work is to construct a digital twin of a specific subject. Therefore, we focused this study on one patient. The clinical data were gathered from a 78 yo female heart failure (HF) patient with atrial fibrillation (AF) selected for cardiac resynchronisation therapy device upgrade from RV pacing to biventricular pacing. Right ventricular (RV) pacing was therefore the baseline rhythm, with a QRS duration of 200 ms measured from a 12-lead electrocardiogram (ECG). The patient underwent ECG-gated CT prior to the upgrade procedure, providing ten CT frames over a cardiac cycle. During the upgrade procedure, the LV pressure was invasively recorded through a pressure wire. After the ectopic baseline beats were removed, 30 beats were left to characterise the patient’s LV haemodynamics. The LV end-diastolic pressure (EDP) and peak in pressure were 2.8±3.2 mmHg and 124.2±7.4 mmHg, while the basic cycle length and the LV pressure systolic duration were 854±9 ms and 467±15 ms, respectively.

Four-chamber electromechanics framework

Four-chamber heart geometry.

The four-chamber heart geometry was generated from the end-diastolic computed tomography (CT). The pipeline to segment and generate 1mm linear tetrahedral four-chamber meshes was described previously [2, 3]. The atria were refined with the resample algorithm from meshtool [13] to have at least 3 elements across the wall thickness to reduce locking effects. The ventricles were assigned with a myofibre orientation using Bayer’s rule-base method [14] (Fig 1, bottom right), with the fibre and sheet angles at the endocardium and epicardium set to be +60° and -60° [4], and -65° and +25° [14], respectively. Atrial myofibre orientation was assigned by computing universal atrial coordinates on the atria and mapping an ex-vivo diffusion tensor MRI dataset onto the endocardial and the epicardial surfaces (Fig 1, top right) [15, 16]. The transmural fibre orientation was then set to be the endocardial and the epicardial orientation for elements below and above 50% of the wall thickness, respectively.

thumbnail
Fig 1. Four-chamber heart geometry.

Left: patient-specific four-chamber heart mesh. Right: refines atria, atrial and ventricular transmural myofiber orientation from endocardium (red) to epicardium (blue).

https://doi.org/10.1371/journal.pcbi.1011257.g001

Electromechanics simulation framework.

The electrical activation of the heart was simulated with a reaction-Eikonal model without diffusion [17]. The Eikonal model in Eq (1) solves for local activation times ta(x) at node location x given V(x) containing the squared local conduction velocities (CV) in the fibres, sheet and normal to sheet directions, and sites of initial activation Γ, which activate at a prescribed time t0: (1)

Atrial and ventricular myocardium were represented as transversely isotropic conductive regions and were assigned CV in the fibre direction (CVf,V and CVf,A) and an anisotropy ratio (kft,V and kft,A), respectively. The remaining regions were passive. To represent fast endocardial activation due to the His–Purkinje system, we defined a one-element thick endocardial layer extending up to 70% in the apico-basal direction in the ventricles [18, 19] with faster conduction velocity (CV) compared to the rest of ventricular myocardium of a factor kFEC (Fig 2, right). Equivalently, to account for the Bachmann bundle, we defined a region between the LA and the RA with fast CV compared to the rest of the atrial myocardium of a factor kBB (Fig 2, left) [16]. To fully control the atrioventricular (AV) delay, we defined a passive region along the AV plane to insulate the atria from the ventricles. Atrial activation was then initiated at the location of the RA lead, while ventricular activation was initiated at the RV lead location with a delay defined by the AV delay, included as a free parameter in the simulator (AVdelay). The RA and RV lead locations were selected by segmenting the pacemaker leads from the CT image by thresholding the image intensity.

thumbnail
Fig 2. Electrophysiology simulations.

Left: atria with the region representing the Bachman bundle (red) and the atrial activation site (light-blue circle). Right: ventricles with the fast endocardial conduction layer (blue) and the ventricular activation site (orange circle).

https://doi.org/10.1371/journal.pcbi.1011257.g002

The action potential of ventricular and atrial myocytes was simulated with the ToR-ORd model with dynamic intracellular cloride [20, 21] and the Courtemanche model [22], respectively. At each point in the domain x, a foot current inducing the initial increase in the transmembrane potential Vm is imposed as stimulus to locally activate the cell membrane, at a local activation time ta(x) computed with the Eikonal model [17]. The intracellular calcium transient computed by the ionic model was then provided as input to the Land contraction model [23] to compute the active tension transient in the atria and ventricles. For simplicity, we assumed that active contraction occurred only in the fibre direction.

In the cell electrophysiology models, we selected model parameters that have a biophysical interpretation but can be, at the same time, considered uncertain [24]. Ion channel conductances, pumps and exchangers leading to transmembrane currents depend on protein expression. Therefore, these can be assumed to be different between individual cells and, since they relate to channel biophysics, they can be considered to be epistemic uncertainty. On the other hand, the natural variability in other parameters such as gating kinetics of ion channels, which we did not consider, can be attributed to aleatoric uncertainty, which is irreducible [25]. Based on this, we initially included all ion channel conductances and all calcium kinetics parameters (Table 1). All 17 parameters for the Land model were included in the initial analysis for both atria and ventricles. Tables 1 and 2 summarise the parameters we considered and their meaning.

Prior to the whole organ simulations, the ToR-ORd-Land and the Courtemanche-Land cell models were ran for 500 beats at a basic cycle length of 854.0 ms, consistent with the patient’s heart rate, to bring the models to a steady state. The cell models were then initialised at a steady state in the whole heart simulations.

Passive material properties of atria and ventricular myocardium were represented with a transversely isotropic Guccione model [26]: (2) where J is the determinant of the deformation gradient, E represents the Cauchy-Green strain tensor and f, s and n are the fibre, sheet and normal to sheet directions. Parameters a, bf, bft and bt are the stiffness parameters, and κ = 1000 kPa is the bulk modulus, penalising volume change and therefore enforcing near incompressibility [27, 28]. All material stiffness parameters for the atria and the ventricles were initially included as free simulator parameters. Passive material properties of all other tissues were represented with a Neo-Hookean model, with the parameter stiffness set according to previous studies [2, 3].

Mechanics boundary conditions.

We represented the effect of the pericardium on the heart with normal springs, as described in [2, 29]. The spring stiffness kperi was included as a free simulator parameter. This value was scaled on the ventricles according to a map derived from motion data [3], to constrain the motion of the apex but not the base, allowing for physiological AV plane downward displacement during ventricular systole. A similar analysis on the atria, described in [30], showed that the roof of the atria moved the least, while the regions around the AV plane moved the most as they were stretched down by the contracting ventricles. We therefore defined a scaling map on the atria to include this constraint in the model, by assigning maximum penalty to the roof of the atria and zero penalty towards the AV plane (Fig 3A). In addition, we applied omni-directional springs to the right inferior and superior pulmonary veins and at the superior vena cava rings. The stiffness of these springs was fixed to 1.0 kPa/μm and was not considered as a free parameter as this was only imposed to constrain the motion during the unloading of the mesh, performed in the absence of the pericardium.

thumbnail
Fig 3. Mechanics boundary conditions.

A Penalty map scaling the normal spring stiffness for the effect of the pericardium. B Afterload and preload boundary conditions represented with CircAdapt. Symbols and abbreviations: p = pressure, R = resistance, q = flow across a valve, LV = left ventricle, RV = right ventricle =, LA = left atrium, RA = right atrium, Ao = aorta, Pa = pulmonary artery, Ve = veins, PVe = pulmonary veins, sys = systemic, pulm = pulmonary, MV = mitral valve, TV = tricuspid valve, AV = aortic valve, PV = pulmonary valve.

https://doi.org/10.1371/journal.pcbi.1011257.g003

The three-dimensional (3D) four-chamber electromechanics model was coupled with CircAdapt [31, 32] (Fig 3B), a closed-loop zero-dimensional (0D) model representing the following components of the circulatory system: aorta, pulmonary artery, veins, systemic and pulmonary peripheral resistances, the four cardiac valves (aortic, pulmonary, mitral and tricuspid) and flows across the pulmonary veins into the LA and across the systemic veins into the RA. The monolithic 3D solid–0D fluid coupling method is described in detail by Augustin et al. [29]. Briefly, the pressures of the LA, LV, RA and RV were included as additional unknowns to the monolithic scheme, and the following equations are added to the mechanics equilibrium equations: where V3D and V0D are the volumes of the cavity computed from the deforming 3D mesh and predicted by the 0D model, respectively, t is the time and u is the displacement field. Table 3 lists all CircAdapt parameters that were included in the initial analysis. In this study, CircAdapt adaptation rules were not applied.

The ventricles of the end-diastolic mesh were unloaded from an end-diastolic LV and RV pressure, while the atria were not unloaded, under the assumption that the active tension in the atrial myocardium balances the pressure [23]. The measured mean LV end-diastolic pressure was only 2.8 mmHg, which was not enough to achieve physiological end-diastolic strains. This may reflect an offset error or drift of the pressure catheter measurements. To account for this potential artifact, we introduced an additional free offset parameter called EDPshift,LV to shift the LV EDP during unloading. This not only allowed the simulations to achieve higher end-diastolic strains, but it also accounted for uncertainty on the available clinical measurements. The RV end-diastolic pressure was not available, so we added it as a free simulator parameter (EDPunload,RV). The unloading was performed with a backward displacement method [33]. During the unloading, we did not apply boundary conditions for the effect of the pericardium. Then, prior to the start of the 3D-0D coupled simulation, we reloaded the ventricles to retrieve the end-diastolic mesh while the atrial pressure was initialised at 0 mmHg. The simulation was then started at end-diastole and the pericardium boundary conditions were turned on. To minimise the effect of these initial conditions, we ran all simulations for 5 beats, to reach a near-to-steady-state behaviour.

Patient-specific parameters.

We used the available clinical data to estimate 17 simulator parameters, therefore reducing the computational cost of the final GSA (light-blue rows, Tables 1 and 3). First, a CT motion tracking algorithm [34] was applied to track the motion of the patient’s heart, and the displacement field was used to deform the patient-specific end-diastolic mesh and derive an LV volume transient. This provided the patient’s LV end-diastolic and end-systolic volumes, and therefore an LV stroke volume of 70 mL.

The 30 baseline LV pressure beats were used to compute an average basic cycle length (BCL) of 854 ms, corresponding to a heart rate (HR) of 70 bpm. The ratio between the SV and the BCL provided the systemic flow qref = 82 mL/s in the CircAdapt model. The wall volumes of the LV and RV free walls, septum, LA and RA were computed from the tetrahedral mesh. The mitral valve open orifice area was computed by selecting points along the mitral valve leaflets from the end-diastolic CT, a sphere was fitted to the selected points and the radius was measured assuming a circular orifice shape (consistently with CircAdapt). The aortic valve open orifice area was measured with the same approach, by selecting points on the aortic valve leaflets on the end-systolic CT frame. The leaflets of the tricuspid and pulmonary valves were not visible on the CT, so their open orifice area was assumed to be the same as the mitral and the aortic valve, respectively, again consistent with the CircAdapt model. The systemic and pulmonary orifices areas were computed as the average of the cross-sectional area of the superior and inferior vena cava, and of all pulmonary veins, respectively, all computed from the end-diastolic CT. The reference area for the aorta, the pulmonary artery, the systemic veins and the pulmonary veins tubes were set to be the same as the corresponding adjacent valve. The aortic and the pulmonary artery wall area were then computed assuming a circular cross-section (consistent with CircAdapt) and a wall thickness of 2 mm [35]. The systemic and the pulmonary veins were assigned a wall thickness of 1 mm, because the veins are thinner due to the presence of less smooth muscles compared to arteries [36] as they operate at lower pressures. The blue rows in Table 3 show the values of the parameters that were computed from the patient-specific data and mesh.

Numerical methods.

As shown in Tables 13, the final analysis included 45 parameters. In order to achieve accurate GPE accuracy, we selected N = 500 samples for GPE training. Due to simulator complexity, this requires a lot of computational resources. We therefore selected the numerical tolerances in the mechanics simulations to reduce as much as possible the computational cost of the electromechanical simulations.

We assumed that 5 beats were enough to achieve a near steady state solution. The first three beats were ran by setting the number of Newton iterations to 1. As shown by Augustin et al. [29], this approach brings the simulation closer to a steady state before solving non-linear mechanics more accurately with more Newton iterations. The last two beats were ran by setting the number of Newton iterations to 2 to have a better approximation of the stretch rate for the cell model. We also increased the tolerance for the solution of the linearised system to 10-4 for all beats. In S1 File, we show that these settings have a limited effect on the pressure-volume dynamics simulated by the model, with differences in the pressure and volume features always below 3%, while allowing up to 3 times speedup in the simulation time. Finally, the time step for the ionic models and the mechanics were set to 0.02 ms and 1.0 ms, respectively.

All simulations were run with the cardiac arrhythmia research package (CARP) [37] on a supercomputer on 512 cores. Details about the software implementation is provided in [17, 29].

Emulators, global sensitivity analysis and history matching

GPEs, GSA and HM methods used in this study were based on Longobardi et al. [11], using code available on GitHub (https://github.com/stelong/Historia.git). Briefly, a GPE provides a statistical model for how a scalar simulator output f(x) varies as a function of D input parameters x = (x1, …, xD). Under our GPE, the simulator outputs are modelled as being jointly Gaussian with prior mean and covariance (3) (4) where k(x, x′) is a positive-definite kernel, and β0, …, βD are the degrees of freedom of the GPE, which need to be estimated. In this case, we chose the exponentiated quadratic kernel, also known as squared exponential or the Gaussian kernel, because it offers a parsimonious model even in situations where the data is limited and noisy [43]. In S8 File, we show that GPE accuracy was not affected when using Matérn rather than an exponentiated quadratic kernel. We trained one GPE for each scalar simulator output, using a training set consisting of simulator outputs computed for a variety of different input parameter configurations as described in [11].

GPE accuracy was evaluated by performing a 5-fold cross-validation. For each split, we computed the coefficient of determination R2 and the independent standard error (ISE) between a vector of outputs Y and the posterior predictions of the GPE with expected value and standard deviation and as: (5) (6) where RSS and TSS are the sum of squared residuals and the total sum of squares, respectively. R2 provides an error on the point-wise estimate of the GPE, with values close to 1 indicating a low error between predictions and observations. The ISE evaluates the uncertainty quantification of the GPEs, and represents the number of data points that fall within two standard deviations away from the average prediction. As it is presented as a percentage over all evaluated points, an ISE close to 100% indicates that, for most points, the distance between the GPE prediction and the corresponding observation is within the GPE uncertainty. Once the GPEs were validated, we trained an additional GPE for each output using all available simulations that was then used for GSA and HM.

Using the GPEs to predict the simulator outputs, we performed a Sobol variance based GSA to quantify the importance of each parameter [44]. We computed the total effects ST using the Saltelli method with the SALib Python library. A Saltelli sampling [44, 45] was used to generate a convenient structure of the samples to efficiently approximate the sensitivity indices as described in S7 File. When computing the sensitivity indices of each output, we sampled the posterior distribution of each GPE with N = 1000 samples, computed the total effect for each sample and taken the mean of the total effects. To an extent, this allowed us to account for GPE uncertainty in the GSA.

We used the total effects to rank the parameters from most to least important by computing the maximum total effect across all outputs. The maximum total effects for each parameter were normalised to that they would sum up to 1 (representing 100% of output variance), and the parameters that collectively explained 90% of variation in simulator output were included in the analysis, while the others were fixed to values provided in Tables 13, since the simulator outputs were not sensitive to how these parameters are specified.

In order to make sure that the ToR-ORd-Land, Courtemanche-Land and the Eikonal models operated within physiological ranges, and to avoid wasting computational effort in the electromechanics simulator, we used HM. HM with GPEs is thoroughly explained in Longobardi et al. [11], and Fig 4 shows a schematic of our approach. First, the parameter space is explored with a Latin hypercube design. This provides a set of samples in the parameter space where the simulations are performed. These samples and outputs are then used to train a GPE. A test parameter set of 100,000 samples is then constructed with Latin hypercube sampling, and the GPE is evaluated at each sample. For each test sample, we evaluate the implausibility measure I(x) as: (7) (8) where fi is the GPE trained for the ith output. Scalars μi and σi are the measured target value and standard deviation for the ith output, available from the literature or the clinical data. By defining a threshold Ith on the implausibility measure I(x), the test samples are split in plausible (I(x) below the threshold) and implausible (I(x) above the threshold). From the implausible region, we then select a number Nsimul samples where the simulator is evaluated, and we construct a new test parameter set Ntest = 100000 for the next HM iteration (or wave). The GPEs are re-trained by adding the new simulated samples to the training set, and evaluated on the test set constructed at the previous wave. For each wave, we report the percentage of non-implausible points, the mean and maximum I and the mean and maximum ratio of GPE variance prediction over the measured variance σ2 (Vratio). The latter quantity gives a measure of how uncertain the GPEs are, relative to the uncertainty on the data we are trying to match. We set the final implausibility threshold to 3, based on Pukelsheim’s 3-sigma rule [46].

thumbnail
Fig 4. History matching.

Orange, blue and purple boxes indicate Latin hypercube sampling, simulator evaluation and GPE training, respectively. Blue and purple X indicate points where the simulator or the GPEs are evaluated.

https://doi.org/10.1371/journal.pcbi.1011257.g004

Sensitivity analysis and history matching on simulator components

Tables 13 list all parameters we initially considered. Excluding the values we could derive from the available clinical data, the total number of parameters was 117. This means that we would need to run over 1000 simulations for GPE training, which constitutes a considerable computational cost. To overcome this issue, we used GPEs, GSA and HM methods described above on the different simulator components to exclude unimportant parameters prior to performing the fully coupled simulator runs. Fig 5 summarises the analysis performed on the different sub-models, and below we describe it in detail.

thumbnail
Fig 5. Summary of sub-models analysis.

The diagram shows a summary of the analysis performed on the sub-models before performing the GSA on the four-chamber electromechanics model.

https://doi.org/10.1371/journal.pcbi.1011257.g005

In S2 File, we used GSA and parameter ranking to investigate which parameters affected the calcium transient as, in our simulator, the intracellular calcium is the only coupling variable between the ionic model and the Land model at the cellular level. This process reduced the number of parameters of the ToR-ORd model from 29 to 10. These 10 parameters were then combined with the 17 Land model parameters to identify which ionic and cellular mechanics parameters affected the active tension transient the most, as the cellular active tension drives contraction at the whole organ scale. The GSA on the coupled ToR-ORd-Land model led to a final set of 13 parameters to be considered in the fully coupled simulator. We applied a similar approach to the Courtemanche-Land model (S3 File) to identify 3 out of the 18 Courtemanche parameters and 10 out of the 17 Land model parameters as those affecting cellular atrial active tension the most. For the electrophysiology model at the tissue level, we performed a GSA on activation times only to exclude tissue parameters that did not affect the total atrial and ventricular activation time (S4 File). Due to the potential role fibres play in cardiac activation, we repeated the analysis presented in S4 File with ventricular transmural fibre orientation set to 50°/ − 50° and 70°/ − 70° ventricular fibre direction. There is high uncertainty on the size of the fast conducting regions in the atria and the ventricles (Fig 2). To show that our analysis on the electrophysiology model was independent of this modelling choice, we repeated the analysis with different sizes of the Bachmann bundle area in the atria and of the fast endocardial conduction layer in the ventricles (S10 File). Both analysis showed that the parameters we excluded following the GSA on the tissue electrophysiology model were independent of ventricular fibre direction (S9 File) and the size of the fast conducting regions (S10 File).

We performed a GSA on the passive mechanics model during inflation alone to exclude unimportant stiffness parameters (S5 File). We finally used the CircAdapt ODE model to study the effect of the circulatory system parameters on pressure-volume relationships of all four-chambers (S6 File), and to discard the ones that affected the simulated dynamics the least. The GSA on different simulator components allowed us to systematically reduce the number of free parameters from 117 to 45 prior to running the fully coupled simulator, and to reduce the number of simulations required for GPE training to a more feasible number. The final parameters are represented by the white rows in Table (Tables 13).

In addition to running a GSA on each simulator component, we used HM to restrict the important parameters to areas in the parameter space where either the outputs were physiological, or where the mechanics simulations were unlikely to be unstable and diverge, therefore leading to significant waste of computational resources. We constrained the ToR-ORd-Land model so that it provided physiological calcium transients based on the literature [47, 48] (S2 File), and where peak, rest and duration of the isometric tension transient were within ranges of measured values in mammals [49, 50] and allowed the simulations to achieve LV pressure peak and duration that were consistent with clinical data. Similarly, the calcium transient simulated with the Courtemanche model was constrained to be consistent with the literature [47, 48, 51] (S3 File), while the isometric active tension was bound to be physiological [9, 50, 5259]. We finally applied HM on the electrophysiology Eikonal model alone to achieve physiological total activation times of the atria and the ventricles, as the simulated activation times are independent of the mechanics simulation (S4 File).

In S11 File, we provide a detailed explanation of how we constructed the training samples for the GPEs Briefly, the last HM wave on the ToR-ORd-Land model, the Courtemanche-Land model and the electrophysiology tissue model provided N = 90906, N = 148527 and N = 99723 physiological samples, respectively. The bound for the parameters included in the HM was therefore implicitly defined by the plausible regions of the last HM wave (Tables 1 and 2, second column marked as HM). The bound for all the other parameters is provided in the second column of Tables 13. The AVdelay was set between 100 and 200 ms, consistently with ranges reported for CRT patients [38]. The range for ventricular stiffness parameters was set to ±50% from values reported in [9], to account for high variability in myocardium stiffness reported in the literature. The bounds for the bulk atrial stiffness a were increased based on higher collagen density in the atria compared to the ventricles [41]. Similarly, the range for the scaling parameters for the RV bulk stiffness and the reference tension was chosen to account for higher RV collagen content and lower myocyte density compared to the LV [41]. The length and the stiffness of the aorta were set to ±25% from their default value in CircAdapt [32], while the pulmonary and the systemic resistance factors were set to be between 1.0 and 4.0 to allow for physiological peak in pressure for the ventricles. The pericardium stiffness was bound between 0.5 and 2.0 kPa/mm, consistent with previous studies [4]. The shift on the LV end-diastolic pressure and the RV end-diastolic pressure were set according to measurements in AF patients in sinus rhythm [42]. For these parameters, we performed a Latin hypercube sampling with N = 90906 samples. These were then combined with the plausible regions from the HM on the ToR-ORd-Land model, the Courtemanche-Land model and the electrophysiology tissue model. Having identified the plausible region, we constructed a space filling design in this area on which to run the fully coupled electromechanics simulator and then to train the final GPEs. Because the plausible region had a complicated geometry, we used the function psa_select from the python library diversipy to extract N = 500 uniformly distributed samples.

We trained GPEs to predict the following outputs:

  1. LV/RV end-diastolic volume (EDV)
  2. LV/RV end-diastolic pressure (EDP)
  3. LV/RV end-systolic volume (ESV)
  4. LV/RV peak in pressure (pmax)
  5. LV/RV maximum pressure derivative (dp/dtmax)
  6. LV/RV minimum pressure derivative (dp/dtmin)
  7. LA/RA end-diastolic volume prior to atrial contraction (EDV)
  8. LA/RA end-systolic volume at the end of atrial contraction (ESV)
  9. LA/RA maximum volume during the v-wave (EDVvwave)
  10. LA/RA peak in pressure during atrial contraction (pmax).

The trained GPEs were then used to predict the model outputs and to run a GSA on the fully-coupled electromechanics framework, to investigate how the 45 parameters we considered affected whole heart function. During the GSA, we ensured that the GPEs were not evaluated outside their training space by generating a base sequence with Nbase = 2000 for each simulator sub-component. The GPEs trained to predict total ventricular and atrial activation times (S4 File), atrial and ventricular calcium and active tension transient features (S2 and S3 Files) were evaluated to exclude samples that provided unphysiological activation times, calcium and active tension. The base sequence was resampled for each simulator sub-component with an increasing initial number of samples until >2000 plausible samples were found. A Saltelli sampling was then generated, the samples were screened again with the GPEs for each simulator sub-component to ensure all samples provided physiological outputs. We finally quantified the total effects by using only the samples that were not excluded during the screening process. This allowed us to quantify the effect of simulator parameters on whole heart function while remaining within the GPE emulator training space.

The code to perform the GSA and the HM is provided at https://github.com/MarinaStrocchi/Strocchi_etal_2023_GSA. At this link, we also provide an example on how to train the GPEs, perform a GSA and rank the parameters for the ToRORd ionic model (S2 File).

Results

Fig 6 shows the pressure-volume loops for the LV (top), RV, LA and RA (bottom) that were simulated with the four-chamber electromechanics framework. The orange and the blue dots represent the simulated LV peak in pressure and LV pressure transient duration for each sample, respectively. The horizontal and the vertical lines indicate the average measured peak in pressure and pressure transient duration available from clinical data. The simulations achieve the patient’s peak in pressure and result in physiological pressure transient duration, thanks to the HM run on the ToR-ORd-Land model (S2 File). Although some simulations predict a low and unphysiological LV peak in pressure, we included these samples in the GPE training set because these simulations still provided information about the effect of the simulator input parameters and output, while parameter inference was outside the scope of this study. Out of 500 samples, 42 and 51 simulations failed during the unloading and the coupled simulation, respectively. This provided 407 simulations we extracted the pressure and volume features from to train the GPEs. Table 4 shows the mean R2 and the ISE scores for the GPEs traines for each output. The last column reports the mean scores across all five splits. For all outputs, the mean R2 and ISE were above 0.76 and 82.0, respectively, indicating that the GPEs provide an accurate estimate for all outputs.

thumbnail
Fig 6. Simulated pressure-volume transients.

Top: LV pressure-volume transients and pressure traces simulated with the fully-coupled simulator (black), with the orange and the blue circles representing the resulting LV peak pressure and transient duration. The horizontal and the vertical lines show the clinically measured average LV peak pressure and transient duration, respectively, to show that the simulator achieves physiological values. Bottom: RV, LA and RA simulated pressure-volume transients.

https://doi.org/10.1371/journal.pcbi.1011257.g006

thumbnail
Table 4. GPEs performance.

Average R2 score and ISE over the five cross-validation splits, reported for each output.

https://doi.org/10.1371/journal.pcbi.1011257.t004

The 407 model evaluations used for GPE training took an average of 6 h and 20 minutes per job on 512 cores, corresponding to a mean of 3247 core-hours per evaluation and 1321420 core-hours in total for all simulations. For the GSA, we performed 94000 GPE evaluations. Given the computational cost of the whole-heart simulations, it would not be feasible to perform these many model evaluations. Once trained, the GPEs require only ∼1s to estimate the outputs, allowing us to perform the GSA at a treatable computational cost, reduced from an estimate of >305 million core-hours to ∼1.3 million core-hours, due to running the simulations for GPE training. Therefore, thanks to the GPEs, we were able to achieve more than a 300 fold decrease in computational cost.

Fig 7 shows the heatmap for the total effect of all input parameters (y-axis) over all outputs (x-axis), with dark colors showing high interaction between inputs and outputs. Below, we provide a detailed analysis of the sensitivity analysis on the ventricular and atrial diastolic and systolic properties. We also show barplots of the total effects of important parameters for each output, multiplied by the sign of the linear regression coefficients βi from Eq (4) to give information about positive and negative interactions. In S12 File, we show that the standard deviation of the total effects computed from 1000 different samples of the posterior distribution of the GPEs is small, indicating that GPE uncertainty has negligible effects on the sensitivity indices we reported.

thumbnail
Fig 7. Global sensitivity analysis.

Total effects of input simulator parameters (y-axis) on model outputs (x-axis) normalised between 0 and 1 for each output.

https://doi.org/10.1371/journal.pcbi.1011257.g007

Ventricular diastolic and atrial systolic properties

Fig 8 shows the total order effects on the simulated end-diastolic volume and pressure of the ventricles (first two columns). Blue and red bars represent simulator parameters that are needed to explain 90% of output variation, with positive and negative effects on the output, respectively. Gray bars show parameters needed to explain up to 95% of output variation, while all other parameters are not displayed. The EDV of the ventricles was negatively affected by the unloading pressure of the LV (EDPLV,shift) and the RV (EDPunload,RV), and by the systemic (Rsys) and the pulmonary (Rpulm) resistances. A higher unloading pressure for the ventricles caused higher strains and therefore higher rest tension, that in turn reduced the EDV simulated at the limit cycle. The systemic resistance increased the LV EDV and the LV EDP. This was caused by increased LA ESV (and therefore decreased LA ejection) in response to increased systemic resistance. Similarly, the RV EDV was positively affected by the pulmonary resistance. The parameters that had a negative effect on the LV EDV had a positive effect on the RV EDV and vice-versa. This was caused by the interaction between the ventricles, modulated by the presence of the pericardium and by the intra-ventricular septum. When the LV EDV is bigger, the septum bulges towards the RV, leading to smaller RV EDV and therefore opposite effects of simulator parameters on the LV EDV and the RV EDV. Finally, passive material properties of the ventricles (aV and bt,V) positively affected the EDV of the LV. Although this might be counter-intuitive, decreased stiffness causes smaller unloaded volumes and higher initial strains at the beginning of the simulations, leading to higher ventricular rest tension and therefore decreased LV EDV.

thumbnail
Fig 8. The effect of simulator parameters on ventricular diastole and atrial systole.

Barplots of the total order effects of the input parameters, ranked from most to least important, for ventricular EDV, EDP and atrial ESV and peak pressure. The red and blue bars indicate parameters that explain >90% of the output variance, the gray bars indicate parameters that explain up to 95% of variance, while the rest of parameters are not displayed. Red and blue indicate positive and negative interactions, respectively. When this resulted in more than 10 parameters, the number of parameters was limited to 10 or to the number of parameters explaining 90% of output variance.

https://doi.org/10.1371/journal.pcbi.1011257.g008

The systemic and pulmonary peripheral resistances also affected atrial systolic peak in pressure (Fig 8, third and fourth columns). The LA peak in pressure was decreased by higher pulmonary resistance, caused by smaller LA preload (negative interaction between LA EDV in Fig 9). Similarly, increased systemic resistance led to smaller RA peak in pressure. The systemic resistance had also a significant positive impact on the LA peak in pressure. The RA pressure was also driven up by the pulmonary resistance, but with smaller interaction compared to the LA pressure and the systemic resistance. The unloading pressure (EDPunload,RV) and the reference tension scaling factor (Tref,LvRv) of the RV had minor positive effects on the LA peak in pressure. In addition to the systemic resistance, the RA peak in pressure was affected by cellular ionic and active tension parameters for velocity dependence (ϕA, Aeff,A), cross-bridge kinetics (rs,A, rw,A) and calcium handling (ca50,A, TRPN50,A and gCaL). Atrial stiffness (aA and bf,A) and atrial conduction velocities (CVf,A and kBB) also played a minor but significant role in RA peak pressure, with stiffer and slower myocardium causing lower RA pressure values.

thumbnail
Fig 9. The effect of simulator parameters on ventricular systole and atrial diastole.

Barplots of the total order effects of the input parameters, ranked from most to least important, for ventricular peak pressure, ESV and atrial maximum volume during the v-wave. The red and blue bars indicate parameters that explain >90% of the output variance, the gray bars indicate parameters that explain up to 95% of variance, while the rest of parameters are not displayed. Red and blue indicate positive and negative interactions, respectively. When this resulted in more than 10 parameters, the number of parameters was limited to 10 or to the number of parameters explaining 90% of output variance.

https://doi.org/10.1371/journal.pcbi.1011257.g009

The unloading pressures of the LV and RV affected the EDV of the ventricles the most, as they determined the end-diastolic strains and therefore the ventricular rest active tension. The systemic and pulmonary resistances had significant effects on the EDV as well, due to their effect on atrial ejection, while the AV delay mostly affected the EDP. Parameters for atrial active tension timing, CV in the ventricular endocardium and in the Bachmann bundle showed minor effects on the EDP of the ventricles, due to their effect on timing of RV initial contraction and atrial contraction rise and decay, respectively. Finally, the effect of parameters on the LV EDV was opposite to the effect on the RV EDV due to ventricular interaction modulated by the pericardium.

Ventricular systolic properties

The systemic and the pulmonary resistances affected the peak in pressure of the LV and the RV the most (Fig 9, top). Increased systemic and pulmonary resistances increase the afterload of the LV and the RV, therefore leading to higher LV and RV peak in pressure. The negative effect of the pulmonary resistance on the LV EDV caused smaller ventricular load, leading to a negative effect on the LV peak in pressure. Similarly, the RV peak in pressure was negatively affected by the systemic resistance. Higher bulk passive stiffness of the ventricles (aV) decreased the peak in pressure and increased ventricular ESV of both ventricles due to less significant deformation during active contraction. Parameters for ventricular cellular active tension also had a significant effect on the peak in pressure of the ventricles. Degree of velocity dependence (Aeff,V), cross-bridges kinetics parameter (rs,V, TRPN50 and μV) and the reference active tension (Tref,V) had significant effects on the peak in pressure, as they affected the timing and extent of cellular peak in cellular active tension, as shown in S2 File. Furthermore, the RV peak in pressure was positively affected by the LV-RV scaling of the reference tension Tref,LvRv.

LV ESV was positively affected by the systemic resistance and the unloading pressure of the RV. On the other hand, the unloading pressure of the LV negatively impacted LV ESV. The RV ESV was similarly affected by the unloading pressure of the RV and the LV, and by the pulmonary resistance. Increased peripheral resistance downstream to the LV and the RV led to higher ventricular afteload, therefore causing higher ESV. The passive stiffness of the ventricles aV and bt,V had minor positive effects on the ESV of both ventricles, therefore leading to smaller LV and RV ejection.

Ventricular pressure rise and decay rate were strongly affected by the systemic and the pulmonary resistances, but also by ventricular active tension parameters (Fig 10). Increased systemic resistance led to higher LV peak in pressure and therefore faster rise and decay in LV pressure. Similarly, RV rise and decay in pressure was strongly affected by the pulmonary resistance due to its effect on the RV peak in pressure. Parameters for ventricular cross-bridge kinetics (rs,V, rw,V, ku,V) affected both LV and RV maximum pressure derivative due to their effect on cellular active tension rising time and maximum tension derivative, as we showed in S2 File. Increased ventricular passive stiffness (aV) caused slower LV and RV pressure rise in both ventricles. In addition, LV pressure rise was impacted by the unloading pressures of the LV and the RV due to their effect on end-diastolic strains, by the rate of transition from blocked to unblocked cross-bridge binding sites (ku,V) and by ventricular velocity dependence (Aeff,V). The reference tension of the ventricles (Tref) and the scaling for the RV reference tension (Tref,LvRv) had an effect on the maximum pressure derivative of the LV and the RV, respectively, due to their effect on cellular active tension peak and timing. Apart from velocity dependence and ventricular stiffness, all important parameters for the LV and RV rise in pressure also impacted the LV and RV minimum pressure derivative. In addition, ventricular minimum pressure derivative was affected by the maximum troponin concentration () of the ToR-ORd model, due to its effect on calcium relaxation (S2 File), and by nTRPN,V and nTm,V, as they affected ventricular active tension decay at the cellular level.

thumbnail
Fig 10. The effect of simulator parameters on ventricular pressure rise and decay rate.

Barplots of the total order effects of the input parameters, ranked from most to least important, for ventricular maximum and minimum derivatives. The red and blue bars indicate parameters that explain >90% of the output variance, the gray bars indicate parameters that explain up to 95% of variance, while the rest of parameters are not displayed. Red and blue indicate positive and negative interactions, respectively. When this resulted in more than 10 parameters, the number of parameters was limited to 10 or to the number of parameters explaining 90% of output variance.

https://doi.org/10.1371/journal.pcbi.1011257.g010

The systolic dynamics of the ventricles are strictly related to atrial filling dynamics due to venous return. RV output was decreased by high pulmonary resistance. This in turn caused smaller LA maximum volume during the v-wave, e.g. smaller LA venous return. Similarly, due to decreased LV output, systemic venous return to the RA got smaller when the systemic resistance was high. Atrial stiffness parameters aA, bf,A and bt,A negatively affected venous return to the atria due to decreased atrial compliance. Finally, the unloading pressure and the passive stiffness of the ventricles aV had minor effects on venous return to the atria, due to their effect on ventricular output.

The systemic and the pulmonary resistances had the highest impact on ventricular systolic dynamics because they determine ventricular afterload. However, ventricular calcium and active tension parameters had significant impact on systolic function, in particular on pressure rise and decay. Atrial filling during the v-wave was affected by the peripheral resistances, but also by atrial stiffness parameters.

Discussion

In this study, we presented the first GSA on a four-chamber electromechanics simulator coupled with the circulatory system and accounting for atrial and ventricular contraction triggered by the calcium transient, as well as the effect of the pericardium. GPEs and GSA were used to identify important parameters for tissue electrophysiology, passive mechanics, the circulatory system and cellular dynamics. This analysis allowed us to reduce the number of parameters from 117 to 45. HM was then used to restrict the parameters for the cellular electromechanics simulator and tissue electrophysiology to areas in the parameter space where calcium, tension and total ventricular and atrial activation times were physiological. This provided 407 successful four-chamber electromechanics simulations to train the GPEs to predict features for whole organ function, and run a GSA on the fully coupled simulator. We showed that our analysis allows us to link cellular dynamics for calcium and tension to whole organ function. This framework can be used to improve our understanding of cardiac physiology and pathophysiology, and ultimately to provide novel insights into patient-specific treatment planning.

The effect of increased vascular resistance on heart function

Our results show that systemic and pulmonary peripheral resistance have significant effects on ventricular systolic function. In agreement with our analysis, Linde et al. [60] measured elevated systemic vascular resistance in hypertensive patients, contributing to elevated systolic pressure. Furthermore, drug studies reported on the efficacy of beta-blockers with vasodilatory activity for hypertension treatment, because the drug-induced decrease in the systemic resistance drives the LV systolic pressure down [61]. Increased systemic resistance was not only related to increased LV pressures, but also to decreased LV ejection in patients without known history of cardiovascular diseases [62], consistent with our analysis. Similarly, pulmonary hypertension patients have been reported to have increased pulmonary vascular resistance [63], therefore leading to increased RV pressure and reduced RV SV compared to normal subjects.

Not only systolic but also diastolic ventricular properties are affected by peripheral resistances. Harshaw et al. [64] explored the effect of a drug-induced decrease in systemic resistance in patients with mitral valve dysfunctions, showing a positive relationship between LV EDP and peripheral resistance, similar to our analysis. Further, patients with cirrhosis were reported to have increased left heart volumes and decreased right heart volumes due to reduced systemic resistance compared to controls, in agreement with our study [65]. RV EDV is elevated in pulmonary hypertensive patients [63] due to increased pulmonary vascular resistance, and was shown to decrease as a consequence of lower-body suction, aimed at reducing RV afterload through reduced peripheral resistance [66]. Several studies report on inverse effects of left and right heart volumes in response to increased resistance, with LV EDV decreasing with higher pulmonary resistance and RV EDV decreasing with higher systemic resistance [6567]. This has been attributed to diastolic ventricular interaction exerted by the pericardium [66], which causes an inverse relation between the LV EDV and the RV EDV, consistent with our results.

The systemic and the pulmonary resistances connect the left and the right sides of the hearts. Therefore, changes in the pulmonary and the systemic resistances induce alterations in the atria as well as in the ventricles. We have found a positive and negative linear correlation between the LA volumes and the systemic and the pulmonary resistances, respectively. Consistently with these findings, Marston et al. [68] reported a negative relation between pulmonary vascular resistance and LA volumes in patients with pulmonary hypertension. Furthermore, changes in pulmonary resistance following blood clots removal led to increased LA volumes and LV EDV. In addition, also consistent with our findings, in hypertensive patients, where the systemic resistance is higher than in normal subjects, the LA was dilated [69].

Increased or decreased systemic and pulmonary vascular resistances can occur in a wide range of pathologies, such as hypertension or cirrhosis. Vascular resistance has important effects on both the left and the right sides of the heart due to diastolic ventricular interaction mediated by the pericardium, and due to the interaction between atria and ventricles, mediated by the circulatory system. By coupling the three-dimensional mechanics framework with a closed-loop model for blood circulation and by including the effect of the pericardium on the heart, our simulator is able to account for these complex interactions and simulate how all four chambers respond to altered peripheral resistances.

Links between cell and whole organ function

Reduced contractility and consequent decreased atrial ejection during atrial contraction is a well recognized consequence of AF, even after cardioversion [70]. However, over the years, several studies have proposed different explanations for this mechanism. Active tension measurements in myocytes isolated from the right atrial appendage in AF patients have shown reduction of cellular contractile force compared to control, mainly caused by reduced the L-type calcium current [71, 72]. In agreement with these findings, our results showed increased atrial systolic volumes in response to lower L-type calcium ion channel conductance. Conduction slowing caused by loss of cell-to-cell coupling and/or by the presence of fibrosis might also contribute to depressed atrial contractility in AF patients [70, 73], consistent with the correlation between atrial conduction velocity and atrial end-systolic volume we found in our analysis. Finally, increased atrial stiffness, reported in AF patients [70] but also in other cardiac diseases, including hypertrophic cardiomyopathy [74], leads to impaired atrial compliance and smaller atrial end-diastolic volume, in agreement with our findings.

Efficiency of contraction can be measured using pressure biomarkers, such as peak in pressure or maximum pressure derivative, while pressure decay provides an indication of velocity of relaxation. Pressure relaxation was reported to be slower in diastolic HF patients [75], while systolic HF causes slower pressure rise due to sub-optimal contraction. Using our analysis, we were able to link these whole organ biomarkers to cellular function. We have shown that ventricular peak in pressure and pressure derivatives were significantly affected by tension velocity dependence, reference tension and cross-bridge kinetics parameters. Cellular measurements on isolated failing myocytes in animals and human tissue preparations showed that tension velocity dependence [76], cross-bridge cycling rate [77] and force development [78] were altered in failing subjects. Furthermore, diseased myocytes were reported to have slower calcium transient relaxation [48, 78], which, according to our analysis, might lead to slower pressure relaxation at the whole organ level.

Cardiac pathophysiology is complex, and can span across different length scales, from subcellular processes involving ionic channel remodelling and cross-bridge kinetics up to changes at the whole organ level. Although the parameter ranges we chose do not necessarily capture pathological variability caused by the cardiac diseases mentioned above, we have shown that our analysis can be used to link cellular function to biomarkers measured at the whole organ level. This demonstrates that multi-scale electromechanics models have the potential to improve our understanding of the pathophysology of a wide range of atrial and/or ventricular disorders.

Limitations

Although our study constitutes the first GSA on a four-chamber electromechanics model, accounting for multi-scale processes spanning from the sub-cellular to the whole organ level, it has some limitations.

Due to model complexity and to reduce computational costs, we did not run a fully converged Newton solution by constraining the maximum number of Newton iterations to 2. This introduced small errors in the simulated pressure and volume features below 3%, but it also allowed to speed up the simulations by up to 3 times. This reduced the computational cost of our GSA, making it possible to run more than 400 simulations and train GPEs that could accurately provide output predictions as function of 45 model parameters.

Our four-chamber electromechanics framework does not account for mechano-electrical feedback, assuming that mechanical deformation does not affect cellular and tissue electrophysiology. However, there are known mechanisms in which stretch induces changes in the electrical activity of the heart through, for instance, stretch-activated ion channels [79]. These mechanisms might play a role in a some pathologies, such as arrhythmogenesis initiation [80]. Therefore, depending on the application, additional feedback mechanisms might need to be included in the model, to provide a more accurate representation of the patient’s heart.

In our analysis, we only performed four-chamber simulations on one patient-specific anatomy due to their extensive computational cost. While the anatomy might impact the sensitivity analysis at the whole organ level, the analysis ran on the circulatory model, and the ventricular and atrial ionic and contraction models is independent of the geometry and therefore the conclusions for these models remain valid for other patients. In future, the anatomy could be accounted for in the GSA by performing a principal component analysis on a cohort of geometries and adding the principal vectors as additional parameters in the GSA, as in Rodero et al. [81]. In our analysis, we also fixed the fibre direction of atrial and ventricular myocardium, despite the high uncertainty in fibre direction reported in the literature. In S9 File, we showed that ventricular myofibre orientation has limited effects on the tissue electrophysiology GSA, passive mechanics and electromechanics model outputs, consistent with previous computational studies [82, 83]. However, we did not repeat the whole study with different fibre orientation, and we did not investigate the effect atrial myofibres as their arrangement is more complex and less well-established than that of ventricular fibre orientation. When more measurements of atrial myofibre orientation are available, it will be important to establish their effect on cardiac function.

To retrieve the stress-free configuration of the heart, we unloaded the ventricles but not the atria. This prevented excessive deformation in the atria during the unloading procedure, which makes the mechanics simulations more likely to diverge. Furthermore, we discarded diastolic residual active tension, even though this might play a role in relaxation dynamics. In future, residual active tension and atrial pressures could be accounted for in the unloading procedure [84], in order to have a more accurate estimation of the stress-free geometry.

In S2 and S3 Files, we performed GSA and HM on the ToR-ORd and Courtemanche models alone and coupled with the Land model to identify unimportant parameters for atrial and ventricular calcium and the active tension transients. The cell simulations used to train the GPEs in these instances were run with a basic cycle length of 1000 ms, corresponding to 1 Hz pacing, as opposed to using the basic cycle length derived from the clinical data. This was done to match the literature data we used as target values for the calcium transient features. Additionally, this makes the findings at the cellular level more general and applicable to other patients.

The GSA results we presented remain true within the parameter bounds and model outputs we included in the analysis. We only considered a subset of the 117 parameters that, within an expected physiological variability, had the largest impact on specific model outputs. Extreme changes in parameters, for example pathological or pharmacological inhibition of ion channels outside our assumed physiological variability, may result in other parameters outside of the 45 identified as being important in our analysis. Nevertheless, additional parameters could be re-introduced after the fitting process to investigate a particular cardiac pathology/treatment of interest.

In this study, we performed a GSA on five sub-models (S1S5 Files) to exclude parameters that did not affect sub-model outputs, under the assumption that these parameters would also be unimportant for whole-organ biomarkers. This could be done thanks to our knowledge of how the different sub-models are coupled within the four-chamber electromechanics framework. While this approach might fail to capture interactions between the parameters we excluded from the sub-models, this also allowed us to reduce the computational cost required to generate the simulations for emulators training and therefore the GSA. Nevertheless, since the parameters we excluded had small, if any, effects on sub-model outputs, we expect the combined effects between unimportant parameters of different sub-models to be small as well. Furthermore, in our GSA and HM, we did not consider any potential correlation between different outputs. Training a GPE for each output separately provided us the flexibility to fit the hyperparameters for each GPE independently. This could have been accounted for by training a multi-output GPE instead [85], where the GPE degrees of freedom are estimated for all outputs at the same time rather than independently. Furthermore, HM with a multi-variate GPE would have potentially helped in ruling out more implausible areas of the parameter space.

The GSA we presented does not account for the uncertainty of the multi-scale model, since quantifying model error is a more challenging task than learning the parameters of the model itself. In [86], Kennedy and O’Hagan showed how to account for model discrepancy, also referred to as model inadequacy, by adding it as an additional source of uncertainty. While the Kennedy and O’Hagan approach is simple and easy to apply, challenges arise during model parameter inference, since the observed data may be attributed either to the model itself or the model error, or both, meaning that the model parameters can no longer be consistently identified [87]. An elegant resolution involving orthogonal Gaussian processes was proposed in [88]. This approach is applicable on a linear model but, for non-linear models like the multi-scale four-chamber model presented in this study, retaining consistency of parameter estimates while accounting for model error remains an open research problem. Furthermore, although we sampled the GPEs posterior distribution and computed the average total effects from these samples, the sensitivity indices in our GSA were computed without directly accounting for GPE uncertainty. In future studies, the uncertainty of the GPE prediction could be accounted for by using alternative approximations of the sensitivity indices that directly include GPE standard deviation as well as the expected value.

Despite its limitation, our GSA allows us to link cellular processes with whole organ function, and has the potential to provide novel insights into patient pathophysiology and response to treatment that could never be tested in-vivo.

Conclusion

Our four-chamber electromechanics simulator is able to account for LV-RV and AV mechanical interaction, modulated by the presence of the pericardium and the coupling with a closed-loop model for the circulatory system. The GSA we presented allowed us to make considerations about how cellular dynamics translate into altered whole organ function. Thanks to the wide range of dynamics we can simulate, this analysis could potentially be applied to investigate a wide range of pathologies affecting the atria and the ventricles, and the consequences these have on all other chambers.

Supporting information

S1 File. Numerical tests.

We compared the model solution obtained with our numerical settings and a fully converged solution.

https://doi.org/10.1371/journal.pcbi.1011257.s001

(PDF)

S2 File. ToR-ORd-Land model sensitivity analysis and history matching.

We trained GPEs to predict calcium and active tension transient features simulated by the ToR-ORd model coupled with the Land contraction model, and used them to run a GSA to identify important parameters for these dynamics. Then, we used HM to isolate areas in the parameter space where the ventricular calcium and tension were physiological.

https://doi.org/10.1371/journal.pcbi.1011257.s002

(PDF)

S3 File. Courtemanche-Land model sensitivity analysis and history matching.

We trained GPEs to predict calcium and active tension transient features simulated by the Courtemanche model coupled with the Land contraction model, and used them to run a GSA to identify important parameters for these dynamics. Then, we used HM to isolate areas in the parameter space where the atrial calcium and tension were physiological.

https://doi.org/10.1371/journal.pcbi.1011257.s003

(PDF)

S4 File. Tissue electrophysiology sensitivity analysis and history matching.

We trained GPEs to predict the total atrial and ventricular activation times, and used them to run a GSA to identify important conduction parameters for simulated activation times. Then, we used HM to isolate areas in the parameter space where the total ventricular and atrial activation times were physiological.

https://doi.org/10.1371/journal.pcbi.1011257.s004

(PDF)

S5 File. Passive mechanics sensitivity analysis.

We trained GPEs to predict the maximum volumes for the four chambers and ventricular and atrial fibres strains during a passive inflation, and used them to run a GSA to identify important stiffness parameters.

https://doi.org/10.1371/journal.pcbi.1011257.s005

(PDF)

S6 File. CircAdapt sensitivity analysis.

We trained GPEs to predict a wide range of pressure and volume output features for all four chambers simulated with CircAdapt, and used them to run a GSA to identify important circulatory parameters for pressure and volume dynamics.

https://doi.org/10.1371/journal.pcbi.1011257.s006

(PDF)

S7 File. Saltelli sampling construction.

Detailed explanation of how we constructed the Saltelli samples for the GSA.

https://doi.org/10.1371/journal.pcbi.1011257.s007

(PDF)

S8 File. Gaussian processes emulators kernel comparison.

We investigated how the choice of kernel affected GPE accuracy and GSA. The GPE training and GSA performed for the ToR-ORd model (S2 File) were repeated with a Matérn rather than an exponentiated quadratic kernel and the results between the two cases were compared.

https://doi.org/10.1371/journal.pcbi.1011257.s008

(PDF)

S9 File. The effect of ventricular fibre orientation.

We performed the GSA on the electrophysiology tissue model with two different ventricular fibre orientation to investigate the effect of myofibre arrengement on electrophysiology simulations and GSA. We also repeated a passive inflation and a four-chamber electromechanics simulation to quantify the effect of ventricular fibres on mechanics model outputs.

https://doi.org/10.1371/journal.pcbi.1011257.s009

(PDF)

S10 File. The effect of fast conducting regions size on the electrophysiology model.

We performed the GSA on the electrophysiology tissue model altering the size of the fast endocardial conducting layer of the ventricles and the Bachmann bundle.

https://doi.org/10.1371/journal.pcbi.1011257.s010

(PDF)

S11 File. Non-implausible region extraction.

Schematic of how the non-implausible areas from the HM of the sub-models were used to construct the training samples for the GPEs at the whole organ scale.

https://doi.org/10.1371/journal.pcbi.1011257.s011

(PDF)

S12 File. The effect of emulators uncertainty on sensitivity indices.

We show the mean and standard deviation of the total effects computed for the whole organ sensitivity analysis obtained when sampling the posterior distribution of the emulators.

https://doi.org/10.1371/journal.pcbi.1011257.s012

(PDF)

Acknowledgments

This work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).

References

  1. 1. Niederer SA, Lumens J, Trayanova NA. Computational models in cardiology. Nature Reviews Cardiology. 2019;16(2):100–111. pmid:30361497
  2. 2. Strocchi M, Gsell MA, Augustin CM, Razeghi O, Roney CH, Prassl AJ, et al. Simulating ventricular systolic motion in a four-chamber heart model with spatially varying Robin boundary conditions to model the effect of the pericardium. Journal of Biomechanics. 2020;101:109645. pmid:32014305
  3. 3. Strocchi M, Augustin CM, Gsell MA, Karabelas E, Neic A, Gillette K, et al. A publicly available virtual cohort of four-chamber heart meshes for cardiac electro-mechanics simulations. PLoS One. 2020;15(6):e0235145. pmid:32589679
  4. 4. Pfaller MR, Hörmann JM, Weigl M, Nagler A, Chabiniok R, Bertoglio C, et al. The importance of the pericardium for cardiac biomechanics: From physiology to computational modeling. Biomechanics and modeling in mechanobiology. 2018; p. 1–27. pmid:30535650
  5. 5. Augustin CM, Neic A, Liebmann M, Prassl AJ, Niederer SA, Haase G, et al. Anatomically accurate high resolution modeling of human whole heart electromechanics: a strongly scalable algebraic multigrid solver method for nonlinear deformation. J Comput Phys. 2016;305:622–646. pmid:26819483
  6. 6. Krishnamurthy A, Gonzales MJ, Sturgeon G, Segars WP, McCulloch AD. Biomechanics simulations using cubic Hermite meshes with extraordinary nodes for isogeometric cardiac modeling. Computer aided geometric design. 2016;43:27–38. pmid:27182096
  7. 7. Fritz T, Wieners C, Seemann G, Steen H, Dössel O. Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomech Model Mechanobiol. 2014;13(3):627–641. pmid:23990017
  8. 8. Gerach T, Schuler S, Fröhlich J, Lindner L, Kovacheva E, Moss R, et al. Electro-mechanical whole-heart digital twins: a fully coupled multi-physics approach. Mathematics. 2021;9(11):1247.
  9. 9. Land S, Niederer SA. Influence of atrial contraction dynamics on cardiac function. Int J Numer Method Biomed Eng. 2018;34(3):e2931–e2931. pmid:28990354
  10. 10. Salvador M, Dede L, Manzoni A. Non intrusive reduced order modeling of parametrized PDEs by kernel POD and neural networks. Computers & Mathematics with Applications. 2021;104:1–13.
  11. 11. Longobardi S, Lewalle A, Coveney S, Sjaastad I, Espe EK, Louch WE, et al. Predicting left ventricular contractile function via Gaussian process emulation in aortic-banded rats. Philosophical Transactions of the Royal Society A. 2020;378(2173):20190334. pmid:32448071
  12. 12. Galappaththige S, Gray RA, Costa CM, Niederer S, Pathmanathan P. Credibility assessment of patient-specific computational modeling using patient-specific cardiac modeling as an exemplar. PLoS computational biology. 2022;18(10):e1010541. pmid:36215228
  13. 13. Neic A, Gsell MA, Karabelas E, Prassl AJ, Plank G. Automating image-based mesh generation and manipulation tasks in cardiac modeling workflows using meshtool. SoftwareX. 2020;11:100454. pmid:32607406
  14. 14. Bayer JD, Blake RC, Plank G, Trayanova NA. A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann Biomed Eng. 2012;40(10):2243–2254. pmid:22648575
  15. 15. Labarthe S, Bayer J, Coudière Y, Henry J, Cochet H, Jais P, et al. A bilayer model of human atria: mathematical background, construction, and assessment. Europace. 2014;16(suppl4):iv21–iv29. pmid:25362166
  16. 16. Roney CH, Pashaei A, Meo M, Dubois R, Boyle PM, Trayanova NA, et al. Universal atrial coordinates applied to visualisation, registration and construction of patient specific meshes. Medical image analysis. 2019;55:65–75. pmid:31026761
  17. 17. Neic A, Campos FO, Prassl AJ, Niederer SA, Bishop MJ, Vigmond EJ, et al. Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model. J Comput Phys. 2017;. pmid:28819329
  18. 18. Lee AW, Nguyen UC, Razeghi O, Gould J, Sidhu BS, Sieniewicz B, et al. A rule-based method for predicting the electrical activation of the heart with cardiac resynchronization therapy from non-invasive clinical data. Medical Image Analysis. 2019;57:197–213. pmid:31326854
  19. 19. Strocchi M, Lee AW, Neic A, Bouyssier J, Gillette K, Plank G, et al. His Bundle and Left Bundle Pacing with Optimised Atrio-ventricular Delay Achieve Superior Electrical Synchrony over Endocardial and Epicardial Pacing in Left Bundle Branch Block Patients. Heart Rhythm. 2020;. pmid:32603781
  20. 20. Tomek J, Bueno-Orovio A, Passini E, Zhou X, Minchole A, Britton O, et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. Elife. 2019;8:e48890. pmid:31868580
  21. 21. Tomek J, Bueno-Orovio A, Rodriguez B. ToR-ORd-dynCl: an update of the ToR-ORd model of human ventricular cardiomyocyte with dynamic intracellular chloride. BioRxiv. 2020;.
  22. 22. Courtemanche M, Ramirez RJ, Nattel S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. American Journal of Physiology-Heart and Circulatory Physiology. 1998;275(1):H301–H321. pmid:9688927
  23. 23. Land S, Park-Holohan SJ, Smith NP, dos Remedios CG, Kentish JC, Niederer SA. A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. J Mol Cell Cardiol. 2017;106:68–83. pmid:28392437
  24. 24. Coveney S, Clayton RH. Sensitivity and uncertainty analysis of two human atrial cardiac cell models using Gaussian process emulators. Frontiers in Physiology. 2020;11:364. pmid:32390867
  25. 25. Mirams GR, Pathmanathan P, Gray RA, Challenor P, Clayton RH. Uncertainty and variability in computational and mathematical models of cardiac physiology. The Journal of physiology. 2016;594(23):6833–6847. pmid:26990229
  26. 26. Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. Journal of biomechanics. 1995;28(10):1167–1177. pmid:8550635
  27. 27. Flory PJ. Thermodynamic relations for high elastic materials. Transactions of the Faraday Society. 1961;57:829–838.
  28. 28. Ogden RW. Nearly isochoric elastic deformations: application to rubberlike solids. J Mech Phys Solids. 1978;26(1):37–57.
  29. 29. Augustin CM, Gsell MA, Karabelas E, Willemen E, Prinzen FW, Lumens J, et al. A computationally efficient physiologically comprehensive 3D–0D closed-loop model of the heart and circulation. Computer methods in applied mechanics and engineering. 2021;386:114092. pmid:34630765
  30. 30. Strocchi M, Augustin CM, Gsell MA, Karabelas E, Neic A, Gillette K, et al. The Effect of Ventricular Myofibre Orientation on Atrial Dynamics. In: International Conference on Functional Imaging and Modeling of the Heart. Springer; 2021. p. 659–670.
  31. 31. Arts T, Delhaas T, Bovendeerd P, Verbeek X, Prinzen FW. Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model. American Journal of Physiology-Heart and Circulatory Physiology. 2005;57(4):H1943–H1943. pmid:15550528
  32. 32. Walmsley J, Arts T, Derval N, Bordachar P, Cochet H, Ploux S, et al. Fast simulation of mechanical heterogeneity in the electrically asynchronous heart using the multipatch module. PLoS computational biology. 2015;11(7):e1004284–e1004284. pmid:26204520
  33. 33. Bols J, Degroote J, Trachet B, Verhegghe B, Segers P, Vierendeels J. A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels. J Comput Appl Math. 2013;246:10–17.
  34. 34. Shi W, Jantsch M, Aljabar P, Pizarro L, Bai W, Wang H, et al. Temporal sparse free-form deformations. Medical image analysis. 2013;17(7):779–789. pmid:23743085
  35. 35. Mensel B, Kühn JP, Schneider T, Quadrat A, Hegenscheid K. Mean thoracic aortic wall thickness determination by cine MRI with steady-state free precession: validation with dark blood imaging. Acad Radiol. 2013;20(8):1004–1008. pmid:23830606
  36. 36. Bacakova L, Travnickova M, Filova E, Matějka R, Stepanovska J, Musilkova J, et al. The role of vascular smooth muscle cells in the physiology and pathophysiology of blood vessels. Muscle Cell and Tissue-Current Status of Research Field. 2018;1:13.
  37. 37. Vigmond EJ, Hughes M, Plank G, Leon LJ. Computational tools for modeling electrical activity in cardiac tissue. J Electrocardiol. 2003;36:69–74. pmid:14716595
  38. 38. Hyde ER, Behar JM, Crozier A, Claridge S, Jackson T, Sohal M, et al. Improvement of right ventricular hemodynamics with left ventricular endocardial pacing during cardiac resynchronization therapy. Pacing Clin Electrophysiol. 2016;39(6):531–541. pmid:27001004
  39. 39. O’Hara T, Virág L, Varró A, Rudy Y. Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS computational biology. 2011;7(5). pmid:21637795
  40. 40. Nasopoulou A, Shetty A, Lee J, Nordsletten D, Rinaldi CA, Lamata P, et al. Improved identifiability of myocardial material parameters by an energy-based cost function. Biomech Model Mechanobiol. 2017;16(3):971–988. pmid:28188386
  41. 41. Oken DE, Boucek RJ. Quantitation of collagen in human myocardium. Circulation research. 1957;5(4):357–361. pmid:13447177
  42. 42. Alboni P, SCARFŒ S, Fuca G, Paparella N, Yannacopulu P. Hemodynamics of idiopathic paroxysmal atrial fibrillation. Pacing and Clinical Electrophysiology. 1995;18(5):980–985. pmid:7659571
  43. 43. Rasmussen CE, Williams CK, et al. Gaussian processes for machine learning. vol. 1. Springer; 2006.
  44. 44. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer physics communications. 2010;181(2):259–270.
  45. 45. Saltelli A. Making best use of model evaluations to compute sensitivity indices. Computer physics communications. 2002;145(2):280–297.
  46. 46. Pukelsheim F. The three sigma rule. The American Statistician. 1994;48(2):88–91.
  47. 47. Coppini R, Ferrantini C, Yao L, Fan P, Del Lungo M, Stillitano F, et al. Late sodium current inhibition reverses electromechanical dysfunction in human hypertrophic cardiomyopathy. Circulation. 2013;127(5):575–584. pmid:23271797
  48. 48. Piacentino V III, Weber CR, Chen X, Weisser-Thomas J, Margulies KB, Bers DM, et al. Cellular basis of abnormal calcium transients of failing human ventricular myocytes. Circulation research. 2003;92(6):651–658.
  49. 49. Backx PH, Gao W, Azan-Backx MD, Marban E. The relationship between contractile force and intracellular [Ca2+] in intact rat cardiac trabeculae. The Journal of General Physiology. 1995;105(1):1–19. pmid:7730787
  50. 50. Land S, Niederer SA, Aronsen JM, Espe EKS, Zhang L, Louch WE, et al. An analysis of deformation-dependent electromechanical coupling in the mouse heart. The Journal of physiology. 2012;590(18):4553–4569. pmid:22615436
  51. 51. Brixius K, Pietsch M, Hoischen S, Muller-Ehmsen J, Schwinger RH. Effect of inotropic interventions on contraction and Ca2+ transients in the human heart. Journal of Applied Physiology. 1997;83(2):652–660. pmid:9262464
  52. 52. Niederer S, Hunter P, Smith N. A quantitative analysis of cardiac myocyte relaxation: a simulation study. Biophysical journal. 2006;90(5):1697–1722. pmid:16339881
  53. 53. Allen D, Kurihara S. The effects of muscle length on intracellular calcium transients in mammalian cardiac muscle. The Journal of physiology. 1982;327(1):79–94. pmid:7120151
  54. 54. Kentish JC. Combined inhibitory actions of acidosis and phosphate on maximum force production in rat skinned cardiac muscle. Pflügers Archiv. 1991;419(3):310–318. pmid:1745606
  55. 55. Hinken AC, McDonald KS. Inorganic phosphate speeds loaded shortening in rat skinned cardiac myocytes. American Journal of Physiology-Cell Physiology. 2004;287(2):C500–C507. pmid:15084471
  56. 56. Wolska BM, Vijayan K, Arteaga GM, Konhilas JP, Phillips RM, Kim R, et al. Expression of slow skeletal troponin I in adult transgenic mouse heart muscle reduces the force decline observed during acidic conditions. The Journal of physiology. 2001;536(3):863–870. pmid:11691878
  57. 57. Hibberd M, Jewell B. Calcium-and length-dependent force production in rat ventricular muscle. The Journal of physiology. 1982;329(1):527–540. pmid:7143258
  58. 58. Ebus J, Papp Z, Zaremba R, Stienen G. Effects of MgATP on ATP utilization and force under normal and simulated ischaemic conditions in rat cardiac trabeculae. Pflügers Archiv. 2001;443(1):102–111. pmid:11692273
  59. 59. Papp Z, Szabó Á, Barends JP, Stienen G. The mechanism of the force enhancement by MgADP under simulated ischaemic conditions in rat cardiac myocytes. The Journal of physiology. 2002;543(1):177–189. pmid:12181290
  60. 60. Linde T, Sandhagen B, Hägg A, Mörlin C, Wikström B, Danielson BG. Blood viscosity and peripheral vascular resistance in patients with untreated essential hypertension. Journal of hypertension. 1993;11(7):731–736. pmid:8228192
  61. 61. Brett S, Forte P, Chowienczyk P, Benjamin N, Ritter J. Comparison of the effects of nebivolol and bisoprolol on systemic vascular resistance in patients with essential hypertension. Clinical drug investigation. 2002;22(6):355–359.
  62. 62. Mandry D, Girerd N, Lamiral Z, Huttin O, Filippetti L, Micard E, et al. Relationship Between Left Ventricular Ejection Fraction Variation and Systemic Vascular Resistance: A Prospective Cardiovascular Magnetic Resonance Study. Frontiers in Cardiovascular Medicine. 2021;8. pmid:35004914
  63. 63. Claessen G, La Gerche A, Dymarkowski S, Claus P, Delcroix M, Heidbuchel H. Pulmonary vascular and right ventricular reserve in patients with normalized resting hemodynamics after pulmonary endarterectomy. Journal of the American Heart Association. 2015;4(3):e001602. pmid:25801760
  64. 64. HARSHAW CW, GROSSMAN W, MUNRO AB, MCLAURIN LP. Reduced systemic vascular resistance as therapy for severe mitral regurgitation of valvular origin. Annals of internal medicine. 1975;83(3):312–316. pmid:1180426
  65. 65. Møller S, Sondergaard L, Mogelvang J, Henriksen O, Henriksen JH. Decreased right heart blood volume determined by magnetic resonance imaging: evidence of central underfilling in cirrhosis. Hepatology. 1995;22(2):472–478. pmid:7635415
  66. 66. Atherton JJ, Moore TD, Lele SS, Thomson HL, Galbraith AJ, Belenkie I, et al. Diastolic ventricular interaction in chronic heart failure. The Lancet. 1997;349(9067):1720–1724. pmid:9193381
  67. 67. Cecconi M, Johnston E, Rhodes A. What role does the right side of the heart play in circulation? Critical care. 2006;10(3):1–7. pmid:17164017
  68. 68. Marston NA, Auger WR, Madani MM, Kimura BJ, Strachan G, Raisinghani AB, et al. Assessment of left atrial volume before and after pulmonary thromboendarterectomy in chronic thromboembolic pulmonary hypertension. Cardiovascular ultrasound. 2014;12(1):1–5. pmid:25109313
  69. 69. Eshoo S, Ross DL, Thomas L. Impact of mild hypertension on left atrial size and function. Circulation: Cardiovascular Imaging. 2009;2(2):93–99. pmid:19808574
  70. 70. Heijman J, Linz D, Schotten U. Dynamics of atrial fibrillation mechanisms and comorbidities. Annual review of physiology. 2021;83:83–106. pmid:33064962
  71. 71. Schotten U, Ausma J, Stellbrink C, Sabatschus I, Vogel M, Frechen D, et al. Cellular mechanisms of depressed atrial contractility in patients with chronic atrial fibrillation. Circulation. 2001;103(5):691–698. pmid:11156881
  72. 72. Schotten U, Duytschaever M, Ausma J, Eijsbouts S, Neuberger HR, Allessie M. Electrical and contractile remodeling during the first days of atrial fibrillation go hand in hand. Circulation. 2003;107(10):1433–1439. pmid:12642366
  73. 73. Veeraraghavan R, Gourdie RG, Poelzing S. Mechanisms of cardiac conduction: a history of revisions. American Journal of Physiology-Heart and Circulatory Physiology. 2014;306(5):H619–H627. pmid:24414064
  74. 74. Sanada H, Shimizu M, Sugihara N, Shimizu K, Ino H, Takeda R. Increased left atrial chamber stiffness in hypertrophic cardiomyopathy. Heart. 1993;69(1):31–35. pmid:8457391
  75. 75. Zile MR, Baicu CF, Gaasch WH. Diastolic heart failure—abnormalities in active relaxation and passive stiffness of the left ventricle. New England Journal of Medicine. 2004;350(19):1953–1959. pmid:15128895
  76. 76. Capasso JM, Fitzpatrick D, Anversa P. Cellular mechanisms of ventricular failure: myocyte kinetics and geometry with age. American Journal of Physiology-Heart and Circulatory Physiology. 1992;262(6):H1770–H1781. pmid:1621835
  77. 77. Hajjar RJ, Gwathmey JK. Cross-bridge dynamics in human ventricular myocardium. Regulation of contractility in the failing heart. Circulation. 1992;86(6):1819–1826. pmid:1451254
  78. 78. Pérez NG, Hashimoto K, McCune S, Altschuld RA, Marbán E. Origin of contractile dysfunction in heart failure: calcium cycling versus myofilaments. Circulation. 1999;99(8):1077–1083. pmid:10051303
  79. 79. Kohl P, Hunter P, Noble D. Stretch-induced changes in heart rate and rhythm: clinical observations, experiments and mathematical models. Progress in biophysics and molecular biology. 1999;71(1):91–138. pmid:10070213
  80. 80. Peyronnet R, Nerbonne JM, Kohl P. Cardiac mechano-gated ion channels and arrhythmias. Circulation research. 2016;118(2):311–329. pmid:26838316
  81. 81. Rodero C, Strocchi M, Marciniak M, Longobardi S, Whitaker J, O’Neill MD, et al. Linking statistical shape models and simulated function in the healthy adult human heart. PLoS computational biology. 2021;17(4):e1008851. pmid:33857152
  82. 82. Niederer SA, Lamata P, Plank G, Chinchapatnam P, Ginks M, Rhode K, et al. Analyses of the redistribution of work following cardiac resynchronisation therapy in a patient specific model. PLoS ONE. 2012;7(8). pmid:22952697
  83. 83. He J, Pertsov AM, Cherry EM, Fenton FH, Roney CH, Niederer SA, et al. Fiber Organization has Little Effect on Electrical Activation Patterns during Focal Arrhythmias in the Left Atrium. IEEE Transactions on Biomedical Engineering. 2022;.
  84. 84. Xi J, Lamata P, Niederer S, Land S, Shi W, Zhuang X, et al. The estimation of patient-specific cardiac diastolic functions from clinical measurements. Medical image analysis. 2013;17(2):133–146. pmid:23153619
  85. 85. Conti S, Gosling JP, Oakley JE, O’Hagan A. Gaussian process emulation of dynamic computer codes. Biometrika. 2009;96(3):663–676.
  86. 86. Kennedy MC, O’Hagan A. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2001;63(3):425–464.
  87. 87. Brynjarsdóttir J, O’Hagan A. Learning about physical parameters: The importance of model discrepancy. Inverse problems. 2014;30(11):114007.
  88. 88. Plumlee M, Joseph VR. Orthogonal Gaussian process models. Statistica Sinica. 2018; p. 601–619.