Modeling of the contrast-enhanced perfusion test in liver based on the multi-compartment flow in porous media

Modeling blood perfusion in tissues belongs to the most challenging problems in the biomechanical and biomedical research. The topic is of interest especially in liver surgery (Plantefèv et al, 2016), brain neurosurgery, or neurology (Dai et al, 2016; Reichold et al, 2009; Reichold, 2011). Although the functionality of the two related organs, liver and brain, is completely different, there are some similar aspects for both. In contrast with the myocardium, another highly perfused tissue, for both the organs the phenomenon of deformation can be neglected unless special concerns are in focus, see e.g.  Peterlik et al (2012). It is desirable to find accurately parts of the tissue with an insufficient blood supply, to localize anomalies in the blood microcirculation, and to quantify locally the perfusion efficiency. It is worth to note that also the kidney belongs to highly perfused tissues with its important role in the blood circulation system, but we do not tackle this issue in the present paper.

In most perfused organs, the blood is transported through complex hierarchical branching systems of vessels with various diameters which diminish consecutively at each bifurcation, beginning from arteries down to capillaries. The ascending system of veins is arranged in the reversed sense, starting at the capillary level, where the blood participates in the metabolism process, delivering the oxygen to the tissue and collecting the waste products. The blood circulation in the tissue can be represented by two penetrating, but mutually separated vascular trees connected only at the capillary level. In reality, collateralization vessels exist which make shunts between the two trees, but this phenomenon is usually negligible for the global hierarchically built system.

Although the modeling approach is presented in a rather general setting, in this paper, we focus on the liver tissue which has a special structure. The blood flows in the whole organ from two separated vascular trees, one connected to the hepatic artery (HA) and the other to the portal vein (PV), whereby the venous supply dominates. The two branching systems bifurcate to about 20 generations (Debbaut et al, 2014). The blood is filtered in the tissue functional units, called hepatic units, or lobules, which typically are considered as hexagonal prisms. At their corner edges, the hepatic artery and portal vein trees terminate and the blood flows through hepatic capillaries constituting the sinusoids to the so-called central veins, the terminal branches of the draining hepatic vein (HV) located in the center of each lobule. In the proposed model, for simplicity we neglect the HA, as its blood supply is minor in the comparison with the one by the portal vein¹¹1The perfusion model is intended to estimate the volume regeneration capacity of the liver parenchyma (Brůha et al, 2015). In this context, it is important to capture the supply through the portal vein. Concerning the perfusion CT simulation, we focus on the second stage when the contrast bolus arrives at the PV inlet so that the first stage associated with the HA is over..

The steadily growing body of published works related to the modeling of liver perfusion shows its importance for the medical research. Obviously, relevant computational models which can be used for complex simulations involving various parameters require detailed information about the liver structure. However, there is a pertaining difficulty in describing the vascular trees bridging several scales. The current standard computed tomography (CT), magnetic resonance imaging (MRI) and related image processing techniques provide resolution of the vascular trees down to millimeters, which is not enough to capture details on about 10 underlying tree generations. Therefore, various methods and algorithms have been proposed to generate artificial, but anatomically correct vascular architectures. These are coupled with the authentic patient-specific reconstructed geometry to obtain complete vascular trees with a required resolution of microvessels. The most popular approach is based on the constrained constructive optimization (CCO), see e.g. Georg et al (2010), Schwen and Preusser (2012). However, there are more advanced approaches based on the metabolism driven processes of the vascular tree reconstruction (Schneider et al, 2012).

