Other conditions that characterize the milieu mammalian articular chondrocyte would also be expected to have very significant electrophysiological consequences, mainly by altering Em. Perhaps the main one of these in chondrocytes is the significant elevation of intracellular [Na+]i levels to (perhaps) as much as 20 mM as opposed to values of 8–10 mM in most other mammalian cells. This elevated [Na+]i may partly explain the somewhat atypical pHi values, since intrinsic Na+/H+ exchange activity and thus pH regulation would be altered. What is likely more important, [Na+]i levels in the 15–30 mM range would strongly regulate (activate) the electrogenic Na+/K+ pump. The resulting increase in outward electrogenic pump current (see **Figure 7A**) would provide a significant hyperpolarizing influence to the chondrocyte Em. This possibility is illustrated

FIGURE 10 | Demonstration of the effects of *in silico* changes in the voltage-dependence for activation of the delayed rectifier K+ current, IK−DR, Panels (A,B), illustrate the altered activation curves and I-V relationships, respectively. Panel (A) shows the derived activation curve for IK−DR (black line) at baseline and after shifting its midpoint ±10 mV to begin an analysis of the effects of changes in surface potential on chondrocyte Em. Panel (B) shows the (black line computed chondrocyte I–V curves for IK−DR at baseline and after the activation of IK−DR is shifted ±10 mV. Panel (C) shows the corresponding changes in resting membrane potential in the baseline model and after [Na+] i . The red line in (C) was computed assuming that intracellular [Na+] i is 15 mM. The resulting augmentation/stimulation of the Na+/K+ pump provides an electrogenic outward current that can significantly hyperpolarize Em (see Discussion section) at all membrane potentials.

by the black traces in **Figure 10C**. Note that when [Na+]i is increased (8.6 to 15 mM) the extra outward Na+/K+ pump current can significantly hyperpolarize the chondrocyte Em. After key principles of chondrocyte Ca2+ transport are further understood, it may be possible to revise and improve related aspects of our chondrocyte model as has been done for the PC 12 cell (Duman et al., 2008). One reason for working toward these changes/improvements is the possibility of gaining new insights into [Ca2+]i regulated apoptosis and autophagy (Harr and Distelhorst, 2010) as well as senescence (Mobasheri et al., 1998).

# Ca2+-Influx and Ca2+-Dependent Currents

As noted, activation of TRPV4 channels would be expected to result in a significant influx of Ca2+ and Na+, and a related depolarization, when the chondrocyte membrane potential is negative to approximately −20 mV, the reversal potential (see **Figures 8**, **9**) for ion flux through these channels. However, any such TRP channel-induced depolarization would be transient, and relatively small. This is because depolarization and/ or Ca2+ influx would be expected to be quickly followed by a K+ efflux. This would result in repolarization back to Em or even beyond (a hyperpolarization). This K+ efflux/outward current would be due to a combination of: (i) the depolarizationinduced activation of the delayed rectifier K+ currents, IK−DR, and/or (ii) turning on of the BK and/or other Ca2+-activated K + current (cf. Horrigan and Aldrich, 2002; Magleby, 2003; Sun et al., 2009). In other nonexcitable cells, even a very small Ca2+ influx can trigger a large Ca2+ release from the ER, thus producing a significant increase in intracellular Ca2+, [Ca2+]i . Moreover, the extent to which [Ca2+]i in the ER changes can trigger a secondary but significant net Ca2+ influx mediated by ICRAC channels expressed in the surface membrane (c.f. Shaw et al., 2013). In a human chondrocyte-derived cell line, Funabashi et al. (2010b) have shown that histamine can strongly enhance a Ca2+-activated K+ current, and thus hyperpolarize the chondrocyte membrane potential. This relatively hyperpolarized transmembrane voltage, Em, is maintained until [Ca2+]i and/or Ca2+-dependent signaling mechanisms reset to "resting" values. This Ca2+ dependent hyperpolarization may also contribute to chondrocyte differentiation (Muramatsu et al., 2007).

# Connexin-Mediated Current Flow and Electrotonic Interactions

As mentioned in the Introduction, chondrocytes in articular joints from adult humans function as isolated, single cells. Interestingly, however, Cell Physiological data from early adolescent articular joint preparations suggest that the growth plate of articular joints includes small groups of closely opposed chondrocytes. In these preparations, expression of selected members of the connexin family of intercellular communication proteins (e.g., Cx43) has been detected using standard immunohistochemical approaches (Chi et al., 2004; Mayan et al., 2013a,b; Asmar et al., 2016).

We also note that is Cx43 increased expression in isolated adult chondrocytes can result in release ATP in response to e.g. mechanical perturbations (Millward-Sadler et al., 2004; Knight et al., 2009; Garcia and Knight, 2010). One plausible mechanism for this, is the occurrence of transient openings of these "hemichannels" that consist of either pannexin or connexin subunits (Saez et al., 2003; Garcia and Knight, 2010; Penuela et al., 2013). Given this important functionality, our mathematical model also includes a connexin-mediated conductance (**Figure 4**). However, since this initial model development was focused on simulating adult chondrocyte behavior under a restricted set of physiological conditions, these connexin single-cell channels were shut off (assigned a conductance value of 0 pS) in all simulations that form the basis of this paper. We acknowledge that this choice does **not** allow this model to account for any of the interesting function properties that arise from the so-called connexon-43 hemichannel behavior in chondrocytes (Knight et al., 2009). This "hemichannel activity" in chondrocytes can also be mediated by the pannexin family of integral membrane proteins (Bond et al., 2011; Matta et al., 2015).

# Limitations of the Mathematical Model of the Human Chondrocyte Resting Membrane Potential

The mathematical model that we have developed is an important advance. However, we recognize that it has significant limitations. These include:

# SUMMARY

This mathematical model of chondrocyte electrophysiology provides a reliable platform for integrating and evaluating both recent and well-established experimental data that is relevant to the generation of the resting potential. At a minimum, given that the chondrocyte is in a unique but relatively inaccessible, environment, our model provides new insights into: the biophysical effects of alterations in ionic strength of synovial fluid on ion channel voltage-dependent gating (zeta potential effects). Key elements of the zeta potential working hypothesis can be tested if recent advances in voltagedependent dye methods are implemented to provide a means of separating the various components of the signal that underlies the "transmembrane potential" in small nonexcitable cells (Cohen and Venkatachalam, 2014). It will also be necessary to account for the effects of cyclic stretch on chondrocyte ion channels by defining the main strain dependent alterations in channel gating voltage dependence or kinetics (cf. Hille, 2001; Maleckar et al., 2009). Approaches for detecting and determining the limitations of present patch clamp methods that can bias key electrophysiological data sets and influence their interpretation must continue to be utilized and refined. However, even in its present form this first-order model will continue to be useful for rationalizing and bringing together genomic data from microarray expression profiles, and understanding ion channel/antiporter drug target initiatives.

# AUTHOR CONTRIBUTIONS

RBC carried out all of the experimental work in this paper, and provided valuable input into the development and validation of the mathematical model. MM when supervising Dr. H. Narayanan at Simula Research Inc in Oslo, Norway lead the mathematical modeling component of this study. BV, a Senior Scientist at Glaxo Smith Kline, Philadelphia, USA provided valuable input into experimental design and facilitated transfer of the GSK properiarty cell lines and pro-drugs. WRG lead the project and wrote the manuscript. All authors have read and revised the manuscript, and are in agreement with its main findings and conclusions.

# ACKNOWLEDGMENTS

An Alberta Innovates—Health Solutions (AI-HS) Scientist Award (WRG), an AI-HS Starter Grant, and a Canadian Institutes of Health Research Grant were combined to support this experimental work and related model development in the Giles laboratory at the University of Calgary. We are grateful to the Southern Alberta Tissue Transplant Facility (Dr. R. Krawetz) for supply of human articular joint (knee) tissue. Ms. Colleen Kondo prepared and maintained the human chondrocyte cultures and was responsible for project management in the Giles laboratory. Postdoctoral Fellowship funding for Dr. Harish Narayanan from both the Simula Research Laboratory, Oslo, Norway and the AI-HS is gratefully acknowledged. Dr. Narayanan provided input into early aspects of the computational work in this paper. Glaxo Smith Kline Laboratories in Philadelphia, Pennsylvania supplied the human chondrocyte cell line and TRPV4 agonists/antagonists that were used to generate some of the data.

### REFERENCES

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.00974/full#supplementary-material

potassium channels in human and equine articular chondrocytes. Osteoarthr. Cartil. 15, 1–8. doi: 10.1016/j.joca.2006.06.017

acutely isolated canine articular chondrocytes. J. Physiol. 557, 93–104. doi: 10.1113/jphysiol.2003.058883

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Maleckar, Clark, Votta and Giles. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Multiphysics Biventricular Cardiac Model: Simulations With a Left-Ventricular Assist Device

Azam Ahmad Bakir 1 \*, Amr Al Abed1 , Michael C. Stevens 1,2, Nigel H. Lovell 1 and Socrates Dokos 1

*1 Graduate School of Biomedical Engineering, University of New South Wales, Kensington, NSW, Australia, 2 Innovative Cardiovascular Engineering and Technology Laboratory, Critical Care Research Group, The Prince Charles Hospital, Brisbane, QLD, Australia*

### Edited by:

*Christopher Basciano, Becton Dickinson, United States*

### Reviewed by:

*Ethan Kung, Clemson University, United States Joakim Sundnes, Simula Research Laboratory, Norway*

> \*Correspondence: *Azam Ahmad Bakir<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

Received: *30 April 2018* Accepted: *21 August 2018* Published: *11 September 2018*

### Citation:

*Ahmad Bakir A, Al Abed A, Stevens MC, Lovell NH and Dokos S (2018) A Multiphysics Biventricular Cardiac Model: Simulations With a Left-Ventricular Assist Device. Front. Physiol. 9:1259. doi: 10.3389/fphys.2018.01259*

Computational models have become essential in predicting medical device efficacy prior to clinical studies. To investigate the performance of a left-ventricular assist device (LVAD), a fully-coupled cardiac fluid-electromechanics finite element model was developed, incorporating electrical activation, passive and active myocardial mechanics, as well as blood hemodynamics solved simultaneously in an idealized biventricular geometry. Electrical activation was initiated using a simplified Purkinje network with one-way coupling to the surrounding myocardium. Phenomenological action potential and excitation-contraction equations were adapted to trigger myocardial contraction. Action potential propagation was formulated within a material frame to emulate gap junction-controlled propagation, such that the activation sequence was independent of myocardial deformation. Passive cardiac mechanics were governed by a transverse isotropic hyperelastic constitutive formulation. Blood velocity and pressure were determined by the incompressible Navier-Stokes formulations with a closed-loop Windkessel circuit governing the circulatory load. To investigate heart-LVAD interaction, we reduced the left ventricular (LV) contraction stress to mimic a failing heart, and inserted a LVAD cannula at the LV apex with continuous flow governing the outflow rate. A proportional controller was implemented to determine the pump motor voltage whilst maintaining pump motor speed. Following LVAD insertion, the model revealed a change in the LV pressure-volume loop shape from rectangular to triangular. At higher pump speeds, aortic ejection ceased and the LV decompressed to smaller end diastolic volumes. After multiple cycles, the LV cavity gradually collapsed along with a drop in pump motor current. The model was therefore able to predict ventricular collapse, indicating its utility for future development of control algorithms and pre-clinical testing of LVADs to avoid LV collapse in recipients.

Keywords: computational modeling, multiphysics modeling, heart, LVAD, cardiac, physiome, electromechanics, fluid-structure interaction

# 1. INTRODUCTION

The optimal function of the heart relies on multiphysics phenomena, comprised of electrical activation, as well as myocardial contraction and blood hemodynamics. A deficiency in any of these physics will impact the others. As such, development of cardiac medical devices has to consider these multiphysics interactions during both pre-clinical and clinical studies. Multiphysics cardiac computational models can be a test-bed for early ideas to assess how therapeutic interventions could impact the diseased heart. This is especially beneficial to reduce the risk to patients with rare medical conditions, where the number of subjects that can be recruited for device clinical trials is limited. Results from computational simulations could be used to generate more plausible hypotheses for future testing; thus, helping to reduce development and technology transfer costs and potential clinical adverse reactions.

Numerous studies have modeled the electromechanical interactions in the myocardium with varying levels of complexity (Usyk et al., 2002; Gurev et al., 2011; Quarteroni et al., 2017). Most were coupled to Windkessel-type models to mimic the systemic circulation (Usyk et al., 2002; Gurev et al., 2011). Instead of modeling the spatially varying fluid dynamics, a uniform endocardial pressure load was calculated to constrain the volume during the isovolumic phases. Whilst such an approach may be sufficient to model most cardiac electromechanics phenomena, it limits model usability: for example, to simulate cardiac interaction with intra-cavity implants. In cases with large intraventricular pressure gradients, such as in hypertrophic obstructive cardiomyopathy (Geske et al., 2011), fluid loading on the endocardium may be spatially heterogeneous. This could render the spatially uniform pressure load constraint inaccurate. Furthermore, formation of vortices, a potential measure of cardiac health index (Hong et al., 2008), cannot be simulated in such models.

Fluid-structure interaction (FSI) models of the ventricles are important in studying fluid loading on the myocardium as well as the heart's blood pumping efficiency. For the sake of simplicity, the contribution of electrical activation on wall mechanics was not considered in many ventricular FSI models. Instead, a temporal activation function was used to govern contraction uniformly across the entire left ventricle (LV) (Nordsletten et al., 2011). Therefore, the effect of spatially heterogeneous activation on fluid behavior cannot be investigated in such models. Several studies have performed simulations of cardiac electrophysiology coupled with FSI (Watanabe et al., 2004; Choi et al., 2015). Nevertheless, their approach only utilized one-way coupling for some of the modeled physics, limiting the scope of physics interactions. A fully coupled fluid-electromechanics model can provide added advantages in terms of predictive power and future-expansion capability to test devices and drugs, or to be applied as a surgery planning tool to predict the heart's motion and blood flow profile post-surgery. The review by Quarteroni et al. (2017) provides a detailed summary of various approaches of modeling cardiac electromechanics as well as fluid-structure interactions.

A potential application of cardiac multiphysics models is to study heart-implant interactions such as the left ventricular assist device (LVAD), classified as a Class III device by the United States' Food and Drug Administration (FDA). As it may pose a high risk to patients, the device needs to be rigorously tested prior to marketing approval. Whilst the LVAD works well for patients with the most severe heart disease, excessively high pump flow rate relative to ventricular filling rate can cause the ventricles to collapse. This phenomenon is known as ventricular suction and it represents a common complication (Vollkron et al., 2007). Moreover, increased arrhythmic tendency has been noted among LVAD recipients, which may reduce right ventricular (RV) output resulting in reduced LVAD preload, worsening the suction (Robertson et al., 2017). Furthermore, several LVAD recipients also suffer from post-implantation RV dysfunction and commonly require an additional RV assist device (Neyer et al., 2016). Positioning and design of the LVAD cannula can also affect its performance (Ong et al., 2013). Therefore, LVAD complications are inherently multiphysics in nature and not only fluid-structure interactions, but also ventricular electromechanics as well as adjacent ventricular function. If any of these issues can be predicted beforehand, a strategy can be developed to lessen the risk.

In this study, we develop a multiphysics modeling framework for simulating biventricular electrophysiological behavior along with myocardial mechanics and two-way bloodmyocardium interaction. The framework is implemented in an idealized biventricular structure that embeds a Purkinje fiber network to trigger myocardial electrical activation. We present computational simulations of a healthy heart as well as LVAD intervention in a dilated heart, including the impact of pump speed on both the pump and ventricles, as well as LV collapse in mitral stenosis case. We utilized commercial computational modeling software for easy replication and future model expansion.

# 2. METHODS

We first present our standard framework for formulating a healthy biventricular model based on our preliminary work presented in Bakir et al. (2017). **Figure 1A** summarizes all physics included along with their coupling. We then describe the implementation of heart-LVAD interaction by introducing a failing dilated heart and adding a LVAD model to the standard framework. Unless otherwise stated, descriptions and parameter values are listed in the **Supplementary Material**.

# 2.1. Standard Biventricular Model Formulations

### 2.1.1. Gross Geometry and Microstructure

The geometry was constructed using two truncated ellipsoids merged at the interventricular sulci, as shown in **Figure 2A**. The aortic and mitral valve boundaries were generated using circular and semilunar shapes, respectively, whilst the pulmonary and tricuspid valve boundaries were created using circular and elliptical shapes, respectively.

model as well as the cannula structure implanted in the dilated heart geometry.

To obtain a realistic activation sequence, an idealized Purkinje fiber tree structure (**Figures 2C,D**) was constructed within the myocardium to mimic the gross anatomy of the Purkinje network (Vigmond and Stuyvers, 2015) without restricting meshing requirements. The Purkinje structure begins at the septal base with two major bundles: the left and right bundle branches. The

left bundle branch was further divided into three minor bundles near the basal septal entry, whilst the right bundle branched into anterior and posterior minor bundles. These minor bundles further branched out circumferentially, every 1 cm between the apex and base, forming the Purkinje-myocyte junctions (PMJ). The Purkinje fibers were represented in the geometry as 1D edges embedded in 3D space: and were numerically assigned a radius, r, that was set to decrease linearly from 500 µm at the basal origin to 50 µm at the endpoints of the PMJs based on the measurements by James (1961).

The myocardial microstructure is governed by the fiber, sheet and normal-to-sheet orientations (LeGrice et al., 1997). For simplicity, the following settings and assumptions were implemented:

3. As the apical microstructure is difficult to measure experimentally (LeGrice et al., 1997), and adaptation of above principles could result in a singularity, the microstructure of the apex region, defined by a 1 cm diameter cylinder, was assumed to be isotropic (Nordsletten et al., 2011).

### 2.1.2. Electrophysiology and Gap Junction-based Propagation

Myocardial action potential (AP) formulations Equations (1)– (5) were based on the phenomenological Nash and Panfilov (2004) model, modified to incorporate units into the originally dimensionless form. Variables Vm and Vp represent myocardial and Purkinje membrane potentials respectively, whilst R is an auxiliary recovery variable for both myocardium and Purkinje fiber. Parameter a in Equations (3) and (5) was chosen to vary linearly such that aepi is 0.12 and aendo is 0.07 to allow the myocardium to relax in the opposite direction of its activation (Glukhov et al., 2010).

$$\beta\_{sv} \left( C\_m \frac{\partial V\_m}{\partial t} + i\_{ion} + i\_{pr} \right) = \nabla\_X \cdot \left( \sigma \nabla\_X V\_m \right) \tag{1}$$

$$r\left(C\_m \frac{\partial V\_p}{\partial t} + i\_{ion}\right) = \nabla\_X \cdot \left(\frac{r^2}{2\rho\_i} \nabla\_X V\_p\right) \tag{2}$$

$$\begin{split} i\_{ion} &= k\_1 k\_2 (V\_m - B) \left( \left[ \frac{V\_m - B}{A} \right] - a \right) (V\_m - 1) \\ &+ k\_2 R (V\_m - B) \end{split} \tag{3}$$

$$i\_{\mathbb{P}^r} = R\_{\mathbb{P}^{mj}} (V\_m - V\_{\mathbb{P}}) \tag{4}$$

$$\begin{split} \frac{\partial R}{\partial t} &= \left( \varepsilon\_0 + \frac{\mu\_1 R}{\left[ \frac{V\_m - B}{A} \right] + \mu\_2} \right) \\ &\times \left( -R - k\_1 \left[ \frac{V\_m - B}{A} \right] \left( \left[ \frac{V\_m - B}{A} \right] - a - 1 \right) \right) \end{split} \tag{5}$$

Membrane potential dynamics and AP propagation within the myocardium were governed by Equation (1), and within the Purkinje fibers by Equation (2), based on the modified cable equation developed by Alqahtani et al. (2017) to take into account the effect of varying fiber radius. Equation (3) represents the transmembrane ionic current flow while Equation (5) governs AP recovery in both the myocardium and Purkinje fiber domains. For the Purkinje fibers, parameters a and k2 were adjusted to achieve a three-fold increase in conduction velocity and a 2-fold increase in AP upstroke velocity relative to their respective values in the myocardium (Vigmond and Stuyvers, 2015). The Purkinje current source, ipr, was supplied to the myocardium only at the the PMJs to trigger myocardial activation (Equation 4).

The ∇X operator in Equations (1) and (2) is defined within the material frame (X, Y, Z) such that ∇X = ( ∂ ∂X , ∂ ∂Y , ∂ ∂Z ) T , corresponding to a gap junction-based constitutive formulation described in an earlier study (Bakir and Dokos, 2015). Except at the apex, the electrical conductivity tensor, σ , was set to exhibit an anisotropy ratio of 4:2:1 along the fiber, Fˆ, sheet, Sˆ, and normal-to-sheet, Nˆ directions respectively, in accordance with Hooks et al. (2007). The conductivity tensor, σ , was determined from conductivities defined in the local fiber (σf ), sheet (σs) and normal-to-sheet (σn) directions according to:

$$
\sigma = \sigma\_f(\hat{\mathcal{F}} \otimes \hat{\mathcal{F}}) + \sigma\_s(\hat{\mathcal{S}} \otimes \hat{\mathcal{S}}) + \sigma\_n(\hat{\mathcal{N}} \otimes \hat{\mathcal{N}}) \tag{6}
$$

The value of σf was set such that the epicardial activation time matched reported ex vivo experimental findings in human heart (Durrer et al., 1970). The electrical conductivity of the isotropic apex was assumed to be σf . Zero flux boundary conditions were implemented for all myocardial boundaries unless otherwise noted.

### 2.1.3. Active Stress

The simplified form of active stress proposed by Nash and Panfilov (2004), and further modified by Göktepe and Kuhl (2010), was chosen to generate mechanical contraction triggered by the AP. Active stress, Ta, was represented by the following phenomenological ordinary differential equation (ODE):

$$\frac{\partial T\_a}{\partial t} = \epsilon (V\_m) \left( k\_{Ta} \left[ \frac{V\_m - B}{A} \right] - T\_a \right) \tag{7}$$

where kTa controls the maximum magnitude of Ta. A, and B are fixed parameters introduced here to give units to the originally dimensionless formulation. The delay function, ǫ(Vm), was given by Göktepe and Kuhl (2010) as follows:

$$\epsilon(V\_m) = \epsilon\_0 + (\epsilon\_\infty - \epsilon\_0) \exp(-\exp(-\xi \left(V\_m - V\_{threshold}\right))) \tag{8}$$

where ǫ0, ǫ∞, ξ and Vthreshold are also fixed parameters. The active stress was set such that the myocardial membrane potential, but not the Purkinje fiber, triggers active stress generation.

### 2.1.4. Myocardial Mechanics

The myocardium was represented by the transverse isotropic myocardial hyperelastic formulation of Holzapfel and Ogden (2009), as detailed in Equations (9)–(11), which can simulate the myocardial biaxial tension response with four parameters. A nearly-incompressible constraint was enforced through the adaptation of a volumetric strain energy function, ψvol in Equation (12), based on the formulation of Doll and Schweizerhof (1999):

$$
\psi = \psi\_{isotropic} + \psi\_{fiber} + \psi\_{vol} \tag{9}
$$

$$
\psi\_{isotropic} = \frac{a\_i}{2b\_i} \exp(b(I\_1 - 3))\tag{10}
$$

$$
\psi\_{fiber} = \frac{a\_f}{2b\_f} \exp(b\_f(I\_{4f} - 1)^2 - 1) \tag{11}
$$

$$
\psi\_{\rm vol} = \frac{\kappa (f - 1) \ln(f)}{2} \tag{12}
$$

I1 in Equation (10) denotes the first invariant of the isochoric elastic right Cauchy-Green tensor, **C**, whilst I4f = Fˆ · (**C**Fˆ). J in Equation (12) represents the determinant of the deformation gradient tensor, **F**. In order to ensure material stability, I4f was set to zero when I4f < 0, assuming that myocardial fibers do not contribute significantly to passive mechanics during compression (Holzapfel and Ogden, 2009). The apical strain energy function was solely given by ψisotropic due to the isotropic apex assumption.

The equation of motion governing myocardial deformation is given by Equation (13):

$$
\rho\_s \frac{\partial^2 \vec{u}\_s^\star}{\partial t^2} + \alpha \rho\_s \frac{\partial \vec{u}\_s^\star}{\partial t} = \nabla\_X \cdot \left( \mathbf{FT} + \beta \frac{\partial \mathbf{FT}}{\partial t} \right) \tag{13}
$$

where #»us is the myocardial displacement. Equation (13) also contains Rayleigh damping, represented by the first order time derivative of displacement and stress, to eliminate unphysiological oscillations during fluid loading. Rayleigh damping parameters were adapted from Fritz et al. (2014) such that α = 100 s −1 and β = 0.01 s.

In order to couple the electrical and mechanical formulations, Ta was added to the 2nd Piola Kirchhoff stress, **T**, along the fiber, sheet, and normal-to-sheet directions, as shown in Equation (14) such that the magnitude of active stress in the sheet and normalto-sheet directions was 40% of its value along the fiber direction (Usyk et al., 2002). In the isotropic apex, an isotropic active stress of Ta was assumed.

$$\mathbf{T} = \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{E}} + T\_a(\hat{\mathbf{F}} \otimes \hat{\mathbf{F}}) + 0.4 \ T\_a(\hat{\mathbf{S}} \otimes \hat{\mathbf{S}}) + 0.4 \ T\_a(\hat{\mathbf{N}} \otimes \hat{\mathbf{N}}) \tag{14}$$

