A Projection-Domain Low-Count Quantitative SPECT Method for α-Particle-Emitting Radiopharmaceutical Therapy

Single-photon emission-computed tomography (SPECT) provides a mechanism to estimate regional isotope uptake in lesions and at-risk organs after administration of α-particle-emitting radiopharmaceutical therapies (α-RPTs). However, this estimation task is challenging due to the complex emission spectra, the very low number of detected counts (~20 times lower than in conventional SPECT), the impact of stray-radiation-related noise at these low counts, and the multiple image-degrading processes in SPECT. The conventional reconstruction-based quantification methods are observed to be erroneous for α-RPT SPECT. To address these challenges, we developed a low-count quantitative SPECT (LC-QSPECT) method that directly estimates the regional activity uptake from the projection data (obviating the reconstruction step), compensates for stray-radiation-related noise, and accounts for the radioisotope and SPECT physics, including the isotope spectra, scatter, attenuation, and collimator–detector response, using a Monte Carlo-based approach. The method was validated in the context of 3-D SPECT with 223Ra, a commonly used radionuclide for α-RPT. Validation was performed using both realistic simulation studies, including a virtual clinical trial, and synthetic and 3-D-printed anthropomorphic physical-phantom studies. Across all studies, the LC-QSPECT method yielded reliable regional-uptake estimates and outperformed the conventional ordered subset expectation-maximization (OSEM)-based reconstruction and geometric transfer matrix (GTM)-based post-reconstruction partial-volume compensation methods. Furthermore, the method yielded reliable uptake across different lesion sizes, contrasts, and different levels of intralesion heterogeneity. Additionally, the variance of the estimated uptake approached the Cramér–Rao bound-defined theoretical limit. In conclusion, the proposed LC-QSPECT method demonstrated the ability to perform reliable quantification for α-RPT SPECT.


II. THEORY
Consider a SPECT system imaging a radioisotope distribution f(r), where r = (x, y, z) denote the spatial 3-D coordinates. Denote the measured projection data by the M-dimensional vector g. Assume that the object being imaged and the projection data lie in the Hilbert space of square-integrable functions, denoted by L 2 ℝ 3 and the Hilbert space of Euclidean vectors, denoted by E M , respectively. Then, the SPECT system, denoted by the operator ℋ, is a transformation from L 2 ℝ 3 to E M . In SPECT with α-RPTs, the stray-radiation-related noise occupies a substantial portion of the measured counts due to the very low count levels. We model this noise as Poisson distributed with the same mean ψ for all projection bins. Let Ψ be an M-dimensional vector, with each element equaling ψ, that denotes the mean of the stray-radiation-related noise across all M projection bins. Denote the entire noise in the imaging system by the M-dimensional random vector n. Then, the projection data g are Poisson distributed with mean ℋf + Ψ. Thus, the imaging system equation is given as follows: g = ℋf + Ψ + n . (1) Our objective is to estimate the regional uptake within a set of VOIs. Mathematically, we first define a 3-D VOI function ϕ k VOI (r), where ϕ k VOI (r) = 1, if r lies within the k th VOI 0, otherwise. (2) Denote λ as the K-dimensional vector of regional uptake, where λ k is given as follows: λ k = ∫ d 3 rf(r)ϕ k VOI (r) ∫ d 3 rϕ k VOI (r) .
Our objective is to estimate λ.

A. Reconstruction-Based Quantification Methods
The conventional procedure to estimate λ is to first reconstruct the activity uptake distribution over a voxelized grid, and then estimate the activity uptake in a discretized version of the VOI as defined in (2). Mathematically, the activity uptake distribution is described using a voxel basis function, denoted by ϕ n VOX (r), as follows: f vox (r) = ∑ n = 1 N θ n ϕ n vox (r) .
Multiple algorithms are available to estimate the coefficients θ n . These algorithms yield an estimate of θ n , denoted by θ n . To estimate λ k , we first define a discretized mask matrix M.
One procedure to define the elements of this matrix is given as follows: M n, k = 1, if a majority of the n th voxel lies inside the k th VOI . 0, otherwise. (5) Then, the estimate of λ k obtained from the reconstructed image, denoted by λ k recon , is given as follows: λ k recon = 1 N ∑ n = 1 N θ n M n, k . (6) This procedure to estimate λ has several issues. First, a large number of voxels need to be estimated during reconstruction, leading to a highly ill-posed problem, especially when the number of counts is low. This may, in turn, lead to a biased estimate [26]. A second issue is the bias introduced due to partial-volume effects (PVEs) [27]. PVEs arise due to two distinct phenomena namely, finite system resolution and tissue-fraction effects. When defining M n,k , an element of this matrix is 1 when a majority of this voxel is within the VOI. Therefore, this does not account for tissue fraction effects, causing bias when estimating λ k recon using (6). The third issue is the activity inside a voxel θ n is fundamentally not estimable [28]. Next, by the data-processing inequality, the process of reconstruction can only lead to information loss [17], [29]. Finally, these RBQ approaches are often based on maximum-likelihood expectation maximization (MLEM) [30] or ordered subset expectation maximization (OSEM) [31]. However, at low counts, these methods have limited precision and deviate from the theoretically lowest possible Cramér-Rao bound (CRB) [26]. All these issues serve as sources of error in the estimated regional activity. Thus, as previous studies have reported, even highly fine-tuned versions of these methods yield unreliable estimates of regional uptake [13]- [15].

