Simulation of the Perfusion of Contrast Agent Used in Cardiac Magnetic Resonance: A Step Toward Non-invasive Cardiac Perfusion Quantification

This work presents a new mathematical model to describe cardiac perfusion in the myocardium as acquired by cardiac magnetic resonance (CMR) perfusion exams. The combination of first pass (or contrast-enhanced CMR) and late enhancement CMR is a widely used non-invasive exam that can identify abnormal perfused regions of the heart via the use of a contrast agent (CA). The exam provides important information to the diagnosis, management, and prognosis of ischemia and infarct: perfusion on different regions, the status of microvascular structures, the presence of fibrosis, and the relative volume of extracellular space. This information is obtained by inferring the spatiotemporal dynamics of the contrast in the myocardial tissue from the acquired images. The evaluation of these physiological parameters plays an important role in the assessment of myocardial viability. However, the nature of cardiac physiology poses great challenges in the estimation of these parameters. Briefly, these are currently estimated qualitatively via visual inspection of images and comparison of relative brightness between different regions of the heart. Therefore, there is a great urge for techniques that can help to quantify cardiac perfusion. In this work, we propose a new mathematical model based on multidomain flow in porous media. The model is based on a system of partial differential equations. Darcy's law is used to obtain the pressure and velocity distribution. CA dynamics is described by reaction-diffusion-advection equations in the intravascular space and in the interstitial space. The interaction of fibrosis and the CA is also considered. The new model treats the domains as anisotropic media and imposes a closed loop of intravascular flow, which is necessary to reproduce the recirculation of the CA. The model parameters were adjusted to reproduce clinical data. In addition, the model was used to simulate different scenarios: normal perfusion; endocardial ischemia due to stenosis in a coronary artery in the epicardium; and myocardial infarct. Therefore, the computational model was able to correlate anatomical features, stenosis and the presence of fibrosis, with functional ones, cardiac perfusion. Altogether, the results suggest that the model can support the process of non-invasive cardiac perfusion quantification.

1. Introduction

Cardiovascular diseases have a close relation with high blood pressure, smoking, diabetes, lack of exercise, obesity, high blood cholesterol, poor diet, and excessive alcohol (Mendis et al., 2011; Mehta et al., 2015). Many of these conditions can affect myocardial perfusion (MP), which is the distribution process of oxygen and nutrients to the cardiac tissue. Among all the injuries related to a compromised cardiac perfusion that can evolve inside the heart, one that occurs frequently is the obstruction of a coronary artery, or stenosis. This leads to perfusion problems, the coronary artery diseases (CAD), also referred to as ischemic heart diseases (IHD) or simply as ischemia. Basically, the stenosis reduces the supply of oxygen and nutrients (perfusion defect) in the region of the myocardium that should be irrigated by the coronary artery in question. This condition may lead to more severe cases, like an infarct, which can provoke irreversible damages to the cardiac tissue, including death. There are several exams used to evaluate myocardial perfusion in clinical practice. Two that are widely used are computer tomography (CT) and the magnetic resonance imaging (MRI). Both consists in the administration in the patient of a contrast agent (CA) which, in the images generated by the chosen protocol, assume a specific contrast on poor perfused regions. The information that can be inferred with the acquired images has a fundamental role in the diagnosis, management, and prognosis of ischemia and infarct. In this work we present a new mathematical model to describe cardiac perfusion in the myocardium as acquired by cardiac magnetic resonance perfusion exams. Our results suggest that the model is an valuable tool that can support the process of non-invasive cardiac perfusion quantification. The next section gives more detailed explanations of MP and of contrast-enhanced MRI, which is the focus of this work.

1.1. Background

For the heart to have the necessary strength to pump blood throughout the body, it needs a constant supply of oxygen and nutrients. This process can be summarized as follows: the blood enters into the myocardium through the epicardial coronary arteries. As Figure 1 summarizes, the coronary arteries branches in smaller arteries, the arterioles. These in turn, branches into even smaller vessels. This process continues down to the capillaries, which are the smallest vessels in terms of diameter. It is on this level where the leakage of oxygen and nutrients from intra to extravascular media happens. After that, the capillaries increase in size and become veins. The venous blood follows the opposite direction and arrives at the venules followed by the coronary veins, these last vessels leave the epicardium to the lungs to be filtrated. Previous studies have shown that the blood-tissue exchange on the capillary level may be a complex phenomenon. For instance, in Bassingthwaighte et al. (1989) and Dash and Bassingthwaighte (2006), a multidomain model that accounts for the blood flow in capillary, endothelium, interstitium and parenchymal cells was developed, that includes the fluid exchange between each domain. However, the present work will not take into account such level of details.

FIGURE 1

Figure 1. Simplified scheme representing the relation between arteries and veins, in addition with the direction of blood flow. The coronary arteries get smaller and become arterioles, which get smaller again and become arterial capillaries. After the leak of oxygen and nutrients, the arterial capillaries now are called venous capillaries, which get bigger, becoming venules, and finally coronary veins.

In order to evaluate MP, several exams can be used. One that provides high quality images is computer tomography (CT). However, CT uses a harmful and allergenic substance as contrast agent (Brenner and Hall, 2007; Kroeker et al., 2011). Another option is contrast-enhanced MRI (Nikolaou et al., 2011). This one consists in the administration of a dose of Gadolinium with diethylenetriaminepentacetate (Gd-DTPA) in the patient that helps to identify patterns in the myocardial tissue. Contrast-enhanced MRI has different protocols. We can mention, for example, the first pass (FP) (Gerber et al., 2008) and the late-enhancement (LE) (Pop et al., 2013). As the name says, the FP uses the first passage of the CA through the circulatory system to generate the images, in which the regions that are well perfused by the CA acquire a characteristic glow, whereas dark regions indicate problems in perfusion. These problems can be related for instance with ischemia or infarct regions. The FP evaluates the myocardium for 50 s. On the other hand, LE enables the differentiation between a ischemic region, potentially still viable, and infarct. This exam evaluates the myocardium 10 min after the CA administration, and opposite to the FP, the infarct area is the one that gets brighter. Therefore, both FP and LE can be used as complementary MRI perfusion exams.

Nevertheless, the evaluation of the images acquired by the exams is a tedious and error-prone process performed by highly specialized physicians. To support and improve this process different tools that seek a more quantitative evaluation of myocardial perfusion are under development. In Chiribiri et al. (2013) and Alves et al. (2017), such process is studied using phantoms specifically designed to simulate contrast-enhanced MRI exams. In Aquaro et al. (2013), flow quantification is approached by calculating the maximum upslope of the myocardial signal. In Jerosch-Herold (2010) and Zarinabad et al. (2012), the quantification is given in terms of signal deconvolution techniques (Fermi deconvolution). In addition, several mathematical models have been developed with this goal. In Michler et al. (2013), a model is presented that captures pressure and velocity profiles of the blood in 3D representation of the heart. In Alves et al. (2016), the authors performed simulations considering regions with low permeability, showing the behavior of the CA in different scenarios. In Cookson et al. (2014), for the purpose of reproducing the exchange of CA between intra and extravascular spaces, the authors take into account a 2-domain formulation. Although clearly substantial progress has been made, there is still room for further improvements.

Another invasive exam that indirectly evaluates perfusion is the Fractional Flow Reserve (FFR) (Lee et al., 2017), which gives information related to the severity of a coronary stenosis. It is important to mention that there are softwares that can estimate FFR by combining mathematical models (also PDEs) and images acquired by CT exams. These softwares are available for clinical use (Hlatky et al., 2015).

In this work, we propose a similar tool that also uses mathematical models based on PDEs and images, in this case from contrast-enhanced MRI exams. The model simulates the blood flow in a short-axis slice of the myocardium (cardiac muscle). Simulations of the model provides the dynamics of CA in the initial stage, as captured by first pass, and after many minutes of initial CA injection, as captured by late enhancement. In order to describe myocardial blood flow, Darcy approach for fluid in porous media is used, coupled to advection-diffusion-reaction equations to evaluate the CA dynamics. Clinical data were used to adjust the parameters of the model. Some specific contributions of this work are: (1) the simulation of CA dynamics in a domain that represents the short axis of the myocardium; (2) the development and test of a mathematical formulation that takes into account three domains (intravascular, extravascular, and fibrosis); (3) a model that contains a closed loop of intravascular flow, which is necessary to reproduce the recirculation of the CA; (4) the model parameters were adjusted to reproduce clinical data; and (5) analysis of the CA dynamics and its quantification for scenarios of normal perfusion, endocardial ischemia, and myocardial infarct.

