Embedded Computational Heart Model for External Ventricular Assist Device Investigations

Purpose External cardiac assist devices are based on a promising and simple concept for treating heart failure, but they are surprisingly difficult to design. Thus, a structured approach combining experiments with computer-based optimization is essential. The latter provides the motivation for the work presented in this paper. Methods We present a computational modeling framework for realistic representation of the heart’s tissue structure, electrophysiology and actuation. The passive heart tissue is described by a nonlinear anisotropic material law, considering fiber and sheetlet directions. For muscle contraction, an orthotropic active-strain model is employed, initiated by a periodically propagating electrical potential. The model allows for boundary conditions at the epicardium accounting for external assist devices, and it is coupled to a circulation network providing appropriate pressure boundary conditions inside the ventricles. Results Simulated results from an unsupported healthy and a pathological heart model are presented and reproduce accurate deformations compared to phenomenological measurements. Moreover, cardiac output and ventricular pressure signals are in good agreement too. By investigating the impact of applying an exemplary external actuation to the pathological heart model, it shows that cardiac patches can restore a healthy blood flow. Conclusion We demonstrate that the devised computational modeling framework is capable of predicting characteristic trends (e.g. apex shortening, wall thickening and apex twisting) of a healthy heart, and that it can be used to study pathological hearts and external activation thereof.


INTRODUCTION
Advanced heart failure is currently treated most successfully by heart transplantation. However, donor organ availability is very limited. Thus, over the past decades, ventricular assist devices (VADs) have become the standard alternative medical therapy. 25 Usually, turbodynamic flow pumps support or replace left ventricular functions; and in severe cases, both heart chambers. Some patients receive a VAD only as a bridge to heart transplant, whereas for some patients it is their destination therapy. Over the years, many products have been invented; some work in a pulsatile manner, but most of them provide a more or less continuous flow. What all of the nowadays implanted devices have in common is that they remove blood from a heart chamber and deliver it into the ascending great vessels. Downsides of existing VADs are the surgical opening of anatomical structures of a closed system, direct contact with the blood stream and unnatural temporal actuation. To avoid these shortcomings, other concepts have been suggested, e.g. to use an inflating balloon in the right ventricle 15 or the intrapericardial cavity, 11,17,21 two contracting rings around the heart, 33 inflatable diaphragms on the heart's surface, 20,31,48 curvature inverting devices, 27,28,30 a twisting system 7 and the criscone device. 28 Here, an external mechanical actuation concept, inspired by open cardiac massage and previous work, 16,41 is considered. Similar as the surgeon's palms, two opposed force patches periodically compress the heart and thus enhance the circulation. Before such devices can be used even for short term applications, such as during heart surgeries, numerous challenges need to be addressed. In addition to the difficulties in achieving desired flow rates and pressure differences, the consequences of attaching force patches directly on the epicardium are unclear. As illustrated in Fig. 1, epicardial force patches may lead to undesired deformations of the papillary muscles and heart valves, tissue fatigue, squeezing of the coronary arteries underneath the patches, and critical stress peaks. Understanding how these factors influence cardiac function are important to assess the efficacy of the device.
The main focus of this paper lies on the development of a computational framework suited to investigate external application of force patches while gaining a better understanding of cardiac mechanics. Investigating the impact of VADs on the heart function using experiments can be difficult, if not impossible, and expensive. On the other hand, computer simulations of the heart allow to study a variety of pathologies more easily. For example, one can use computational heart models reflecting coronary artery diseases, electrophysiological abnormalities or pressure-volume anomalies. Moreover, a computational model allows for optimization studies; objectives are e.g. cardiac output, overall contraction pattern, stress peaks or the right ventricle's vulnerability. By adjusting force patches as functions of space and time, the desired flow rates and pressures can be imposed. The center piece of the computational framework presented in this paper is based on the finite element library LifeV, which implements a detailed model of the heart structure accounting for the highly anisotropic muscle, mechanical and electrical behavior. A Holzapfel-Ogden model 8,14 is employed, and required fiber-and sheetlet directions are generated according to Refs. 43,52. Rigid external force patches are modeled using time-dependent non-homogeneous Dirichlet boundary conditions. Muscle activation, which is modeled by an orthotropic active strain formulation 2,29,38,43,44 and causes the cardiomyocytes to contract, is triggered by an electrical wave. 3 Although similar frameworks have been published, 9,12,13,34,37,47,50 activation has not been described by active strain nor has the effect of force patches been studied. Finally, in order to calculate blood pressure in the two ventricles, which is required for properly predicting deformations, the finite element model is coupled to a simplified hydraulic network accounting for pulmonary and systemic circulations. 6,46 For the heart valves a novel pressure gradient dependent model is proposed.
The paper is organized as follows. In Embedded Heart Model and Governing Equations section we describe the governing equations of the individual submodels for electrophysiology, activation, structure mechanics and the lumped parameter circulation. The computational framework, that is, the numerical solution algorithms and coupling thereof, is presented in Methods: Computational Framework section. Numerical results are shown in Validation and Results section; first of a healthy heart and then of one suffering from coronary heart disease. Moreover, simulation results with the pathological heart model ð Þ and f 2 X; t ð Þ attached to the outer heart surface for actuation. The three-dimensional human heart model is from Ref. [1]. assisted by two force patches and a parameter sensitivity study are presented. The paper closes with a discussion and conclusions in Discussion and Conclusion section.