B. Proposed Method
To address the above-described issues with RBQ methods, we recognize that our objective is to estimate the mean uptake within certain regions. Thus, we directly represent the object f(r) in terms of the VOI basis functions. These VOI basis functions are given by ϕ k VOI (r) as defined in (2). The activity distribution is then represented in terms of these basis functions as follows: f VOI (r) = ∑ k = 1 K λ k ϕ k VOI (r) (7) where, if the activity inside each VOI is constant, then f VOI (r) = f(r). Inserting this definition for f(r) in (1) yields the following expression for the mth element of the vector g g m = ∫ ℎ m (r)f(r)d 3 r + ψ + n m = ∑ k = 1 K λ k ∫ ℎ m (r)ϕ k VOI (r)d 3 r + ψ + n m . (8) This can be written in vector form as follows: g = Hλ + Ψ + n (9) where H is the M × K-dimensional system matrix with elements given by Given the measured projection g, to estimate λ, we maximize the probability of occurrence of the measured data. Denote Pr(x) as the probability of a discrete random variable x. Then, the probability of the measured projection data is given as follows: where we have used the fact that the measured data across the different bins are independent. Now, the measured data g m are Poisson distributed with mean (Hλ) m + ψ. Thus This gives the likelihood of the measured data g. To estimate λ, we maximize the logarithm of the likelihood of λ given g λ = arg max λ ln[Pr(g | λ)] . (13) To maximize this log likelihood, we follow the same process as used to derive the conventional MLEM technique [17]. Briefly, we differentiate the log likelihood with respect to the elements of λ and equate that to 0 to find the point at which the log likelihood is maximized. This yields the following iterative equation to estimate λ k : where λ k (t) denotes the estimate of λ k at the tth iteration. We refer to this procedure as the LC-QSPECT method. Our approach advances on existing methods that directly quantify from projection data [18] by providing the ability to model stray-radiation-related noise. As we see later, this ability to model stray-radiation-related noise plays a key role in the task of reliable quantification. Furthermore, the system matrix H models all key image-degrading processes in α-particle SPECT using a Monte Carlo (MC)-based approach, which further improves the performance of this technique on the task of quantification.
Our approach alleviates the issues outlined earlier with RBQ approaches. Typically, the number of VOIs K is less than the number of voxels N, so the problem is less ill-posed.
Besides, the method is less sensitive to PVEs since we define the boundaries of VOIs before estimating the regional uptake. In particular, the tissue-fraction effects are minimized since there is no voxelization. Furthermore, while it is true that the mean VOI activity λ k could also be inestimable, but since the VOI is generally larger than a voxel, the estimation bias is lower. We also directly estimate the regional uptake from the projection data, thus avoiding any reconstruction-related information loss. Finally, as our results later show, the method yields estimates with a precision that is close to the CRB.

A. Implementation
Implementing the proposed LC-QSPECT method required obtaining the elements of the system matrix H mk , as we see from in (14). We obtained these elements using an MCbased approach. More specifically, we used SIMIND, a well-validated MC-based simulation software [32], [33] to model the isotope emission and the SPECT system. Next, for a given patient, we obtained the definition of the VOIs. These can be obtained, for example, by segmenting the CT that is acquired along with the SPECT. We assigned unit uptake to the VOI and zero uptake elsewhere. We assumed that the attenuation map of the patient is available. Projection data for this activity and attenuation map were generated by simulating more than 100 million photons for each VOI. The simulations modeled all relevant imagedegrading processes in SPECT, including attenuation, scatter, collimator response, septal penetration, and scatter, characteristic X-ray from both the α-emitting isotope and the lead in the collimator, finite energy and position resolution of the detector, and the backscatter in the detector. Scaling the projection data according to the acquisition time of the projections yielded the corresponding columns of the system matrix. While the bremsstrahlung was not modeled in this study, this does not significantly affect the accuracy of the simulation, as we show in Section IV-A.
Next, the LC-QSPECT method required obtaining the mean of stray-radiation-related noise, i.e., ψ in (14). We estimated this experimentally from a planar blank scan acquired on the SPECT system for over 10 min. Averaging the projection bin counts in this scan yielded the mean background counts. This was then scaled to the acquisition time to estimate the mean stray-radiation-related noise.
The computed system matrix and mean stray-radiation-related noise were used in (14) to estimate the regional uptake directly from the projection data. As the system matrix modeled all relevant image-degrading processes, these processes were automatically compensated during quantification.

B. General Evaluation Framework
Evaluating the performance of the LC-QSPECT method on the estimation task of regionaluptake quantification required a setup where the ground-truth regional uptake was known.
For this purpose, we conducted realistic simulations, including a virtual clinical trial (VCT), and physical-phantom studies. We describe these evaluations in detail in Section III-C. In this section, we describe the methods to which we compared our method and the figures of merit.

