Model-derived input function

PET imaging of cerebral glucose metabolism with [18F]FDG is typically based on the irreversible two-tissue compartment model [3], in which [18F]FDG enters the tissue via glucose transporters defined by an influx rate constant \(_\) (in mL/g/min). Once in the first compartment (i.e., the free pool), [18F]FDG can either return to the blood pool at an efflux rate defined by \(_\) (in min−1) or be phosphorylated at a rate defined by \(_\) (in min−1) and subsequently trapped in the metabolic pool. These processes are defined by the following two differential equations:

$$\frac_\left(t\right)=__\left(t\right)-\left(_+_\right)_\left(t\right)$$

(1)

$$\frac_\left(t\right)=__\left(t\right)$$

(2)

where \(_\left(t\right)\) represents the activity concentration in the free pool and \(_\left(t\right)\) in the metabolic pool. Eq. (3) provides the solution to Eqs. (1) and (2) for the total measured activity (i.e., \(_\left(t\right)=\left(1-_\right)\left(_\left(t\right)+_\left(t\right)\right)+__\left(t\right)\)), which includes an additional term \(_\) (in mL/g) to account for blood-borne [18F]FDG activity [26]:

$$_\left(t\right)=\left(1-_\right)\frac_}_}\left(_+_^_t}\right)*_\left(t\right)+__\left(t\right)$$

(3)

where \(*\) represents the convolution operation; \(_\left(t\right)\) and \(_\left(t\right)\) represent the plasma and whole-blood activity concentration of [18F]FDG, respectively; and \(_=_+_\).

The previous equations can be combined and rearranged to derive an expression for the input function, i.e., the MDIF (see Appendix for full derivation). By selecting the WB TAC (\(_\left(t\right)\)) to define the MDIF, considering it has the highest SNR obtainable in a dynamic PET brain image, and assuming \(_\left(t\right)=R_\left(t\right)\) (where \(R\) is the blood-to-plasma ratio), the MDIF is given by:

$$_^\left(t\right)=\frac_}\left[_\left(t\right)+\frac_}_}_\left(\frac_-_}_-_}^_t}-\frac_-_}_-_}^_t}\right)*_\left(t\right)\right]$$

(4)

where \(_}=a\mp \sqrt^-b}\), with \(a=\frac\left(\frac_}_}_+_\right)\) and \(b=\frac_}_}__\).

By employing this approach, the input function is defined by the four model parameters (\(_, _,_\) and \(_\)) that characterize the WB TAC. Eq. (4) was developed for [18F]FDG given it obeys the irreversible 2TCM and follows two key assumptions. First, [18F]FDG does not require metabolite correction as the tracer and its metabolites are trapped in the cell and clearance of [18F]FDG-6-phosphate is slow for acquisition times of up to 60 min [27]. Second, the blood-to-plasma ratio was assumed to be a constant value for simplification purposes [28], although this assumption might not be valid for non-primates [29]. The appendix provides the general solution for reversible radiotracers. We are currently developing a variation of the MDIF method for tracers that require metabolite correction, but this is beyond the scope of the current work.

MDIF SIME implementation

The SIME procedure [17] requires TACs from multiple volumes-of-interest (VOIs) to derive the whole-brain kinetic parameters that define the MDIF (Eq. (4)). These TACs can be from anatomical VOIs or grouped based on functional similarity, as previously suggested by Wong et al., to identify TACs with distinct kinetics [18]. We implemented a k-means clustering algorithm (explained below) to group TACs by functional similarity. After initial investigation, we observed that 3 clusters of each tissue type (i.e., grey [GM] and white matter [WM]) were sufficient to derive the MDIF, with little-to-no improvement observed by using more than 6 TACs (results not shown). A diagram outlining the MDIF SIME implementation is shown in Additional file 1: Fig. S1.

The parameters needed to compute the MDIF using Eq. (4) were obtained by simultaneously fitting six TACs to the irreversible 2TCM solution (Eq. (3)) using a non-linear least squares (NLLS) fitting routine. The algorithm runs until the cost function, comprising of the residual sum of squares (RSS; in Eq. (5)), is minimized and the WB model parameters used to define the MDIF is obtained.

$$RSS=\sum_^\sum_^_^\left(_\right)-_\left(_\right)\right)}^+\sum_^_^-_\right)}^$$

(5)