Mechanical Response of the Myocardium
Let X 0 & R 3 be the undeformed structure configuration of the heart with coordinates denoted by X, and x being the current position of the deformed body X & R 3 . The undeformed boundary @X 0 is divided into C D with Dirichlet boundary conditions, C N with Neumann boundary conditions and C R with Robin boundary conditions. Inertial and body forces (due to gravity) are negligible compared to other forces 51 and thus are ignored here. The resulting quasi-static boundary value problem reads where u ¼ x À X is the displacement, P the first Piola-Kirchhoff stress tensor, m the outward surface normal in the reference configuration and u D and p N ¼ p m Dirichlet and Neumann boundary conditions, respectively. Note that r is operated with respect to X in the reference configuration. The myocardium is considered hyperelastic, and therefore its mechanical response can be fully characterized by a scalar-valued strain-energy function W. We consider an orthotropic description of the myocardium given by the Holzapfel-Ogden constitutive law, 14 which leads to for the passive strain-energy function with the four invariants which are functions of the right Cauchy-Green tensor C ¼ F T F and the fiber and sheetlet directions f and s, respectively. The constants a, b, a i and b i i 2 ff; s; fsg ð Þwere found by fitting experimental data. 5 The stress tensor is defined by the relation P ¼ @W @F , in which F ¼ r x ¼ I þ r u is the deformation gradient tensor.
Because numerical experiments have shown that the active-strain formulation is capable of reproducing appropriate contraction patterns, we follow such an approach. 29,43 Assuming the deformation gradient admits a multiplicative decomposition, F ¼ F E F A , in which F E is the passive, elastic deformation and F A is an active deformation described by the expression in which I is the identity tensor. The active-strain functions c f , c s and c n provide information on the contractile strains. Notice that f 0 , s 0 and n 0 are orthogonal directions in the reference configuration; see Fig. 2. Given the decomposition of the deformation gradient tensor into an active and elastic component, the active strain formulation assumes that the strain energy can be expressed by Eq. (2.2), where the invariants I 4f , I 4s , I 8fs and I 1 have to be replaced by 4Þ and ð2:6Þ ð2:7Þ respectively.

Active Strain Dynamics
In Eq. (2.3), the active part of the deformation gradient F A is determined by the three active-strain functions c f , c s and c n . Appropriate definitions for these have been derived in Ref. [43]. Thus, we assume the active strain in fiber direction c f satisfies where c is the calcium concentration coming from electrophysiology, I 4f is the fiber elongation from mechanics, l A and a A are constants and R FL is the sarcomere force-length relationship, which can be found in Ref. [43]. The two other active strain functions c s and c n are determined by imposing volume conservation, i.e., detðF A Þ ¼ 1, and are functions of c f and an orthotropic parameter j, i.e.,

Electrophysiology
Before any contraction of the heart muscle starts, the sinoatrial node (SA-node) initiates a wave of electrical depolarization propagating from the right atrium to the left atrium and to the atrioventricular node (AVnode). The AV-node delays the signal by approximately 0.12 s, before it propagates through the ventricular muscles. From the AV node, the electric signal travels through the fast cardiac conduction system to quickly activate the endorcardial surface in multiple locations and ensure a uniform ventricular contraction. To reproduce the signal after the AV-node, we use the monodomain equations combined with the minimal model for human myocardial action potential. 3 In the reference domain X 0 the monodomain equations read in which V is the transmembrane potential, the vector w represents the gating variables, c is the vector of ionic concentrations, I app is an externally applied density current, I ion and r link the two equations, G is the spatial (in the deformed configuration) anisotropic conductivity tensor, depending on the fiber and sheetlet directions f 0 and s 0 , v and C m are constants, F is the deformation gradient and J ¼ detðFÞ is its determinant.