1) Methods Compared:
We compared the performance of the LC-QSPECT method with two widely used RBQ methods.

a) OSEM-based method:
Here, the activity maps were first reconstructed using an OSEM-based approach. This approach, implemented using the customizable and advanced software for tomographic reconstruction (CASToR) [34] software, compensated for attenuation, scatter, collimator-detector response, and stray-radiation-related noise. Scatter was compensated using the triple-energy-window (TEW) method [35]. To minimize noise amplification that could be caused by using the TEW method, we applied pre-reconstruction Butterworth filters to the photopeak and scatter-window projections [36]. These filters were optimized by minimizing the normalized root-mean-square error (NRMSE) between the true and estimated regional-uptake values. The cutoff frequencies of the filters applied to photopeak and scatter-window projections were optimized by a 2-D grid search and the optimized frequencies were found to be 0.15 and 0.05 cycle/pixel, respectively. The optimized orders were found to be 8 for both filters. We also compensated the strayradiation-related noise with an additive term in the iteration equation of the OSEM-based method. This additive term was similar to that in the proposed method and is as described in Section III-A. The reconstructed-image dimensions were 128 × 128 × 91, with a voxel side length of 4.418 mm. We also optimized the number of iterations and subsets based on the NRMSE between the true and estimated regional uptake. The optimized number of iterations and subsets were found to be 20 and 6, respectively, which was consistent with that reported in [15]. From the reconstructed image, the uptake in different VOIs was calculated.
b) GTM-based method: PVEs are known to degrade quantification accuracy in SPECT [27]. The LC-QSPECT method implicitly assumes constant uptake within each VOI. Under this assumption, PVEs can also be compensated post-reconstruction. A widely used approach for this purpose is the geometric transfer matrix (GTM)-based method [37]. Thus, we also compared our approach to this method. The elements of the GTM were obtained from the projections and reconstructions of the VOIs. These projections were obtained using the process described in Section III-A and the reconstruction was done using the OSEM-based method as described above. The rest of the implementation of this method was as described in [37].

2) Figures of Merit:
We evaluated the accuracy, precision, and overall error of the LC-QSPECT, OSEM, and GTM-based methods on the task of estimating regional uptake. In most of our experiments, we generate multiple instances of projection data for a single phantom, where each instance corresponds to a separate noise realization. Denote the total number of noise realizations by R. Denote the true and estimated activity uptake in the kth VOI for the rth noise realization by λ rk and λ rk , respectively. In these experiments, the accuracy of the estimated uptake was quantified using the normalized bias (NB), which, for the kth VOI, is given as follows: The precision of the estimated uptake was quantified using the normalized standard deviation (NSD), which, for the kth VOI, is given as follows: Finally, the overall error in estimating the uptake was quantified by the NRMSE. For the kth VOI, this is given as follows: To evaluate the performance of the methods over populations, we used the ensemble NB and ensemble NRMSE. Denote the number of samples in the population by S. Denote the true and estimated activity uptake in the kth VOI for the sth sample by λ sk and λ sk , respectively. The ensemble NB for the kth VOI is given as follows: The ensemble NRMSE for the kth VOI is given as follows: Next, in cases where we had just a single noise realization, we used normalized error to quantify performance. This is defined as the difference between the true and estimated uptake values, normalized by the true uptake value.
Finally, we also computed the CRB, which is the minimum variance that can be achieved by an unbiased estimator, as a benchmark to compare the precision of the activity estimated using the proposed method. The CRB is given by the diagonal elements of the inverse of the Fisher information matrix for the estimated parameter [17]. Denote this matrix by F. Denote λ k and λ k as the true and estimated activity uptake, respectively, in the kth VOI of a targeted digital phantom. Then, the elements of this matrix are given as follows: . (20) Substituting (12) in (20) yields

