Sensitivity analysis of a strongly-coupled human-based electromechanical cardiac model: Effect of mechanical parameters on physiologically relevant biomarkers.

The human heart beats as a result of multiscale nonlinear dynamics coupling subcellular to whole organ processes, achieving electrophysiologically-driven mechanical contraction. Computational cardiac modelling and simulation have achieved a great degree of maturity, both in terms of mathematical models of underlying biophysical processes and the development of simulation software. In this study, we present the detailed description of a human-based physiologically-based, and fully-coupled ventricular electromechanical modelling and simulation framework, and a sensitivity analysis focused on its mechanical properties. The biophysical detail of the model, from ionic to whole-organ, is crucial to enable future simulations of disease and drug action. Key novelties include the coupling of state-of-the-art human-based electrophysiology membrane kinetics, excitation-contraction and active contraction models, and the incorporation of a pre-stress model to allow for pre-stressing and pre-loading the ventricles in a dynamical regime. Through high performance computing simulations, we demonstrate that 50% to 200% - 1000% variations in key parameters result in changes in clinically-relevant mechanical biomarkers ranging from diseased to healthy values in clinical studies. Furthermore mechanical biomarkers are primarily affected by only one or two parameters. Specifically, ejection fraction is dominated by the scaling parameter of the active tension model and its scaling parameter in the normal direction ( k ort 2 ); the end systolic pressure is dominated by the pressure at which the ejection phase is triggered ( P ej ) and the compliance of the Windkessel fluid model ( C ); and the longitudinal fractional shortening is dominated by the fibre angle ( ϕ ) and k ort 2 . The wall thickening does not seem to be clearly dominated by any of the considered input parameters. In summary, this study presents in detail the description and implementation of a human-based coupled electromechanical modelling and simulation framework, and a high performance computing study on the sensitivity of mechanical biomarkers to key model parameters. The tools and knowledge generated enable future investigations into disease and drug action on human ventricles.

Introduction

Scope. Cardiac disease is the leading cause of death worldwide, affecting millions of people every year [1]. Disease generally affects the human heart through complex mechanisms, involving from subcellular processes such as ionic currents, to whole-organ properties such as tissue structure; and impairing, in one way or another, the heart’s primary role of pumping blood through mechanical contraction and relaxation driven by electrophysiological activity. The great complexity of the pathophysiological mechanisms involved in cardiac electromechanical function is difficult to assess in clinical and experimental settings, due to ethical and practical limitations, and also limited access to human tissue. Novel human-based approaches are therefore required to accelerate very much needed improvements in diagnosis and treatment [2].

Computational multiscale modelling and simulations provide a powerful framework to dissect the highly nonlinear processes involved in cardiac electromechanical disease. In cardiac science, computational modelling and simulation has been an active area of research, at the forefront of computational biomedicine [3]. For the last five decades, cell electrophysiology, excitation–contraction coupling and active tension models have been developed in close iteration with experiments. Detailed models and simulations of human ventricular electrophysiology have replicated experimental and clinical recordings in a range of healthy and disease conditions, and also under drug action (see for instance [3], [4], [5]). Their maturity has triggered interest and impact beyond academia, such as the adoption of the state-of-the-art O’Hara–Rudy (ORd) model [6] for industrial and regulatory purposes, within the CiPA initiative sponsored by the US Food and Drug Administration [7]. Studies have shown prediction of drug-induced clinical arrhythmic risk in human with 89% accuracy while data obtained from previously conducted animal studies for similar datasets showed up to 75% accuracy [8]. Furthermore, recently a new model of human-based ventricular active contraction was proposed and evaluated by Land et al. [9] using experimental recordings obtained in human tissue preparations. However, the coupling and incorporation of these two models in the framework of organ-level ventricular simulations has not been achieved yet.

Progress has also been impressive in enabling simulations considering multi-physics coupling including electrophysiology and mechanics for whole heart dynamics (see for instance [10], [11], [12], [13], [14], [15], amongst others). Some of these studies have also considered fully coupled electro-mechanico-fluidic ventricular dynamics, where an explicit representation of the fluid via the Navier–Stokes equations is considered [13], [14], [16]. Building on previous work mainly from [14], [17], [18], [19], [20], in this study we describe a human-based physiologically-detailed, and fully-coupled ventricular electromechanical model together with a suitable simulation framework and a sensitivity analysis focusing on its mechanical properties. One of the novelties of this formulation is the level of detail considered (from ionic to whole-organ), which is crucial to eventually enable simulations of disease and drug action. Other key additions include coupled state-of-the-art human-based electrophysiology membrane kinetics, excitation–contraction and active contraction models [9], [21], and the implementation of a dedicated model that allows us to consider pre-stressing and pre-loading of the ventricles in a dynamical regime. This overcomes some of the limitations associated with one or more of the following approaches: quasi-static instead of dynamic formulations (i.e. ignoring inertial effects, which may indeed be relevant in cardiac electromechanics [22], [23]), simplified models for cell electrophysiology and/or active tension [24], [25], the lack of pre-stress and/or lack of representation of the four phases of the cardiac cycle [18], [26], and lack of an appropriate representation of epicardial mechanical boundary conditions [27], [28]. Moreover, to enable the reproducibility of the results, we also include a full description of the numerical schemes and their computational implementation, giving sufficient detail for each model component.

Furthermore, through high performance computing (HPC), we conduct a sensitivity analysis in which we quantify variations in simulated mechanical biomarkers caused by varying key parameters in the human ventricular electromechanical model. Whereas previous studies have investigated the variability generated by electrophysiological parameters [29], [30], the role of variability in mechanical parameters on mechanical biomarkers in the human heart has not been systematically assessed yet. After recalling classical approaches to deal with uncertainty in parameter values, we perform HPC simulations varying key mechanical parameters within appropriate ranges, we then quantify physiologically-relevant biomarkers, and we finally validate our results against clinically-relevant values.