In section 2 the multi-domain mathematical model is presented. Section 3 describes the numerical methods. Section 4 presents the results obtained with the simulations. Section 5 presents a general discussion of the results and the limitations of the model. Finally, section 5.3 presents the conclusions.

2. Multi-Domain Model for the Blood Circulation and the Contrast Dynamics

In this work, the blood flow in the myocardium is considered as a single-phase flow in porous media. Basically, a porous media is a solid filled with empty spaces, which are connecting each other. The fraction of empty spaces in the total sample volume, which is called porosity, is given by

$\begin{array}{cc}\varphi =\frac{V_p}{V_t},& (1)\end{array}$

where V_p is the volume of the pores and V_t is the total volume of the sample. This kind of analysis is widely used in oil extraction (Taber, 1969) and aquifer contamination by organic compounds (Abriola and Pinder, 1985).

In order to better understand the complex phenomenon of myocardial perfusion and, in particular, how MRI images are related to ischemia or infarct, previous works have considered the myocardial arterial tree as a porous media (Bassingthwaighte et al., 1989; Michler et al., 2013; Cookson et al., 2014; Alves et al., 2015, 2016). This approach represents a major simplification. As presented in Figure 1 the vessel network is discrete and involves non-trivial topological features. Therefore, a lot of information is lost when we represent it by a continuous porous media. Nevertheless, this simplified model is still able to provide useful insights about the phenomena of myocardial perfusion. For instance, in Michler et al. (2013), Darcy's law for flow in porous media was used to model the pressure field and blood velocity in the intravascular tree that perfuses the myocardium. In addition, the same work proposes the use of three different compartments, each one dealing with a diameter range of the vascular tree. In addition to this study, in Alves et al. (2015, 2016), a second equation was proposed, this was a diffusion-advection PDE used to describe the CA dynamics in scenarios of both healthy and diseased tissues. For this purpose, a small region in the domain was set with a lower permeability. In Cookson et al. (2014), a third equation was added to capture the CA dynamics in two different domains, intravascular and extravascular. In this work, we will present a new mathematical model that is able to reproduce three important scenarios: perfusion on a healthy tissue, on an ischemic region, and on an infarct. The last two can happen due to a low pressure difference between epicardial and endocardial surfaces, which can be caused by a stenosis in a coronary artery. Compared to the aforementioned previous work, the new model is more realistic:

(I) It reproduces the Gd-DTPA recirculation: the CA substance is partially washed out after the first pass in the myocardium. However, as can be seen in Lima et al. (1995) and Daly and Kwong (2013), the blood with CA returns to perfuse the heart in cycles. The new arrival of this CA wave implies in a small increase of the SI (Signal Intensity), which is captured by the new model here presented;

(II) Here, our domain is not a simple rectangle as those used before in Cookson et al. (2014) and Alves et al. (2016). In this work, we use a realist short-axis cut of the cardiac left ventricle. Therefore, the simulation results are closer to the images actually acquired by an MRI exam;

(III) The CA is confined to the interstitial (extracellular) space, since it is known that it does not to enter into the cardiac cells (intracellular space);

(IV) When the scenario of infarct is simulated, we added a sortion term in the equations to capture how CA is trapped inside the fibrosis network present in the scars (injured portion of the interstitium).

Cardiac perfusion happens together with cardiac contraction, see Figure 2 for an example of the different phases of cardiac contraction as captured by cine MRI exam. However, for contrast-enhanced MRI physicians use a single phase, usually end-diastole, and analyze CA dynamics comparing these images acquired on different heart beats. Therefore, the sequence of contrast-enhanced MR images is analyzed on a non-contracting heart. It is a simplification, but serves well the purpose as most of the perfusion takes place during diastole. In this work, we will use the same hypothesis and consider CA dynamics on a non-contracting heart. Therefore, the variables of interest are myocardial pressure and CA concentration. The blood fluid is taken as incompressible, homogeneous, and Newtonian.

FIGURE 2

Figure 2. Six MR images at different phases of cardiac cycle (it is possible to observe the movement of the myocardium). Extracted from Koch et al. (2011).

2.1. Mathematical Model

Figure 3 shows a framework of a multidomain model for cardiac perfusion of the left ventricle. Each domain is modeled by a porous media where the variables of interest are pressure, flow, and CA concentration. On the left side of the figure we have the domains that represent the intravascular flow in the heart: arterial domain (top), capillary domain (middle), and venous domain (bottom). On the right side of the figure we have the extravascular domain in the heart which is subdivided in interstitial space and fibrotic network (shaded area associate to an infarction). This multidomain model captures the flux between the domains as part of a closed circuit. There is an influx through the coronary arteries on the epicadial surface of the arterial domain. The blood follows to the capillary domain, where there is the exchange of plasma and CA with the extracellular domain. The flow continues from the capillary domain to the venous one, from where the coronary veins leave the heart. This outflux enters a recirculation circuit, which is a simplification of the pulmonary and systemic circulation processes.

FIGURE 3

Figure 3. Multidomain model for cardiac perfusion of the left ventricle. Each domain is modeled by a porous media. On the left side of the figure we have the domains that represent the intravascular flow in the heart: arterial domain (top), capillary domain (middle), and venous domain (bottom). On the right side of the figure we have the extravascular domain in the heart which is subdivided in interstitial space and fibrotic network (shaded area associate to an infarction). A recirculation circuit, which is a simplification of the pulmonary and systemic circulation processes, enforces a closed circuit.

This multidomain model can be further extended or simplified, depending on the main goals of a particular study. For instance, it could be used to reproduce the transmural perfusion and respective pressure gradients from subepicardium to subendocardium in the arterial domain as well as the outward perfusion, i.e., from subendocariaum to subepicardium, in the case of the venous domain. In addition, by setting an appropriate anisotropic permeability tensor of the vascular domain, it may be used to reproduce the preferential perfusion direction of the micro-vascular system, which follows the orientations of cardiac fibers. In this work, we have used a simplified version of the described multidomain approach by combining the arterial and capillary domains into one intravascular domain and by neglecting the venous domain. This simplification is presented in Figure 4A. As our results will show, this simplified model is quite successful in reproducing several aspects of CA dynamics in myocardial perfusion. In particular, it allowed the study of the influence that a coronary obstruction has in the pressure distribution and cardiac perfusion, including scenarios of ischemia and infarction.

FIGURE 4

Figure 4. The simplified two-domain model (A), used to simulate healthy and acute ischemia cases; and the three-domain model extended by a fibrosis domain (B), used to simulate an infarct region.

To obtain a transmural pressure gradient, boundary conditions are prescribed, which generates a pressure gradient and therefore a blood flow from epi- to endo-cardium. This part of the model can be summarized as follows:

$\begin{array}{cc}\nabla ·\overrightarrow{v}=0\text{in}{\Omega }_i,& (2)\end{array}$

$\begin{array}{cc}\overrightarrow{v}=-\mathbf{K}\nabla p\text{in}{\Omega }_i,& (3)\end{array}$

where K is the permeability tensor; p is the pressure; $\overrightarrow{v}$ is the Darcy velocity and Ω_i is the intravascular domain. The boundary conditions for Equation 2, written in terms of pressure, are given by Dirichlet's boundary condition:

$\begin{array}{cc}p=p_o\text{in}{\text{Γ}}_o,& (4)\end{array}$

$\begin{array}{cc}p=p_i\text{in}{\text{Γ}}_i,& (5)\end{array}$

where Γ_o and Γ_i are epicardial and endocardial boundaries, respectively. For the CA dynamics we considered advective ($\overrightarrow{v}C$) and diffusive (−D∇C) components for both the intravascular and extravascular domains. Thus, this system is given by

$\begin{array}{cc}\frac{\partial (\varphi C_i)}{\partial t}+\nabla ·\overrightarrow{v}{C}_i-\varphi \nabla ·(D_i\nabla C_i)+f=0\text{in}{\Omega }_i,& (6)\end{array}$

$\begin{array}{l}\frac{\partial ((1-\varphi )\lambda C_e)}{\partial t}-(1-\varphi )\lambda \nabla ·(D_e\nabla C_e)-f+(1-\varphi )\lambda k_e{C}_e+g\\ =0\text{in}{\Omega }_e,& (7)\end{array}$

where f models the transfer between the two domains and is given by

