1 Introduction

Positron Emission Tomography (PET) is a powerful molecular imaging technique which can be used to quantify the density of a target of interest in vivo, using a targeted radiotracer. However, full quantification requires a long, dynamic image acquisition to cover the delivery, binding and washout of the tracer, to facilitate PK modelling and the extraction of parameters of interest. Since an acquisition time of 60 min or more is not clinically feasible, instead a single static 10 min acquisition is performed and the standardised uptake value ratio (SUVR) is calculated to estimate the relative uptake of the tracer. However, since this estimate is derived from a single time point, variations in blood flow, which influence the delivery and washout of the tracer, cannot be accounted for. Consequently SUVR measures become biased, confounding results in longitudinal studies where, in conditions such as Alzheimer’s Disease (AD), there will be both changes in blood flow and the abundance of the biological target [1].

A framework to perform PK modelling on 30 min of PET data by incorporating blood flow information from simultaneously acquired MRI was proposed in [2]. This approach halves the acquisition time by assuming that the relationship between cerebral blood flow (CBF) measured using arterial spin labelling (ASL) MRI and the PET tracer delivery parameter \(R_{1}\) can be approximated using a global linear regression derived from regional average values. However, image artefacts, limitations in image acquisition and assumptions in the ASL model used to calculate CBF mean that the errors within the CBF estimates vary across the brain. This will alter the relationship between ASL-CBF and PET tracer delivery [3]. Furthermore, due to low SNR and the resulting noise in the CBF estimates, performing linear regression for each region is non-trivial.

In this paper we propose to synthesise PET tracer delivery \(R_{1}\) maps from ASL-CBF maps and structural T1 data using information propagation from a database of subjects with ASL-CBF, T1 and PET-\(R_{1}\) data. Local similarity between the unseen candidate data and the database is used to propagate database PET-\(R_{1}\) information into the candidate subject space. The use of a local similarity metric ensures that small-scale variation can be accounted for, whilst the use of a multi-modal database ensures that the method is robust to artefacts and uncertainties in the CBF estimation. This allows the technique to be extended to voxel-wise analysis which is required to detect local changes which differentiate between forms of disease and allow early diagnosis.

The methodology was evaluated in 32 subjects participating in a study of ageing/preclinical AD where the tracer target was amyloid-\(\beta \), as quantified by the binding potential (\(BP_{ND}\)). Our approach gives a significantly improved estimation of \(R_{1}\) on both regional and voxel-wise scales compared to [2]. This translates into significantly lower errors in \(BP_{ND}\) estimates, such that amyloid burden can be accurately quantified using a single 30 min PET/MRI acquisition.

2 Methods

2.1 MRI-PET Database Construction and \(R_{1}\) Synthesis

For each of the N subjects forming the database, the MRI data (including structural T1 and ASL-CBF) are affinely registered into PET-\(R_{1}\) space. The input database thus consists of ASL-CBF and T1 data, and the output database contains the corresponding subject PET-\(R_{1}\) maps, as shown in Fig. 1. For \(R_{1}\)-map synthesis, the multi-contrast approach from [4] is used. The ASL-T1 pairs in the input database are registered to the candidate data using an affine registration to initialise the non-rigid registration [5]. The convolution based local normalised cross correlation (LNCC) is then calculated between each registered ASL-T1 pair in the database and the candidate ASL-T1 pair. This local similarity metric is then summed across the 2 channels, and is used to weight the contribution of each propagated \(R_{1}\)-map from the output database to the synthetic \(R_{1}\)-map [6].

Fig. 1.
figure 1

(a) Multi-channel registration transforms input database images into the candidate subject space. The local similarity between the images is calculated, (b) and is used to weight the propagation of \(R_{1}\) information from the output database into candidate space. (c) The synthesised \(R_{1}\) map is combined with dynamic PET data and PK modelling is used to derive the target density map.

2.2 PET Quantification Using the SRTM

