Single-voxel delay map from long-axial field-of-view PET scans

1 Introduction

Dynamic positron emission tomography (PET) scans are widely used as a non-invasive technique to estimate physiological parameters of different tissues, such as the blood perfusion of organs using [¹⁵O]H₂O (1–3) or the metabolism of glucose in the brain or tumors using [¹⁸F]FDG (4). Kinetic parameters of interest can be calculated by modeling the tissue time activity curve (TAC) as a response to an input curve or input function, typically employed in the form of compartmental models (5, 6).

Prior to kinetic modeling and particularly for the short-lived isotopes it is essential to correct for relative transport delay between the input function, i.e., the tracer concentration in the blood as a function of time, and the tissue response curve, which are typically acquired from two different anatomical locations in the body. On a clinical standard size PET scanner of approximately 25 cm and if the organ of interest cannot be in the same field of view as a large arterial space, the input function is usually measured by sampling blood from an arterial cannulation—also called an arterial input function (AIF)—while on whole-body scanners the input function is often derived from a segmentation of the aorta—called an image-derived input function (IDIF)—that does not require any invasive measurements (7, 8).

Long-axial field-of-view (LAFOV) PET scanners are capable of capturing much more of the available signal due to the extended coverage, enabling dynamic scans with faster frame rates of around one second, reducing motion artifacts and increasing the accuracy of kinetic modeling and generation of parametric images (9–11). The increased sensitivity of these scanners also allows for lowering the dose of radiotracer administered to the patient, ultimately lessening the radiation exposure experienced by the patient (12). In addition, the high spatial resolution of LAFOV scanners reduces partial volume effects and allow for a more accurate image-derived input function (7, 8, 13).

Modeling the kinetics and physiological parameters of an organ often involves using a single mean TAC for the entire organ. This approach improves the signal-to-noise ratio and the robustness of the subsequent fit. However, this method has a drawback as it potentially yields less accurate regional modeling results by neglecting the heterogeneity within the organ (14). To address this issue and analyze the variability of kinetic parameters across different regions of the organ, single-voxel modeling can be employed. Despite its advantages, single-voxel modeling faces significant challenges, primarily due to the low signal-to-noise ratio. The advancement of LAFOV scanners, which benefit from an axial FOV of 1 m or more, has mitigated this problem. These scanners offer enhanced sensitivity and improved time resolution, making the kinetic modeling of single voxels more reliable and accurate (15).

Another aspect to consider is the kinetic modeling of organs where regions differ significantly in terms of delay time between the input function and the tissue curve. For instance, in the case of the brain in patients with single-sided carotid stenosis, employing an average delay based on the mean organ TAC may lead to inaccuracies, since a regional delay or even voxel-based delay, which considers the characteristics of smaller, albeit noisier, regions, might be more appropriate to capture the heterogeneity in tracer arrival times (16).

PET radiotracers are designed to target distinct biological functions, and thus their uptake, retention and clearance from tissue differ depending on the specific tissue they pass. In terms of compartmental modeling, [¹⁵O]H₂O is described well by a one-tissue compartment model as it diffuses freely between blood and tissue, approximating an extraction fraction (EF) of one, while [¹⁸F]FDG is better described by a two-tissue model as it can enter into a metabolized state, and in addition, other tracers can also exhibit a number of other properties, which in turn alter the specific shape of the tissue TAC and the exact model used (17).

The conventional approach of delay correcting mean tissue TACs to an input function may therefore not be adequate due to both the noisy nature of single voxels with high time resolution as well as the different characteristics of radiotracers. In this article we explore a new method of estimating the delay of these tissue TACs which involves fitting the cumulatively summed TAC to increase the robustness of the fit. Additionally, we examine the performance of this algorithm across a range of tracers, each with distinct noise characteristics.

To assess the performance, we will compare our proposed method (18) to two commonly applied methods, namely cross-correlation and a one-tissue compartment model including delay (19). Herein, we perform an evaluation on both synthetically produced TACs across a range of delays and levels of noise. Whole body delay maps are calculated for each method for three tracers covering a large range of EFs.

2 Methods

2.1 Data acquisition

Dynamic PET scans from 15 patients were used in the development and testing of the algorithm, divided equally (n = 5) among each of the tracers [¹⁵O]H₂O, [¹⁸F]FDG and [⁶⁴Cu]Cu-DOTATATE. The patients were scanned on a Siemens Biograph Vision Quadra PET/CT scanner with an axial coverage of 106 cm with a scan start approximately 10 s prior to tracer bolus injection. Data were reconstructed using a maximum ring difference (MRD) of 85 in a 440 × 440 voxels matrix resulting in a 1.65 mm in-plane voxel size. Raw list mode data were subsequently reconstructed using an ordered subset expectation maximum (OSEM, 4 it., 5 subsets) algorithm with scatter, point-spread function, time-of-flight and CT-based attenuation correction applied. Reconstruction filter used was a 3D Gaussian filter with 2.00 mm FWHM for [¹⁵O]H₂O and an all-pass for [¹⁸F]FDG and [⁶⁴Cu]Cu-DOTATATE. For [¹⁵O]H₂O, the frame durations were 40 frames × 1 s + 5 frames × 4 s + 6 frames × 10 s + 3 frames × 20 s. For [¹⁸F]FDG and [⁶⁴Cu]Cu-DOTATATE, frame durations were 40 frames × 1 s + 10 frames × 5 s + 15 frames × 10 s + 6 frames × 60 s. For analysis of TACs, the frame mid times are used.