$\begin{array}{cc}f\{\begin{array}{l}P\times (C_i-C_e),if{C}_i>C_e,\\ 0,\text{otherwise},\end{array}& (8)\end{array}$

where P is the endothelial permeability. This parameter is known to vary for different pathologies (Bassingthwaighte et al., 1989). Section 5 presents the sensitivity of it. ϕ is porosity, i.e., the fraction of intravascular domain, and thus 1−ϕ is the fraction of the extravascular domain. C_i is the concentration of CA in the intravascular domain (Ω_i), C_e is the concentration of CA in the interstitial space (Ω_e). D_i and D_e are diffusion coefficients. −k_eC_e models the flux from the interstitial space to the venous system. As mentioned before, the extravascular space can be divided in interstitial (extracellular) and intracellular portions. We assume that the CA does not enter in the myocytes. Thus, λ represents the fraction of extravascular (1−ϕ) space that is due to the interstitium.

As mentioned before, in the case of chronic infarct, scars may develop by the increase of fibrosis to replace dead myocytes. Therefore, when performing the two protocols of contrast-enhanced MRI perfusion, two different situations may happen: in first pass, because it takes around 50 s, there is not enough time for the CA to reach the fibrotic area. Thus, on this exam, the images provided by the MRI will indicate a dark color in the injured area, whereas the remain will be bright. On the other hand, in late enhancement, there is enough time for the CA to reach the area. In addition, the wash out of CA is delayed, since the network of fibrosis will behave as a trap for CA. The described phenomenon, of fluid, particles or substances been attached to a “solid phase,” in physics, is called sortion. Therefore, the dynamics of the CA in the fibrosis is given by

$\begin{array}{cc}\frac{\partial ((1-\varphi )\lambda {\lambda }_f{C}_f)}{\partial t}+(1-\varphi )\lambda {\lambda }_f{k}_f{C}_f-g=0,\text{in}{\Omega }_f,& (9)\end{array}$

where C_f is the concentration of CA in the fibrosis domain (Ω_f) and k_fC_f models the flux from the fibrotic network to the venous system. g reflects the flow of CA from interstitium to the fibrotic network. It is nonzero only when we simulate the infarct scenario. It is given by

$\begin{array}{cc}g=(1-\varphi )\lambda {\lambda }_f{k}_{ef}{C}_e,& (10)\end{array}$

where k_ef is the rate in which the contrast goes from the interstitium to the fibrosis. λ_f is the fraction of the interstitium with fibrosis. See Figure 4B.

For the outer boundary condition (Γ_o), a prescribed flux ($\overrightarrow{v}Q(t)$) is used, applying a transient Gaussian formula:

$\begin{array}{cc}Q(t)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{1}{2}\left(\frac{t-t_{peak}}{\sigma }\right)}+X(t,x,y),& (11)\end{array}$

where σ² reflects the variance of the contrast agent infusion, t_peak is the mean of the Gaussian which sets the time value for the peak of the gadolinium intensity. X(t, x, y) is a term that must be added to simulate the CA cyclic behavior observed in experiments for the first pass and late enhancement MRI perfusion (Lima et al., 1995; Daly and Kwong, 2013). In the next section we will present more details about this term.

Homogeneous Neumann boundary conditions are used in Γ_i for both equations:

$\begin{array}{cc}{D}_i\nabla C_i·\overrightarrow{n}=0\text{in}{\text{Γ}}_i,& (12)\end{array}$

$\begin{array}{cc}{D}_e\nabla C_e·\overrightarrow{n}=0\text{in}{\text{Γ}}_i,& (13)\end{array}$

where $\overrightarrow{\mathbf{\text{n}}}$ is the normal vector. The same is assumed for Γ_o in Equation 7:

$\begin{array}{cc}{D}_e\nabla C_e·\overrightarrow{n}=0\text{in}{\text{Γ}}_o.& (14)\end{array}$

As initial conditions, we have:

$\begin{array}{cc}{C}_i(x,0)=C_{i_0}(x)\text{in}\Omega ,& (15)\end{array}$

$\begin{array}{cc}{C}_e(x,0)=C_{e_0}(x)\text{in}\Omega .& (16)\end{array}$

2.2. 1D Model to Simulate the CA Cyclic Behavior

Like Equation 6, that simulates the contrast agent dynamics in the myocardium, a similar model, however simplified to 1D, was used to capture the recirculation of contrast. Figure 4A uses a line to represent the recirculation. This line connects the outflux and influx via the variables C_L(t) and C_R(t), respectively.

Basically, the values C_L(t) and C_R(t) are the amount of contrast at the beginning and at the end of the recirculation, respectively. The amount C_L(t) is given by

$\begin{array}{cc}{C}_L(t)=\frac{{\displaystyle {\int }_{{\Omega }_i}{C}_i}d{\Omega }_i}{\left|{\Omega }_i\right|}.& (17)\end{array}$

The value given by Equation 17 represents the average of the contrast in the intravascular medium. It is considered as a Dirichlet condition in the left boundary of the 1D domain. Thus, for a given velocity, the average of contrast travels in the domain through the advection-diffusion-reaction equation

$\begin{array}{cc}\frac{\partial C_{out}}{\partial t}+\nabla ·{\overrightarrow{v}}_{out}{C}_{out}-\nabla ·(D_{out}\nabla C_{out})+k{C}_{out}=0\text{in}{\Omega }_{out}.& (18)\end{array}$

The amount obtained in the right boundary of the 1D domain, C_R(t), is taken as the recirculation term X(t, x, y) in Equation 11. Physiologically, a portion of the contrast is eliminated in the recirculation, specially by the kidneys. This is modeled by the reaction term kC_out. This phenomenon explains the decrease of CA in the subsequent passes in the myocardium. Although the recirculation is more complex and involves the dispersion of the CA in the whole body in different quantities, this simple 1D model captures the CA elimination and its cyclic behavior in the MRI perfusion exam, as presented in section 4.

3. Numerical Methods and Experiments

For the discretization, the FVM (Finite Volume Method) formulation was used, which consists in the evaluation of influx and outflux of a control volume around each node of the mesh. The FVM applied to Equation 2 (an elliptic PDE) gives us a linear system with one equation for each node. The control volumes were taken as squares of size equal to h. This system was solved by the iterative method of Jacobi (Quarteroni et al., 2010), with an error tolerance (stop criteria) set to 10⁻⁸. Once the pressure field is known, Equation 3 is used to obtain the velocity in the faces of the control volumes. The velocity is used in the second term of the left side in Equation 6 and the discretization of the diffusive term follows that of the pressure equation. For the diffusive term, the central differences scheme works very well, without any major problems. But for the advective term a more sophisticated scheme must be used to ensure accuracy and convergence of the numerical method. For the numerical approximation of this term, the third-order polynomial upwind scheme (TOPUS) was implemented. For more details about this upwind scheme, see Ferreira et al. (2012) and Alves et al. (2016). We used the explicit Euler method in the temporal part of the advection-diffusion-reaction equations. All the mathematical calculations in each step of the discretization can be found in Alves et al. (2016).

3.1. Numerical Simulations

We simulated the CA dynamics on a short-axis slice of a human ventricle taken from an open repository (Koch et al., 2011) (see Figure 5). During the segmentation process, we have included the papillary muscles as part of the myocardium (Gommans et al., 2016). However, we note that the results presented in this work are not affect by this particular short-axis or segmentation. Similar results can be obtained with different short-axis and different segmentation schemes that do not include the papillary muscles.

FIGURE 5

Figure 5. Performed simulations: (A) healthy, (B) ischemic, and (C) infarct. (D) Shows the regions of interest for the evaluation of CA signal intensity. When simulating ischemia and infarct, we sum up the concentration of CA in these regions of interest, for each time instant. The colored circle that surrounds the epicardium represents the pressure boundary condition for each scenario. Note the different minimum values of the two color bars.

We discretized the segmented image using a total of 10,000 finite volumes, i.e., the resolution used for the simulations was Δx = Δy = 0.5mm. The MRI image used had a resolution of 1.5 × 1.5 mm² (Koch et al., 2011). Therefore, the resolution used for the simulations was three times higher than the one from the MRI image. However, as we have not processed the segmented image with any smoothing technique, our simulation domain presented in Figure 5 is also pixel-based, with non-smooth boundaries.