C. Evaluation Using Realistic Simulation Studies
We conducted realistic simulation studies in the context of imaging patients with bone metastases of prostate cancer treated with 223 Ra dichloride, a widely administered United States (U.S.) Food and Drug Administration (FDA)-approved α-RPT indicated for the treatment of patients with castration-resistant prostate cancer [38], [39]. A major site of these osseous metastases is the pelvis, which is adjacent to the anatomical location of the prostate. Thus, we focused on the pelvic region.
To generate a realistic patient population for our study, digital 3-D activity and attenuation maps of the pelvic region were generated using the extended cardiac-torso (XCAT) [40] phantom (representative slice shown in Fig. 1). To simulate continuous activity distribution, the activity map had a high resolution of 512 × 512 along the axial dimensions and 364 slices along the depth dimension. The side length of the voxels was 1.105 mm. Our analysis of clinical SPECT data of patients administered 223 Ra dichloride therapy suggested that there were three primary sites of uptake in the pelvic region: 1) lesion (indicated by the arrow); 2) bone; and 3) gut (through which the majority of administered activity is cleared). The rest of the region in the patient typically had the same low uptake. We refer to the VOI corresponding to the rest of the patient as the background. This led to a total of four VOIs. The variability in anatomies and regional uptake within these phantoms were simulated based on clinical data, as described later (Section III-C3).
Next, a GE Optima 640 SPECT system with a high energy general purpose (HEGP) collimator was simulated using SIMIND. The scintillation detector in the system had an intrinsic spatial resolution of 3.9 mm and an energy resolution of 9.8% at 140 keV, where the resolutions were quantified in terms of full width at half maximum (FWHM). Photons were acquired at 60 angular positions spaced uniformly over 360°. The photopeak window was set as 85 keV ± 20% [41]. The isotope emission and all relevant image-degrading processes in SPECT were simulated. The mean of the stray-radiation-related noise in each projection bin, ψ, was determined as described in Section III-A. The value of stray-radiation-related counts in each projection bin was individually sampled from a Poisson distribution with a mean equal to ψ and was added to the corresponding projection bin in the simulated projections. The generated projections had around 5000 counts per axial slice, simulating a clinically realistic low-count α-RPT SPECT acquisition. We used the procedure described in Section III-A and III-B to quantify the uptake from the simulated projection data using the proposed, OSEM, and GTM-based methods. Then, the methods were evaluated on the task of quantifying the regional uptake using the figures of merit defined in Section III-B2.

1) Validating the SPECT Simulation:
We conducted our simulations with SIMIND. While SIMIND has been validated for multiple SPECT studies [33], [42], we further validated the accuracy of our simulation approach in the context of α-RPT SPECT. For this purpose, we compared the projection data obtained with our simulation approach to that obtained on an actual scanner. More specifically, we considered a NEMA phantom that was scanned on a GE Optima 640 SPECT system with the HEGP collimator using the procedure as described in more detail in Section III-D. We then modeled this acquisition using our simulation approach. For our simulation, the activity map was designed to simulate the known 223 Ra activity concentrations filled in the physical phantom and the attenuation map of the NEMA phantom was derived from the CT scans. Then, we generated the simulated projections using SIMIND, modeling the same acquisition process as described later in Section III-D1. The profiles of the projection data at four angular positions spaced uniformly over 360° obtained with the scanner and that with the simulation approach were compared directly, without any normalization. The match of these profiles, as we see later in the results section, provided evidence of the accuracy of our simulations.

2) Evaluating Convergence of the LC-QSPECT Method:
For this purpose, we generated five 3-D XCAT-based phantoms with dimensions similar to average patient size. Each phantom had a lesion of a different size but at the same location in the pelvis. The lesion diameters varied between 15 and 35 mm. The activity uptake in the four VOIs was in the ratio of 2:5:25:20 in the background, bone, gut, and lesion regions. These ratios were derived based on our analysis of clinical data. Projection data corresponding to these phantoms were generated as described in Section III-C. We applied the LC-QSPECT method to the generated projections and the error in the estimated lesion uptake after each iteration was computed. Six hundred iterations were performed for each phantom.

3) Evaluating Performance in Simulated Clinical Scenario Using Virtual
Clinical Trial: VCTs are an emerging evaluation paradigm that provide the ability to rigorously and objectively evaluate the performance of new imaging technology in simulated clinical scenarios that model patient population variability [43]. In our VCT, we simulated 50 digital 3-D male patients with different anatomies using the XCAT phantom. As in [40], the heights and weights of the 50 patients were sampled from a Gaussian distribution with the mean equal to the height and weight of a 50th percentile male U.S. adult and a 10% standard deviation. Next, based on clinically derived parameters [12], the lesion diameter was sampled from a Gaussian distribution with a mean of 33.75 mm and a standard deviation of 12.64 mm. Subsequently, as described above, the activity uptake in the four VOIs was sampled from a normal distribution with a clinically derived mean uptake ratio of 2:5:25:20 in the background, bone, gut, and lesion regions. Projection data corresponding to these patients were generated as described in Section III-C. The LC-QSPECT, OSEM, and GTM-based methods were applied to these data, thus comparing the performance of these methods for this realistic population. We evaluated the accuracy and overall error of the regional-uptake estimates yielded by the LC-QSPECT method across all these regions for this clinically realistic patient population.

4) Comparing Precision of Estimated Regional Uptake Using LC-QSPECT
Method With Cramér-Rao Bound: From our results, we observed that the LC-QSPECT method was yielding approximately unbiased estimates of the regional activity uptake. For an unbiased estimator, the minimum variance that can be obtained is given by the CRB. Thus, we compared the precision of the estimated regional uptake with this bound.
To conduct this study, we generated 50 noise realizations for a 3-D male patient simulated using XCAT. The LC-QSPECT method was used to estimate the mean regional activity in each of the 50 noise realizations. The NB and NSD of these estimates were computed. The NSD of estimates in each VOI was compared to the square root of the CRB.