The simplified reference tissue model (SRTM) [7] is used to quantify the PET data, using the linearised formulation as in (1), where a set of basis functions for \(C_{R}(t) \otimes e^{-\theta t}\) are pre-calculated over a physiologically plausible range of \(\theta \) [8]. The tracer-target interaction is modelled using one tissue compartment where the tracer concentration in the reference region \(C_R(t)\) is the input function, and the target tissue concentration is \(C_T(t)\). The model contains 3 parameters, \(R_{1}\) (the rate of delivery in the target tissue relative to reference tissue), \(k_2\) (the transfer rate constant from target tissue to blood), and the parameter of interest \(BP_{ND}\) (the binding potential which is related to target density and consequently amyloid-\(\beta \) burden). The parameters are estimated using curve fitting to \(C_R(t)\) and \(C_T(t)\) measured from \(t=0\) at tracer injection over a sufficient duration.

$$\begin{aligned} \begin{aligned} C_{T}(t) = R_{1}C_{R}(t)+\phi C_{R}(t) \otimes e^{-\theta t} \\ \text {where} \quad \phi =k_{2}-R_{1}{k_{2}}/{(1+BP_{ND})},\quad \theta = {k_{2}}/{(1+BP_{ND})} \end{aligned} \end{aligned}$$
(1)

Cerebellar grey matter is used as the reference region for \(C_R(t)\) as it is assumed to be devoid of amyloid-\(\beta \) [9]. \(BP_{ND}\), \(R_{1}\) and \(k_2\) were calculated from dynamic PET data acquired from \(t=0\) to \(t=60\) (0:60) min as the gold standard.

2.3 CBF Estimation from ASL MRI

Cerebral blood flow (CBF) maps were estimated from pseudo-continuous ASL (PCASL) data and saturation recovery images at 3 recovery times (1, 2, 4 s) [10]. The proton density, \(S_0\), was estimated by fitting the saturation recovery images for \([\text {T1}, S_{0}]\). The parameter values used were 0.9 ml/g for the plasma/tissue partition coefficient, 1650 ms for the blood T1, and 0.85 for labelling efficiency.

2.4 SRTM with Incomplete PET Scan and CBF

Population-Based Extrapolation of Reference Input C\(_{\varvec{R}}\)(t). To calculate the convolution term in (1), \(C_R(t)\) for \(t\in [0, t_e]\), where \(t=0\) at injection and \(t={t_e}\) at the end of the acquisition, is required. When the PET acquisition time is reduced to \(t\in [t_s, t_e]\), \(C_R(t)\) for \(t\in [0, t_s]\) must be estimated. This is performed as in [2], where a set of subjects with \(C_R(t)\) for \(t\in [0, t_e]\) are used to generate a mean population \(C_R(t)\), which is then scaled to the available \(C_R(t)\) \(t\in [t_s, t_e]\) to estimate a subject specific \(C_R(t)\) for \(t\in [0, t_e]\).

SRTM with CBF-Derived \({\varvec{R}}_\mathbf{1}\) and Extrapolated C\(_{\varvec{R}}\)(t). (1) is re-written \(C_{T}(t) - R_{1}C_{R}(t)= \phi C_{R}(t) \otimes e^{-\theta t}\), to group the measured terms (\(C_{T}(t)\), \(C_{R}(t)\) for \(t\in [t_s, t_e]\)) with the CBF-derived \(R_1\). Basis functions were generated using a range of \(\theta \) values and the extrapolated \(C_{R}\) for \(t\in [0, t_e]\), and \(\phi \) is estimated using a least squares fit of the data. The combination of \(\phi \) and \(\theta \) which minimise the least squares error in the fit are then used to calculate \(BP_{ND}\) and \(k_{2}\).

3 Experiments and Results

