*2.5. Image Post-Processing*

The reconstructed dynamic PET data underwent two additional image correction processes. First, we conducted a frame-by-frame 3D PET image registration to minimize the between-frame misalignment due to the involuntary patient motion during the PET acquisition. The last frame was used as the reference for the frame-by-frame registration. Second, partial volume correction was performed for all PET datasets with the geometric transfer matrix method [28] provided by the PETPVC toolbox [29]. The anatomical maps were derived from the segmentation results using Freesurfer, which performs subcortical region segmentation over the T1-weighted scans [30]. All image processing methods were implemented in MATLAB (version 2020a, Mathworks, Inc., Natick, MA, USA) The Freesurfer-derived segmentation maps of the prior step were also applied to the PET dynamic data and used to extract the tissue time-activity curves of the regions of interest. We chose the following nine regions as target regions of evaluation: putamen, caudate, thalamus, hippocampus, frontal cortex, temporal cortex, occipital cortex, parietal cortex, and cerebellum.

#### *2.6. IDIF Extraction Procedure*

The IDIF extraction procedure involved three steps: an image segmentation process to extract arterial voxels, a factorization process to separate the blood signal from the tissue signal, and a scaling process to set the separated blood signal back to the accurate activity concentration units. In the first step, we aimed to segmen<sup>t</sup> the voxels that are most likely to be within the carotid arteries. For each subject, we first took the Freesurfer-derived segmentation maps to determine the lowest slice of the cerebellum and only segmen<sup>t</sup> between this slice and the overall lowest slice in the field of view. For each of those axial slices, we took the first minute of dynamic frames and searched on the left side for the voxel with the highest intensity over these six frames. Assuming an approximated diameter of 6 mm for the carotid arteries [31], this voxel was then dilated with a five-by-five diamondshaped structural element to form a segmented mask. The same operation was repeated for the right side of the same image slice to complete the carotid segmentation for this slice. All the selected slices underwent the same segmentation procedure to complete the segmentation for carotid arterial blood voxels. With the partial volume effect, the intensity of the segmented voxels that resemble the arterial blood activities is in fact a mixture of the arterial blood activities and the surrounding tissues.

To extract the arterial blood activities, it is assumed that all the surrounding tissue of the selected voxels can be approximated as a single tissue type and shares the same tracer uptake kinetics. Accordingly, the activity concentration of a specific voxel *i* can be expressed as:

$$\mathcal{C}\_{i,PET}(t) = \mathfrak{a}\_i \mathcal{C}\_{AIF}(t) + (1 - \mathfrak{a}\_i)\mathcal{C}\_{TISSILE}(t) \tag{1}$$

where *CPET* represents the PET-measured activity, *CAIF* represents the arterial blood activity, and *CT ISSUE* represents the surrounding tissue activity. *αi* is the voxel-dependent mixing fraction of the arterial blood for the specific voxel.

Assuming there is a total of *n* segmented voxels and m PET image frames for the dynamic study, an n-by-m matrix A can be formed by combining *Ci*,*PET* from all the m frames and n segmented voxels. With the goal to extract the underlying arterial input function *CAIF* and tissue activity *CT ISSUE*, we developed a novel method that decomposes the matrix A into a 2-by-m matrix H that contains the two components of *CAIF* and *CT ISSUE* and an n-by-2 weighting matrix W so that the difference between A and W\*H can be minimized. Unlike other blind matrix factorization methods, our matrix factorization method is based on physiological models and shares similar concepts with guided matrix factorization [32] and knowledge-driven matrix factorization [33] by incorporating the prior knowledge of underlying factors during the matrix factorization process. During this model-based matrix factorization (MBMF), it is assumed that *CAIF* and *CT ISSUE* can both be modeled and parameterized. Accordingly, the factorization process becomes an optimization problem that estimates the underlying parameters for *CAIF* and *CT ISSUE*, instead of a blind and direct search for the time activity curves of *CAIF* and *CT ISSUE*. In this study, we used the 7-parameter input function model developed by Feng et al. [22,34] as:

$$\mathcal{C}\_{AIF}(t) = (A\_1(t-\tau) - A\_2 - A\_3)e^{-\lambda\_1(t-\tau)} + A\_2e^{-\lambda\_2(t-\tau)} + A\_3e^{-\lambda\_3(t-\tau)}\tag{2}$$

And the tissue time-activity function *CT ISSUE* was modeled as a two-tissue compartment model output as [35]:

$$\mathbb{C}\_{TISSIZE}(t) = \mathbb{C}\_{AIF}(t) \otimes \frac{K\_1}{(B\_2 - B\_1)} \left[ (k\_3 + k\_4 - B\_1)e^{-B\_1t} + (B\_2 - k\_3 - k\_4)e^{-B\_2t} \right] \tag{3}$$

where

$$B\_{1,2} = \frac{1}{2} [(k\_2 + k\_3 + k\_4) \mp \sqrt{\left(k\_2 + k\_3 + k\_4\right)^2 - 4k\_2k\_4}]\tag{4}$$