5) Evaluating Performance as Function of Lesion Size and Contrast:
To evaluate the sensitivity of the method to lesion size, we used the same simulation setup as described in Section III-C2. Fifty noise realizations were generated for each of the five phantoms, where each phantom had a different lesion size. To evaluate the sensitivity of the method to lesion contrast, we recognize that the lesion is present within the bone. We generated six XCAT-based phantoms with average patient size and a pelvic lesion of diameter 33.75 mm [12]. Each phantom had a different lesion-to-bone uptake ratio (LBUR), ranging from 1:1 to 6:1. The uptake in the background, bone, and gut regions was in the ratio of 2:5:25. Fifty noise realizations were generated for each phantom and used to evaluate the method performance as a function of lesion contrast.

6) Evaluating Effect of Spatial Intralesion Heterogeneity:
In contrast to the OSEM-based method, the LC-QSPECT method assumes that the activity uptake in the VOI is homogeneous. This assumption may not hold in clinical settings [44], [45]. Thus, we evaluated the impact of spatial intralesion heterogeneity on the performance of the LC-QSPECT method.
We generated five phantoms with average patient size and a lesion of diameter 33.75 mm [12] in the pelvis. Each phantom had a different amount of spatial intralesion heterogeneity.
To simulate intralesion heterogeneity, we modeled the uptake in the lesion as a 3-D lumpy model [46]. Denote the support of the lesion in the object space by s(r). Then, the lesion activity uptake, denoted by f l (r), is given as follows: where P denotes the total number of lumps, and c p , a p , and σ p denote the center, magnitude, and width of the pth lump function, respectively. Different levels of heterogeneity were simulated by varying the values of P, c p , a p , and σ p . The heterogeneity was characterized by entropy, where a higher value of entropy refers to more heterogeneity in uptake. We used Shannon entropy [47], a commonly used metric to quantify entropy. As described in [48], we first calculated the histogram of the activity map of the lesion region. The histogram had 256 bins {b w , w = 0, 1, … , 255}. Denote v as the total number of voxels in the lesion region. The normalized histogram {B w , w = 0, 1, … , 255} was calculated using the expression B w = b w /v. Then, the entropy, denoted by E, was computed as follows [48]: The generated lesions and the corresponding entropy are shown in Fig. 2. All lesions had the same mean activity uptake. The mean activity uptake in the background, bone, gut, and lesion region was in the ratio of 2:5:25:20. We generated 50 noise realizations for each phantom. Using these data, we evaluated the accuracy and precision of the proposed method as a function of different levels of spatial intralesion heterogeneity.

7) Evaluating Effect of Compensating for Stray-Radiation-Related Noise:
A key feature of the proposed method is the ability to compensate for stray-radiation-related noise. To evaluate the impact of compensating for this noise on quantification performance in α-RPT, a comparative test was performed using the VCT setup. The LC-QSPECT method was modified to not compensate for stray-radiation-related noise by setting the term ψ = 0 in (14). This was compared to the proposed LC-QSPECT method that compensated for this noise.

D. Evaluation Using Physical-Phantom Studies
Evaluation with physical-phantom studies quantifies the performance of the method with real scanner data. We conducted this study with two phantoms: 1) NEMA phantom (Data Spectrum ™ , USA) [ Fig. 3(a)] and 2) 3-D printed vertebrae phantom [ Fig. 3(b)]. The NEMAphantom study was conducted to evaluate the performance of the LC-QSPECT method for different lesion sizes, with the spheres in the phantom simulating lesions. The vertebraephantom study was conducted to simulate the imaging of a lesion within the spine bone in the thoracic region. Details on the designing, printing, and preparation of this phantom for the experiments are provided in the supplementary material.

1) Phantom Scanning:
We filled both the phantoms with clinically relevant 223 Ra activity concentrations as described in the supplementary material. Next, we scanned the NEMA phantoms on a GE Optima 640 SPECT/CT system with a HEGP collimator and the vertebrae phantom in the same system with a medium energy general purpose (MEGP) collimator, with the goal of evaluating the robustness of our method for different collimator configurations. We placed each phantom at the center of the field of view of the γ camera.
The photopeak and scatter windows were set as 85 keV ± 20% and 57 keV ± 18%, respectively. Scans were acquired at 60 angular positions spaced uniformly over 360°. The acquisition time at each angular position was set to 30 s, as in clinical studies. The size of the projection at each angular position was 128 × 128 pixels, where the pixel side length was 4.4 mm. A body-contour orbit was used to improve resolution. Corresponding low-dose CT scans were also acquired for each phantom (120 kVp, 10 mA, and 512 × 512). The axial pixel spacing of the CT images was 0.98 mm and the spacing between slices was 5.0 mm. The CT and SPECT scans were registered.

