NiftyPAD - Novel Python Package for Quantitative Analysis of Dynamic PET Data

Current PET datasets are becoming larger, thereby increasing the demand for fast and reproducible processing pipelines. This paper presents a freely available, open source, Python-based software package called NiftyPAD, for versatile analyses of static, full or dual-time window dynamic brain PET data. The key novelties of NiftyPAD are the analyses of dual-time window scans with reference input processing, pharmacokinetic modelling with shortened PET acquisitions through the incorporation of arterial spin labelling (ASL)-derived relative perfusion measures, as well as optional PET data-based motion correction. Results obtained with NiftyPAD were compared with the well-established software packages PPET and QModeling for a range of kinetic models. Clinical data from eight subjects scanned with four different amyloid tracers were used to validate the computational performance. NiftyPAD achieved \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2>0.999$$\end{document}R2>0.999 correlation with PPET, with absolute difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{-2}$$\end{document}∼10-2 for linearised Logan and MRTM2 methods, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2>0.999999$$\end{document}R2>0.999999 correlation with QModeling, with absolute difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{-4}$$\end{document}∼10-4 for basis function based SRTM and SRTM2 models. For the recently published SRTM ASL method, which is unavailable in existing software packages, high correlations with negligible bias were observed with the full scan SRTM in terms of non-displaceable binding potential (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2=0.96$$\end{document}R2=0.96), indicating reliable model implementation in NiftyPAD. Together, these findings illustrate that NiftyPAD is versatile, flexible, and produces comparable results with established software packages for quantification of dynamic PET data. It is freely available (https://github.com/AMYPAD/NiftyPAD), and allows for multi-platform usage. The modular setup makes adding new functionalities easy, and the package is lightweight with minimal dependencies, making it easy to use and integrate into existing processing pipelines. Supplementary Information The online version contains supplementary material available at 10.1007/s12021-022-09616-0.


Introduction
As the biomedical research community is moving towards earlier disease detection, it has become increasingly relevant to capture subtle changes in biological processes Guo et al. (2020), Cohen and Klunk (2014). A technique that allows for both visualising and quantifying these processes in vivo is Positron Emission Tomography (PET) Guo et al. (2020), Cohen and Klunk (2014). PET scans can be used to obtain highly accurate, quantitative measurements, provided that a suitable scanning protocol is used. To this end, dynamic or more recently established dual-time window scanning protocols may be required, because of their ability to take into account physiological changes Bullich et al. (2018), Lammertsma (2017), Heeman et al. (2019), van Berckel et al. (2013). These protocols require longer scanning time, complex processing and pharmacokinetic analyses.

3
As the size of PET datasets continue to increase, so will the demand for efficient processing pipelines that provide high throughput Rabinovici et al. (2021), Ferrucci (2018), LaMontagne et al. (2019).
In addition to existing commercial software packages, such as PMOD (PMOD Technologies, Zurich, Switzerland), publicly available software packages for pharmacokinetic analysis of dynamic brain PET data have been introduced, including COMKAT Muzic and Cornelius (2001), PPET Boellaard et al. (2006), MIAKAT Gunn et al. (2016) (previously available), QModeling López-González et al. (2019), Magia Karjalainen et al. (2020b) and kinfitr Matheson (2019). PPET is written in IDL (L3Harris Geospatial Solutions, Inc.). COMKAT, MIAKAT, QModeling and Magia are written and run in MATLAB (The MathWorks, Inc.). Both IDL and MATLAB are not freely available. kinfitr is written in the freely available R language. In this work we present NiftyPAD, where PAD stands for package for quantitative analysis of dynamic PET data. It is written in Python, which is completely open source, and designed to increase the capacity of a pharmacokinetic modelling software package. An important feature NiftyPAD provides is the ability to analyse PET data acquired in a dual-time window protocol Heeman et al. (2019), to support the field given the growing number of studies that acquire early PET data after tracer injection in addition to the late PET scan LopesAlves et al. (2020), Cecchin et al. (2017), Son et al. (2020). A dual-time window protocol is useful for quantitative studies using tracers with slow kinetics because of their long dynamic acquisition times (up to 130 min). To decrease the scanning burden for the participant, a break is inserted in the middle of the scan during which the participant can leave the scanner and rest. The break also allows for interleaved acquisition, in which the scan of a second participant can be initiated during the break of the first participant. To allow for kinetic modelling, these dual-time window or dual-phase data require appropriate interpolation of the reference tissue curve, which has been rarely implemented in existing software packages Funck et al. (2018), Karjalainen et al. (2020a). Another novel feature of Nifty-PAD is pharmacokinetic modelling through incorporation of ASL-derived relative perfusion measures for simultaneous PET-MR scans Scott et al. (2019). This feature circumvents the need for an early PET scan by combining a simultaneously acquired ASL scan with a static PET scan to allow for kinetic modelling.
The growing number of dynamic datasets and their complex analyses, demonstrate the need for software tools that facilitate straightforward and automated analyses of PET data. In addition, combining these aspects in an opensource software tool will improve the reproducibility of the results by allowing others to use the same algorithm, whilst maintaining the possibility of tailored software solutions. Therefore, NiftyPAD was designed to support several important features which are not available in other existing software packages for kinetic modeling: • Freely available (Python-based) and open-source, allowing for full transparency • Analysis of static, dynamic and dual-time window PET data • Pharmacokinetic modelling with the incorporation of arterial spin labelling (ASL)-derived relative perfusion measures for simultaneous PET-MR scans • Motion correction, through a built-in kinetics based realignment function, and options for excluding PET frames with large motion-misaligned attenuation correction in kinetic analysis The present paper discusses features that are incorporated in the NiftyPAD software package, with focus on the modelling aspects. Implementation of NiftyPAD was evaluated by comparing its numerical results of estimating kinetic parameters with those obtained by using the established software packages PPET Boellaard et al. (2006) and QModeling López-González et al. (2019). Clinical PET data from various amyloid tracers were used for the evaluation. A selection of core models were assessed given the model availability in the existing software packages. Figure 1 summarises the features of NiftyPAD and its workflow. NiftyPAD provides a group of reference-based kinetic models to generate parametric images or regional kinetic parameters. For quantification, NiftyPAD requires the user to provide a PET scan, data for obtaining a reference tissue time-activity curve (TAC), a region of interest (ROI) template in case of regional analyses, and modelling settings for the selected kinetic models. More specifically, reference tissue input processing is implemented for interpolating missing reference tissue data points of PET data acquired according to a dual-time window protocol and to improve accuracy of the reference tissue curve in case of noise or motion. Other implemented features include support for weighting the temporal data points and options for motion correction.

User Input
As shown in Fig. 1, NiftyPAD requires a user to provide: 1) PET imaging data (static, dynamic or dual-time window) for parametric analysis, in combination with frame start and end times in seconds, 2) a ROI template (parcellation) for extracting regional TACs, or user-pre-defined TACs for regional analyses, 3) reference tissue input data, provided directly by the user or using NiftyPAD to extract for a named reference region, and 4) model settings to select the quantification methods and configurations of model parameters.
With regard to the PET scan, most commonly used imaging data formats are supported, such as dicom, nifti, analyze and ecat (see the NiBabel Python library for a full list of supported file types, https:// doi. org/ 10. 5281/ zenodo. 42955 21). In addition, fitting settings (e.g. the number of basis functions and the applied boundary conditions) can be specified.

Input Processing
NiftyPAD allows for processing of the reference tissue input data. This can be used to improve accuracy in case of noise or motion, or to interpolate missing data in the specific case of dual-time window scans where the dynamic data are acquired in separate sessions.
The following four methods for reference input processing are provided in NiftyPAD:

Linear Interpolation
The reference input is interpolated using linear functions, this is usually suitable for standard scan protocols without dual-time window.

Cubic Interpolation
The reference input is interpolated using cubic functions, when linear functions are not suitable for interpolating between two time points in the reference data.

Exponential Interpolation
The decay part of the reference input is interpolated using an exponential function: where C R (t) is the tracer concentration in the reference region, a 0 , a 1 , b 1 are the unknown parameters to be determined by the reference input data, t s and t e are the starting and ending times for the exponential decay phase chosen by the user.

Feng's Plasma Input Model + 1TC
In addition to cubic and exponential interpolations, NiftyPAD also provides a compartmental model-based method for reference input processing. The tracer behaviour in the reference tissue is described by a one-tissue-compartment (1TC) model, with the operational function where C R (t) and C p (t) are tracer concentrations in the reference tissue and arterial plasma, respectively, a 3 and b 3 are the unknown parameters to be determined by the reference input data. The plasma input function C p (t) is described by Feng's plasma input model Feng et al. (1994) as where a 0 , a 1 , a 2 , b 0 , b 1 , b 2 are unknown parameters determined together with a 3 and b 3 by fitting Equation 2 to the reference input data.
Examples of using cubic and exponential interpolations, and the Feng's input model + 1TC method on clinical data are shown in "Comparison with Established Softwares".

Pharmacokinetic Quantification
The current version of NiftyPAD focuses on a group of reference tissue input based methods that do not require invasive arterial blood sampling or processing, which are widely applied to a range of neuro studies LaMontagne et al. (2019), Golla et al. (2016), Mertens et al. (2020).
The code for all methods is available at https:// github. com/ AMYPAD/ Nifty PAD/ blob/ master/ nifty pad/ models. py. (1) The kinetic parameters of interest are estimated by numerical optimisation, using weighted (unweighted if weights are not given) least squares as the objective function. For quantification of both region of interest (ROI) and voxel data, the following linear and nonlinear models are implemented in NiftyPAD:

SRTM and SRTM2
The simplified reference tissue model SRTM Lammertsma and Hume (1996), and SRTM2 Wu and Carson (2002) with pre-defined tissue-to-plasma clearance k ′ 2 of the reference tissue are implemented with the operational functions and respectively. Where C T (t) and C R (t) are tracer concentrations in the target and reference tissue compartment, respectively, R 1 is the relative tracer delivery parameter ( K 1 /K ′ 1 ), k 2 is the rate constant for transfer from free to plasma compartment, k 2a is the apparent rate constant for transfer from the specific compartment to plasma and t is the time in minutes. The kinetic parameters of interest R 1 , k 2 , and k 2a = k 2 ∕(1 + BP ND ) are estimated by nonlinear optimisation. The non-displaceable binding potential, BP ND , is then derived by BP ND = k 2 ∕k 2a − 1 for SRTM, and

SRTM Basis and SRTM2 Basis
The linearisation of SRTM with basis functions is implemented based on Gunn et al. (1997), with the operational function where B i are the number of basis functions and the basis functions are calculated as The pre-defined k ′ 2 version SRTM2 with basis functions is implemented based on Wu and Carson (2002), with the operational function where the basis functions are calculated as Details on solving the linearised problem and calculations of the kinetic parameters of interest R 1 , k 2 and BP ND can be found in Gunn et al. (1997) and Wu and Carson (2002).
SRTM ASL NiftyPAD features a recently developed model SRTM ASL Scott et al. (2019) for analysing PET data acquired by a simultaneous PET-MR scanner where arterial spin labelling (ASL) is available to provide the perfusion information and derive the relative influx rate R 1 . In SRTM ASL, R 1 is derived from ASL MRI data by using Equations 3 and 4 in Scott et al. (2019). With pre-determined R 1 , the SRTM basis model can be re-written as where C d T (t) is a dummy variable determined before fitting. is then solved with the pre-defined basis functions B i (t) , which follows the same definition as in Eq. 7. The kinetic parameters of interest k 2 and BP ND are then derived as in Scott et al. (2019).

Logan and Logan2
The graphical analysis methods with reference input Logan plot, and Logan2 with pre-defined k ′ 2 are both implemented based on Logan et al. (1996), by and respectively, where int ′ stands for the intercept in the linear regression, DVR for the distribution volume ratio, and binding potential BP ND = DVR − 1.

MRTM and MRTM2
The Ichise's Multilinear Reference Tissue Model MRTM and MRTM2 are both implemented based on Ichise et al. (2003). The operational equations for MRTM and MRTM2 with pre-defined k ′ 2 are and (9) respectively, where the coefficients 1 , 2 and 3 are determined by linear regression. The kinetic parameter of interest BP ND is then calculated by BP ND = −( 1 ∕ 2 + 1).

SUVr
Standardised uptake value ratio (SUVr), the most commonly used semi-quantitative method is calculated as the ratio of target to reference tissue activity for a pre-defined time window.

Weighting Schemes
To improve the accuracy of kinetic analysis, weighting schemes are commonly used when fitting a kinetic model to dynamic PET data. Introduction of weighting schemes can improve model fitting especially with noisy PET data, or data acquired using a dual-time window protocol. Nifty-PAD therefore provides a number of weighting schemes to facilitate optimal fitting Yaqub et al. (2006), and also allows the user to input pre-determined weights.

Motion Correction
NiftyPAD features an integrated motion correction module, which can be used with PET imaging data. The module consists of two features: 1) exclude PET frames that suffer from severe motion and cause a mismatch in attenuation correction, and 2) correct for motion between PET frames by performing a groupwise registration that minimises the errors in the fits of the pharmacokinetic modelling. This approach is implemented based on Jiao et al. (2014). Note that for dual-time windows scans, this motion correction approach can be applied before aligning the early and late frames. However the co-registration of the early and late frames using external structural images can greatly improve the motion estimation and correction.