Rather than modeling the flow on the complex branching system, many publications have been devoted to detailed studies of microcirclation. Perfusion characteristics of blood flow through hepatic microcapillaries of human lobular sinusoids were analyzed by Debbaut et al (2012) using tissue samples processed by corrosion casting and the micro-CT imaging; anisotropic tensors of the sinusoidal permeability were calculated by a simple averaging of the flow field obtained using standard CFD simulations. Blood circulation in the liver lobule was considered in Bonfiglio et al (2010) using both Newtonian and shear-thinning flow models. The liver tissue is able to adapt rapidly to modified hydraulic conditions; in Ricken et al (2010), the sinusoidal permeability has been studied in terms of a vascular remodeling process induced by outflow obstruction. Sinusoidal vascular canals were assumed to align with the blood pressure gradient. In Ricken et al (2015), using the general framework of a two-scale approach and the theory of porous media, the basic coupled transport-reaction mechanisms between blood perfusion and the hepatic cell metabolism were described at the level of lobule. Lymphatic flow production in the liver lobules arranged in periodic lattices was treated in Siggers et al (2014). A similar approach has been used in Lettmann and Hardtke-Wolenski (2014) to study microcirculation driven transport of oxygen and cytokine in the tissue. Besides the blood flow in healthy liver itself, diseased tissue deserve attention from medical point of view, but also as a challenge for mathematical modeling; there is a considerable number of issues related to modeling pathophysiological conditions, such as the fluid transport through vascular tumors (Pozrikidis, 2010; Shipley and Chapman, 2010; Mescam et al, 2010).

Since the straightforward mathematical description of blood flow on the complete vascular trees featured by several hierarchies using the Navier-Stokes equations leads to a prohibitively complicated numerical model, some convenient model reduction is needed. We adhere to the general approach proposed in Formaggia et al (2009) which consists in decomposing the tree into a higher (HH) and a lower (LH) hierarchy, see Fig. 1.

The HH tree is constituted by larger vessels, so that the flow can be described in the 3D geometry using the Navier-Stokes equations. Alternatively the vessels can be replaced with line segments, so that the flow model is obtained using the 1D description of flow on such branching system of line segments.

To describe the blood flow in the LH tree, we resort to a type of the Darcy flow model which is defined in the 3D bulk of the liver parenchyma. Since the LH part involves vessels associated with different hierarchies, several compartments are introduced and Darcy model is used within each of them. The flow between different compartments is described in terms of the difference between the compartmental pressure at a given location; details are given below.

The simulation of blood perfusion in tissues, such as the hepatic or cerebral ones, belongs to problems typically requiring a sort of multiscale modeling approach (D’Angelo, 2007). The term “multiscale” is employed in a slightly different manner than in Formaggia et al (2009), where the “geometrical multiscale” modeling was introduced in the context of the cardiovascular system. The difficulty of blood perfusion modeling arises from the nature of the flow in the perfusion trees incorporating blood vessels of very different diameters. In liver, the blood is conveyed to the organ by the portal vein which has the diameter of about $D=1\,\textrm{cm}$. The perfusion trees have several hierarchies associated with bifurcations, where vessel diameters reduce. The tissue parenchyma is arranged in hexagonal structures (lobules) with the characteristic size of 1.5 mm where the blood is transported through the capillaries with the diameter of about $d=10\,\mu\textrm{m}$ from the portal compartment to the hepatic one. Thus, the overall scale change characterized by the ratio of the largest and smallest diameters of the vascular trees is approximately $d/D=0.001$.

The model which we develop should allow also for simulation of the transport of the contrast fluid (the tracer) during the contrast-enhanced dynamic perfusion test. The aim is to provide a computational feedback which would enable us to analyze more accurately the CT scans obtained from the standard dynamic perfusion test of a patient (Materne et al, 2000). The methods currently being used are based on the deconvolution, or maximum slope techniques. They provide some local integral characteristics like the blood flow, blood volume, mean transition times and others for each voxel of the tissue (Koh et al, 2006; Fieselmann et al, 2011). Our approach should provide an alternative and more detailed interpretation of the measured CT scans. The proposed strategy announced first in our conference paper (Jonášová et al, 2014) is based on the following tasks:

• •

  to describe flows in the perfusion trees including the parenchyma by combining suitable models which are associated with different hierarchies and which use the most information about the structures and geometry, involving only few undetermined parameters;