In order to assess the model in different scenarios, three main simulations were performed: the first is a healthy situation, with normal values of pressure in epi- and endo-cardium. The second simulation tries to reproduce an ischemia. For this purpose, a drop of pressure is considered in the epicardium, which emulates a coronary stenosis. For the third simulation, we consider a more severe stenosis and a region with fibrosis to reproduce an infarct. The values for the prescribed pressures in the non-ischemic case are taken from Chapelle et al. (2010). For the diseased simulations, the values for the drop of pressure in a location of the epicardium (that simulates a coronary artery) were taken from works that study the FFR technique (Pijls et al., 2000; Sant'Anna et al., 2007; Lee et al., 2017). Here we present the features of each of these scenarios (see Figure 5):

• 1st (normal): p(x, y) = 0.0kPa, ∀(x, y) ∈ Γ_i, and p(x, y) = 2.0kPa, ∀(x, y) ∈ Γ_o, see Figure 5A;

• 2nd (ischemia): In a small region of the right-hand side of Γ_o pressure is reduced (to simulate a 30% stenosis), p(x, y) = 1.4kPa. The exact location is presented by the arrow in Figure 5B. In the rest of Γ_o, pressure is set to its normal value, p(x, y) = 2.0kPa. See the circular color bar around the epicardial surfaces that presents pressure values at the epicardial surface in Figure 5B. At the endocardial surface p(x, y) = 0.0kPa, ∀(x, y) ∈ Γ_i;

• 3rd (infarct): In a small region of the right-hand side of Γ_o pressure is reduced (to simulate a 80% stenosis), p(x, y) = 0.4kPa. The exact location is presented by the arrow in Figure 5C. In the rest of Γ_o, pressure is set to its normal value, p(x, y) = 2.0kPa. See the circular color bar around the epicardial surfaces that presents pressure values at the epicardial surface in Figure 5C. At the endocardial surface p(x, y) = 0.0kPa, ∀(x, y) ∈ Γ_i. In addition, this scenario also considers the presence of a fibrotic region near the endocardial surface and opposed to the simulated stenosis. This fibrotic region is called diseased, and its location and size are presented in Figure 5D.

The 3rd scenario has more particularities than just the pressure gradient. As mentioned before, the CA is confined to the interstitial space, and the measurements of this volume have several particularities (Al-Wakeel-Marquard et al., 2017). For the purpose of this work, we use λ = 0.25 (Arai, 2011) for the healthy and ischemic scenarios. In the case of an infarct region, many cells are dead and thus the fraction related to the interstitium is big, and we set λ = 1.0. In the same way, λ_f, the fraction of fibrosis, is zero for the healthy and ischemic cases, whereas for infarct λ_f = 0.50. Finally, for the infarct scenario the wash out of CA is delayed, since it is trapped in the fibrotic network. Therefore, the rate of CA from fibrosis to the venous system is k_f = 0.0007 s⁻¹, much slower than for the cases without fibrosis, where k_e = 0.003 s⁻¹.

The total porosity ϕ was assumed to be constant and equal to 0.10 (Cookson et al., 2014). The Gaussian formula (Equation 11) was applied only at the epicardium, Γ_o, with σ = 7.0 and t_peak = 22 s. The simulation time was set to 600 s, which is the average time used for the LE exam (Knowles et al., 2008; Arai, 2011). It is known that the exchange of oxygen, nutrients and CA from capillaries to the tissue is a phenomenon that involves different components of the biological system, such as the capillary, endothelium, interstitium, and parenchymal cells) (Bassingthwaighte et al., 1989). However, here we fix the endothelial permeability P of Equation 8 to P = 0.03 s⁻¹, in the three scenarios. Nevertheless, section 4 presents a sensitivity analysis of this parameter. With respect to the 1D model, the parameter values were: $\overrightarrow{v}=1.1$ mm s⁻¹ (velocity), $D_{out}=0.5\times 10^{-1}$ mm² s⁻¹ (diffusion coefficient), k = 0.015 s⁻¹ (elimination of contrast in the recirculation). Table 1 presents a summary of the parameter values for each scenario.

TABLE 1

Table 1. Parameters used in each scenario.

3.2. Anisotropy and Heterogeneity of the Intravascular Domain

The subendocardial permeability is, in average, two times higher than the subepicardial one (Smith et al., 2014). In addition, the preferential direction of permeability, in the short-axis, is parallel to the epicardial and endocardial surfaces (Smith et al., 2014). Therefore, in order to take these features into account, a permeability tensor K(x, y) was generated ∀(x, y) ∈ Ω_i. First we generate a transmural gradient of permeability by solving the Laplacian equation, Δw = 0 ∈ Γ_i, with w = 0 ∈ Γ_o and w = 1 ∈ Γ_i, i.e., using Dirichlet boundary condition at the epicardial and endocardial boundaries, respectively. This gives us values between 0 and 1 for each node in the simulated domain. Next, we compute the pair of vectors {e_l, e_t}, where e_l is the first and e_t the second preferential direction of permeability by making e_t = ∇w, and ensuring that e_l is perpendicular to e_t. With these two vectors calculated for each point (x, y) we can write the permeability tensor K(x, y) as follows:

$\begin{array}{cc}K(x,y)=K_l(x,y)I+(K_t(x,y)-K_l(x,y))e_l(x,y)e_l{\left(x,y\right)}^t,& (19)\end{array}$

where K_l is the permeability value along e_l and K_t is the permeability along e_t (transversal to e_l). These values can be calculated using the weight w:

$\begin{array}{cc}Kl(x,y)=K1(1-w(x,y))+2.0K1w(x,y),& (20)\end{array}$

$\begin{array}{cc}Kt(x,y)=K2(1-w(x,y))+2.0K2w(x,y),& (21)\end{array}$

where K₁ and K₂ are the permeabilities in outer boundary, with K₁ approximately 10.8 times higher than K₂ (Smith et al., 2014); and the multiplication by 2.0 in the equation reflects that the subendocardial permeability is, in average, two times higher than the subepicardial one (Smith et al., 2014).

Table 1 presents a summary of the parameters used for each considered scenario in this work.

3.3. Computing

The simulator is an in-house code implemented in C. Simulations were executed on an Intel(R) Core(TM) i7 3 GHz 8 Gb. Each simulation of the late enhancement protocol (10 min of perfusion) takes approximately 708 s. This first version of the code was developed as a sequential one, i.e., parallel programming was not used, but it is planned for future works.

4. Results

Section 4.1 presents the pressure profile obtained from the Darcy approach (Equation 2). In the sequence, we present qualitative results of the CA perfusion on a short axis. Both, first pass and late enhancement are considered for the three different scenarios (healthy, ischemic, and infarction). In section 4.2 we present quantitative results of the model, and compare the signal intensities of the CA in the first pass and in the late enhancement exams to clinical data from the literature (Knowles et al., 2008; Daly and Kwong, 2013).

4.1. Spatio-Temporal Dynamics of CA Perfusion in the Heart

Figure 6 presents the numerical solutions obtained for the pressure distribution for the three considered scenarios: (a) healthy situation; (b) for the case where a stenosis of 30% is simulated; and (c) for a stenosis of 80%. In the real phenomenon, the pressure profile changes during the cardiac cycle. As mentioned before, we use the same hypothesis as the clinicians and evaluate the dynamics of perfusion at the same instant of each cardiac cycle, i.e., heart contraction is not taken into account. The calculated pressure and velocity will be used in the transport equation (Equation 6). From the figure we observe that the stenosis decreases the trasmural gradient of pressure.

FIGURE 6

Figure 6. Pressure field of the (A) healthy (B) ischemic and (C) infarction scenario (in kPa). It is possible to observe that the pressure is higher near the epicardium area and lower near endocardium, as expected. The Dirichlet boundary condition mentioned in section 3.1 generates a pressure gradient. Using the Darcy's law we obtain a velocity field that is used in the advective term of Equation 6.

Figure 7 shows the final result, in terms of qualitative images of the myocardial short axis for the three simulated scenarios. For each scenario, we shown the concentration in the end of both first pass and late enhancement exams. The simulation domain was adapted from a real short axis MRI image (Koch et al., 2011). The model reproduces very well the qualitative behavior, for both first pass and late enhancement exams. During first pass, in a healthy heart the CA is homogeneously distributed, whereas in the cases of ischemia and infarction the affected region is revealed by low localized contrast values. At the time of late-enhancement, the CA attaches to fibrosis in the infarct region. This is reflected by the high localized contrast values at late-enhancement.

FIGURE 7

Figure 7. Qualitative images in the end of first pass (50 s) and late enhancement (600 s) exams. In the first scenario (healthy), the CA is homogeneously distributed, while in the second (ischemia), the diseased region is revealed by low localized contrast values. At the time of late-enhancement, the CA gets stuck in the infarction region as reflected by the high localized contrast values. All the simulations reproduced well the characteristics seen by the clinicians in the considered scenarios.

