Use of dynamic reconstruction for parametric Patlak imaging in dynamic whole body PET

Dynamic whole body (DWB) PET acquisition protocols enable the use of whole body parametric imaging for clinical applications. In FDG imaging, accurate parametric images of Patlak K i can be complementary to regular standardised uptake value images and improve on current applications or enable new ones. In this study we consider DWB protocols implemented on clinical scanners with a limited axial field of view with the use of multiple whole body sweeps. These protocols result in temporal gaps in the dynamic data which produce noisier and potentially more biased parametric images, compared to single bed (SB) dynamic protocols. Dynamic reconstruction using the Patlak model has been previously proposed to overcome these limits and shown improved DWB parametric images of K i . In this work, we propose and make use of a spectral analysis based model for dynamic reconstruction and parametric imaging of Patlak K i . Both dynamic reconstruction methods were evaluated for DWB FDG protocols and compared against 3D reconstruction based parametric imaging from SB dynamic protocols. This work was conducted on simulated data and results were tested against real FDG dynamic data. We showed that dynamic reconstruction can achieve levels of parametric image noise and bias comparable to 3D reconstruction in SB dynamic studies, with the spectral model offering additional flexibility and further reduction of image noise. Comparisons were also made between step and shoot and continuous bed motion (CBM) protocols, which showed that CBM can achieve lower parametric image noise due to reduced acquisition temporal gaps. Finally, our results showed that dynamic reconstruction improved VOI parametric mean estimates but did not result to fully converged values before resulting in undesirable levels of noise. Additional regularisation methods need to be considered for DWB protocols to ensure both accurate quantification and acceptable noise levels for clinical applications.