• •

  to reconstruct the organ (liver) shape and the blood vessels geometries up to a certain hierarchy using the image segmentation techniques (based on the CT “static” data) (Ryba et al, 2013; Jiřík, 2013–2015) possibly supplemented by the CCO approach employed also in Lukeš et al (2014);

• •

  to identify selected parameters of the perfusion model (like permeability), so that the simulated dynamic CT examination well approximates the measured data (Rohan et al, 2015a).

Such a “tuned” model would enable us to analyze the blood flow in particular compartments and to predict effects of intended medical treatment, like resection of a part of the tissue.

In the proposed strategy, several difficulties arise: on one hand, the model should reflect the microstructure and fit the complex geometry of the perfusion trees, on the other hand, it should be parametrized using just not too many parameters, to prevent ill conditioning of the “tuning” step involving simulations of the steady blood flow and the dynamic tracer transport. Therefore, the tuning parameters must reflect some important features of the perfusion system.

The model of the perfusion and of the contrast fluid transport described in this paper consists of the following two parts:

1. 1.

   The “inlet” and “outlet” trees which are associated with the portal and hepatic vein networks²²2We have in mind the liver perfusion modeling, so that we also adapt the notations to fit with this particular application. (in the case of the liver; recall that we have neglected the hepatic artery in the model). These trees are denoted by $\mathcal{T}_{P}$ and $\mathcal{T}_{H}$ (label $P$ stands for the portal vein system, whereas label $H$ stands for the hepatic vein system). They are formed by vessel segments and junctions representing the bifurcations. The steady flow on these trees is described by a 1D model which is introduced in Section 2.3, see also Jonášová et al (2014). The model is relevant for higher hierarchies of the vascular trees, down to a certain size of the vessels, so that the terminal branches³³3or leaves in the graph theory terminology of $\mathcal{T}_{P}$ and $\mathcal{T}_{H}$ may have diameters about 1 mm.

2. 2.

   The tissue parenchyma including lower hierarchies of the “inlet” and “outlet” trees. The flow in this part of the vascular system is governed by the multi-compartment Darcy flow model in a porous medium. At any point of the liver occupying domain $\Omega\subset\mathbb{R}^{3}$, several (at least two) different pressure values are defined which are associated with different compartments, see Section 2.2. There are point sources and sinks defined in $\Omega$, where the 1D trees are connected to the parenchyma model.

In this section we explain how to simulate the dynamic perfusion tests which are used as a principal method to assess blood flow in highly perfused organs, e.g.  in the liver, or brain. It is based on the computed tomography scans which provide information on the tissue density. This quantity is proportional to the local concentration of the contrast fluid (the tracer) which is dissolved in the blood. Its relative content is expressed by the saturation $S$. This quantity is defined for each of the individual compartments of the parenchyma treated by the Darcy flow model and also for all branches of the upper level vascular trees, as explained in Section 2.3.

Transport equations for resolving the saturations in all the compartments can be derived using the pre-computed perfusion velocities, the mass conservation law, and taking into account fluid exchange between the compartments. Thus, we obtain a system of hyperbolic equations for resolving the saturations. Then the tissue contrast is defined locally as the weighted sum of all the saturations; the weights are given by the volume fractions.

The model introduced in the preceding sections has been implemented in our non-commercial codes. In this section, we illustrate the model response using numerical examples with real geometry of the human liver, but with the perfusion trees of the portal and hepatic veins generated using the CCO method (Georg et al, 2010). Before presenting particular results which serve as a proof concept, we explain a flowchart of the currently developed computational modeling tool.

CT scans of human liver processed by the LISA (LIver Surgery Analyser) software (Jiřík, 2013–2015) are the starting point to the numerical simulation of liver perfusion. The LISA code includes the semi-automated segmentation method based on the graph-cut algorithm and returns the geometrical model of liver parenchyma that can be transformed into a volumetric finite element mesh. The code is also used to identify the position and diameter of the portal and hepatic veins at the point of entry into the liver parenchyma; this data is necessary to generate two artificial vascular trees associated with the portal and hepatic veins, which in part are employed directly in the simulation of the upper-level flow using the 1D flow model described in Section 2.3, but also to establish permeabilities ${\bm{K}}^{i}$ and perfusion coefficients $G_{i}^{j}$ appearing in the multicompartment Darcy flow model. The hepatic artery is not considered in our liver model because of its small contribution to the overall blood perfusion through the liver tissue. A constrained constructive optimization approach is used for the creation of vascular structures reflecting main physiological principles, see Georg et al (2010), and the geometry of a particular liver.

