Raw data from 12 patients with CKD-MBD enrolled in our previous study of method implementation (2) were reanalyzed for reproducibility. For the optimized analyses in this study, one obese patient was excluded due to extremely poor counting statistics resulting in outlying data in all series with small VOI volumes ≤1 ml.

In the original study, 5-ml venous blood samples were collected at −5, 30, 40, 50, 60, and 90 min after injection and prepared for well counting. The well-counter and PET/CT scanner were cross-calibrated as previously described by Vrist et al. (2). Whole blood and plasma data from our original study were reused for this study.

To convert measured activity from image-derived whole blood to plasma activity curves, plasma to whole blood activity ratios (PWR) were calculated for each of the samples.

In addition, all plasma values were transformed using the natural logarithm function. The slope and intercept of the resulting line at 40–60 and 40–90 mpi (plasma exponentials) were then determined by linear regression analysis and extrapolated back to the time for the peak. The plasma exponentials were used as a substitution for the corresponding image exponentials for the construction of the various input functions as described below. The 90-min samples were included for two reasons: 1.

Interpolation of plasma data for calculation of K_i at the time of whole-body (WB) data acquisition, instead of extrapolation.

The original PET/CT images were acquired on a Siemens Biograph mCT-4R 64 slice PET/CT scanner with a 22-cm axial field of view (FOV). The participants were positioned with the heart and the thoracic vertebrae Th7–Th10 centered in the FOV. After an intravenous bolus injection of 150 MBq ¹⁸F-NaF flushed with 20 ml isotonic saline, a 60-min list-mode dynamic scan was acquired immediately followed by a WB scan from the middle of the femur to the vertex of the skull acquired in 6–7 FOVs of 3 min per bed position.

The original PET images for dynamic analysis were re-binned into 50-time frames: 20 × 3 s, 12 × 5 s, 4 × 30 s, and 14 × 240 s. The reconstruction of PET scans used filtered back-projection with a Gaussian filter of 5 mm and a matrix size of 256 × 256 (3.2 mm × 3.2 mm × 1.4 mm). All dynamic images were automatically decay-corrected to the study injection time (study reference time). Image data from the WB scan were automatically decay-corrected to the start of the WB scan requiring additional decay correction to the study injection time for comparison with dynamic data.

Low-dose CT scans were performed, and the images reconstructed in three utilization-dependent series: (1) attenuation correction; (2) localization and identification of the thoracic vertebrae in the dynamic scan; and (3) localization of the relevant bone regions in the WB scan.

PMOD® version 4.206 software (PMOD Technologies LLC, Switzerland) was used for the non-linear regression analysis of the dynamic data and analysis of the static WB data.

For the reanalysis of the original data, all VOIs were constructed as described in the original study (2), whereas for the new optimization studies the vertebral VOIs (Th7/8–Th10/11) were drawn using a circular region of interest (ROI) in a single slice of 1-voxel thickness (1.4 mm) with a diameter of 4 voxels (12.8 mm), which was centered in the spongious bone and propagated through six slices resulting in a cylinder with a volume of approximately 1 ml in each of the four vertebrae (Figure 1A). Care was taken to avoid cortical bone and areas of obvious abnormal bone turnover (e.g., in a compressed vertebrae or for near lying large osteophytes). The data from the four individual VOIs were then unified to one combined VOI using the Union function in the PMOD VOI tools to improve counting statistics.

These combined VOIs were kept constant for each subject dataset and used throughout all studies of the various input functions described below.

The VOIs for the IDIFs were constructed in five different ways. In the two studies using hot contouring of activity in the LV, the LV was delimited using a box VOI (approximately 40 mm × 40 mm × 40 mm). Within this VOI, a hot contour was drawn using PMOD's contour tool with the cutoff values described below.

The VOI definitions for these basic (unmodified) input functions were: 1.

For reanalysis with original VOI definition (LV-Orig): a hot contour delineating 50%–70% of the max value in the box (Figure 1D). VOI volume: 50.8 ± 25.2.

All contours were visually inspected for overlap with myocardial background and having a distance of at least two voxels (∼6.4 mm) from the myocardial background in all planes.

The basic input functions were derived from the activity data in the various input VOIs and converted to plasma activity by multiplying the whole blood activity with the PWR in all frames.

