TOWARDS MICROSTRUCTURE BASED TISSUE PERFUSION RECONSTRUCTION FROM CT USING MULTISCALE MODELING

The paper deals with modeling of blood perfusion and simulation of dynamic CT investigation. The flow can be characterized at several scales for which different models are used. We focus on two levels: flow in larger branching vessels is described using a simple “1D” model based of the Bernoulli equation. This model is coupled through point sources/sinks with a “0D” model describing multicompartment flows in tissue parenchyma. We propose also a model of homogenized layer which can better reflect arrangement of the microvessels; the homogenization approach will be used to describe flows in the parenchyma by a two-scale model which provides the effective parameters dependent on the microstructure. The research is motivated by modeling liver perfusion which should enable an improved analysis of CT scans. For this purpose we describe a dynamic transport of the contrast fluid at levels of the “1D” and “0D” models, so that time-space distribution of the so-called tissue density can be computed and compared with the measured data obtained form the CT.


Introduction
Problem of tissue perfusion reconstruction is quite challenging for the biomechanical and biomedical research.It is desirable to find accurately parts of the tissue with an insufficient blood supply, to localize anomalies in the blood micro-circulation and to quantify locally the perfusion efficiency.Such information is need namely to asses physiological functionality of the brain tissue and trends to pathologies.For the liver segmentation, it is important to identify parts of the whole organ which are supplied with blood by a selected vein branch.Our effort is to improve modeling techniques for processing the standard CT investigations.
Microstructurally-oriented and computationally efficient models allow for simulations of tissue perfusion and transport of the contrast fluid whose instantaneous volume concentration in the parenchyma corresponds to the density registered by the CT.To cope with the complexity of the perfusion trees, we use a combination of 3D and 1D fluid dynamics for large and medium vessels, respectively, whereas Darcy flow "0D" models are used to describe blood redistribution in the parenchyma and smaller vessels.In this paper we present a simplified model where only a 1D model is employed to describe flow on the branching medium vessels.
We developed several multiscale and multicompartment models based on homogenization [8,10,11];they involve effective parameters computed using solutions of local "microscopic" problems in RVEs where the geometry and topology of different vessels is specified; each compartment is associated with one pressure field, their local differences reveal the amount of the tissue perfusion.To describe the hierarchical arrangement of the perfused tissue, we proposed a model of homogenized perfusion in the 3-compartment medium constituted by several transversely periodic layers with dual porosities [9].Thus, a 3D volume can be replaced by a number of 2D "homogenized layers" coupled by conditions governing the fluid exchange between them.The 1D flow model on the branching vessel network is coupled with the homogenized 0D flow model of the tissue perfusion by a non-compound solving procedure -the sheared data (pressures and fluxes) are being exchanged between the two models on point and line interfaces defined in the 3D domain.Once the perfusion fields are computed, the tracer redistribution can be simulated using a transport model.The local tracer concentrations in different compartments are averaged to obtain the density maps comparable with the CT outputs.
The paper is organized, as follows: First, in Section 2 we explain conception of a hierarchical flow model in the context of the liver perfusion simulation.For modeling at the whole organ level we use a lumped model with multiple sectors and ad hoc defined permeabilities and perfusion parameters, as discussed in Section 2.1.In Section 3 we report on the two-scale model of homogenized perfused layer; it is outlined how this model can be used to describe perfusion of liver lobules at the lowermost level where the blood is transported between the portal and hepatic veins.In Section 4 we describe a model of the transport of the tracer through the branching network of portal and hepatic veins