An image-derived input function (IDIF) was derived from a TotalSegmentator (20) segmented aorta subsequently transferred to the PET scan. The entire aorta was divided into four segments: the ascending part of the aorta, the aortic arch, the proximal descending part of the aorta and the distal descending part. The proximal descending part of the aorta was used to derive the IDIF across all tracers as the average time activity curve within a 1 ml volume of interest positioned lengthwise inside to avoid partial volume effects.

The choice of tracers included ones that are often used in dynamic PET imaging to measure blood flow and rate of metabolis, [¹⁵O]H₂O and [¹⁸F]FDG, while [⁶⁴Cu]Cu-DOTATATE was included to test the performance and limitations of the proposed model on a tracer with very low signal-to-noise ratio and first pass extraction fraction.

The project was approved by the Departmental Review Board Rigshospitalet, University of Copenhagen on 17 September 2021. All patients provided written informed consent prior to inclusion.

2.2 Delay map algorithm

The algorithm estimates the onset times for the input function and a tissue time activity curve independently and subtracts the two to estimate the delay between the input and the tissue ($\mathrm{\Delta }T=T_{\mathrm{tissue}}-T_{\mathrm{input}}$). Both the input function and the tissue TAC go through the same processing as described below.

Due to the inherent noisy characteristics of single voxels, several approaches were employed to sensibly handle the noise. First the TAC was summed cumulatively. This eliminates most of the approximately Gaussian distributed noise. Secondly, plateaus were detected and removed in case of noise in the early part of the TAC before fitting the summed TAC for up to 10 different ranges to find the best matching fit within a set of criteria. In case of no acceptable fits, the average TAC of a 5 × 5 × 5 cube around the failed voxel was calculated and the algorithm was rerun on these average TACs only. This cube will have a volume of 0.56 ml for the [¹⁵O]H₂O scans and 0.68 ml for the [¹⁸F]FDG and [⁶⁴Cu]Cu-DOTATATE scans.

2.2.1 Fitting the summed TAC

For ease of fitting, a simple fitting function with few parameters is proposed, which was empirically found to fit different types of TACs well:

$\begin{array}{c}{A}_{\mathrm{tissue}}(t;a,b,c)={\displaystyle \frac{a}{1+\exp \left({\displaystyle \frac{c-t}{b}}\right)},}\end{array}(1)$

Where a is the amplitude or saturation level when t is large, i.e., $t\gg c$; b indicates how rapid the function transitions from 0 to the saturation level, a, and c is the center point around which the function transitions. Here $b=0$ means that the change happens instantly.

This function fits the TACs from zero activity at the beginning and up to a saturation level. However, as observed in Supplementary Figure S1, time activity curves for single voxels can be inherently dominated by noise, which in turn complicates the fitting process and the robustness of fitting directly to them. Thus, to reduce the effect of the noise and increase the robustness of the fit, the TACs are summed cumulatively, and fitted by the integral of $A_{\mathrm{tissue}}(t)$ (see Figure 1):

$A_{\mathrm{sum}}(t;a,b,c)={\int }_0^t{A}_{\mathrm{tissue}}(t^{\mathrm{\prime }})d{t}^{\mathrm{\prime }}=ab\ln \left[\exp \left({\displaystyle \frac{c-t}{b}}\right)+1\right]+a(t-c),(2)$

where a, b and c were defined in Equation (1). The onset time being determined from the fitting parameters as

$\begin{array}{c}T=c-3.5\cdot b.\end{array}(3)$

As c is the central point of transition from flat to a rising curve, the onset point lies before this time. To compensate, a correction term which depends on the tightness of the transition, b, is added. The constant value of 3.5 has been chosen empirically to provide a uniform appearance across the parameter maps depicted in Figures 3A,B.

Figure 1

Figure 1. Summing the TACs cumulatively eliminates most of the noise and hence increases robustness of the fit using Equation (2).

Figure 3