We used the trust-region-reflective optimization algorithm within MATLAB's 'fmincon' function to perform the numerical optimization for the *CAIF*(*t*) and *CT ISSUE*(*t*). In each iteration of the parameter optimization, the parameter set of ∅ = {*<sup>τ</sup>*, *A*1, *A*2, *A*3, *λ*1, *λ*2, *λ*3, *K*1, *k*2, *k*3, *k*4 } at the current iteration was applied to Equations (2) and (3) and then used to form matrix H. The weighting matrix W at this iteration was derived from MATLAB's linear system solver that minimizes the least-square errors. The optimization then searched for the optimal solution for the parameter set ∅ by minimizing:

$$\mathcal{D} = \operatorname{argmin}\_{\{\mathcal{Q}\}} ||WH - A||^2 \tag{5}$$

After the convergence of MBMF, optimal ∅ was used to determine the extracted and normalized functions *CNorm*,*AIF* and *CNorm*,*<sup>T</sup> ISSUE* that represented the decomposed arterial blood and tissue components, respectively. Those functions were denoted with 'Norm' because they were in a normalized form after the MBMF extraction, where ∑*m j*=1 *CNorm*,*AIF*,*<sup>j</sup>* equaled to one (*j* is the frame index). The last step was to scale the extracted *CNorm*,*AIF* from an arbitrary unit to the correct physical magnitude of activity units as *CAIF*. To scale *CNorm*,*AIF* to the correct units of activity concentration, we related *CNorm*,*AIF* to the true AIF *CAIF* by a scaling factor *sAIF* as:

$$\mathbf{C}\_{AIF}(t) = \mathbf{s}\_{AIF}\mathbf{C}\_{Norm,AIF}(t) \tag{6}$$

Similarly, *CNorm*,*<sup>T</sup> ISSUE* was related to the true *CT ISSUE* by a scaling factor *sT ISSUE*. The individual activity of a specific voxel *i* was estimated by the MBMF extraction as:

$$
\dot{\mathbb{C}}\_{i,PET}(t) = w\_{i, AIF} \mathbb{C}\_{Norm, AIF}(t) + w\_{i, TISSIZE} \mathbb{C}\_{Norm, TISSLE}(t) \tag{7}
$$

*wi*,*AIF* and *wi*,*<sup>T</sup> ISSUE* were the MBMF-estimated mixing fractions for the *CNorm*,*AIF* and *CNorm*,*<sup>T</sup> ISSUE*, respectively. Note *<sup>C</sup><sup>i</sup>*,*PET* is not identical to *Ci*,*PET* since it is a weighted summation of the extracted and estimated functions *CNorm*,*AIF* and *CNorm*,*<sup>T</sup> ISSUE*.

We further derived:

$$\tilde{\mathbf{C}}\_{i,PET}(t) = \frac{\omega\_{i, AIF}}{s\_{AIF}} \cdot \mathbf{C}\_{AIF}(t) + \frac{\omega\_{i, TIISILE}}{s\_{TIISILE}} \cdot \mathbf{C}\_{TISSIIE}(t) \tag{8}$$

From Equation (1), it was assumed that the mixing fractions of the AIF and tissue activity should sum up to one under ideal situation. Therefore, the optimal values of *sAIF* and *sT ISSUE* were optimized by minimizing:

$$\left\| \left. \right\|\_{AIF\prime} \left. \right\|\_{TISSIE} = \operatorname\*{argmin}\_{\left( \left. \right|\_{AIF\prime} \left. \right|\_{TISSIE} \right)} \sum\_{i=1}^{n} \left( \frac{\omega\_{i,AIF}}{s\_{AIF}} + \frac{\omega\_{i,TITS}}{s\_{TISSIE}} - 1 \right)^{2} \tag{9}$$

where *<sup>s</sup>*<sup>ˆ</sup>*AIF* and *<sup>s</sup>*<sup>ˆ</sup>*T ISSUE* denoted the optimized scaling factors for extracted AIF and tissue activity, respectively. n denoted the total number of segmented voxels. We used the trustregion-reflective optimization algorithm with MATLAB's 'fmincon' function to perform the numerical optimization for minimizing Equation (5). The estimated *<sup>s</sup>*<sup>ˆ</sup>*AIF* was then plugged into Equation (6) to scale the IDIF into the correct physical units.

## *2.7. Metabolite Correction*

Since the IDIF method extracts the whole-blood activity and cannot separate the metabolite signal from the unmetabolized tracer, we used a population-based approach to convert the extracted IDIF to a metabolite-corrected plasma time-activity curve [36,37]. With the metabolite data and plasma time-activity curves measured in the blood-sampling cohort, the individual parent fraction was multiplied to the plasma-to-whole blood activity fraction. The individual composite fraction was averaged for each time point at the 5th, 15th, 30th, and 60th minutes and then fitted to a single exponential function with constant [37,38]. As a result, the corrected IDIF is expressed as:

$$\mathcal{C}\_{IDIF,MCPC}(t) = \\$\_{AIF}\mathcal{C}\_{N,AIF}(t) (1 - 0.29(1 - e^{-0.03t})) \tag{10}$$

where *t* is in the unit of minutes, and *MCPC* denotes metabolite-corrected plasma concentration of activity.