Fiber and Sheetlet Directions
While DT-MRI data for reconstructing fiber fields are available, information on the collagen sheets is usually unavailable. Here, both fields are reconstructed using a rule-based approach. 8,43,52 Assuming the sheetlet field s 0 to be divergence-and curl-free, s 0 can be derived from a scalar potential / as s 0 ¼ r /. Taking the divergence of the equation, a Laplaceequation for the potential / is obtained, that is, with Dirichlet boundary conditions / ¼ g on C D and Neumann boundary conditions @/ @m ¼ h on C N . Having solved Eq. (2.11) inside the undeformed myocardium for /, the sheetlet vector is s 0 ¼ r /=jjr /jj and points from the endocardium toward the epicardium. Next, consider the left ventricular centerline k, pointing from the apex to the base. Its projection on the plane orthogonal to s 0 is Figure 2 shows an enlarged section of the heart tissue with fiber-and sheeltlet vectors f and s, respectively; n ¼ f Â s. Note that the direction of the fibers depend on the layer. A preliminary, horizontal fiber field is orthogonal to both s 0 as well as k p and readsf 0 More details about building R s 0 can be found in Ref. [43], where it also is shown that computationally obtained fiber fields are in good agreement with dissections. 45

Circulation
The body circulation is represented by the circuit depicted in Fig. 3, which is adapted from Ref. [6,19,46]. The driving force of the circulation is exerted by both the right and left heart chambers. This force is a boundary condition in Eq. (2.1) and can be divided into pressure, viscous and inertial forces. Because forces due to inertia and wall shear stress inside the heart are much smaller than the pressure force, 36 we neglect their contribution. Moreover, because in both ventricles the dynamic blood pressure is much smaller than the static pressure variations during a heart beat, we assumed a spatially constant pressure field. 9,10,24 Ventricles are drawn as variable capacitors and stand for the three-dimensional finite element model of the heart. Both the pulmonary and systemic circulations include three resistances, three capacitances and an inertance. Each of the heart valves is modeled by a resistance and a diode, and both atria are modeled as a capacitance. The dependent variables are pressure values p at the nodes, from which the flow rates can be calculated. Lumped parameter models as used here obey Kirchhoff's circuit laws and are by construction conservative. Thus, the sum of all flow rates q iÀj from node i to its connected nodes j is zero, i.e., X j q iÀj ¼ 0: ð2:12Þ The resistance R between two nodes i and j is the ratio of pressure difference p i À p j and flow rate q iÀj , which results in Inertia in both aorta and pulmonary artery is modeled by the inertences L sa and L pa , respectively, i.e., the corresponding pressure difference is a function of temporal flow rate change according to Blood vessel elasticity is modeled by compliances C between a certain node i and the ground, from which one obtains To account for the heart valves, a new diode model is proposed here. In an ideal setting, diodes are open for a positive pressure gradient and closed otherwise. However, heart valves are not ideal diodes and therefore it is suggested here to model them as that is, as a compartment consisting of a resistance R and a diode D ¼ D p i ; p j ; t À Á 2 0; 1 ½ . In other words, a valve is described as a temporally varying resistance, and the three simple intermediate differential equations Note, however, that the focus of this study is not on valve models. The resistance of the outflow valves additionally depends on the ventricular pressure values. 6 For example, the resistance of the aortic valve is R ao ¼ k lv p lv , where k lv is a constant parameter. Moreover, when the heart valves open, ventricular volume decreases and thus the contraction is no more considered isometric. An instantaneous ventricular pressure p lv;inst is obtained as the difference between the isometric pressure p lv and viscous losses, 6 namely, p lv;inst ¼ p lv À R ao q ao . The same is done for the instantaneous pressure in the right ventricle.
The link between the circuit and the mechanics is established via pressure and volume of the two nodes representing the ventricles and will be explained in more detail later. Changes of left and right ventricular volumes are described as dV l dt ¼ q mi À q ao and ð2:21Þ where q mi , q ao , q tr and q pu denote the flow rates through the four heart valves.