Hierarchical model of perfusion
The simulation of blood perfusion in tissues like liver or brain belongs to problems that necessitate a sort of multiscale modeling to be used, cf.[4];the term "multiscale" is employed in a slightly different manner than in [5],where the "geometrical multiscale" modeling was introduced in the context of the cardiovascular system.Obviously, the difficulty of the problem consists in the nature of the flow on branching structures incorporating blood vessels of very different diameters.For instance, the blood is conveyed to the liver by the portal vein with its diameter of D=1.3 cm.At the tissue parenchyma level, the blood channels are formed in the lobular hexagonal structures with the characteristic size of 1.5 mm, so that the sinus through which the blood flows from the portal to the hepatic compartments is formed by microvessels with the diameter of about d=10 µm.Moreover the scale change characterized by the ratio d/D = 0.001 is continuous, the perfusion trees have several hierarchies distinguished by the number of "new" branches and reduced vessel diameters.
The model which we develop should allow also for simulation of the transport of the contrast fluid (the tracer) during the dynamic perfusion test.The aim is to provide a computational feedback which would enable to analyze more accurately the CT scans obtained from the standard dynamic perfusion test of a patient, cf.[7].Nowadays the methods based on the deconvolution, or maximum slope techniques provide some characteristics like the blood flow, blood volume, transition times and others for each voxel of the tissue, cf.[6].Our approach should provide an alternative and more detailed interpretation of the measured CT scans.The proposed strategy is based on the following tasks: • to describe all involved hierarchies of the perfusion trees including the parenchyma by suitable models involving 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); • to tune (by solving some optimization problems) selected parameters of the parenchyma model (like permeability), so that the simulated dynamic CT test is as close as possible to the measured data.
Such a "tuned" model would enable to analyze the blood flow in particular compartments and to predict effects of intended medical treatment, like resection of a part of the tissue.The proposed strategy contains several hurdles; 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 parameterized just using not too many parameters to prevent ill conditioning of the "tuning" step involving simulations of the steady blood flow and the dynamic tracer transport.For this purpose, the tuning parameters must reflect some important features of the perfusion system.
In this paper we describe a simplified model on which we test the modeling of the perfusion and tracer transport.It consists of the following parts: 1.The "inlet" and "outlet" trees which (in the case of the liver) described the portal and hepatic veins, respectively.These trees are formed by pipes and junctions representing the branching.For simplicity we consider incompressible, inviscid steady 1D flow, as described in Section 2.2.
2. The parenchyma "0D" model is constituted by equations governing a multicompartment Darcy flow in a 3D porous medium involving a pressure field associated with each compartment.Thus, at any point of the organ domain Ω ⊂ R 3 several (at least two) different pressure values are defined which are associated with different compartments, see Section 2.1.There are point sources and sinks defined in Ω where the "1D" trees are connected to the "0D" model.It is worth noting that our definition of a "0D" model is different from the one used in [5] We intend to describe the blood flow in the parenchyma using a more accurate model based on homogenization of the tissue microstructure.In Section 3 we report on a model which can be adapted to describe flows at the level of lobules forming a periodic structure.

Macroscopic lumped modeling of parenchyma
We approximate the perfusion at the level of tissue parenchyma using a lumped "macroscopic" model describing parallel flows in multiple compartments, cf.[2].The model involves pressures {p i } i , i = 1, . . ., ī which must satisfy the following equations: where K i is the local permeability of the i-th compartment network and G i j is the perfusion coefficient related to compartments i, j, so that J i j describes the amount of fluid which flows from i to j (drainage flux; obviously J j i = −J i j ) and f i is the local source/sink flux of the i-th compartment.Equations (1) are supplemented by the non-penetration conditions, In our numerical tests we considered just two compartments, ī = 2.The local source/sink functions for each compartment is defined using the Dirac distributions, see (4).They establish the interface between the "0D" and "1D" models, as explained below.

Flow on branching networks
Let T ({J j } j { e } e ) a branching tree formed by pipes e and junctions J j .By J 0 we denote the input junction, whereas the terminal branches end by junctions Ĵk through which they are connected with the parenchyma.Any junction J j = {e} of T joins several vessels (pipes) e , although the input and terminal junctions are just one-element sets.Further by n we denote the number of all terminal junctions.
The simplest possible model of flow on T is presented by the following n+1 equations where (3) 1 are the Bernoulli equations and (3) 2 is the mass conservation; by A k we denote the cross-sections of terminal and input branches, i.e. of the associated vessels (we use a model with constant cross-sections along each vessel which is a relevant simplifying assumption with no effect on the contrast fluid transport, as will be explained in Section 4).Now the problem is to compute the input pressure p 0 and the terminal velocities {w k } k , k = 1, . . ., n for a given input velocity w 0 and terminal pressures {p k } k .
Although this model neglects the dissipation and may lead to unaccurate results in general, in this paper we use it just for its simplicity, since the main focus is in modeling the parenchyma perfusion.Anyway, it can be extended by some energy loss terms related to the tree geometry.