Data. Imaging data were collected from 32 cognitively normal subjects participating in Insight 46, a neuroimaging sub-study of the MRC National Survey of Health and Development, who underwent simultaneous PET and multi-modal MRI on a Siemens Biograph mMR 3T PET/MRI scanner [11]. List mode PET data were acquired for 60 min following intravenous injection of [18F]florbetapir, an amyloid-\(\beta \) targeting radiotracer. Structural 3D T1- and T2-weighted MR images were used to estimate the attenuation map [4]. Dynamic PET data were binned into \(15\,\mathrm{s}\,\times \,4\), \(30\,\mathrm{s}\,\times \,8\), \(60\,\mathrm{s}\,\times \,9\), \(180\,\mathrm{s}\times \,2\), \(300\,\mathrm{s}\,\times \,8\) time frames, and reconstructed into \(2\,\times \,2\,\times \,2\) mm voxels using in-house software with corrections for dead-time, attenuation, scatter, randoms and normalisation [12]. PCASL ASL data were acquired using a 3D GRASE readout at \(3.75\,\times \,3.75\,\times \,4\) mm and reconstructed to \(1.88\,\times \,1.88\,\times \,4\) mm resolution voxels. 10 control-label pairs were acquired with a pulse duration and post labelling delay of 1800 ms. For regional analysis and reference region delineation, T1 data were parcellated [6] and propagated into PET space.

Validation. Leave one out analysis was used such that \(N-1\) subjects were used in the database, for linear regression between ASL-CBF and PET-\(R_{1}\), and for extrapolation of the reference region \(C_{R}\). The comparison of methods was performed for both \(R_{1}\) and \(BP_{ND}\) estimation using the mean absolute error (\(\text {MAE}={1/v}\sum _{v}|I_{v}^{est}-I_{v}^{GS}|\) where I is intensity, v is the number of voxels or regions, GS is the gold standard and est is the estimated value). Statistical tests on MAE were performed using the paired Wilcoxon signed rank test. For SUVR estimation, 50:60 min of PET data were normalised by the mean reference region intensity, as this time window gave the lowest error compared to the gold standard. For PK modelling using reduced acquisition time, 30:60 min of PET data were used. PK modelling using 30:60 min of data and the gold standard ‘true’ \(R_{1}\) was calculated as the ideal case where \(R_{1}\) is estimated exactly.

3.1 Regional Analysis of \(R_{1}\) and \(BP_{ND}\) Estimation

Regional analysis was performed by averaging across 17 regions: 6 cortical grey matter regions, accumbens, amygdala, brainstem, caudate, cerebral white matter, hippocampus, pallidum, putamen, thalamus, cerebellar white matter and cerebellar grey (reference region), with right and left hemispheres combined.

Fig. 2.
figure 2

Regional errors in \(R_{1}\) and \(BP_{ND}\), outliers excluded for visualisation.

\({\varvec{R}}_\mathbf{1}\) Estimation Accuracy. Figure 2 left shows the errors in the \(R_{1}\) estimation compared to the gold standard, demonstrating that the proposed synthesis method gives improved estimates over those obtained using linear regression on the ASL values directly. For the proposed synthesised \(R_{1}\) the MAE is significantly reduced (\(p<0.001\)), from 0.08668 using the ASL regression method to 0.06326 using the synthesised \(R_{1}\). This reflects the flexibility of the synthesis method to account for regional differences in the relationship between CBF and \(R_{1}\).

Binding Potential \(\varvec{(BP_{ND})}\) Estimation Accuracy. The errors in the estimation of \(BP_{ND}\) using fixed \(R_{1}\) values and 30 min of PET data are shown in Fig. 2 right. This illustrates the known positive bias in the SUVR compared to the gold standard and the reduction in errors when PK modelling is used. For the fit using the proposed synthesised \(R_1\) map, the error is reduced compared to linear regression on the ASL values. This is not reflected in the MAE where the proposed method and the ASL regression method achieve an MAE of 0.07495 and 0.07501 respectively. However, the distributions are statistically significantly different (\(p=0.02\)). Both \(R_{1}\) estimation techniques perform significantly better than SUVR (\(p<0.001\)), and produce similar errors to those obtained using the gold standard \(R_{1}\) (\(\text {MAE}=0.06145\)).