The mathematical notation defined in this section largely follows the one used in [13], [31], adopting tensor notation throughout this study. Unless otherwise specified, scalars are denoted with Greek or Latin italic characters (e.g.  or , respectively); vectors, or first-order tensors, are denoted by Latin bold lower-case characters (e.g. ); second-order tensors are denoted with Greek or Latin bold upper-case characters (e.g.  or , respectively); and fourth-order tensors are denoted by Latin double-barred upper-case characters (e.g. ). Reference and spatial coordinates are respectively indicated with and . Reference and spatial configurations of a body are respectively denoted with and , with and being their corresponding boundaries. Subscripts are used for naming purposes and, when relevant, they are followed, after a space, by a time step, iteration number, or additional naming subscripts; superscripts are also used for naming purposes. Additional notation details are given in Appendix A.

Materials and methods

Geometry of the left ventricular model

A simplified left ventricular geometry is modelled as an ellipsoid, truncated at the base. The geometry was obtained from the cardiac mechanics benchmark in [32], and isotropically scaled to obtain an end-diastolic volume (EDV) similar to that of the human left ventricle (EDV 159.43 ml) [33]. The finite element (FE) mesh is composed of linear tetrahedra, and has a maximum edge length of 0.025 cm, resulting in 33M elements and 5.6M nodes; the used step size ranges from 0.001 to 0.01 ms. Such mesh size is needed for the convergence of numerical algorithms as discussed below.

Description of the governing equations

The main components of the human-based ventricular electromechanical model are: (1) the anisotropic description of the mechanical behaviour of ventricular muscle (Fig. 2.1 Blue “blob”), (2) the electrical propagation model (Fig. 2.1 Green “blob”), (3) the cell electrophysiology model including calcium dynamics (Fig. 2.1 Yellow “blob”), and (4) the active-tension model coupled to the excitation–contraction model (Fig. 2.1 Red “blob”). As in [13], the present subsection describes physiological models and coupling mechanisms for each of these phenomena, primarily focusing on their mathematical structure. A complete description of the specific way of coupling human-based electrophysiology and active tension models is one of the novelties of this study, and the results are detailed and thoroughly discussed in Section ??.

Fig. 2.1

Diagram of the human ventricular electromechanical modelling and simulation framework, which includes: (blue) nonlinear solid mechanics, (green) electrical propagation, (yellow) ventricular cell electrophysiology, and (red) active-tension plus excitation-contraction. The arrows represent the different coupling mechanisms and their corresponding directions, with the quantities written in red being exchanged between the different components of the model. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Ventricular mechanical activity

The end-diastolic configuration of the left ventricle is considered as the region at , also referred to as reference, or material, configuration, with the following boundaries: epicardium and left ventricular endocardium (Fig. 2.2). The motion of the ventricles is defined through the deformation operator which maps the position of every location in the ventricles onto the spatial, or current, configuration for a time . The deformation gradient tensor field is defined as , and maps tangent vectors from the reference configuration onto the spatial configuration.

Fig. 2.2

Left ventricular geometry, including local fibre architecture and domain boundaries. The centreline of the ventricular cavity is aligned with the -axis.

The mechanical behaviour of ventricular tissue is characterised by the balance of linear momentum in a total Lagrangian formulation, considering inertial contributions, but ignoring volumetric forces (e.g. gravitational effects or body loads), such that where is the (initial) density of the considered material in the reference configuration, is the field of displacements, is the deformation gradient , is the second Piola–Kirchhoff stress, is the pressure applied to the endocardial surface and obtained through the Windkessel model described below, is the stiffness corresponding to the elastic spring boundary condition applied to the epicardium, and is the outward normal defined on the whole material boundary . Note that in Eq. (2.1c) , where could be seen as a linear operator which projects to the subspace spanned by and then scales it by ; this would act as a spring which is positioned normally to the epicardial surface , and with constant stiffness .

The considered constitutive law for the passive mechanical behaviour of ventricular tissue is a nearly, or quasi-, incompressible version of Holzapfel and Ogden [34]; justified by in-vivo measurements of the material properties of the myocardium [35], [36], it is assumed that the ventricular tissue admits a small degree of compressibility. The strain energy density function defining the orthotropic and hyperelastic continuum is where its volumetric, and isochoric isotropic/anisotropic contributions are respectively defined as and , , and ; and are the fibre and sheet directions in the reference configuration. These, together with the normal direction in the reference configuration form a right-handed orthonormal set of basis vectors (note the over-bar, which emphasises the difference between the normal to the fibre-sheet plane and the normal to a boundary of the considered domain, as previously specified). The parameters have stress dimensions, whereas are dimensionless parameters; all considered positive, and being the bulk modulus. The term is the th invariant of the right Cauchy–Green strain tensor , defined as From these specifications, the precise form of the passive contribution to the second Piola–Kirchhoff stress tensor can be obtained as Equation (2.1b) describes the effect of blood flow on the left ventricle. This was implemented through a state machine with four (five counting the initiation) phases, one for each of the cardiac phases during a heartbeat: isovolumetric contraction, ejection, isovolumetric relaxation, and filling; specified as follows. For additional details on the Windkessel model, see e.g. [24].

Initiation. This phase is needed in order to bring the system to the end of diastole. The ventricular pressure applied to the endocardium is brought up to healthy left ventricular end-diastolic pressure linearly throughout the phase. In order to keep the volume constant during this phase, a pre-stress is assumed to follow the fibre direction and is initialised at with a relatively small value and is further adjusted during this phase via the relation where and can be considered as the inverse of penalty terms for volume difference and volume rate, respectively. The last summand is used as a stabilisation term in order to avoid oscillations of the ventricular pressure due to inertial effects. The pre-stress is kept constant for the rest of simulation, at the converged value obtained in the last time step of this phase.

Isovolumetric contraction. At this phase the aortic valve is closed, and the ventricular pressure increases due to contraction of the ventricles caused by electrical activation, while keeping the volume constant. The ventricular pressure required to maintain the condition of constant ventricular volume is obtained via where and can be considered as the inverse of penalty terms for volume difference and volume rate, respectively. Again, the last summand is used as stabilisation against spurious oscillations of the ventricular pressure, which might occur due to inertial effects.