4.2. Signal Intensity of Contrast Agent

The variable of interest in this section is the signal intensity of the contrast agent in specific regions of interest (ROI). To compare the CA dynamics of the healthy and pathological cases, two regions of interest in the myocardium were marked as presented in Figure 5D: one captures the injured region and the other is on a remote healthy region. Each value of the time series correspond to the sum of the CA concentration in one specific ROI at a specific time.

Since we are setting a reduced pressure on the right-hand side of the epicardial surface to simulate a stenosis, we placed the diseased ROI as a subendocardial region opposite to the stenosis point. The remote ROI was placed on the subendocardial left-hand side of the short axis (see Figure 5D). These are common choices of ROIs to compare a diseased versus a healthy region (Hsu et al., 2006; Schuster et al., 2015). In our simulations, the remote ROI has a total area of 70 mm², while the diseased ROI has 55 mm². These sizes and locations of ROIs are commonly reported in the literature (Hsu et al., 2012; Pop et al., 2013; Schuster et al., 2015). For the third scenario that simulates an infarct, a fibrotic region was modeled at the same location of the diseased ROI.

Section 4.2.1 presents the results for the ischemic scenario. Section 4.2.2 presents the results of the infarction. Is important to remark that the fibrosis mentioned before was placed only in the diseased region, with a particular volume fraction, as discussed in section 2.1.

4.2.1. First Pass

Figure 8A shows the results obtained with respect to the intensity of CA signal during a first pass MRI perfusion exam in the two considered regions of interest. These results are in agreement with those found in the literature (Nagel et al., 2003; Jenson et al., 2013; Wang et al., 2015; Bakir et al., 2016). In particular, the experimental curves presented in Figure 8A are taken from Daly and Kwong (2013), where the first pass of healthy region is compared to a ischemic one. It is possible to observe the second pass of the CA, i.e., the remaining amount of it that was not eliminated by the kidneys and recirculates. Figure 8B presents how each domain, intra and extravascular, contributes to the total SI of CA over time for the case of ischemia. These results are also in agreement with the literature (Bogaert et al., 2005).

FIGURE 8

Figure 8. (A) Comparison between SI of the CA during time observed in Daly and Kwong (2013) and the numerical results of this work in some ischemic area (red) and remote and healthy area (blue). The oscillatory decay is due to the Gd-DTPA recirculation. (B) The contribution of intra- and extravascular SI of the CA in the ischemic area: (purple) total, (red) intravascular, and (black) extravascular concentrations vs. time.

4.2.2. Late Enhancement

After the first pass detects a region with poor perfusion, late-enhancement is performed to study myocardial viability, i.e., to differentiate between ischemia, but yet viable tissue, or infarct. This is possible because fibrosis is present in chronic infarct and the CA is trapped in the fibrotic network, taking longer to be washed out of the myocardium. If the perfusion problem is not so severe, like under ischemia, the CA is washed out at a similar rate as a remote/healthy region. The results shown in Figure 9 are in agreement with the late enhancement dynamics (Arai, 2011) and reproduces well the clinical data taken from Knowles et al. (2008). Figure 9 presents the contribution of both intra- and extravascular spaces in the CA dynamics in the infarct region. The intravascular curve shows that the CA in this domain is quickly washed out if compared to the CA in the extravascular space.

FIGURE 9

Figure 9. (A) SI of the contrast agent in scenarios remote/healthy (red) and infarction (black). The yellow curve presents the clinical data taken from Knowles et al. (2008). The CA remains in the infarct area because it attaches to fibrosis, and thus takes longer to be washed out. (B) Contribution of both intra and extravascular spaces in the CA dynamics in the infarct region. The intravascular curve (purple) shows that the CA in this domain is quickly washed out, while the CA remains in the extravascular space (red curve) after 2 min.

5. Discussion

As we present a new model for cardiac perfusion with focus on the application of contrast-enhanced MRI, this section begins by presenting sensitivity tests of the proposed models. This will contribute to the discussion regarding the importance of our main hypothesis, the positive outcomes of our model as well as its limitations.

5.1. Sensitivity of the Parameters

This section is divided in three parts. Section 5.1.1 shows the importance of the recirculation coupling in the model (Equation 18). In section 5.1.2, an effort to reproduce clinical data of LE for the infarction scenario without taking into account the domain of fibrosis is presented. In section 5.1.3, we demonstrate that the coupling of the model with a fibrotic domain and the inclusion of a sortion term provide a better alternative to model the experimental and clinical data.

5.1.1. Recirculation

Figure 10 is presented to reinforce the idea that the recirculation coupling is important to obtain results close to those obtained in clinical exams. The red curve is the clinical information from Daly and Kwong (2013) (same as in Figure 8A), and due to the lack of data, it stops after 50 s. The figure also shows the dynamics of the CA with the recirculation model (black, same as in Figure 8A) and without it (yellow) for 200 s of perfusion. As presented, if one does not take into account the recirculation model and its parameters (${\overrightarrow{v}}_{out}$, D_out and k from Equation 18), the oscillatory decay observed in the clinical and experimental data is not captured.

FIGURE 10

Figure 10. Simulations regarding ischemic scenarios: the first (black), shows the model coupled with the recirculation part (Equation 18) and the second (yellow), without it. The red curve shows the experimental data from Daly and Kwong (2013).

The whole process of recirculation is quite complex: the kidneys, for instance, play an important role in the physiological process of elimination of CA; the constant flux velocity assumed by the model is in fact a simplification, since the flux varies and depends on the location if the body, the structure of the vessels, viscosity, and many other features. Therefore, other variables and parameters could be added in this part of the model. However, it is not in the scope of this work to go deep into this topic. The solution presented here is taking into account the blood cycle in the whole human body as a single-phase flow in a simple 1D domain, which provides good qualitative agreement with experimental data. For future works, it would be interesting to study this phenomenon and look for other solutions.

5.1.2. Modeling Perfusion Without a Mathematical Description of Fibrosis

As mentioned before, without taking into account (Equation 9), that describes fibrosis, the model captures well the scenario of ischemia and is able to reproduce clinical data (Daly and Kwong, 2013). However, here we show that without the mathematical description of fibrosis the model can not reproduce the dynamics of infaction, mainly those found in the time-scale of the late-enhancement exam.

Figure 11 presents the sensitivity analysis of the model without the coupling to Equation (9). In other words, it shows an attempt of reproduce the clinical data of infarction from Knowles et al. (2008) using just a part of the proposed mathematical model. All curves in Figure 11 refer to the dynamics of the CA in the diseased region according to Figure 5D. The parameters were only varied is this ROI, and kept constant in the other regions. In addition, when varying one parameter, all the others are fixed according to the infarction column in Table 1.

FIGURE 11

Figure 11. Sensitivity analysis of the parameters using the model two domains: (A) λ—the fraction of extravascular media that corresponds to the interstitium; (B) P—the rate of CA that goes from intravascular to extravascular space; (C) Pressure difference between epicardium and endocardium in the infarcted region; (D) k_e—the rate of CA going from extravascular media to the venous domain (a term of decay). When varying one parameter, all the others are fixed according to the infarction column in Table 1.

Figure 11A presents the sensitivity of λ (fraction of interstitium in the extravascular domain). For high values of λ there is more space for the CA. Thus, more CA is hold in the interstitial domain: time-to-peak, peak value and value at 600 s increase. These are all in the correct direction to reproduce the data in Figure 11D (black and yellow). However, even setting λ = 1, i.e., to its maximum value, the values of Figure 11A are far from those of Figure 11D (black and yellow). Figure 11B shows the sensitivity of P, which models the rate of CA exchange between intra and extravascular domains. For instance, P = 0 s⁻¹ indicates that there is no contrast agent flowing from one domain to the other, which means that it is confined to the intravascular domain. Hence, in this particular case, the CA is quickly washed out. Increasing P, the CA goes faster to the extravascular domain, hence it takes longer to be washed out to the venous system, but it is still completely washed out at 600 s. Figure 11C shows how the pressure gradient, which is related to the severity of coronary stenosis, influences the signal intensity of the CA in the diseased region. For example, in the case of an ischemia (this one is related with Figure 5B), the obstruction of 30% leads to a pressure gradient of up to 1.45 kPa in the right portion of the myocardium. In the worst scenario shown here, an obstruction of 80% leads to a pressure gradient of up to 0.45 kPa between epi- and endo-cardium (related with Figure 5C). As expected, the pressure gradient is directly proportional to the rate with which the CA reaches the diseased area. In Figure 11D, two different values for k_e (transfer of CA from the extravascular media to the venous system) are used. As expected, decreasing the value by two order of magnitude, the CA takes longer to be washed out. However, once gain, the curves obtained are far from the goal 11D (black and yellow).

