A comprehensive mathematical model for cardiac perfusion

The aim of this paper is to introduce a new mathematical model that simulates myocardial blood perfusion that accounts for multiscale and multiphysics features. Our model incorporates cardiac electrophysiology, active and passive mechanics, hemodynamics, valve modeling, and a multicompartment Darcy model of perfusion. We consider a fully coupled electromechanical model of the left heart that provides input for a fully coupled Navier–Stokes–Darcy model for myocardial perfusion. The fluid dynamics problem is modeled in a left heart geometry that includes large epicardial coronaries, while the multicompartment Darcy model is set in a biventricular myocardium. Using a realistic and detailed cardiac geometry, our simulations demonstrate the biophysical fidelity of our model in describing cardiac perfusion. Specifically, we successfully validate the model reliability by comparing in-silico coronary flow rates and average myocardial blood flow with clinically established values ranges reported in relevant literature. Additionally, we investigate the impact of a regurgitant aortic valve on myocardial perfusion, and our results indicate a reduction in myocardial perfusion due to blood flow taken away by the left ventricle during diastole. To the best of our knowledge, our work represents the first instance where electromechanics, hemodynamics, and perfusion are integrated into a single computational framework.


Introduction
Myocardial perfusion is the process by which oxygenated blood is delivered through the coronary arteries to the heart muscle or myocardium, enabling its oxygenation and metabolism.The microvasculature of the myocardium is responsible for facilitating the exchange of oxygen and nutrients with the blood.However, when the coronary circulation is obstructed due to factors such as arterial stenosis or cardiac pathologies like aortic regurgitation and arrhythmias, the blood supply to the cardiac muscle may be limited.This restricted blood flow can lead to ischemia and potentially trigger a myocardial infarction, commonly known as a heart attack [1].
Stress myocardial computed tomography perfusion (stress-CTP) is a method for quantitatively assessing myocardial blood perfusion through myocardial blood flow maps (MBF), obtained by exposing patients to additional radiation (with respect to standard angiography) and administering an intravenous stressor during a CT scan.In-silico computational models [2][3][4] can provide valuable insights into physiological processes and enable the simulation of virtual scenarios under multiple pathological conditions, making them useful for studying e.g.coronary by-passes [5] and ventricular hypertrophy [6].However, the development of a comprehensive model of myocardial perfusion requires accounting for the complex interactions among multiple physical processes, including the coexistence of multiple spatial scales in the coronary circulation.The coronary arterial tree can be subdivided into epicardial coronary arteries (large coronaries) and intramural vessels (arterioles, venules and microvasculature) [7].From a modeling point of view, the blood flow in the large epicardial vessels can be described using full 3D fluid dynamics or fluidstructure-interaction equations [7][8][9][10][11][12], or geometrically reduced hemodynamics model, as 1D models [13][14][15].Differently, below a threshold length scale, the blood flow in the myocardium can be represented as in a porous medium [16], thanks to Darcy or multicompartment Darcy models [4,13,[17][18][19][20].The integration of these models yields a coupled mathematical problem, featuring dynamic and kinematic conditions at the interface between large coronaries and microvasculature [17,18].
Figure 1 displays the intricate processes involved in myocardial perfusion, which result from the interplay of various physical phenomena and scales, including electrophysiology, mechanical activation, tissue mechanics, cardiac hemodynamics, and valve dynamics.In this paper, we propose for the first time a novel mathematical model that unifies these different aspects within a single framework.Our mathematical model includes core models for electrophysiology, active and passive mechanics, blood fluid dynamics in the left atrium, ventricle, and aorta, mitral and aortic valve dynamics, and myocardial blood perfusion.To partially decouple the problem, we use a fully-coupled electromechanical model to trigger a fully coupled Navier-Stokes -multi-compartment Darcy perfusion model.To the best of our knowledge, this work represents the first attempt in the literature to integrate electromechanics, fluid mechanics, and perfusion into a single computational framework.
Our computational model provides physiological coronary flow rates and myocardial blood flow maps for the healthy case, as demonstrated by our results.In addition, we successfully simulate a severe aortic valve regurgitation, which can cause reduced oxygen delivery to the myocardium due to steal of coronary flow during diastole.
Our novel integrated model is mathematically sound and physiologically accurate, as it does not require any assumptions about boundary conditions at the inlet sections of large coronary arteries -as commonly done for instance in refs.[18,[21][22][23][24] -and features a detailed 3D electro-mechano-fluid model to provide precise inputs for myocardial perfusion.Our computational model enables the simulation of the effects of various pathologies on perfusion, as demonstrated in our study on aortic regurgitation.Our work is a significant advancement towards the realization of an integrated model of the whole human heart function, which would enable in-depth investigations of physiological and pathological perfusion scenarios, including the myocardial ischemia resulting from a coronary artery occlusion.

Methods
To describe the methodology that we develop for the multiphysics simulation of cardiac perfusion, we first introduce our mathematical model, then we give details on numerical methods, software libraries, and computational setup.

Mathematical model
For the mathematical model, we consider the time domain (0, T ) and three different spatial domains: • The left heart solid domain Ω s , comprising the ventricle and atrial myocardium, and by the aortic vessel wall.In Ω s , we define the electromechanical problem.We consider a Lagrangian framework set in the reference unloaded configuration Ω s .
• The left heart fluid domain Ω f , comprising the left ventricle and atrium chambers, together with the aorta and the epicardial coronaries.In Ω f , we define the fluid dynamics problem.Ω f is a time dependent domain, but we omit the subscript t to keep the notation simpler.
• The biventricular geometry Ω p , that we model as if it were a porous medium, where we set our perfusion model.We assume Ω p to be non deformable [17][18][19].
We give a graphical representation of each domain in Figure 2, top.Notice that we ignore fluid dynamics in the right heart since coronaries originate from the left heart.Hence, there is not a direct feedback of the right heart hemodynamics on myocardial perfusion.Accordingly, also the electromechanical simulation has been performed in the left heart solely.In the following, we describe each physical problem occurring in the different domains, and provide details on the coupling conditions.The overall multiphysics model is sketched in Figure 2, bottom.

Electromechanics
To model the electric and mechanical activity of the heart, several mathematical and numerical models have been proposed in the literature [25][26][27][28][29][30][31][32].We consider the model presented in [33,34] which is set in the left heart solid domain Ω s shown in Figure 2, top.For the recovery of the reference configuration [26,[35][36][37], we refer specifically to the algorithm presented in [33].We reconstruct cardiac fibers with the Laplace-Dirichlet Rule-Based Methods [38][39][40], using the algorithms for ventricles and atria presented in [40] and [41], respectively.We assume that the active mechanics triggered by electrophysiology is present only in the left ventricle Ω s LV .Conversely, since our focus is on the dynamics downstream the aortic valve, we treat the atrial tissue as an electrically passive material.This is a simplification that has been commonly adopted in other electromechanics [31,[42][43][44] and electro-mechano-fluid [45][46][47] models of the heart.
We model electrophysiology by the evolution of the transmembrane potential v in the left ventricle via the monodomain equation [48].We denote the electrophysiology problem in compact form as where w and z are the gating variables and ionic concentrations, respectively.Note that the monodomain equation is augmented with mechano-electric feedbacks [49][50][51], as highlighted by the dependence from the solid displacement d s .We couple Equation (1) with the ten Tusscher and Panfilov ionic model [52], that we denote in short as We model the active contractile force [53] by means of the biophysically detailed RDQ20 activation model [54], which accounts for the force-sarcomere length relationship and the force-velocity relationship thanks to fiber strain-rate feedback, which we deem essential to faithfully predict blood fluxes and velocities in the CFD simulation [34,55].Denoting by s the state variables related to the active contractile force and by SL the sarcomere length, which depends on the displacement d s , we express the activation model in compact notation as where [Ca 2+ ] i represents the intracellular calcium concentration stored in the vector function z.Following [54], Equation (3) allows then to compute the active contractile force T act (s, SL).
For the structural problem, we consider the elastodynamic equation, in the unknown d s , in which the first Piola-Kirchhoff stress tensor is split into a passive term (depending on d s only) and an active term (depending on d s and T act ).For the passive part, we use the Usyk anistropic strain energy function [56].In short, we denote the structural problem as equipped with the following boundary conditions: generalized Robin boundary conditions [33] to model the action of the pericardium, homogeneous Dirichlet boundary conditions (i.e.no displacement) on the rings of the pulmonary veins and homogenous Neumann boundary conditions (i.e.no stress) on the ring of the ascending aorta.Furthermore, for simplicity, we set homogenous Dirichlet boundary conditions on the epicardial valvular ring.On the endocardium, we set the fluid pressure as described in the following paragraph.
In this work, we consider a one-way coupling between the electromechanics and the 3D fluid dynamics problems [57][58][59][60] (see below).Specifically, the 3D electromechanics problem is solved off-line, prior to the 3D fluid dynamics-perfusion problem.However, in order to account for feedback of the fluid on the electromechanical model, we strongly couple the 3D electromechanics with a 0D lumped parameter model of the circulation [33,61,62].Thus, on the endocardium, we enforce the continuity of the 0D fluid -3D solid cavity pressures and cavity volumes.Accordingly, p cavity and V cavity are the vectors collecting the pressures and volumes of the left atrium, left ventricle and ascending aorta.We denote the circulation model as where c is the state vector that includes pressures, volumes and fluxes in different compartments.Particularly, the pressure p cavity acts as a Lagrangian multiplier to enforce the volumetric costraint V cavity (d s ) = V 0D cavity [33].

The fluid geometry and fluid dynamics models
Let Ω f ⊂ R 3 be the fluid dynamics domain (that is the region occupied by the fluid) in its reference configuration (see Figure 2).The fluid domain in the current configuration is being the domain displacement (for the sake of brevity, we omit the subscript t from the fluid domain and its boundaries).The latter is computed by solving a Laplace equation in ∂ Ω f × (0, T ) with Dirichlet boundary conditions on the physical wall: d f = d f w , with d f w equal to the electromechanical displacement d s restricted on the endocardium, and zero on the coronaries walls (suitably smoothly merged, see [60,63] for further details).We compute the fluid domain velocity by u ALE = ∂d f ∂t .We compactly denote the fluid geometry problem as To model blood flows in the left heart and large epicardial coronaries, we consider the Navier-Stokes equations expressed in Arbitrary Lagrangian Eulerian (ALE) framework [64].We set our fluid dynamics problem in the domain Ω f , delimited by These boundaries represent the pulmonary veins sections, ascending aorta section, coronary outlet sections and endocardial wall, respectively (see Figure 2, top).In particular, we consider J coronary outlet sections: We denote by u f and p f the fluid velocity and pressure, respectively.We model blood as if it were a Newtonian fluid with constant density ρ = 1.06 • 10 3 kg/m 3 and dynamic viscosity µ = 3.5 • 10 −3 kg/(m s), with the total stress tensor defined as σ(u f , p f ) = −p f I + µ(∇u f + ∇ T u f ).Moreover, we account for the presence of mitral and aortic valve by means of the Resistive Immersed Implicit Surface (RIIS) method [65,66].We consider both valves as immersed surfaces in the fluid domain.Each valve is characterized by a resistance R k and by a parameter ε k , representing the half-thickness of the valve leaflets.We introduce in the momentum balance a resistive term R(u f , u ALE ) that enforces the kinematic mismatch between the relative fluid velocity and the one of the valve.We refer to refs.[6,60,65,67] for further details on this method and for the definition of R. We let the valves open and close instantaneously, at the initial and final times of isovolumetric phases (computed from the electromechanical simulation).The fluid dynamics model reads: where At the wall, we prescribe the ALE velocity (computed in Equation ( 6)).On the coronary outlets, we prescribe Robin boundary conditions, where p j c is the pressure that arises from the coupling condition with the multi-comparment Darcy model and α j are conductances, with j = 1, . . ., J [18].Moreover, on Γ f c , we also assume null tangential tractions, being τ i two tangential vectors [18].Furthermore, on Γ f pv and Γ f aa , we prescribe Neumann boundary conditions.Specifically, we set a constant physiological pressure equal to 10 mmHg on the inlet pulmonary vein sections.On the outlet section of the ascending aorta, we prescribe the systemic arterial pressure resulting from the 3D-0D eletromechanical simulation.The fluid dynamics model is also equipped with a zero velocity initial condition.

The multi-compartment Darcy model
To model blood perfusion, we consider a multi-compartment Darcy model in the biventricular myocardial domain Ω p (see Figure 2, top).This model allows us to describe several length scales featuring the myocardium and its microvasculature as a porous medium made of different compartments [4,18,19,68].Specifically, we consider the three compartments Darcy equations [18,19] in the unknown u p i , p p i , representing the Darcy velocity and pore pressure, respectively, with i = 1, 2, 3: K i is the permeability tensor, g i a volumetric sink term and the coefficients β i,k are the inter-compartment pressure-coupling coefficients.Following ref. [18], g 1 is provided by epicardial blood hemodynamics (i.e. by the coupling condition with the Navier-Stokes problem, see below) and g 2 = 0 since the second compartment does not exchange mass with the outside.Furthermore, to account for the effect of the cardiac contraction on perfusion -still avoiding the use of a poromechanical model [20] -we propose g 3 to surrogate the reservoir effect of the coronary bed by setting where γ is a suitable coefficient and the new contribution p bed (t) is a function of the left ventricular pressure p LV (t).The latter is obtained from the 3D-0D electromechanical problem (1)-( 2)-( 3)-( 4)- (5).In (9), a 1 and a 2 are calibrated in order to match physiological fluxes.The biventricular domain Ω p is partitioned into J non-overlapping perfusion regions, such that each epicardial vessel feeds only one portion [18].For the estimation of parameters K i , β i,k , with i, k = 1, 2, 3, and for the strategy employed to partion Ω p , we refer to ref. [18].

Coupling conditions
In this section, we describe the coupling conditions enforced to match the different physics.In Figure 2 bottom, we sketch the overall multiphysics model and we highlight the coupling conditions.For the fully coupled electromechanical model, we refer the reader entirely to ref. [33].
For the coupling between electromechanics and cardiac hemodynamics, we consider the following kinematic condition: where d s is defined on the atrial and ventricular endocardium and on the endothelium of the ascending aorta.We recall that, for simplicity, we set null displacement on the coronaries wall.
For the fully coupled Navier-Stokes -Darcy model, the coupling conditions read [18]: where Equation (11a) and Equation (11b) enforce the balance of internal forces and mass conservation, respectively.χ Ω pj is the characteristic function of the j-th partition [18].

Computational setup
We consider a realistic cardiac geometry provided by the Zygote solid 3D heart model [69]: an anatomically CAD model representing an average healthy human heart reconstructed from high-resolution CT scan data.We generate three meshes for the electromechanics, fluid dynamics and multicompartment Darcy problems with vmtk [70], using the methods and tools discussed in [63].Details on the generated meshes are provided in Table 1 and displayed in Figure 3a.Note that the CFD mesh is refined near the valves to accurately capture them with the RIIS method [6,60,65].Immersed valves in their open and closed configurations are displayed in Figure 3b.Notice also that we used the same mesh for electrophysiology and mechanics, with a value of the mesh size which is tipically considered too coarse to accurately resolve the traveling electrical front [48,71].However, we suitably increase the conductivities to compensate for the use of a coarse electrophysiological mesh [32,72,73].We believe that this choice is adequate for our purposes since the electromechanics simulation has the sole purpose to provide the endocardial displacement for the CFD problem.
In Figure 3c, we report the perfusion regions of the biventricular geometry: one for each terminal vessel.For the complete setup of the multicompartment Darcy model, and for the preprocessing methods used to generate the perfusion regions, we refer the interested reader to ref. [18].To surrogate the reservoir effect of the coronary bed (see Equation ( 9)), we choose a 1 = 0.4 and a 2 = 1500 Pa, which corresponds to a coronary bed pressure in the range [14.2, 61.4] mmHg.We discretize our mathematical models in space by the Finite Element (FE) method.We use linear FEs for electrophysiology and mechanics.The fluid dynamics problem is solved with linear FEs with VMS-LES stabilization [74,75], acting also as a turbulence model to account for possible transition-to-turbulence effects [76].The convective term is treated semi-implicitly.The multicompartment Darcy model, solved for the pressures, is discretized with linear FEs.As temporal advancing scheme, we use Backward Difference Formula (BDF), with the time-step sizes listed in Table 1.For additional details on numerics, we refer to refs.[18,33,60] for the electromechanics, fluid dynamics, and perfusion models, respectively.
To efficiently solve the coupled problem, we first solve the electromechanical problem using a Segregated-Intergrid-Staggered method [33,34].We pick the displacement on the fifth heart cycle -once the solution has approached a period limit cycle in terms of pressure and volume transients -and we use it as unidirectional input [60] (one-way) for the fully coupled fluid dynamics -multicompartment Darcy problem.The electromechanical displacement is linearly projected onto the CFD wall mesh.To solve the fluid dynamics -multicompartment Darcy problem, we use an implicit method with an iterative splitting strategy, i.e. we subiterate discretizations of ( 7) and ( 8) until convergence [18].We start our simulation at the end of the filling phase.We simulate two heartbeats and we report the solution of the second cycle to cancel the influence of a non-physical null velocity initial condition.
We solve the multiphysics problem in life x [77], a high-performance C++ FE library developed within the iHEART project, mainly focused on cardiac simulations, and based on the deal.II finite element core [78][79][80].Numerical simulations are run in parallel on the GALILEO100 supercomputer (528 computing nodes each 2 x CPU Intel CascadeLake 8260, with 24 cores each, 2.4 GHz, 384GB RAM) at the CINECA supercomputing center, using 288 cores.

Results
We start our analysis from a physiological simulation of a coupled electromechanicsblood dynamics -myocardial perfusion obtained by means of the proposed multiphysics model (Test I).In Figure 4, we report results from this test.Concerning electromechanics (Figure 4(a)), we show the calcium concentration, along with the displacement magnitude when the ventricle contracts.We display the intracardiac hemodynamics during filling and ejection in Figure 4(b), by reporting the volume rendering of velocity magnitude and pressure on the boundary of Ω f .Notice that the model can faithfully predict the formation of the clockwise jet in the left ventricle during filling, which redirects the blood in the aortic root during systole [81,82].Considering the cardiac chambers only, we find larger velocities during ejection, compared to the filling phase.Conversely, focusing on the coronaries only (Figure 4(b), bottom), we notice that blood is faster during the filling with respect to the ejection phase; accordingly, a larger pressure drop is also computed in the coronary tree during ventricular diastole, allowing the blood to accelerate and to perfuse the cardiac muscle.In Figure 4(c), we report the multicompartment Darcy's pressures during the filling stage.
Figure 5 shows quantitative results of the electromechanical simulation in Test I. Consistently with clinical ranges from literature [83][84][85], we compute the left ventricular stroke volume, ejection fraction, and peak pressure (the latter coming from the 0D hemodynamic model) equal to 83.0 ml, 54.2%, and 125.4 mmHg, respectively (see Figures 5a and 5b, where pressure-volume loop and volume in time are represented).From Figure 5b, it is possible to distinguish isovolumetric contraction, systolic ejection, isovolumetric relaxation, and diastolic filling phases.We show transients of the Navier-Stokes -Darcy simulation in Figure 6.We report the flow rate computed at the aortic section and the total flow rate in the pulmonary veins in Figure 6c: we compute a peak aortic flow equal to 562.0 ml/s -consistently with physiological ranges [86] -and the peak total flux in the pulmonary veins is 267.2 ml/s.In Figure 6d, we show the total fluxes computed at the outlets of the Left Coronary Artery (LCA) and Right Coronary Artery (RCA).Our mathematical model faithfully predicts a peak blood flow rate at the beginning of the filling phase (diastole), resulting from the myocardium relaxation after the systolic contraction.Our finding is consistent with clinical evidences [87]; furthermore, as also experimentally measured in [88], the flux in the LCA is larger than the one in the RCA.Pressures in the fluid dynamics domain are reported in Figure 6e.We obtain a peak systolic arterial pressure of 103.3 mmHg and a minimum diastolic pressure equal to 80.3 mmHg: both results are physiologically consistent [89].In Figure 6e, we also show the coronary pressure by averaging the average pressure in each coronary outlet: we get similar LCA and RCA  pressures.Figure 6f shows the pressure in the three Darcy's compartments.As expected, we have a decreasing pressure going from one compartment to the following one, and comparable values during the systolic peak, due to the contraction of the muscle and the consequent partial obstruction of vessels.
We aim now to study the effects that a valvular pathology has on myocardial perfusion.This allows us to explore the capabilities of the mathematical model in simulating also pathological scenarios.To this aim, we consider Test II, where the case of Aortic Regurgitation (AR) is considered.This pathology consists of a leaking of the aortic valve leaflets causing the blood to flow from the aorta to the left ventricle during the filling stage.To model the leaking of the aortic valve, we replace the "closed" physiological configuration PH used in Test I by the regurgitant configuration AR used in Test II, as we display in Figure 7a.We obtain the AR configuration by introducing a regurgitant orifice which is about the 4.5% of the aortic annulus section.Furthermore, since AR is associated with an increased systolic and a decreased diastolic aortic pressure [90], we modify the systemic arterial pressure prescribed on Γ f aa accordingly.In fact, we increase and decrease the pressure by 20% in systole and diastole, respectively (see Figure 7c).Figure 7b shows the volume rendering of the velocity magnitude in the AR case.During the filling stage, we observe reverse blood flow from the aorta to the left ventricle, yielding a mix of blood between the mitral and AR jets.In Figure 7d, we compare the coronary flowrates against time in the PH and AR cases.The diastolic flowrate decreases in the AR case, with a peak reduction of 24.8%.This trend is also confirmed by Figure 7e, where we show the velocity glyphs in the coronary tree at the diastolic peak: in the AR case, we measure much lower velocities.Differently, during ejection, we observe that the AR case produces a slight increase of the coronary flow (Figure 7d) due to a larger systemic arterial pressure than in the PH case.To better assess the consequences of this pathology in terms of myocardial perfusion, we quantify the amount of blood inside the microvascolature.Accordingly, we compute the Myocardial Blood Flow (MBF) as: MBF represents the amount of blood reaching the third compartment, i.e.where oxygen and nutrients are exchanged at a capillary level.This value is normalized over 100 ml of cardiac tissue and the factor 60 s/min allows to express the perfusion rate in minutes.Figure 7f shows a comparison between the PH and AR cases in terms of MBF at the diastolic peak.Overall, in both cases, we can observe a heterogeneous distribution of the MBF due to different resistance of the vessels inside each perfusion region, provided by the heterogeneous parameters of the Darcy model [18].More quantitatively, in the PH case, we compute an average MBF equal to 87.5 ml/min/100 ml.Our result is consistent with clinical studies, which measure a normal MBF at rest from 57.6 to 96.1 ml/min/100 ml [91].Differently, the pathological case is characterized by a much lower perfusion: at the diastolic peak we measure 68.2 ml/min/100 ml.Thus, the ventricular reverse flow due to a regurgitant aortic valve is responsible for a steal of coronary flow, and hence an abnormal and impaired myocardial perfusion.

Discussion
In this paper, we proposed for the first time a computational model to simulate myocardial perfusion accounting for the interaction of different cardiac physical processes.Our model comprises 3D electrophysiology, active and passive mechanics, blood dynamics, and myocardial perfusion, and it was successfully applied to both a healthy and an aortic regurgitant scenarios.By carrying out simulations on a realistic cardiac geometry, we showed that the model faithfully predicted electromechanics and blood dynamics quantities as previously shown also in [33,55,60].Moreover, as a new outcome of this work, we were able to predict cardiac perfusion in both physiological and aortic regurgitant cases by means of a comprehensive cardiac function model.The inclusion of the whole left heart function and geometry in our simulations allowed us to avoid any arbitrary prescription of the fluid pressure and velocities at the inlet of epicardial coronaries.Indeed, in previous perfusion models, due to the absence of any electromechanical and fluid dynamics model in the left ventricle, it was necessary to prescribe transients of flowrates or pressures at such sections [18,[21][22][23][24].Moreover, it is well known that vascular resistance increases during systole because the contraction of the myocardium compresses the intramyocardial coronary arteries [7], producing a peak flowrate during diastole.In our simulations, we can correctly capture this phenomenon without prescribing any data on the inlet of coronaries, but thanks to the interplay of different features that we included in the model, as discussed in what follows.• By including the entire left heart geometry and modeling its motion through the use of an electromechanics model, we were able to achieve a physiological blood velocity pattern in the ascending aorta.
• By modeling the aortic valve during the ejection phase, we were able to simulate partial obstruction of the coronary ostia when the valve is open.This approach resulted in a reduction of the systolic coronary flow rate, which is consistent with clinical evidence [92].Our simulations indicated that neglecting the modeling of the aortic valve in its open configuration (i.e. by setting a null resistance during systole) leads to the computation of larger and non-physiological flow rates.For the sake of brevity, we do not include the results of this case.
• The contraction of cardiomyocytes is a well-known cause of impediment to systolic coronary flow [93,94].To account for this effect in our perfusion model, we introduced a novel time-dependent coronary bed pressure (refer to Equation 9) that emulates an increase in vascular resistance [92], thereby enabling the simulation of a diastolic coronary flow rate.
Furthermore, the development of a multiphysics mathematical model allowed us to investigate that a regurgitant aortic valve produces a reduction of coronary flow during diastole, by redirecting the blood in the left ventricle, as highlighted in Figures 7b and  7d.The main consequence of this aspect is a reduced perfusion of the myocardium during diastole, accordingly with clinical evidence [95] and quantified by the computation of a reduced MBF at the diastolic peak (see Figure 7f).Furthermore, we faithfully captured also a slight increase of the epicardial coronary flow during ejection with respect to the physiological case, as described in [95,96] (see Figure 7d, blue lines during ejection).This is due to a larger aortic pressure during systole, with respect to the physiological case.In addition, the clinically meaningful outcome of the regurgitant simulation permitted us the prove the robustness of the proposed model with respect to geometric variations.
The study presented in this work features some limitations.We highlight that setting a null displacement on the epicardial valvular ring in the electromechanical simulation is not physiological.Indeed, the base of the ventricle should be free to move up and down during the cardiac cycle.Accordingly, the coronaries should follow this motion.Moreover, they should be compliant, whereas we have assumed here that they are rigid and static.
We noticed that the coronary pressure during ejection found by our simulation in subject PH is smaller if compared with standard physiological values and other computational studies [8,97].In particular, from Figure 6c, the LCA and RCA pressures should be about 25 mmHg greater during ejection.We believe that this is due to the open aortic valve configuration which is representative and not obtained by an FSI simulation.This limitation may result in an excessive occlusion of the coronary ostia and then in an augmented resistance during systole, which provokes a decrease of the computed coronary pressure.

Figure 1 :
Figure 1: The perfusion is the result of complex interactions among different models.

Figure 2 :
Figure 2: Top: Computational domains; Bottom: Coupling of the different physics.

Figure 3 :
Figure 3: (a) The three computational meshes for the multiphysics problem.(b) Aortic and mitral valves in the open and closed configuration.(c) Perfusion regions in the biventricular geometry.

Figure 4 :
Figure 4: Results from a physiological simulation.a) electromechanics of the left heart: calcium concentration during ventricular depolarization and displacement magnitude during ventricular contraction.b) left heart hemodynamics: on the top, volume rendering of velocity magnitude and pressure during filling and ejection phases; on the bottom, focus on the epicardial coronary arteries.c) myocardial perfusion: Darcy pressure in three different compartments during filling.Test I.

Figure 6 :
Figure 6: Results from physiological 3D Navier-Stokes -3D Darcy simulation in a representative heartbeat: (a) flow rates in pulmonary veins and aortic outlet section; (b) flow rates in epicardial coronary arteries; (c) average pressure in aortic outlet sections and coronary outlets; (d) pressure in the three Darcy's compartments.Test I.

Figure 7 :
Figure 7: Simulation of Aortic Regurgitation (AR): (a) aortic valve in ventricular diastole modeled with the RIIS method, comparison between physiological (PH) and AR cases; (b) volume rendering of velocity magnitude during ventricular filling; (c) aortic pressure prescribed on the ascending aorta outlet section in the PH and AR cases; (d) coronary flowrates over time, comparison between PH and AR cases; (e) velocity during diastole in the coronary tree, comparison between PH and AR cases; (f) Mean Blood Flow at diastolic peak in PH and AR cases.Test II and comparisons with PH case.

Table 1 :
Simulation Physics Mesh size [mm] Mesh details and time step sizes for the electromechanics (EM) and computational fluid dynamics (CFD) -Darcy simulations.Electrophysiology-ionic E − I, force generation A and mechanics problems M are solved on the same left heart mesh.C is the 0D circulation problem.F and D are the fluid dynamics and Darcy problems, respectively.