2) Regional-Uptake Estimation:
The NEMA phantom had seven VOIs, including the six spheres and the background. The vertebrae phantom had two VOIs, namely, the lesion chamber and the background. VOI definitions were obtained by manually segmenting fused CT scans and OSEM reconstructed SPECT images of both phantoms. The segmented VOI masks had a 5 mm distance between adjacent slices, which was the same as the CT scans. We performed a cubic interpolation along the depth dimension to generate VOI masks with a 1 mm distance between adjacent slices, yielding the definition of VOI boundaries along the depth dimension. Using these VOIs, the system matrices for the LC-QSPECT method were generated using the process as described in Section III-A. From the projection data, the regional uptake was estimated using the LC-QSPECT method. In addition, we compared the performance of the LC-QSPECT method on the task of estimating the regional uptake of the NEMA phantom with and without compensating for the stray-radiation-related noise.
Next, the activity image was reconstructed on the clinical workstation XELERIS (General Electric, USA) using the scanner-based OSEM technique with the following parameters: two iterations, ten subsets, and a Butterworth filter with cutoff frequency of 0.48 cycle/cm and order of 10. Attenuation and scatter were compensated using the CT-based attenuation map and the dual-energy window (DEW) scatter-compensation method, respectively, [49]. From the reconstructed image, the regional uptake was estimated, yielding the output with the OSEM-based method. Finally, the GTM-based method was applied to the reconstructed image, using the process described in Section III-B1. We computed the normalized difference between the estimated and true regional uptake, termed as the normalized error, for all three methods.

A. Validating the SPECT Simulation Procedure
Projections of the NEMA phantom from the simulated and physical SPECT systems are shown in Fig. 4. We compared the projections at four angular positions spaced uniformly over 360°. For each angular position, the profiles along the dashed line in both projections are also shown. For each point in the profile, the number of counts was calculated by averaging among eight adjacent pixels on both sides of the dashed line. The averaging was performed to reduce the noise-related variation of the profile. We observe that the profiles of the simulated projection match those acquired on the scanner at all angular positions. This provided evidence for the accuracy of the SPECT system simulation procedure used in this manuscript.

B. Convergence of the Proposed Technique
The normalized error in the estimated lesion uptake as a function of the iteration number of the LC-QSPECT method is shown in Fig. 5. We observed that the method converged after 256 iterations for all five lesion diameters. Thus, we chose 256 as the number of iterations for the LC-QSPECT method for subsequent experiments.

1) Virtual Clinical Trial:
The absolute ensemble NB and ensemble NRMSE of the estimated regional uptake in the lesion, gut, bone, and background (BKGD) in the VCT setup are shown in Fig. 6. Also provided is a violin plot that shows the distribution of the normalized error between the true and estimated uptake for all 50 patients using the proposed method. We observe that for all the regions, the LC-QSPECT method consistently outperformed the OSEM and GTM-based methods, based on both accuracy and overall error. The LC-QSPECT method had at least three times lower ensemble bias in the estimated lesion uptake compared to the OSEM and GTM-based methods. Furthermore, the absolute ensemble NB obtained with the LC-QSPECT method for all the regions was consistently lower than 3.4%. Finally, for 90% of the simulated patients, the normalized error in estimating the lesion uptake with the proposed method was within ± 20%. Table I shows the NB of the estimated regional uptake from 50 noise realizations for a representative patient using the LC-QSPECT method.

2) Comparing Precision of Estimated Regional Uptake Using LC-QSPECT Method With Cramér-Rao Bound:
The NB values are all close to 0%, with the maximum value being below 4%. Based on these and other results as presented later, we observe that the LC-QSPECT method is approximately unbiased. Thus, in Fig. 7, the NSD of the corresponding estimated regional uptake using the LC-QSPECT method was compared with that obtained by taking the square root of the normalized CRB, the lower bound of variance for any unbiased estimator. We observe that for all the regions, the proposed method yielded an NSD very close to that obtained from the CRB.

3) Performance as Function of Lesion Size and Contrast:
The absolute NB, NSD, and NRMSE of the estimated lesion uptake as a function of the lesion diameter using the different methods are shown in Fig. 8. The NSD of the estimates is also compared with the square root of the CRB. The LC-QSPECT method consistently yielded close to zero bias. Furthermore, the NSD of the estimates obtained with this method was close to that derived from the CRB for all lesion diameters. Additionally, the method yielded the lowest NRMSE and consistently outperformed both the OSEM and GTM-based methods for all lesion sizes.
The absolute NB and NSD of the estimated lesion uptake as a function of the LBUR using the LC-QSPECT, OSEM, and GTM-based methods are shown in Fig. 9. We also compared the NSD of the estimates with that derived from the CRB. Again, the LC-QSPECT method consistently yielded almost zero bias and NSD close to the square root of CRB for all LBUR values. Furthermore, the method consistently outperformed both the RBQ methods.

4) Impact of Spatial Intralesion Heterogeneity:
The absolute NB and NSD of the estimated lesion uptake for different degrees of spatial intralesion heterogeneity, as quantified with the entropy parameter, using the LC-QSPECT, OSEM, and GTM-based methods are shown in Fig. 10. We observe that the LC-QSPECT method consistently yielded close to zero bias for all degrees of heterogeneity. Furthermore, even the precision of the estimated mean uptake was not significantly impacted by the intralesion heterogeneity.
Finally, based on both the accuracy and precision, the proposed method outperformed both the RBQ methods.