All basic input functions – LV-Orig, LV-New, AO-Fix, and AO-Peak with the exception of the LV-Fix – were corrected for PVE and Bg using the recovery coefficient RC_ß as described previously (2, 3):(2)RIDIF(t)=RCß⋅CA(t)+(1−RCß)⋅CBg(3)suchthatRCß=(RIDIF(t)−CBg(t))/(CA(t)–CBg(t)),where R_IDIF is activity measured in VOIs in either LV or AO, C_Bg is background activity, and C_A is the “true” activity in arterial blood, which after 30 min equals the activity in venous blood (3).

For the comparison and estimation of the influence of background activity, all input data were also corrected using the simpler calibration factor RC_CF (1) without correction for spill-over from background activity:(4)RCCF=RLV(t)/(CLV(t).This equation can be used with VOIs with small volumes ≤1 ml, e.g., LV-Fix placed so far from surrounding background activity (>25 mm) that spill-over from background activity mostly can be ignored. An illustrative phantom measurement has shown that activity spill-over beyond this distance is rather constant and less than 8.2% (see Supplementary Section 1.5).

After correction, basic input functions were given identifiers by adding either -ß or -CF to the geometry nomenclature—e.g., “LV-VOI-Orig-ß” and so on.

In all the basic image curves the terminal exponentials from 40 to 60 mpi were replaced by an exponential calculated from the plasma samples using either 40–60 mpi plasma exponentials or 40–90 mpi plasma exponentials.

All image values of the IDIFs were transformed using the natural logarithm function. The slope and intercept for the resulting line from 40 to 60 mpi (“terminal exponential”) was then determined by linear regression analysis and extrapolated back to the time for the peak.

The values of this line were subtracted from the values of the entire logarithmic curve to obtain the residual curve of the initial “fast image exponentials.”

The input curves were then reconstituted in two ways: 1.

The logarithmically transformed plasma curve derived from the plasma samples (40–60 mpi and 40–90 mpi) were added to the logarithmically transformed residual curve and then exponentially retransformed as in our original study (2)—“the multiplicative logarithmic” (“Log”) method with the given identifiers “LV-Orig-Pl-40-90-Log (t₆₀), and so on.

K_i ml/min⁻¹·ml⁻¹ was calculated as the mean value of four thoracic vertebrae (Th7–Th10 or Th8–Th11). No correction for delay was made as we found the delay to be of only a few seconds, and attempts at correcting the very noisy data in the first few acquisition frames failed to make the data more consistent.

The PMOD® version 4.206 software (PMOD Technologies LLC, Switzerland) was used to perform a two-tissue compartment dynamic NLR analysis of ¹⁸F-NaF-turnover as described by Hawkins et al. (8). The exchange of ¹⁸F-NaF between the compartments—plasma, extravascular, and bone—is described by the kinetic parameters K₁-k₄ and the parameter for regional bone plasma clearance K_i is defined as(5)Ki=K1k3k2+k3This NLR method was further analyzed for two values of the fractional blood volume vB, which was either fixed at 0.05 (vB-Fix) or determined by PMOD as a free-fitted (vB-Free) parameter.

Assuming the efflux rate constant k₄ to be negligibly small (k₄ ≈ 0 min⁻¹), the Patlak graphical analysis (9) provides a simpler alternative analysis method for measuring K_i as described by equation 6 (5):(6)CB(T)CPl(T)=Ki∫0T⁡CPl(t)dtCPl(T)+V0This equation approximates a straight-line fit with K_i as the slope where C_B and C_Pl are the respective concentrations of tracer bound in bone and freely diffusible in plasma at each time point (t). V₀ is the intercept of the ordinate and represents the apparent volume of distribution.

K_i was calculated from the 60-min dynamic PET/CT scan using a bone TAC and each of the various IDIF modifications or selected semi-population input functions as described in Section 2.14. All Patlak results generated by the various input functions were compared to the basic corrected IDIF within each series and to the corresponding PMOD results.

For the future analysis of static WB scans with the “best” input function as defined below, we derived an optimized three-exponential SPIF, as previously described by Blake et al. (6). The SPIF was derived from a population residual curve (PopRes) and then added to the terminal exponential derived from the plasma samples (described above under “Blood samples”), where the PopRes was derived from the corrected basic IDIFs scaled to a reference dose of 100 MBq, as used in the study by Puri et al. (1):(7)SPIF=PopRes⊗Inj⋅dose100⊕(PlasmaExp)The residual curve represents the sum of the early fast exponentials and was derived by subtracting the terminal exponential (all data values ≥40 mpi) from the entire image-derived curve.

All residual curves were adjusted so that the times of peak count rate for all curves were coincident with the most frequent unadjusted peak time (16.5 s). The residual curves were then averaged to define the PopRes for combination with the plasma exponential to make the SPIF.

A mathematical model for the PopRes was fitted and is described in Supplementary Section 1.4.

As in the previous study (2), the static scan analysis was performed using a modified Patlak analysis with only two data points (4, 6). The first (0,V₀) was obtained from either the original dynamic function or the reconstructed SPIF using the individual patient's blood samples. The second time point was obtained as the start time for the WB scan and the activity in the vertebrae (C_B(t), C_Pl(t)) at that time using the same VOIs as for the dynamic scan but adjusted for proper alignment. The values of the static scan data were corrected for decay to the time of injection/start of the dynamic scan for comparison with the dynamic data.

The K_i values were then calculated as the slope of the line between (0,V₀) and (C_B(t), C_Pl(t)) (4–6).

The results are presented as mean ± standard deviation (SD) but as the majority of the datasets obtained using the different input functions were not supposed to be normally distributed, visual representation of the data are presented in box-and-whisker plots showing: Max, 75% quartile, Median, 25% quartile, and Min values as well as the difference (the box) between the 75% quartile and the 25% quartile (the interquartile range (IQR)).

Differences-between-method results were evaluated using Bland–Altman plots (10) showing the mean differences between the corresponding data points and the upper and lower 95% confidence limits.

Correlations between K_i values obtained using different analysis models were calculated using Pearson's correlation coefficient but were not used as selection criteria when choosing the best agreement between the methods/parameter choices.

The percentage coefficient of variation of the PopRes curves was obtained as the ratio between the SD and the mean of the PopRes curve.

The PMOD χ² test was used to evaluate the model fit of the input functions to the applied model.

The original and reanalyzed data were compared using a paired, two-tailed t-test.

Comparison of the reanalyzed and original data (LV-Orig β series) using the original method definitions are summarized in Supplementary Section 1.2. The reanalyzed and original data (PWR, RC_ß, AUCs of derived input curves, NLR-K_i and Patlak-K_i results) were all comparable.

The mean inter-observer difference for the NLR-ß-K_i values were 0.0006 ± 0.0052, r = 0.91, and p = 0.72 (NS) and −0.0018 ± 0.0042. r = 0.89, and p = 0.17 (NS) for the Patlak-ß-K_i values.

The VOI volumes for input VOIs, Bg, and Bone with their corresponding RCs for the basic input IDIFs are shown in Table 1. The LV-Orig VOI has the lowest RC_ß (0.73 ± 0.17) compared with the AO-Fix VOI's RC_ß (0.99 ± 0.09). In contrast, the AO-Peak VOIs had RCs of 1.73 ± 0.62. The RC_CFs were generally higher as spill-over was not included.

All median and range data supporting the following observational results for the uncorrected and corrected Basic IDIF AUCs can be found in the Supplementary Section 1.3 (Figure S1A,B) box plots:

Correction of the basic IDIFs with either RC_ß or RC_CF tended to shift the AUCs toward higher AUC values (Supplementary Section 1.3, Figure S1B) but with less change in AUC values for smaller VOIs, the highest values being 377.53 ± 66.88 kBq·min and the lowest 293.57 ± 66.09 kBq min. The AO-Fix curves were almost unchanged while the AO-Peak curves showed lower values (Δ-AUC-mean: −187.1 ± 122.21 kBq min) with an unacceptably wide 95% CI of −426.6 to 52.5 kBq min.

The effect of substituting the final image exponential in the corrected basic IDIFs with plasma exponentials is shown in the box plots (Supplementary Figure S2A–D).

IDIFs that were modified using the exponential additive method generally had a lower IQR and range for the 40–90 mpi data analysis compared with the corresponding 40–60 mpi analysis curves (Δ-IQR: −33.6–1.8; range: −7.7 to −19.7). However, the mean AUCs for the LV input functions were not significantly different from the corresponding basic input functions as shown in Supplementary Section 1.3 (Table S3). This also applied to the SPIFs reconstructed from the LV-New-ß series. The small differences in AUCs between the SPIFs and their corresponding IDIFs were not significant.

In Figure 2A, the model of our PopRes was compared with data from the observed, optimized PopRes curve. The visual fit was very good and use of the curve results in comparable K_i results, as described in the following subsections. The mean difference was −0.40 ± 2.79 with a 95% CI of −5.87 to 5.06. The AUC_1800sec was 4,078 kBq s.

The mathematical best-fit curve for our PopRes is presented in Supplementary Section 1.4.

As a test of inter-observer reproducibility, all original image data (2) were reanalyzed using the original method parameters (Section 3.1). The observed slight difference (NS) in the plasma-to-whole-blood ratio (PWR-reanalyzed = 1.16 ± 0.02. vs. PWR-Orig = 1.17 ± 0.03) lies within the range of values presented in Figure 3 in Vrist et al. (2) and is most likely due to the more selective exclusion of plasma outliers in the reanalysis in the present study.

The RC_ß correction factor used in the present study (0.72 ± 0.17) was comparable to that of the original study (0.69 ± 0.15).

The AUC of the basic input curves in this study is 374 ± 59 kBq min, which is comparable with the 353 ± 59 kBq min result from the original study, with the resulting K_i values of the NLR and Patlak analyses differing slightly, but with an acceptable inter-observer difference (NLR results: 1.3%; Patlak results: 6.3%).

Reproducible values between the original data analysis of this and the original study (2) are a prerequisite for being able to attribute subsequently observed differences in kinetic parameters as being due to this study's use of “optimized” input functions while keeping all other parameters (PWR, Bg-VOIs, and Bone-VOIs) almost constant in all series.

This assumption is supported by previous repeatability studies on the use of ¹⁸F-NaF-PET to evaluate SUV as a measure for uptake in cancer patients (11, 12), as well as ¹⁸F-NaF-PET being previously used successfully for studies on changes in K_i values before and after treatment of osteoporosis (4, 7, 13).

Comparing the VOI volumes and their corresponding RCs, it is obvious there is a correlation between these parameters, as the RCs are derived from the activity in the VOIs and compared with the actual blood sample activities (Table 1): the smaller the VOI, the closer to 1 the correction factor, except for the AO-Peak VOI in which the factors are much higher than 1. This correlation has been described previously in cardiac FDG-PET studies (14) and can be ascribed to a combination of PVE and an increasing mean activity (kBq/ml) in the VOIs as the smaller volumes are placed, more selectively, over areas with the highest activity.

The RC_ß corrects the activities measured in the input VOI for both PVE and Bg (3), whereas the RC_CF only corrects for PVE. The RC_ß (0.99 ± 0.09) and RC_CF (0.99 ± 0.09) are identical for the AO-Fix-ß VOIs, indicating the background in and around the aortic wall is negligible, at least at this level of the thoracic aorta. These factors are comparable to the coefficients reported in a study by Puri et al. (15). However, for the AO-Peak VOIs, the correction factors are unrealistically high (1.73 ± 0.63) and result in a much larger variance in the derived K_i results when compared with the other corrected input functions.

As the correction factors were primarily lower than 1, the general effect on the uncorrected input functions was to increase the AUCs as shown in Supplementary Figure S1. The closer the correction factor is to 1, as in the AO-Fix curves, the lesser the shift in values.

The exception to this was the AO-Peak curves which, due to the correction factors of 1.73 ± 0.63, resulted in smaller AUCs and larger variance. The most likely reason for the poor AO-Peak VOI performance is due to its inherent poor counting statistics caused by the small VOI volume of only a few voxels in combination with the amount of activity injected and the short time resolution in the initial dynamic acquisition bins.

The general effect of substituting the final image exponentials with plasma exponentials was smaller IQRs, as seen in the box plots in Supplementary Figure S2. However, the AUC mean differences were not significantly different across the series of input functions (Supplementary Section 1.3, Table S3).

The substitution methods producing AUC results (kBq/ml min) closest to the basic corrected IDIFs were the LV-New-Pl-40-90-Log (Δ-Mean: −1.17, 95% Cl −17.7 to 15.3) and the LV-New-40-90-Exp (Δ-Mean: 3.2, 95% Cl −9.7 to 16.2). For comparison, the corresponding values for the AO-Fix-40-90-Exp were much higher (Δ-Mean: −11.0, 95% Cl −44.4 to 22.5; p = 0.04).