Ejection. When the ventricular pressure surpasses the arterial pressure, the aortic valve opens and the ejection phase begins, eventually leading to a reduction in ventricular volume. To model the blood pressure of the systemic circulation system during the ejection phase, the following two element Windkessel model is used (see e.g. [27]) where defines the compliance, and defines the impedance of the systemic circulation system. Here stands for arterial pressure.

Isovolumetric relaxation. When the ventricular flow reverses, i.e. , the isovolumetric relaxation phase begins. In this phase, the ventricular pressure decreases while keeping the ventricular volume constant at the end-systolic volume (ESV). The pressure is therefore obtained analogously to the isovolumetric contraction phase from (2.3).

Filling. The isovolumetric relaxation phase ends when the pressure drops below a specified threshold. Blood flows from the atria to the ventricle while the ventricular volume returns to its initial value (end-diastolic volume). In this phase, the ventricular pressure is modelled with the following decay equation where is the decay constant.

The left ventricular volume is calculated throughout the cardiac cycle by using the divergence theorem and the assumption that the volume is constrained by a flat lid located on the basal plane [13], leading to where . Because of the orientation of the reference state of the ventricle and of its centreline, here we consider that on the basal plane.

Electrical propagation model

Simulation of propagation of electrical excitation was achieved through the monodomain equations [37] in the form: where is the transmembrane potential, are the gating variables that regulate the transmembrane currents, are the ionic concentrations inside of the cell, is the right-hand side of the system of ODEs corresponding to the generic state variable vector , is the surface to volume ratio, is the membrane capacitance per unit area, are the ionic currents, and is the applied current triggering the electrical depolarisation acting on specific anatomical locations. The orthotropic tensor of local conductivities in the reference configuration is defined as where , and are, respectively, the conductivities in the fibre, sheet and normal directions.

Human-based ventricular cellular electromechanical model

The electrophysiological activity at the cellular level is encoded in (2.6b), (2.6c) and in the term of (2.6a). Eqs. (2.6b), (2.6c) describe the membrane kinetics as well as the local transport of ions within cells. In this study, the modified version by [21] of the human ventricular cell electrophysiology model by O’Hara et al. [6] was used. It consists of a system of ODEs with 41 state variables, and it exhibits a highly stiff behaviour, mainly due to the activation of sodium ionic current, which causes a sudden spike in the transmembrane potential . This is the human model that shows closer agreement with experimental recordings in a range of conditions and that has been adopted as the basis for the CiPA initiative for drug assessment, sponsored by the Food and Drug Association (FDA) for the in-silico [7].

Coupling of the electrophysiology to the contractile machinery of the cardiomyocytes is represented through the human-based model of excitation–contraction and active tension by [9]. The model is based on sets of ODEs describing the local dynamics of a vector of state variables , which represents the contractile mechanisms in cardiomyocytes (see a detailed discussion about these models in [38]). The corresponding system of ODEs reads as (note the dependencies on stretch and stretch-rate) where is the stretch in the fibre direction. The Land et al. [9] model consists of six state variables : and are variables associated with the crossbridge binding, respectively being the post-powerstroke and pre-powerstroke states, represents the fraction of troponin C units with calcium bound to their regulatory binding site, represents the fraction of blocked myosin binding sites on actin, and the state variables dictate the pre-powerstroke and post-powerstroke distortion in a distortion-decay model.

The right-hand side terms defining in (2.8) assume the following form where the state variable-dependent parameters are defined as When the solution of the system of ODEs (2.8) is achieved, the expression for the active tension in the fibre direction can be retrieved as where and the remaining parameters are constant.

Equations (2.6c), (2.8) are coupled bidirectionally, as specified by (2.8) and such that (2.6c) now becomes Particularly, the calcium bound to troponin needed by the cell electrophysiology model is obtained from the Land active tension model, which means that depends on as (a more complete derivation, following [6], and a specification of model constants are shown in Appendix B) where and . On the other hand, the calcium concentration needed by the active tension model is obtained from the cell electrophysiology model, i.e.  depends on specifically through the evolution of the calcium bound to troponin [9], in the form

Description of the electromechanical coupling mechanisms

The electrophysiological and mechanical activities of the heart at the ventricular level are bidirectionally coupled through active stress (from electrical diffusion to solid deformation), the effect of deformation on diffusivity and stretch-activated ionic currents (the two latter are also known as mechano-electric feedback, i.e. from solid deformation to electrical diffusion).

Active stress

As an active stress formulation is used, the stress term in (2.1a) is additively decomposed into passive and active terms, such that the balance of linear momentum reads Here the active Second Piola–Kirchhoff stress tensor is defined as where and are activation parameters acting in the sheet and normal directions, respectively. These constants are taken as in [39], [40], [41], in which an orthotropic active stress tensor is assumed, considering that mechanical activation occurs differently in each local direction. The specific values of the weighting activation parameters can contribute to achieving physiological regimes of wall thickening and ejection fraction. It is also important to point out that in (2.11) the pre-stress given by (2.2) appears only in the fibre direction.

Mechano-electric feedback

Mechano-electric feedback generally comprises of two coupling submechanisms in which the mechanics of cardiac tissue affect the electrophysiological activity. By adding stretch-activated ionic currents to the spatial formulation of (2.6a) and then pulling back to the reference configuration (see e.g. [42]), one obtains The orthotropic tensor of local conductivities (2.7) is now defined as where and are the stretches in the sheet and normal directions, respectively. The term in (2.12) denotes stretch-activated ionic currents defined as where and are, respectively, the reversal potential and the maximal conductance of the channels [43], [44]. Note that this term only becomes active when the tissue is stretched in the fibre direction, i.e. . In this study is set to 0, regardless of the strain state.

It is worthy to note a few aspects of the chosen formulation (2.12)–(2.13). The pulled-back counterpart of the conductivity tensor defined on the spatial domain (2.13) is which assumes that the conductivity coefficients have been measured in the spatial configuration. On the other hand, one could assume conductivity coefficients measured in the reference configuration, which would lead to the conductivity tensor in the reference configuration As one can observe, the equivalences between the latter conductivity coefficients in (2.15), (2.16) are which means that if the conductivity coefficients are fixed in one configuration, they vary in the other configuration according to the corresponding strain distribution. In particular, the formulation chosen for this study (2.12)–(2.13) assumes conductivity coefficients measured in the spatial configuration; this means that, for instance, under an isovolumetric extension in the fibre direction the conductivity coefficient in the reference configuration would be reduced by a factor . A more exhaustive description of the effect of deformation on diffusivity, its different formulations and implications can be seen in [13].