5) Impact of Compensating for Stray-Radiation-Related Noise:
The ensemble NRMSE in estimating regional uptake using the LC-QSPECT method with and without compensating for the stray-radiation-related noise in the VCT is shown in Table II. We observe that compensating for the stray-radiation-related noise led to significantly more reliable regional-uptake estimates in this clinically realistic α-RPT SPECT setup.

1) NEMA-Phantom Study:
The normalized absolute error in estimating the regional uptake as a function of the sphere diameter is shown in Fig. 11. We observe that the LC-QSPECT method consistently outperformed both the RBQ methods. Next, the normalized error of the LC-QSPECT method with and without compensating for the stray-radiationrelated noise is shown in Table III. We observe that the values of absolute normalized error in the largest three spheres of the NEMA phantom are at least three times higher than that of the LC-QSPECT that compensates for the stray-radiation-related noise.

2) Vertebrae-Phantom Study:
The normalized absolute error in estimating the lesion uptake using the LC-QSPECT, OSEM, and GTM-based methods is shown in Table IV. Again, the LC-QSPECT method significantly outperformed the other two approaches. Furthermore, the LC-QSPECT method yielded an error as low as 2.86% in these challenging real-world conditions.

V. DISCUSSION
We have designed, implemented, and evaluated an LC-QSPECT approach to quantify regional activities from low-count SPECT data for α-RPTs. Our results demonstrate the efficacy of the LC-QSPECT method for α-RPTs. In Fig. 6, we observed that the LC-QSPECT method yielded averaged absolute ensemble NB and ensemble NRMSE values as low as 1.6% and 5.2%, respectively. In contrast, the OSEM-based method yielded averaged absolute ensemble NB and ensemble NRMSE of 28.2% and 32.1%, respectively, and the GTM-based method yielded averaged absolute ensemble NB and ensemble NRMSE of 12.6% and 18.7%, respectively. These observations for the OSEM-based method are consistent with previous literature [13]- [15] and confirm the limitation of this approach for quantitative SPECT with α-RPTs. Furthermore, our results demonstrate the reliability and superiority of the proposed method for quantitative SPECT with α-RPTs.

