Introduction

The quantitative determination of hepatic blood flow has the potential to provide important information in the assessment and follow-up of liver disorders, which are almost invariably accompanied by abnormalities in organ perfusion, representing a prognostic indicator and responding to disease amelioration [4, 10, 19, 21, 28–31, 35]. Positron emission tomography (PET) and 15O-labeled water (H2 15O) enable to assess hepatic perfusion quantitatively [29, 30, 35], as based on tracer kinetic modeling, requiring the notion of the time variation of radiotracer concentrations in the liver tissue and in the blood entering the organ (input function).

The liver is characterized by a dual blood supply, comprising the hepatic artery and the portal vein, draining venous blood from the gastrointestinal tract. Thus, in the modeling of PET data from liver, two blood time–activity curves are required to represent the input function (dual input function [DIF]). However, blood withdrawal from a peripheral artery [8, 9, 16, 17, 26, 33] is not always successful and risk-free, and it requires careful correction in time delay between the sampling site and the tissue. More importantly, the portal vein cannot be accessed from any peripheral site, making its blood collection impractical in humans.

The aim of the present study was to develop a new technique to estimate the two components of the DIF non-invasively from dynamic H2 15O PET images. The present method was characterized by use of a model input function to create a tissue model function, which was used to simultaneously fit multiple tissue curves from PET image. The parameter obtained in the input function model reproduced the input function. Computer simulation studies were performed to examine the magnitude of potential biases in the parameter estimates caused by the inherent assumptions and to optimize the extraction of multiple tissue curves from the image. The present investigation was conducted in pigs because the comparison between measured and estimated values necessitated deep catheterization and invasive Doppler flow measurements.

Materials and methods

Theory and computation of non-invasive DIF

A model function was created to shape the input function according to the dose of tracer, administration process, body weight, and physiological state in each subject [18]. The model function introduced is

$$\begin{aligned} & \begin{array}{*{20}c} {C_{{\text{A}}} {\left( t \right)} = 0.}{{\left( {t < t_{1} } \right)}} \\ { = \frac{A}{{K_{{\text{e}}} ^{2} {\left( {1 + \alpha } \right)}^{2} }}{\left( {1 - \exp {\left( {K_{{\text{e}}} {\left( {1 + \alpha } \right)}{\left( {t_{1} - t} \right)}} \right)}} \right)}}{{\left( {t_{1} \leqslant t \leqslant t_{2} } \right)}} \\ { = \frac{A}{{K_{{\text{e}}} ^{2} {\left( {1 + \alpha } \right)}^{2} }}{\left( {\exp {\left( {K_{{\text{e}}} {\left( {1 + \alpha } \right)}{\left( {t_{1} - t_{2} } \right)}} \right)} + \exp {\left( {K_{{\text{e}}} {\left( {1 + \alpha } \right)}{\left( {t_{2} - t} \right)}} \right)} - 2 \cdot \exp {\left( {K_{{\text{e}}} {\left( {1 + \alpha } \right)}{\left( {t_{1} - t} \right)}} \right)}} \right)}}{{\left( {t > t_{2} } \right)}} \\ \end{array} {\text{ }} \\ & {\text{ }}\\ & {\text{}} \\ \end{aligned} $$
(1)

Details of the model function are given in the Appendix. Briefly, A indicates the height, and t 1 and t 2 − t 1 indicate the appearance time of tracer and administration duration, respectively. K e (ml/min) and K i(= αK e) (ml/min) represent the tracer bidirectional diffusion rates between arterial blood and whole body interstitial spaces, respectively.

The portal vein blood model function was generated by introducing the gut compartment model [29–31, 35], that is, a single compartment model between arterial blood and gut compartment, assuming no difference in appearance time between arterial and portal blood (or delay time of portal input), with diffusion rate k g in the gut system as

$$C_{{\text{P}}} {\left( t \right)} = k_{{\text{g}}} C_{{\text{A}}} {\left( t \right)} \otimes e^{{ - k_{{\text{g}}} \cdot t}}$$
(2)

Using these arterial and portal input model functions, the tissue response function can be expressed by assuming a single tissue compartment model [29–30, 35] and that tracers in arterial and portal blood were well mixed before exchange with liver tissue as

$$C_{{{\text{TIS}}}} {\left( t \right)} = {\left( {f_{a} C_{{\text{A}}} {\left( t \right)} + f_{{\text{p}}} C_{{\text{P}}} {\left( t \right)}} \right)} \otimes e^{{ - k_{2} \cdot t}}$$
(3)