Weak formulation

The first step towards obtaining a numerical solution of the proposed nonlinear coupled problem is to define an adequate weak formulation of Eqs. (2.10), (2.1b), (2.1c), (2.12), (2.6b), (2.6c), (2.8), posed in the reference configuration. Assuming sufficient regularity for initial conditions and source terms, for a fixed time, one proceeds to contract the governing equations by sufficiently regular test functions, to integrate by parts the divergence terms, and to use appropriately the boundary conditions to obtain the following variational problem:

For , find (where the power denotes the dimension of the vector of each species ), such that where the functional spaces of the test functions in (2.17a)–(2.17b), (2.17c)–(2.17e) are respectively specified as and note that these coincide with the trial spaces in (2.17a), (2.17b) since no part of the boundary assumes essential boundary conditions.

The problem under consideration is highly nonlinear and multiscale, and it features intricate interdependencies between submodels. Obtaining approximate solutions based on fully monolithic approaches is therefore not feasible as the number of state variables of the physiologically-relevant cell electrophysiology models can be large (frequently these models also feature Heaviside functions, which compromise the well-definiteness of Jacobian operators). Instead, a semi-implicit partitioned approach is proposed, and is made precise in what follows. The use of segregated approaches is additionally advantageous as the propagation of electrical depolarisation consists of a moving wave with a steep wavefront (with a width of fractions of millimetre, and with a duration of the upstroke of the action potential of 2 ms), while the tissue deformation does not exhibit such a particularly local behaviour [45], [46]. This would allow, for instance, the use of different meshes: a finer mesh for electrical propagation and a coarser mesh for solid mechanics, although this has not been done in this study.

The detailed description of the numerical implementation of the model, including the finite element discretisation, the time discretisation, operator splitting and subsequent linearisation can be found in Appendix D, Appendix E.

Software and simulation platform

The software used in this study is Alya, developed at the Barcelona Supercomputing Center (BSC) [14], [17], [18], [20], [47]. Alya is a multi-physics code that works in parallel, specially designed to efficiently simulate strongly coupled problems in large-scale HPC platforms. Electrical diffusion and mechanical deformation are both solved in the same code instance, in a staggered manner. Parallelisation is based on a hybrid message passing interface (MPI)/OpenMP paradigm, in which the mesh is partitioned at a pre-processing stage by using METIS [48]. With respect to its scalability, an example can be found in [19], where the largest simulation so far of a cardiac electromechanical model of a biventricular geometry is performed, with a mesh consisting of 3.4 billion tetrahedra, reaching an almost linear scalability up to 100,000 cores in NCSA’s Blue Waters. Verification tests have been conducted as part of previous studies [20], [47].

The 3D ventricular simulations were run on the supercomputer MareNostrum IV (BSC), by using 2400 MPI processes (equivalent to 2400 computing cores), and would take from 2 to 40 hours (including the initiation phase), depending on the considered time step. The assigned computational time was enough to ensure that both isovolumetric contraction and ejection phases were completed for a single heartbeat. This enables calculation of the considered clinically-relevant mechanical biomarkers.

Initiation of the simulations. The human-based coupled electromechanical cellular model was iterated in a single instance before starting the FE simulation itself. The pacing was chosen as 70 bpm during 1000 beats, the latter being, according to previous experience by the authors, enough to bring a cellular model of similar characteristics to a “periodic steady-state” [21]. Stretch and stretch rate were respectively set to 1 and 0 during the initiation of the cellular model. The results of the state variables of the cell electrophysiology model were assigned to each node of the FE mesh as their initial values ; likewise, the results of the state variables of the active tension model were assigned to each integration point of the FE mesh (as their initial values ). The activation sequence of the reaction–diffusion starts with the activation of the endocardium and then progresses by diffusion to the rest of the myocardium. A current stimulus was applied for 2 ms, which enabled activation of the ventricle.

Model parameters. Myocardial tissue was modelled as an orthotropic material for both mechanical deformation and electrical diffusion. The corresponding diffusivity values were 0.001 cm, 0.0005 cm and 0.00025 cm for, respectively, the fibre, sheet and normal directions [49]. Fibre directions are parallel to the basal plane in the middle of the cardiac wall, and vary linearly up to forming an angle with the basal plane at the endocardium/epicardium [50], and tangent to the circumferential direction. Sheet directions are parallel to the vectors normal to the endocardial/epicardial surface, and do not vary throughout the cardiac wall. Other parameter values will be detailed in Section 2.6.

Sensitivity analysis

Choice of sensitivity parameters

The parameters chosen for the sensitivity analysis were the following: stiffness of the epicardium Robin boundary condition , end-diastolic pressure , pressure at which the ejection phase is triggered , compliance of the 2-elements Windkessel model , resistance of the 2-elements Windkessel model , orthotropic activation parameter in the normal direction , active stress scaling parameter of the active tension model , the angle of fibre distribution , bulk modulus , and the linear parameter (i.e. that with stress magnitude) of the passive strain energy density function in the sheet direction . These parameters were chosen amongst the whole set of mechanical parameters due to a combination of the following factors: (1) uncertainty of their values; (2) feasibility of the study; (3) their effect on the mechanical biomarkers in cellular and tissue sensitivity studies. With respect to passive properties, simulations were carried out to investigate which of the mechanical parameters of the strain energy density function had a substantial impact on the mechanical biomarkers. For this, the so-called linear parameters (those with stress dimension) in the Holzapfel–Ogden strain energy density function [34] were considered. Variations of 50%–200% in single parameters with respect to the baseline values were tested, and only was found to result in variations larger than 5% in EF. Furthermore, in the cellular active contraction model, was identified as the most critical parameter for active tension, and therefore was considered in the sensitivity analysis.

General approach