3.2 Voxel-Wise Analysis of \(R_{1}\) and \(BP_{ND}\) Estimation

Voxel-wise analysis was performed on all voxels contained within the subject brain mask. MAE was calculated for each subject then averaged across subjects.

\({\varvec{R}}_\mathbf{1}\) Estimation Accuracy. Figure 3 left shows the reduction in \(R_{1}\) estimate error using \(R_{1}\) synthesis compared to the ASL regression technique. This is reflected in the significantly reduced MAE, which is 0.21175 for ASL regression as opposed to 0.15979 for \(R_{1}\) synthesis (\(p<0.001\)). Figure 3 right shows that the synthesised \(R_{1}\) maps are better able to capture local differences in tracer delivery than those generated using ASL regression. This is due to the shallow slope obtained from the linear regression which reduces the range of the ASL-CBF maps.

Fig. 3.
figure 3

Voxel-wise \(R_{1}\): Left- errors (outliers not shown), Right- example maps

Binding Potential \(\varvec{(BP_{ND})}\) Estimation Accuracy. Figure 4 left shows the mean error in \(BP_{ND}\) estimates compared to the gold standard and as for the regional analysis, SUVR shows a positive bias. The ASL regression technique also has a small positive bias which is not present for the synthesised and true \(R_{1}\) values. Again, the MAE is significantly lower for the synthesised \(R_{1}\) compared to the ASL-regression method. The mean errors using the synthesised \(R_{1}\) and gold standard \(R_{1}\) appear similarly distributed, however the MAE is significantly lower when using the gold standard, where \(\text {MAE}=0.1593\) compared to 0.17687 for synthesised \(R_{1}\) (\(p<0.001\)). Figure 4 right demonstrates that the proposed method is a good approximation of the \(BP_{ND}\) using the gold standard \(R_{1}\), however all of the PK modelling techniques using 30 min of data show corruption due to the noise in the PET data on a voxel level which is enhanced when halving the acquisition time.

Fig. 4.
figure 4

Voxel-wise \(BP_{ND}\) results: Left- Errors in \(BP_{ND}\) (outliers excluded for display) Right- Example \(BP_{ND}\) maps overlaid on T1 images for visualisation

4 Discussion and Conclusion

In this paper we present a novel methodology to synthesise PET tracer delivery, \(R_{1}\), maps from a database of ASL-CBF and T1 images to significantly improve the quantification of PET data using a clinically feasible acquisition time. The results were compared to a technique proposed in [2], where \(R_{1}\) values were derived directly from the ASL data by linear regression of a population of subjects with ASL-CBF and \(R_{1}\) maps. For regional analysis, the proposed synthesised \(R_{1}\) estimation performed significantly better than the ASL regression method, and the performance in the estimation of \(BP_{ND}\) following PK modelling on 30:60 min of PET data was also significantly better. The analysis was further extended to a voxel level where the higher accuracy and resolution of the synthesised \(R_{1}\) map led to a significant improvement in \(BP_{ND}\) estimation.

This work demonstrates that \(BP_{ND}\) quantification is possible on a voxel level for a 30 min PET/MR acquisition. However, fitting noisy voxel-wise data is susceptible to errors which are increased for reduced acquisition times. In future work we will apply the SRTM with a spatial constraint to improve the stability of voxel-wise fitting [13]. Furthermore, we will consider that the images within the database are imperfect, i.e. the ASL-CBF maps are susceptible to artefacts, and the \(R_{1}\) maps are voxel-wise fits to noise corrupted data. This information may be used to modify the similarity estimation between the candidate and the database to increase robustness in the presence of artefacts.