In summary, using the formulation with only two domains (intra- and extra-vascular), it was not possible to reproduce the observed scenario of CA dynamics for a case of infarction. However, in Cookson et al. (2014), a formulation with two domains was used and some results were similar to those observed for infarction. To achieve this, however, the authors had to change the permeability of the media (K in Equation 3) in specific regions. Their purpose was to evaluate the CA dynamics in different scenarios, without specifying any pathology. For the cases of healthy, ischemic or infarct myocardium, the three scenarios treated here, the changes used in the permeability K would be difficult to justify. Thus, for the simulation of an infarction it seems to be important to model fibrosis, as we propose in our new formulation.

5.1.3. Modeling Perfusion With a Mathematical Description of Fibrosis

Figure 12 shows the dynamics of the CA in the infarct region taking into account the third domain and Equation 9. Figure 12A shows the dynamics of the CA for different values of λ (interstitial fraction of the extravascular media). For the healthy regions, its values is around 0.25 (Knowles et al., 2008; Arai, 2011). Once again, for all the tests, when varying one parameter, all the others are fixed according to the infarction column in Table 1. As in Figure 11A, for high values of λ there is more space for the CA. Thus, more CA is hold in the interstitial domain. However, using this model, time-to-peak, peak value, and value at 600 s are all much more sensitive to λ. And therefore, few adjustments of the parameter are enough to fit different experimental values. Figure 12B presents the sensitivity of k_f, the rate in which the CA goes from fibrosis to the venous system. The higher its value, the faster the CA is washed out of the myocardium. If k_f = 0.0 s⁻¹, no CA is leaving the fibrotic region, hence a plateau appears. The best value that reproduces the clinical data is $k_f=7\times 10^{-4}$ s⁻¹. Figure 12C presents the sensitivity of P. The value that best fitted the clinical data of infarction is P = 0.025 s⁻¹. Once again, if P = 0.0 s⁻¹, there is no communication between intra and extravascular domains, hence the CA is quickly washed out. Figure 12D shows results for different values of k_e (decay rate of CA). It reaches saturation at 7 × 10⁻⁵ s⁻¹. In summary, the presented formulation that models a third domain for fibrosis has proved to be a good alternative for the simulation of the dynamics of CA in the case of infarction. In addition, the curves generated with this model were more sensitive to many of its parameters. Therefore, a few adjustments of some parameters are enough to fit different experimental values.

FIGURE 12

Figure 12. Sensitivity analysis of the parameters in the infarcted scenario using the model with three domains: (A) λ—the fraction of interstitium in the infarcted area; (B) k_f—the exchange rate of CA going from interstitium to the venous system; (C) λ_f—the fraction of interstitium that has fibrosis (scars); (D) decay of CA from the fibrosis. When varying one parameter, all the others are fixed according to the infarction column in Table 1.

5.2. Limitations and Future Works

For low doses of contrast agents, there is an approximately linear relationship between the SI of the CA and its concentration (Donahue et al., 1997). For cardiovascular imaging, the dosage is not low enough for this linearity assumption to apply (Gerber et al., 2008). To overcome this issue, several researchers have presented different solutions (Cernicanu and Axel, 2006; Engblom et al., 2017; Kellman et al., 2017). In this work, the presented models focus on the dynamics of CA concentration. In future works, we will adopt the aforementioned solutions to convert CA concentration to SI and vice-versa.

As mentioned before, in the series of images provided by the contrast enhanced MRI perfusion, the clinicians take the same image, or phase, at different cardiac cycles to analyze the perfusion of the CA. Doing so, the short axis images of the left ventricle have the same shape over time. Hence it is reasonable to assume a static configuration of pressure. However, this assumption, which does not take into account poroelastic effects, has its shortcomings. To improve this model, a general mathematical model that couples porous flows and hyperelastic formulation was developed in Chapelle and Moireau (2014). Furthermore, the same authors used a similar formulation in order to evaluate myocardial perfusion (Chapelle et al., 2010). In future works, the model presented here will be extended in this direction, in order to obtain more reliable results and to evaluate the myocardium at different scenarios and for different pathologies.

It is important to mention that the pressure and flow distributions in the intravascular tree are more complex than those obtained by our model. The fine and complex structure with its connectivity, branches and varying radius are important features that are not present in our continuum porous media model. Therefore, in the near future we will replace our continuum approach for the intravascular domain by a detailed arterial tree model. In order to generate arterial trees, we will adopt algorithms based on the CCO (Constrained Constructive Optimization) method (Karch et al., 1999; Blanco et al., 2013; Brito et al., 2017; Meneses et al., 2017). This will allow us to model stenosis in a more precise way, by directly altering the radius of a particular branch of the coronary arteries.

We have chosen to use a segmentation of the short axis that includes the papillary muscles. Indeed, recent papers show that these muscles have impact on the analysis of MR images (Gommans et al., 2016). However, this controversy is out of the scope of our work. The segmentation choice was arbitrary and our model and implementation is general enough to use any domain, short-axis, long-axis, with or without papillary muscles. The same applies for the location and shapes of the modeled infarct region.

A very important feature here presented is the simulation of fibrosis. To reproduce late-enhancement images for the case of infaction it was fundamental to couple the third domain for fibrosis to the two-domain model (intra- and extra-vascular domains). However, this brings more parameters to be adjusted. In future works, we plan to study how the numbers of parameters of the models affect the correct description of patient specific perfusion data.

A similar model to the two-domain model (intra- and extra-vascular domains) presented here has been successfully used in the quantification of myocardial perfusion with cardiac MRI (Engblom et al., 2017; Kellman et al., 2017; Brown et al., 2018). It is a simplification of the model presented in Bassingthwaighte et al. (1989) and was called the BTEX model. This simplified model was shown to be very useful to speed up the inverse problem of parameter estimation to reproduce patient-specific CA dynamics. Here, we will show the simplification steps needed to convert our two-domain model into the BTEX model. First, the pressure equation is not taken into account, and tissue perfusion is characterized by a single scalar parameter F which replaces our vector $\overrightarrow{v}=-\mathbf{\text{K}}\nabla p$ in Equation (6). Another simplification of BTEX is that it is used as a set of independent 1D two-domain models. The perfusion F is independently estimated for each MRI pixel via a 1D model. This approach allows a great simplification for the problem of parameter estimation. In this work, we have presented 2D models that can be easily extended to full 3D models of the heart. Finally, a very important feature here presented was the simulation of fibrosis. To reproduce late-enhancement images for the case of infarction, it was fundamental to couple the third domain for fibrosis to the two-domain model. In summary, the new models proposed in this work have some potential advantages for the process of cardiac perfusion quantification. However, the models have more equations and parameters than the simplified models that are currently used. As a consequence, parameter estimation via the solution of an associated inverse problem will be a difficult task. One solution is to use additional non-invasive information to adjust the model, such as T1 and T2 mapping from MRI exams and FFR (pressure and flux) estimations from CT exams. As future work, we plan to study how these larger number of parameters and equations may affect the correct estimation of patient-specific cardiac perfusion by performing a clinical study similar to those presented in Engblom et al. (2017), Kellman et al. (2017), and Brown et al. (2018).

5.3. Conclusion

We presented a mathematical model that describes the blood and contrast perfusion in a short axis cut of the left ventricle of the heart. The simulation results are in agreement with the observed phenomena, specifically with the MRI images obtained by the non-invasive first pass and late enhancement exams and with clinical data on the time dependence of the contrast agent (CA) in different situations. This model takes into account the flow of Gadolinium-based MRI CA (Gd-DTPA) in the arterial vascular domain, its leakage to the interstitial domain and its retention in fibrotic tissue (when simulating a scenario of infarction). Having in mind the substantial information that quantitative methods can provide for the physicians, the mathematical model proposed here and its future extensions can be considered as a first step to develop better tools to understand the complexity of the spatio-temporal CA dynamics and the relation between clinical data acquired via MRI and the underlying cardiac perfusion.

Author Contributions