Figs. 7, 8(b), and 9(b)
show that the NSD of the regional-uptake estimates obtained using the LC-QSPECT method approached the lowest theoretical limit as defined by the CRB. This is an important finding as quantifying uptake precisely is a major challenge in α-RPT SPECT [15]. Furthermore, our results in Figs. 6(a), 8(a), and 9(a), and Table I indicate that the LC-QSPECT method is approximately unbiased. These findings, while empirical, are theoretically consistent as the LC-QSPECT method is an ML estimator, and when an efficient estimator exists, the ML estimator is efficient [17]. Thus, our results suggest that the LC-QSPECT method may be the optimal estimator in terms of bias and variance properties. In Fig. 8(b), we observe that for certain lesion sizes the RBQ methods yielded a lower NSD than the proposed method and that derived from CRB. However, we also observe that the RBQ methods are biased for all the lesion sizes considered in this experiment, so the comparison of the variance of these methods with CRB is not meaningful. Furthermore, we see that even for those lesion sizes, the proposed method yielded an improved overall accuracy and precision as quantified by the NRMSE.
The efficacy of the LC-QSPECT method was also observed for different lesion sizes (Fig. 8) and different LBURs (Fig. 9). The bias in the estimated uptake was close to zero for all the considered lesion sizes and LBUR values. In particular, the low bias at different lesion sizes demonstrates that the LC-QSPECT method is relatively insensitive to PVEs. PVEs are a major source of error in quantitative SPECT. Note that even though the GTM-based method is designed exclusively to compensate for PVEs, the method still yielded large values of bias for small lesions. In contrast, the proposed method was approximately unbiased (Fig.  8). This result thus demonstrates an important advantage of this method compared to OSEM and GTM-based methods. Similarly, the proposed method was approximately unbiased for different LBUR values, demonstrating the accuracy of the method for different signal contrast values.
While the LC-QSPECT method assumes that the activity uptake in each VOI is homogeneous, the results in Fig. 10 show that the intralesion heterogeneity considered in this study does not substantially affect the reliability of the estimated uptake. Furthermore, very importantly, even when such heterogeneity is present, the proposed method outperforms the OSEM and GTM-based approaches. Thus, even when the assumption of uniform uptake in the lesion is violated to some extent, these results suggest that the LC-QSPECT method can still be a reliable and superior option.
In the physical-phantom studies, we again observed that the LC-QSPECT method significantly outperformed both the RBQ methods ( Fig. 11 and Table IV). These results were consistent with the simulation-study results. More importantly, they demonstrate the feasibility of the proposed method in real-world settings. The LC-QSPECT method requires the definitions of VOIs, and the physical-phantom studies show that these can be obtained from segmenting fused CT and reconstructed SPECT images. CT scans are typically acquired in conjunction with the SPECT scans for attenuation compensation, and often using SPECT/CT systems.
compensating for this noise resulted in over 180% NRMSE. We observed that the proposed LC-QPSECT method effectively compensated for this stray-radiation-related noise. Similar results were observed in the physical-phantom study where compensating for the strayradiation-related noise improved the estimation accuracy of the proposed method for the four larger spheres in NEMA phantom, as shown in Table III. For the smallest two spheres, however, the absolute values of normalized error without stray-radiation-related noise compensation using the LC-QSPECT method are slightly lower than those with this compensation. Here, we note that for the LC-QSPECT method, VOIs are defined from manual segmentation of the low-dose CT scans. These VOIs may have errors. For the small sphere sizes, even a small misdefinition of the VOI may cause a large error. Furthermore, these results are only for a single noise realization and the impact of stochastic noise-related errors is high for smaller spheres. Thus, it is likely that, for these smaller spheres, the overestimation caused due to the stray-radiation-related noise was canceled by the underestimation due to the noise-related stochastic error and the potential misdefinition of the VOI maps. However, the average values of absolute normalized error among all spheres in the NEMA phantom using the LC-QSPECT with and without compensating for the stray-radiation-related noise were 29.2% and 46.3%, respectively. This demonstrated that overall, the method yielded much improved performance.
Another important feature of the LC-QSPECT method is the use of a MC simulation-based approach to generate the system matrix. This approach yields highly accurate modeling of SPECT physics. The MC approach is computationally feasible because the number of VOIs is typically quite small. Thus, this matrix can be precomputed and stored. In this study, the system matrix of each patient took less than 50 min and 30 MB for generation and storage, respectively. In clinical applications, the CT scans of the patient could be acquired first, then the system matrix can be generated simultaneously when the SPECT scan is acquired. Furthermore, the estimation process is fast; and 256 iterations required < 30 s on a standard desktop computer equipped with an Intel Core i7-10700K CPU with 16 cores and 32-GB RAM. This is unlike developing such an approach for OSEM-based methods, for which a similar system matrix may require up to 30 TB of memory. Thus, the proposed method provides a mechanism for highly accurate MC-based system modeling, which is not possible with RBQ methods. This is another advantage of the proposed method.
A limitation of the LC-QSPECT method is that the method quantifies regional uptake rather than voxel values. Thus, dosimetry estimates are confined to the mean absorbed dose across the region. This may then preclude the estimation of a dose-volume histogram within tumors and normal organs, which is used in the external-beam radiation therapy. However, given the errors with voxel-based reconstruction methods, the estimation of dose-volume histograms may be infeasible at these low counts. In contrast, our method shows that the mean regional uptake can be estimated reliably. Another point to note is that the exchange of (purported) resolution from conventional voxel-based to this region-based method is of limited consequence to α-RPT dosimetry as no imaging system can resolve isotope distribution at the sub-100 μm scale.
The LC-QSPECT method, similar to the RBQ methods, requires reliable VOI definitions. While these VOI definitions could be obtained from the CT scans that are acquired for quantitative SPECT, there is a possibility that the SPECT and CT scans may be misaligned. Studying the effect of this misalignment is an important research direction. Methods to register SPECT and CT scans [50], [51] would help address this limitation. Another area of future research is the clinical validation of these techniques. One issue with such validation is the absence of ground-truth quantitative values. To address this issue, no-gold-standard evaluation techniques are being developed [52], including in the context of evaluating quantitative SPECT reconstruction methods [53]. This may provide a mechanism to clinically evaluate the proposed methods.

VI. CONCLUSION
An LC-QSPECT method was proposed for quantitative SPECT with α-RPTs, where the number of detected γ-ray photons is very small. The method yielded reliable (accurate and precise) values of regional uptake and outperformed the conventional OSEM and GTM-based methods, as evaluated in the context of α-RPT with 223 Ra. The method yielded reliable activity estimates in a VCT simulating imaging of patients with bone metastasis who were administered with 223 Ra dichloride therapy. Additionally, the method yielded reliable estimates of uptake for different lesion sizes and different LBUR values. Furthermore, the method yielded reliable lesion uptake estimates for different degrees of intralesion heterogeneity. The method was observed to be approximately unbiased and yield a standard deviation close to that defined by the CRB, indicating that the method may be an efficient estimator. Evaluation with physical-phantom studies on SPECT/CT scanners using NEMA and an anthropomorphic vertebrae phantom reinforced this reliable quantification performance and provided evidence for the feasibility of the approach in practical settings. Overall, the results provide strong evidence for further evaluation and application of this method to perform quantitative SPECT for α-RPTs.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. Digital (a) activity map and (b) attenuation map for the pelvic region generated using the anthropomorphic XCAT phantom. Projections of the NEMA phantom and the profiles along the yellow-dashed lines in the projections acquired at the four angular positions using the scanner and simulated SPECT system.  Normalized error in the estimated activity uptake in the lesion region with five different diameters as a function of the iteration number of the LC-QSPECT method.  Comparing the NSD of the regional activity estimates obtained from the LC-QSPECT method with the lower bound for the NSD as defined by the CRB.    (a) Absolute NB and (b) NSD of the estimated lesion uptake as a function of different amounts of intralesion heterogeneity as quantified using the entropy parameter. Normalized absolute error of regional-uptake estimates of the NEMA phantom using the LC-QSPECT, GTM, and OSEM-based methods.