The sensitivity analysis was divided into two parts: a one-at-a-time approach and a partial rank correlation coefficient (PRCC) sampling-based global approach. One-at-a-time approaches involve variation of a single parameter while keeping the rest fixed. Sampling-based approaches involve the generation and exploration of a mapping from the model parameters, sampled in this study by Latin Hypercube Sampling (LHS) [51], to the desired mechanical biomarkers (known in the uncertainty quantification and sensitivity analysis literature as “quantities of interest”) [52]. Then, the importance of the individual parameters with respect to the uncertainty in the mechanical biomarkers is assessed by the analysis of the PRCC, as proposed in [53], [54]. Uniform probability density distributions were assumed for the considered input parameters, as parameter distributions are unknown [55]. PRCC-based sensitivity analysis studies are not new in the context of cardiac electrophysiology [56], [57], [58].

Correlation provides a measure of the relationship between model input parameters and output biomarkers. A correlation coefficient (CC) between an input parameter and a biomarker is where The CC varies from −1 to +1, where a positive value would indicate that and tend to increase, or decrease, together and a negative value would indicate the opposite; a value of 0 indicates the presence of no linear relationship. The partial correlation coefficient (PCC) between and can be defined as follows. First, the regression models are constructed. The results of these linear regressions are used to defined the new variables and . Partial correlation characterises the linear relationships between an input parameter and a biomarker after a correction for the linear effects on of the remaining input parameters is made. The PCC between and are defined as the CC (2.18) between and , i.e. . Similarly to PCC, partial rank correlation performs correlation on rank-transformed data, which means that and are replaced by their corresponding ranks (e.g. the smallest value is assigned the rank of 1, the following largest value is assigned the rank of 2, up to the largest value, which is assigned the rank of ). PRCC is a robust sensitivity coefficient for nonlinear but monotonic relationships between and , as long as little to no correlation exists between the input parameters [53], [59], [60].

A one-at-a-time approach was considered for the following two reasons, (1) as a preliminary exploration of the parameter space, with the purpose of approximately assessing the potential monotonicity and nonlinearity of the mechanical biomarkers with respect to each individual model parameter; (2) to probe suitable ranges of model parameters for the subsequent PRCC-based global sensitivity analysis. Thus simulations varying one parameter at a time informed the global approach for the sensitivity analysis using LHS. The one-at-a-time approach considered ten values for each of the parameters, while keeping the others constant. Therefore, a total number of samples equal to 100, ten samples per parameter, was considered. These values for the parameters were obtained as where and are, respectively, the lower and upper bounds of the parameter ranges. A summary of the parameters included in the sensitivity analysis and the ranges of considered values are shown in Table 2.1; the baseline values of these parameters are shown in Table 2.2. These ranges of values generally extended from 50% of the baseline value to 200%–1000% of the baseline value, with the exception of , whose values were completely unknown; , whose maximum value was obtained from [40], [41] and not increased further due to the stability of the numerical treatment (a higher value of this parameter could compromise the positive-definiteness of the tangent operator, hampering the convergence of the Newton–Raphson scheme); and the fibre angle , whose values were directly obtained from the estimations made in [50].

Table 2.1

Parameters included in the sensitivity analysis, their ranges of values, units, (sub)model where they belong, the domain to which the corresponding (sub)model is applied to, and reference(s) from where the ranges of values are inferred.

Parameter	Range of values	Units	(Sub)Model	Applied to	Ref.
k_epi	0–1,000,000	Ba cm−1	Robin	∂Ω_0epi	N/A
P_0	7.5–30	mmHg	Windkessel	∂Ω_0endo	[24], [27]
P_ej	33.75–135	mmHg	Windkessel	∂Ω_0endo	[24]
C	0.1–0.4	ml mmHg−1	Windkessel	∂Ω_0endo	[27]
R	75–750	mmHg ms ml−1	Windkessel	∂Ω_0endo	[27]
k_ort2	0.0–0.6	N/A	Orthotropic activation	Ω_0	[40]
T_ref	60–1200	kPa	Active stress	Ω_0	[9]
ϕ	45–90	deg	Fibre distribution	Ω_0	[50]
K	500–5000	kPa	Passive stress	Ω_0	[27]
a_s	512.8–7692	Pa	Passive stress	Ω_0	[27]

Table 2.2

Values of the baseline parameters.

Parameter	k_epi (Ba cm−1)	P_0 (mmHg)	P_ej (mmHg)	C (ml⋅mmHg−1)	K (kPa)
Value	200,000	20	67.5	0.2	3333
Parameter	R (mmHg⋅ms⋅ml−1)	k_ort2	T_ref (kPa)	ϕ (deg)	a_s (Pa)
Value	750	0	120	60	2564

The sampling of the parameter space for the global sensitivity analysis considered LHS, with the limits defined by the ranges of values in Table 2.1. A number of samples equal to 10 times the number of parameters were also considered, leading to 100 samples (200 in total). This number of LHS samples was deemed suitable as there are studies in the literature that used a number of LHS samples of the same order for a much higher number of independent parameters [61], [62].

Choice of mechanical biomarkers

The present electromechanical simulations aim to reproduce clinically-relevant ventricular activity (and to then be further used for simulations under diseased conditions or drug action). A suitable set of outputs or mechanical biomarkers was chosen. These included ejection fraction (EF), which represents the amount of blood that is pumped during a heartbeat end-systolic pressure (ESP), which is the maximum ventricular pressure achieved in the left ventricle during a heartbeat; longitudinal fractional shortening (LFS) [63], which is a fractional version of the relative displacement between the endocardial apex and the base. LFS is evaluated as where and are, respectively, the endocardial apico-basal distance at the end of the ejection phase and the initial (or end-diastolic) endocardial apico-basal distance (i.e. the “length” of the cavity of the truncated ellipsoid); and wall thickening (WT) [64], which is the fractional version of the relative displacement between points in the endocardium and epicardium that are at the same position relative to the apico-basal axis where and are, respectively, the cardiac wall thickness evaluated at the end of the ejection phase and the initial cardiac wall thickness. and are evaluated in this study via the displacements of two points which are at the same angle relative to the apico-basal axis and at the same height relative to the endocardial apex, at the beginning of the simulation, one of them is on the endocardium and the other on the epicardium; the thickness values were obtained as the difference of the distance of these points to the apico-basal axis.

