Extracting a respiratory signal from raw dynamic PET data that contain tracer kinetics

Data driven gating (DDG) methods provide an alternative to hardware based respiratory gating for PET imaging. Several existing DDG approaches obtain a respiratory signal by observing the change in PET-counts within specific regions of acquired PET data. Currently, these methods do not allow for tracer kinetics which can interfere with the respiratory signal and introduce error. In this work, we produced a DDG method for dynamic PET studies that exhibit tracer kinetics. Our method is based on an existing approach that uses frequency-domain analysis to locate regions within raw PET data that are subject to respiratory motion. In the new approach, an optimised non-stationary short-time Fourier transform was used to create a time-varying 4D map of motion affected regions. Additional processing was required to ensure that the relationship between the sign of the respiratory signal and the physical direction of movement remained consistent for each temporal segment of the 4D map. The change in PET-counts within the 4D map during the PET acquisition was then used to generate a respiratory curve. Using 26 min dynamic cardiac NH3 PET acquisitions which included a hardware derived respiratory measurement, we show that tracer kinetics can severely degrade the respiratory signal generated by the original DDG method. In some cases, the transition of tracer from the liver to the lungs caused the respiratory signal to invert. The new approach successfully compensated for tracer kinetics and improved the correlation between the data-driven and hardware based signals. On average, good correlation was maintained throughout the PET acquisitions.