In accordance with the modeling principles described above, the calculation of blood perfusion and contrast fluid distribution in the liver can be divided into two stages. The first one includes the solution of 1D flow model (11)-(12) coupled with the multicompartment Darcy model (10). The resulting velocities associated with vessel segments of the 1D flow model and velocity fields ${\bm{w}}^{i}$ in compartments $i$ are then employed in the second stage in which the compartmental saturations $S^{i}$ yielding the total concentration $C$ are computed. The 1D flow model is implemented in the Python language, the multicompartment Darcy flow model is discretized using the FE method in software SfePy – Simple Finite Elements in Python, see Cimrman (2014). The model describing the tracer transport is implemented using the finite volume method in the Matlab system. The overview of the simulation process is illustrated in Fig. 4.

To demonstrate features of the hierarchical model described in this paper, a three compartment model is considered. The first compartment represents the vascular network of the portal vein consisting of veins with diameter in the range of about $10^{-3}-10^{-4}$ m. The same diameter range is chosen also for the third compartment which stands for the system of the hepatic vein. Compartment 2 involves small veins of diameters smaller than $10^{-4}$ m including the hepatic capillaries of the lobular sinusoids, and it may therefore be perceived as a filtration part of the perfusion model. The generated artificial vascular trees representing the vascular networks in the liver are depicted in Fig. 5, their properties are summarized in Tab. 1. The higher hierarchy trees, where the blood flow is approximated by the 1D model, see Fig. 1, are obtained by taking the vessel segments with the Horton–Strahler (HS) number (Georg et al, 2010) greater or equal to 6. The branching complexity of the generated portal vascular tree characterized by the HS number is illustrated in Fig. 6.

To determine parameters ${\bm{K}}^{1}$, ${\bm{K}}^{3}$ and $G_{1}^{2}$, $G_{2}^{3}$, as well as the porosities $\phi^{1}$, $\phi^{3}$ in each finite element, a volume averaging is performed over the corresponding vascular tree for a given range of the vessel diameters. The anisotropic permeabilities ${\bm{K}}^{1}$, ${\bm{K}}^{3}$ of compartments 1 and 3 are evaluated using the expression derived by Huge and van Campen (Huyghe and Campen, 1995), whereas the coupling perfusion coefficients $G_{1}^{2}$, $G_{2}^{3}$ are calculated with the help of the theory presented in Hyde et al (2013). Values of the model parameters at two selected points are listed in Tab. 2. The isotropic permeability in compartment 2 is taken as $K^{2}=2\times 10^{-14}$ and the uniform porosity as $\phi^{2}=0.15$, both values are chosen with respect to the parameters identified in Debbaut et al (2012). The domain of the liver is discretized by a finite element mesh consisting of 11 019 tetrahedrons and 2 442 nodes. The rather coarse mesh was chosen due to large computational cost associated with the tracer transport calculation.

The 1D flow model is connected to the multicompartment model through 45 junctions, 34 of them are the sources located in compartment 1 and remaining 11 are the sinks located in compartment 3. The upper-level vascular trees, terminal points of which correspond to the sources and sinks, are shown in Fig. 7 and Fig. 8 (middle). The input blood velocity to the portal vein is taken as $v_{in}=0.25\,{\textrm{m}}\cdot{\textrm{s}}^{-1}$ and the output blood pressure prescribed in the hepatic vein is $p_{out}=10^{3}$ Pa.