METHODS: COMPUTATIONAL FRAMEWORK
The overall framework has been implemented by using the parallel finite element library LifeV (http:// www.lifev.org). Initial data and model parameters can be found in the Tables 3, 4 and 5.

Finite Element Discretization for Electromechanics
Equation (2.1) is solved with a Galerkin finite element method. Tetrahedral elements with quadratic FIGURE 3. Lumped parameter model of the blood flow in the human body. The pulmonary circulation, where blood is oxygenated, is drawn at the top. At the bottom, the systemic circulation is shown, which delivers the oxygen to the body's tissue. This model is an adaption from Ref. [6,19]. The meaning of the symbols is shown in Table 1.
shape functions are employed for all electromechanical fields and a parallel version of the generalized minimal residual method (GMRES) is used to deal with the resulting system of linear equations. Solution of the non-linear problem is achieved with Newton's method.
Incompressibility of the heart tissue is enforced by decomposing the strain energy function into an isochoric and a volumetric part, i.e., where the volumetric penalty function W V I 3 ð Þ is added to Eq. (2.2) and only depends on the invariant used in this work. For the term W A I 1 ; I 4;f ; I 4;s ; I 8;fs À Á in Eq. (3.1) to be isochoric, all four invariants in the argument list must be isochoric as well. However, significant changes only concern the first invariant I 1 , 42 which is replaced in Eq. (2.2) by Coefficients for the isochoric part can be found in Table 2.
Active strain evolution in Eq. (2.8) is solved explicitly in time and its dependency on the prestretch is removed by a Taylor-expansion, 43 i.e., Active strain measures c f , c s and c n are spatially discretized in the same manner as the structural mechanics; however, they are scalar fields. Time discretization of Eq. (3.2) is done by a forward Euler scheme. Active strain model coefficients are a A ¼ À6:25 lM À2 , l A ¼ 5000 slM À2 , c 0 ¼ 0:2155 and j ¼ 4. Equation (2.10) is solved by a semi-explicit method. 38 Therefore, the gating variables w are up-dated explicitly by a forward Euler scheme first. The transmembrane potential V is determined in a second step by ionic current interpolation (ICI). 32 Spatial discretization is the same as for activation and the mechanics. Ionic cell model (Minimal Model) parameters are found in Ref. [3] for the myocardium and monodomain model parameters as well as electrophysiology initial data in Table 2.

Ordinary Differential Equations for Circulatory Model
The hydraulic circuit depicted in Fig. 3 is described by Eqs. (2.12)- (2.20). Note that the circuit model is not an integral part of the finite element library LifeV, i.e., it had to be implemented separately. Coupling with the heart model is described in ''Coupling Mechanics and Circulation'' section. By combining Eqs. (2.12)-(2.16), the system of ordinary differential equations is obtained, where the vector of unknowns, consists of a vector b, which stores nodal pressure values, and of a vector q storing flow rates between neighboring nodes. Time discretization is done by a backward Euler scheme leading to with the unknown solution x nþ1 at the new time t nþ1 , the previously computed solution x n at the old time t n and the time step size s c ¼ t nþ1 À t n . Considering Eqs. (2.17)-(2.20), which determine whether a heart valve is open or closed, it can be seen that the valve state D is a function of the pressure difference across the heart valve. This implies that matrices M and A as well as the vector u depend on x. In order to retain an implicit scheme, a fix-point iteration algorithm is employed, which reads where k is the step number, . At the initial step k ¼ 0, x nþ1;0 is set to x n . Fix-point iterations are carried out as long as jjx nþ1;kþ1 À x nþ1;k jj 2 jjx nþ1;k jj 2 >e c ; ð3:6Þ where e c ¼ 10 À6 is the desired relative change from one to the next iteration step. Parameters for various electric/hydraulic elements can be found in Table 3.

BIOMEDICAL ENGINEERING SOCIETY
At the beginning of every iteration step of Eq.  Table 4.
Once D is determined for all four heart valves at step k, M k , A k and u k are assembled and the linear system is solved for x nþ1;kþ1 . Initial conditions for the four valve states, pressure and flow rate values can be found in Table 5.

Coupling Mechanics and Circulation
Electrophysiology, activation, mechanics and circulation require different temporal resolutions. To capture the fast dynamics of the steep propagation front a small timestep for electrophysiology and activation is used (s e ¼ 0:05 ms), whereas for the mechanics and the circulation larger timesteps are possible (s m ¼ s c ¼ 1 ms). Local load steps for the mechanics are performed in between two timesteps (s ls ¼ 0:5 ms). The temporal evolution is outlined by Algorithm 3.1.
Coupling between the mechanics and the body circulation is achieved via ventricular pressure and volume. 19 Therefore, the volume of the finite element mechanics model is conserved as well. Based on our experience, only an implicit coupling is stable during all phases of a heartbeat. An implicit scheme is obtained by applying the same pressure boundary conditions p ¼ ½p lv p rv T inside the left and right ventricles on both the mechanical model and the body circulation. By minimizing jjRjj 1 ¼ jjV fe À V circ jj 1 where V ¼ ½V lv V rv T contains both ventricular volumes, the pressure boundary conditions are updated through Newton's method. Thus, the Jacobian matrix reads where the derivatives of V fe and V circ are determined by pressure perturbation. For the results presented here, we have altered pressure values by 10 À3 mmHg. Perturbation of the mechanical model is based on the previously computed linearization of PðuÞ. Eventually, Newton iterations with step k are performed as long as jjRjj 1 >e n is larger than a defined maximum error of e n ¼ 10 À5 mL .   Algorithm 3.1 explains the temporal solution algorithm structure of the framework.

VALIDATION AND RESULTS
All simulation results presented here have been performed on the Euler cluster at ETH Zurich. They were executed on eight Hewlett-Packard BL460c Gen8 nodes, each with 128 GB Ram and two 12-core Intel Xeon E52697v2 processors (2.7 Ghz nominal, 3-3.5 GHz peak); that is, in total on 192 cores. With 160 k elements it took 16.43 hours to simulate one heart beat of 800 ms.

Validation with Healthy Heart Model
Validation of the whole framework can be done in different ways, for example by comparing contraction patterns from imaging data or by comparing with electrophysiological and flow measurements. However, the range of measurements among individuals is large and therefore our framework shall reproduce common observations, e.g. by how much the distance between base and apex shortens, the wall thickens and the apex twists.
Most heart geometry measurements are noisy and vary greatly in size. Thus, this framework uses idealized geometrical data, which are very detailed. 1 The geometry is prepared for simulations according to Fig. 4. As seen in the first state of Fig. 4, raw heart data are detailed and contain many features not required for our studies. From the first to the second state, features such as the great vessels, both the atria, coronary arteries, heart valves and papillary muscles are removed. The final spatial domain, consisting only of the strong ventricular muscles, can be prepared for finite element discretization; it is divided into tetrahedra using GMSH (www.gmsh.info). Three differently sized meshes with 160 k, 250 k and 500 k elements are deployed. While the solid mechanics solutions obtained with 160 k elements hardly differ from those obtained with higher resolution, the electrophysiological problem remains under-resolved with 160 k elements. However, at low resolution the coefficients r f and r s of the conductivity tensor G can be adjusted such that the propagation speed of the electric signal is in very good agreement with reality. Therefore, the 160 k mesh was used for the subsequently presented results. Transition from the mesh to fiber and sheetlet fields (state 3 to 4 in Fig. 4) is based on the algorithm described in ''Fiber and Sheetlet Directions'' section. Fiber angles vary between À 60°to + 60°from both endocardia to the septum and epicardium.
Insulation boundary conditions are applied to the electrical potential, the ionic model and the activation. Boundary conditions for the mechanics are set as follows. On both endocardial surfaces pressure boundary conditions, i.e., P Á m ¼ p m, are applied, where p comes from coupling with the circulation. At the base, around the four heart valve openings, where the great vessels are attached to the heart, Robin boundary conditions, i.e., P Á m ¼ a u with a ¼ a b ¼ 750 mmHg cm À1 , are  enforced. In order to imitate surrounding tissue such as the lungs, Robin boundary conditions with a ¼ a e ¼ 0:75 mmHg cm À1 are employed on the epicardium as well.
Before the actual simulation starts, the heart is preloaded with initial pressure values, that is, 5:7 mmHg inside the left ventricle and 2:4 mmHg in the right one. This preloading requires about 50 load steps causing a volume increase of 25 mL and 23 mL in the left and right ventricles, respectively. Every heartbeat lasts 800 ms and is initiated by an electric wave traveling through the domain; it spreads within 60 ms through the entire ventricles. Subsequently, we initiate the signal in a simplified fashion by an applied electric current I app (see Eq. (2.10)), which is applied in the lower third of both endocardial surfaces. Figure 5 presents the evolution of the potential V, the normalized calcium concentration c and the activation level c f at a spot inside the left myocardium at half height. At other coordinates the curves would qualitatively be similar, but shifted in time. Results from the lumped parameter model are shown at the top of Fig. 6. Depending on the initial conditions of the model, the results change from one heartbeat to the next one, until they have converged. In our case, this takes usually two to three heartbeats. On the left side of Fig. 6, temporal pressure and flow rate evolution indicates a qualitatively physiological blood flow behavior in comparison to measurements. 4 In addition, flow rate overshoots are visible as a result of the valve modeling and the phase of refilling through the mitral valve is rather short. On the right side of the figure, p-V loops of both ventricles are drawn. The total power output of the heart model is about 1.3 W, which is in good agreement with measurements. 18 In order to validate the contraction pattern of the heart model with phenomenological measurements, Fig. 6 displays a sequence of slices through the heart such that both ventricles are seen from inside. The color indicates the calcium concentration spreading through the domain. It can be seen how the heart shortens and how the wall thickness increases, as reported by others. 23,39 Three phenomenological values, which are appropriate for validation of the ven-tricular contraction pattern, are the base-apex shortening, wall thickening and base-apex twist.
The plot on the left of Fig. 7 depicts apex shortening, apex twist and wall thickening as function of time. The right side of Fig. 7 exhibits a cut through the heart before and during contraction, highlighting apex shortening and wall thickening.

Results with Pathological Heart Model
We consider a pathological heart suffering from coronary heart disease. Thus, the myocardium's ability to contract during systole is locally reduced, and Eq. (2.8) is adjusted by setting the activation level c f to zero in a specified region. Figure 8a depicts in red the damaged region of 3 cm radius imposed on the anterior left ventricle. Figure 8b shows a cut through both heart chambers and the myocardial infarction is marked dark red. Black lines indicate the perimeter of a systolic healthy heart. It can be observed how the affected tissue bends outward.
Myocardial infarction not only changes the contraction pattern, but also the cardiac output. Figure 8c makes this change visible by comparing the p-V loops of a healthy heart with those of diseased ones. Since the infarction region is in the left ventricle, the right ventricle's stroke work hardly changes for all three diseased cases. Left ventricular stroke work for an infarction radius of 3.5 cm decreases by almost 50%, whereas for a radius of 2.5 cm, the reduction is only 20%.

Results with External Actuation
Force patches (see Fig. 1) are modeled by Dirichlet boundary conditions for Eq. (2.1). For simplicity, two patches with spatially uniform displacement vectors FIGURE 7. Left: temporal evolution of apex shortening, wall thickening and apex twist, which peak at about 15 %, 45 % and 12°, respectively. Right: cut through the heart, highlighting apex shortening and wall thickening.
ð Þ and f 2 X; t ð Þ are introduced on the circular boundaries C D 1 and C D 2 , respectively, which are defined withing the anterior and posterior epicardial surfaces of the left ventricle. Both patches are rigid (i.e., they contract as one), have a radius of r ¼ 2 cm and their unit motion vectors n point toward each other, that is, n 1 ¼ Àn 2 . Their displacement vectors f i i 2 f1; 2g ð Þallow for a smooth variation between timesteps and are given by wheref i is the peak displacement magnitude,t i the peak time, T i the impact duration and n i the unit vector in motion direction. Both patches move synchronously such that T ¼ T 1 ¼ T 2 andt ¼t 1 ¼t 2 . Their peak displacement magnitudes will be identical as well, that is,f 1 ¼f 2 . For the two patches presented above, parameterst and T are chosen such that the patch motion is synchronous to natural heart activa- tion, i.e.,t ¼ 250 ms and T ¼ 300 ms. Figure 9 shows the patch displacement jfjðtÞ for a heart beat duration of 800 ms,t ¼ 250 ms, T ¼ 300 ms andf ¼ 0:9 cm. Further, Fig. 10 visualizes the heart's deformation when the patches meet their peak displacements. In this demonstration, Dirichlet boundary conditions for patches are only applied in perpendicular directions to the centerline of the left ventricle k, i.e., sliding of the tissue underneath the patches is allowed in direction of k. Therefore the heart muscle can still properly contract along k. Under these conditions, to regain a healthy stroke work, we need to impose a peak displacement of magnitudef ¼ 0:9 cm; the corresponding p-V loops are shown at the bottom of Fig. 11. Left ventricular results in Fig. 11 are promising regarding stroke work, i.e., a healthy p-V-loop is recovered. However, an unexpected pressure increase can be detected in the right ventricle. It emerges from the heart's geometry, where the right ventricle is wrapped to some extent around the left one. The top of Fig. 11 presents the deformation of the assisted case during one heart beat. In contrast to Fig. 6, the deformation pattern is most different in the area where patches are attached.

VAD Sensitivity Analysis
In this section we vary the four patch parametersf, t, T and r and analyze their influence on the p-V loop. From the top of Fig. 12 it can be seen, how the p-V-Loop changes with the displacement magnitude, i.e., forf ¼ 0:7, 0:8, 0:9, 1:0 and 1:1 cm. Despite the small changes, a considerable difference in stroke work and pressure values can be observed. Next, the sensitivity of the p-V-loop is analyzed by choosing different peak timest (see bottom of Fig. 12). While the left heart chamber is not much affected, the right ventricular peak pressure and stroke volume decrease with increasingt. Thus, overloading the right ventricle may be reduced by adjusting the peak timet appropriately. Here, the patches on the epicardium have circular shape, which may not be optimal. Nevertheless, the patch radius is varied in order to study the influence on the p-V-loop. The top left plot of Fig. 13 demonstrates that a larger patch radius yields more stroke work by the left ventricle. A similar response is observed for the right ventricle.
Eventually, the period T in the temporal patch activation function, i.e., in Eq. (4.1), is changed. From the bottom of Fig. 13 hardly any effect can be observed for the left ventricle, whereas the right peak pressure slightly increases for longer durations T.

DISCUSSION AND CONCLUSION
The objective of this paper is to propose a modeling framework for the human heart, which allows for investigations of external actuation via force patches attached to the epicardium. The framework consists of a detailed FEM model for the structural mechanics, electrophysiology and actuation, which is coupled to a simplified hydraulic network model representing pulmonary and systemic blood flow circulations. For the latter we introduced a novel valve model, which accounts for a pressure gradient dependent dynamics with typical opening and closing times.
A nonlinear anisotropic Holzapfel-Ogden material law accounting for fiber and sheetlet orientations, and an active strain model responding to the signal of a  sophisticated electrophysiological model, lead to very realistic deformations. This has been demonstrated for a healthy heart model. We showed that the computed base-apex shortening, apex-twist and wall thickening are within the expected ranges. The twist first decreases by about 3°and as soon as the outflow valves open, it increases to about 12°. This is in very good agreement with measurements. 22,40,49 Wall thickening first decreases by about 12%, before it increases up to 45%, as reported similarly in Ref. [23,39]. The apex shortens from the beginning of activation and the shortening peaks at 15%, which is consistent with observations. 23 Moreover, predicted pressure and flow rate evolutions are in good agreement with measurements. Of course, this sort of validation is only conclusive to some extent. However, if a model matches these observations, it is likely to reproduce a good overall contraction pattern.
Simulation studies with a pathological model of a heart suffering from myocardial infarction show (as expected) a dramatic reduction of the induced flow rates and significantly different deformations of the heart muscle. Computational results with external actuation using two opposed force patches applied on the epicardium of the diseased heart predict left ventricular pressure and volume changes which are in very good agreement with those of the healthy heart. The right ventricular pressure, however, rises above the desired peak value with this simple actuation scenario. Sensitivity studies with different patch sizes, peak displacements of the patches, peak displacement times, and different impact durations show that the right ventricular power output can be controlled. The gained knowledge can be used to adjust such patch configurations. It is emphasized that a healthy p-V-loop is not the only objective. For example, tissue stress and fatigue, as well as valve and papillary muscle deforma- FIGURE 11. Top: left ventricular deformation during one heartbeat of a diseased model, which is assisted by two opposite patches. The contraction pattern is different compared to healthy conditions (Fig. 6), particularly in the patch regions. The calcium spread is visualized by the color. Bottom: p-V loops for a healthy, a pathological and an assisted pathological heart are presented. The diseased heart has an infarction radius of 3 cm. In order to obtain the same power output as the healthy heart, the pathological heart has to be assisted with a peak displacement of 0.9 cm for both patches.
tion have to be taken into account. And of course it is crucial that the interface between patches and epicardium does not damage the tissue. However, in line with the focus of the paper, this study only serves to demonstrate the potential of the presented computational framework to perform optimization studies, and no attempt has been made yet to find a better combination of parameters, patch geometry and motion. Nevertheless, already the presented results are encouraging for future external VAD developments, and they demonstrate the usefulness of this computational model.
The main limitations of the modeling framework are uncertainties regarding geometry, distribution of material properties and of other parameters. They have been chosen based to the best of our knowledge relying on literature data. Especially among pathological hearts the variation is huge, and it is a separate challenge to retrieve complete descriptions. All presented simulation studies used the same basic heart geometry, and only parameters regarding patches and their actuation have been varied. Further sensitivity studies would be of interest and are planed in the future. There is a demand of personalizing the modeling of human pathophysiology, 26,35 and it is important to notice that the devised framework can be applied to investigate patient-specific or animal hearts, provided geometry, material properties and all other relevant parameters are available. Another limitation is the simplified flow model. Therefore it is not directly possible to take into account effects of damaged heart valves or the impact of force patches on them. At this point, a further shortcoming is the force patch boundary condition, which only accounts for the force component normal to the epicardium, i.e., tangential forces are ignored. Further, the problem that patches may slide off the slippery surface cannot be addressed, interactions between patches and coronary vessels are ignored, and no conclusions regarding long term tissue damage introduced by the force patches are possible. Also, at this point, independent actuation is implemented; in the future synchronization with the electrical pulses is envisioned.
Despite several limitations, the framework is of significant value for future investigations and optimization of external VADs with force patches, as demonstrated by the computational studies in this paper. In addition, it can be used to study different disease scenarios and/or various assist device designs. On one hand, it can be expected that it will help to reduce the number of required in vivo experiments, and on the other hand it allows for investigations which are not even possible in vivo. Due to its flexibility, the presented framework has the potential to become an integral part of patient-specific applications. Extraction of geometrical and material information, estimation of further parameters, and uncertainty assessment would be other components of such a much bigger framework.
Potential improvements of the presented computational framework include local modifications of the conductivity tensor in pathologic regions, a completely incompressible treatment of the elasticity by a mixed formulation approach, and atrial contraction by an additional compliance in the hydraulic circuit. Further, a resolved flow field in the two ventricles would be of interest in some cases, e.g. if heart valve functions become part of the investigations. And after all, there exists potential to improve the computational efficiency, which is of particular interest for parameter and optimization studies.
Regarding device development, further sensitivity studies would be of great interest, and optimization not only of patch location, patch displacement peak, displacement peak time and patch size, but also of patch shape and different actuation patterns. In addition to p-V curves, other objectives will be considered for optimization. These will include constraints for patch induced stress and deformation in the valve regions. FIGURE 13. As presumed, a larger patch radius leads to more power output by both ventricles (see top plots). Therefore, patch peak displacement could be reduced by using larger patches, which again results in lower stress in the tissue. The two plots at the bottom indicate that hardly any effect occurs in the left heart chamber when the duration T is altered. In the right ventricle, however, the peak pressure decreases with shorter durations T. As a reference, healthy p-V loops are shown by gray areas.
Alongside computational studies, also ex-and in vivo experiments are planned. The intermediate goal is a device for emergency situations during open heart surgery.

ACKNOWLEDGMENTS
The authors gratefully acknowledge Simone Deparis and Volkmar Falk for their valuable suggestions and insights. Further, the authors also thank the reviewers; their constructive comments helped to significantly improve the paper.

AUTHOR CONTRIBUTIONS
TK: model development and implementation, code optimization, numerical studies, analysis of results, manuscript writing. SR: development and implementation of original heart model, manuscript revision. SV: validation plan review, result evaluation & interpretation, device design, manuscript revision. SD: validation plan review, result evaluation, defining clinical application limits, manuscript revision. PJ: supervision, modeling advice and problem solving, result evaluation, funding, manuscript drafting and revision.

AVAILABILITY OF DATA AND MATERIAL
Computer simulation results fully available.

CONFLICT OF INTEREST
Author Thomas Kummer, Author Simone Rossi, Author Stijn Vandenberghe, Author Stefanos Demertzis and Author Patrick Jenny declare that they have no conflict of interest.

CODE AVAILABILITY
Custom code implemented within the finite element library LifeV (http://www.lifev.org).

CONSENT FOR PUBLICATION
Not applicable.

CONSENT TO PARTICIPATE
Not applicable.

ETHICS APPROVAL
Not applicable.

FUNDING
Open access funding provided by Swiss Federal Institute of Technology Zurich. Self funding by Prof. Patrick Jenny, Institute of Fluid Dynamics, ETH Zurich.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://crea tivecommons.org/licenses/by/4.0/.