where \(_\left(t\right)\) is the [18F]FDG-PET activity concentration of the jth VOI, \(_^\left(t\right)\) is the corresponding model estimate, \(n\) is the number of VOIs (not including the WB TAC), and \(T\) is the number of time-frames. Given the need for scalers (or anchors) to accurately obtain an input function with SIME [15], our procedure included two late time-point blood samples to act as anchors. \(_\) is the kth blood sample activity concentration, and \(_^\) is given by Eq. (4) at the same time as the blood sample.

The fitting procedure was performed in MATLAB (The MathWorks Inc., R2023a) using the optimization routine lsqnonlin and four parameters per TAC were included in the fitting routine (\(_\), \(_\), \(_\), and \(_\), for a total of 28 parameters). The upper bounds were set to 0.2 mL/g/min, 0.4 min−1, 0.2 min−1, and 0.10 mL/g for \(_\), \(_\), \(_\), and \(_\), respectively (based on literature values, see Additional file 1: Table S1). All lower bounds were set to 0.01. A constraint was added to the nonlinear optimization function implemented in this study to ensure the distribution volume (= \(_/\left(_+_\right)\)) is less than unity (i.e., \(_<_+_\)), which was based on literature values (Additional file 1: Table S1).

k-means clustering

For the animal and human studies, VOIs were generated by clustering the TACs based on functional similarity by implementing a k-means clustering algorithm (Euclidean distance, 500 iterations, 10 replicates) that used the k-medoids MATLAB function kmediods [30]. For clustering only, PET images were denoised by applying the 3D highly constrained backprojection (HYPR3D) method [31], in which a Gaussian filter kernel of standard deviation 3 voxels was used (equivalent to 6.3 × 6.3 × 6.1 mm3 for the human PET data used in this study). Although HYPR3D introduces bias to the dynamic signal, which is known to affect quantification [32], this bias is not expected to influence the k-means clustering algorithm. Instead, denoising with HYPR3D improves the clustering of similar TACs. Additional k-means clustering implementation details are given in the appropriate sections below.

IDIF SIME implementation

For comparison to our MDIF SIME method, the original SIME approach based on the work by Feng et al. and Wong et al. [17, 18] was implemented to extract an IDIF by estimating its 7 parameters (Eq. (6)), while simultaneously fitting six VOI TACs (WB TAC not included) using a NLLS fitting routine. The same cost function from the MDIF SIME was used (Eq. (5)), except for an additional weight (= 10) included in the blood samples portion of Eq. (5). All fitting was performed in MATLAB using the optimization routine lsqnonlin and four parameters per TAC were included in the fitting routine (\(_\), \(_\), \(_\), and \(_\); same constraint and bounds as described above) in addition to the 7 parameters required for the IDIF (Eq. (6)): \(_\), \(_\), with i = 1, 2, 3, and a delay term \(\delta\) (total of 31 parameters). Upper bounds were 4000 kBq/mL/min, 100 kBq/mL, 50 kBq/mL, 25 min−1, 1 min−1, and 0.1 min−1, for \(_\), \(_\), \(_\), \(_\), \(_\), and \(_\), respectively. Lower bounds were set to zero. The delay term (\(\delta\)) was set to vary within ± 10 s of an initial delay computed a priori based on the rise of the WB TAC.

$$_^\left(t\right)=(_\left(t-\delta \right)-_-_)^_\left(t-\delta \right)}+_^_\left(t-\delta \right)}+_^_\left(t-\delta \right)}$$

(6)

Standalone fitting routine

To extract microparameters from anatomical VOIs, each TAC was fit to the irreversible 2TCM solution (Eq. (3)) using the MATLAB optimization routine lsqnonlin. Four parameters were included in the weighted NLLS (WNLLS) fitting routine (\(_\), \(_\), \(_\), and \(_\)). Upper bounds were set to 0.5 mL/g/min, 0.5 min−1, 0.2 min−1, and 1 mL/g for \(_\), \(_\), \(_\), and \(_\), respectively. All lower bounds were set to zero. The cost function was defined as a weighted RSS (WRSS), given by \(WRSS=\sum_^__^\left(_\right)-_\left(_\right)\right)}^\). Weights for the fitting (\(_\)) were defined as the inverse of variance of the PET measurement error (i.e., the ratio between the time-frame duration and the WB activity concentration, normalized to the maximum weight) [33]. A brain density of 1.05 mL/g was used throughout the analyses.