Introduction
Positron emission tomography (PET) imaging is well known and established in clinical applications and pathways, with an important role towards the delivery of precision medicine (Subramaniam 2017). The established clinical practices rely on static imaging after a certain uptake period and semi-quantitative measures, such as the standardised uptake value (SUV). But these measures are vulnerable to many unknown factors that can vary between PET examinations, such as body composition, retention clearance, inconsistencies in uptake time and imaging practices (Boellaard 2011). On the other hand, dynamic PET imaging can be used to fully characterise underlying tracer kinetics and provide fully quantitative measures that could overcome many limitations of current static imaging practices and enable the use of PET for new applications in clinical practice (Lammertsma 2017, Meikle et al 2020, Dimitrakopoulou-Strauss et al 2021. Current clinical scanners are limited in coverage by their axial field of view (A-FOV), with values ranging from 15 to 26 cm (Vandenberghe et al 2020). This is sufficient for single organ dynamic studies but cannot directly provide synchronous wholebody coverage, which is essential for some clinical applications such as tumour staging in oncology. In practice for static imaging whole-body coverage is achieved using multiple bed positions at different axial locations to provide the desired axial coverage (Schubert et al 1996), or alternatively via continuous bed motion (CBM) during the acquisition (Panin et al 2014). Recently scanners with increased A-FOV have been developed , Siegel et al 2020, even with nearly 2 m long A-FOV which provides total-body coverage (Cherry et al 2018). But these scanners are still not widely adopted in the clinic. Using similar methods as in static whole-body imaging, dynamic whole body (DWB) protocols have been developed using multiple bed positions and repeated whole-body passes (Karakatsanis et al 2011, 2013, Rahmim et al 2019. These types of acquisition protocols have also been incorporated into clinical products (Hu et al 2020), and it has been shown that their use in clinical practice is feasible (Fahrni et al 2019, Dias et al 2020. The immediate effect of transition from single bed (SB) dynamic studies to multi-bed dynamic studies is the introduction of temporal gaps in the acquired data of any given bed position. These are introduced at each bed position by the time spent on imaging other bed positions and by scanner system delays due to the time required to move the bed to the next position and prepare for the next acquisition. These gaps cause a significant reduction in the sensitivity of the acquisition, with fewer total counts collected for each axial location when compared to SB dynamic acquisitions. Furthermore, estimation of fast temporal changes in tracer uptake is compromised as the early time points of the acquisition are not fully sampled for all beds. Finally, the established clinical protocols that make use of image derived input function (IDIF) to ease integration in clinical practice further sacrifice imaging time in the studys early phase, which is spent in acquiring fast frames over a SB location centred over the heart and the aorta (Hu et al 2020).
The generation of parametric images from dynamic data requires the fitting of the dynamic model of interest on time activity curves (TAC) for every voxel in the image. Due to the poor statistics and high noise associated with TAC measurements at the voxel level, in particular for DWB acquisitions, parametric image estimates can be heavily corrupted by noise and potentially biased. The use of direct dynamic reconstruction has been proposed to improve on this task by making use of dynamic models directly in the reconstruction. These techniques allow for more accurate modelling of the noise from the raw PET data in the generation process of parametric images and can improve parametric image noise and reduce bias (Reader and Verhaeghe 2014). For DWB acquisitions specifically, it has been shown using simulated and real data that direct dynamic reconstruction provides reduced noise, bias and improved parametric image contrast when compared to postreconstruction parametric imaging (Karakatsanis et al 2016). So far, dynamic reconstruction has been used in DWB applications using the Patlak model and a generalised Patlak model for K i parametric imaging, using only frame data after sufficient time from injection to satisfy the steady state conditions of the Patlak model.
In this work, we evaluate the performance of dynamic reconstruction algorithms for DWB fluorodeoxyglucose ([ 18 F]FDG) PET imaging for various dynamic reconstruction methods and different DWB acquisition protocols. The novelty in this work is the application of a spectral analysis based dynamic reconstruction on DWB data for Patlak K i parametric imaging, which enables the use of all acquired DWB data from injection in the reconstruction. The evaluation is based on simulations of a unique bed position for a SB and various multi-bed dynamic acquisition protocols and results are illustrated in a real dynamic FDG PET study. An additional novelty in this work is the comparison of DWB results against SB dynamic data to assess the degradation effect of the temporal gaps, which are introduced in the transition from SB studies to multi-bed DWB studies, and to assess the potential of dynamic reconstruction in counteracting these effects. Our work also evaluates the use of different optimization methods for dynamic reconstruction in the setup of DWB imaging.
In detail, we evaluate for WB Patlak K i parametric imaging (1) the benefits of using direct Patlak dynamic reconstruction in DWB protocols against SB dynamic protocols and indirect parametric imaging from regular 3D reconstruction, (2) the use of the Spectral analysis dynamic model (Cunningham and Jones 1993) in dynamic reconstruction for indirect parametric K i imaging, (3) the use of two different optimization algorithms for dynamic reconstruction and finally (4) the impact of different DWB acquisition strategies.

Simulated acquisition protocols
A SB dynamic protocol (continuous in time with no temporal gaps) and three DWB acquisition protocols of five bed positions (with temporal gaps) were considered for the simulation study, using the geometry characteristics of the GE Signa PET/MR scanner (Grant et al 2016). With the provided 25 cm A-FOV per bed and a bed overlap of 3.34 cm, axial coverage of 110.3 cm can be achieved with five bed positions. The relatively small overlap (approximately half compared to routine clinical protocols) was selected to reduce the number of beds and subsequently the acquisition temporal gaps. A total study duration of 60 min was used in the design of all protocols, including an initial SB dynamic phase of 3 min centred over the aorta to mimic requirements for IDIF estimation.
• The first DWB protocol (DWB-1) considers a step and shoot (S&S) acquisition using measured system delays from an experimental DWB protocol on the Signa PET/MR, resulting in a delay of 6 s between adjacent bed positions and 36 s between whole-body sweeps, for a total of 8 whole-body sweeps in the duration of the study.
As the system has no dedicated DWB protocol, a large fraction of the system delays between whole-body sweeps is caused by data transfers and internal processes.
• The second DWB protocol (DWB-2) differs from DWB-1 by use of system delays that mimic a CBM acquisition of the same length, with no delays between adjacent bed positions and 12 s delay between wholebody sweeps to account for bed speed and acquisition overscan (used to obtain reasonable axial sensitivity at the edges of the acquisition (Panin et al 2014)). This protocol results in 9 whole-body sweeps. In this protocol it was assumed that the CBM acquisition is dedicated and optimized for DWB acquisitions. The physical time required for translation of the bed between whole-body sweeps was deduced by measurements of table movements on the Signa PET/MR scanner. In this work, the actual simulation for DWB-2 did not make use of CBM within the simulation process but made use of the accurate timings that reflect the reduction of delays and framing achievable with CBM on the same system geometry.
• The third DWB protocol (DWB-3) replicates the timing properties of DWB-2 for CBM acquisition but utilises a bi-directional acquisition that reduces delays between sweeps to the time spent for the over-scan. This acquisition motion provides more sweeps in the same study duration but results in non-uniform sampling. This protocol results in 10 whole-body sweeps.
In this study we assumed that individual bed positions of the DWB studies are to be reconstructed and processed independently and the result K i parametric maps to be combined after reconstruction for the creation of WB parametric K i maps, as proposed by Karakatsanis et al (2016). Due to practical limitations in reconstruction time, the simulation study was focused on a SB position from the assumed DWB protocols of five bed positions. The second bed position from the top was selected, centred over the upper chest as seen in figure 1. Different bed positions will result in different sampling of frames, resulting in an offset in sampling for the S&S and uni-directional CBM DWB protocols, but most importantly resulting in sampling uniformity differences for bi-directional CBM acquisition. These sampling uniformity differences are eliminated at the central bed position and maximised at the edge positions (first and last bed positions) of the assumed 5-bed DWB protocol. The choice of the second bed position was made as a compromise between the central and edge bed positions, to include non-uniformity effects in the analysis of the DWB-3 protocol results. The PET data simulations were conducted solely for this bed position using the framing of the DWB protocols described above. The exact framing information is available in the supplementary material.

Digital phantom and analytical simulation
The Zubal brain phantom (Zubal et al 1994) was chosen for the simulations of PET [ 18 F]FDG data, even though its anatomy doesn't correspond to the anatomy that would be found in the axial location of the simulated bed position. The choice of the phantom was made to incorporate higher complexity structures than those offered from common lung/chest phantoms. Furthermore, the placement of the brain phantom in the centre of the FOV resulted in no area of the phantom falling in the overlapping regions of the DWB acquisition protocols. Therefore, the simulation of adjacent bed positions was not necessary for this study. The behaviour of structures falling in the overlapping regions of DWB studies as well as the treatment of the overlap data during and after reconstruction is yet another subject for investigation. The Zubal brain phantom was segmented into 19 unique regions and a non-reversible two tissue compartment model was assigned uniformly to each region to simulate realistic FDG kinetics, with K 1 , k 2 , k 3 and V b values drawn from the literature and a real measured input function. A selection of the simulated kinetic parameters is provided in the supplementary material. An analytical simulator was used to generate raw PET sinogram data (Stute et al 2015). The simulations included attenuation and detector resolution effects, scattered and random coincidences, while Poisson noise was added to the sinogram data. The same attenuation coefficients and random and scattered coincidences analytical distributions used in the simulation (that is, the analytical distributions that were used in the simulation before adding Poisson noise) were used in the reconstruction, resulting in unbiased corrections. In addition, the exact input function used within the simulation process was used for all dynamic reconstructions. The simulation did not include time of flight (TOF) information in the data. Fifty different noise realisations were simulated for each DWB protocol. For the SB protocol the number of noise realisations was reduced to twenty for practical reasons, as it involves a much higher number of frames resulting in substantially longer reconstruction time.

Reconstruction and kinetic modelling
Both 3D and dynamic iterative reconstructions were used in this study. All reconstructions were run within an OSEM framework, for 40 iterations and 28 subsets, using the open source fully quantitative reconstruction platform CASToR (Merlin et al 2018). All reconstructions were performed with a voxel size of 2.2 mm × 2.2 mm × 2.8 mm and included resolution modelling as well as corrections for attenuation and for random and scattered coincidences (generated from the simulation). Dynamic reconstructions are made by combining dynamic models that describe the tracer kinetics with the tomographic reconstruction process. When the dynamic model of interest is used this technique results in direct reconstruction of parametric images of interest and allows for accurate modelling of the raw data noise in the estimation process (Carson and Lange 1985, Matthews et al 1995, Kamasak et al 2005, Wang et al 2008. The use of generic dynamic models can also be made for dynamic reconstruction, where dynamic models impose temporal regularisation in the frame activity estimation process. Post-reconstruction (indirect) estimation of parametric images can then be made and indirectly benefit from the use of dynamic reconstruction Verhaeghe 2014, Novosad andReader 2016). Linear dynamic models can be directly implemented within the system matrix (Matthews et al 1995, Wang et al 2008, Reader and Verhaeghe 2014 but result in algorithms with substantially slower overall convergence properties. Instead a nested optimization framework Qi 2010, Matthews et al 2010) can be used, which decouples the dynamic model fitting process on image space data from the tomographic update process over the raw PET data. This allows for multiple nested sub-iterations of dynamic model fitting to be run within each tomographic iteration of the reconstruction process, resulting in convergence acceleration and reasonable computing time requirements. The separation of the two processes allows for various optimization algorithms to be implemented in the nested dynamic model fitting process (Matthews et al 2010). One particular method of interest is the non-negative least squares (NNLS) algorithm (Lawson and Hanson 1995), which enforces non-negativity and is commonly used in post-reconstruction kinetic modelling. The NNLS algorithm is self-terminating and for dynamic reconstruction it has been shown that a single execution of nested NNLS results in parametric images with similar root mean square (rms) error to that of 15 iterations of nested MLEM (Matthews et al 2010). Therefore the use of nested NNLS optimization has the potential for reduced overall reconstruction time. In this work we made use of both nested MLEM and NNLS optimizations to compare their performance and reconstruction time requirements for their implementation within CASToR. The nested MLEM optimization was used with 20 sub-iterations of the dynamic model fitting process after each subset of the OSEM tomographic update process, which has been found to be the optimal number of sub-iterations in a previous similar study (Karakatsanis et al 2016). Hereafter we will refer to nested dynamic reconstruction simply as 4D reconstruction.
In this work we evaluated two different linear dynamic models for 4D reconstructions, the Patlak model and the spectral analysis model.
• 4D Patlak: The Patlak model describes the activity in tissue C where C P is the activity concentration in arterial blood plasma at time point t, K i is the steady state trapping rate and V α is the apparent volume of distribution. The Patlak model is valid once steady state conditions have been reached (denoted as t ss ). For a PET measurement the observed activity is where conventionally it is assumed that the blood fraction V B is small ( 0.05) in most tissues. If we assume for FDG that the total blood activity concentration C B is proportional to C P with C B = rC P , define the Patlak slope θ 1 = (1 − V B )K i and the Patlak intercept θ 2 = V α + rV B , then the observed activity of acquisition frame f between time points t start and t end is modelled according to the Patlak model as is the observed activity map. Using this representation a linear model of two basis functions can be constructed, which when fitted to TAC data provides parametric images θ = [θ 1 , θ 2 ]. 4D dynamic reconstruction with the Patlak model directly results in parametric images of θ 1 and θ 2 . In our study 4D Patlak reconstruction was applied using frame data after the first 15 min, from where we assumed steady state conditions (t 15 min ss = ). It is important to note that a limitation of the Patlak model is that the estimated K i from the Patlak slope θ 1 is susceptible to systematic errors in its estimation and can deviate from the true underlying K i ( = K 1 k 3 /(k 2 + k 3 )). In addition, V B is not necessarily known a priori and Patlak analysis can not distinguish between K i and (1 − V B ). In this study we use the Patlak slope θ 1 as the K i value of interest for parametric imaging.
• 4D spectral: The spectral analysis method has been proposed by Cunningham and Jones (1993) for the analysis of TAC data from dynamic PET studies. The method makes use of a pre-defined set of exponential functions, convolved with the input function, to describe the measured TAC data and therefore relies on few assumptions. It can be used to analyse and compare data, while also allow for the deduction of information of the underlying kinetic behaviour (Gunn et al 2002). It was first used in dynamic PET for parametric imaging by Meikle et al (1998) and later as a generic temporal regularisation method within dynamic reconstruction in various works (Reader et al 2007, Matthews et al 2010, Wang and Qi 2010. In the latter, the resulting temporally regularised activity data per frame can be used for post-reconstruction parametric imaging. This process provides improved estimates, compared to post-reconstruction analysis from regular 3D frame reconstructions. This has been shown previously for post-reconstruction analysis using a two-tissue compartmental model (Novosad and Reader 2016). In this work, we propose a novel use of the spectral model in dynamic reconstruction for regularisation of activity frame data, followed by post-reconstruction Patlak analysis for parametric K i imaging. We will refer to this dynamic reconstruction method as the 4D Spectral reconstruction. The expected benefit of this practice is that the spectral model can make use of all frame data within the reconstruction, allowing all the dynamic acquisition statistics to be indirectly used within postreconstruction Patlak analysis. By contrast, the Patlak model based dynamic reconstruction is limited to frame data after steady state conditions have been reached. In the proposed practice, the 4D Spectral reconstruction for parametric K i imaging can be regarded as indirect dynamic reconstruction.
For a measurement within an acquisition frame f between time points t start and t end , the observed PET activity can be described according to the spectral analysis method with M + 1 number of parameters, decay rates β and their positive weights f as Assuming C B is proportional to C P then the parametric map f M is proportional to the blood fraction V B , while for irreversible kinetics the decay rate β 0 → 0 and the parametric map f 0 describes tracer trapping. Parametric maps f 1 ...f M−1 describe the exchange between compartments, with decay rates β 1 ...β M−1 chosen to be logarithmically spaced within a range of values that covers the expected underlying kinetics. In our tests we used 3 different sets of numbers of basis functions (M + 1 = 17, 9 and 6), with β 1 ...β M−1 logarithmically spaced within the range of 3-0.001 min 1 -. Although post-reconstruction spectral analysis makes use of hundreds of basis functions, this number is reduced when the model is incorporated in dynamic reconstruction to reduce the number of estimated parameters and variance of estimates. The choice is commonly set to less than 100 (Reader et al 2007, Novosad andReader 2016) and in this study it was chosen to be close to the number of frames per bed position of the DWB PET data. We made use of three sets of basis, to evaluate the performance of the model when the number is lower, similar or greater than the number of available frames per bed position in the DWB data. Unlike the Patlak model, the spectral analysis model is valid from the start of the acquisition and by default was applied to all time frames. The parameters f b of the spectral model have physiological meaning and in combination can be used to derive macro-parameters maps such as K 1 and either K i or V D , depending on the irreversible or reversible kinetic behaviour (Gunn et al 2002). However, this derivation implies that the acquisition starts at the injection time, which is not the case for DWB protocols, except for the bed position corresponding to the initial dynamic phase used for IDIF estimation. In this sense, the use of 4D Spectral reconstruction for parametric K i imaging can be regarded as indirect dynamic reconstruction.
As highlighted above the spectral analysis model makes use of all frame data, while the Patlak model uses data after t ss . In order to make a closer comparison between the two models for 4D reconstruction, an additional comparison was made using only data after t ss (reconstructions labelled with t > t ss ).
With the exception of 4D Patlak reconstructions that directly output parametric images of K i , all other 4D and 3D reconstruction activity maps were fitted post-reconstruction with the Patlak model at the voxel level, using the ordinary least squares (OLS) optimization algorithm, to generate parametric K i images. For all 4D reconstructions and post-reconstruction fitting processes the true input function was used.
The reconstructionʼs namings and parameters are summarised in tables 1 and 2.

Evaluation metrics
The reconstructed and generated parametric K i images were evaluated across noise realisations for voxel based and VOI based metrics. We define θ j,n as the image K i value for voxel j in noise realisation n, θ VOI,n the VOI K i mean value and VOI GT q the ground truth value. The following voxel-based metrics of bias and coefficient of variation (COV) were calculated, including their rms spatial average within VOIs using equations (7) The true simulated K i value (K i = K 1 k 3 /(k 2 + k 3 )) was used as the ground truth value VOI GT q . Because Patlak analysis provides different fits on the simulated TACs depending on the DWB protocol framing, noiseless TAC Patlak fits were conducted to estimate these differences. These values are provided in table 4, for comparison against the findings presented in results section 3.4 on the comparison of DWB protocols.
The cortex and an eroded thalamus VOIs, defined from the respective regions of the Zubal brain phantom, were evaluated in the analysis. The thalamus VOI was eroded by 2 voxels in order to be less susceptible to partial volume effects from the surrounding white matter structure. By contrast the cortex VOI was included in its exact shape, which is subject to partial volume effects, to show the behaviour of these effects with different reconstruction algorithms. Details of the VOIs are provided in table 3.

Real data
A SB dynamic examination centred over the lungs region was used to test the performance of the evaluated algorithms and to compare results against the simulation findings. Approval for the retrospective use of the real patient data was obtained for this study. The data had been collected with approval from a local ethics committee. The original dataset was acquired on a Signa PET/MR, starting at the injection of 177 MBq of FDG tracer to the patient, for a duration of 1 hour. The imaged patient had been diagnosed with a non small cell lung cancer at the left lung. The raw list-mode dataset was retrospectively reprocessed (replayed) to create two new datasets. One dataset using the framing of the simulated SB study and one dataset using the framing of the simulated DWB-1 study (DWB) including temporal gaps. The TOF information available with the original dataset from the Signa PET/MR was not used in the reconstruction process. An IDIF from the ascending aorta was measured on activity image data from 3D reconstruction and used for 4D reconstruction and postreconstruction analysis. The two datasets were reconstructed using the same 3D and 4D dynamic reconstruction algorithms that were used in the simulation study using identical parameters. No respiratory motion correction or gating was applied to the data. Similar to the simulation study, post-reconstruction Patlak analysis at the voxel level was performed with OLS to generate parametric images of K i . VOIs were drawn over the tumour, the tumourʼs background (left lung) and the liver, as shown in figure 2, to compare between reconstructions and against the findings of the simulation study. Using these the contrast to noise ratio (CNR) was estimated according to where SD VOI is the spatial standard deviation of a VOI. Similarly, CNR was calculated using the simulation study data to enable direct comparison with the real data. In this case, the eroded thalamus VOI was used as the target region and the white matter as the background for both contrast and noise estimation.

Comparison between SB and DWB protocol data
The VOI and voxel based metrics comparing 3D and 4D Patlak reconstructions for the SB and DWB-1 protocols are shown in figure 3. For both voxel based and VOI based metrics the 3D reconstruction followed by postreconstruction Patlak fitting using DWB data resulted in higher CoV values, compared to 3D reconstruction of SB data at matched bias. For the first few iterations the 3D reconstructions of both datasets resulted in similar bias values, while further iterations resulted in a wider range of bias values for the DWB data compared to SB data within 40 iterations.
The use of 4D Patlak reconstruction on DWB data produced results with lower CoV on both evaluated metrics and VOIs, compared to the 3D reconstruction of the same data at matched bias, and a shorter range of bias values within 40 iterations. For the VOI metrics, CoV values of the 4D Patlak reconstruction on both evaluated regions approach those of 3D reconstruction of SB data. Furthermore, the 4D Patlak reconstruction of DWB data resulted in eroded thalamus bias values that evolved towards a steady value of positive bias, at approximately iteration 12, after which further iterations resulted in small step changes towards lower bias. For the voxel metrics, similar behaviour is seen on early iterations of 4D Patlak reconstruction on DWB data for the CoV, with values approaching those of 3D reconstruction of SB data. But at further iterations the CoV for the 4D Patlak reconstruction in both VOIs surpasses values from 3D reconstruction on DWB data. On the eroded thalamus this was the case beyond iteration 24, while for the cortex from iteration 22 and beyond. These results show that there is a risk of increasing parametric image noise, greater than that of 3D reconstructions, when the 4D reconstruction is run at high iterations to achieve more favourable and stable mean VOI behaviour.
The use of 4D Patlak reconstruction with SB data showed similar effects on behaviour for CoV and bias on both metrics, compared to 3D reconstruction of SB data, and resulted in the lowest CoV values for these comparisons.

Comparison between 4D dynamic reconstructions
The VOI and voxel based metrics are shown in figure 4 for comparison of 4D Patlak and 4D Spectral reconstructions of DWB and SB data. On both metrics and for both regions the use of Spectral reconstruction with 6 basis functions provided the lowest CoV values at matched bias compared to other 4D reconstructions, on DWB and SB data respectively. However 4D Spectral reconstruction with 6 basis also provided the highest respective bias values in the eroded thalamus. On the cortex the difference on bias metrics was relatively small between all 4D reconstructions of respective DWB and SB data.
The 4D Spectral reconstruction with 9 and 17 basis functions resulted in similar bias and CoV values at both regions on DWB data. Their use resulted in lower CoV compared to 4D Patlak reconstruction and 3D reconstruction of SB data, but higher compared to 4D Spectral using 6 basis. Nonetheless, at the eroded thalamus use of 9 and 17 basis provided improved bias values at matched CoV when compared to the use of 6 basis, closer to bias values from 4D Patlak reconstruction. On SB data, a greater separation in rms Bias is seen between 4D Spectral reconstructions with 9 and 17 basis functions. When the 4D Spectral reconstructions of  DWB data were provided with the same frame data as the 4D Patlak reconstructions (4 frames with t > t ss ), instead of all data (8 frames for DWB-1), it resulted in a noticeable increase of the CoV values with very close noise versus bias trade-offs between 6 and 9 basis. Their trade-off curves got closer to the one of the 4D Patlak reconstruction, with lower bias values on both metrics but higher rms CoV compared to the 4D Patlak reconstruction.

Comparing between nested optimizations in 4D reconstruction
Results of 4D Spectral and 4D Patlak reconstructions of DWB data using MLEM and NNLS nested optimization are shown in figure 5. An additional figure comparing the two nested optimization methods on SB data is provided in the supplementary material. A clear difference in behaviour is seen going from MLEM to NNLS from early iterations, with 4D reconstructions using nested NNLS optimization resulting in higher CoV values at matched bias compared to the respective 4D reconstruction using nested MLEM (with 20 nested sub-iterations). At the same time, NNLS nested optimization often resulted in a slight reduction in bias at matched CoV values.
No difference was seen in convergence properties such as convergence speed between the two nested optimization options. Nevertheless, the use of a single run of NNLS optimization in each nested loop, instead of 20 nested MLEM iterations, resulted in a notable reduction of overall reconstruction time. The average reconstruction time of DWB-1 data using the two methods on a computer with a 16-core 2.20 GHz processor and 96 GB of RAM memory is shown in table 5.

Comparison between DWB protocols
A comparison of 4D Patlak and 4D Spectral reconstructions between the three simulated DWB protocols is made in figure 6. For VOI based metrics on both regions and 4D reconstructions, data from all three DWB protocols resulted in bias values within the expected range of the underlying systematic bias in Patlak analysis. These biases are caused by the sampling differences of the DWB protocols, which result in the values shown in table 4 for noiseless TAC data. Differences in CoV at matched bias values are more profound in VOI metrics of the cortex region, where data from protocol DWB-2 resulted in the lowest CoV values and DWB-1 and DWB-3 data resulted in closer CoV values. This level of reduction in CoV was not seen on the eroded thalamus, neither on VOI or voxel based metrics. In the eroded thalamus the differences of DWB protocols on CoV at matched bias were smaller and their ordering was mixed between the two 4D reconstructions.

Comparison with real data
The comparison of reconstructions using a real FDG dataset, reprocessed and reconstructed with the SB and DWB-1 protocol framings, is made in figure 7. Results on the real data showed similar evolution of CNR with increasing iterations for 4D and 3D reconstructions respectively. The 4D Spectral reconstruction using 6 basis functions provided the highest CNR values throughout all iterations, followed by the 4D Spectral reconstruction using 9 basis functions and 4D Patlak. Overall, CNR of all 4D reconstructions of DWB data was higher than that of 3D reconstruction of SB data and 3D reconstruction of DWB data. The liver SD versus tumour VOI average trade-off curves showed close behaviour between 4D reconstructions, resulting to SD values in the first 15-17 iterations of 4D reconstructions which were lower compared to 3D reconstruction of DWB data at matched tumour mean values. Compared to 3D reconstruction of SB data, all 4D reconstructions of DWB data resulted in higher SD values (for matched tumour mean values where comparison is possible). Parametric K i images from 3D and 4D reconstructions at iterations with matched liver SD values of approximately 3.1·10 −3 min −1 are shown in figure 8. The K i images show the degradation of the tumourʼs and heartʼs high-intensity definition in the transition from SB to DWB with 3D reconstruction based parametric imaging. The 4D reconstructions of DWB data improve on these definitions with 4D Spectral reconstruction providing a similar definition on the heartʼs structure to that of SB data with 3D reconstruction. The remaining differences in the heartʼs definition can be also attributed to differences in convergence between the presented iterations of each respective reconstruction.
For comparison against the simulation study, metrics of CNR and SD versus VOI average trade-off curves from simulation DWB-1 are shown in figure 9. Similar to the real data, values of CNR are highest for 4D Spectral reconstruction using 6 basis functions, followed by 4D Spectral reconstruction using 9 basis functions. Contrary to the real data, 4D Patlak reconstruction provided relatively lower CNR values and closer to the values of 3D reconstruction from SB data. The evolution of CNR with increasing iterations was similar to the real data, with maximum values attained at approximately 3 to 5 iterations for all reconstructions. Furthermore, the simulation data SD versus VOI average trade-off curves showed a larger separation of 4D Spectral reconstructions, compared to the real data, with the reconstruction of 6 basis providing lower SD values compared to 3D and 4D Patlak reconstruction of DWB data, but at higher bias values from other 4D reconstructions. Finally, a difference in the bias paths with increasing iteration between the real and simulated data can be seen that is most prominent for 4D reconstructions. These differences can be attributed to differences in the VOIs used to generate these curves, between the real data and the simulated data, as well as differences in their surrounding structures. Nevertheless, in both cases the curves show similar behaviour on bias and SD step changes with increasing iteration, with bias step changes reducing noticeably after approximately 10 iterations while further iterations result in small step changes of bias and SD. Parametric K i images of 3D and selected 4D reconstructions at matched rms CoV values, of approximately 32% as measured at the eroded thalamus VOI, are shown in figure 10 for a single noise replicate along with images of mean bias over noise replicates. Additional sample K i images for all evaluated reconstruction algorithms are provided in the supplementary material. The single replicate images show that structures of the thalamus seen in 3D reconstruction of SB data are better resolved in DWB when using the 4D Spectral reconstruction. The images of bias show similar behaviour in the thalamus over reconstructions and demonstrate the partial volume effects at the cortex region.

Discussion
Our simulation study shows that the dynamic reconstruction of DWB FDG data resulted in a substantial reduction of Patlak K i image noise and more favourable convergence behaviour, compared to 3D reconstruction based parametric imaging. These results, limited to a single level of noise, are in agreement with the findings of Karakatsanis et al (2016). Moreover, we directly compared against an SB dynamic protocol, processed with 3D reconstruction, and showed that comparable values of parametric image noise and bias can be achieved with DWB protocols by the use of a dynamic reconstruction.
In this study, a SB location was used in the evaluations of reconstructions and comparisons between SB and DWB protocols. Although this study is limited to a single axial location, we anticipate the results to hold for FDG Patlak K i parametric imaging of other bed positions on both S&S and uni-directional CBM DWB protocols. We assume this because all bed positions of these DWB protocols are sampled uniformly with time and differ only by a time offset in their sampling. This fact, along with the fact that FDG Patlak K i parametric imaging relies on relatively slow dynamic changes after steady state conditions have been reached, make the assumption of generalisation to other bed positions reasonable. By contrast, bi-directional CBM introduces non-uniform  frame sampling whose characteristics vary significantly per bed position. In this case we cannot safely generalise the findings of this SB study to other bed positions of bi-directional CBM, which necessitate simulation of all bed positions. Furthermore, the results cannot be extended to parametric imaging of other kinetic parameters which might be more sensitive to fast dynamic changes and hence timing differences between bed positions and differences due to acquisition temporal gaps. For this study we have assumed that individual bed positions of the DWB studies are to be reconstructed independently and the result parametric K i images to be merged after reconstruction into a WB K i map. This method initially suggested by Karakatsanis et al (2016) combines the reconstructed images, including the overlapping information, in parametric space. Although this technique is practical in terms of implementation,   it may not be making best use of acquired statistics at the overlap region. Ideally, the combined bed statistics and timing information at the overlapping regions should be considered directly in the reconstruction process, either using slice dependent timing information for CBM DWB reconstruction (Hu et al 2020) or with synchronous use of all bed dynamic information from S&S DWB within a direct multi-bed DWB reconstruction (Chalampalakis et al 2020). Further research into the behaviour of the overlapping regions within reconstruction and into methods for use of multi-bed DWB data within the reconstruction process is needed to evaluate these techniques for parametric imaging over the entire WB FOV. The simulation study is conducted for a number of noise realisations, that is limited by the computational requirements of the reconstructions. This limitation results to some uncertainty in the evaluation metrics that is not explored in this study. Furthermore, differences in the number of noise realisations and noise properties of the different simulated protocols and reconstruction algorithms can result in different levels of uncertainly, most notably for rms Bias which had shown to be sensitive to remaining noise on Bias images over replicates. Nevertheless, comparisons against the use of fewer replicates in our results showed that these effects did not affect the qualitative findings of our results on the ranking of different reconstructions between SB and DWB protocol data. Additional work is required to estimate the uncertainty of the used metrics and a practical number  of replicates at which this uncertainly is minimised for all simulated protocols and reconstructions over all evaluated OSEM iterations.
The choice of the iteration number to terminate a 4D reconstruction algorithm is not evident. For a VOIbased analysis, the convergence of the mean K i value in the cortex or the eroded thalamus was not seen in the range of the 40 evaluated OSEM iterations, in particular for a 3D reconstruction. This behaviour was also observed for a 4D reconstruction algorithm, but to a lesser extend. A high number of iterations of 4D reconstruction algorithms provided more stable VOI mean values, but at risk of resulting to higher parametric image noise than that of a 3D reconstructions on the same DWB data. The results obtained with one real dataset showed similar behaviour, with mean VOI values continuing to slightly increase even after 40 OSEM iterations and with 4D reconstruction parametric image noise surpassing that of 3D reconstruction at late iterations. This example illustrates that the relative aspects of 4D to 3D comparisons with simulated data for the tested DWB and SB protocols have the capacity to translate to studies with different levels of noise. Overall, the risks of excessive parametric image noise and under-converged K i values will be lesser for 4D based reconstruction methods than for a 3D reconstruction, which demonstrated considerably more instability with increasing iterations. To ensure convergence of the K i values while suppressing the increase of parametric image noise, further regularisation techniques can be used with methods such as 4D MAP reconstruction (Wang et al 2008, Reader andVerhaeghe 2014) or kernel 4D dynamic reconstruction (Novosad andReader 2016, Gong et al 2018).
Our nested optimization tests using NNLS instead of multiple MLEM sub-iterations did not provide any differences in the acceleration of the convergence and showed comparable behaviour to a previous study on the use of NNLS with the spectral model (Matthews et al 2010). NNLS did provide computing acceleration by a factor of approximately two for our datasets, but resulted in an increase of the parametric image noise compared to MLEM sub-iterations for similar bias characteristics. Equivalent or higher acceleration could be potentially achieved if the nested MLEM optimization was conducted in graphical processing units instead of the CPUs.
In this work, we evaluated the use of an indirect dynamic reconstruction method based on a generic 4D Spectral reconstruction algorithm followed by a post-reconstruction Patlak model fitting. The genericity of the spectral model allows for flexibility in modelling dynamic processes that do not necessarily fall under the idealised behaviour of the kinetic model of interest. In this simulation study, we were limited to irreversible FDG kinetics that can be sufficiently described by the Patlak model. In this case, 4D Spectral reconstruction making use of the full dynamic data outperformed the direct Patlak reconstruction in terms of parametric image noise, while maintaining similar bias behaviour. The benefit of the 4D spectral reconstruction was less obvious when fewer frames were used in reconstructions using t > t ss , indicating that its favourable behaviour was mostly due to the use of more temporal frames than the Patlak reconstruction. In real FDG studies, it can be desirable to account for reversible FDG kinetics and reduce the bias of the estimated macro-parameters arising from poorly modelled kinetics. The spectral model can allow for more complex compartmental modelling with no strong prior knowledge or enforcement of a specific model. As such it can account for more complex kinetic behaviours, including reversibility of tracer, in the reconstruction process and allow for post-reconstruction exploratory modelling to identify the best model to describe and present the data. Moreover, for DWB studies where not all body regions and organs will necessarily be adequately described by a single dynamic model of interest, the proposed indirect method can allow for the assignment of different kinetic models in different regions of the body to ensure appropriate representation of the dynamic tracer behaviour. Depending on the availability of early frame data, the fitted spectral model can be used to directly estimate K 1 (Meikle et al 1998, Matthews et al 2010, while post-reconstruction micro-parameter estimation could be performed for potential uses in clinical applications (Novosad andReader 2016, Zaker et al 2020) and indirectly take advantage of the 4D reconstruction temporal regularisation. An important parameter to configure for the spectral model is the number of basis functions. Contrary to post-reconstruction spectral analysis where hundreds of basis functions are used to finely sample the space of kinetic exchange rates, a smaller number of basis is desirable in reconstruction to favour reduced image noise. In some cases of our findings the lowest number of basis functions used (6 basis) resulted in higher bias values which indicates less than adequate modelling of the underlying kinetics, compared to reconstructions with more basis functions and to 4D Patlak reconstruction. However, this was not the case when fewer frames were used in reconstructions using t > t ss data. These findings indicate that the selection of the number of basis functions is important not only for controlling the produced image noise but also for controlling bias by adequately modelling the kinetics behaviour in reconstruction. For the higher numbers of basis functions, with 9 and 17, almost identical behaviour was seen on the DWB-1 dataset (of 8 frames). Overall on the choice of the number of basis functions, results indicate a greater risk in image bias when using a too small number of basis functions, and a lesser risk in image noise when using more basis functions than strictly needed to properly model the underlying kinetics. In any case the number of basis needs to be tuned for every DWB protocol, depending on the number of frames within the dataset and the range of underlying kinetics as well as the level of noise in the PET data.
In this study, the investigation between S&S and CBM DWB protocols was limited to aspects of sampling frequency and uniformity within the total examination time. Our results showed small differences in parametric image bias but a noticeable reduction in parametric image noise when utilising CBM acquisition with uniform sampling. Overall differences were inline with previous findings in comparison of S&S and CBM DWB imaging on a real data study using different metrics (Karakatsanis et al 2015). Furthermore, beyond the aspects of reduced acquisition delays and higher sampling frequency, CBM acquisition has other desirable properties for DWB acquisitions as outlined previously (Karakatsanis et al 2015). The most important aspect is the result uniform axial sensitivity profile at any choice of acquisition speed. That can be of importance in DWB parametric imaging where multiple regions of interest are expected in the effective FOV. In our study we have not considered this aspect for the CBM protocols and we did not examine regions in the overlap range of the S&S protocol. But the observed improvements related to reduced delays in acquisition coupled with uniformity of axial sensitivity favour the choice of CBM over S&S protocols. We investigated further potential reductions in system delays by allowing for non-uniform frame sampling using bi-directional CBM. In that case we did not see the same effects as in the transition from S&S to CBM. But our results on bi-directional CBM are limited to the specific framing of the evaluated protocol design which offered more total frames but resulted in fewer total counts compared to the other protocols. As discussed above, the findings of our study are limited to a single axial bed location. Since non-uniform sampling of the bi-directional CBM protocol is strongly affected by axial location we cannot generalise the results of the comparison between the two CBM protocols over the WB. Additional tests are required on the exploitation of the flexibility offered by bi-directional CBM to assess other potential benefits against uni-directional CBM.

Conclusion
4D dynamic reconstruction is necessary in DWB parametric imaging to achieve accurate and stable quantification. For FDG Patlak K i parametric imaging we have shown results of direct Patlak dynamic reconstruction with noise and bias values that were comparable to 3D reconstruction based parametric imaging from SB dynamic studies. In this work, we proposed the use of an indirect method for DWB parametric imaging, based on the spectral analysis model. This more flexible approach allows for complex kinetic modelling to be used during reconstruction for temporal regularisation, with minimal assumptions on the underlying kinetics. In Patlak K i parametric imaging this method outperformed the direct Patlak approach, by making use of all the acquired data for temporal regularisation from which post-reconstruction parametric imaging benefitted by further reduction of noise compared to the Patlak approach. Furthermore, the spectral model approach can be used for more complex post-reconstruction modelling, for example in parametric imaging of FDG microparameters. Finally, we investigated the impact of various acquisition modes (for CBM and S&S) resulting in different temporal sampling of the data. Benefits of reduced delays and increased acquisition statistics were partially seen in reduced parametric image noise for the CBM protocol with uni-directional axial sampling. By contrast, CBM using bi-directional motion resulted in parametric image noise levels that were similar to the S&S protocol. Further investigation is required to assess the potential benefits from bi-directional CBM against unidirectional CBM and the effects of non-uniform sampling over the entire FOV of the DWB protocols.
Overall, the use of 4D dynamic reconstruction for DWB parametric imaging offers desirable properties that enable the transition from SB dynamic studies and common 3D reconstruction parametric imaging practices without loss of image quality and with additional benefits for accuracy of parametric images. Potential applications of DWB parametric imaging are expected to rely on the quantification of images and so there should be no compromise between parametric image accuracy and image noise. Our results showed that 4D reconstructions need to be sufficiently iterated to ensure accurate quantification, with potential for improvement in maintaining low parametric image noise by use of additional regularisation methods.