The steady state results computed for the three compartment perfusion model described above are presented in Fig. 8. This figure shows the distribution of pressure $p^{i}$ in each of the three compartments, $i=1,2,3$ (1 – portal, 2 – filtration, 3 – hepatic) as well as the corresponding velocity vectors ${\bm{w}}^{i}$, which well illustrate an interplay between the filling and draining functionality of the portal and hepatic compartments, respectively. From the surface and volumetric visualizations of the computed pressure fields (left, middle), it can be noticed that the portal (top) and hepatic (bottom) pressures depend strongly on the structure of the vascular trees, especially their upper-level hierarchies; this is in accordance with observations made on real human liver.

To get a better idea on how the different perfusion properties of the three compartments influence the final distribution of blood within the liver and also to simulate a real dynamic perfusion test, a simulation of the contrast fluid transport is carried out using the model presented in Section 3. For this purpose, a given external source saturation $S_{\textrm{in}}(t)$ is prescribed in the form of a time bolus defined at the root segment of the portal vein system, Fig. 9,

where $\bar{S}=0.4$ and $T=2\,\textrm{s}$. In Fig. 11, the obtained portal $S_{pv}$, lobular $S_{lob}$ and hepatic $S_{hv}$ saturations and total apparent concentration $C$ are shown in two selected sections of the liver, positions of which are outlined in Fig. 10 and denoted as planes A and B. Due to the time-dependent character of the tracer transport simulation, the results in Fig. 11 are visualized at three time instants ($t_{1}=1.17\,\textrm{s}$, $t_{2}=1.95\,\textrm{s}$, and $t_{3}=8\,\textrm{s}$). Note that the first two time instants reflect the time at which most of the portal $S_{pv}$ and lobular $S_{lob}$ saturations attain their maximum, see graphs in Fig. 12. The top graphs show the time development of the portal $S_{pv}$, hepatic $S_{hv}$ and lobular $S_{lob}$ saturations at points $\hat{\textrm{A}}$ and $\hat{\textrm{B}}$ located in planes A and B (Fig. 10), while the bottom graph depicts the evolution of the total apparent concentration $C$ at the same points.

The saturation isocontours presented in Fig. 11 clearly demonstrate the transport of the contrast fluid through each of the three considered liver compartments. In other words, the most apparent saturation of the respective compartment occurs successively—first in the portal vein system at the time $t_{1}=1.17\,\textrm{s}$, followed by the lobular one at $t_{2}=1.95\,\textrm{s}$ and then concluded in the hepatic vein system at $t_{3}=8\,\textrm{s}$. This dynamics of the perfusion can also be recognized from the varying grey levels of the total apparent concentration $C$ in Fig. 11.

To get another, more time continuous view of the tracer transport through each of the three compartments, we refer to the graphs in Fig. 12; from there it is possible to observe a relatively fast “filtration” of the tracer through the portal compartment ($S_{pv}$ reaches its maximum shortly after the saturation bolus have entered the compartment; moreover, the shape of the bolus remains almost unchanged over time, see Fig. 12 (top left). On contrary, the filtration compartment, which is characterized by a relatively high porosity ($\phi^{2}=0.15$) compared to the ones computed in the filling/draining compartments ($\phi^{1}$ and $\phi^{3}\approx 10^{-4}$, see Tab. 2), fulfills its function as a filtration system, in which the saturation $S_{lob}$ decreases only slowly over time, see Fig. 12 (top right). As a consequence, the perfusion properties of the filtration system influence the time development of the hepatic saturation $S_{hv}$ which, compared to the portal saturation $S_{pv}$, changes much more gradually depending on tracer transfer from the filtration system.

In the paper, we proposed the multicompartment model of the tissue blood perfusion and of the superimposed time-space distribution of the contrast fluid. This latter feature enables to establish a modeling feedback which is necessary to tune the model parameters according to information available due to standard clinical CT examination combined with essential flow parameters and basic structural and morphological data of a specific perfused tissue. Although the model is intended and currently being developed for modeling the liver perfusion, it can be adapted for other applications, namely the cerebral perfusion; this issue is in our focus for future work.