Simulations

A theoretical AIF (\(_\left(t\right)\)) was generated by Eq. (6) using \(_\) = 850 (in arbitrary units [a.u.]/min), \(_\) = 22 and \(_\) = 21 (in a.u.), \(_\) = 4, \(_\) = 0.12 and \(_\) = 0.01 (in min−1), and \(\delta\) = 0. These constants were obtained from Feng et al. and are based on experimental AIF measurements from a human study [34]. The theoretical AIF was used in the irreversible 2TCM solution (Eq. (3)) to generate simulated regional TACs for six theoretical VOIs (Table 1). The microparameters shown in Table 1 were chosen from a variety of studies [9, 35,36,37] and reflect the need for TACs with distinct kinetics for the SIME approach [18, 38]. The net clearance rate constant of [18F]FDG (\(_\) = \(__/_\), in mL/g/min), is also shown in Table 1. A theoretical WB TAC was obtained by using the average microparameters from the six VOIs. Fig. 1 shows the theoretical AIF alongside simulated TACs. These curves and microparameters were considered the ground-truth for the simulations. Lastly, simulations used anchors extracted from the theoretical AIF at 28.5 and 53.5 min timepoints, as indicated by the arrows in Fig. 1A.

Table 1 Microparameters used to define regional TACs for six theoretical VOIs included in the simulationsFig. 1

Theoretical A AIF and B TACs used in the simulations. Curves were interpolated to the time frames used in the human experiments (symbols). The arrows represent the scalers used in the SIME simulations. AIF: Arterial input function; TAC: Time-activity curve; VOI: Volume-of-interest; WB: Whole brain

To confirm the accuracy of the MDIF SIME routine, errors due to estimating the WB microparameters (Eq. (4)) with the SIME routine were evaluated. The MDIF was generated as described above and used to fit the six VOI TACs (Table 1) to the irreversible 2TCM solution (Eq. (3)) in the standalone fitting routine. Then, the best-fit estimates of the four parameters (\(_\), \(_\), \(_\), and \(_\)) were compared to their input values.

Animal study

The MDIFs were assessed in comparison to AIFs using retrospective data from animal experiments [12] collected at the Lawson Health Research Institute. All experiments were conducted according to the regulations of the Canadian Council on Animal Care and approved by the Animal Care Committee at Western University. Data from juvenile Duroc pigs were collected on a 3 T hybrid PET/MR scanner (Biograph mMR, Siemens Heathineers, Erlangen, Bavaria, Germany) using a 12-channel PET-compatible receiver head coil. The five animals included in this study (weight, 21 ± 2 kg; injected [18F]FDG activity, 90 ± 20 MBq; blood glucose level, 5.0 ± 1.6 mmol/L) were scanned under an anaesthetic combination of 1–3% isoflurane and 6–25 mL/kg/h intravenous infusion of propofol. [18F]FDG was administered via a cephalic vein immediately followed by a 60-min acquisition of PET data in list-mode. Arterial blood was continuously withdrawn from a femoral artery using an automated MR-compatible system (Swisstrace GmbH, Menzingen, Canton of Zug, Switzerland) at a sampling rate of 4 mL/min for the first 5 min, 1 mL/min for the next 5 min, and 0.5 mL/min for the remaining 50 min. Arterial blood samples were processed and converted into AIFs via rebinning to match the reconstructed PET time frames using pSample (PMOD Technologies LLC). Each measured AIF was corrected for dispersion and delay [39].

Dynamic PET images were reconstructed offline with the Siemens e7 tools into 51 time-frames (15 × 2 s, 6 × 5 s, 8 × 15 s, 4 × 30 s, 5 × 60 s, 5 × 120 s, 8 × 300 s) using a 3-dimensional ordinary Poisson ordered subset expectation maximization (3D-OP-OSEM) method (3 iterations, 21 subsets) with corrections for decay, random coincidences, dead-time, detector normalization, data rebinning, attenuation, and scatter. Images were reconstructed into a 344 × 344 × 127 matrix with voxel size of 0.8 × 0.8 × 2.0 mm and smoothed with a 3D gaussian filter of 4 mm. Simultaneously to the PET acquisition, T1-weighted MR images were acquired (magnetization-prepared rapid gradient-echo sequence [MPRAGE]; repetition/echo times [TR/TE], 2000/2.98 ms; inversion time, 900 ms; field-of-view [FoV], 256 × 256 mm2; isotropic voxel size, 1 mm3; flip angle [α], 9°; 176 slices). At the end of the experiment, the animals were euthanized according to the animal care guidelines and transported to a computed tomography (CT) scanner to obtain a post-mortem CT-based attenuation correction map.