Evaluation
Performance of the main functions provided in NiftyPAD was evaluated using clinical data. Here, we present the results of 1) processing the reference tissue input data for dual-time window scans with various methods, 2) comparing the kinetic parameters computed by NiftyPAD and established software packages, 3) applying the ASL SRTM method to a dynamic dataset, and 4) generating parametric images.

Performance on Input Processing
For PET scans acquired using dual-time window protocols, appropriate reference tissue input interpolation of missing data-points between the two parts of a scan is an essential step before performing kinetic analysis, as linear interpolation or simply concatenating the two parts of the scan does not follow the kinetics in the reference data. Figure 2 shows an example of the interpolated reference tissue TACs obtained by using cubic interpolation, exponential interpolation and Feng's plasma input + 1TC model, which are available in the Nifty-PAD software package and described in "Input Processing".

Comparison with Established Softwares
As the core function of NiftyPAD is performing kinetic analysis on dynamic PET data, we evaluated the implementation of the kinetic models SRTM basis, SRTM2 basis, Logan and MRTM2, by comparing the kinetic parameters computed by NiftyPAD with those computed by the two established softwares, PPET Boellaard et al. (2006) and QModeling López-González et al. (2019). Given the model availability and implementation similarity, PPET was used for the graphical models Logan and MRTM2, and QModeling was used for the basis function based methods SRTM basis and SRTM2 basis. We did not use PPET for the SRTM models as the computation of the basis functions and the solution of the linear problem were implemented differently in PPET.
Clinical data from eight subjects scanned with four different amyloid tracers were used for this evaluation. The subjects include A -positives and A -negatives. Four were from the AMYPAD Prognostic and Natural History Study (PNHS) in which two were scanned with [ 18 F]flutemetamol and two with Feng's plasma input + 1TC time (sec) Activity concentration (Bq/mL) C Fig. 2 Interpolation of reference tissue TAC from one [ 18 F]florbetaben scan with dual-time window acquisition, using A cubic interpolation, B exponential interpolation, and C Feng's plasma input + 1TC model. Reference tissue input processing serves as an essential first step before kinetic analysis of dual-time window PET scans [ 11 C]PiB dataset ), Ossenkoppele et al. (2012. All participants provided written informed consent before participating in the study. The local Medical Ethics Review Committee approved the study protocol and scans were acquired at the Amsterdam UMC, location VUmc. The full demographic information of the subjects can be found in Table 1 in Supplementary Materials. Detailed acquisition and processing of these scans are described in LopesAlves et al. (2020), , Verfaillie et al. (2021). Notably, the [ 18 F]flutemetamol and [ 18 F]florbetaben scans followed the dual-time window protocol. These scans were not included in the comparison of NiftyPAD and QModeling, as the latter one does not support data acquired according to a dual-time window protocol.
For each scan, brain parcellation was applied to extract the time-activity curves (TACs) from anatomical regions of interest. Thirty regional TACs were selected from each subject with a range of volume sizes and binding levels. These regional TACs were used to quantitatively compare the kinetic parameter values estimated by NiftyPAD and the established softwares. Correlation and Bland-Altman plots of the kinetic parameters obtained by using Nifty-PAD and the reference software are shown in Figs. 3,4,5,6,7 and 8. Comparing with PPET for the graphical models, the absolute differences were in the order of 0 − 10 −2 for BP ND . This can be due to the difference in precision between the programming languages, as PPET is written in IDL and has by default single precision, whereas NiftyPAD is written in Python and has by default double precision. Also the implementation of the linear regression used in the graphical models can results in differences between these two packages. Comparing with QModeling for the basis function based SRTM models, the absolute differences were in the order of 0-10 −4 for BP ND and R 1 . We did not find differences in implementation between NiftyPAD and QModeling of these two SRTM models. The differences are considered to be caused by the difference in handling the floating point calculations between Python and MAT-LAB which QModeling uses.

Application of SRTM ASL to Regional Data
In order to show the effectiveness of the SRTM ASL model for estimating BP ND from a shorter acquisition, the model was applied to TACs from one dynamic [ 18 F]flutemetamol data set that had been acquired using a dual-time window protocol with an acquisition phase from 0-30 minutes and 90-110 minutes post injection (A -negative, male, 63.0 years old). The scan was acquired using a Siemens mMR Biograph and, simultaneously, a T1 MPRAGE Sagittal MRI was acquired. The T1-weighted MR image was segmented into 11 different regions (white matter, cortex, cerebellar white  Data points correspond to different brain regions from each subject. The dashed lines are the line of identity on the left, and the mean difference on the right matter, thalamus, caudate, putamen, pallidum, brain stem, hippocampus, amygdala and accumbens) using Freesurfer version 6.0.0 Fischl (2012), and TAC data were extracted from the PET scan for each region. Note that for the purpose of code validation, we did not use the R 1 derived from ASL MRI as in Scott et al. (2019), for the reason that the noise in ASL will propagate into the derived R 1 and the evaluation of the SRTM ASL implementation will be confounded. Instead, SRTM basis was applied to the full PET data to determine R 1 and then, SRTM ASL was applied to the 90-110 minutes data with the determined R 1 as an input parameter varied per region. The BP ND values from SRTM ASL and SRTM basis were compared using linear regression analyses.
Excellent correlation and negligible bias was observed between the SRTM basis and SRTM ASL methods (Fig. 9). The BP ND values from SRTM ASL and SRTM basis are not identical due to the difference in the length of scan involved in analysis, in this case 0-110 minutes for SRTM basis and 90-110 minutes for SRTM ASL.

Parametric Image Generation
One scan of dynamic [ 11 C]PiB data of 90 minutes duration was used to illustrate the parametric image generation. The scan belongs to a patient with Alzheimer's disease dementia (A -positive, male, 60 years old) . Pre-processing was done as described previously  and Logan reference model was used to generate parametric images with cerebellar grey matter as reference tissue. BP ND images derived using NiftyPAD and PPET are shown in Fig. 10.

Discussion
In this paper we present NiftyPAD, a software package for performing kinetic analysis of dynamic PET images from clinical studies. It is written in Python, is open source and freely available, thus providing full transparency. The software is multi-platform, stand-alone, has minimal dependencies, all of which are standard Python packages, and is therefore easy to be extended and integrated into any PET data processing and analysis pipeline. The key novelties of NiftyPAD includes the analyses of dual-time window scans with reference input processing, and pharmacokinetic modelling with shortened PET acquisitions through the incorporation of arterial spin labelling (ASL)-derived relative perfusion measures. Results produced by NiftyPAD were compared with those produced by two established software packages PPET and QModeling. NiftyPAD achieved reliable results and we noticed that differences in the programming languages used, and especially differences in the implementation of the numerical operations underlying the kinetic models can result in differences in the final outcome measures. For data consistency, we recommend to use the same kinetic analysis software to reduce the variability when comparing kinetic parameters.
Compared with currently available pharmacokinetic modelling software packages, a major novelty of NiftyPAD is the ability to analyse PET data acquired in a dual-time window protocol. This feature is of great importance to the field, given the growing number of studies that acquire early PET data at tracer injection in addition to a late, static scan. More specifically, to allow for kinetic modelling, these dualtime window data require interpolation of the missing data points of the reference tissue curve, which can be done via reference input processing. Of course, as each tracer has distinct kinetics, the optimal interpolation scheme requires validation per tracer. Another distinct feature of NiftyPAD is pharmacokinetic modelling through incorporation of ASL-derived relative perfusion measures for simultaneous PET-MR scans. This feature circumvents the need for an early PET scan by combining relative perfusion measures from an ASL scan with a late PET scan to allow for kinetic modelling. However, as this model has only been applied to [ 18 F]florbetapir data in a proof-of-concept study, its applicability to other (amyloid) tracers is still warranted Scott et al. (2019). Finally, given that patient motion occurs both in clinical and research practise, NiftyPAD offers two distinct features to limit the effect of motion on quantification.
The current version of the NiftyPAD software package can be extended in several directions. For example, plasma input-based models, which are of particular interest for tracer-validation studies can be added to broaden the analysis capacity. At present, documentation, demos/tutorials and graphical user interface are under development in order to facilitate the use of NiftyPAD in routine clinical studies. Lastly, future versions will focus more on quality assessment of the generated output, which is currently also under development.

Conclusion
NiftyPAD is a freely available, open-source software program for quantitative PET analyses. It features unique modules that allow for analysing data acquired in a dual-time window protocol, optional kinetic based motion correction and incorporation of ASL derived relative perfusion measures into the pharmacokinetic modelling routine. Furthermore, the software has been shown to provide accurate estimates of pharmacokinetic parameters and can be successfully applied to PET data on a regional basis as well as voxel-by-voxel.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.