Coupled "1D" -"0D" model
The "1D" model is coupled with the "0D" model through the terminal junctions which specify the sources and sinks f i for all the considered compartments.For the i-th compartment saturated by the tree T we define where δ(x − xk ) is the Dirac distribution at point xk ∈ Ω which is associated with the terminal junction k.Condition (4) 2 couples the pressures at the terminal junctions of the "1D" model with the pressure fields in Ω.In practice, we use an approximation of δ(x− xk ) which is based upon the finite element discretization.
We conclude by a simple iterative algorithm used to compute the perfusion pressure and velocity fields.Since there are two trees T , as illustrated in Fig. 5, we shall label the corresponding solutions by indices P and H, denoting quantities associated with the portal and hepatic veins.For given values wP 0 and pH 0 , i.e. the velocity in the portal (inlet) vein and the pressure in hepatic (outlet) vein, the computation proceeds by repeating the following steps: 1. Set all interface velocities and pressures to zero, namely {p P k } = 0 and {w H k } = 0. Set i = 0 and τ = 1/N , for a given N ∈ N. 5. Use the conditions (4) to update the interface variables for the next iteration i + 1.

For the new iteration
6. Go to step 2, unless a steady state is obtained.

Model of perfusion homogenized double-porous layers
Our aim is to develop a microstructurally-oriented and computationally efficient modeling tool which will allow for simulations of tissue perfusion.Here we focus on modelling the tissue parenchyma.At the level of small vessels and microvessels, the perfusion can be described using the Darcy flow in double porous structure consisting of 3 compartments: two mutually disconnected channels (small arteries and veins) and the matrix (microvessels and capillaries), represented as the dual porosity, where the permeability is decreasing with the scale parameter -the size of the microstructure.
We use techniques based on the asymptotic analysis by the periodic unfolding method [3] and treat a specific microstructure topology which gives rise to several macroscopic pressure fields, like in [14],cf.[10], satisfying specific homogenized equations.The layer-wise decomposition is a new feature in the context of homogenized models, up to our knowledge.A 3D layered structure occupying domain Ω H ⊂ R 3 can be replaced by a finite number of 2D "homogenized layers" Γ 0 ⊂ R 2 coupled by conditions governing the fluid exchange between them, as proposed in [12].

Macroscopic equation for single layer
The homogenized problem for pressures p A and p B , associated with two channels A and B, describes 2D parallel flows in homogenized layer Γ 0 ⊂ R 2 , cf. [14].Each channel system forms a connected domain (so, we assume at least a small co-lateralization of vessels in the perfusion tree).

Microstructure and the reference cell.
The periodic microstructure is generated by the reference periodic cell, Y = Ξ × I z , where the channels A and B, respectively, whereas Y M represents the dual porosity (a network of capillaries).The structure is periodic w.r.t.coordinates y := (y 1 , y 2 ) (this establishes the notion of the Y-periodicity which is refered to below), the transversal coordinate is denoted by z.The upper and lower boundary segments for D = A, B, M .The channels have inlet / outlet branches entering through the layer faces Γ + and Γ − .The inlet / outlet "surfaces" For instance, in Fig. 1 each of the two depicted channels A, B has three such surfaces.

Macroscopic equations.
Two coupled "macroscopic" equations (one for A and one for B) involve the following homogenized coefficients: permeabilities (K ab ) A,B of the channels, the transmission G and drainage (S a ) A,B,k (for channel branches k ∈ J D ) coefficients.They govern the fluid redistribution between the two channel systems A and B.
Perfusion in the homogenized layer occupying the domain Γ 0 is governed by two equations, each per one channel system (in the paper the summation convention applies for indices a, b ∈ {1, 2} of the coordinates x = (x 1 , x 2 )) where fluxes ĝ+/− , ḡD and gk D must be given such that the following solvability conditions hold: Above F A+/− , c hA are constants (the summation w.r.t.repeated indices a, b applies).The term G p A − p B , evaluated at point x ∈ Γ 0 , expresses the amount of fluid (blood) perfused through the matrix (the dual porosity) between sectors A and B. The details are reported in [9].A problem involving ( 6) and ( 7) is supplemented by boundary conditions on ∂Γ 0 representing the "side" boundary of the layer after the problem dimension has been reduced from 3D to 2D.As the result of homogenization, the original "no-flow" (Neumann type) condition in 3D is replaced by the following condition