For the k-means clustering algorithm, a semi-automatic procedure was used to define a VOI encompassing the brain in each anatomical slice. The final WB VOI, which was the composite of all the slices, was used to generate seven clusters with the k-means clustering approach as described above. The cluster containing mostly non-brain regions (i.e., blood vessels and cerebrospinal fluid voxels) was excluded from the SIME algorithm. The SIME approach incorporated arterial samples collected at 20–30 min and 40–60 min post-injection as anchors. Although a constant blood-to-plasma ratio may not be valid for the porcine model used in this study [29], we assumed \(R\) = 1 throughout the analysis of the animal data.

Human study

Retrospective data from 18 neurologically healthy volunteers (44 ± 15 years, 77 ± 17 kg, 9 M/9F; average injected activity, 180 ± 40 MBq [range 130–260 MBq]; average injected activity per body weight, 2.43 ± 0.39 MBq/kg [range 1.80–3.70 MBq/kg]; blood glucose level, 5.0 ± 0.7 mmol/L [range 4.1–6.8 mmol/L]) collected at the Lawson Health Research Institute were used to assess the feasibility of obtaining MDIFs from dynamic [18F]FDG PET data. The study was approved by the Western University Health Sciences Research Ethics Board and was conducted in accordance with the Declaration of Helsinki ethical standards. Participants provided written informed consent in compliance with the Tri-Council Policy Statement of Ethical Conduct for Research Involving Humans.

Scanning was performed on a 3 T hybrid PET/MR scanner (Biograph mMR) using a 16-channel PET-compatible coil (12- and 4-channel head and neck coils, respectively). Each participant had their head immobilized during the 60-min list-mode PET acquisition performed following the [18F]FDG bolus injection, which was immediately followed by saline flush. PET data were reconstructed offline with the Siemens e7 tools into 50 time-frames (15 × 2 s, 6 × 5 s, 8 × 15 s, 4 × 30 s, 5 × 60 s, 3 × 120 s, 8 × 300 s, and 1 × 240 s) using the iterative reconstruction algorithm 3D-OP-OSEM (3 iterations, 21 subsets, 3D gaussian filter of 2 mm3, and a zoom factor of 2). Dynamic PET images were reconstructed into a 172 × 172 × 127 matrix with voxel size of 2.1 × 2.1 × 2.0 mm (FoV, 359 × 359 × 258 mm3). Corrections for decay, random coincidences, dead-time, detector normalization, data rebinning, attenuation, and scatter were performed. MR-based attenuation correction was performed with a vendor-provided ultrashort echo time MRI sequence (TR, 11.94 ms; TE, 0.07 and 2.46 ms; α, 10º; FoV, 300 × 300 × 300 mm3; voxel size, 1.6 × 1.6 × 1.6 mm3). Motion correction was applied to the time frames 24–50 (corresponding to 1.5 to 60 min) by realigning each frame to the mean image using statistical parametric mapping (v12, SPM12; https://www.fil.ion.ucl.ac.uk/spm/software/spm12/). Each session also included acquiring an MPRAGE image (TR/TE, 2400/2.25 ms; α, 8°; FoV, 205 × 205 mm2; 240 slices; voxel size, 0.8 × 0.8 × 0.8 mm3), which was coregistered to the PET space using SPM12.

The MPRAGE image was segmented into six tissue classes using SPM12. A WB mask was created by combining GM and WM tissue probability maps (80% threshold; cerebrospinal fluid voxels excluded). Five clusters for each tissue class (i.e., GM and WM) were generated by implementing the k-means clustering algorithm as described above. Out of the 5 clusters, the 3 largest ones were chosen for the SIME algorithms. Two venous samples collected during dynamic PET imaging (one 20–30 min post-injection and the other 45–60 min post-injection) were used as scalers when arterial-to-venous [18F]FDG equilibrium was assumed [6, 40]. The blood-to-plasma ratio was assumed to be ~ 0.9 (\(R\) = 0.9) to derive the MDIF from the human data [28].

Variance of estimated \(}}_}}\)