JA and RdS conceived the mathematical model. JA implemented the code and performed the simulations. JA and RdQ developed the numerical approximations. MB provided clinical data to validate the model. All authors evaluated the results and reviewed the manuscript.

Funding

The authors would like to acknowledge the financial support provided by FAPEMIG, CAPES, UFJF, FINEP, CNPq, and PTB.

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.

Acknowledgments

The authors are grateful to Tobias Schäffter for the fruitful discussions.

References

Abriola, L. M., and Pinder, G. F. (1985). A multiphase approach to the modeling of porous media contamination by organic compounds: 1. Equation development. Water Resour. Res. 21, 11–18. doi: 10.1029/WR021i001p00011

CrossRef Full Text | Google Scholar

Alves, J., de Queiroz, R., and dos Santos, R. (2017). Comparing myocardium perfusion data acquired by a MRI-phantom and a mathematical model. IFMBE Proc. 60, 308–311. doi: 10.1007/978-981-10-4086-3_78

CrossRef Full Text | Google Scholar

Alves, J. R., Queiroz, R. A. B., and Santos, R. W. (2016). Simulation of cardiac perfusion by contrast in the myocardium using a formulation of flow in porous media. J. Comput. Appl. Math. 295, 13–24. doi: 10.1016/j.cam.2015.04.008

CrossRef Full Text | Google Scholar

Alves, J. R., Santos, R. W., Novaes, G. M., and Queiroz, R. A. B. (2015). Simulating the perfusion of contrast in cardiac transmural wall. IFMBE Proc. 49, 837–840. doi: 10.1007/978-3-319-13117-7_213

CrossRef Full Text | Google Scholar

Al-Wakeel-Marquard, N., Rastin, S., Muench, F., Yilmaz, S., Berger, F., Kuehne, T., et al. (2017). Cardiac t1 mapping in congenital heart disease: bolus vs. infusion protocols for measurements of myocardial extracellular volume fraction. Int. J. Cardiovasc. Imaging 33, 1961–1968. doi: 10.1007/s10554-017-1191-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Aquaro, G. D., Todiere, G., Di Bella, G., Guiducci, L., Pingitore, A., and Lionetti, V. (2013). A fast and effective method of quantifying myocardial perfusion by magnetic resonance imaging. Int. J. Cardiovasc. Imaging 29, 1313–1324. doi: 10.1007/s10554-013-0220-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Arai, A. E. (2011). The cardiac magnetic resonance (CMR) approach to assessing myocardial viability. J. Nucl. Cardiol. 18, 1095–1102. doi: 10.1007/s12350-011-9441-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bakir, M., Wei, J., Nelson, M. D., Mehta, P. K., Haftbaradaran, A., Jones, E., et al. (2016). Cardiac magnetic resonance imaging for myocardial perfusion and diastolic function - reference control values for women. Cardiovasc. Diagn. Ther. 6, 78–86. doi: 10.3978/j.issn.2223-3652.2015.09.03

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassingthwaighte, J., Wang, C., and Chan, I. (1989). Blood-tissue exchange via transport and transformation by capillary endothelial cells. Circul. Res. 65, 997–1020. 78-86.

PubMed Abstract | Google Scholar

Blanco, P. J., Queiroz, R. A. B., and Feijóo, R. A. (2013). A computational approach to generate concurrent arterial networks in vascular territories. Int. J. Numer. Methods Biomed. Eng. 29, 601–614. doi: 10.1002/cnm.2547

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogaert, J., Dymarkowski, S., and Taylor, A. M. (2005). Clinical Cardiac MRI. Heidelberg: Taylor & Francis.

PubMed Abstract | Google Scholar

Brenner, D. J., and Hall, E. J. (2007). Computed tomography - an increasing source of radiation exposure. N. Engl. J. Med. 357, 2277–2284. doi: 10.1056/NEJMra072149

PubMed Abstract | CrossRef Full Text | Google Scholar

Brito, P. F., Meneses, L. D. M., Rocha, B. M., Santos, R. W., and Queiroz, R. A. B. (2017). Construction of arterial networks considering the fahraeus-lindqvist effect. IFMBE Proc. 60, 277–280. doi: 10.1007/978-981-10-4086-3_70

CrossRef Full Text | Google Scholar

Brown, L. A. E., Onciul, S. C., Broadbent, D. A., Johnson, K., Fent, G. J., Foley, J. R. J., et al. (2018). Fully automated, inline quantification of myocardial blood flow with cardiovascular magnetic resonance: repeatability of measurements in healthy subjects. J. Cardiovasc. Magn. Reson. 20:48. doi: 10.1186/s12968-018-0462-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Cernicanu, A., and Axel, L. (2006). Theory-based signal calibration with single-point t1 measurements for first-pass quantitative perfusion mri studies. Acad. Radiol. 13, 686–693. doi: 10.1016/j.acra.2006.02.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapelle, D., Gerbeau, J.-F., Sainte-Marie, J., and Vignon-Clementel, I. E. (2010). A poroelastic model valid in large strains with applications to perfusion in cardiac modeling. Comput. Mech. 46, 91–101. doi: 10.1007/s00466-009-0452-x

CrossRef Full Text | Google Scholar

Chapelle, D., and Moireau, P. (2014). General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes. Eur. J. Mech. B Fluids 46, 82–96. doi: 10.1016/j.euromechflu.2014.02.009

CrossRef Full Text | Google Scholar

Chiribiri, A., Schuster, A., Ishida, M., Hautvast, G., Zarinabad, N., Morton, G., et al. (2013). Perfusion phantom: An efficient and reproducible method to simulate myocardial first-pass perfusion measurements with cardiovascular magnetic resonance. Magn. Reson. Med. 69, 698–707. doi: 10.1002/mrm.24299

PubMed Abstract | CrossRef Full Text | Google Scholar

Cookson, A. N., Lee, J., Michler, C., Chabiniok, R., Hyde, E., Nordsletten, D., et al. (2014). A spatially-distributed computational model to quantify behaviour of contrast agents in MR perfusion imaging. Med. Image Anal. 18, 1200–1216. doi: 10.1016/j.media.2014.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Daly, C., and Kwong, R. Y. (2013). Cardiac MRI for myocardial ischemia. Method. DeBakey Cardiovasc. J. 9, 123–131. doi: 10.14797/mdcj-9-3-123

PubMed Abstract | CrossRef Full Text | Google Scholar

Dash, R. K., and Bassingthwaighte, J. B. (2006). Simultaneous blood–tissue exchange of oxygen, carbon dioxide, bicarbonate, and hydrogen ion. Ann. Biomed. Eng. 34, 1129–1148. doi: 10.1007/s10439-005-9066-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Donahue, K. M., Weisskoff, R. M., and Burstein, D. (1997). Water diffusion and exchange as they influence contrast enhancement. J. Magn. Reson. Imaging 7, 102–110. doi: 10.1002/jmri.1880070114

PubMed Abstract | CrossRef Full Text | Google Scholar

Engblom, H., Xue, H., Akil, S., Carlsson, M., Hindorf, C., Oddstig, J., et al. (2017). Fully quantitative cardiovascular magnetic resonance myocardial perfusion ready for clinical use: a comparison between cardiovascular magnetic resonance imaging and positron emission tomography. J. Cardiovasc. Magn. Reson. 19:78. doi: 10.1186/s12968-017-0388-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, V. G., Queiroz, R. A. B., Lima, G. A. B., Cuenca, R. G., Oishi, C. M., Azevedo, J. L. F., et al. (2012). A bounded upwinding scheme for computing convection-dominated transport problems. Comput. Fluids 57, 208–224. doi: 10.1016/j.compfluid.2011.12.021

CrossRef Full Text | Google Scholar