**E** in Equation (14) is the Green-Lagrange strain tensor whilst ψ is the myocardial strain-energy function.

For simplicity, the ventricular base, including the valves, was assumed to be fixed to avoid translation and rotation. The endocardial boundary was also subjected to pressure and viscous stress from the blood fluid, discussed further below. The remaining myocardial boundaries were allowed to move freely.

### 2.1.5. Blood Hemodynamics and Circulation

The fluid velocity, #»vf , and pressure, p, in the blood cavity were determined using the incompressible Newtonian Navier-Stokes formulation (Equation 15) and the continuity equation (Equation 16).

$$\rho\_f \frac{\partial \vec{\boldsymbol{\nu}}\_f^\bullet}{\partial t} + \rho\_f (\vec{\boldsymbol{\nu}}\_f^\bullet \cdot \nabla\_\mathbf{x}) \vec{\boldsymbol{\nu}}\_f^\bullet = \nabla\_\mathbf{x} \cdot \left( -p\mathbf{I} + \mu\_f (\nabla\_\mathbf{x} \vec{\boldsymbol{\nu}}\_f^\bullet + (\nabla\_\mathbf{x} \vec{\boldsymbol{\nu}}\_f^\bullet)^T) \right) \tag{15}$$

$$\nabla\_{\mathbf{x}} \cdot \vec{v\_f} = 0 \tag{16}$$