Microscopic problems.
The homogenized coefficients are evaluated using the characteristic responses (also called the "corrector basis functions") which solve microscopic problems defined in subsectors Y D , D = A, B, M of the decomposed cell Y .By ∂ # Y ⊂ ∂Y we denote the "periodic" boundary of Y .There are four groups of the microscopic problems.
where K = (K ij ) is the permeability in the channel D.
The solution to ( 9) is the characteristic response of the pressure in channel Y D for the imposed "unit" pressure gradient, whereby no flow through the external boundary ∂Y D \∂ # Y D is admitted.Corrector functions π b D determine K D ab , see [9].2nd micro-problem.Find p g,D (Y-periodic functions in y 1 , y 2 ) such that for D = A, B and given gk For solvability of (10), fluxes gk D must satisfy (7) 1 .Solution p g,D presents the pressure response in channel Y D to external prescribed fluxes gk D .The characteristic response to "unit" fluxes are given by corrector functions γ k D which satisfy the decomposition and determine coefficients S D,k a .For more details on computing γ k D and S D,k a we refer the reader to [9].
3rd micro-problem.Find η A (Y-periodic function in y 1 , y 2 ) such that where κ = (κ ij ) is the permeability in the dual porosity Y M .Thus, for a given scale ε 0 , the "true" permeability is ε 2 0 κ.The solution of (12), i.e. the pressure in the dual porosity, is the characteristic response to the unit pressure prescribed on walls of Y A , whereas zero pressure is prescribed on walls of Y B .No flux is admitted through the external face ∂ ± Y M .Corrector function η A is employed to compute the homogenized coefficient G, see [9].
4th micro-problem.Find γ +/− (Y-periodic functions in y 1 , y 2 ) such that Given data for the simulation of perfusion in the double-porous layer.

Numerical illustration
We consider a square representing the mid-plane segment of the associated layer.The periodic structure depicted in Fig. 1, which is used in our computations, is generated by the cell Y where both Y A and Y B are channels with 3 external branches; note that both channels form connected domains in the entirety of the layer.It is important to note that the homogenized layer model with the double porosity has a real meaning for a defined scale ε > 0 which also is related to the layer thickness.To use the homogenized model established by (6), the fluxes describing channel branch saturations {g k A }, {g k B } and external flows through matrix interface ĝ+ , ĝ− must be given such that (7) are satisfied.Consider a "real structure" with given fluxes {g l,ε A }, {g k,ε B } and ĝ±,ε , see Fig. 2. Eq. ( 7) 1 holds if where the last definition expresses proportionality of the matrix fluxes to the smallness of the dual permeability, as discussed in [9], see also [12].Fluxes displayed in Fig. 2 were introduced two macroscopic points: such that (7) 2 holds, whereby ĝ±,ε = 0, thus, the dual porosity (the matrix part) was insulated on the upper and lower edges of the layer.In Fig. (3) the distribution of both macroscopic pressures in the homogenized layer is depicted and the pressure field in the dual porosity of the microstructure is depicted for two different macroscopic positions.
The multiscale modeling based on the homogenization method is implemented in our in-house developed code SfePy, see [1].

Adaptation of the model for modeling of the liver parenchyma
As explained in the introduction, the homogenized layer model was proposed to cope with perfusion on complex branching structures whereby the fluid can flow from one to another through the dual porosity.The idea was to decompose a 3D volume into "parallel" layers having periodic, but mutually different structures.In the context of the liver parenchyma, the discussed model of one homogenized layer have a special meaning.On the microscale level the tissue is formed by an almost periodic structure generated by lobules, see Fig. 4; they are arranged as hexagonal structures with a limited thickness.The two channel systems A and B represent the portal and central hepatic veins, whereas the dual porosity (the matrix) can represent the lobule sinus.Other compartments (the hepatic artery and the bile ducts) are neglected.
There are two possible methods how to employ the perfusion layer model within the ad hoc macroscopic layer model.The first one is more coherent with the original idea of decomposing a 3D volume into a number of homogenized layers, which means the dimension reduction "3D to N times 2D".The second one consists in proposing a lumping scheme which associates the structure of the two-compartment layer model ( 6) to with the analogical volume perfusion model ( 1) defined with respect toa coordinate system reflecting locally the tissue orthotropy: namely the pressure gradients (and associated fluxes) in the transversal directions must be established for a layer with its thickness proportional to the scale ε > 0, as discussed above in Section 3.2.To develop such a lumping scheme is a challenging task which may be of interest for many other applications.