Gerber, B. L., Raman, S. V., Nayak, K., Epstein, F. H., Ferreira, P., Axel, L., et al. (2008). Myocardial first-pass perfusion cardiovascular magnetic resonance: history, theory, and current state of the art. J. Cardiovasc. Magn. Reson. 10:18. doi: 10.1186/1532-429X-10-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Gommans, D. H. F., Bakker, J., Cramer, G. E., Verheugt, F. W. A., Brouwer, M. A., and Kofflard, M. J. M. (2016). Impact of the papillary muscles on cardiac magnetic resonance image analysis of important left ventricular parameters in hypertrophic cardiomyopathy. Netherlands Heart J. 24, 326–331. doi: 10.1007/s12471-016-0805-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hlatky, M. A., De Bruyne, B., Pontone, G., Patel, M. R., Norgaard, B. L., Byrne, R. A., et al. (2015). Quality-of-life and economic outcomes of assessing fractional flow reserve with computed tomography angiography: Platform. J. Am. Coll. Cardiol. 66, 2315–2323. doi: 10.1016/j.jacc.2015.09.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsu, L.-Y., Groves, D. W., Aletras, A. H., Kellman, P., and Arai, A. E. (2012). A quantitative pixel-wise measurement of myocardial blood flow by contrast-enhanced first-pass cmr perfusion imaging: microsphere validation in dogs and feasibility study in humans. Cardiovasc. Imaging 5, 154–166. doi: 10.1016/j.jcmg.2011.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsu, L.-Y., Rhoads, K. L., Holly, J. E., Kellman, P., Aletras, A. H., and Arai, A. E. (2006). Quantitative myocardial perfusion analysis with a dual-bolus contrast-enhanced first-pass mri technique in humans. J. Magn. Reson. Imaging 23, 315–322. doi: 10.1002/jmri.20502

PubMed Abstract | CrossRef Full Text | Google Scholar

Jenson, S., Craig, S., York, G., and Steel, K. (2013). Review of cardiac MR perfusion imaging. J. Am. Osteopath. Coll. Radiol. 2, 2–7. Available online at: https://www.jaocr.org/articles/review-of-cardiac-mr-perfusion-imaging

Google Scholar

Jerosch-Herold, M. (2010). Quantification of myocardial perfusion by cardiovascular magnetic resonance. J. Cardiovasc. Magn. Reson. 12:57. doi: 10.1186/1532-429X-12-57

PubMed Abstract | CrossRef Full Text | Google Scholar

Karch, R., Neumann, F., Neumann, M., and Schreiner, W. (1999). A three-dimensional model for arterial tree representation, generated by constrained constructive optimization. Comput. Biol. Med. 29, 19–38. doi: 10.1016/S0010-4825(98)00045-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kellman, P., Hansen, M. S., Nielles-Vallespin, S., Nickander, J., Themudo, R., Ugander, M., et al. (2017). Myocardial perfusion cardiovascular magnetic resonance: optimized dual sequence and reconstruction for quantification. J. Cardiovasc. Magn. Reson. 19:43. doi: 10.1186/s12968-017-0355-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Knowles, B. R., Batchelor, P. G., Parish, V., Ginks, M., Plein, S., Razavi, R., et al. (2008). Pharmacokinetic modeling of delayed gadolinium enhancement in the myocardium. Magn. Reson. Med. 60, 1524–1530. doi: 10.1002/mrm.21767

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, H., Bousseljot, R.-D., Kosch, O., Jahnke, C., Paetsch, I., Fleck, E., et al. (2011). A reference dataset for verifying numerical electrophysiological heart models. Biomed. Eng. Online 10:11. doi: 10.1186/1475-925X-10-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroeker, K. I., Lam, S., Birchall, I., and Fedorak, R. N. (2011). Patients with IBD are exposed to high levels of ionizing radiation through CT scan diagnostic imaging: a five-year study. J. Clin. Gastroenterol. 45, 34–39. doi: 10.1097/MCG.0b013e3181e5d1c5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, K. E., Lee, S. H., Shin, E.-S., and Shim, E. B. (2017). A vessel length-based method to compute coronary fractional flow reserve from optical coherence tomography images. Biomed. Eng. Online 16:83. doi: 10.1186/s12938-017-0365-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Lima, J. A., Judd, R. M., Bazille, A., Schulman, S. P., Atalar, E., and Zerhouni, E. A. (1995). Regional heterogeneity of human myocardial infarcts demonstrated by contrast-enhanced MRI. Circulation 92, 1117–1125. doi: 10.1161/01.CIR.92.5.1117

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, P. K., Wei, J., and Wenger, N. K. (2015). Ischemic heart disease in women: a focus on risk factors. Trends Cardiovasc. Med. 25, 140–151. doi: 10.1016/j.tcm.2014.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Mendis, S., Puska, P., and Norrving, B. (2011). Global Atlas on Cardiovascular Disease Prevention and Control. Geneva: World Health Organization.

Meneses, L. D. M., Brito, P. F., Rocha, B. M., Santos, R. W., and Queiroz, R. A. B. (2017). Construction of arterial networks considering a power law with exponent dependent on bifurcation level. IFMBE Proc. 60, 545–548. doi: 10.1007/978-981-10-4086-3_137

CrossRef Full Text | Google Scholar

Michler, C., Cookson, A. N., Chabiniok, R., Hyde, E., Lee, J., Sinclair, M., et al. (2013). A computationally efficient framework for the simulation of cardiac perfusion using a multi-compartment Darcy porous-media flow model. Int. J. Numer. Methods Biomed. Eng. 29, 217–232. doi: 10.1002/cnm.2520

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagel, E., Klein, C., Paetsch, I., Hettwer, S., Schnackenburg, B., Wegscheider, L., et al. (2003). Magnetic resonance perfusion measurements for the noninvasive detection of coronary artery disease. Circulation 108, 432–437. doi: 10.1161/01.CIR.0000080915.35024.A9

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikolaou, K., Alkadhi, H., Bamberg, F., Leschka, S., and Wintersperger, B. J. (2011). MRI and CT in the diagnosis of coronary artery disease: indications and applications. Insights Imaging 2, 9–24. doi: 10.1007/s13244-010-0049-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Pijls, N. H., De Bruyne, B., Bech, G. J. W., Liistro, F., Heyndrickx, G. R., Bonnier, H. J., et al. (2000). Coronary pressure measurement to assess the hemodynamic significance of serial stenoses within one coronary artery. Circulation 102, 2371–2377. doi: 10.1161/01.CIR.102.19.2371

PubMed Abstract | CrossRef Full Text | Google Scholar

Pop, M., Ghugre, N., Ramanan, V., Morikawa, L., Stanisz, G., Dick, A., Wright, G., et al. (2013). Quantification of fibrosis in infarcted swine hearts by ex vivo late gadolinium-enhancement and diffusion-weighted MRI methods. Phys. Med. Biol. 58, 5009–5028. doi: 10.1088/0031-9155/58/15/5009

PubMed Abstract | CrossRef Full Text | Google Scholar

Quarteroni, A., Sacco, R., and Saleri, F. (2010). Numerical Mathematics, Vol. 37. Heidelberg: Springer Science & Business Media.

Google Scholar

Sant'Anna, F. M., Silva, E. E., Batista, L. A., Ventura, F. M., Barrozo, C. A. M., and Pijls, N. H. (2007). Influence of routine assessment of fractional flow reserve on decision making during coronary interventions. Am. J. Cardiol. 99, 504–508. doi: 10.1016/j.amjcard.2006.09.097

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuster, A., Sinclair, M., Zarinabad, N., Ishida, M., van den Wijngaard, J. P., Paul, M., et al. (2015). A quantitative high resolution voxel-wise assessment of myocardial blood flow from contrast-enhanced first-pass magnetic resonance perfusion imaging: microsphere validation in a magnetic resonance compatible free beating explanted pig heart model. Eur. Heart J. Cardiovasc. Imaging 16, 1082–1092. doi: 10.1093/ehjci/jev023

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A. F., Shipley, R. J., Lee, J., Sands, G. B., LeGrice, I. J., and Smith, N. P. (2014). Transmural variation and anisotropy of microvascular flow conductivity in the rat myocardium. Ann. Biomed. Eng. 42, 1966–1977. doi: 10.1007/s10439-014-1028-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Taber, J. (1969). Dynamic and static forces required to remove a discontinuous oil phase from porous media containing both oil and water. Soc. Petrol. Eng. J. 9, 3–12. doi: 10.2118/2098-PA

CrossRef Full Text | Google Scholar

Wang, J., Xiang, B., Lin, H.-y., Liu, H.-y., Freed, D., Arora, R. C., et al. (2015). Collateral circulation formation determines the characteristic profiles of contrast-enhanced MRI in the infarcted myocardium of pigs. Acta Pharmacol. Sin. 36, 463–472. doi: 10.1038/aps.2014.158

PubMed Abstract | CrossRef Full Text | Google Scholar

Zarinabad, N., Chiribiri, A., Hautvast, G. L., Ishida, M., Schuster, A., Cvetkovic, Z., et al. (2012). Voxel-wise quantification of myocardial perfusion by cardiac magnetic resonance. Feasibility and methods comparison. Magn. Reson. Med. 68, 1994–2004. doi: 10.1002/mrm.24195

PubMed Abstract | CrossRef Full Text | Google Scholar