Introduction
In PET imaging, respiratory gating can be used to divide an acquisition into a series of time frames, which in sequence, represent one cycle of respiratory motion. Gated PET images contain reduced motion effects and can be assessed directly, registered and combined, or selectively combined to use gated data only in regions which move (Kesner et al 2013). While hardware based methods remain the standard approach for gating and are supported by most manufacturers of PET cameras, they require additional equipment and are prone to error (Liu et al 2009, Wells et al 2010. Data-driven gating (DDG) provides an appealing alternative as this approach extracts the respiratory signal directly from acquired PET data, requires no additional hardware, and usually does not require any modification to the PET acquisition protocol.
Several data-driven respiratory gating methods have been developed for PET imaging, and many of these rely on the principle that the PET counts within a respiratory-motion affected region of a PET acquisition will vary over time in correlation with the respiratory-motion itself (Visvikis et al 2003, He et al 2007, Kesner et al 2007, Büther et al 2009, Schleyer et al 2009, Kesner and Kuntner 2010, Thielemans et al 2011. In some applications, dynamic PET imaging is used to measure tracer kinetics to provide quantitative information un-obtainable in static imaging. Dynamic PET, acquired from the time of tracer injection, is required for absolute myocardial perfusion analysis, and correcting for respiratory motion using respiratory gated PET has been shown to improve quantitative accuracy (Pourmoghaddas et al 2013).
In dynamic PET imaging, we expect that the spatial distribution of tracer will vary during the acquisition. This can negatively impact data-driven gating where it is assumed that regions affected by respiratory motion are not transient. Here, we used dynamic NH 3 cardiac PET data to demonstrate how a DDG method that automatically identifies motion-affected regions fails in the presence of tracer kinetics. We then describe a method that analyses the PET data as a series of optimised temporal windows to allow for tracer kinetics.

Methods
The DDG method used is referred to as the Spectral Analysis Method (SAM), which uses spectral analysis to produce a static 3D mask that identifies respiratory-motion affected regions in raw PET data (Schleyer et al 2009(Schleyer et al , 2011. In this process, the fast-Fourier transform (FFT) of each sinogram location is calculated in the temporal domain over the entire duration of the scan unlisted to spatially down-sampled, smoothed, and decay corrected short duration sinograms (here we use 500 ms), g(t). The sinogram locations that demonstrate a high spectral magnitude within a pre-defined frequency window that represents respiration (here we use 0.1-0.4 Hz) are included in the mask.
In the SAM approach, a respiratory signal for each time point, r(t), is then obtained by taking the weighted sum of values in the sinograms, where w k is 0 if sinogram bin k is not in the mask. The choice of the weights, w k , is important. As illustrated in figure 1, the contents of a region of interest (ROI) defined over the leading edge of a moving object cancel out the contents of an ROI defined over the object's opposing trailing edge if both ROIs are weighted equally. Therefore, to let r(t) represent the size of motion, SAM assigns a negative weight to one of the ROIs. However, as described in Schleyer et al (2009), the weight values only indicate if edges are opposing or equal and it is not known if edges are trailing or leading. It follows that the absolute direction of motion is not defined by r(t), i.e. the relationship between the sign of the weight values (and of the respiratory signal) and physical direction is arbitrary.
An overview of the SAM method is illustrated in figure 2.
As the FFT is calculated over the entire duration of the scan, a static mask is produced which does not allow for a change in tracer distribution. An alteration in tracer distribution within the mask will contribute to the overall respiratory signal. Additionally, tracer moves within the body during the dynamic acquisition leading to a time varying tracer distribution in different organs, so the regions in the sinograms where the motion can be detected will also move. To overcome this, we used a short-time Fourier transform (STFT)  SAM method uses an FFT over the entire acquisition duration to produce a single 3D mask. The KRG method uses STFT windows to produce a 4D mask to allow for tracer kinetics. The sign of each mask must be chosen to ensure consistency.
to analyse the data as a series of contiguous temporal segments so that a single mask is generated for each segment. This approach, referred to as Kinetic Respiratory Gating (KRG), is illustrated in figure 2.
In selecting the temporal width of the STFT window, a trade-off exists between the amount of data available to the spectral analysis (a wider window produces more sample points and less noise) and the temporal resolution of the mask (a wider window decreases sensitivity to tracer kinetics); we require the widest window possible without compromising the sensitivity to tracer kinetics. However, the rate of tracer kinetics within the body will not be constant throughout the acquisition. Specifically, immediately after injection the tracer is distributed rapidly in blood before it transfers more gradually into various tissues, typically at a rate which decreases over time as equilibrium is approached. As such, we would expect that the optimal STFT window size will also vary during the acquisition according to the rate at which the tracer distribution changes. For this reason, KRG varies the width of the STFT window depending on the rate at which tracer distribution changes, and this is estimated using a sinogram-derived measure, y′(t), described below. In the separate optimisation process described below, patient data was used to determine the best window width to use during specific ranges of tracer-change rate.
The tracer-change, δ(t), was approximated by calculating the sum of differences between each 500 ms sinogram in the decay corrected series g(t) and the last 500 ms sinogram, g (final), where h is a 12 s box-car filter, and bins refers to all sinogram bin locations. To produce a smoothed and monotonic approximation of tracer-change, a dual-exponential, y(t), was fitted to the section of δ(t) from its peak to the end of acquisition, and the rate of tracer-change was approximated with the gradient function y′(t). This function does not correctly describe the tracer-change rate resulting from the first-pass of tracer distribution seen prior to the curve peak (see figure 5), however as we show in the results this is not required as only one stationary-window width produced the highest correlation before the curve peak for all patients; one STFT window-width was required for all time-points prior to δ(t) peak in all cases. 53 NH 3 PET studies were used to demonstrate the KRG method. Studies were acquired in 2D mode (with the septa in the field of view) for 26 min from the time of injection. For comparison, a hardware based respiratory signal was acquired with a Varian RPM device (Varian Medical Systems, Palo Alto, CA) for all PET acquisitions. Images were acquired on a GE Discovery VCT PET/CT scanner as a dynamic, with the following frames as typically used for cardiac tracer kinetic analysis: 12 × 10 s; 6 × 20 s; 2 × 60 s; 1 × 20 min. Data driven and hardware derived signals were compared over the entire duration of the acquisition, and over each of these dynamic frames.
The two processes of STFT optimisation and correcting for inconsistent mask-sign are described below.

STFT width optimisation
This process was performed as a training step to determine the best STFT windows. Once completed, KRG operates independently of any hardware based signal. Based on the assumption (described above) that the best window size will largely depend on the rate at which tracer moves in the body, the optimal STFT windows were defined as a function of y′(t).
Each data-set was first processed with static STFT widths of 10,14,16,18,20,40,60,120,180,300,780 and 1560 s (where 1560 s is the full dataset and is equivalent to SAM). For each window width, correlation between the data-derived and RPM derived respiratory signals was calculated as a function of time 3 . By averaging over all patients, a mean correlation time-series was produced for each window width, and these were plotted against the mean y′(t) (also averaged over all patients). The highest correlation with RPM was not produced by one window width for all of y′(t); different stationary widths produced the highest correlation for different sections of y′(t). These widths, and the corresponding ranges of y′(t) where they produced the best correlation, parameterised the optimal non-stationary STFT width.

Sign correction
A 3D phase weighted mask was created at each contiguous STFT window using the SAM weighting values described in Schleyer et al (2009). As described above, the (overall) sign of the weight values and therefore of the respiratory signal are arbitrary. As each STFT window is evaluated independently, a process was required to ensure that the relationship between respiratory signal sign and the direction of respiratory motion remained consistent for all STFT windows.
This was achieved by inverting the mask sign as necessary, in order to minimise the difference between neighbouring masks. For STFT window number n in the contiguous series of N masks (corresponding to the series of N STFT windows) a sign of α n was assigned to mask w n that minimised where d is the neighbourhood size (the number of previous masks to consider, d = 5 in this work), and α ν is the sign for mask w ν . Following this, the 4D mask, created by stacking the 3D masks, was smoothed in the temporal domain with a Gaussian kernel width of 3 temporal indices. Finally, a segment of respiratory curve, r n (t), was obtained for each contiguous STFT window by applying the corresponding 3D mask to the sinogram data acquired during the STFT temporal window defined by times t n 1 and t n 2 , and summing over all sinogram bins, (4)

Overall respiratory signal
An overall respiratory signal was created by z-score normalising the respiratory curve from each segment (by subtracting the mean and dividing by the standard deviation), and combining all segments in sequence.

Signal comparison
Before comparison, complete data-driven and RPM curves were straightened by subtracting the low frequency content and then z-score normalised. Data-driven and RPM curves were compared by calculating the correlation between the overall signals, and between the signals during each of the dynamic frames as defined by the acquisition protocol. All RPM curves were interpolated to match the temporal sampling of the data-driven curves.

Results
An example NH 3 patient is shown in figure 3. Good correlation can be seen between the SAM and RPM during the initial 300 s of the acquisition, after which the correlation becomes negative and the SAM respiratory signal is inverted. This corresponds to a reversal in activity gradient across the liver and lungs. Good correlation is seen between RPM and KRG signals throughout the acquisition.

STFT window optimisation
The mean correlation (averaged over all patients) between each stationary STFT derived respiratory signal and the RPM signal is plotted against the mean gradient of tracer-change, y′(t) in figure 4. The 20 s and 40 s STFT windows are shown to produce the highest correlation across the entire range of y′(t). Notably, the 40 s width produced the highest correlation in the range y′(t) > − 0.0294, and the 20 s window in the range y′(t) < − 0.0294, and these were chosen as the optimal STFT window parameters. The 40 s and 20 s window widths produced a statistically significantly different correlation (p < 0.05) in both of these ranges. Figure 5 indicates the two STFT window widths for an example patient, showing the tracer change, δ(t), and the fitted dual-exponent, y(t). Over all patients, the mean time-point where y′(t) = − 0.0294 was 42 s (±8.7 s), and this was later than the time at δ(t) peak in every case (on average, 28 s later). The remainder of the results are using the optimised STFT parameters. Figure 6 shows an example patient where the correlation between STFT-based data-driven signal (without sign-correction) and the RPM signal becomes inverted during the acquisition. DDG resp. curve

RPM SAM KRG
The sign of the correlation remains positive throughout the entire acquisition after correcting the sign of the masks. After applying sign-correction, a sign-constant RPM correlation was found for all dynamic frames, for all patients.

Signal comparison
Left-ventricle (LV) activity was measured by defining a region of interest on non attenuation corrected reconstructed dynamic PET images. On average, LV activity peaked at 35 s. As expected, early dynamic frames with little or no activity in the field of view produced DDG traces that correlated poorly with the RPM traces. This is exemplified in figure 7 which shows the first four 10 s frames for a patient. The first 2 frames demonstrate low correlation, however 0% and 0.6% of the peak LV activity were present in these frames. Initial frames with very low activity are not expected to contribute to kinetic image analysis, so frames before the peak with an LV activity of less than 5% of the peak were excluded from the evaluation. The mean correlation (for the entire scan) with the RPM signal for all patients was 0.5 (±0.2) for SAM, and 0.75 (±0.08) for KRG. The mean (over all patients) SAM and KRG correlation with RPM is shown for individual dynamic frames in figure 8 (including the relative LV count-rate averaged over all patients). The mean correlation with RPM was above 0.7 for all dynamic frames except for the second frame (10-20 s) where the mean correlation was 0.59. All first frame images (0-10 s) contained an LV activity of less than 5% peak and were excluded.

Discussion
Tracer kinetics can severely degrade data-driven respiratory signals. In some cases the transfer of tracer across respiratory affected boundaries (such as between the liver and lungs) can cause an inversion of the activity gradient, so that the relationship between derived respiratory  signal sign and direction of motion also inverts. Inverting the respiratory signal where the activity gradient inverts could partly resolve this, however as seen in figure 3 the respiratory signal gradually degrades over a period of approximately 100 s, so applying a signal inversion can not completely restore accuracy during this phase. Additionally, the time when the activity gradient inverts is not known without assessing the image data. By using a short-time Fourier transform to generate a series of masks and combining the respiratory signals from each segment, these issues were overcome. Importantly, however, the sign of the mask produced at each STFT segment must be corrected to ensure consistency throughout the acquisition.
Data-driven gating methods derive a respiratory signal from regions of activity that move with respiration. As very low activity may be present in the field of view immediately following injection, it can be expected that data-driven methods may not recover an accurate respiratory signal during the first few seconds of the acquisition. Our results confirmed the correlation between KRG and RPM signals was typically lower for the first 2 or 3 dynamic frames (20-40 s in total), but remained generally high for the rest of the acquisition. Reduced correlation is likely to indicate less effective respiratory gating, however we assume that this is of no consequence as these frames contain little or no activity (less than 5% of the LV peak); early frames with LV activity greater than 5% demonstrated high correlation (0.59 for frame 20-30 s and 0.73 for frame 30-40 s).
Optimised non-stationary STFT windows were used in this work, and this allows for the movement of tracer distribution typically seen in dynamic images acquired from the time of injection. As the best STFT window size is expected to vary with different rates of tracer change, the non-stationary window width was assigned according to the rate of tracer-change which was estimated directly from sinograms. This allows for inter-patient differences in tracer kinetics and does not rely on strict adherence to imaging protocol. However, the non-stationary STFT temporal window size was only optimised for the NH 3 cardiac perfusion PET studies in this work, and further optimisation may be require to apply the method to tracers other than NH 3 .
We expected that shorter STFT windows would be required for the initial stages of acquisition when tracer change is typically the fastest, and longer STFT windows later in the acquisition when the rate of tracer-change reduces. However, 40 s windows were found to be optimal for the first 42 s of the acquisition (on average), and 20 s windows for the remainder of the acquisition. A likely explanation for this unexpected result is that the very low activity in the field of view during the initial 20 s of imaging was insufficient for the spectral analysis to produce a meaningful mask. The longer 40 s window contained adequate activity to form a Figure 7. An example of the first 4 10 s dynamic frames. KRG and SAM traces appear to match the RPM trace closely from approximately 12 s. The first 2 frames demonstrate low correlation between DDG and RPM tracers, however 0% and 0.6% LV relative activity was present in these frames. mask, and although tracer-change rate is likely to be high in this early phase this was of no consequence to the respiratory signal as an inversion of activity-gradient was not seen in the early blood-flow data.
In dynamic PET imaging, early frames which represent the initial blood-pool distribution of the injected tracer are generally used to obtain an image-derived input function. Respiratory motion can contribute to partial volume effects where the input function is obtained from the LV (Du et al 2013), and respiratory gating Rubidium-82 PET has been shown improve the accuracy of absolute myocardial perfusion quantitation (Pourmoghaddas et al 2013). However, the effect of gating on other applications of quantitative dynamic imaging are currently not known, and the approach presented offers the potential to retrospectively respiratory gate dynamic data acquired without hardware-based measurements of respiration. Phase matched attenuation correction data remains a requirement for respiratory motion correction, specifically where airtissue boundaries are present such as for cardiac imaging (Chin et al 2003, Meunier et al 2006, Martinez-Möller et al 2007, Gould et al 2007, McQuaid and Hutton 2008. 4D CT is frequently used for this, however pseudo gated 4D AC data can be generated by applying PET-derived motion fields to static 3D CT images (Dawood et al 2006, Fayad et al 2008, Wells et al 2010, McQuaid et al 2011, which has been shown to be robust (Schleyer et al 2013).
To allow for tracer kinetics, a data-driven respiratory signal was obtained by analysing the acquisition as a series of temporal windows and then combining the respiratory signal from each window. Several other existing data-driven gating methods are likely to also fail in the presence of tracer kinetics, and the windowed approach used in this work could potentially also be used for these methods , He et al 2007, Büther et al 2010, Kesner and Kuntner 2010, Thielemans et al 2011.  Computation time was relatively fast. Using Matlab 2013a (Mathworks, Natick, MA) in single threaded mode running on a 3.3 GHz processor, unlisting the listmode file into the 500 ms series took approximately 70 s, the KRG method took approximately 11 s, and the SAM method approximately 5 s. All images were acquired in 2D mode in this work, and while the KRG method was not tested with 3D mode dynamic data we have previously found that the SAM approach works well with 3D mode static data (Schleyer et al 2011).
A caveat of the evaluation presented in this paper is the use of an RPM hardware-based signal as the gold standard respiratory signal. This may introduce error, as the chest wall position may not be a true surrogate for internal-structure respiratory motion (Ozhasoglu and Murphy 2002, Hoisak et al 2004, Fayad et al 2011, i.e. a correlation of less than 1 does not necessarily indicate less accurate gating. In some cases a time-varying phase difference between external and internal can exist (Ozhasoglu and Murphy 2002), and in figure 7 a slight phase shift can be noted between the RPM and DDG signals between 10-20 s. While we found a generally high correlation between the DDG and RPM signals, consideration should be given to this limitation.

Conclusion
Data-driven gating methods that extract a respiratory signal by observing the counts within motion-affected regions can fail in the presence of tracer-kinetics. This was successfully overcome by analysing the PET acquisition as a series of temporal windows with a non-stationary STFT, and ensuring a consistent relationship between signal sign and direction of motion for all windows. The non-stationary STFT was optimised for NH 3 cardiac perfusion PET images, such that the STFT windows were defined according to an acquisition-specific estimate of tracer change which was derived directly from sinograms. 53 cardiac perfusion PET images were tested, and on average, good correlation between the data-driven respiratory signal and a hardware based signal was found.