Convected contrast fluid in the porous materials
In this section we explain how to simulate the dynamic perfusion tests which are used as principal method to assess blood flow in the brain or in the lever.It is based on the CT (computed tomography) investigation which provides scans of the tissue density.This quantity is proportional to the local concentration of the contrast fluid (the tracer) transported being dissolved in the blood.We assume that the tracer does not diffuses in the solution, so it is only convected.Its content in the solution is expressed by the saturation S.

Tracer redistribution in parenchyma
Let S i be the tracer saturation of the i-the compartment, so that c i = φ i S i (no summation) is the tracer partial concentration associated with i.Clearly, S i ∈ [0, 1]; this boundedness must be guaranteed by the transport equations (the conservation law).The total apparent concentration (the grey level) is then given as The local conservation in a domain ω ⊂ Ω for the i-th compartment is expressed, as follows: where S in is the external source saturation, f i + > 0 is the positive part (flow-in) of f i (f i − is the out-flow, vice versa) and the Z switches: From ( 16) we deduce the following problem: given {w i } i and {p i } i , for a given initial conditions where Instead of the switch Z we may introduce corresponding index sets: Further, by introducing the "true mean velocities" v i = (φ i ) −1 w i , we can rewrite (18) 1 , as follows: where Dt is the material derivative w.r.t. to v i .