where k 2 is defined as (f a + f p)/V L, and V L (ml/g) is the distribution volume of water between blood and tissue. In the present study, V L was fixed to 0.7 ml/g, which was suggested to fix in a sensitivity analysis by Ziegler et al. [35] and was obtained as 0.71 ± 0.03 ml/g for same subjects in our preliminary evaluation using measured blood input functions. Including a blood volume term into this equation, the model function for liver tracer concentrations, as measured by PET (C PET), can be expressed as

$$C_{{{\text{PET}}}} {\left( t \right)} = {\left( {1 - V_{0} } \right)}{\left( {f_{{\text{a}}} C_{{\text{A}}} {\left( t \right)} + f_{{\text{p}}} C_{{\text{P}}} {\left( t \right)}} \right)} \otimes e^{{ - k_{2} \cdot t}} + V_{0} C_{{{\text{input}}}} {\left( t \right)}$$
(4)

where C input(t) is defined as

$$C_{{{\text{input}}}} {\left( t \right)} = r_{{\text{a}}} C_{{\text{A}}} {\left( t \right)} + r_{{\text{p}}} C_{{\text{P}}} {\left( t \right)}$$
(5)

where r a and r p are arterial (f a ml/min/g) and portal vein blood flow (f p ml/min/g) ratios to total hepatic flow, i.e., r a  = f a/(f a + f p) and r p = f p/(f a + f p) [35]. The flow chart to estimate input functions in this procedure is simplified in Fig. 1. Multiple tissue time–activity curves (TAC) from liver image were used to estimate the input functions. First, the model function in Eq. 4 was individually fitted to tissue TACs, assuming that k g in Eq. 2 is constant by a non-linear fitting method (variable-metric method in the PAW environment: version 2.13/08 [http://wwwasd.web.cern.ch/wwwasd/paw/]), and the set of seven parameters of A, t 1, t 2, K e(1 + α), f a, f p, and V 0 in Eqs. 1 and 4 was obtained for each tissue TAC. Then, means and standard deviations of t 1, t 2, and r a(=f a/(f a + f p)) were calculated, and the tissue TACs with values of t 1 or t 2 > 1 standard deviation of respective means were excluded to avoid the potential influence of TACs outside the liver. In the second step, assuming that all parts of the liver share the same input functions, values of t 1, t 2, and r a were fixed to their means, and the other two parameters (A and K e(1 + α)) were estimated by minimizing the following equation:

$$S^{2} = {\sum\limits_{i \in Tis} {{\sum\limits_k {{\left( {C^{{i,k}}_{{{\text{PET}}}} - {\left( {{\left( {1 - V^{i}_{0} } \right)}{\left( {f^{i}_{{\text{a}}} C_{{\text{A}}} {\left( t \right)} + f^{i}_{{\text{p}}} C_{{\text{P}}} {\left( t \right)}} \right)} \otimes e^{{ - k^{i}_{2} \cdot t}} + V^{i}_{0} C_{{{\text{input}}}} {\left( t \right)}} \right)}} \right)}^{2} } }} }$$
(6)

where C PET i ,k is the tissue TAC for kth frame in ith tissue region of interest, t is the corresponding time of kth frame, and f a i, f p i( = f a(1 − r a)/r a) and V 0 i are values of arterial and portal vein blood flows and of blood volume for ith tissue, respectively. In this procedure, S 2 was minimized by the grid search method to avoid dependency of initial guess, where S 2 was calculated for 1,000 discrete values of both A and K e(1 + α) between ranges of three standard deviations from respective mean values, omitting the negative value. In this procedure, for a given input function, i.e., given A and K e(1 + α), f a and V 0 for each TAC were computed by the grid search method, with acceptable ranges of 0–100 ml/min/g and 0–1 ml/ml, and steps of 1 ml/min/g and 0.01 ml/ml, respectively, and then substituted in Eq. 6. Finally, the image-based input function was obtained by substituting the estimated parameters into Eq. 1.

Fig. 1
figure 1

A schematic diagram of the procedure to estimate the input functions using multiple tissue TACs. Step 1 The model function (Eq. 4) was individually fitted to N tissue time–activity curves (TAC). Then, means and standard deviations of t 1, t 2, and r a were calculated, and the tissue TACs with values of t 1 or t 2 > 1 standard deviation of respective means were excluded (indicated as N′ TACs) to avoid the potential influence of TACs outside the liver. In the second step, assuming that all parts of the liver share the same input functions, values of t 1, t 2, and r a were fixed to their means, and the other two parameters (A and K e(1 + α)) were estimated by minimizing Eq. 6 by the grid search method. Finally, the image-based input function was obtained by substituting the estimated parameters into Eq. 1

Full size image

Simulation study

The present method for generating portal vein input assumes that the diffusion rate in the gut system, k g, is a fix constant, and there is no time delay between portal and arterial blood. It is not a priori known how these assumed factors degrade the accuracy of estimated DIF and flow. Moreover, tissue TACs from PET images convey some degree of noise, and the accuracy of the estimated input function might depend on either the degree of noise, or the applied number of tissue TACs, or both. A simulation study was designed to reveal the influence of the above elements on the accuracy of the current method.

To this purpose, we selected one arterial curve from one of the present experiments. First, a portal input curve was created by assuming k g = 0.5/min, corresponding to the estimated mean in all animals. The combination of these arterial and portal vein curves was treated as the ‘true DIF’. In the present experimental study, the average of activity concentrations in an area of the summed image was distributed with a 20% range around the mean for the whole liver, and this percentage was independent of the size of the selected areas in regions >50 pixels. This supports the assumption that flow values in the liver distribute around a 20% range around a mean of arterial flow of 15 ml/min/100 g [22]. Thus, by assuming ten values of f a as 13, 13.5, 14, 14.5, 15, 15, 15.5, 16, 16.5, and 17 ml/min/100 g, and ratio r( = f p/f a) = 6 [22], one set of ten hepatic tissue TACs was generated from the true DIF using Eq. 6.

The propagation of an error in k g and delay time to blood flow estimation was simulated. The sequence of steps in this procedure is simplified in Fig. 2a and b. For k g, simulated portal input curves were created from the selected arterial curve by changing k g from 0.35 to 0.65/min for error simulation in k g, and combinations of the arterial and simulated portal vein curves were used as the simulated DIF. Sets of tissue TACs were generated from these simulated DIFs, with the same assumptions of f a and r as given above. In turn, each set for each k g was used to back-estimate DIF (arterial and portal components), fixing k g as 0.5/min in this process as presented above. Finally, f a and f p were calculated from estimated DIFs for each k g by the Gauss–Newton non-linear fitting method in the interactive modeling and data analysis system called PyBLD (http://homepage2.nifty.com/peco/pybld/pybld.html) [5] using Eq. 4. For delay time, simulated portal input curves were created from the selected arterial curve by shifting the time from 0 to 10 s, and combinations of the arterial and simulated portal vein curves were used as the simulated DIF. Sets of tissue TACs were generated as above. In turn, each set for each delay time was used to back-estimate DIF fixing time delay as 0.0 s. Finally, f a and f p were calculated from estimated DIFs for each delay time. Mean of percent difference between computed and assumed (‘true’) flow values are presented as a function of k g and delay time.

The influence of noise versus number of TACs on the accuracy of the method was explored. As shown by Edward et al. [7], as the noise on tissue TACs increased, the standard deviation of uptake ratio of tracer increased; as more regions were used, the standard deviation tended to decrease. However, if the number of TACs is larger, the noise on tissue is also large and vice versa. Our simulation was intended to reveal an optimal number of tissue TACs to be extracted from the whole region of the liver. The procedure is summarized in Fig. 3. First, tissue TACs with noise were generated as follows: Gaussian noise at peak was imposed on the set of ten hepatic tissue TACs generated above. Two levels of noise were introduced, corresponding to 10% and 20% of counts at the level of the peak and 10% and 20% each of the square root of counts at the other points. This procedure was repeated 100 times and 100 sets of noisy tissue TACs, embracing a total of 1,000 pixels obtained. Next, the ith set of tissue TACs in kth frame with f a defined as \(C_{{\text{fa}}} ^{i,k} \) were summed for same f a as

$$\overline{C} ^{k}_{{{\text{fa}}}} = \frac{1}{{N_{{\text{T}}} }}{\sum\limits_i^{N_{{\text{T}}} } {C^{{i,k}}_{{{\text{fa}}}} } }$$
(7)

where N T indicates the summed number of tissue TACs and corresponds to the summed number of pixels. N T were set to 5, 10, 20, 50, 100, and 200, corresponding to a number of tissue TACs (N tis) of 200, 100, 50, 20, 10, and 5, respectively. Here, when N T was 200, the 100 tissue TACs were summed as N T = 100 and additionally combinations of f a = 13 and 13.5, 14 and 14.5, 15 and 15, 15.5 and 16, and 16.5 and 17 ml/min/100 g were summed. For each N tis and each level of noise, DIF was estimated, as described. Then, arterial and portal vein blood flow values were computed as above, using estimated DIF and tissue TACs with f a of 15 ml/100 g/min. This procedure was repeated 100 times, and the bias and deviation in values of arterial and portal vein flow results were calculated. Their bias and deviation was presented as a function of N tis.

Fig. 2
figure 2

Schematic diagram of the procedure to analyze error sensitivity in hepatic arterial (f a) and portal flow (f p) values against assumed k g (a) and time delay (b). Portal input curves were created by changing the value of k g from 0.35 to 0.65/min in (a) and by shifting the time from 0 to 10 s in (b), respectively, and combinations of the arterial (C A ) and simulated portal (C P ) curves were used as the simulated dual input functions (DIF). Sets of tissue time–activity curves (TAC) were generated from these simulated DIFs by assuming ten values of f a from 13 to 17 ml/min/100 g. In turn, each set of tissue TACs was used to back-estimate DIF fixing k g as 0.5/min and time delay as 0.0 s. Finally, f a and f p were calculated from estimated DIFs for each k g and delay time

Full size image
Fig. 3
figure 3

Schematic diagram of the procedure to analyze statistical accuracy of hepatic arterial (f a) and portal flow (f p) values against noise on tissue curves. First, tissue time–activity curves (TAC) with noise were generated by imposing Gaussian noise on the set of ten hepatic tissue TACs. This procedure was repeated 100 times, and 100 sets of noisy tissue TACs were obtained. Next, the N T (=5, 10, 20, 50, 100, and 200) sets of tissue TACs with the same flow value were summed. For each N T, dual input function (DIF) was estimated. Then, arterial (f a ) and portal blood flow (f b ) values were computed using estimated DIF and tissue TACs. This procedure was repeated 100 times

Full size image

Experimental study

PET experiment

Fourteen pigs under anesthesia with weight 30.0 ± 1.1 kg were studied. Data on glucose metabolism in these animals have been previously reported [13, 14]. Animals were deprived of food on the day prior to the study at 5:00 pm. Anesthesia was induced with ketamine (1.0 g) into neck muscles and maintained by ketamine and pancuronium (total of 1.5 g and 40 mg, respectively) administered intravenously during the experiment. Animals were intubated through a tracheostomy, and their respiration was controlled by a ventilator providing oxygen and normal room air (regulated ventilation, 16 breaths per minute). Catheters were inserted into the carotid artery for arterial blood sampling and the femoral vein for administration of H2 15O. Splanchnic vessels were accessed by sub-costal incision; after dissection of the hepato-gastric ligament, purse string sutures were allocated to allow catheter insertion via a small incision in the portal vein. A catheter was inserted directly in the portal vein for portal vein blood sampling. Ultrasound-based flow-probes (Medi-Stim Butterfly Flowmeter, Medi-Stim AS) were placed around the portal vein and hepatic artery to determine blood velocity in each vessel. The diameter of the hepatic artery and portal vein were measured off-line from B-mode ultrasound images acquired using an Acuson Sequoia 512 mainframe with a 13-MHz B-mode linear array transducer. The area of the vessel was calculated assuming circular shape. Then, blood flow was obtained for each vessel during the PET scans. The surgical access was closed, and the distal catheter extremities were secured to the abdominal surface to avoid tip displacement. The animals were then transported to the PET center for tracer administration, liver imaging, and blood sampling. Vital signs, blood pressure, and heart rate were monitored throughout the study.

PET acquisition was carried out in 2D mode using an ECAT 931-08/12 scanner (CTI Inc, Knoxville, TN, USA) with a 10.5-cm axial field of view and a resolution of 6.7 mm (axial) × 6.5 mm (in-plane) full width at half maximum. After transmission scan for attenuation correction, the dynamic scan was started after the injection of H2 15O (274 MBq, 30-s bolus injection), consisting of 20 frames with gradually increasing individual durations (6 × 5, 6 × 15, and 8 × 30 s).

During PET scanning, blood was withdrawn continuously from the carotid artery and portal vein through catheters (1.4 mm in inner diameter; length of tube was 900 mm to the detector and 60 mm in the detector sensitive region) by using a peristaltic pump (Scanditronix, Uppsala, Sweden) with a withdraw speed of 6 ml/min. Radioactivity concentrations in blood were measured with a BGO coincidence monitor system. The detectors had been cross-calibrated to the PET scanner via ion chamber [26].

At the end of the experimental period, animals were sacrificed by potassium chloride injection and anesthetic overdose, the abdominal cavity was rapidly accessed, and the whole liver was explanted and weighed and its volume was measured by water displacement; liver density was calculated as the ratio of organ weight-to-volume to derive the ultrasound-based flow to PET-equivalent unit (i.e., flow per unit of tissue volume).

The protocol was reviewed and approved by the Ethical Committee for Animal Experiments of the University of Turku.

Data processing

Dynamic sinogram data were corrected for dead time in each frame in addition to detector normalization. Tomographic images were reconstructed from corrected sinogram data by the median root prior reconstruction algorithm with 150 iterations and Bayesian coefficient of 0.3 [1]. Attenuation correction was applied with transmission data. A reconstructed image had 128 × 128 × 15 matrix size with a pixel size of 2.4 mm × 2.4 mm and 6.7 mm with 20 frames.

Measured arterial and portal vein blood TACs were corrected for physical decay and dispersion [11] as τ = 2.5 s, which was experimentally obtained and usually applied in our center. The arterial TAC corrected for decay and dispersion was then corrected for delay by fitting to a whole-liver tissue TAC [12]. The arterial curve obtained, C a(t), was used as the measured arterial input function. Then, the portal vein curve, corrected for dispersion (τ = 2.5 s) and delay with the same delay time for arterial TAC, C P(t), was fitted according to the following equation:

$$C_{\text{P}} \left( t \right) = k_g C_{\text{a}} \left( {t + \Delta t_{\text{p}} } \right) \otimes e^{ - k_g \cdot \left( {t + \Delta t_{\text{p}} } \right)} $$
(8)

to obtain k g and to account for the appearance time (Δt p seconds) via the gut system. Obtained measured curves were directly fitted with Eqs. 1 and 2 to examine adequacy for a usage of model functions.

A region of interest (ROI) was placed on the whole region of the liver in a summed image and subsequently divided plane-by-plane into sub-regions of 700 pixels each, corresponding to 11–22 sub-regions. Sub-regions were created by extracting pixels firstly from horizontal then vertical directions inside the whole ROI in each slice. Each sub-region consisted of a single area with the same number of pixels. Tissue TACs in the sub-regions were extracted from dynamic images. Then, DIF was estimated according to the procedure introduced above. In the first step, initial values and boundary conditions for the non-linear fitting (PAW environment) for each parameter were 20,000 between 0.0000002 and 200,000,000 Bq/ml for A, 5 between 2 and 20 ml/min for K e(1 + α), 1 between −10 and 100 s for t 1, 20 between 1 to 60 s for t 2 − t 1, 20 between 1 and 100 ml/min/g for f a, 100 between 1 and 400 ml/min/g for f b, and 0.05 between 0 and 1 ml/ml for V 0. In the second step, S 2 value in Eq. 6 was minimized, and the image-based input function was obtained. Areas under the curves (AUC) for measured and image-based inputs were calculated for 0 to 180 s. Their percent difference was calculated.

Perfusion values f a and f p were calculated by non-linear Gauss–Newton fitting method (PyBLD environment). Results obtained with the new technique were compared with (a) those obtained with the measured input function and (b) the ones from our independent reference method, i.e., ultrasonography, after their normalization to the organ volume to derive PET-equivalent units.

Statistical analysis

Data are shown individually or as mean ± SD. The Student’s paired t test was used for intra-individual comparisons of flow values. Regression analyses were performed according to standard techniques. A p < 0.05 was considered to be significant. Differences between the flow values were calculated as (f X  − f Y )/f Y , where f X and f Y are flow values from the non-invasive method and from the measured input or ultrasonography, respectively, and plotted in Bland–Altman plot [3].

Results

Simulation study

The biases in values of arterial, portal vein, and total blood perfusion due to a fixed k g and delay time are presented in Fig. 4a and b as a function of the value of k g and delay time, respectively. The error in total flow results did not exceed 10% for a ≤20% (i.e., 0.4–0.6 min−1) difference between the fixed and the assumed (true) k g and for a <10-s time delay.

Fig. 4
figure 4

Error in values of arterial (f a ), portal vein (f p ), and total (f a +f p ) blood flow propagated from error in k g (a) and delay time (b)

Full size image

The influence of noise and number of tissue TACs, i.e., the bias and deviation on both arterial and portal blood flow values, showed to be minimal for a number of tissue TACs of 10 to 20 at both noise levels (Fig. 5). As shown in Fig. 5, if the number of tissue TACs is increased, noise on each curve for input estimation becomes larger. On the other hand, a smaller number of tissue TAC corresponds to less information from tissue TAC in terms of variation of flow values. This result suggested that the optimal number of tissue TACs to be applied to preserve accuracy is in the above range, which is independent of the two noise levels. Among the five parameter composing the model input functions, the three parameters t 1, t 2, and r a were determined with same accuracy, i.e., both the difference and deviation in those values were less than 1 s for t 1 and t 2 and 5% for r a, respectively, for the noise level of 10%, independent of the number of tissue TACs. Bias and deviation of the remaining two parameters A and K e(1 + α) depended on the number of tissue TACs following the same tendency as the bias and deviation on blood flow values, as described above.

Fig. 5
figure 5

Bias (left) and deviation (right) in the arterial and portal vein blood flow values as a function of the number of time–activity curves applied to the estimation of the input function

Full size image

Experimental study

Reconstructed images are shown in Fig. 6, together with divided sub-regions. In the first step of our procedure, the obtained value of t 1 from TACs extracted from sub-regions overlapping the vena cava (e.g., lower sub-region at upper left side image in Fig. 6) was 10 to 13 s earlier than the mean, and these TACs were omitted from further processing. The estimated t 2 − t 1 was 27 ± 3 s, which was similar to the tracer administration duration.

Fig. 6
figure 6

View of liver H2 15O PET images in four slices and sub-regions (solid line). The small area with high activity levels on the mid-right and mid-left side of the image corresponds to the vena cava and aorta, respectively

Full size image

Figure 7 shows the curves of the model arterial and portal input functions (Eqs. 1 and 2) directly fitted to measured curves. The model functions for those were superimposable to measured curves, although both modeled curves slightly overestimated at the late times. This result suggested that the model function was almost adequate to use for the estimation of input.

Fig. 7
figure 7

Time–activity curves representing the arterial (broken line) and portal (solid line) model input functions (Eqs. 1 and 2) in comparison with the measured arterial (black circles) and portal (open circles) input functions

Full size image

The mean ± SD of k g was 0.497 ± 0.153 ml/min/g and that of Δt p was 0.7 ± 5.1 s obtained by fitting the portal TAC using arterial TAC by Eq. 8.

Estimated, image-derived arterial and hepatic input functions were almost superimposable to the measured curves (Fig. 8). The mean ± SD and range of difference of AUCs were −3.15 ± 8.73% ranging from −13.5% to 17.9% and 1.47 ± 8.87% ranging from −13.5% to 10.2% for arterial and portal input functions, respectively. The coefficient of variation of the estimated flow ratio between artery and portal vein in the first step across sub-regions was 26 ± 9%. The mean ± SD of that ratio across subjects was 0.15 ± 0.07 and those from ultrasonography was 0.16 ± 0.06, and paired t test showed no significant difference between them. This suggests supporting the assumption that the ratio between arterial and portal input defined in Eq. 5 relates to the flow values.

Fig. 8
figure 8

Estimated arterial (red line) and portal vein (blue line) input functions from PET images and their comparison with measured arterial (plot in light blue) and portal input (plot in pink) functions

Full size image

The Bland–Altman plot between values of hepatic arterial, portal, and total perfusion, as estimated by using the image-derived versus the measured blood curves, is shown in Fig. 9. This plot demonstrates a small overestimation by image-derived method with a bias of 0.01 and 0.07 ml/min/g for arterial and portal flow, respectively, and that 0.08 ml/min/g for total flow. Respective regression lines were the following: y = 0.00 + 1.09 x (r = 0.97, p < 0.001), y = 0.05 + 1.02 x (r = 0.87, p < 0.001), and y = 0.02 + 1.06 x (r = 0.90, p < 0.001). Paired t test showed no significant difference between the methods. Differences were −6.8 ± 20.0%, −4.9 ± 14.3%, and −5.8 ± 15.6% for arterial, portal, and total blood flow values, respectively.

Fig. 9
figure 9

a Bland–Altman plot for arterial (square), portal (triangle) and b total hepatic blood flow differences between measured and image-derived input functions

Full size image

The Bland–Altman plot between values of hepatic arterial, portal, and total perfusion, as estimated by using the current method versus ultrasonography, is given in Fig. 10. This plot demonstrates an overestimation by image-derived method with a bias of 0.02 and 0.22 ml/min/g for arterial and portal flow, respectively, and that 0.24 ml/min/g for total flow. Respective regression lines were the following: y = 0.06 + 0.69 x (r = 0.69, p = 0.12), y = 0.41 + 0.98 x (r = 0.54, p = 0.025), and y = 0.24 + 0.97 x (r = 0.60, p = 0.022). Again, paired t test showed no significant difference between the methods. Differences were 3.6 ± 52.0%, 15.5 ± 31.3%, and 16.9 ± 33.0% for values of arterial, portal, and total blood flow, respectively.

Fig. 10
figure 10

a Bland–Altman plot for arterial (square), portal (triangle), and b total hepatic blood flow differences between ultrasonography and kinetic modeling using image-derived input functions

Full size image

The total flow values ranged from 0.5 to 2 ml/min/g in the animals (Figs. 9 and 10). However, only two out of 14 showed smaller values of 0.5 ml/min/g (i.e., approximately 500 ml in the whole organ), which is still physiologically reasonable, while the great majority clustered between 1 and 2 ml/min/g.

Discussion

In the current work, we developed and validated a method to estimate the two components of the hepatic dual input function from liver H2 15O PET images and quantify hepatic perfusion. Computer simulations were used to evaluate the influence of assumptions, noise in raw data, and number and size of the regions of interest to be used in the analysis. After demonstrating that k g can be assumed within a 20% range by introducing a negligible error in perfusion estimates and that 10–20 regional time–activity curves appear optimal, the method was validated experimentally by showing its coherence with measured blood tracer levels and with liver perfusion results obtained by an independent technique.

The current approach estimated the hepatic arterial and portal input functions from multiple tissue curves to calculate respective and total organ perfusion. A high degree of overlap and tight correlations were observed between the estimated input functions and those obtained during online blood sampling/counting. Consequently, calculated flow values were consistent between the methods. Alternative to the present procedure, a ROI-based input extraction from PET images has been used for the carotid artery in [11C]flumazenil brain studies [27], abdominal artery for kidney blood flow quantification with H2 15O [15], and aorta for cardiac 18F-FDG metabolism [32], and for tumor blood using H2 15O [34]. In these approaches, ROIs are drawn in visible vessels, and the partial volume effect must be taken into account by testing different ROI sizes or by thresholding the pixels inside arterial ROIs; the need for partial volume correction remains a necessary limitation. Closer to the current analysis, Edward et al. applied multiple tissue curves to estimate quantitative kinetic parameters in the brain [7] and well reproduced the input function for H15 2O. However, their formula did not take into account the radioactivity from the blood component inside the tissue ROI, and the validity of their method may not be directly extrapolated to the liver because of the large proportion of blood, which is typically ranging between 0.27 and 0.40 ml/ml [23] in this organ. The present method showed that the height of the estimated input is almost doubled if the blood volume is not included in the formula and if the arterial volume contributes 10% of radioactivity in the tissue TACs in our preliminary study (data not shown). The shape of an arterial input function from multiple tissue TACs has been well reproduced in brain 18FDG or [11C]MPDX studies by using an independent component analysis-based method (extraction of the plasma TAC using independent component analysis [24, 25] still requiring one arterial blood sample, and the combination of the latter and the current techniques may be of further simplification and deserves investigation since it would entail neither a model function nor direct blood measurements.

One advantage of introducing a model function was to shape the model input function by imposing constrains to the parameters range, allowing to overcome noise problems caused by limited scan duration and short half-life of 15O. The present approach may be applicable to a study group including subjects with hepatic disorders as far as measurement conditions are equivalent and the shape of the input function can be expected to be similar, though the validity of the present method was tested in normal animals. We expect no relevant limitation in the extension of the assumptions concerning the shape to other species and in a majority of hepatic conditions. A drawback in the use of a model function, however, is that the feasibility is unknown for a group in which the shape of input functions could be extremely different or cannot be expressed by the present model function. This is not a commonly expected case. In this situation, the present method would require to, and may still be adapted, by introducing group-specific parameter constraints or a modified model function. The present model function was created by assuming the model, namely, tracer bidirectional diffusion to whole body as in differential Eq. 9. The solution was derived as Eq. 10, and the model function was modified to avoid the one order term of t, which would complicate calculations in the following procedures, i.e., model function for portal input and for tissue response functions. This modification could deteriorate the physiological mean of parameters, such as K e and K i; however, the input functions obtained in the present study using this modified equation well reproduced the shape of measured inputs. The modified model function and derived portal model function seemed to be superimposable to measured blood TACs, although there were slight, few-second systematic misalignments in the peak of arterial blood and overestimations at the late phase. This suggested that the error in the position of the peak and in late phase in the estimated input function against measured ones (Figs. 7 and 8) is due to a limitation in the description of the model function.

The present estimation procedure followed two steps, as designed to fit tissue curves individually, and then simultaneously. The first step allowed careful exclusion of tissue TACs showing t 1 or t 2 values over one standard deviation from the mean to eliminate the influence of radioactivity outside the liver region. In fact, in the experimental procedure, H2 15O was injected in the femoral vein, draining into the vena cava, which is not distant from the liver. Other adjacent high-perfusion organs include the kidneys. The influence of ROIs drawn in proximity of these regions was not included in the model. Thus, special attention was paid at excluding confounding tissue TACs by examining t 1 and t 2. In the above examples, the tracer was expected to show an early peak in case of an anatomical overlap with the vena cava, and the extracted TAC covering this region was omitted in this step. The second step was introduced to facilitate the achievement of the convergence by fixing the values of t 1, t 2, and r a to their calculated means (as obtained above) to estimate the remaining two parameters. Generally, if many parameters are estimated in a fitting procedure such as in the present method, there could be many local minima, and uniqueness of parameter solution might not be guaranteed. As shown in the simulation study, the three parameters t 1, t 2, and r a were estimated independent of the number of tissue TACs; however, the remaining two parameters A and K e(1 + α) were dependent on that. This suggests that correlation among parameters due to their numerosity could not be prevented. However, the shapes of input functions were reproduced, and flow values were consistent with other two methods, i.e., those computed from measured inputs and from ultrasonography. Thus, the correlation among parameters did not seem to affect the estimation of flow values, although further study is required for optimization.

We used a fixed value of k g to represent the diffusion rate of water between arterial blood and the gut compartment in the estimation of the portal input. The deviation in this rate constant was about 26% in the current study group. The simulation analysis showed that values within 20% of the assumed true k g number corresponded to a propagated error of 10% in the final estimation of hepatic perfusion. The value of k g used in our final computations is in accordance with the recently reviewed concept that [20] in mammals, the general biological rate (uptake ratio) varies approximately in proportion to the 3/4 power of body size and, given a body mass of ∼60 kg, k g, which is the uptake rate of water in the gut system, can be predicted to fall around 0.45 min−1 in humans. This number is consistent with a mean figure of 0.5 min−1, as obtained in this study, suggesting that the present assumption could be implemented to obtain the liver input function in humans. We also assumed a time delay of portal input to be zero against the arterial one. The deviation in this time was about 0.7 s in the current study group. The simulation analysis showed that an error in this value within 10 s corresponded to a propagated error of less than 10% in the final estimation of hepatic perfusion. Of further strength, a close agreement was shown between estimated and measured blood activity curves and estimated and Doppler-determined liver flow results. The larger difference of the latter result against the former result might be due to the model assumptions in flow calculations, as well as in the assumption of circular shape when estimating the area of the arterial and portal vessels by ultrasonography and in the accuracy of ultrasonography data (from multiple measurements of flow data, coefficient of variation was 13 ± 5% for portal flow and 18 ± 10% for hepatic arterial flow with this study [data not shown]). In this study, the flow values were calculated assuming the dual input, single compartment model [2, 29, 30, 35]. Altogether, the above observations support the use of a fixed k g and the current model in the fully non-invasive quantification of liver perfusion.

The validation of the current approach, as obtained in this study, is especially valuable in the liver for multiple reasons. First, the inaccessibility of the portal vein prevents its direct blood sampling in humans. Arterial blood can be obtained [8, 9, 16, 17, 26, 33], but blood counting requires corrections for dispersion, delay between target organ and sampling device, and cross calibration between PET scanner and radioactivity counter, which are all potential sources of errors, in the same magnitude as that expected with the current method. Second, liver perfusion can be compromised both as consequence and cause of hepatic disease and is considered a prognostic indicator and useful marker during progression or treatment follow-up [6, 22]. Third, the possibility to distinctly quantify portal and arterial perfusion is important because their reciprocal compensation may be masked once only if total hepatic flow is measured.

The present simulation study allowed to establish that the optimal number of tissue TACs for DIF estimation was 10 to 20, independent of the noise levels, among the ones selected in this investigation. As pointed out by Edward et al. [7], as the noise on tissue TACs increased, the standard deviation of uptake ratio of tracer increased. Also, they suggested that the standard deviation tended to decrease when more regions were used. The present study intended to investigate the optimal number of tissue TACs from the whole region of liver. The noise in the liver can be minimized by placing a ROI to cover the whole organ and subsequently dividing it in a number of sub-regions corresponding to 10–20 under the conditions of the current experiments. The present results may depend on the reconstruction method. However, as far as the PET image is calculated quantitatively and the distribution of flow values in the extracted TACs is in the same order of magnitude as the present study, the results of optimization in this study would be applicable because those two conditions were assumed in the present simulation study. We assumed that the ratio of blood flow between the hepatic artery and the portal vein was uniform in the whole organ, as supported by an extended literature on the healthy liver and on a majority of metabolic disorders involving the organ. Conversely, the quantification of flow in hepatic tumors in which perfusion from arterial blood is predominant may be best approximated by simplifying the procedure to a single input or by fitting the relative vascular (arterial and portal) contributions as additional parameters in the model. The current procedure was validated for the determination of liver perfusion with H2 15O PET data. Required conditions were a model function to describe the input function and a kinetic model for tracer exchange between blood and tissue. In theory, the present method might be adapted to other tracers and organs if tracer kinetics in the tissue can be described with a model function.

In conclusion, our results demonstrate that arterial and portal vein concentrations of labeled water can be estimated directly from tissue time–activity curves obtained through dynamic H2 15O PET imaging. The calculated hepatic arterial, portal, and total perfusion values using estimated or measured input functions were similar and consistent with ultrasonography measurements.