Coefficient of variation (CV) of \(_\) estimates obtained by fitting each cluster TAC to the 2TCM solution was computed as \(CV\left(\%\right)=100\sqrt^}/_\), where the variance (\(^\)) was obtained from the error propagation rule: \(^=\left[\begin\frac_}_}& \frac_}_}& \frac_}_}\end\right]\gamma \left(_,_,_\right)\frac_}_}& \frac_}_}& \frac_}_}\end\right]}^\), where the Hessian (covariance) matrix (\(\gamma\)) was approximated as \(\gamma =^J\right)}^WRSS/T\) (\(J\) is the Jacobian matrix from the MATLAB function lsqnonlin and \(T\) the number of time frames).

Goodness of fit

Following the derivation of the MDIF and IDIF with their respective SIME algorithms, goodness of fit of each input function was assessed via the WRSS from the standalone fitting routine. For this, each input function was used to fit the GM and WM TACs, and WRSS was calculated for the first 3 min of data, as well as for the entire 60 min of data. In addition, fitting to the irreversible 2TCM solution was evaluated qualitatively for the GM and WM TACs of one representative subject, as well as for TACs from two randomly selected GM and WM voxels.

Regional measurements

Anatomical MPRAGE images were normalized to the Montreal Neurological Institute (MNI; McGill University, Montreal, Quebec, Canada) space using SPM12 and the same deformation field was applied to the dynamic PET data. Average TACs were extracted for GM and WM (99% threshold; subject space), as well as twelve VOIs (frontal, occipital, parietal, and temporal lobes, as well as insula, cingulate, hippocampus, precuneus, caudate nucleus, putamen, thalamus, and cerebellum; MNI space) using the automated anatomical labelling atlas (Wake Forest University PickAtlas, http://fmri.wfubmc.edu/cms/software). The mask generated for the twelve VOIs was multiplied by the GM mask (80% threshold; MNI space) to extract the central portion of the VOIs, avoiding PVE at the boundaries.

Patlak analysis (\(t>^\) = 20 min) was performed for each VOI, from which the slope (\(_\); in mL/100 g/min) was converted to CMRGlu (in µmol/100 g/min) by assuming a lumped constant of 0.52 [35] and incorporating measurements of the blood concentration of glucose (in µmol/mL).

Voxelwise parametric images of microparameters were generated by incorporating the MDIF into the variational Bayesian approach proposed by Castellaro et al. [33]. Specifically for this step, the k-means clustering algorithm described above was used to generate a functional atlas containing 40 clusters for each subject: 20 GM, 10 WM and 10 non-brain clusters, the latter including cerebrospinal fluid and voxels surrounding the brain. In the variational Bayesian fitting routine, each cluster TAC was fit to the irreversible 2TCM solution, with its rate constants used as priors for all intra-cluster voxels assuming a Gaussian distribution of possible values. Parametric maps of CMRGlu were computed from \(_\) images (= \(__/\left(_+_\right)\)) [35]. Resulting parametric images (i.e., CMRGlu, \(_\), \(_\), \(_\), and \(_\)) were normalized to the MNI space using SPM12.

Statistics

Area under the curve (AUC) was used to assess the similarity between input functions. Percent error was calculated as \(100\left(m-\widehat\right)/\widehat\), where \(m\) and \(\widehat\) are the observed and expected values, respectively. Percent difference between two observed measurements \(_\) and \(_\) was computed as \(100\left|_-_\right|/\overline\), where \(\overline\) is the average measurement. A repeated measures two-way ANOVA was used to evaluate if estimates were affected by input function (MDIF vs. IDIF) and VOI; multiple comparisons tests were used to evaluate differences between MDIF and IDIF estimates. Linear regression was used to compare MDIF and IDIF results, from which the line-of-best fit was obtained alongside the 95% confidence intervals (CIs) for both the slope and intercept. Correlation was assessed by means of the Pearson correlation coefficient (\(\rho\)). All datasets were found to be normally distributed with normal Q-Q plots. Statistical tests were performed using GraphPad Prism (version 9, GraphPad Software, San Diego, California USA, www.graphpad.com). Statistical significance was defined by \(\alpha\) = 0.05. Measurements are expressed in terms of mean ± one standard deviation alongside the CIs in square brackets when relevant.