∇x in Equations (15) and (16) describes the spatial derivative in the spatial frame such that ∇x = ( ∂ ∂x , ∂ ∂y , ∂ ∂z ) T .

No-slip boundary conditions were applied at the ventricular base surface except at the inlets and outlets. FSI was implemented as a strongly coupled two-way framework where (1) the fluid velocity at the endocardium was set to be equal to the endocardial surface velocity (Equation 17), and (2) the fluid total stress, **s**, a sum of pressure and viscous stress, was applied along the normal of the endocardial surface (Equation 18):

$$\left. \vec{v\_f} \right|\_{endo} = \left. \frac{\partial \vec{u\_s}}{\partial t} \right|\_{end} \tag{17}$$

$$\mathbf{s} = \hat{n}|\_{endo} \cdot (-p\mathbf{I} + \mu\_f (\nabla\_{\mathbf{x}} \vec{\nu}\_f^{\star} + (\nabla\_{\mathbf{x}} \vec{\nu}\_f^{\star})^T)) \tag{18}$$

To couple **s** to the solid formulation in Equation (13), a transformation from the spatial frame (i.e. the deformed frame where the Navier-Stokes equations are solved) to the material frame (i.e. the reference or material frame where the solid mechanics are solved) is needed, as shown in Equation (19), where dv and dV are mesh element scale factors in the spatial and material frame respectively.

$$\mathbf{S} = \mathbf{s} \cdot \frac{d\nu}{dV} \tag{19}$$

Here, **S** denotes the traction from the fluid applied to the endocardial wall in material frame coordinates. The scale factors simply scale the local mesh element coordinates to the spatial and material frame coordinates.

The fluid boundary deformation was governed by the myocardial structural deformation, #»us , whilst deformation of the mesh within the fluid domain was governed by a mesh smoothing formulation (Equation 20). The hyperelastic mesh smoothing algorithm in Equation (20) was solved to determine the minimum mesh deformation energy, ψ**mesh**, within the blood domain, f :

$$\psi\_{m\text{sch}} = \frac{1}{2} \int\_{\Omega\_f} \mathbf{C}\_1 (I\_{1,m\text{sch}} - \mathbf{3}) + \mathbf{C}\_2 (I\_{1,m\text{sch}} - \mathbf{3})^2$$

$$+ \kappa\_{m\text{sch}} (I\_{m\text{sch}} - \mathbf{1})^2 \, dV\_{\Omega\_f} \tag{20}$$

The invariants, Jmesh and I1,mesh, were given by Equations (21) and (22):

$$J\_{mesh} = \det(\nabla\_{\lambda} \mathbf{x}) \tag{21}$$

$$I\_{1,msh} = (I\_{msh})^{-2/3} \operatorname{tr}((\nabla\_X \mathbf{x})^T \nabla\_X \mathbf{x}) \tag{22}$$

where C1 and κmesh are artificial shear and bulk moduli, both set to a value of 1, whilst ∇Xx = ∂x ∂X , where **x** represents the spatial frame coordinates and **X** the material frame coordinates. The non-linear stiffening parameter C2 was set to 0 for the healthy simulation and 100 for the LVAD simulations to cope with greater mesh deformation subsequent to cannula inclusion. Boundary conditions for Equation (20) were set to equal the displacement of the endocardial wall, #»us , and held fixed at the ventricular base.

The finite element model was coupled to a closed-loop Windkessel circulation (**Figures 1C,D**), based on the Kerckhoffs et al. (2007) circuit, modified to include an unstressed volume parameter in the capacitances, and to reflect healthy human circulation. Laminar flow rate conditions were implemented at the inlet and outlet boundaries. Each flow rate was determined from the difference between the average boundary pressure, taken at each of the inlet and outlet boundary of the finite element model, and the Windkessel circulatory pressure. The flow rate variables, described in the **Supplementary Material**, were applied as follows: Qmi was applied to the mitral inlet, Qao to the aortic outlet, Qtri to the tricuspid inlet, and Qpa to the pulmonary artery outlet.

These laminar inflow and outflow boundary conditions were implemented by appending a fictitious tube to each boundary, and applying:

$$[L\_{\rm ext} \nabla\_{\mathbf{x},t} \cdot [-p\mathbf{I} + \mu(\nabla\_{\mathbf{x},t} \vec{v}\_f^\star + (\nabla\_{\mathbf{x},t} \vec{v}\_f^\star)^T)] = -p\_{\rm ext} \vec{\eta}^\star \tag{23}$$

where ∇x,t is the tangential derivative in spatial frame at the inlet or outlet surfaces. This equation essentially determines the velocity applied at the boundary, assuming the inlet and outlet boundaries are connected to a tube of Lext external length. A constraint, pext, was applied to ensure that the flow rate through the tube matched the assigned flow rate from the above. To ensure a laminar profile and prevent back-flow when the valve is closed, Lext needs to be sufficiently large to eliminate numerical backflow was set to 20 m in our model. It should be noted that the Lext does not represent an actual physical tube length, but a value chosen to satisfy the conditions.

### 2.1.6. Simulation Protocol

1. To begin, the Windkessel circuit was decoupled from the heart and the LV and RV chambers were inflated until they reached their end diastolic volumes, Vlv,∞ and Vrv,∞, respectively. The filling rate for each ventricle was calculated using Equation (24) based on a proportional controller:

$$Q\_{\rm fill} = -k\_{\rm flow}(V\_{l\nu} - V\_{l\nu,\infty})\tag{24}$$

where kflow is a fixed parameter set to 200,000 s−1 .

2. Using the filled ventricle, the Windkessel circuit was connected and the model was solved using the initial Windkessel volumes listed in the **Supplementary Material**, which were tuned to be near their steady-state cycle' volumes. The model was run for three steady cycles. A pulse stimulus of 2 ms duration at 60 beats per minute (bpm) was applied at the basal entry of both left and right bundles of the Purkinje network.

### 2.2. LVAD Simulations

A left ventricular assist device (LVAD) model was attached to the biventricular model via a cannula.

### 2.2.1. Heart Failure Modifications

Prior to adding the LVAD, the biventricular model was modified to reflect dilated cardiomyopathy (DCM), a common etiology of LVAD recipients. The geometry was modified by increasing the LV short axis and reducing its wall thickness as shown in **Figures 2E,F** (Mohiaddin, 1995; Dandel et al., 2008). This increased the LV chamber size whilst the RV chamber size was minimally affected. The heart rate was increased to 80 bpm (Mohiaddin, 1995; Dandel et al., 2008) and the LV contraction strength was reduced by setting the kTa parameter of the standard model to half its value, in accordance with isometric twitch measurements in isolated epicardial tissue obtained from DCM patients (Hasenfuss et al., 1992). Other settings described in the standard model were left unchanged.

### 2.2.2. LVAD Model

The LVAD was represented by an ODE model that relates the hydraulic state experienced by the pump and pump motor current. The LVAD model, developed by Lim et al. (2010), consists of three components:

1. Motor windings electrical equation, describing the current, I, needed to maintain the pump motor speed, ω:

$$V\_{pump} = -2k\_{\text{e}}\,\omega + R\_{pump}\,I + L\frac{dI}{dt} \tag{25}$$

The voltage, Vpump, was controlled using a proportional controller, which maintains the pump speed at the desired speed, ωset:

$$V\_{pump} = k\_{pump} \left(\omega - \alpha\_{set}\right) \tag{26}$$

2. Electromagnetic torque transfer equation, relating the electromagnetic torque, Te, produced by the pump to the input flow rate, Q (L min−1 ):

$$T\_e = 3k\_e I = J\frac{d\boldsymbol{\alpha}}{dt} + a\_p \, Q^2 \, \boldsymbol{\alpha} + b\_p \, Q \, \boldsymbol{\alpha}^2 + c\_p \, \boldsymbol{\alpha} + d\_p \, \boldsymbol{\alpha}^3 \tag{27}$$

3. Pump hydraulic equation, relating the differential pressure, 1P, produced by the pump to the input flow rate and pump speed:

$$
\Delta P = e\_{\mathcal{P}} + f\_{\mathcal{P}} \, \mathcal{Q}^3 + \mathcal{g}\_{\mathcal{P}} \, \mathcal{o}^2 \tag{28}
$$

Pump inflow rate was obtained via the equation:

$$\frac{dQ}{dt} = \frac{\Delta P - (P\_{as} - P\_{lv,annual}) - (R\_{in} + R\_{out}) \text{ Q}}{L\_{in} + L\_{out}} \tag{29}$$

Variable inflow and outflow resistances, Rin and Rout, were set to be proportional to the pump flow rate to account for the turbulence effect in the pump (Lim et al., 2010), and are given by

$$R\_{in} + R\_{out} = k\_r \ Q \tag{30}$$

The LVAD equations were linked to the ventricular model via the numerical difference between LV pressure averaged at the cannula outlet boundary, Plv,cannula, and the systemic arterial pressure, Pas of the Windkessel circulation as shown in Equation (29). This difference is the pressure head of the pump, 1P, minus the pressure loss across the cannula. Equations (25)–(29) describe the relationships that determines the overall pump flow rate, Q, using this pressure difference.

### 2.2.3. Cannula Design

A cylindrical generic pump cannula, 1.6 cm inner diameter and 0.2 cm wall thickness, was inserted at the LV apex of the biventricular geometry as shown in **Figure 1C**. The cannula extends only a quarter of the LV apicobasal distance as suggested by Ong et al. (2013). The cannula was set to be a linear elastic silicon material (E = 170 GPa, ν = 0.28, ρs = 2,329 kg m−3 ), was held fixed in space at its external parts outside the LV.

Contact modeling between the LV endocardium and the cannula is necessary to prevent mesh collapse when surface contact occurs. In essence, the contact formulation applies loads that are inversely proportional to the distance between contacting surfaces when the surfaces move closer than a certain threshold. A wall distance formulation (Fares and Schröder, 2002), described in the **Supplementary Material**, was applied to calculate the distance between the endocardial surface and contact destination surfaces (cannula and ventricular base). The ventricular base of the fluid domain was included as a contact destination, due to the possibility of the collapsing wall contacting the basal surface, inducing mesh collapse.

A boundary load Pcontact, given by Equation (31), was applied along the normal direction of the endocardial surface to counter its movement toward the contact destination surfaces:

$$P\_{\text{contact}} = \begin{cases} k\_{\text{contact}} (\text{G}\_{\text{contact}} - \text{G}\_{\text{contact,lim}}) & \text{G}\_{\text{contact,lim}} \le \text{G}\_{\text{contact}}\\ 0 & \text{otherwise} \end{cases} \tag{31}$$

Gcontact is the reciprocal of distance between surfaces, whilst Gcontact,lim is the set threshold (150 m−1 ). kcontact was set to 103 N mm−1 , which was sufficient to prevent the wall contact. To enhance the blood domain mesh deformation with the added cannula, parameter C2 in Equation (20) was set to 100, introducing a non-linear stiffening to the mesh smoothing formulation.

### 2.2.4. Simulation Protocol

### 2.3. Mesh Settings

The mesh was constructed using a combination of tetrahedral, pyramidal and prism elements. A fluid boundary layer of two elements was placed at the fluid domain side of the fluid-structure interface. The outlet and inlet surface meshes were refined by setting their edges to an average element size of 0.1 cm. The resulting normal mesh is shown in **Figure 1C**. Mesh element size was determined as the length of the element's largest edge. The average mesh size for the myocardium was 0.52 cm whilst in the fluid domain it was 0.4 cm, which resulted in 22661 myocardial mesh elements and 44078 blood domain elements. A mesh convergence analysis was conducted on the standard model by reducing the average element size to 0.35 and 0.31 cm for the myocardium and fluid domains respectively, defined as the finer mesh setting from here on. The LVAD simulation was computed using a similar setting to the normal mesh to the standard healthy heart simulation, but with an average element size of 0.21 cm within the cannula. This resulted in 28005 myocardial mesh elements and 103569 blood domain elements.

Mesh convergence analysis was performed for the biventricular model using the following quantities:

1. Characteristic quantities, Cq

$$C\_q = \frac{1}{N} \sum\_{n=1}^{N} \frac{1}{M} (\sum\_{m=1}^{M} \zeta\_{mn}) \tag{32}$$

2. Relative root mean square, %RMS

$$\% \overline{RMS} = 100 (\frac{\sqrt{\frac{1}{M} \sum\_{m=1}^{M} (\zeta\_{m,normalmech} - \zeta\_{m,refinedmech})^2}}{\sqrt{\frac{1}{M} \sum\_{m=1}^{M} (\zeta\_{m,refinedmech})^2}}\tag{33}$$

where the ζmn quantities represent the local electrical activation time, myocardial displacement and fluid pressure at the times of peak LV ejection and peak LV filling rate at a set of 40x40x40 material datapoints generated across the biventricular structure. The magnitudes of the first derivatives of Vm, uEs , vEf and p were extracted, and electrical activation time was determined by finding the time of maximum first derivative of Vm at each point.

### 2.4. Solver Settings

The physics involved in the model are summarized in **Figure 1A**. In brief, the myocardial AP (Equation 1) was triggered by the Purkinje current (Equation 2). AP activation triggers the generation of active stress in Equation (7). This active stress was in turn coupled to the transverse isotropic hyperelastic myocardial formulation by directly adding the active stress tensor (Equation 14) into the 2nd Piola-Kirchoff stress in Equation (14). The latter was determined from the derivative of the strain energy function with respect to the strain components, as indicated in Equation (14). The electromechanics formulations (Equations 1– 14) are implemented only in the myocardial domain and solved in the material frame, to consider for gap junction-controlled AP propagation as described in Bakir and Dokos (2015).

The fluid physics (Equations 15–16) was solved in the spatial frame within a moving mesh. The fluid-solid coupling occurs at the endocardial surfaces (Equations 17–19), with the fluid velocity set to the endocardial wall velocity and fluid total stress, **S**, applied as a boundary load for the myocardial mechanics. Deformation of the mesh in the blood domain was governed by the hyperelastic formulation (Equation 20), with boundary conditions given by the endocardial displacement (Equation 13) while the basal plane was held fixed.

We utilized COMSOL Multiphysics 5.2a (COMSOL A.B., Sweden) for generating the geometry and mesh as well as solving all formulations. Fully-coupled PARDISO direct linear solver with an automatic damped Newton method were implemented, where all formulations were assembled in a single matrix and solved simultaneously together with preordering algorithm set to nested dissection multithreaded. Streamline diffusion and crosswind diffusion artificial stabilizations of the Navier-Stokes equations were applied, which enabled same-order velocitypressure numerical integration and eased the requirement for a very fine fluid domain mesh. We adopted P1-P1 linear Lagrange elements for the Navier-Stokes equations and quadratic Lagrange elements for the membrane potential potential and recovery, R, variable as well as for all the solid mechanics formulations. All formulations were solved simultaneously and all variables were updated at each time step using a second order Backward Differentiation Formula (BDF) adaptive time-stepping method with an event detection algorithm implemented to ensure the BDF solver did not miss the stimulus events. Simulation results were output at 2 ms time steps.

# 3. RESULTS

### 3.1. Standard Biventricle Model 3.1.1. Overall Multiphysics Behavior

Electric activation, myocardial motion, and blood velocity vector streamlines are shown in the snapshots of model behavior in **Figure 3** and in the **Supplementary Animation 1**. The snapshots were taken at various phases in the third cycle of the steady cycle model.

The model started from the end-diastolic volume (EDV) phase, where vortices from the previous filling phase still persist in both ventricles. The stimulus began at the basal entries of the Purkinje fibers. As the myocardium began to electrically activate, a torsional motion was noted. The blood flow vortices began to diminish and the flow was gradually directed toward the outlets until ejection was initiated. As the heart reached total electrical activation, blood flow streamlines were oriented toward the outlets in near straight lines. The apex also moved toward the base during contraction. The straight streamlines gradually diminished as the heart entered isovolumic relaxation.

AP relaxation first occurred in the epicardium, whilst the endocardium repolarized last: the opposite of the activation direction. The filling was initiated even before the heart contraction fully relaxed. As the ventricles were being filled, blood vortices formed near the inlets. The vortices grew larger and moved apically. LV vortex shape was somewhat symmetrical whilst the RV displayed a more elongated and crescent-like shape due to the RV geometry. At the end of the cycle, the ventricles returned back to their starting EDV shape with larger vortex presents within the central region of the ventricular cavities. Smaller vortices can also be noted in the lateral regions nearing the myocardial wall, near the inlets.

### 3.1.2. Electrical Activation

Following external stimulation of the basal entries of the Purkinje network, the AP propagated through the branching fibers toward the apex. The Purkinje fiber network took 34 ms to be wholly activated. As the AP traveled through the PMJ edges the myocardium was excited. The site of first myocardial activation occurred at the septal region 20 ms after the Purkinje base was stimulated.

LV activation was initiated at three endocardial sites that merged into a single wavefront, which propagated radially toward the epicardium before finally reaching the base (**Figures 4A,B**). Epicardial breakthrough occurred first in the anterior and posterior interventricular sulci. These breakthroughs formed two wavefronts that propagated toward the unactivated regions in the RV. These wavefronts moved in a near tangential fashion at the free wall of the RV toward the base. Overall, it took 94 ms for the whole structure to be activated, with the RV free wall base the last activated region. On the other hand, the repolarization sequence was the opposite of the activation sequence such that it was initiated at the epicardium rather than the endocardium.

The action potential and active stress profiles are displayed in **Figures 4C–H**. The APD90 measured at the endocardium and epicardium was 324 and 266 ms, respectively. The presence of this APD gradient was sufficient to generate a recovery sequence opposite to that of the activation sequence. The Purkinje AP morphology exhibited an upstroke velocity twice that of the myocardium and the APD was longer (**Figure 4D**). At the Purkinje-myocyte junction (PMJ), the myocardial AP appeared similar to the Purkinje AP morphology. The time to peak active stress was 192 and 222 ms whilst the time to half relaxation was 142 and 156 ms in the epicardium and endocardium, respectively. These measures were extracted from three material points, located mid-apicobasal distance of the LV free wall, septum and RV free wall.

### 3.1.3. Mechanical Variables

The fiber, sheet, and normal-to-sheet components of the 2nd Piola-Kirchoff stress and Green-Lagrange strain tensors were extracted from three material points, located at the epicardium of the LV free wall, septum and RV free wall. These points were taken at midway between the base and apex.

To assess regional work, stress-strain loops were generated (**Figure 5**) as in Russell et al. (2012). The fiber direction was the dominant direction of deformation and stress. The LV free wall stress-strain loop was the largest compared to the those measured at the septum or RV free wall, indicating greater regional work in the LV free wall relative to the other regions. The stress strain loops for the sheet and normal-to-sheet components showed more skewed and distorted loops.

Torsional motion was assessed at three points in the myocardium: LV free wall, septum and RV free wall. These material points were positioned near the apical region at three quarters of the apicobasal distance. All points were placed at the epicardial surfaces. Torsion was mostly predominant in the LV free wall, with 16.5◦ counterclockwise motion when viewed upward from the apex, as shown in **Figure 6**. The torsion was substantially less in the RV free wall and septal region with magnitude of 5◦ . During the cardiac cycle, an apicobasal displacement of 1.2 cm was also noted in the model.

Supplementary Animations 1, 2 respectively.

activation. (B) Electrical activation time map in the standard model. Comparison of (C) LV free wall myocardium endocardial and epicardial, (D) Purkinje fiber, and (E) Purkinje-myocyte junction action potential profiles. (F) Active stress profile of the myocardium taken at the same location as (C). The time derivative of membrane potential of the (G) Purkinje fiber and (H) Purkinje-myocyte junction, taken at the same locations as in (D,E).

### 3.1.4. Hemodynamics

**Table 1** lists the hemodynamic quantities measured at the third cycle of the biventricular model, compared against the normal range of healthy humans.

The pressure-volume (PV) loop plots in **Figure 3** depict the commonly observed rectangular shape for these loops. This indicates that the four cardiac cycle phases were simulated: isovolumic contraction, ejection, isovolumic relaxation, and filling. As the model reached stable cycles, the stroke volume of both ventricles equalized, as shown in the PV loop and **Table 1**. At maximum ejection rate, aortic velocity was 3.18 m s−1 , whilst the pulmonary velocity was 2.94 m s−1 , as shown in **Figure 3**. The mitral velocity was 1.04 m s−1 during peak filling rate, whereas the tricuspid inlet velocity was 1.07 m s−1 .

During peak ejection rate, a pressure gradient existed predominantly at the outlets (**Figure 7**). The pressure difference

FIGURE 5 | Stress vs. strain plots taken from the (A) fiber, (B) sheet, and (C) normal-to-sheet components of the 2nd Piola-Kirchoff stress and Green-Lagrange strain tensors. The results are taken from the standard healthy model at the third cycle (*t* = 2,000–3,000 ms).

between the outlet and the inner cavity region were 31.3 mmHg and 28.9 mmHg in the LV and RV, respectively. During the filling phase, LV and RV pressure distributions showed more dispersed patterns, with the lowest pressure at the vortex centers. The LV, RV, aortic and pulmonary arterial pressure waveforms rose and dropped at nearly the same time (**Figure 6**). The time to peak pressure for the LV and RV also appeared to occur at a similar time, indicating a synchronous behavior. Similarly, peak filling rates occurred in near synchrony, where the time of ejection start and end occurred at the roughly the same time (outflow rates plot **Figure 6**).

The mean kinetic energy (KE) of blood in the ventricles was computed using a similar method to Carlsson et al. (2012):

$$\overline{\rm KE} = \rho\_f \int |\vec{\nu\_f}|^2 \, dV\_{\Omega\_f} \tag{34}$$

where ρf is the blood density, vf blood velocity, and dVf the ventricular fluid cavity. KE of both ventricles is shown in **Figure 6C**. KE was near zero prior to the start of the first contraction as the blood velocity magnitude was near zero, and was largest during the ejection phase of the cardiac cycle with the LV showing the larger magnitude of KE. During isovolumic relaxation, KE dropped and rose again during the subsequent filling phase. As filling phase ceased, KE gradually declined but never dropped to near zero. Consequently, the maximum KE of the subsequent cycle's ejection phase was higher than the first cycle's.

The average vorticity, |ω|, was also measured by averaging the vorticity magnitude across the entire ventricular cavities. Vorticity, ωE was obtained via the the curl operator on the blood velocity vectors:

$$
\vec{w} = \nabla \times \vec{\nu\_f} \tag{35}
$$

The |ω| waveform showed two main peaks in both LV and RV (**Figure 6D**): (1) during the ejection phase and (2) during the filling phase. The filling phase exhibited a larger magnitude of |ω|, which is also qualitatively visible in the snapshots of **Figure 3** with the blood swirling motion.

Following ventricular contraction, |ω| increased momentarily before dropping significantly, but then increased again during the filling phase. Overall, the RV and LV |ω| showed comparable magnitudes and waveform phases throughout the cardiac cycle, indicating synchronized behavior.

### 3.2. Mesh Convergence

Mesh sensitivity analysis showed good convergence (< 5%) for global characteristic quantities, Cq, of the mechanical and fluid dynamics, as summarized in **Table 2**. A more local measure, %RMS comparing normal and finer meshes, revealed a good convergence (< 5%) for local electrical activation time, myocardial displacement, and fluid pressure at the time of peak LV ejection and peak LV filling rate. Whilst the fluid velocity showed a good %RMS convergence at the time of peak LV ejection between the two meshes, the same measure was slightly higher at the time of LV filling rate.

# 3.3. LVAD Simulations

The model simulated the general impact of LVAD on the failing heart as well as, the impact of pump speed, reproducing ventricular collapse, when mitral stenosis was introduced.

### 3.3.1. Heart Failure Without LVAD

By reducing LV contractility and inducing heart dilation, the LV ejection fraction (EF) dropped to 16.5% and the LV end-diastolic volume (EDV) increased to 299 ml. EF was also reduced in the healthy RV to 40%. The LV pressure, along with aortic pressure, dropped to a peak systolic pressure of merely 89 mmHg, whilst the RV peak systolic pressure reached 32 mmHg.

The qualitative behavior is shown in the snapshots of **Figure 3** and the **Supplementary Animation**. Overall, minimal change was noted in the electrical activation sequence. The most obvious LV changes are a weaker contraction, blood vortices persisting even during systole, and a slower peak aortic ejection. On the other hand, the RV showed minimal changes relative to the standard model.

### 3.3.2. Low LVAD Pump Speed

Snapshots of the fluid streamlines in **Figure 8** during systole show that the aortic ejection persisted with the LVAD, albeit at lower velocity magnitude. Due to the presence of two outlets during systole, a low speed vortex region was observed in the LV free wall region. The smallest LV volume occurred at about the same time as the initiation of ventricular filling. As the LV was being filled, a fluid vortex formed at the LV free wall region and traveled toward the apex. Toward the end of the filling phase, the vortex separated into two main vortices: one near the basal LV free wall and one near the cannula. The RV showed typical systolic and diastolic streamlines as was seen in the healthy model. Throughout the cycles, the LV endocardium was well out of contact with the cannula.

Throughout the three cycles, the PV loops overlapped, indicating stability (**Figure 9A**). It can be seen that aortic ejection was preserved at this pump speed, although lower than the pre-LVAD setting. The aortic pressure peaked at 95 mmHg and reached a minimum of 86 mmHg, displaying less pulsatility (**Figure 10**). The LV pressure waveform observed at the LV base and the pressure waveform at the cannula showed synchronous characteristics, and were nearly identical. This was also illustrated in the LV pressure distribution of **Figure 11** where the LV pressure is lowest at the cannula. The LV pressure gradients were highest during the systolic phase (≈ 9 mmHg) and lowest at diastole (≈ 2 mmHg).

The LVAD model also provides measures of pump electrical current, pump flow rate and speed. With the proportional controller applied, the pump speed was nearly constant with very insignificant (< 1%) drops during systolic contraction. The pump motor current showed a bump to 0.57 A from a baseline of 0.35 A, coinciding with increased pump flow rate to a peak of 6.2 L min−1 from a baseline of 2.2 L min−1 . Qualitatively, it can be seen that the pump motor current and flow rate waveforms coincided with the LV pressure waveform due to direct correlation with 1P in Equations (25)–(28).

### 3.3.3. Sudden Increment of Pump Speed

A sudden increase of pump speed to 3,150 RPM immediately resulted in an increased aortic pressure, causing a cessation of aortic ejection in the subsequent cycles (**Figure 10**). Both peak LV and RV pressures also gradually dropped from cycle to cycle as the high pump speed was maintained. Moreover, LV pressure gradient between the cannula and ventricular base increased, as illustrated in the pressure snapshots of **Figure 11**. At the end of the third cycle, the pressure became negative in the cannula region (**Figure 11**).

As shown in **Figure 10**, the increase in speed was accompanied by an immediate increase in the pump flow rate from the previous baseline of 2.2 L min−1 to a new baseline of 7 L min−1 . During systole, the pump flow rate peaked at ≈9 L min−1 as opposed to 6.2 L min−1 at the lower, 2,100 RPM, speed. Consequently, the

FIGURE 8 | Snapshots of the biventricular model with LVAD, depicting the membrane potential, myocardial movement and fluid velocity streamlines at the times of peak LV pressure, smallest LV volume, and end-diastole. The snapshots are taken at the third cycle of stable 2,100 RPM pump speed (Top), fourth cycle following speed increment to 3,150 RPM (Middle), and third cycle of mitral stenosis (Bottom). Streamline density is proportional to the blood velocity magnitude. The animations of these snapshots are available in Supplementary Animation 3 for the 2,100 RPM setting, Supplementary Animation 4 for speed increment setting and Supplementary Animation 5 for mitral stenosis case.

FIGURE 10 | (A) The LV, aortic, RV and pulmonary arterial pressure waveforms from the third cycle of stable 2,100 RPM pump speed, speed increment to 3,150 RPM, and mitral stenosis. (B,C) The aortic, mitral, pulmonary and tricuspid flow rate waveforms (B) along with ventricular volume waveforms (C) with similar LVAD settings as the top row. The pump current (D), pump flow rate (E), and pump speed (F) waveforms from the third cycle of stable 2,100 RPM pump speed, speed increment to 3,150 RPM, and mitral stenosis.

pump current also elevated to a baseline of 1.2 A with a maximum current of 1.5 A.

With the cessation of aortic ejection, the blood flow streamlines during systole were directed toward the cannula at a higher velocity, as shown in the snapshots in **Figure 8**. Slight swirling was also observed in the region surrounding the cannula. Similar to the low speed setting, the smallest LV volume occurred at about the same time as the mitral filling phase. At end-diastole, both ventricles exhibited fluid vortices. The cannula and the LV endocardium were not in contact throughout the simulation (**Figure 8**).

As the cycles went on, the LV pressure-volume loops gradually became more triangular (**Figure 9B**). The LV EDV also gradually reduced, along with overall LV blood pressure characteristics. The increased speed also appeared to improve RV function, with enhanced ejection fraction and lower peak RV pressure. The decongested LV also resulted in a small, yet gradual, wall collapse as quantified by the decrease in distance between the LV free wall and septum (**Figure 12**).

The stress-strain loops measured at the epicardial sites along the mid-apicobasal distance, LV and RV free walls and septum, suggest a potential change in regional work in the ventricles. The fiber components of the stress-strain loop showed a gradual shift to the left as the speed increased, especially for the LV free wall (**Figure 13**). On the contrary, minimal change was noted in the RV stress-strain loops.

### 3.3.4. Mitral Stenosis With LVAD

In the first cycle with mitral stenosis, the systolic measures did not appear to be affected by the stenosis. During diastole, the introduction of mitral stenosis severely impeded the inflow into the LV, resulting in an immediate imbalance between the pump flow rate and input flow rate into the pump. As such, the LV collapsed as the cycles continued.

The major impact of mitral stenosis appeared from the second cycle onwards, in which the ejected blood was not replenished in the LV for the subsequent cycles, resulting in a lower LV EDV in subsequent cycles. Consequently, the LV pressure generated from

FIGURE 11 | Snapshots of the biventricular model with LVAD depicting the LV blood pressure distribution at the times of peak LV pressure, smallest LV volume, and end-diastole. The snapshots are taken at the third cycle of stable 2,100 RPM pump speed (Top), fourth cycle following speed increment to 3,150 RPM (Middle), and third cycle of mitral stenosis (Bottom).

the second cycle onwards was significantly lower (52.9 mmHg) than prior to mitral stenosis (86.6 mmHg). The LV pressure dropped even further in the third cycle to 30 mmHg, nearly similar to the RV pressure, as shown in **Figure 9C**.

The collapsing LV also resulted in the LV PV loop to be smaller than the RV (**Figure 9**). The LV PV loop also shifted downwards and leftward toward the negative pressure region. The RV PV loop was also affected, as it shifted to the right, since the RV started to dilate with a slight increase in peak RV pressure. In addition, the fiber direction stress-strain loop, exhibited in **Figure 13**, also shrank significantly for the LV free wall and septum across all strain and stress tensor components.

The snapshots in **Figure 8** showed an obvious shift in the septum orientation toward the LV. Mitral flow was noted to initiate earlier, and last for a longer duration (**Figure 10**). The smallest LV volume occurred at an earlier phase of the cardiac cycle, compared to previous settings. During end-diastole, the fluid vortices were notably absent in the LV, and only present in the RV. The snapshots of LV blood pressure distribution in **Figure 11** also revealed greater negative pressure, with the filling phase showing an entirely negative pressure throughout the whole chamber. With the ventricular wall collapsing into the LV, it can be seen that the endocardium was in contact with the cannula throughout the cardiac cycles. The distance between LV free wall and septum was reduced significantly in the presence of mitral stenosis (**Figure 12**).

The reduction in LV pressure also manifested in a reduced peak pump flow rate and pump current (**Figure 10**).

FIGURE 12 | Distance between the LV free wall endocardium and LV septal endocardium midway between the apex and the base. The lower 2D plots depict the cross-section of the ventricles taken at midlevel along the apicobasal direction for all LVAD simulation settings. The red surface indicates the myocardial cross-section at time of smallest LV whilst the black outline indicates the myocardial cross-section at time *t* = 0. "X"; denotes the two points used to calculate the wall to wall distance.

Nevertheless, the baseline pump flow rate and current were maintained at approximately 6 L min−1 and 1.3 A, respectively. The waveforms also became less pulsatile with the onset of ventricular collapse.

# 4. DISCUSSION

In this study, we developed a cardiac multiphysics modeling framework for biventricular electrical, mechanical and fluid dynamics physics in a healthy heart. The modeling framework was then applied to simulate a failing heart with LVAD support. This framework can be applied in future pre-clinical investigations of device performance under various pathological conditions.

### 4.1. Standard Biventricular Model 4.1.1. Electrical Characteristics

The simulated LV activation pattern is nearly radial and the RV activation pattern is nearly tangential ending at the basal region, whilst the interventricular septum is predominantly activated in a left to right sequence. The first epicardial breakthrough occurred in the RV due to its thinner wall structure. As such, the model's electrical activation reproduced normal human electrical activation sequence (Durrer et al., 1970). Furthermore, the total duration of electrical activation, 94 ms, falls within the range of the QRS duration in healthy humans (80-120 ms), which corresponds to the ventricular activation time (Guyton and Hall, 2006). The model also reproduced the transmural APD gradient observed in the healthy myocardium. In addition, AP repolarization began from the epicardium and ended at the endocardium, the opposite direction to the activation. Such behavior has been recorded in ventricular tissue strips from rejected donor hearts by Glukhov et al. (2010).

The Purkinje fiber network implementation enabled simulation of a healthy activation sequence, and can be applied in future studies to simulate diseased activation such as bundle branch block by simply disabling select branches of the Purkinje tree. The Purkinje fiber APD is longer than that of the nearby contractile myocardium to prevent retrograde activation (Vigmond and Stuyvers, 2015). Although current flow from the myocardium to the Purkinje fibers was not included in this study, retrograde activation can be considered in the future. For example to simulate ectopic beat re-entry, an additional current could be added to the Purkinje AP model of Equation (2) (Vigmond and Stuyvers, 2015).

AP propagation physics was formulated in the material frame to emulate gap junction controlled propagation as described in our previous work (Bakir and Dokos, 2015). Since the resistivity of the gap junctions is higher than the intracellular space of cardiomyocytes, the myocardial substrate can be treated as network of interconnected resistors whose total electrical resistance remains unchanged as the structure deforms. This phenomenon has been experimentally observed in a number of animal species by Penefsky and Hoffman (1963), who measured papillary muscle conduction velocity while they were stretched. For this reason, formulating the electrical propagation in the material frame represents a more suitable approach. Although, the use of material frame-based formulation may have an insignificant effect on rhythmic activation patterns, we have shown that in more complex electrical activation scenarios, such as arrhythmias, the effect is prominent (Bakir and Dokos, 2015).

### 4.1.2. Mechanical Variables

The predicted apex to base shortening of 1.2 cm is within the range measured in healthy humans (Alam et al., 1990). The torsion was predominantly within the LV due to its larger muscle mass, as opposed to the RV. The LV wall torsion magnitude was within the range reported by Henson et al. (2000). On the other hand, RV torsion measures have been under-reported in the literature. Nonetheless, it is expected that the RV exhibits lesser twist as its thinner wall accommodates lesser circumferential fibers and more oblique fibers (Pettersen et al., 2007). The RV also restricted the septal motion generating lesser torsion in the septum as opposed to the LV free wall. Thus, the common healthy ventricular mechanics were observed in our biventricular model.

Stress-strain loops are indicators of regional ventricular work, which may be correlated to myocardial oxygen consumption and potential structural remodeling (Russell et al., 2012). Our model predicted a relatively large and box-shaped loop for the LV free wall, indicating a greater work load on this section of the myocardium to eject blood, compared to those at the LV septum. This is in agreement with a clinical study that reported regional work is larger in the LV free wall compared to LV septum (Russell et al., 2012). The predicted RV free wall stressstrain loop was smaller than the LV free wall. This is due to the lower afterload pressure imposed on the RV, even though both ventricles produced a similar cardiac output.

### 4.1.3. Hemodynamics

Global hemodynamics measures predicted by our model were similar to normal quantities reported in healthy humans. The simulated LV ejection fraction and stroke volume were on the lower end of normal ranges, possibly due to the smaller amount of myocardial thickening exhibited by our model. We speculate that this in turn was due to our assumption of transverse isotropy. A multiscale model by Washio et al. (2013) includes tissue level structure in its myocardial mechanics formulations, enabling modeling of cleavage planes. This facilitates sheet sliding, which helps to eject an additional 7.5 ml of blood compared to modeling without this level of detail. However, inclusion of realistic microstructure will hugely increase computational workload. Our fixed constraint at the base and lack of geometrical outflow tract may also impede further motion of the ventricles, thus rendering a lesser amount of myocardium usable for volume displacement, and therefore blood ejection.

Our simulations predicted a slightly higher pressure gradient within the LV (≈ 31 mmHg) than typical values reported for healthy humans (< 30 mmHg), even though the outlet boundary area is within the healthy human range (Westaby et al., 1984; Geske et al., 2011). Nevertheless, our predicted value does not exceed the severe obstruction threshold of 50 mmHg defined by Geske et al. (2011). As a consequence, the predicted peak velocity of inlet and outlet flows are also higher than normal (Mowat et al., 1983). We hypothesize that this disparity is due to our simplified generic structure with a flat ventricular base, whilst in reality the structure is slightly angled and funnel-like, which may help direct the flow (Greenbaum et al., 1981). This issue could be resolved when the framework introduced here is applied to an anatomically realistic structure for future patient-specific simulations.

The filling phase FSI has been of considerable interest in cardiac mechanics research, particularly the formation of a vortex ring and its possible link to cardiac efficiency (Hong et al., 2008). In our healthy simulation, a vortex ring formed in each of the ventricles during the early filling phase, growing in size as they approached the apex. As the filling ceased, the vortex ring moved to the basal region near the outlets. This simulated LV behavior matched the description reported in the echocardiographic study of Hong et al. (2008), which attributed the vortex to the occurrence of a shear layer between the high speed mitral jet and the low speed blood movement in the cavity, causing the mitral jet to roll up, forming a vortex. These investigators reasoned that the vortex helps to maintain kinetic energy of the blood prior to ejection.

Blood KE and |ω| have been proposed as potential markers for predicting cardiac dysfunction (Carlsson et al., 2012). The systolic and diastolic peaks of KE have been calculated clinically from 4D flow MRI data (Carlsson et al., 2012). Whilst our model's RV KE is similar in behavior to that reported by Carlsson et al. (2012), this is not the case for the LV. The clinical measurement revealed the LV KE is larger during diastole than during systole, which is not seen in our simulations. Carlsson et al. (2012) argues that two mechanisms contribute to the diastolic KE: elastic recoil of the myocardium and the displacement of the ventricular base into the position previously occupied by the atria. Since the atrial structure is not incorporated in our model, this latter effect cannot be simulated.

The KE never dropped to zero, which can be attributed to the vortex present in the chambers during diastole. The vortex has been suggested to aid energy and momentum transfer in the fluid, and the vorticity magnitude has been proposed as a measure to predict LV diastolic dysfunction caused by ventricular interdependency in cases of RV failure (Schäfer et al., 2016). Both RV and LV |ω| magnitude was largest during the filling phase as expected, due to the presence of rotating flow observed in **Figure 6**. A slight bump was observed during systole, which can be caused by the fluid being pushed against the flat fixed basal boundaries. It should be noted that the atrial contraction will generate an additional vortex, affecting |ω| (Schäfer et al., 2016).

In comparison to the passive filling vorticity measured by Schäfer et al. (2016) in healthy human LVs (≈ 30 1/s), LV |ω| simulated by the standard model was higher. We hypothesize that this discrepancy is due to filtering out the larger vorticity magnitude (> 70 1/s) during clinical image analysis, in addition to the lack of internal ventricular structures such as valves and papillary muscles in our model geometry, which could affect the vortex strength. Although both KE and |ω| did not entirely replicate realistic measurements, the details presented here could provide basic insights on how disease states affect blood energetics. More accurate measures would be obtained should a more realistic structure be implemented.

# 4.2. LVAD Simulation

We simulated the success of the LVAD in restoring the aortic pressure in a failing heart. An immediate reduction in the native aortic flow rate was also observed. Moreover, our computational simulations demonstrated that an imbalance of pump flow rate and LV filling rate can result in ventricular collapse, a phenomenon known as suction, which can impede the pump's performance. As such, our modeling framework can be applied as a simulation tool to optimize LVAD design, as well as predict and investigate any clinical risks for various heart conditions prior to clinical studies, considering LVAD is designated as high risk (Class III) device by the US FDA.

By simply adjusting the geometry and myocardial contractile strength parameter, kTa, the standard model was modified to simulate typical characteristics of LVAD recipients. The diseased heart model's LV EF (16.5%), LV internal diameter (7.4 cm) and LV EDV (299 ml) are within the ranges observed in severe dilated heart failure patients (Dandel et al., 2008). Therefore, the standard model can be easily modified to simulate behaviors of relevant diseases.

### 4.2.1. Impact of LVAD on Hemodynamics

Computational electrical-FSI simulations offer advantages over most imaging modalities to investigate fluid velocity profiles in the presence of LVADs. This is due to their superior resolution and no risk of image artifacts brought by the metallic implants (Carr et al., 2010). Analysis of the fluid streamlines can be beneficial to determine risk of thrombosis associated with cannula design and placement, as was performed in the 2D modeling study of Ong et al. (2013).

Following LVAD placement, our model displayed an immediate increase in aortic pressure due to the additional input on the systemic circulation by the LVAD outflow. As such, higher LV pressure was needed to enable opening of the aortic outlet. However, with a weaker contractility, the LV was not capable of generating sufficient pressure to overcome the restored aortic pressure. This explains the reduction in the aortic flow rate and aortic ejection duration. With further increase in pump speed, the aortic pressure rose more, resulting in a complete cessation of aortic ejection. Aortic valve closure is a known complication in many LVAD recipients as it may result in valve fusion, thus requiring speed adjustment as a prevention (Rose et al., 2000). Therefore, using computational simulations under patient-specific settings, it is possible to determine whether a safe window of pump speed operation exists for a certain ventricular state prior to clinical implantation.

Our simulations predicted separation of pressure waveforms measured at the LV base and cannula as the pump was maintained at a higher pump speed. Such behavior has been observed in a canine experiment where separation of pressure waveforms was associated with suction or overpumping characteristics (Salamonsen et al., 2015). Following implantation, the normally rectangular LV pressure-volume loop transformed into a more triangular shape. This was due to the loss of isovolumic phases, since the LVAD model is based on a continuous flow pump. The triangular shape became more prominent at a higher pump speed as the LV pressure generation capability was negated by the stronger LVAD flow rate.

Under reduced LV filling induced by mitral stenosis, the LV pressure-volume loop was severely depressed because the pump flow rate was larger than the inflow rate supplied by the left atrium. Therefore, the LVAD emptied the LV, causing it to collapse. Consequently, native heart pumping capability was impaired, as predicted by the end-systolic pressure volume relation (ESPVR), which linearly relates the maximum pressure that can be generated with the LV volume. During LV collapse following mitral occlusion, the LV is held at an extremely low volume, so, according to the ESPVR it is expected that the generated LV pressure is low (Guyton and Hall, 2006).

Increasing the LVAD pump speed appeared to result in diminished diastolic vortex formation, as shown by the fluid streamlines in **Figure 8**. This is further aggravated in the collapsed LV induced by mitral stenosis. We speculate that this change can be attributed to the LV geometry becoming less ellipsoidal, which could reduce space for the vortex to form. This effect is especially exaggerated in the collapsed LV.

### 4.2.2. Impact of LVAD on Myocardial Mechanics

Following the increase in LVAD pump speed, our simulations revealed that the stress level experienced by the LV wall gradually dropped from cycle to cycle. The stress-strain loops assessed along the fiber direction exhibited a leftward trend, and the local strain was reduced, which are the result of LV decongestion by the LVAD.

On the contrary, our model predicted only a small drop in RV stress level following LVAD intervention, whilst the RV strain components were unchanged. This is similar to predictions by the modeling study of Sack et al. (2016) that RV stresses remain unaltered following LVAD intervention. The improvement in the RV PV loop under LVAD support, shown by our simulation, may only be due to the alteration in septal stresses and strains due to the unloaded LV. As such, only a slight PV improvement was noted and the RV free wall mechanical work remained unchanged.

When LVAD support was simulated in the presence of mitral stenosis, the LV was substantially drained and all LV stress-strain tensor components reduced significantly, indicating a severe reduction of LV mechanical regional work. The LV chamber lost its ellipsoidal geometry and its diameter was reduced significantly. At the apex, mechanical contact was observed between the endocardium and the cannula. If such a suction event is not alleviated, myocardial injury may occur and could trigger arrhythmia (Vollkron et al., 2007).

### 4.2.3. Behavior of LVAD Variables

In our simulations, pump motor current and flow rate waveforms approximately followed the LV pressure waveform. During contraction, pump flow rate increases due to reduced pressure difference between the aorta and the LV that needs to be overcome by the pump. The increase in pump flow results from the increased hydraulic loading on the pump impeller, as noted by the torque transfer formulation in Equation (27) (Lim et al., 2010). Therefore, the pump motor current has to increase as well to generate more torque to overcome the increased hydraulic loading.

At higher pump speed, the simulation predicted loss of pulsatility (**Figure 10**), which coincided with the loss of peak ventricular pressure. Since the pump produces stronger torque at higher speed, the increased hydraulic loading can be easily overcome by the pump without requiring an extensive increase in pump current during systole. The collapsed LV imposed greater resistance to the pump flow rate, causing the pump flow pulsatility to be severely weakened.

Pump current waveforms have been used as a measure to detect suction events in the heart under LVAD support (Yuhki et al., 1999). As the pump current waveform shows good agreement with the LV pressure waveform, this measure can be used to detect suction events, heart rate and aortic valve closure. The pump current is a non-invasive measure that is readily available in the pump. Furthermore, it does not suffer baseline drifting, typically experienced by pressure sensors (Troughton et al., 2011). Our biventricular model can be applied as a simulation tool to aid design and test better controllers capable of detecting these events and providing counter measures, such as lowering the pump speed or alerting medical personnel (Stevens et al., 2014; Robertson et al., 2017), before the patient situation gets worse.

Overall, we observed a stable pump speed throughout the simulation, with the pump speed was nearly similar to the target speed, ωset as shown in **Figure 10F**. Whilst Lim et al. (2010) employed a proportional-integral (PI) controller to maintain the speed, we did not model the integral component. Nevertheless, we did not observe any major discrepancy in the pump speed relative to ωset, indicating a sufficient controller function for the simulation.

### 4.2.4. Impact of LVAD on the RV

Our simulations predicted a small improvement in the RV function with LVAD support, in particular a stroke volume increase and a drop in peak RV pressure to a healthy level. Further improvements were observed with an increase in pump speed. These improvements may be attributed to increased RV filling subsequent to LVAD-enhanced LV output, as well as the direct result of reduced afterload (Morgan et al., 2013). As the LV decompressed, LV filling improved and subsequently reduced the pulmonary arterial pressure. With reduced afterload, RV ejection fraction is expected to improve according to the end-systolic pressure volume relation (Maughan et al., 1984).

On the contrary, when mitral stenosis was induced in the simulation, RV dilation and increased peak RV pressure reappeared. The mitral stenosis essentially increased RV afterload, preventing more blood from leaving the pulmonary circulation. As such, RV pressure increased to counter the increased afterload. With the baseline pump flow rate barely changed and the increased afterload impeding RV outflow, the RV began to dilate. This helped push the septum further into the LV chamber as the LV collapsed.

Clinically there is considerable interest in the impact of LVADs on the RV due to the significant pool of LVAD recipients developing RV failure (Neyer et al., 2016). Nonetheless, a clinical observation by Morgan et al. (2013) noted that LVAD actually improved RV function in some patients. RV failure may be caused by an undetected RV disease that is masked by more significant LV failure. LVAD support improves the venous return to the RV; however this may aggravate the undetected RV disease as a weaker RV will be unable to cope with the improved preload. This could be investigated in future with this model to optimize LVAD settings prior to implantation.

To our knowledge, previous finite element models of heart-LVAD interaction did not consider the biventricular interaction. For example, the study by McCormick et al. (2013) expands the LV FSI model of Nordsletten et al. (2011) by adding an immersed cannula structure, but does not implement any electrophysiology physics. Other studies by Sack et al. (2016) and Heikhmakhtiar et al. (2017) simply modeled the LVAD contribution by adding a constant flow from the heart and into the aortic circuit of the four-chamber Living Heart Model (Baillargeon et al., 2014) and the Gurev et al. (2011) biventricular model, respectively. Nevertheless, the cannula itself and the fluid dynamics were not modeled in these studies.

The framework presented here offers a number of advantages over previous heart-LVAD models. In particular the capability to simulate

Considering the LVAD is classified as high risk class III medical device by the FDA, we believe computational modeling via the framework introduced here can help improve device safety and efficacy prior to pre-clinical studies and clinical testing in patients with severe heart disease. Furthermore, some devices, such as the pediatric Berlin Heart LVAD, have also been approved for rare medical cases through the Humanitarian Device Exemption (HDE) rule, which relaxes the clinical efficacy requirement due to the lack of test subjects (Almond et al., 2011). Computational simulations can help fill this gap and provide a higher degree of confidence among clinicians and regulatory bodies to approve the usage of this high risk device in this cohort of patients.

With all major cardiac physics present in this model, simulations can be easily expanded for other treatments and diseases. The inclusion of basic fluid hemodynamics allows the model to simulate internal flow characteristics of the heart with various implants such as artificial valves, which cannot be presented in electromechanical simulations such as Kerckhoffs et al. (2009). On the other hand, FSI simulations with spatiallyuniform contraction stress, such as Krittian et al. (2010), cannot be extended to include electrical abnormalities. Inclusion of a closed-loop circulation and biventricular structure also provide the capability to study biventricular interaction, which cannot be simulated in the previous LV fluid-electromechanics studies (Watanabe et al., 2004). As such, the model presented here provides a framework for more intensive and extensive cardiac simulations.

# 4.3. Limitations of Current Modeling Framework

A morphologically-realistic geometry was not implemented as the main aim of this study was to demonstrate the model's capability to simulate cardiac multiphysics phenomena. An idealized geometry eases the computational load, making it more suitable for framework development. We consider this study a stepping stone for future simulations based on patient-specific anatomies, extracted from imaging modalities such as computed tomography (CT).

In terms of electrophysiology, the model presently does not include several mechanofeedback components such as the stretch-activated currents and stretch effect on the cellular membrane's caveolae. Stretch-activated channels have been proposed to be capable of triggering electrical activation via mechanical stretch alone (Sachs, 2010). However, more characterization needs to be performed to better understand the channel kinetics. The effect of stretch-activated channels can be easily added to the source term in Equation (1) in the future. Stretch may also reduce the number of caveolae structures, leading to a reduction in membrane area and hence reduced membrane capacitance. This effect has been found to slow electric conduction in the cardiac myocyte (Pfeiffer et al., 2014). Thus, this suggests that the membrane capacitance must be made stretch-dependent if such characteristics are to be modeled. Nevertheless, it should be noted that the addition of these stretch-dependencies may increase the non-linearity of the electrical physics. It is likely that these additions may not alter the predictions of most organ-level simulations as the electrical to mechanical interaction is the more crucial coupling.

The use of phenomenological AP and active stress formulations provided considerable simplification, since it reduced the number of variables required to be solved relative to biophysically-based models. This also enabled the use of larger time step steps, for more efficient computation. Mechanical dependencies of the contraction strength may be necessary for specific future studies to simulate the Frank-Starling mechanism, which was not considered in the present model. The mechanism enables the heart to pump at a greater ejection fraction should it be filled with a larger amount of blood (Guyton and Hall, 2006). Nevertheless in failing heart, the Frank-Starling mechanism is known to be absent (Schwinger et al., 1994). With mechanical dependencies, numerical instability may arise especially if a segregated-type solver is implemented (Niederer and Smith, 2008). While a fully coupled-type solvers as adapted here and elsewhere (Göktepe and Kuhl, 2010)are more robust, they tend to require larger computer processing power and memory. As such, additional tuning of the solver settings may be necessary.

Passive myocardial mechanical properties were simply assumed to be a transverse isotropic hyperelastic, whilst in reality, the myocardium exhibits orthotrophic properties (Dokos et al., 2002). Nonetheless, the transverse isotropic hyperelastic material should be able to replicate the uniaxial and biaxial mechanical testing by Demer and Yin (1983). However, if realistic shear characteristics are desired, the orthotrophic terms of Holzapfel and Ogden (2009) can be simply added to the existing strainenergy function.

Viscoelasticity is a common property of myocardial tissue, contributed by the collagenous structure and intercellular fluids, helping to damp out vibration, due to the blood flow (Cansiz et al., 2017; Quarteroni et al., 2017). In the model presented here, a Rayleigh damping component was added to damp out oscillations, where the Rayleigh damping parameters were set proportional to the mass and stiffness terms in the equation of motion. In future, incorporation of the viscoelastic components would be a more accurate way to damp out such oscillation without the need for Rayleigh damping.

Active atrial contribution to the filling phase was not considered in the present study. The atrial contraction is expected to produce an additional peak in the filling flow rate waveform, which forms the "A" wave commonly observed in Doppler imaging of the flow waveforms (Sohn et al., 1997). This atrial contraction will form an additional vortex following the early passive filling (E-wave), and will thus affect the blood kinetic energy and average vorticity plots Hong et al. (2008). The "A" waveform could be simulated by adding an active atrial function in the Windkessel atrial compliances to simulate their physiological pressure-volume relations. In addition, filling kinetic energy is dependent on the basal displacement into the spatial location, previously occupied by the atria (Carlsson et al., 2012). However, incorporating this effect into the simulation would require structural modeling of the atria as well.

Valve structure will influence the intra-ventricular blood flow profile. Furthermore, the angled position of the valves may also ease flow resistance during filling and ejection. Among LVAD recipients, tricuspid regurgitation has also been noted, in particular in those who developed RV failure, since the RV may not be able to cope with the restored systemic venous return (Hayek et al., 2014). Furthermore, aortic valve fusion is a known issue, experienced among LVAD recipients (Rose et al., 2000). Consequently, incorporating the valve structure mechanics into the multiphysics models may be necessary in future due to the high prevalence of valve failure in many ventricular diseases. Inclusion of valvular structure will likely require contact modeling, as was performed between the endocardium and cannula, as well as employing a re-meshing algorithm to overcome severe mesh distortion. Although we used laminar Navier-Stokes formulations, the turbulent form may need to be considered in future, especially in the presence of structural valves. Turbulence models have rarely been considered in multiphysics simulations, since they are more computationally expensive, mainly applied in pure fluid simulation studies (Chnafa et al., 2015).

The model mesh we have used was selected as a compromise between global convergence and computational cost. Fluid dynamics, especially during the diastolic phase, can involve more complex flow patterns, due to vortex formation and the presence of low velocity regions; these may require a finer mesh to be accurately modeled. With the present mesh setting, a converged solution was obtained for blood flow systolic characteristics, but the convergence slightly worsens during the isovolumic phases and diastolic filling. Nevertheless, qualitative and global measures such as vortex shape and flow waveforms can be sufficiently obtained with present mesh settings. In addition, changes to action potential parameters such as σ and k2, resulting in increased conduction velocity, will require finer meshing to achieve a converged conduction velocity and to prevent numerical convergence failure at large element size. Therefore, future mesh settings need to be tailored to the requirements of the study.

# 5. CONCLUSION

The cardiac multiphysics framework developed was applied in a biventricular structure with an idealized Purkinje fiber network utilizing a modified form of the standard cable equation, taking into account the change in Purkinje fiber radius. The Windkessel circulation was expanded by including a closedloop circulation which linked both ventricular chambers. The model managed to simulate healthy LV and RV electrical activation sequences, mechanical behavior as well as global hemodynamics. This highlights the model's capability for future biventricular modeling work. The model was also tested to simulate LVAD support and impact of the pump on ventricular function. The model was shown to be capable of predicting changes in pump variables following changes in the state of the heart, possibly aiding pump controller design in future. The fluid-electromechanics model enables study of the link between electrical activation, myocardial mechanics and blood hemodynamics characteristics, allowing holistic simulations of new therapeutic approaches in their ideation stage, aiding to speed up technology transfer to pre-clinical and clinical trial stages.

# AUTHOR CONTRIBUTIONS

AAB carried out the computational modeling. AAB, AA, and SD developed the framework for the standard model. AAB and AA drafted the manuscript. AAB, AA, MS, NL, and SD developed the framework for LVAD simulations, conceived of and coordinated the study, reviewed the manuscript and gave final approval for publication.

# FUNDING

AAB is funded by a University of New South Wales (UNSW Sydney) Postgraduate Research Scholarship.

# ACKNOWLEDGMENTS

This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. We would also like to thank the High Performance Computing unit at UNSW for providing computational resources and assistance.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.01259/full#supplementary-material

# REFERENCES

ventricular pressure–strain loop area: a non-invasive index of myocardial work. Eur. Heart J. 33, 724–733. doi: 10.1093/eurheartj/ehs016

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Ahmad Bakir, Al Abed, Stevens, Lovell and Dokos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Advancing Regulatory Science With Computational Modeling for Medical Devices at the FDA's Office of Science and Engineering Laboratories

### Tina M. Morrison\*, Pras Pathmanathan, Mariam Adwan and Edward Margerrison

*Office of Science and Engineering Laboratories, Center for Devices and Radiological Health, U.S. Food and Drug Administration, Silver Spring, MD, United States*

### Edited by:

*Markus Reiterer, Medtronic, United States*

### Reviewed by:

*Marco Viceconti, Insigneo Institute of in silico Medicine, United Kingdom Adam Himes, Medtronic, United States*

> \*Correspondence: *Tina M. Morrison*

<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Translational Medicine, a section of the journal Frontiers in Medicine*

Received: *30 April 2018* Accepted: *08 August 2018* Published: *25 September 2018*

### Citation:

*Morrison TM, Pathmanathan P, Adwan M and Margerrison E (2018) Advancing Regulatory Science With Computational Modeling for Medical Devices at the FDA's Office of Science and Engineering Laboratories. Front. Med. 5:241. doi: 10.3389/fmed.2018.00241* Protecting and promoting public health is the mission of the U.S. Food and Drug Administration (FDA). FDA's Center for Devices and Radiological Health (CDRH), which regulates medical devices marketed in the U.S., envisions itself as the world's leader in medical device innovation and regulatory science–the development of new methods, standards, and approaches to assess the safety, efficacy, quality, and performance of medical devices. Traditionally, bench testing, animal studies, and clinical trials have been the main sources of evidence for getting medical devices on the market in the U.S. In recent years, however, computational modeling has become an increasingly powerful tool for evaluating medical devices, complementing bench, animal and clinical methods. Moreover, computational modeling methods are increasingly being used within software platforms, serving as clinical decision support tools, and are being embedded in medical devices. Because of its reach and huge potential, computational modeling has been identified as a priority by CDRH, and indeed by FDA's leadership. Therefore, the Office of Science and Engineering Laboratories (OSEL)—the research arm of CDRH—has committed significant resources to transforming computational modeling from a valuable scientific tool to a valuable regulatory tool, and developing mechanisms to rely more on digital evidence in place of other evidence. This article introduces the role of computational modeling for medical devices, describes OSEL's ongoing research, and overviews how evidence from computational modeling (i.e., digital evidence) has been used in regulatory submissions by industry to CDRH in recent years. It concludes by discussing the potential future role for computational modeling and digital evidence in medical devices.

Keywords: medical devices, computational modeling, regulatory science, virtual patients, virtual clinical trials, FDA

# INTRODUCTION

The mission of the U. S. Food and Drug Administration (FDA) is to protect and promote public health, and it does so by ensuring the safety, effectiveness and security of FDA-regulated products1 These products include, but are not limited to, medical devices, drugs for humans and animals, and biological products such as vaccines and the blood supply, each of which are managed by separate Centers within the Agency. The FDA accomplishes its mission by performing pre-market clearance, approval and post-market monitoring of the performance and safety of products, enforcing, and ensuring compliance to manufacturing processes and quality control, and conducting regulatory science research. The latter, although less well-known in the scientific community, is fundamental to support sciencebased regulatory decision-making by FDA. Regulatory science encompasses the development of new methods, standards, and approaches to assess the safety, efficacy, quality, and performance of FDA-regulated products and products under development. Each Center in the FDA is committed to advancing these efforts, which have accelerated the product development pathway and regulatory review cycle so that new, innovative products can be made available to the American public.

The FDA faces many challenges (1), such as new and evolving public health threats; rapid scientific breakthroughs and emerging technologies resulting in novel products that may raise unique testing and safety issues; globalization of public health, science, manufacturing and supply chains; and providing timely, accurate and useful consumer information in an age of information overload. To enable the Agency to meet today's public health needs and to be fully prepared for the challenges and opportunities of tomorrow, FDA leadership developed a strategic plan identifying nine target areas, stating that investment in these areas is essential to mission success2 Of those nine, four priority areas identified an important role for computational modeling3 , see **Table 1**. These priorities also have relevant aspects related to medical devices4 . regulated by the Center for Devices and Radiological Health (CDRH)5 , as mentioned by the FDA Commissioner in a blog posted in July 20176 .

CDRH's mission goes beyond protecting and public health; with a vision to be the world's leader in medical device innovation, they provide consumers, patients, their caregivers, and providers with understandable and accessible science-based information about the products it oversees, and facilitate innovation by advancing regulatory science. Science-based regulatory decisions are made with evidence collected from four different models: animal, bench, computational7 , and human (i.e., clinical trials), see **Figure 1A**. While each model has its advantages and limitations for evaluating different aspects of medical device performance (3), computational modeling is a promising one for supporting the future of medical devices and healthcare. FDA's Office of Science and Engineering Laboratories (OSEL) has committed significant resources for transforming computational modeling from a valuable scientific tool to a valuable regulatory tool because of its potential for significant cost-savings in evaluating medical devices, simulating performance under scenarios that may not be possible with human use or that could more effectively be evaluated with simulation.

OSEL has a unique role in medical device regulation serving as the research arm for CDRH. OSEL's expertise spans a variety of scientific, engineering, and mathematical disciplines8 , with a diverse group of 130 full time scientists and engineers (supported by numerous post-doctoral fellows and interns) that provide expert support internally to the regulatory teams and externally to industry, clinical and the scientific communities. They conduct cutting-edge research, ensure readiness for emerging and innovative medical technologies, develop evaluation strategies and testing standards, create accessible and understandable public health information, deliver timely decisions for products across their life cycle, and readily share data and engage with stakeholders to advance regulatory science. The growing area of computational modeling is fully supported by OSEL and CDRH through research and development of methods and tools, serving as expert consultants by reviewing and assessing computational modeling submitted by medical device companies, i.e., sponsors, in regulatory submissions, and publicly sharing computational modeling that supports regulatory decision-making. Moreover, OSEL houses a high-performance computing center which supports scientific computing needs for CDRH and other Centers across FDA9

Because OSEL scientists have different roles to support CDRH, we designed and conducted a 35-question survey to better understand computational modeling in research and regulatory domains, including goals and objectives on the use of computational modeling in research, reliance of evidence from computational modeling and simulation (i.e., digital evidence) in regulatory submissions, and opportunities for the future with computational modeling and simulation. Thirty-six OSEL scientists with direct involvement in computational modeling projects and initiatives were

1https://www.fda.gov

2http://www.fda.gov/ScienceResearch/SpecialTopics/RegulatoryScience/

3Computational modeling is the process of representing a real-world system by means of a computer and then running the simulation by implementing a numerical scheme.

4The FDA definition of medical devices excludes drugs, which achieve their effects through chemical action within or on the body. The vast majority of medical device manufacturers are classified by the U.S. Department of Commerce in five industries: x-ray and electromedical equipment, surgical and medical instruments, surgical appliances and supplies, dental equipment, and ophthalmic goods (2) 5https://www.fda.gov/AboutFDA/CentersOffices/

OfficeofMedicalProductsandTobacco/CDRH/CDRHReports/ucm274152.htm 6https://blogs.fda.gov/FDAvoice/index.php, post July 7, 2017.

7Evidence from algorithms, and computer-based modeling and simulation will be referred to as digital evidence.

8https://www.fda.gov/AboutFDA/CentersOffices/

OfficeofMedicalProductsandTobacco/CDRH/CDRHOffices/ucm115989.htm 9https://scl-wiki-01.fda.gov/wiki/index.php/Main\_Page

TABLE 1 | In 2011, FDA identified an important role for computational modeling in its strategic priorities.

*FDA's Center for Devices and Radiological Health also published a special report on regulatory science*5 *to align with these efforts. The table highlights the priority areas on the left, the medical device relevance in the middle, and the proposed methods and approaches on the right.*

*aNote that (Q)SAR models are classification models that relate the structure of a chemical to its activity, i.e., quantitative structure activity relationship.*

interviewed, their responses were transcribed in a database and then shared with the scientists for fact checking. This perspective will present the results from that survey and highlight the different roles that computational modeling has and can play in medical devices, and discuss the potential future for digital evidence and simulation in medical devices.

# OVERVIEW OF COMPUTATIONAL MODELING FOR MEDICAL DEVICES

Computational modeling can be used to simulate and better understand medical devices in several ways, as depicted in **Figure 1B**. Starting with the upper row, the simplest and most common implementation of computational modeling for medical devices is simply to simulate the device under a variety of conditions that mimic some aspect of the clinical or use environment to investigate some aspect of the device's performance. Computational modeling applications also include simulating the anatomy or serving as computational human phantoms for medical imaging systems or as a platform for assessing implanted devices; simulating physiology, such as electrophysiology during arrhythmias in the heart, or of pancreatic function; simulating chemical toxicology (using (Q)SAR models), which can support our ability to understand whether or not compounds released from medical devices, such as from dyes and coatings, are harmful; simulating the additive manufacturing process to optimize a 3D-printed product or simulating the substrate on which a 3D-printed product will be manufactured. The applications in the upper row of **Figure 1B** typically support design or evaluation of a physical, medical device. Other applications, depicted in the lower row, include computational algorithms embedded in a medical device or serving as the medical device, i.e., software as a medical device (4). An example of the former is embedded control algorithms in glucose monitors, which have the potential for advancing modern artificial pancreas systems10 used in glucose regulation for patients with diabetes (5). The models for the artificial pancreas have been used to replace in vivo animal studies to initiate clinical studies for these closed-loop devices (6). An example of computational modeling as a medical device is the use of personalized simulation to indicate whether a patient is a candidate for a medical device or a pharmaceutical, for example, to simulate an invasive clinical procedure or dosage effect to predict an outcome before the therapy is selected.

Computational modeling can also be used to simulate treatment outcomes. Statistical models have long been used to simulate clinical trial design and interpret results. An evolving concept is that of "virtual patients," and new statistical models to augment clinical trial design with virtual patients to predict treatment outcomes (7, 8). It is important to note that a "virtual patient" is not necessarily a digitized patient; it is an approach that allows previously collected evidence (such as digital evidence or other historical clinical evidence typically referred to as "external evidence") to inform the collection of new evidence from a clinical trial using Bayesian methodologies. Thus, computational modeling can enable a pathway to expose fewer patients to experimental therapies by relying on other sources of evidence. It can offer an opportunity to address questions that we cannot address clinically due to financial or ethical considerations, and investigate aspects of device performance in many more clinically-relevant cases (hundreds of

10https://www.fda.gov/MedicalDevices/ProductsandMedicalProcedures/ HomeHealthandConsumer/ConsumerProducts/ArtificialPancreas/default. htm

computational, and human (i.e., clinical trials). (B) Computational modeling has the potential to transform medical device design and evaluation in several ways. The upper row consists of applications that typically support the design or evaluation of the physical device. The lower row represents other applications, such as those embedded in a device or simulation as a medical device. Moreover, computational modeling can also simulate treatment outcomes or simulate the clinical trial for imaging systems. Lastly, it can play a critical role in the development of data-driven models from real-world data. See the text for more details.

thousands as compared to hundreds). Computational modeling can facilitate the exploration of using a medical device in populations that cannot be investigated clinically, such as in patients with rare diseases or pediatric patients, without harm. Computational modeling has also enabled the complete "in silico" simulation of clinical trials for medical imaging systems. By this we mean the implementation of different computational models to simulate the entire clinical evaluation of an imaging system, creating a "virtual clinical trial," where no patients are physically exposed to the imaging system–more on this later. Lastly, knowledge-base tools can be harnessed to develop data-driven models from big data sources, such as real-world data, and employ deep learning methods to gain relevant insights about medical device use and performance. Computational modeling for medical devices has a broad scope impacting many facets of the product lifecycle, and scientists from OSEL are leading or closely collaborating with leaders in the field in each of the aforementioned categories.

# COMPUTATIONAL MODELING RESEARCH

The research conducted in OSEL is directly motivated by regulatory needs, scientific questions arising from the review of regulatory submissions and anticipating future direction of industry needs through technology forecasting. Moreover, the vision of CDRH is to ensure patients in the U.S. have access to innovative medical devices first in the world, and computational modeling is one tool to support faster more efficient regulatory approvals without sacrificing patient safety or the confidence in regulatory decisions. Some companies have stated that the cost for clinical trials may soon outpace revenue (9), and industry will therefore need other relevant and reliable data sources for demonstrating safety and effectiveness of medical devices; computational modeling is a practical and viable method for gathering clinical information to augment clinical trials (10). More details on this are provided in the closing section.

There is a broad range of modeling disciplines that OSEL scientists are using in their medical device-driven research, including photon transport, fluid dynamics, heat transfer, electromagnetism, solid mechanics, acoustics and optics, along with anatomical, physiological, and mechanistic modeling. Other include (Q)SAR models for assessing molecular carcinogenicity (11), deep learning methods and artificial intelligence for analyzing and synthesizing real-world data. Within this diverse range, OSEL has been advancing different areas of computational modeling for medical devices. The following examples provide a glimpse of the many computational modeling applications in OSEL.

Scientists have developed computational models of medical devices for investigating a specific approach or consideration about the medical device. It is important to emphasize that the computational studies below are not of a specific manufacturer's device, but of generic devices where the study results have broad impact in that device domain and are translatable to other domains. These models include implantable cardiovascular stents for assessing different methods to calculate fatigue safety factor (12); heart valves implanted with non-circular configurations (13) to assess the impact on stresses and strains; inferior vena cava filters to demonstrate a new method for computing embolus transport (14); hip implants for evaluating the impact of the design on contact mechanics (15); radiofrequency coils for MRI systems (16, 17) to investigate the design parameters on the electromagnetic field; surgical facemasks (18) for evaluating aerosol leakage of different designs; blood pump (19) for assessing the ability to predict hemolysis using computational fluid dynamics; and electrical stimulation of implanted lead wires (20) to investigate local heating. They have also developed new methods for simulating photon transport of x-ray emitters (21) and compressive sensing for imaging systems (22). Another computational effort was the development of a complex constitutive models for absorbable polymers used in medical implants (23).

Computational models of anatomy or physiology include improved drug delivery in the cornea with ultrasound energy (24); physiological models of heart cells (25), renal circulation (26), hemodynamic responses to blood volume perturbations (27), left bundle branch block (28), gas dynamics in the retina (29), coupled electrical and mechanical activity in the heart (30); energy absorption in patients with deep-brain stimulators (31– 34), breast tissue expanders (35), in pregnant women and fetus during MRI exams (36); subthalamic nucleus (37), the breast (38), cancellous bone (39), the head (40) and whole body models (41, 42).

A part of OSEL's mission is to improve CDRH's ability to evaluate medical devices and support the regulatory approval of innovative medical devices more efficiently without sacrificing safety. Therefore, some of the research efforts involve demonstrating through examples that the output from computational modeling is a viable source of regulatory-grade evidence, i.e., sufficiently-credible digital evidence that can support regulatory applications (43–45); developing frameworks (46, 47) and metrics (48) for assessing the trustworthiness of models, and studying workflows for creating reproducible models11 (49), and identifying considerations for computational patient models for autonomous medical devices (50).

Other computational tools to assess specific aspects of device performance or safety that industry can employ include a simulator for high-intensity focused ultrasound (HIFU) beams and heating effects (51, 52), benchmarks models for computational fluid dynamics (19), patient-specific workflows for assessing clot trapping efficiency in IVC filters (53), surrogate models for predicting device-specific and speciesspecific hemolysis (Craven et al., under review), optical-thermal light-tissue interactions for photoacoustic breast imaging (54), and an online app for assessing the safety of color additives (55).

Additional efforts are pushing the state of the art of simulation for medical devices, including fluid-structure-interaction of deformable blood clots, computational human phantoms for active implants (37), lesion insertion and image reconstruction (56), computational patient models for closed-loop control devices (27), whole-heart modeling for electrophysiology devices (30), computational modeling for determining hemolysis levels in patients supported by blood-circulating medical devices (57), evaluating exposure risk from nickel leaching devices (58) and risk assessment for framing policy and deciding on the stockpile of personal protective equipment for wide-spread outbreaks or virus epidemics (59).

Lastly, as previously mentioned, one team is developing and validating a framework for streamlining the market entry of imaging systems relying solely on simulation in place of clinical trials. The VICTRE project (virtual imaging of clinical trials for regulatory evaluation) approach involves simulating the anatomical structure of the breast (with or without a neoplastic lesion), the radiological transmission (i.e., imaging

11NIBIB. R01EB024573, https://simtk.org/projects/kneehub

system) and reconstructed images, and the clinical reader studies. By simulating each component of the clinical trial process, there is the potential for minimizing the need for clinical trials and thus the regulatory review of imaging systems (60). Note that the VICTRE project uses statistical analysis tools that evaluate the diagnostic performance of radiologists (virtual or human). An important aspect of the statistical analyses is that they account for radiologist variability and case variability. Such analyses are not trivial and have been developed by OSEL scientists (61–64) with validation based on sophisticated simulation tools (65).

OSEL scientists share their models and data with the public to facilitate the use and broader adoption of these computational tools and approaches. For example, anatomical models of the head (the MIDA model) and whole body models (the Virtual Family) can be downloaded from the IT'IS Foundation's website12, 13 The experimental and computational data from an FDA-led multi-laboratory study for fluid dynamics on generic medical devices can be found here;14 the simulation and statistical tools for the VICTRE project are here15 and here (66, 67); the HIFU simulator here16 Other software applications are being shared through Github17, such as the design of a generic inferior vena cava filter18 and a risk assessment tool for assessing color additives19

### COMPUTATIONAL MODELING IN REGULATORY SUBMISSIONS

The OSEL scientists also serve as expert consultants on regulatory submissions. The review and decision about a medical device regulatory application requires a team of experts led by the regulatory offices in CDRH. OSEL scientists serve as specialized, technical experts on the regulatory teams. More than 2500 consulting reviews were completed by all scientists in OSEL in 2017, and about 500 were completed by the 36 scientists interviewed for this perspective. Of the 500 consults performed by the scientists surveyed, 220 (44%) included computational modeling and digital evidence in the submission. (Note that therefore 9% of all expert consults performed by OSEL for the regulatory offices in 2017 involved computational modeling). Of these 220, the submission-type breakdown is as follows: 36% were for premarket 510(k) notifications (for moderate risk devices), 25% for clinical trial applications, 24% for pre-submissions, 13% for premarket approval applications (for high risk devices), with only a handful for other submission types. With respect

13The Virtual Family: https://www.itis.ethz.ch/virtual-population/virtualpopulation/overview/

to medical areas, the largest number of consults in 2017 were for neurological devices, followed by cardiovascular and orthopedic devices, imaging systems, and surgical devices.

The survey results indicated that the primary use of computational modeling in regulatory submissions was to identify the appropriate bench testing configurations, such as worst-case or clinically challenging conditions, for cardiovascular, orthopedic, and surgical implants. The second most common use of computational modeling was to provide evidence supporting the safety assessment of patients with and without implanted devices when exposed to the radiofrequency (RF) fields of an MR system. A noteworthy example of the latter is the recent clearance20 of the first 7 Tesla MRI system (Siemens Magnetom)21, where the Virtual Family (41) and the MIDA head model (40) were used to predict aspects of safety and effectiveness of the new system. Examples of the former include RF safety evaluation for patients with implanted electrically passive (e.g., joint replacement, stents) or electrically active devices (e.g., neurostimulators, pacemakers, cochlear implants). Other modeling examples include therapeutic ultrasound systems where simulation results of the ultrasound energy delivered to in vivo locations have been used in regulatory submissions as justification for system parameters, or the recent clearance of Compressed Sensing GRASP-VIBE R to support high-resolution dynamic abdominal imaging under free-breathing. From the 510(k) summary22, "A comparison of the functionality was performed between the new feature and the device feature by detailed simulations with a numerical [computational] phantom".

In general, computational modeling can be part of a regulatory submission in two ways. The first is when simulation results serve as supporting (digital) evidence in a marketing application for a medical device. The second is when simulation is a medical device, such as for clinical decision support; this is "software as a medical device." Virtually all consults regarding computational modeling were of the former; the latter, with just a handful of submissions, is a new growth area for CDRH, especially with the release of the FDA guidance that describes the clinical evaluation for these software application, and the new program area on Digital Health Technologies23 Two examples of software as a medical device that have received FDA clearance used patientspecific computational models generated from CT imaging data to non-invasively predict clinically-relevant quantities for treatment selection. Heartflow R generates a personalized 3D model of the patient's coronary arteries and simulates blood flow to predict fractional flow reserve24 The CardioInsight R Mapping System generates a personalized model of the patient's heart and torso, then simulates the electrical activity on the

12The MIDA Model: https://www.itis.ethz.ch/virtual-population/regionalhuman-models/mida-model/,

14https://nciphub.org/wiki/FDA\_CFD

15https://github.com/didsr/victre

16https://www.mathworks.com/matlabcentral/fileexchange/30886-high-

intensity-focused-ultrasound-simulator.

17https://github.com/didsr

18https://github.com/kenaycock/Generic-IVC-Filter

19https://dsaylor.github.io/CHRIS/

20510(k) Summary: https://www.accessdata.fda.gov/cdrh\_docs/pdf17/K170840. pdf

21https://usa.healthcare.siemens.com/news/magnetomterrafdaclearance.html

22510(k) Summary: https://www.accessdata.fda.gov/cdrh\_docs/pdf17/K173617. pdf

23https://www.fda.gov/medicaldevices/digitalhealth/

24https://www.accessdata.fda.gov/cdrh\_docs/pdf15/K152733.pdf

heart surface from body surface potential recordings25 For more information on these and other patient-specific cardiovascular models see (68).