Transport on branching network
We consider a branching network consisting of pipes and junctions.For such structure we can derive the transport (advection) equations.Let =]x 0 , x 1 [, x 0 , x 1 ∈ R be a pipe.Thus, by x we refer to the axial coordinate along the oriented(!)pipe with the end-points x 0 , x 1 , while by X ∈ R 3 we mean the spatial positions associated with x.We consider a velocity w(x) and cross-section A(x) given at any x ∈ , which satisfy the mass conservation (By Q we denote the flux in the pipe.) Positivity of the convection velocity w is established in the context of the orientation of the pipe, i.e.T = x 1 x 0 1/w(x) dx is the transition time of the steady flow in the pipe.
It is now easy to derive the following equation for transport of the tracer, where S(x, t) is the local instantaneous saturation; possible forms of the same equation are: At the pipe ends we consider the boundary conditions: Transition times.We consider given saturations S 0 (t) and S 1 (t) at the end-points of pipe e , see (23).Eq. ( 22) 3 can be written using the material derivative as hence S(x 1 , t 1 ) = S(x 0 , t 0 ) where the transition time x 0 (w(x)) −1 dx.
Mixing and transport through junctions.We consider junctions (X j , J j ) connecting pipes { k } k∈J j , where X j ∈ R 3 is the j-th junction spatial position and J j is the set of indices of pipes connected at the junction.At any junction, a unique saturation Sj is computed using an obvious conservation law.The mean junction saturation Sj satisfies so that velocities w e in pipes e define v e depending on the oriented network topology.We can call J j + the index set of sources and J j − the index set of sinks.
Assembling equations of the transport on the network Due to (24) and knowledge of the transition times, the state of the transport is described by the junction saturations { Sj (t)} j .Equations of the model are assembled, as follows: 1.For e with a given w e define the source junction i = i(e), such that e ∈ J − , i.e. v i e < 0. It means, that S(x, t) for x ∈ e is driven by Si (t).
2. The sink junction j of the pipe e, where e ∈ J + , is associated with the general form equation ( 25), where S e (t) := Si (t − T e ), where i = i(e).
The junction equations (26) can be evaluated for discretized time interval, i.e. for t ∈ {t n } n where t n = t 0 +n∆t.Obviously, for a given T e , the saturation at t n −T e ∈ [t p , t p+1 ] is approximated using the average of values at t p and t p+1 .

Numerical examples
Preliminary testing of the above described multiscale modeling approach was performed for a simplified model of pig liver, Fig. 5 (left); domain Ω is the cube into which two perfusion trees handled by the "1D" model penetrate; they are described schematically in Fig. 5(right), whereby the lengths and average diameters of the main hepatic veins and their segments listed in Table 1 are based on a sonography examination of an experimental pig done at the Faculty of Medicine in Pilsen, Charles University in Prague.Effects of the hepatic artery were not considered in the present example.
The blood flow was computed using the algorithm described in Section 2.3.The "0D" model was treated numerically using the finite element method based on the weak formulation of problem ( 1) and ( 2) with the coupling conditions (4) related to the "1D" flow model.In Fig. 6 we display pressure fields and the associated perfusion velocities computed for the two compartments V 1 and V 2 which are coupled with the "1D" trees T P and T H , respectively.To simulate propagation of the tracer, we consider external source saturation given in the form of a time bolus defined at the input of the portal vein tree, as follows: where S = 0.4 and T = 2s.The porosity of the liver parenchym φ i , i = 1, 2 is kept constant in both systems and set equal to 0.8.The numerical solution of the tracer redistribution in the parenchyma is based on Eq. ( 18), whose spatial discretization is carried out using the cellcentered finite volume method formulated for unstructured hexahedral grids in combination with the upwind scheme.The time integration employs the two-stage Runge-Kutta scheme of second order accuracy.In Fig. 8 we display maps of the two saturations S 1 , S 2 and of the tracer concentration C evaluated using (15) at selected times.The maps are located in a central section oriented in the horizontal plane.

Conclusion
The model of perfusion proposed in the paper describes the blood flow at several scales; this modeling approach should provide a more accurate way of the perfusion quantification when compared to the classical conventional compartment models, cf.[7] The "0D" model, though describing complex flows on branching vessel networks embedded in the 3D body, is based on the Darcy flow in multiple overlapping and mutually connected compartments.The flow between them due to the pressure difference correspond to the blood filtration in the parenchyma.We illustrated a simulation of the perfusion and tracer transport on a simple example with an artificially generated geometry of the blood vessels.The tracer concentration can be compared with measurements obtained from the standard CT perfusion investigations.We plan to study the perfusion in realistic geometries.Besides various simplifications of the modeling, the most serious problem to be solved is to determine the model parameters.For this we intend to formulate and solve an inverse problem: for given maps of concentrations (from measurements), compute an optimal set of parameters which provide the best fit in terms of the computed concentrations.As the next step we intend to develop a "lumping scheme" which allows to use the homogenized model of the perfused layer to replace the ad hoc macroscopic "0D" model.The model can also be extended to capture effects of deformation.In [2] we considered an ad hoc macroscopic multicompartment model; its structure can be interpreted as a lumped version of the homogenization based perfusion models, e.g.[10,11], cf.[13], where the theory is explained in detail for a one compartment double-porosity medium.

Figure 1 .
Figure 1.A periodic structure of the layer with the highlighted reference cell comprising two channel systems A and B.
1st micro-problem.Find π b D (Y-periodic functions in y 1 , y 2 ) such that for b = 1, 2 and the two channels D ) presents the characteristic pressure response in the dual porosity to unit fluxes through the external face ∂ ± Y M : to compute γ + , unit flux on ∂ + Y M and zero flux on ∂ − Y M apply, the opposite conditions hold for γ − .Corrector functions γ +/− determine coefficients F

Figure 3 .
Figure 3.An illustration of the perfusion simulation in a layer.The last picture indicates positions of two points for which the pressure p M and perfusion velocity fields were evaluated.The arrows indicate the perfusion velocities corresponding to the displayed pressure fields.
e if e ends at junction j , v j e = −w e if e begins at junction j , J

Figure 5 .Flow in V 1 Flow in V 2 Figure 6 .
Figure 5. Simplified model of pig liver (left) and the 1D model of main hepatic vessels (right) Geometrical parameters of trees T P and T H are given in Tab. 1.
3. Let a junction j is coupled directly with junctions {k +e } e with e ∈ J j (we consider two junctions (k + e , k − e ), i.e. "source-sink" couple, for each junction.The junction equation at time t is

Table 1 .
Parameters of the 1D vessel model.For the position of the segments A -I, see Fig.5hepatic portal vein T P hepatic veins T H