The physiological ranges of these mechanical biomarkers for healthy subjects are shown in Table 2.3.

Table 2.3

Mechanical biomarker values for healthy subjects. Values are shown as either full ranges, or 95% confidence interval; if the reference showed mean SD, these were converted to 95% confidence intervals through . Experimental data distributions are assumed to be normal to allow for such extrapolation.

Mechanical biomarker	Values	Units	Ref.
EF	48–69	%	[33]
ESP	100–174	mmHg	[65]
LFS	13–21	%	[66]
WT	18–100	%	[67]

Results

Effect of varying one-at-a-time model parameters on mechanical biomarkers

A typical pressure–volume curve showing a complete heartbeat is shown in Fig. 3.1b; different snapshots of ventricular morphology during such heartbeat are shown (Fig. 3.1a). A visual representation of WT and LFS is also shown to aid the reader in their visual representation (Fig. 3.1c). As expected, the volume during isovolumetric contraction remains constant with the corresponding pressure build up. A reduction in volume and a peak in pressure during the ejection phase follow. When the ventricle starts relaxing and the volume flow reverses, isovolumetric relaxation begins, which implies constant volume and a further reduction in pressure while the myocardium continues repolarising. The final phase, filling phase leads to the slow recovery of the initial volume while intraventricular pressure decreases.

Fig. 3.1