Figure 3. Parameter maps for the proposed model (left column), cross-correlation (middle column) and one-tissue compartment model (right column) showing the residual delay (difference between the actual delay and estimated delay) for synthetic TACs created from a one-tissue compartmental model using an IDIF of the tracer [¹⁵O]H₂O as $C_A(t)$ in Equation (5), across combinations of $K_1$ (horizontal axis) and $k_2$ (vertical axis), linearly spaced in the range $1\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }\le K_1\le 350\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ and $0.05\mathrm{mi}{\mathrm{n}}^{-1}\le k_2\le 3\mathrm{mi}{\mathrm{n}}^{-1}$. Shown here for two levels of noise (low noise: A, C, and E; and high noise: B, D, and F) and for an actual delay between input function and tissue TAC of 0 s, as the parameter maps do not change significantly between different actual delays.

While summing the TACs, the frame durations [the dt' in Equation (2)] must also be multiplied to each TAC data point to avoid abrupt changes at the points where the frame timing changes (see Supplementary Figure S2 for a constant activity curve sampled at varying time intervals with and without the correction). The summed TAC is thus:

$\begin{array}{r}TA{C}_{\mathrm{sum}}(t)=\sum _{i=0}^t(TAC(i)\cdot (\mathrm{FrameMidTime}(i+1)-\mathrm{FrameMidTime}(i))).\end{array}(4)$

2.2.2 Simple plateau detection

Figure 2 shows an example TAC with noise spikes, and how these can create steps and plateaus when the TAC is summed cumulatively. These steps bias the fitting function and might set the onset time at the first step, where the second step would be the arterial pass. To mitigate this, a simple plateau detection algorithm is utilized to detect when 3 or more data points lie within a band of 0.1% of the max value of the first two minutes of the summed TAC. The threshold is chosen to be relative to the TAC to accommodate TACs of different amplitudes and can be fine-tuned for an either more or less sensitive detection of plateaus. Only the first two minutes are checked, to prevent the tail of the TAC to falsely be detected as a plateau. When a plateau is detected, everything before that is set to zero and every following point is lowered by the value of the plateau.

Figure 2

Figure 2. Example TAC from a voxel in the arm of injection (green), cumulatively summed TAC (blue), IDIF (yellow). Possible partial volume effects in the early part of the TAC creates double spikes as the tracer passes tissue twice (venous pass and a later arterial pass). This creates plateaus in the summed TAC.

2.2.3 Parameter fitting ranges

The fitting function, Equation (2), tends to the straight-line asymptote: $A_{\mathrm{sum}}(t\gg c)\to a(t-c)$, however, the summed TACs usually fall off after a few minutes. To obtain the best possible fit to the upslope, up to 10 different ranges are fitted. The endpoints of the fitting ranges are selected based on when the summed TAC starts separating from the asymptote of the fitting equation. This is detected by computing where the absolute values of the second derivative of the summed TAC are highest. This way, both cases where either the summed TAC falls off below the fit asymptote or rises above it (for example due to a high uptake rate combined with a low clearance rate) are considered.

These points are identified through simple peak detection (21) by comparing adjacent points, and the top 10 peaks with the largest peak heights are chosen as the endpoints.

The fitting parameters were restricted to $0.1\frac{\mathrm{Bq}}{\mathrm{ml}\cdot \mathrm{s}}\le a\le 2\times {10}^7\frac{\mathrm{Bq}}{\mathrm{ml}\cdot \mathrm{s}}$; $0.01s\le b\le 10s$; and $0s\le c\le 60s$, which were set broad enough to accommodate all types of summed TACs but more importantly provides a starting point for the fitting function.

2.2.4 Filtering the fits and choosing the best ones

After fitting one range of the summed TAC, as described above, the fit is tested against a set of criteria as follows. After checking for any fitting errors returned by the fitting function, the goodness of fit (GoF: $R^2$) is calculated and evaluated against a lower limit of 0.8. Next, the onset time from Equation (3) is calculated, which must be positive and lie before the endpoint of the fitting range (as found in the previous subsection).

The top half of the fits with the highest $R^2$ are collected and the median of their fitting parameters are calculated. From these, the onset time of the TAC is calculated using Equation (3).

2.3 Validation against other methods

The proposed method is compared against simple cross-correlation and a one-tissue compartment model with incorporated delay (19), to determine the strengths and weaknesses of each. For cross-correlation, the point of highest correlation is found, and for the one-tissue compartment model, an “endtime” of 180 s was used while all other parameters were kept at their default setting [as given in (19)].

2.4 Synthetic TACs

Synthetic TACs ($C_T(t)$) were generated by applying a one-tissue compartment model,

$\begin{array}{c}{\displaystyle \frac{d{C}_T(t)}{dt}=K_1{C}_A(t)-k_2{C}_T(t),}\end{array}(5)$

to an input function ($C_A(t)$) for a range of combinations of rate constants ($K_1$ and $k_2$), delays and added Gaussian distributed noise, as well as for input functions from different tracers. The three methods described above were tested on all the synthetic TACs to evaluate how well each method were able to estimate the known delay between the input function and the synthetic TAC.

Delay between tissue TAC and input function was added by moving the tissue TAC and interpolating to the original time sampling.

Noise was added to the tissue TAC by multiplying each timepoint by a random number from a Gaussian distribution centered at 1 with a standard deviation of $\sigma$. The standard deviation is the parameter used to vary the amount of noise. Any points with a negative value after noise was added were set to 0.

The ranges for each of the parameters used are given in Table 1.

Table 1

Table 1. Ranges of parameters for synthetic tissue TACs.

The standard deviation range of the noise was chosen to match the noise level in the acquired PET data. A standard deviation of $\sigma =0.1$ approximately corresponds to the noise of a few hundred voxels averaged, while a standard deviation of $\sigma =0.4$ roughly corresponds to the noise of a single voxel. We will refer to these two noise levels as “low noise” and “high noise”.

3 Results

3.1 Synthetic TACs

Figure 3 shows the residual delay between the known delay added to the tissue TAC and the delay estimated by each model. A negative residual delay indicates that the tissue TAC is estimated to be delayed to a later time than necessary. Each combination of $K_1$ and $k_2$ corresponds to a specific tissue TAC shape as given by Equation (4), where the left side of the parameter maps represent steadily increasing TACs with low rate of clearance and high uptake, while along the bottom the TACs have low uptake and a high rate of clearance, resulting in low signal TACs. Between these two extremes is the more physiologically relevant area with $30\frac{\mathrm{ml}}{100\mathrm{g}}<\frac{K_1}{k_2}<100\frac{\mathrm{ml}}{100\mathrm{g}}$.

All models show predominantly negative residual delays over the entire parameter map, indicating that the estimated delay was greater than the actual delay. Visually, the proposed model shows an overall consistent residual delay, cross-correlation shows very negative residuals, while the one-tissue compartment model with incorporated delay shows low residual delays in the middle of the parameter map for low noise, although as the noise increases, so does the residual delay. Both the proposed model and cross-correlation shows no dependency on $K_1$.

Figure 4 shows the distributions of residual delays across the whole parameter map for increasing levels of noise.

Figure 4

Figure 4. Histograms of the parameter maps in Figure 3 across a range of noise levels, $0\le \sigma \le 0.5$, downward, for the three models tested: proposed (A), cross-correlation (B), and one-tissue compartment model (C). Note that the x-axis is bounded in a range from −8 s–+8 s, and that the residual delay of cross-correlation mainly lies below this range.

The proposed model displays a moderately consistent residual delay of $-0.75\pm 0.38\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for low noise $(\sigma =0.1)$ to $-0.85\pm 0.88\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for high noise $(\sigma =0.4)$.

Cross-correlation displays a strong negative residual delay across the whole parameter map with mean values ranging from $-12\pm 6\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for low noise $(\sigma =0.1)$ to $-13\pm 6\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for high noise $(\sigma =0.4)$. The residual delay increases rapidly as $k_2$ decreases.

One-tissue compartment model displays a wavy histogram with valleys at whole numbers and smaller distributions in between, peaking at half-integers. It has a mean of $-1.5\pm 1.5\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for low noise $(\sigma =0.1)$ to $-2.9\pm 2.1\frac{\mathrm{ml}}{100\mathrm{g}\cdot \min }$ for high noise $(\sigma =0.4)$. The one-tissue compartment model appears to perform well in a triangular region (indicated by the diagonal lines in Figure 3) bounded by the distribution volume in the range $30\frac{\mathrm{ml}}{100\mathrm{g}}<\frac{K_1}{k_2}<100\frac{\mathrm{ml}}{100\mathrm{g}}$. However, as the noise increases this area also becomes increasingly noisy.

3.2 Whole-body delay maps

Figure 5 displays whole-body delay maps for the three models across three tracers. It is observed that [¹⁵O]H₂O generally gives the least noisy delay maps, whereas [⁶⁴Cu]Cu-DOTATATE gives the most noisy.

Figure 5

Figure 5. Coronal slice of whole-body delay maps for the proposed model (left column/subfigures A, D, and G), cross-correlation (middle column/subfigures B, E, and H) and a one-tissue compartment model with delay included (right column/subfigures C, F, and I), across the tracers [¹⁵O]H₂O (top row), [¹⁸F]FDG (middle row) and [⁶⁴Cu]Cu-DOTATATE (bottom row).

Generally, across all models and tracers, negative delays are observed in the vein of injection, heart and lungs as is expected with the input function derived from the descending part of the aorta. The proposed model additionally estimates negative delays in most of the arm of injection as well as the superior liver, while one-tissue compartment model estimates negative delays in the arm of injection with good separation of the liver. Cross-correlation displays very sharp separation of organs, for example between the lung and liver as well as between the vein of inject and the rest of the arm, however, it generally seems to overestimate the delay.

For the [¹⁸F]FDG delay maps, cross-correlation estimates a mean delay across the tested patients of 100 s ± 60 s in the brain, while the other models estimate a delay of around 4 s ± 4 s.

The one-tissue compartment model displays good separation between organs, particularly the lungs and liver for [¹⁵O]H₂O. For [¹⁸F]FDG, this separation becomes less defined, while the delay in the aorta—the location where the IDIF is derived from—is estimated to be around −2 s ± 2.5 s. For [⁶⁴Cu]Cu-DOTATATE, this model estimates spots of negative delays around the body and most noticeably in the brain.

3.3 Organ delay distributions

Figure 6 shows the delay distributions of selected organs estimated by the proposed model across the three tracers. The organs have been segmented using TotalSegmentator (20).

Figure 6

Figure 6. Delay distributions of selected organ for the tracers [¹⁵O]H₂O (left/subfigure A), [¹⁸F]FDG (center/subfigure B) and [⁶⁴Cu]Cu-DOTATATE (right/subfigure C).

Visually, the distributions for [¹⁵O]H₂O are smooth, while for [¹⁸F]FDG and [⁶⁴Cu]Cu-DOTATATE each organ distribution contains several smaller peaks. These smaller peaks occur at a rate corresponding to a frequency of approximately 1 Hz.

The distributions show the temporal progression of the tracer: the lungs can be observed to be located in between the right and left ventricle after which the tracer enters the aorta and later reaches the other organs. The liver, however, spans the whole time scale, with the first half being the superior part of the liver.

3.4 Effect of voxel-wise delay correction on parametric images

The effect of delay correction can be observed from Figure 7, where subfigure (a) shows the modeled blood flow with a common delay correction applied to all voxels based on the mean time activity curve for the organ (found using the proposed model), which is a common approach to adjust for delays in tracer arrival. Subfigure (b) presents a more refined approach, applying voxel-wise delay correction based on the proposed model. This accounts for the differences in tracer arrival time between voxels and on a regional basis. The relative difference between the two approaches is depicted in subfigure (c), where a difference of around 15%–25% can be observed on the side which is most affected by carotid stenosis.

Figure 7

Figure 7. Parametric images showing the blood flow {K₁, [ml blood/(100 g tissue * min)]} of an axial slice of brain in a patient with carotid stenosis on their left side (the right side in the images) receiving the tracer [¹⁵O]H_sO. (A) is delay corrected using the mean organ time activity curve, while (B) uses voxel-wise delay correction with delay correction, and (C) shows the relative difference between the two.

4 Discussion

The negative delays observed in Figure 5 in the arm and shoulder in the tracer injection arm, as well as in the superior part of the liver for the proposed model, may be attributed to partial volume effects. A large amount of concentrated radiotracer was administered via an antecubital vein injection, which, as it travels to the heart through the venous system, affects the tissue surrounded by the high concentration bolus such as the rest of the arm and the shoulder. High levels of activity are also present in the heart, which, together with the heartbeat and breathing can contribute to too early estimates of the onset times in the superior part of the liver, and thus delay between it and the IDIF, as observed in the delay maps.

The above mentioned TACs can therefore exhibit “double-peaks”; the first for the venous bolus passage that affects it by partial volume or scatter-related effects, and the second peak for the arterial pass. The desired peak we aim to detect here is the second one. These double-peaks in the TACs can partly be resolved by searching for plateaus in the summed TACs, as described in the methods section. If there are long periods with little to no signal between two peaks, for example in the case of the arm or superior liver, the second peak should be selected. The distinction between two peaks is, however, not always visible as they can start to blur together, and sometimes the first peak might be selected in these cases.

Plateau correction had to be performed for each tracer in different amounts across the body. For H₂O, the main concern was in the arm of injection and other close-lying tissue where the initial tracer bolus would create a spike early in the time activity curves that would not correspond to the first arterial pass. For FDG or DOTATATE the main concern was additionally the high amount of noise present.

As mentioned, in case the noise is too great to produce a good fit, the average of a 5 × 5 × 5 cube of voxels would be calculated around the failed voxels and used as a new TAC to run the algorithm on once more. In patients this had to be performed for on average 0.3% of voxels for the [¹⁵O]H₂O scans, 0.03% for the [¹⁸F]FDG scans and 3% for the [⁶⁴Cu]Cu-DOTATATE scans. This would usually be in the extremities of the patient, such as the legs, hands, or extremities in the lower torso. The highest number of failed voxels happened in the DOTATATE scans as expected due to its low first pass extraction and thus low signal-to-noise ratio, however the FDG scans have a lower error rate than the H₂O scans, which could be due to the slightly slower kinetic profile of [¹⁸F]FDG over [¹⁵O]H₂O together with the smaller mean free range of the ${\beta }^{+}$ from ¹⁸F compared to ¹⁸O, giving rise to less partial volume effects and an overall more even TAC.

The one-tissue compartment model with incorporated delay works by fitting a one-tissue compartment model to the tissue TAC with included delay. The tracer kinetics of [¹⁵O]H₂O can likewise be described very well by a one-tissue compartment model, and as such, it is expected to perform well when estimating the delay for this particular tracer. This is especially true for high signal-to-noise ratio TACs such as mean organ TACs, but as the noise increases, this affects the model to a high degree, as can be observed in the parameter maps in Figure 3 and in the associated distributions in Figure 4, where a significant broadening in the delay distribution appears with increasing noise.

The proposed model, however, does not use compartmental modelling to fit the tissue delay but rather detects the uprising slope of the TAC and input function separately. On Figure 3A, this model can be observed to result in an overall even but small residual delay across the whole parameter map. The non-zero systematic residual delay is likely due to fitting the upslopes and calculating the onset times individually for the input function and the tissue TAC. Similar behavior has been observed before using straight line tangents to determine the upslope effectively giving an offset between the two curves, in cases of a theoretical delay of 0 s (22). The proposed model is similar, but incorporates a curve around the onset, which is used to correct for this inherent offset, and is the reason to why the fit parameter b, being the tightness of the curve, is present in the calculation of onset time in Equation (3).

Cross-correlation, by nature, strongly favors co-aligning the peaks of the tissue TAC and the input function. In the case where $k_2$ is small (left side of parameter maps in Figures 3C,D) the tracer is not efficiently cleared from the tissue and a tissue TAC will instead accumulate steadily. This effectively moves its peak towards the end which leads cross-correlation to estimate a very large delay between it and the input function, resulting in a large negative residual delay. On the right side of the parameter maps where $k_2$ is large, the TAC will have a very defined peak. Here the curves are delay corrected such that the peaks align, although the peak of the tissue TAC actually appears a few seconds later than the input function peak, resulting in the negative residual delay, however, slightly more positive than for small $k_2$, c.f. in Figures 3C,D.

In the whole-body delay maps (Figure 5), cross-correlation can be observed to frequently estimate larger delays than the other two models, while also exhibiting clearly defined organs. Both are due to peak-weighing of cross-correlation. For example, the brain in the [¹⁸F]FDG delay map (Figure 5E) has a much larger delay for this model than any of the other delay maps, due to the brain's high glucose metabolism causing large amounts of the tracer to accumulate over time, effectively moving the TAC peak near the end where the correlation between the two curves is highest. This high delay would significantly overestimate the kinetic parameters.

As expected due to the free diffusivity of water in tissue, the one-tissue compartment model with incorporated delay performs very well for [¹⁵O]H₂O. Notably the liver is difficult for the proposed model to reliably estimate due to both partial volume effects and motion effects, however, the one-tissue compartment model estimates delays of about 15–20 s corresponding to the venous part of its blood supply, which is also the main component of its blood supply. Ideally the liver should be modeled in relation to a dual input function, which considers both the arterial and venous parts of its blood supplies. This might influence the estimated delays of the one-tissue compartment model as well as the proposed model more than cross-correlation as they are more sensitive to the shape of the TACs, whereas cross-correlation typically favors aligning the peaks.

Also for [¹⁸F]FDG, the one-tissue compartment model shows bias in delay values. Notably for the aorta, where the input function is also derived from, has an average delay of −2 s. The discrepancy here could stem from the fact that FDG has a lower extraction fraction than H₂O and very different kinetics that match best with a two-tissue compartment model, making it noisier on top of the noise at the single-voxel level.

[⁶⁴Cu]Cu-DOTATATE has a very low extraction fraction, especially in the brain, which increases the noise of the tissue TACs significantly. The one-tissue compartment model shows very noisy characteristics in the brain and lungs. The proposed model performs better with a better discrimination of the aorta and lower delays in the surrounding tissues, similar to the other delay maps. However, especially DOTATATE exhibits large amounts of noise due to its very low extraction fraction leading to many different plateaus in the summed TAC, which is a challenge for the proposed model, especially in the peripheral tissues. Experimenting with the plateau parameters here may give better results in specific cases.

From the delay distributions of the tracers in Figure 6, the shape of the organ histograms can be observed to change slightly going from [¹⁵O]H₂O with an extraction fraction close to 1, through [¹⁸F]FDG, and to [⁶⁴Cu]Cu-DOTATATE with a very low extraction fraction. The first histogram is smoother, which relates to how water is easily diffusible from blood to tissue. Whereas for FDG and DOTATATE, the heartbeat can be observed in each organ's histogram. This is due to their lower extraction fractions, effectively resulting in a higher vascular weighing of the signal.

The proposed model also provides more Gaussian distributed organ delays, compared to cross-correlation and the one-compartment model which exhibit troughs and valleys at integer delays. Features such as the heartbeat and the temporal progression of the tracer are also clearly visible.

The perfusion difference in flow in Figure 7C between a mean organ delay correction method and of a voxel-wise delay correction method using the proposed model is shown in Figure 7C. It is apparent that tissue with longer delay, in cases such as carotid artery occlusion, can be underestimated by 15%–25%. Conversely, the kinetics in tissue with shorter delay times compared to the mean organ TAC delay are potentially being overestimated by a similar amount. In the case of Figure 7, the difference in delay between the affected region and healthy brain tissue was around 5 s ± 2 s.

The scope of this study did not include studying tumors but purely to test the performance of the proposed model and not the physiology. A more targeted and organ specific model could be developed which takes into account the different types of tissue, as in the case of mixture models, where tumors are treated differently from other tissue types (23). However, as the proposed model does not take any a priori decisions on the specific tissue type or the tracer used. It is thus independent of tissue type and tracer and should be broadly applicable. For all modeling approaches in general, delay correction is necessary when estimating kinetic parameters to avoid biased estimates.

In this study the tracer was administered as a bolus with an injection time being as fast as possible and was fully completed within the first 40 sec of scan time. Additionally, the injection was administered a few seconds after the scan started to allow for a few dead frames to be used as a baseline. The proposed model has therefore only been tested on bolus injections without considering variations in injection time.

5 Conclusion

We have developed a model capable of estimating the tracer arrival times in various tissues and organs across different tracers exhibiting different properties, such as their extraction fractions, half-life and kinetics. It is important to note that this model was only tested on data from the Siemens Quadra whole-body PET/CT scanner with high-temporal resolution but can in principle be extended to lower time resolution framing images with appropriate validation. Nonetheless, it performs well for a large range of noise levels and as well as for a wide range of tissue TAC shapes as tested using the synthetically produced one-compartment TACs.

It performs better in high-noise environments and for low extraction fraction tracers, such as [¹⁸F]FDG, than the one-tissue compartment model with incorporated delay. This in turn makes it good for using as a single-voxel delay map and can easily be parallelized for shorter calculation times.

Additionally, the organ delay distributions are normally distributed, which the two other models tested here do not produce, and clearly show temporal features such as the heartbeat and the tracer progression through the body and organs.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: Patient confidentiality. Requests to access these datasets should be directed to thomas.lund.andersen@regionh.dk.

Ethics statement

The studies involving humans were approved by Departmental Review board at Department of Clinical Physiology and Nuclear Medicine, Rigshospitalet, Copenhagen. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

FN: Writing – original draft, Writing – review & editing. UL: Writing – review & editing. HB: Writing – review & editing. CJ: Writing – review & editing. IL: Writing – review & editing. BF: Writing – review & editing. FA: Conceptualization, Supervision, Writing – review & editing. TA: Conceptualization, Supervision, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

We would like to thank the technologists and scanner staff for their work.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnume.2024.1360326/full#supplementary-material

References

1. Iida H, Kanno I, Takahashi A, Miura S, Murakami M, Takahashi K, et al. Measurement of absolute myocardial blood flow with H215O and dynamic positron-emission tomography. AHA Journals: Circulation. (1988) 78:104–15. doi: 10.1161/01.cir.78.1.104

Crossref Full Text | Google Scholar

2. Watabe H, Jino H, Kawachi N, Teramoto N, Hayashi T, Ohta Y, et al. Parametric imaging of myocardial blood flow with 15O-water and PET using the basis function method. J Nucl Med. (2005) 46(7):1219–24.16000292

PubMed Abstract | Google Scholar

3. Larsson HBW, Law I, Andersen TL, Andersen FL, Fischer BM, Vestergaard MB, et al. Brain perfusion estimation by Tikhonov model-free deconvolution in a long axial field of view PET/CT scanner exploring five different PET tracers. Eur J Nucl Med Mol Imaging. (2023) 51(3):707–20. doi: 10.1007/s00259-023-06469-w

PubMed Abstract | Crossref Full Text | Google Scholar

4. Weber WA, Schwaiger M, Avril N. Quantitative assessment of tumor metabolism using FDG-PET imaging. Nucl Med Biol. (2000) 27(7):683–7. doi: 10.1016/S0969-8051(00)00141-4

PubMed Abstract | Crossref Full Text | Google Scholar

5. Bertoldo A, Rizzo G, Veronese M. Deriving physiological information from PET images: from SUV to compartmental modelling. Clin Transl Imaging. (2014) 2:239–51. doi: 10.1007/s40336-014-0067-x

Crossref Full Text | Google Scholar

6. Feng T, Zhao Y, Shi H, Li H, Zhang X, Wang G, et al. Tota-body quantitative parametric imaging of early kinetics of 18F-FDG. J Nucl Med. (2021) 62(5):738–44. doi: 10.2967/jnumed.119.238113

PubMed Abstract | Crossref Full Text | Google Scholar

7. Zanotti-Fregonara P, Chen K, Liow JS, Fujita M, Innis RB. Image-derived input function for brain PET studies: many challenges and few opportunities. J Cereb Blood Flow Metab. (2011) 31(10):1986–98. doi: 10.1038/jcbfm.2011.107

PubMed Abstract | Crossref Full Text | Google Scholar

8. de Geus-Oei L, Visser E, Krabbe P, van Hoorn B, Koenders E, Willemsen A, et al. Comparison of image-derived and arterial input functions for estimating the rate of glucose metabolism in therapy-monitoring 18F-FDG PET studies. J Nucl Med. (2006) 47(6):945–9.16741303

PubMed Abstract | Google Scholar

9. Alberts I, Hünermund J-N, Prenosil G, Mingels C, Bohn KP, Viscione M, et al. Clinical performance of long axial field of view PET/CT: a head-to-head intra-individual comparison of the biograph vision quadra with the biograph vision PET/CT. Eur J Nucl Med Mol Imaging. (2021) 48:2395–404. doi: 10.1007/s00259-021-05282-7

PubMed Abstract | Crossref Full Text | Google Scholar

10. Cherry SR, Jones T, Karp JS, Qi J, Moses WW, Badawi RD. Total-body PET: maximizing sensitivity to create new opportunities for clinical research and patient care. J Nucl Med. (2018) 59(1):3–12. doi: 10.2967/jnumed.116.184028

PubMed Abstract | Crossref Full Text | Google Scholar

11. Badawi RD, Shi H, Hu P, Chen S, Xu T, Price PM, et al. First human imaging studies with the EXPLORER total-body PET scanner. J Nucl Med. (2019) 60(3):299–303. doi: 10.2967/jnumed.119.226498

PubMed Abstract | Crossref Full Text | Google Scholar

12. Filippi L, Dimitrakopoulou-Strauss A, Evangelista L, Schillaci O. Long axial field-of-view PET/CT devices: are we ready for the technological revolution? Expert Rev Med Devices. (2022) 19(10):739–43. doi: 10.1080/17434440.2022.2141111

PubMed Abstract | Crossref Full Text | Google Scholar

13. Sari H, Erlandsson K, Law I, Larsson HB, Ourselin S, Arridge S, et al. Estimation of an image derived input function with MR-defined carotid arteries in FDG-PET human studies using a novel partial volume correction method. J Cereb Blood Flow Metab. (2017) 37(4):1398–409. doi: 10.1177/0271678X16656197

PubMed Abstract | Crossref Full Text | Google Scholar

14. Herholz K, Patlak CS. The influence of tissue heterogeneity on results of fitting nonlinear model equations to regional tracer uptake curves: with an application to compartmental models used in positron emission tomography. J Cereb Blood Flow Metab. (1987) 7(2):214–29. doi: 10.1038/jcbfm.1987.47

PubMed Abstract | Crossref Full Text | Google Scholar

15. Dimitrakopoulou-Strauss A, Pan L, Sachpekidis C. Long axial field of view (LAFOV) PET-CT: implementation in static and dynamic oncological studies. Eur J Nucl Med Mol Imaging. (2023) 60:3354–62. doi: 10.1007/s00259-023-06222-3

Crossref Full Text | Google Scholar

16. Li EJ, Spencer BA, Schmall JP, Abdelhafez Y, Badawi RD, Wang G, et al. Efficient delay correction for total-body PET kinetic modeling using pulse timing methods. J Nucl Med. (2022) 63(8):1266–73. doi: 10.2967/jnumed.121.262968

PubMed Abstract | Crossref Full Text | Google Scholar

17. Gunn RN, Gunn SR, Cunningham VJ. Positron emission tomography compartmental models. J Cereb Blood Flow Metab. (2001) 21(6):635–52. doi: 10.1097/00004647-200106000-00002

PubMed Abstract | Crossref Full Text | Google Scholar

18. Nielsen FB. “Delay Map Code (GitHub),” (December 2023). Available online at: https://github.com/fbnie/Dynamic-PET-delaymap (Accessed December 01, 2023).

19. Turku PET Centre. “TPCCLIB 0.8.0 fit_h2o,” (2023). Available online at: https://vesoik.utugit.fi/petanalysis/tpcclib/doc/fit_h2o.html (Accessed December 01, 2023).

20. Wasserthal J, Breit H-C, Meyer M, Pradella M, Hinck D, Sauter A, et al. Totalsegmentator: robust segmentation of 104 anatomic structures in CT images. Radiol Artif Intell. (2023) 5(5):e230024. doi: 10.1148/ryai.230024

PubMed Abstract | Crossref Full Text | Google Scholar

21. Scipy. “Scipy.signal.find_peaks,” Available online at: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html (Accessed December 2023).

22. Meyer E. Simultaneous correction for tracer arrival delay and dispersion in CBF measurements by the H2–15O autoradiographic method and dynamic PET. J Nucl Med. (1989) 30(6):1069–78.2786948

PubMed Abstract | Google Scholar

23. O’Sullivan F, Muzi M, Mankoff DA, Eary JF, Spence AM, Krohn KA. Voxel-level mapping of tracer kinetics in PET studies: a statistical approach emphasizing tissue life tables. Ann Appl Stat. (2014) 8(2):1065–94. doi: 10.1214/14-aoas732

Crossref Full Text | Google Scholar