(a) Sequence of snapshots of contraction throughout the whole heartbeat, the colour code indicates transmembrane potential (red, fully activated and depolarised, and blue, polarised); (b) pressure–volume loop throughout the whole heartbeat, A–E correspond to the snapshots in (a); (c) a sketch depicting wall thickness and apico-basal distance (which are the absolute counterparts of wall thickening (WT) and longitudinal fractional shortening (LFS)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The values of the mechanical biomarkers for the one-at-a-time variations () of the model parameters () are shown in Fig. 3.2, with the clinically-reported healthy ranges of the mechanical biomarkers shown in green. It can be observed that the relationships between parameters and biomarkers are monotonic (and nonlinear) in almost all of the cases. It can be seen that the parameter which affects the mechanical biomarkers more strongly is (Fig. 3.2b, c, f, h), with a particularly strong effect on EF (Fig. 3.2b). By just varying this parameter to high values, healthy values of EF () are attained. Other relationships shown to be important are (Fig. 3.2a), (Fig. 3.2b), (Fig. 3.2c), (Fig. 3.2c), (Fig. 3.2d), (Fig. 3.2e), (Fig. 3.2f), and (Fig. 3.2f).

Fig. 3.2

Mechanical biomarkers versus one-at-a-time variations of the considered parameters. (a) Ejection fraction versus . (b) Ejection fraction versus . (c) End-systolic pressure versus . (d) End-systolic pressure versus . (e) Longitudinal fractional shortening versus . (f) Longitudinal fractional shortening versus . (g) Wall-thickening versus . (h) Wall-thickening versus . Areas shaded in green represent the healthy ranges of the considered mechanical biomarkers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Results show that increasing increases the haemodynamic load of the left ventricle and delays the transition between isovolumetric contraction and ejection phases, directly affecting the maximum pressure obtained during a cardiac cycle (i.e. directly increasing ESP, Fig. 3.2c). Increasing and has opposite effects. The haemodynamic load of the ventricle is reduced, thus making the ventricle more compliant, as is increased or is decreased; this leads to increased EF (Fig. 3.2a–b) and decreased ESP (Fig. 3.2c–d). By just modifying one of these two parameters, healthy values of ESP can be attained (Fig. 3.2c–d). Increasing leads to a higher stiffness of the pericardium, which is the membrane surrounding the heart; LFS is increased (Fig. 3.2e) because of the restricted motion of the apex and the resulting “sliding”-like contractile motion (i.e. the epicardial surface slides over its initial configuration, Fig. 3.3). The parameter (the angle between the basal plane and the fibre direction in the epicardial/endocardial surface) has a similar effect on LFS. Increasing leads to the active tension distribution in the fibre direction to more closely follow the apico-basal direction in the endocardium and epicardium, causing an increase in apico-basal motion (and thus an increased LFS, Fig. 3.2f). Healthy values of LFS are not achieved by one-at-a-time variations of the considered input parameters (Fig. 3.2e–h). The parameter is the fraction of active tension in the fibre direction which is apportioned to the direction. The direction is normal to the fibre-sheet plane and the sheet direction is parallel to the epicardial/endocardial surface normal. Higher leads to higher tension in the apico-basal direction throughout the whole ventricle, increasing the apico-basal motion in a similar fashion to when increasing (and therefore, increasing LFS, Fig. 3.2f).

Fig. 3.3

(Left) Mechanical deformation of the left ventricular geometry from diastole (transparent grey) to systole (red) for two cases with different mechanical parameter values. Left panel corresponds to low pericardial stiffness () or low fibre angle (), whereas right panel results from high or high . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Global sensitivity analysis results

The combined effect of varying the ten model parameters on the mechanical biomarkers was evaluated via a global PRCC- and LHS-based sampling approach. The histograms showing the probabilistic distribution of the considered mechanical biomarkers are shown in Fig. 3.4. The shaded green area corresponds to the healthy range of the corresponding biomarker, shown in Table 2.3. The grey shaded area corresponds to the low ejection fraction zone (35% is the threshold for implantable cardioverter-defibrillator under certain clinical guidelines [68]). The dashed red lines are mean values for exemplar diseases. As it can be noted, the obtained values for the mechanical biomarkers show a skewed distribution. The peak of the distribution of EF lies on the lower range of healthy values; the peak of the distribution of ESP lies in the middle of the range of healthy values; and the peak of the distribution of LFS lies between the mean of the diseased values and the lower threshold for healthy values. Therefore, the developed numerical model can reproduce a wide variety of outcomes, including those typically considered “diseased”, and those within healthy conditions. The obtained values of WT are . Although this range of values does intersect the healthy range shown in Table 2.3, these values lie on the lower range of healthy values. The healthy ranges of WT in literature have a considerable spread [64], [69], [70]; even within the same study the healthy range can vary enormously due to the considered method of evaluation [70], due to inter-individual variability [67], or even due to intra-individual variability [71]. Such sparsity of WT ranges in the literature does not allow for a straightforward comparison with the results obtained in the present study.

Fig. 3.4

Histograms of simulated mechanical biomarkers compared to clinical values. The green shaded area indicates clinical ranges for healthy patients (Table 2.3). (a) Ejection fraction. The grey shaded area corresponds to low ejection fraction (35% is the threshold for implantable cardioverter-defibrillator indicated by clinical guidelines [68]). The dashed red line indicates the clinical mean value for heart failure patients. [72]. (b) End-systolic pressure. The dashed red line indicates the clinical mean value for dilated cardiomyopathy patients [65]. (c) Longitudinal fractional shortening. The dashed red lines indicate the clinical mean values for patients with coronary artery disease [73], left, and high blood pressure with supranormal ejection fraction [74], right. (d) Wall thickening. The dashed red line corresponds to the clinical mean value for patients with ischaemic heart disease with dyskinetic segments [75]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3.5 shows the PRCC values which are statistically significant (); “plus” or “minus” symbols are included in the significant PRCC, indicating a positive or a negative correlation, respectively. Results are consistent with the sensitivity analysis varying just one parameter at a time (Fig. 3.2). EF is positively correlated with (moderate correlation) and negatively correlated with (moderate correlation); ESP is positively correlated with (strong correlation) and negatively correlated with (moderate correlation); LFS is positively correlated with (weak correlation), (weak correlation), and (strong correlation); WT is negatively correlated with and (both with moderate correlation); all the other statistically significant correlation coefficients are relatively weak in comparison with the mentioned ones, i.e. closer to zero. Fig. 3.5 also suggests that EF, ESP and LFS could be dominated by mainly two parameters each.

Fig. 3.5

Correlation plot showing the partial rank correlation coefficients between each of the varied parameters and each of the simulated mechanical biomarkers. White grids have almost zero correlation coefficients and thus an associated -value lower than the threshold for statistical significance (). “Plus” or “minus” symbols only appear in the statistically significant values of the correlation coefficients and indicate a positive or a negative correlation, respectively.

Fig. 3.6 shows the values of EF, ESP and LFS with respect to their most influential parameters according to the PRCC values in Fig. 3.5. EF is the most important mechanical biomarker, as it is usually employed to classify cardiac function into healthy or non-healthy. In this figure, a distinction between “healthy” cases (EF 40%) and “diseased” cases (EF 40%) is made for graphical purposes, according to clinical guidelines [76]. It can be seen that healthy cases are scattered throughout the whole range of ESP and LFS, which suggests that no significant relationship exists between EF and ESP/LFS (Fig. 3.6c–f). Usually, higher contraction is associated with a higher pressure build-up (i.e. EF is usually associated with higher ESP), which is also suggested by the maximum and minimum ESP values being respectively placed in high and low EF (Fig. 3.6a,b). However, one may see that there exist “healthy” cases with either low or high ESP (Fig. 3.6c, d).

Fig. 3.6

Scatter-plots of the mechanical biomarkers with respect to their most influential model parameters. (a) Ejection fraction versus . (b) Ejection fraction versus . (c) End-systolic pressure versus . (d) End-systolic pressure versus . (e) Longitudinal fractional shortening versus . (f) Longitudinal fractional shortening versus . Black circles denote “healthy” cases (ejection fraction higher than 40%), whereas grey circles denotes “diseased” cases (ejection fraction lower than 40%). Coloured circles denote cases that maximise a mechanical biomarker (red for ejection fraction, green for end-systolic pressure and blue for longitudinal fractional shortening). On the other hand, coloured crosses denote cases that minimise a mechanical biomarker. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

In this study, we present a detailed description and a sensitivity analysis of a biophysically-detailed,bidirectionally-coupled ventricular electromechanical modelling and simulation framework including human-based electrophysiology, excitation–contraction coupling and active tension models. The computational model includes the four main components of a human ventricular electromechanical model: (1) solid mechanics, whose passive behaviour is based on the quasi-incompressible version of the Holzapfel and Ogden model [34], (2) the monodomain model for electrical propagation, (3) human ventricular cell electrophysiology using the modified version of the ORd model from [21], and (4) human active tension model as in [9], coupled to solid mechanics through an active stress framework. Two key mechanisms of mechano-electric feedback are also incorporated including how tissue deformation affects electrical diffusion, and the stretch-activated ionic currents. The main novelties and improvements of our framework, compared to previous published work, can be summarised as follows. We have focused on an integrated approach that couples state of the art biophysically-detailed models of electrophysiology and contractility, usually investigated separately or with simplified formulations [24], [25]. We adopted a dynamic formulation taking into account inertial forces, as it has been shown their relevance in cardiac electromechanics [22]. We implemented a pre-stress model to achieve diastolic pressure, allowing the simulation of the four phases of the cardiac cycle, together with an appropriate representation of epicardial mechanical boundary conditions.

From the simulated pressure–volume relationship which recapitulates the different phases of the cardiac cycle, physiologically-relevant mechanical biomarkers (EF, ESP, LFS, and WT) were extracted in a total of 200 simulations (i.e. 100 for one-at-a-time and 100 for LHS-based approaches). Simulation results using LHS were then used to compute the global sensitivities of relevant mechanical parameters, which have not been rigorously studied in human before. While recent studies [29] have addressed how electrophysiological biomarkers are affected by the incorporation of mechanics into a model of cardiac electrical diffusion under the presence of parameter uncertainty, so far the analysis has been restricted to idealised three-dimensional bars. This has implied that physiologically relevant mechanical biomarkers such as EF could not be assessed. Further comparisons will need to be performed once a more abundant body of other studies are available.

Our sensitivity analysis of model parameters with respect to mechanical biomarkers computed with an ellipsoidal model of the left ventricle of clinically-relevant volume [33] has shown that simulations with the human model and parameter variations yield values for the mechanical biomarkers ranging from healthy to disease scenarios. Varying a single model parameter in the simulations yielded variations in mechanical biomarkers that were nonlinear but monotonic (Fig. 3.2). This justified the use of LHS sampling-based PRCC for the evaluation of sensitivities. It is important to point out the strong influence of on EF, which is attributable to the fact that higher leads to higher active tension, causing higher contraction of the ventricles and thus higher EF. In fact, by looking at Fig. 3.2, one can see that variation of alone yields physiological values of EF. Force production in a cardiac electromechanical model is mainly regulated by the active tension model, which needs a calcium input provided by the cell electrophysiology model; the human-based active tension model used in this study [9] was originally parametrised using very high values of calcium (coming from [77] instead of from a cell electrophysiology model, such as [6], [78]), but the values of calcium provided by [6] are much lower, which is translated into a decrease in force output. Therefore, scaling the active tension up to 7- or 8-fold was needed in order to reach physiological values of EF. An alternative would have been to re-parametrise the coupled cell electrophysiology and active tension model in a previous step; however, a simpler yet as effective approach is scaling , which is that considered in this study. This highlights the importance of calibrating the active tension model to a physiological healthy calcium transient, particularly when coupling it to the cell electrophysiology model.

When considering a global sensitivity analysis approach using LHS, a wide range of variability in mechanical biomarkers was obtained, including values lying on healthy and on diseased ranges (Fig. 3.4). Some of the LHS cases were also lying on the healthy ranges for all of the considered biomarkers. This means that such combinations of model parameters could be used as the baseline values for future studies involving this electromechanical model. The comparison with clinical ranges yielded some important observations. The simulated WT values obtained in this study lie on the lower physiological range (Fig. 3.4), which could potentially be due to the following reasons. The formulation of the computational model may need to be updated to include an active strain formulation with orthotropic activation as argued by [23], [79], in lieu of the active stress formulation used in this study, possibly even with heterogeneous orthotropic activation across the cardiac wall [23]. In our simulations, sheet directions are constant across the cardiac wall, with a sheetlet angle of 0 (i.e. parallel to the direction normal to the endocardial/epicardial surfaces). It has been suggested that while fibre angles vary linearly throughout the cardiac wall, sheetlet angles exhibit a non-monotonic variation, with the following approximate behaviour: going from a certain value at the epicardium, to in the middle of the cardiac wall, to again at the endocardium [80]. Sheetlet angles different than zero might also be needed as results may also depend on their distribution [27], [81]. Furthermore, possibly higher values of alone are needed to achieve higher, more physiological values of WT, which were not used in this study as it could compromise the positive-definiteness of the tangent operator of the Newton–Raphson scheme (the preconditioned conjugate gradients solver used for the linear algebraic systems resulting from the spatial discretisation of the solid mechanics subsystem did not converge for some cases with high , even when using very small time steps). From a purely geometrical point of view, the thickness of the considered ventricle was 1.2 cm, which is on the upper range of human left ventricular wall thickness [82]; this may have also had an impact on WT values.

With respect to the sensitivities resulting from the global approach, the LHS-based approach yields tendencies in the sensitivities similar to those in the one-at-a-time results, in the sense that they generally have the same correlation signs. The PRCC suggest that EF, ESP and LFS are mainly dominated by one or two parameters each: EF is dominated by /, ESP is dominated by /, and LFS is dominated by / (Fig. 3.5). This means that extra effort could be placed into further experimentally characterising these parameters in order to obtain more accurate probability distributions of the considered mechanical biomarkers. On the contrary, WT does not seem to be clearly dominated by any of the considered input parameters, which suggests that it is likely affected by many different factors, such as the distribution of sheet and fibre directions, numerical formulation of contraction (i.e. active stress or active strain), or the parameter(s) that control the maximum level of contraction, among others. The distributions of EF, ESP and LFS with respect to their most influential parameters according to the PRCC values in Fig. 3.5 are seen in Fig. 3.6. One can see the shown correlations in the plots. Furthermore, the classification of cases according to their EF values (into “diseased” and “healthy”) helps in seeing that there is no direct correlation between EF and any of the other mechanical biomarkers, which means that these need to be assessed on a case-by-case basis instead of be inferred from EF values.

Some limitations of this study are worth remarking. The present study includes a sensitivity analysis of the effect of key parameters on mechanical biomarkers, but not all parameters are considered to enable study feasibility. The results can inform the design of further studies focused on geometry, mathematical formulation or boundary conditions, as these can also be important. An active stress formulation was used mainly because it is less invasive than an active strain formulation, and also because most active tension models are tailored to output active stress in the fibre direction, possibly due to experimental evidence at the cell scale [9], [83]. Our results show that simulations do not reproduce the whole range of WT reported in clinical studies. A simplified method was used, in which the thickness was calculated from two points that were initially found at the same height with respect to the base. Needless to say, during contraction this relative position is likely to change and thus introduce some bias in the thickness calculation. Additionally, experimental calculation of WT significantly depends on the evaluation method, with linear methods introducing error because of MRI slice positioning. If individual variation is considered, it can be seen that there is a substantial sparsity in results; for example, evaluation of WT on 9 healthy subjects led to values ranging from 18% to 100% [67]. It is likely that WT values are affected by the choice of volumetric term in the passive model for cardiac tissue and thus further research in this respect is warranted [84], [85]. The initiation of the system was done at the cellular level, and then an initiation phase was used to bring the geometry to an end-diastolic state; however, it is likely that the geometry also needs to be brought to periodic steady-state by repeating a cardiac cycle a few times. Needless to say, performing such thorough initiation needs huge computational resources and would indeed limit the number of LHS samples that could be performed. The geometry used in this study is an ellipsoidal left ventricular geometry, and therefore does not represent a biventricular geometry. Future studies can evaluate the role of anatomical variability in mechanical properties. In terms of the mechano-electric feedback, [86] suggested that stretch-dependent changes in capacitance result in changes in conduction velocity, and are important in order to retrieve an appropriate slow-down in conduction (such as that observed in rabbit ventricles [87]). In our study, we incorporate strain-dependent conductivity tensors as the main mechano-electric feedback mechanism, which is comparable in formulation to that of stretch-dependent space constants, and indeed allowed for conduction slow-down, as explained in Section 2.3.2.