Designing an optimized protocol fordetecting transient dopamine releasebyConnor William James BevingtonB.Sc., The University of Waterloo, 2017A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2019c© Connor William James Bevington 2019The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:Designing an optimized protocol for detecting transient dopamine releasesubmitted by Connor William James Bevington in partial fulfillment of the requirements forthe degree of Master of Sciencein PhysicsExamining Committee:Vesna Sossi, Physics and AstronomySupervisor Arman Rahmim, Physics and AstronomyAdditional Examiner (Second Reader)iiAbstractA single positron emission tomography (PET) scan can be used to detect voxel-level transientdopamine release in response to mid-scan task or stimulus. A standard approach uses F-test(significance) thresholding followed by a cluster size threshold to limit false positives. The Fstatistic is computed from two neurophysiological kinetic models, one able to handle dopamine-induced perturbations to the time-varying radiotracer concentration, and the other unable to doso. Through simulation we first demonstrate that extensive denoising of the dynamic PET images isrequired for this method to have high detection sensitivity, though this often leads to a large clustersize threshold—limiting the detection of smaller cluster sizes—and poorer parametric accuracy.Our simulations also show that voxels near the peripheries of clusters are often rejected—becomingfalse negatives and ultimately distorting the F-distribution of rejected voxels. We then suggest anovel Monte Carlo method that incorporates these two observations into a cost function, allowingerroneously-rejected voxels to be accepted under specified criteria. In simulations, the proposedmethod boosts sensitivity by up to 77% while preserving cluster size threshold, or up to 180%when optimizing for detection sensitivity. A further parametric-based voxelwise thresholding isthen suggested to better estimate the dopamine release dynamics in detected clusters. Finally,we apply the Monte Carlo method coupled with the parametric thresholding approach to a pilotscan from a human gambling study, where additional parametrically-unique clusters are detectedas compared to the current best method.iiiLay SummaryThis work improves upon current detection algorithms searching for dopamine release in the brain,using positron emission tomography (PET) imaging. To acquire PET images tracking dopamine,we inject a molecule that competes with dopamine for binding to receptors in the brain, attached toa radioactive element. The radioactive decay products generate a time-varying signal, allowing usto follow the molecule as it travels throughout the brain. We zoom in on tiny volumes in the brain,and statistically analyze the signal to detect dopamine release. The algorithms developed in thiswork borrow theory used when simulating magnetic materials in thermal equilibrium; though theapplications are rather different, the mathematics are quite similar. Applying these algorithms tosimulated and real data, we demonstrate improvements over the current best methods. These algo-rithms can now be applied to full-scale dopamine release studies, which may improve understandingof brain function in healthy and diseased people.ivPrefaceProfessor Vesna Sossi of the UBC PET group conceptualized the project goal of improving uponcurrent dopamine release detection algorithms. I was responsible for developing the novel theoryof the Monte Carlo method described in Chapter 2. Simulation design was conducted by me, withinput from Dr. Ju-Chieh (Kevin) Cheng and Dr. Ivan S. Klyuzhin of the UBC PET group. I wasresponsible for all analysis and interpretations of the simulations in this work.Part of the work in Section 4.1 and Section 4.2 was presented at the 2018 IEEE Medical ImagingConference, and published in the conference record [1]. The theory presented in Section 2.7 andpart of the results in Section 4.4 and Section 4.6 is included in a manuscript current in preparation[2]. In both cases, I performed all simulations, data analysis, and wrote the manuscripts.Section 4.6 contains data and analysis from a single human pilot scan for a gambling study, usingthe methods developed in this work. The gambling study was designed by Professor Catharine A.Winstanley from the UBC Department of Psychology and Dr. Mariya V. Cherkasova from the UBCFaculty of Medicine, Division of Neurology. I analyzed the data from the pilot scan and interpretedthe results presented in this work.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 PET: fundamental concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 PET physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.1 Positron emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.2 Scanner geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.3 Gamma ray detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.4 Scattering and attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.5 Random coincidences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.6 Other physical limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3.1 Acquiring data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3.2 Analytic reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.3 Iterative reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3.4 Ordered Subsets Expectation Maximization reconstruction . . . . . . . . . . 142.3.5 The system matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.4 Denoising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.1 Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.2 HYPR denoising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.4.3 Incorporating HYPR within image reconstruction . . . . . . . . . . . . . . . 192.5 Kinetic modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.5.1 Compartmental modelling: plasma input models . . . . . . . . . . . . . . . . 222.5.2 Compartmental modeling: reference tissue models . . . . . . . . . . . . . . . 242.5.3 Simplified Reference Tissue Model . . . . . . . . . . . . . . . . . . . . . . . . 25viTable of Contents2.5.4 Linear parametric neurotransmitter PET model . . . . . . . . . . . . . . . . 282.6 Detecting DA release . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.7 A novel Monte Carlo method for detecting DA release . . . . . . . . . . . . . . . . . 362.7.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.7.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.7.3 Optimizing the MC method . . . . . . . . . . . . . . . . . . . . . . . . . . . 393 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1 Generating a simulated phantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1.1 Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1.2 Regional fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1.3 PSF filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.1.4 Forward-projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.1.5 Noise addition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.2 Simulation variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.3 Types of simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3.1 Baseline simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3.2 Random release magnitude simulation . . . . . . . . . . . . . . . . . . . . . . 433.3.3 Cluster release simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.4 Reconstruction and denoising protocols . . . . . . . . . . . . . . . . . . . . . . . . . 463.4.1 Reconstructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.4.2 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.4.3 Comparison to previous noise-modeling methods . . . . . . . . . . . . . . . . 473.5 Analysis workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.5.1 Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.5.2 Weighting schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.5.3 Thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5.4 Determining the cluster size threshold . . . . . . . . . . . . . . . . . . . . . . 513.5.5 Assessing detection sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5.6 Parametric accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.5.7 Comparing DA release curves . . . . . . . . . . . . . . . . . . . . . . . . . . 533.5.8 Parametric thresholding: A more representative picture of release dynamics 533.6 Application to real human data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.1 Baseline simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.1.1 Cluster size thresholds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.1.2 Comparing results using previous noise-modeling methods . . . . . . . . . . 594.2 Random release magnitude simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 604.2.1 Sensitivity as a function of release level . . . . . . . . . . . . . . . . . . . . . 604.2.2 Parametric accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.2.3 Validity of the cluster size threshold . . . . . . . . . . . . . . . . . . . . . . . 724.2.4 Comparing results using previous noise-modeling methods . . . . . . . . . . 734.2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3 Cluster release simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.3.1 Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.3.2 Parametric accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.3.3 Comparing results using previous noise-modeling methods . . . . . . . . . . 91viiTable of Contents4.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 944.4 Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 954.4.1 Choosing the optimal hyperparameter β and thresholding value ccrit . . . . . 964.4.2 Baseline simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 994.4.3 Random release magnitude simulation . . . . . . . . . . . . . . . . . . . . . . 994.4.4 Cluster release simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1044.5 Parametric thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1114.6 Application to a human pilot study . . . . . . . . . . . . . . . . . . . . . . . . . . . 1175 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1205.1 Improved noise-modeling over previous work . . . . . . . . . . . . . . . . . . . . . . 1205.2 Variables in simulation design, image processing, and analysis . . . . . . . . . . . . 1215.3 Baseline simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.4 Random release magnitude simulation results . . . . . . . . . . . . . . . . . . . . . . 1215.5 Cluster release simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1225.6 Take-aways from analyzing the simulations using the standard method . . . . . . . 1235.7 Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1235.8 Parametric thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1245.9 Application to a human pilot study . . . . . . . . . . . . . . . . . . . . . . . . . . . 1245.10 Final thoughts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126viiiList of Tables1.1 Motivating questions for DA release studies . . . . . . . . . . . . . . . . . . . . . . . 23.1 Cluster properties in the cluster release simulation . . . . . . . . . . . . . . . . . . . 453.2 Reconstruction and denoising protocols tested in this work . . . . . . . . . . . . . . . 484.1 CST for a 16 mCi bolus dose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.2 CST for a 20 mCi bolus dose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.3 CST for a 24 mCi bolus dose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.4 Cluster size thresholds used for analyzing random release magnitude and clusterrelease simulations, as well as real human data. . . . . . . . . . . . . . . . . . . . . . 594.5 Global sensitivity in the random release magnitude simulation (16 mCi) . . . . . . . 614.6 Global sensitivity in the random release magnitude simulation (20 mCi) . . . . . . . 614.7 Global sensitivity in the random release magnitude simulation (24 mCi) . . . . . . . 614.8 False positive cluster-free rate in the random release magnitude simulation . . . . . . 734.9 Cluster characteristics in the cluster release simulation . . . . . . . . . . . . . . . . . 774.10 Voxelwise sensitivity in the 16 mCi cluster release simulation (traditional weighting) 784.11 Voxelwise sensitivity in the 16 mCi cluster release simulation (uncorrected weighting) 784.12 Binary50 sensitivity in the 16 mCi cluster release simulation (traditional weighting) 784.13 Binary50 sensitivity in the 16 mCi cluster release simulation (uncorrected weighting) 784.14 Binary sensitivity in the 16 mCi cluster release simulation (traditional weighting) . . 794.15 Binary sensitivity in the 16 mCi cluster release simulation (uncorrected weighting) . 794.16 Comparing voxelwise sensitivity in the cluster release simulation between the stan-dard and MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1044.17 Comparing binary50 sensitivity in the cluster release simulation between the stan-dard and MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1054.18 Comparing binary sensitivity in the cluster release simulation between the standardand MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105ixList of Figures2.1 Coincidence events during a PET scan . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Parameterizing a LOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3 Representation of a point source in sinogram space. . . . . . . . . . . . . . . . . . . . 102.4 Spatial comparison of filtering vs. HYPR post-processing . . . . . . . . . . . . . . . 202.5 Temporal comparison of filtering vs. HYPR post-processing . . . . . . . . . . . . . . 202.6 Visualizing tracer exchange in SRTM . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.7 Modulating parameters in the gamma variate function, h(t) . . . . . . . . . . . . . . 302.8 lp-ntPET TACs for a variety of release patterns . . . . . . . . . . . . . . . . . . . . . 312.9 Examples of time-varying PET-estimated % DA release . . . . . . . . . . . . . . . . 322.10 Reducing DAR(t) to DAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.11 Lookup curves for converting between γ and DAR . . . . . . . . . . . . . . . . . . . 332.12 Distribution of false positive cluster sizes . . . . . . . . . . . . . . . . . . . . . . . . . 352.13 F-distributions for true baseline, true release, and rejected voxels . . . . . . . . . . . 373.1 Comparing noiseless simulated images, with and without PSF filtering . . . . . . . . 423.2 Voxel-level DAR in the random release magnitude simulation . . . . . . . . . . . . . 443.3 Distribution of voxel-level DAR before and after PSF filtering . . . . . . . . . . . . . 443.4 Cluster locations and voxel-level DAR in the cluster release simulation . . . . . . . . 463.5 A relative comparison of weighting schemes . . . . . . . . . . . . . . . . . . . . . . . 504.1 Comparing false positive cluster size distributions in OSEM-5mm-5mm: 20 mCi vs24 mCi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.2 Voxelwise sensitivity as a function of DAR (16 mCi) . . . . . . . . . . . . . . . . . . 624.3 Effect of weighting scheme used in the random release magnitude simulation . . . . . 634.4 The effect of dose level in the random release magnitude simulation . . . . . . . . . . 654.5 Bias in DAR(t) introduced by post-processing in the random release magnitudesimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.6 The effect of post-processing on TACs in the random release magnitude simulation . 664.7 The effect incorrect basis functions on TACs in the random release magnitude sim-ulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.8 Mean timing parameter estimates in the random release magnitude simulation . . . 694.9 MAE in timing parameters for the random release magnitude simulation . . . . . . . 694.10 The effect of dose level on MAE of the timing parameters in the random releasemagnitude simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.11 BP (t) in the random release magnitude simulation . . . . . . . . . . . . . . . . . . . 714.12 DAR(t) in the random release magnitude simulation . . . . . . . . . . . . . . . . . . 714.13 Voxelwise sensitivity as a function of DAR for the Gaussian noise addition method . 734.14 Mean timing parameter estimates in the Gaussian noise simulations . . . . . . . . . 744.15 MAE in timing parameter estimates in the Gaussian noise simulations . . . . . . . . 74xList of Figures4.16 BP (t) for the Gaussian noise addition method in the random release magnitudesimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.17 DAR(t) for the Gaussian noise addition method in the random release magnitudesimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.18 Effect of weighting scheme on detection sensitivity in the cluster release simulation . 804.19 Effect of dose level on detection sensitivity in the cluster release simulation . . . . . 814.20 Bias in DAR(t) introduced by post-processing in the cluster release simulation . . . 834.21 Bias in DAR(t) introduced by not having the corrected basis function in the clusterrelease simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.22 Mean estimates of timing parameters in the cluster release simulation . . . . . . . . 864.23 Effect of dose level on mean timing parameter estimates in the cluster release simulation 864.24 MAE in timing parameters for the cluster release simulation . . . . . . . . . . . . . . 874.25 Effect of dose level on MAE of timing parameters in the cluster release simulation . 874.26 BP (t) estimates in the cluster release simulation . . . . . . . . . . . . . . . . . . . . 884.27 DAR(t) estimates in the cluster release simulation . . . . . . . . . . . . . . . . . . . 894.28 Demonstrating the difference between intra- and extra-striatal bias . . . . . . . . . . 914.29 Comparing voxelwise sensitivity from simulations using Gaussian noise addition . . . 914.30 Comparing binary50 sensitivity from simulations using Gaussian noise addition . . . 924.31 Comparing binary sensitivity from simulations using Gaussian noise addition . . . . 924.32 Mean timing parameter estimates in the cluster release simulation using Gaussiannoise addition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 934.33 Timing parameter MAE in the cluster release simulation using Gaussian noise addition 934.34 Optimizing the MC method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 974.35 Demonstrating MC cost function equilibrium . . . . . . . . . . . . . . . . . . . . . . 984.36 Sensitivity as a function of DAR in the random release magnitude simulation for theMC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.37 Visualizing voxelwise sensitivity improvements in the random release magnitude sim-ulation using the MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.38 Mean estimated timing parameters in the random release magnitude simulation forthe MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.39 Distribution of estimated timing parameters in the random release magnitude sim-ulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.40 MAE in timing parameters for the random release magnitude simulation using theMC method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.41 BP (t) in the random release magnitude simulation using the MC method . . . . . . 1034.42 DAR(t) in the random release magnitude simulation using the MC method . . . . . 1044.43 Visualizing voxelwise sensitivity improvements in the cluster release simulation usingthe MC methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1064.45 Distribution of estimated timing parameters in the cluster release simulation . . . . 1074.44 Mean estimated timing parameters in the cluster release simulation using the MCmethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1074.46 MAE in timing parameters in the cluster release simulation using the MC methods . 1084.47 BP (t) in the cluster release simulation using the MC methods . . . . . . . . . . . . . 1094.48 DAR(t) in the cluster release simulation using the MC methods . . . . . . . . . . . . 1104.49 Distribution of estimated timing parameters in the cluster release simulation afterparametric thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1124.50 BP (t) in the cluster release simulation using the parametric thresholding approach(high release clusters) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113xiList of Figures4.51 BP (t) in the cluster release simulation using the parametric thresholding approach(low release clusters) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1144.52 DAR(t) in the cluster release simulation using the parametric thresholding approach(high release clusters) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1154.53 DAR(t) in the cluster release simulation using the parametric thresholding approach(low release clusters) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1164.54 Analysis results from a human pilot scan . . . . . . . . . . . . . . . . . . . . . . . . . 1184.55 Demonstrating the effect of PTA on real human data . . . . . . . . . . . . . . . . . . 119xiiList of AcronymsAD absolute deviationB/I bolus plus constant infusionBGO bismuth germanateBLUE best linear unbiased estimatorBP binding potentialCDF cumulative distribution functionCOV coefficient of variationCST cluster size thresholdDAR dopamine release levelDA dopamineFBP filtered backprojectionFDG [18F]fluorodeoxyglucoseFOV field of viewFWHM full width at half maximumGSO gadolinium oxyorthosilicateHRRT High Resolution Research TomographHYPR HighlY constrained backPRojectionLOR line of responselp-ntPET linear parametric neurotransmitter PET modelLSO lutetium orthosilicateMAD median absolute deviationMAE mean absolute errorMC Monte CarloMLEM Maximum Likelihood Expectation MaximizationMRI magnetic resonance imagingxiiiList of AcronymsOSEM Ordered Subsets Expectation MaximizationPDF probability distribution functionPET positron emission tomographyPMT photomultiplier tubePSF point spread functionPTA parametric thresholding approachPVE partial volume effectRAC [11C]racloprideROI region of interestSiPM silicon photomultipliersSNR signal-to-noise ratioSRTM Simplified Reference Tissue ModelTAC time activity curveTOF time-of-flightWRSS weighted residual sum of squaresxivAcknowledgementsThis work was funded by an NSERC Discovery grant (240670-13), an NSERC Create grant (448110-2014), and a CIHR grant (PJT-156012). I would like to thank my supervisor, Professor Vesna Sossi,and research mentor, Dr. Ju-Chieh (Kevin) Cheng, for their guidance throughout this project andhigh expectations to bring out my best work. I would also like to thank my parents, Sue and Mark,for raising me in a stimulating environment for curiosity and learning. I would not be where I amtoday if not for your diligence and support during my formative years. Finally, I owe gratitude tomy high school physics teacher, Alastair Paterson, who sparked my interest in physics that burnsever-strong to this day.xvChapter 1IntroductionPositron emission tomography (PET) is an in vivo functional imaging modality that uses radi-olabeled molecules to trace physiological processes in the body. PET is frequently used in neu-roimaging, studying a variety of processes including metabolism, abnormal protein aggregation,and neurotransmitter activity [3].One specific application is the detection of dopamine (DA) release in response to a pharma-ceutical intervention or neuropsychological task [4–7]. PET is able to indirectly image DA releaseusing the radiolabeled compound [11C]raclopride (RAC), comprised of the radioisotope 11C andthe dopamine D2 receptor antagonist raclopride.Being an antagonist, a RAC molecule will compete with DA for binding to D2 receptors, thusfunctioning as a proxy for free D2 receptor density [8, 9]. Being labeled with a β+ emitter, a RACmolecule will also generate local radioactive decays which are subsequently detected by a PETscanner to determine its location in the brain.When DA is released above basal levels in a particular region, there will be increased DA bindingand reduced RAC binding, resulting in a decrease in RAC signal in that region. To characterizethis decrease through neurophysiological models, one must perform dynamic (4D) PET imaging: aseries of 3D images quantifying RAC concentration as a function of time. This provides the basisfor indirect DA release imaging.Traditional protocols involve two separate scans to quantitatively assess DA release, commonlyreferred to as a “double RAC” protocol. One scan images the brain of the participant in thebaseline state—that is without any stimulus that would induce DA release above basal levels. Asecond scan attempts to induce a constant DA release level above the basal level. Mathematically,the time course of concentration during each scan can be compared to arrive at a PET-estimatedpercent dopamine release level (DAR).This traditional protocol suffers from three main issues:• Logistics and safety: the participant must be scheduled for two separate scans, each givingthem radiation exposures on the order of 2-3 mSv.• Scan-to-scan variations: an innumerable set of factors contribute to basal DA levels, in-cluding mood [10], level of fatigue [11], level of hunger [12], and placebo [13]. Therefore it islikely that the basal DA levels will differ in the two scans, confounding DAR.• Steady-state assumption: In the second scan, an assumption is made that the dopaminericsystem is in steady-state, meaning DA concentrations are at a constant value above basallevels. This may be violated [14].Relatively recently, a single scan protocol has been suggested [15–17] to circumvent these issues,which we shall call the mid-scan intervention protocol. Here, a single scan is performed with theparticipant in the baseline scan for the first portion of the scan, to allow for tracer uptake into thebody. Following this, an intervention is performed to induce a transient DA release, which can betraced out by analyzing the dynamic signal.1Chapter 1. IntroductionStudy Question Metric(s)Where is DA release located? detection sensitivityWhat is the strength of DA detection sensitivity, signal bias,release in each voxel? parametric accuracyWhat time does DA release sensitivity, parametric accuracystart/peak/end?What is the smallest detectable cluster size thresholdregion of DA release?Metric descriptionsDetection sensitivity: what percentage of the time is true DArelease detected?Signal bias: what is the error between the true and detectedsignal?Parametric accuracy: what is accuracy in the parameters ofthe neurophysiological models?Cluster size threshold: what is the smallest detectable clustersize that has a 90% probability of not being a false positive?Table 1.1: Likely questions asked by a DA release study, and what metric(s) to use to decide on an optimaldenoising technique.With only a single scan performed, the logistical issues are mitigated and the radiation dosegiven to the participant is halved compared to the double RAC protocol. As well, the scan-to-scan variations are no longer a confounding factor. Most importantly, however, the steady-stateassumption is no longer necessary: in fact, we know the intervention will cause a transient response,consisting of a rise from basal DA levels to some peak level, then an eventual decay back to basallevels.If the dynamic signal is analyzed at the smallest volume possible (i.e. at the voxel-level), one canin theory map out the entire spatiotemporal DA release pattern in the brain. However, voxel-levelsignals are very noisy, and thus necessitate sophisticated denoising techniques to allow accurateand consistent quantification of DA release.The first aim of this work is to assess the effect of a variety of commonly-used denoising tech-niques on the ability to detect voxel-level DA release. This shall be studied through a few lenses,since the goals of particular DA release studies may differ. Some questions that a DA release studymay be interested in addressing are listed in Table 1.1, along with suggested metric(s) that shouldbe used to identify the optimal denoising technique to answer these questions.The second aim of this work suggests a new, Monte Carlo-based approach to modify the existingtransient signal detection algorithm. We shall show that this new algorithm provides very noticeableincreases in detection sensitivity while maintaining parametric accuracy and relatively low falsepositive rates.To address these aims, we first design realistic simulations of RAC scans based on real humandatasets. We then validate the results from the simulated data with a real human scan which ispart of a pilot study involving transient DA release detection.The structure of this work is as follows: In Chapter 2, we shall briefly introduce the physicsbehind PET signal acquisition and image reconstruction, followed by mathematical descriptionsof the denoising operators to be tested. Next, the theory behind neurophysiological modelingis addressed, both in general and pertaining to the specific application of transient DA release2Chapter 1. Introductiondetection.In Chapter 3, we first explain the best current detection algorithm. The implementation of thevarious denoising techniques is then discussed. A description then follows of the various simulationsused to test this algorithm—after the different denoising techniques have been applied—and themetrics to be used to assess their efficacy. Next, a formulation of the novel Monte Carlo method ispresented. We then finish with a description of the human pilot study used to validate the resultsfrom the simulations.The results are then presented from the various simulations for each denoising technique andthe novel Monte Carlo method in Chapter 4, as well as for the single human pilot study scan. Adiscussion alongside the results suggests the optimal denoising technique to use for a particular DArelease study, based on the goals of such a study.3Chapter 2Theory2.1 PET: fundamental conceptsBy way of analogy, PET imaging is similar to mapping out the density of airplanes in the sky bycounting contrails. Instead of a direct measurement of airplane density, this is an indirect measure,making use of a signature left behind by the airplane: the contrails. In the case of PET imaging,we are mapping out radiolabeled molecule densities in the body, and the signature left behind istheir radioactive decays.However, in PET imaging we are often a further step removed from the target of interest. Thisis best explained through another example. Say we wish to image metabolism in the brain. Theprimary energy source for the brain is glucose, and so this would serve as our target. However,instead of radiolabeling a glucose molecule, we instead synthesize a surrogate, deoxyglucose, andsubstitute a hydrogen atom with 18F to form [18F]fluorodeoxyglucose (FDG). Uptake of FDGmolecules by glucose transporters will still occur as if it were a glucose molecule. The radioactivedecays of 18F are then what we detect. The number of decays within a specific region serve as aproxy for FDG uptake, which in turn serves as a proxy for glucose uptake.In order than one does not perturb the physiological process we wish to study, the concentrationof radiolabeled molecules introduced into the body is extremely small—on the order of pmol to nmol.Thus the ratio of proxy molecules to target molecules (in our example: the ratio of FDG moleculesto glucose molecules) is extremely small, and the introduced molecule acts as a tracer that does notperturb the state of the system under investigation. Refining our analogy, PET imaging is morelike counting contrails to determine the density of passengers in the sky: contrails are a proxy forairplane density, which in turn is a proxy for the density of passengers in the sky.2.2 PET physics2.2.1 Positron emissionThe process for PET imaging, from a physics perspective, is as follows: a molecule that serves asa surrogate for a specific physiological process is selected, and then labeled with a radioisotopeto form a radiolabeled molecule, or radiotracer. The radioisotope used must undergo β+ decay,meaning a proton in its nucleus is converted to a neutron, releasing a positron (to conserve charge)and an electron neutrino (to conserve lepton number)p→ n+ e+ + νe (2.1)Following a β+ decay event, the positron will travel a short probabilistic distance characterizedby a quantity known as the positron range (10−2 to 10−1 cm [18]), before annihilating with anelectron. This length defines the theoretical resolution of PET imaging, since even if we couldlocalize the annihilation event with perfect accuracy, we still could only localize the location of thepositron emission to within the positron range, on average.42.2. PET physicsThe electron-positron annihilation, provided it occurs at rest in the frame of the body, willproduce two collinear 511 keV gamma rayse+ + e− → γ + γ (2.2)In reality, however, the annihilation will in general not occur at rest, resulting in some non-collinearity (∼ 0.5◦). This further degrades the resolution of PET imaging.2.2.2 Scanner geometryThe two (approximately) collinear gamma rays then travel towards an array of detectors that formpart of a PET scanner. Generally, the detectors are arranged in parallel rings to form a cylindricalenclosure around the participant being scanned. However, octangonal geometries [19] as well asarrays of detectors conforming to the shape of the head [20–22] have also been used.The job of the detectors is to identify coincidence events. Say that an electron-positron annihi-lation occurs somewhere along an imaginary line connecting a detector pair, and—for simplicity—such an annihilation event occurs at the radial centre of the field of view (FOV). Then the pathstraveled by the two gamma rays will always be equal and length. Since they both travel at thespeed of light, they will arrive at their respective detectors at an identical time. This coincidencedetection defines a line of response (LOR), that is the path along which the annihilation eventoccurred, provided neither gamma ray has been scattered.If the annihilation event occurs away from the radial centre of the FOV, then the distancetraveled by one gamma ray will be shorter than the other, resulting in their arrival times beingoffset by short duration. As well, limitations in detector electronics result in a nonzero temporalresolution in detecting arrival time. Therefore coincidence events are identified within a shortcoincidence window, ∆τ , typically on the order of nanoseconds [18]. A detector pair will accept acoincidence event when a signal characteristic of a 511 keV gamma ray is received in each detectorwithin ∆τ . Such events are referred to as prompts.Without any further information, the location of the annihilation event can only be assignedequal probability along the LOR connecting the detector pair. However, with higher timing reso-lution (on the order of hundreds of picoseconds), one can begin to assign a higher probability toa particular segment of the LOR. This type of imaging is known as time-of-flight (TOF) imaging[23]. In TOF, the LORs are divided into bins each representing a Gaussian probability distributionin position, and prompts are assigned to a particular bin depending on their coincidence timingdifference.Consider the following example: if one gamma ray travels 10 cm to a detector, and its pairedgamma ray travels 20 cm to another detector, then the difference in path length is ∆` = 10 cm,which corresponds to a detection timing difference of ∆t = ∆`/c = 330 ps. if we divide a 30 cmLOR into five bins of 6 cm length, then this translates to a bin width of 200 ps in time. A moreuseful way of describing this would be the central bin represents ∆t ∈ ±[−100, 100] ps, the adjacentbins represent |∆t| ∈ [100, 300] ps, and the outermost bins represent |∆t| ∈ [300, 500] ps. Given∆t = 300 ps in our example, the prompt will then be assigned with the highest probability to theoutermost bin closest to the detector that registered the first prompt.2.2.3 Gamma ray detectionIn most currently-operating PET scanners, the detectors consist of an array of crystals coupled toa photomultiplier tube (PMT) [18]. New generation scanners, as well as PET inserts for magneticresonance imaging (MRI) scanners and hybrid PET/MR scanners, use silicon photomultipliers52.2. PET physics(SiPM) instead of PMTs, due to their functionality within high magnetic field as well as high gainat lower voltage settings than PMTs [24, 25].Broadly, the role of the crystal is to produce scintillation events when absorbing a 511 keVgamma ray. This requires high density and atomic number (and thus high probability of capturinga gamma ray), short crystal dead time (the recovery time for a lattice site between producing ascintillation event and being able to produce another), high energy resolution (for discriminating511 keV events), and high light output (one incidence gamma ray produces several scintillationphotons) [26]. Commonly used crystal materials are bismuth germanate (BGO), lutetium orthosil-icate (LSO), and gadolinium oxyorthosilicate (GSO) [18, 27].Once a gamma ray strikes a crystal, lower energy scintillation photons are generated (generallyin the UV or visible range) and move towards the PMT or SiPM. A single incident scintillationphoton will generate an avalanche of photoelectrons via the photoelectric effect which are detectedas a current spike. The number of photoelectrons produced is proportional to the energy of theincident gamma ray, and thus in principle allows for energy discrimination based on the pulse heightin the signal [18, 28].In an ideal scenario, the system will respond linearly with increasing count rates. Thereforewith double the count rate, we expect the number of events detected by the scanner to also double.However, the effect of dead time (τd) causes a deviation in this expected linear response. The deadtime is a function of the crystal type, detector electronics, and count rate [28].Within the detector electronics, an integration time exists where the circuits integrate over thelight pulse produced from scintillation events. At high count rates, the probability of a detectorreceiving multiple impinging photons within this integration time increases. When this happens,the detector circuitry integrates over multiple pulses. This either exceeds the maximum allowableenergy used to discriminate for 511 keV events (thus all events are lost), or the multiple eventsare treated as a single event (thus the positions of the true events are misplaced) [28]. Anotherpotential contributor to dead time is the “reset time” in detector electronics, where following ascintillation event, a small amount of time must pass before the detector may accept another event.The net effect is that the system responds sublinearly to increasing count rate: at higher truecount rates, the measured count rate will be lower than the ground truth [28]. Corrections for deadtime thus become increasingly important at high count rates. A review of correction methods isbeyond the scope of this work, but a summary is presented in [28]. In general, a combination ofanalytic and measured correction factors are used.2.2.4 Scattering and attenuationAs the gamma rays leave the annihilation site, they must pass through layers of tissue and perhapsbone, and then air, before reaching a detector. A photon has the potential to interact with theelectronic and atomic structures of these impediments, resulting in scattering events.For scattering, the 511 keV photons have the highest cross-section of interaction for the Comptoneffect [18, 28]. In the Compton effect, the photon strikes an outer (quasi-free) electron of an atom.The photon changes path and loses energy as some is passed to the electron, which recoils at anangle to the photon to conserve energy and momentum.This poses several problems for PET imaging: the change in direction can result in an incorrectLOR being assigned to the event, known as a scattered coincidence, again reducing the accuracy ofthe image. The detectors in a PET scanner—in addition to timing windows—also discriminate byenergy [18, 28]. Therefore in large angle scattering events, the energy loss could place the energyof the photon outside this window, preventing the event from being recorded at all. Additionally, ascattered photon initially within the FOV may be directed outside the FOV after scattering, again62.2. PET physicspreventing the event from being detected.In all above cases, we denote the original annihilation event as being attenuated—regardless ofwhether the scattered event is detected along an incorrect LOR, or if it goes completely undetected.The combination of missed events and incorrectly-assigned annihilation event locations, if notaccounted for, will add bias to the final reconstructed image. Thus these effects need to be correctedfor (as much as possible) during the process of image reconstruction.2.2.5 Random coincidencesA high number of counts is imperative to accurately mapping out the radiotracer distribution inPET imaging. Depending on the activity of the tracer during a given time of the scan, ∼ 104 to∼ 106 counts are being recorded each second [18]. A high count rate leads to the possibility of arandom coincidence, where two photons emitted at a similar time but from separate annihilationsites are mistakenly registered as coming from the same annihilation site [18, 28]. This is possiblebecause of the finite coincidence timing window.A closely-related issue is multiple coincidences, where more than two photons are detected withinthe detector timing window. In this case it is impossible to determine which LOR to assign, andso the multiple events that produced these detections are missed, or a single event is kept basedon some arbiration process [18, 28]. Diagrams illustrating true, scatter, random, and multiplecoincidences are presented in Figure 2.1.2.2.6 Other physical limitationsBy design, the radiotracers used in imaging ideally have high activities to produce a high numberof counts per second. However, all radioisotopes have a half-life, and thus their activity decreasesas the scan goes on. For instance, the radiotracer RAC has a half-life of 20.4 minutes, and hence ascan comprises several half-lives. Increasing the initial activity of the radiotracer evidently increasesthe total number of counts in a scan, but also increases the radiation dose delivered to the subjectduring the scan. Regulatory agencies set strict limitations on radiation doses that can be deliveredto subjects volunteering for medical research studies, thus limiting the activities of the radiotracersused in PET imaging.Detector geometry can also pose problems during the imaging process. LORs that pass throughdetectors at more oblique angles have a possibility of producing a scintillation event in adjacentcrystals. As mentioned, some scanner designs also are not composed of circular rings of detectors.For instance, the High Resolution Research Tomograph (HRRT) is composed of octagonal rings.Therefore gamma ray pairs that would require LORs defined through the gaps in detectors at thevertices of the octagon will not be detected.The detectors themselves may also scatter an incoming gamma ray [29]. If the gamma rayis scattered by the initial crystal it strikes to a nearby crystal (where it is subsequently fullystopped), then the event will be slightly mispositioned. If the gamma ray remains within the initialcrystal, but is not completely stopped, then the number of scintillation events will be lower, andhence the detected energy of the gamma ray will be underestimated—potentially leading to eventrejection. While, from a physics perspective, such events may be labeled as scattering events, whenmodeling the imaging process we consider these events under the area of detector normalizationand resolution modeling, since the scattering is occurring within the crystals rather than in theimaging FOV. We elaborate on these issues in Section 2.3.5, when discussing the system matrix inimage reconstruction.72.2. PET physics(a) true coincidence (b) scatter coincidence(c) random coincidence (d) multiple coincidencesFigure 2.1: Various coincidence events that can occur during a PET scan. Annihilation events are indicatedby stars. Solid lines indicate the paths of the gamma ray pairs, and dashed lines indicate the LOR(s) recordedas a result of the events.Finally, the subjects imaged are often human beings, which are not able to remain completelymotionless during the scan. This will introduce blurring in the image without motion correction,since the same annihilation site will map to different voxels in image space as the patient moves.82.3. Image reconstruction2.3 Image reconstruction2.3.1 Acquiring dataWhen two prompts are recorded within the coincidence window ∆τ , an event can be added tothe LOR connecting the two detectors. Mathematically, LORs are represented within a sinogrambased on the mathematics of the Radon transform. Consider first a 2D object f(x) = f(x, y). TheRadon transform maps f(x) to a function S(L) which exists in the space of straight lines, L, in R2p(L) =∫Lf(x)|dx| (2.3)For our application, the set L is not infinite, but rather naturally is the set of all LORs. Asa starting point, consider only the LORs that line in the plane that defines f(x). We can theparametrize a LOR, L(s, θ), in the xy-plane asx cos θ + y sin θ − s = 0 (2.4)where s is the distance of L from the origin and θ is the angle of the normal vector of L makeswith the x-axis (cf. Figure 2.2). Equation (2.3) then becomesp(s, θ) =∫ ∞−∞∫ ∞−∞f(x, y)δ(x cos θ + y sin θ − s)dxdy (2.5)where the delta function makes the integrand only nonzero when integrating along the line L(s, θ).To easily visualize this process, consider a single point with unitary magnitude at a location x =(a, b) with its position vector making an angle α with the x-axis. In this case, the set of nonzeroLORs all intersect this single point, and hence can be thought to “rotate” about the point as θtravels through [−pi/2, pi/2). Mathematically, in sinogram spaces = |x| cos(α− θ) (2.6)where |x| = √a2 + b2. Since α and |x| are fixed for all LORs, as θ travels through [−pi/2, pi/2)a sinusoid is traced out in sinogram space. Therefore, points in R2 translate to sinusoids in asinogram (Figure 2.3).To bring this back into the context of PET imaging, consider a positron-emitting point sourcelocated at x = (a, b) with zero positron range (for simplicity). With each positron annihilation,a pair of nearly collinear gamma rays will be emitted. For those pairs that travel within the2D imaging plane along a line at an angle θi to the x-axis, a value of one is assigned to S(s =|x| cos(α− θi), θ = θi).In the limit of infinite emissions and infinitesimal detector resolution, the full range of θ ∈[−pi/2, pi/2) will be sampled, and thus in sinogram space a full sinusoid will be traced out. However,in a real scan there will be finite emissions. As well, the detectors have nonzero surface area andtherefore put a finite limit on the number of LORs that can be defined. As a result the sinogramis discretized and finitely sampled from its true mathematical form in equation (2.5).Additionally, the objects we are interested in go beyond point sources. A real object consistsof a large—but finite—set of positron-emitters. From the linearity of equation (2.5), the sinogramrecorded during the scan is the superposition of the individual sinograms due to each point source.Each of these sinograms are a discretized and undersampled approximation of their theoreticalfunctional form.The above explanation applies to imaging 2D slices of a true 3D object. Although it is possible toperform 3D PET imaging by recording a series of sinograms for each 2D slice, and then stacking the92.3. Image reconstructionFigure 2.2: Parameterizing a LOR connecting two detectors in the xy-plane. An annihilation event isindicated by a star.(a) (b)Figure 2.3: Demonstrating the transformation of a point source in image space [the star in (a)] to a sinusoidin sinogram space (b). Select LORs are mapped to their corresponding points in sinogram space for clarity.slices together post-imaging, this method leads to low scanner sensitivity [18]. Scanner sensitivityrefers to the proportion of positron emissions that are actually detected by the scanner. If weimage each 2D slice (i.e. in the xy-plane) separately, any gamma ray pairs that have an appreciableazimuthal component (i.e. nonzero projection onto the z-axis) will not be detected.The remedy is to define the LORs—and hence the sinogram—in R3 by adding an azimuthalvariable φ, which describes the angle of the normal vector of the LOR with respect to the z-axis.102.3. Image reconstructionThis allows for nonparallel trajectories by gamma ray pairs with respect to the slices to be detected.2.3.2 Analytic reconstructionThe method of PET imaging requires that we record data as projections. Obviously, the endgoal is to obtain the most accurate representation of the initial object. As mentioned, the detectorgeometry constrains the number of LORs we can define and thus contributes to the finite resolutionof the image to be reconstructed. The smallest volume defined in a reconstructed image is referredto as a voxel. In a real object, there is a large number of positron emitters per voxel. This allowsus to model the initial object to be imaged as continuous: f(x) = f(x, y, z). Thus the process ofdata acquisition by the scanning is the continuous-to-discrete mappingf(x, y, z) 7→ S(sp, θq, φr) (2.7)This operation is referred to as forward-projection. From here, the sinogram data need to bereconstructed via the discrete-to-discrete mappingS(sp, θq, φr) 7→ g(xi, yj , zk) (2.8)This operation is referred to as backprojection. Such a mapping would be produced by thediscretized inverse Radon transform. Ignoring the discrete nature for a moment and returning toR2, we can use the central-slice theorem to determine such an inverse.Let F (u, v) be the 2D Fourier transform of f(x, y), and Fn and F−1n represent the Fouriertransform and inverse Fourier transform operators in n-dimensions, respectively. The central-slicetheorem states that the Fourier transform of the projection of a 2D object f(x, y) is equal to the2D Fourier transform of f(x, y), sliced through the origin parallel to the projection lineF1{p(s, θ′)} = F2{f(x, y)}|θ=θ′ (2.9)The true object can be written in terms off(x, y) = F−12 {F (u, v)} =∫ ∞−∞∫ ∞−∞F (u, v) exp [2pii(ux+ vy)] dudv (2.10)Translating to polar coordinates (u = r cos θ, v = r sin θ =⇒ dudv = rd|r|dθ)f(x, y) =∫ 2pi0dθ∫ ∞0|r|F (r cos θ, r sin θ) exp [2piir(x sin θ + y sin θ)] dr (2.11)From the central-slice theorem, F (r cos θ, r sin θ) = P (r), the Fourier transform of the sinogram.We can also recognize that x sin θ + y sin θ = s, and thus can change the integration limits to beover LORs instead of radial lines: s ∈ (−∞,∞), θ ∈ [−pi/2, pi/2)f(x, y) =∫ pi/2−pi/2dθ[∫ ∞−∞|r|P (r) exp(2piirs)ds]s=x sin θ+y sin θ(2.12)The inner integral is simply the inverse Fourier transform of |r|P (r), which is subsequentlyprojected along the line x cos θ+ y sin θ− s = 0. Methodologically, equation 2.12 states that in thelimit of infinite projections p(s, θ), one can recover the true object, f(x, y) from these projectionsas follows:1. Take the 2D Fourier transform of the sinogram, P (s, θ)112.3. Image reconstruction2. Filter in Fourier space with a ramp filter (i.e. multiply by |r|)3. Take the inverse Fourier transform of this product4. Backproject these inverse Fourier transforms along each projection5. Integrate over all filtered backprojectionsThis algorithm is referred to as filtered backprojection (FBP) [30]. In the limit of infiniteprojections and noiseless data (i.e. infinite counts), FBP can produce an image that is identical tothe original object.In reality, we do not have an infinite set of projections, so the FBP operation is undersampled(especially at high frequencies) and hence produces blurred images. Additionally, real data hasnoise which is dominate at higher frequencies. The ramp filter in Fourier space accentuates highfrequencies, and consequently the noise is amplified. Other filters are often substituted in forthe ramp filter to help suppress noise, but then this deviates from the exact analytic solution ofequation (2.12) [18].2.3.3 Iterative reconstructionWhile analytic reconstructions are attractive since they are an exact solution to the inverse prob-lem and are relatively inexpensive from a computational standpoint, they ignore of a wealth ofinformation that is known about the imaging process, such as the random nature of the decayprocess, attenuation, scattering events, detector sensitivities, and the limited spatial resolution ofthe imaging system.Iterative reconstruction methods move away from an exact solution and approach a solution tothe inverse problem via an iterative algorithm. The basis of most iterative algorithms is the theoryof expectation maximization. The method of Maximum Likelihood Expectation Maximization(MLEM) first proposed for use in PET in 1982 [31] and later refined by in 1984 [32].Mathematical derivationMLEM works as follows. Let the experimentally-measured data be represented by the randomvector Y with a likelihood function g(Y,θ), θ a vector of parameters, and the “complete”, truedata by the random vector X with likelihood f(X, θ). In our discretized imaging system, theembedded discrete likelihood function of the measured data can be recovered through summing thecontributions of the likelihood function in the complete sample spaceg(Y, θ) =∑f(X, θ) (2.13)The idea is to maximize the likelihood (or objective) function of the complete data set, given theexperimentally measured data. It is instead easier to maximize the log of the likelihood function.This process is done in two steps, giving MLEM its name: first the conditional expectation of thelog likelihood given the dataE(ln f(X, θ)|Y, θn) (2.14)then maximize with respect to the estimated parameters θn to yield a new estimate θn+1. In thecontext of PET imaging, we can make use of the stochastic (Poisson) nature of positron decay (andthus event counting) to build our complete data set. Therefore if we let Xij be the number of pairs122.3. Image reconstructionof photons emitted from the jth voxel and contributing to the ith LOR, then Yi, the total numberof counts along the ith LOR isYi =∑j∈IiXij (2.15)where Ii is the subset of voxels that are intersected by the ith LOR. The Xij and Yi both formvectors whose values follow the Poisson distribution. However Yi is what we observe whereas Xij isthe (unknown) exact voxel-by-voxel description of coincidence events. In building a Poisson model,the mean positron emission rate is given by Pijλj , where λj is the activity in voxel j, and P is thesystem matrix, each element indicating the probability that an emission in the jth voxel leads to acount in the ith LOR. It is within this matrix that information about the imaging system, such asthe geometric setup, detector inefficiencies, and patient motion, is encoded.From Poisson statistics, the likelihood function for an element of the complete vector is givenasf(Xij , λj) = e−Pijλj (Pijλj)Xij/Xij ! (2.16)Therefore the log likelihood of the complete vector is given asln f(X, λ) =∑i∑j∈Ii−Pijλj +Xij ln (Pijλj)− lnXij ! (2.17)First we take the conditional expectation given the vector of LOR counts, Y , and the vector ofcurrent intensity estimates, λnE(ln f(X, λ)|Y, λn) =∑i∑j∈Ii−Pijλj +Pijλnj Yi∑j∈IiPijλnjln (Pijλj) +R (2.18)where the last term R has no λj dependence and hence will vanish once the derivative is taken.After taking the derivative with respect to λj , setting equal to zero, and solving for λj → λn+1j ,the final result of the EM algorithm produces the iterative computation ofλn+1j =λnj∑i∈JjPij∑i∈JjPijYi∑j∈IiPijλnj(2.19)where Jj is the subset of LORs that the jth voxel contributes to.Conceptual interpretationsMLEM can be thought of as a two step process: for each voxel, first we forward project thecurrent estimate of activities to yield the expectation of the number of counts along each LOR,E(Yi) =∑j∈Ii Pijλnj . The estimate is then updated, before we backproject to yield the next estimateof the activity in the jth voxel. Forward projection indicates moving from an estimate in the imagedomain to the LOR domain, and backprojection correspondingly moves in the opposite direction.Equation (2.19) can be rearranged asλn+1j = λnj∑i∈JjPij [Yi/E(Yi)]∑i∈JjPij(2.20)132.3. Image reconstructionThis representation allows us to interpret the MLEM algorithm as multiplying the currentestimate of the intensity in the jth voxel, λnj , by a weighted sum of: (measured counts along theith LOR) / (expected counts along the ith LOR), normalized by the weights, where the weights arethe elements of the system matrix, Pij .In essence, the best current estimates for each voxel in image space are forward projected intoLOR space, where the resulting LOR expectation values (based on a Poisson model given thecurrent estimate) are compared to the data. For instance, if these expectation values are largerthan the data (the E(Yi) are exclusively greater than the Yi), then accordingly the weighted sum ofequation (2.20) is less than one and hence the next best estimate will be lower. In a more generalcase, the Yi/E(Yi) will be distributed above and below one. In all cases, the algorithm rewardsmore likely voxel contributions via Pij to determine the update λnj → λn+1j .The MLEM algorithm is guaranteed to improve the estimate with each iteration, such that thealgorithm converges to a global maximum likelihood estimator at λ∞. This can be seen from theconcavity of the log likelihood function, provided there are at least as many LORs as voxels [32].In practice, the initial activity estimate, λ0, is often a uniform value for every voxel. From here,the MLEM algorithm takes around 10-100 iterations to converge to reasonable image quality [33].The convergence in terms of computational time however was troublesome, especially given thecomputational resources available in the early 1980s when the algorithm was suggested, requiringimprovements before widespread clinical implementation.2.3.4 Ordered Subsets Expectation Maximization reconstructionThe path to convergence was reduced by grouping the measured LORs in ordered subsets. Thisyields a modified EM algorithm, the Ordered Subsets Expectation Maximization (OSEM) algorithm[34], only by adding a new indexλn+1,`j =λn,`j∑i∈S`Pij∑i∈S`PijYi∑j∈IiPijλn,`j(2.21)where λn,`j is now the image intensity estimate for the jth voxel based on data from the `th subsetafter n iterations of the OSEM algorithm. Therefore one iteration of OSEM is after the λn+1,`j havebeen calculated for each subset image.Subsets are generally orthogonally partitioned to speed up convergence. The number of subsetschosen is referred to as the OSEM level. It was found that for an equal level of error in EM-and OSEM-reconstructed images, the number of iterations required for OSEM compared to EMwas roughly inversely proportional to the OSEM level [34]. For instance, the mean-squared errorachieved with EM after 32 iterations can be achieved by OSEM with 16 subsets after only 2iterations.This partioning cannot be continued ad infinitium; there is an effective limit to the number ofsubsets before the individual subset images are no longer good representations of the mean values.This limit depends of the amount of data (counts) contained within the sinogram: the less data,the more each subset deviates from the mean values.In this work, the reconstructions of all simulations and real data are performed using OSEMand a recently-developed denoised version of OSEM, HYPR-AU-OSEM [35], to be discussed inSection 2.4.3.142.3. Image reconstruction2.3.5 The system matrixWe now turn our attention to how these algorithms address the physical and technical issuesimpeding a quantitatively accurate reconstruction. This information is contained in the systemmatrix, P ∈ RJ×I , where J is the number of voxels in the FOV and I is the number of LORs. Thematrix is often factorized into several components [36].A detector sensitivity component, Pdet.sens ∈ RI×I , corrects for normalization (correcting for theimperfect efficiency of the detector pairs consisting a LOR). One way the normalization correctioncan be achieved is by performing a long scan with a rotating rod source [18]. The source emitsisotropically, and thus if all detectors had identical efficiency, after a very long time the counts alongeach LOR will converge to a common average. Therefore for differences in efficiency, a correctionfactor is assigned by dividing the average counts for all LORs by the number of counts for a specificLOR.Detector blurring effects—such as inter-crystal scattering—can be accounted for in sinogramspace through Pdet.blur ∈ RI×I , whereas blurring in image space—resulting from effects such as thepositron range—can be accounted for in Ppositron ∈ RJ×J .A geometrical mapping between image space and sinogram space is handled via Pgeom ∈ RI×J .This can be constructed by as follows. For a LOR with length Li, each voxel the LOR passesthrough contains a fraction xj/Li of the total length. Geometrically, the probability that a countalong the ith LOR came from the jth voxel is then just that fraction. If TOF is used, then theprobabilities can be narrowed to a more specific region.Finally, corrections for attenuation are contained within Pattn ∈ RI×I , where each diagonalelement corresponds to the attenuation coefficient at that voxel. Attenuation coefficients can beobtained from a blank scan and a transmission scan [18]. As the name implies, the blank scan isconducted with a rotating rod source, without the patient. A transmission scan with both the rodand the patient is then conducted. For each LOR, with the object in place the photons emittedfrom the rotating rod may be attenuated along any given LOR. Based on the reduced number ofcounts in the transmission scan compared to the blank scan for each LOR, an attenuation map(µ-map) can be constructed by backprojecting these data into image space.Beyond the system matrix, the only additional knowledge from the physics of the system thatremains is scattering and random coincidences. As detectors detect impinging photons, singlesare recorded. If another detector detects a photon within the coincidence window, then a promptis assigned to the corresponding LOR. Of all prompts, only a fraction are true coincidences, ortrues [29]. The remaining prompts are from scatter coincidences (scatters) or random coincidences(randoms).In the past, randoms were commonly estimated using a delayed coincidence method [28, 37],where a timing delay τ ′ was introduced to the timing window, such that all trues and scatters aremissed, thus all prompts recorded are randoms. Consider detectors a and b, which form the LORLab, having singles rates s˙a and s˙b, respectively. The randoms rate along Lab is thenR˙ab = 2τ′s˙as˙b (2.22)This technique has the disadvantage of requiring separate bandwidth for detecting trues andrandoms. Therefore randoms are now commonly estimated from singles [28]. If one integratesequation (2.22), assuming s˙a and s˙b evolve similarly over time, the total number of randoms overa scan duration T detected by Lab isRab = 2τ′sasbT∫0exp(−2λt)dt = ksasb (2.23)152.4. Denoisingwhere sa and sb are the total number of singles recorded in detectors a and b, respectively, through-out the scan, λ is the decay constant of the tracer being used, and k is a constant. Thus from themeasured number of singles over the scan duration in each detector, once can estimate the randomsfor each LOR.Removing scatters is a more difficult task. As mentioned, virtually all 511 keV photons scattervia the Compton effect. The effect of scatter is detrimental to image quality. Consider a 30%energy discrimination in the detectors. From the Compton energy formula, this results in a ratherlarge maximum scattering angleθmax = arccos[1− 511hν(hνhν ′− 1)]= arccos[1−(10.7− 1)]≈ 55◦ (2.24)The cross section for Compton scattering is known through the analytic Klein-Nishina equation.Therefore Monte Carlo or other simulation methods can be employed to model the number ofscatters, given the detector geometry and recorded prompts [38–40]. Simulated methods generallycannot estimate scatter from outside the FOV, but one can tail-fit the sinograms assuming aGaussian distribution to account for this effect [41].2.4 DenoisingDeviations in a PET image from the true object can be due to resolution, bias, and noise. Resolutionissues have been discussed before, and are an unavoidable consequence of the positron range,noncollinearity of positron-electron annihilation, crystal geometry effects, crystal scattering effects,and, in some cases, the finite number of voxels in an image1Since finite resolution is unavoidable, the theoretical best we can do in an imaging operationis approach a version of the true, continuous object, discretized and degraded to the resolution ofthe imaging system, which we shall call the ground truth. Any deviations from the ground truthare referred to as bias. Bias can be noise-induced, or produced during the acquisition process,reconstruction, or in post-processing when applying denoising operations.Noise in PET imaging stems from the random nature of positron emission. β+ decays obeyPoisson statistics. In a Poisson process, if repeated observations of the process over a constant timeinterval ∆T show a mean of N events, then the variance across observations will also be N . Itfollows then that the uncertainty in the number of counts is√N , corresponding to a relative errorof 1/√N .Unlike other medical imaging modalities such as MRI or computed tomography, PET is inher-ently noisy, with the major source of noise arising from the imaging itself rather than electronicnoise in the detection apparatus. The diminishing returns in Poisson relative error also make preci-sion rather difficult in PET imaging: to get to 1% precision in the reading for a particular LOR, wewould require 104 counts. With billions of possible LORs in a typical research scanner, we wouldthen require on the order of 1013 counts per scan, which would require an exorbitant amount ofradiotracer—not feasible in terms of radiation safety towards the participant, nor in an economicsense.2.4.1 FilteringThe inherent noisy and hence low precision nature of PET motivates the need for extensive denoisingprocedures. The current standard for denoising is simply to filter the image. In a filtering operation,1Often, the voxel size chosen is smaller than the intrinsic resolution of the scanner, such that the voxel size doesnot introduce further resolution degradation.162.4. Denoisingthe image data, I(x, y, z, t) are convolved with a kernel F (x, y, z, t) to produce a filtered imageIF (x, y, z, t)2IF (x, y, z, t) = I(x, y, z, t)⊗ F (x, y, z, t) (2.25)This formulation states that one can denoise dynamic (4D) images with a full spatiotemporalkernel. Depending on the application, such a process may be necessary. For instance, in kineticmodeling the precision of the time-course of the radiotracer concentration—known as the timeactivity curve (TAC)— in a particular voxel or region of interest (ROI) is paramount. In dynamicdatasets, the number of counts per frame can be low, and hence the noise is quite high in voxel-levelTACs. Therefore some temporal smoothing may be helpful.However it is more common to simply apply a spatial kernel alone to the images correspondingto each time point (hereafter referred to as frames), or a spatial kernel followed by a temporalkernel. In the latter case, the denoised image is computed in two steps: first to separately generatespatially filtered frames, ISF , which are subsequently temporally filtered to yield a spatiotemporallyfiltered image ISTFISF (x, y, z, t`) = I(x, y, z, t`)⊗ F (x, y, z) (2.26)ISTF (xi, yj , zk, t) = ISF (xi, yj , zk, t)⊗ F (t) (2.27)The functional form of the kernel is arbitrary, but is usually a Gaussian or a simple boxcar. Fora general spatial filter, one can define a kernel of dimensions (2m+ 1)× (2n+ 1)× (2p+ 1) (m, n,p all integers) such that a given voxel in the denoised image is computed asISF (xi, yj , zk) =m∑m′=−mn∑n′=−np∑p′=−pwm′n′p′I(xi+m′ , yj+n′ , zk+p′)m∑m′=−mn∑n′=−np∑p′=−pwm′n′p′(2.28)where wm′n′p′ are the kernel weights. For a simple boxcar, each weight is one. For a Gaussian, theweights arewm′n′p′ = exp[−((m′)22σ2x+(n′)22σ2y+(p′)22σ2z)](2.29)where σ represents the standard deviation in the labeled dimension, in units of voxels. The nor-malization factor in the denominator of equation (2.28) ensures that the mean value of the voxelsinvolved in each step of the convolution operation is preserved, independent of the functional formof the weights.The filtering operation is effective at denoising since in the weighted average of neighbouringvoxels, the true signal will add coherently, whereas the noisy signal will add decoherently. However,filtering fundamentally degrades resolution. Therefore the trade-off in this denoising operation isa loss of high frequency features, such as functional boundaries. Although the mean of the voxelswithin a kernel-sized neighbourhood is preserved, the higher values will be negatively biased by thelower values, and vice versa. Therefore the distribution of voxel values within the neighbourhoodwill approach uniformity.2Hereafter the spatiotemporal variables x, y, z, and t are to be assumed as discrete, with the spatial variablesindexing over the number of voxels in each dimension, and time indexing over the number of 3D images in the dynamic4D dataset172.4. DenoisingThe most obvious example is to consider imaging an axial line of positron-emitting point sources(of uniform activity). The ground truth image would only have a single line of nonzero voxels alongthe axial direction. The noisy image will have a certain variance along the line. Applying a spatialfilter will reduce this variance, but also will assign nonzero values to the voxels adjacent to the line.Since the mean in each convolution step needs to be preserved, the values along the line will benegatively biased. Combined, we can see that the filtering operation causes some activity to “spillout” from the line into the surrounding regions.This “spillover” effect is related to a greater issue referred to as the partial volume effect (PVE).In general, PVE occurs when the structure being imaged is on the order of the resolution of theimaging system—or smaller. As mentioned, the finite resolution of the imaging system meansthe true activity from positron emitters in a voxel-sized volume of the object is mapped to thecorresponding voxel in the image via a weighted average, the weights corresponding to the functionform of the point spread function (PSF) of the scanner. Typically, the PSF is modeled with aGaussian, and hence we can model the resolution degradation and corresponding PVE using thefunctional form of equations (2.28) and (2.29). In other words, the imaging process itself can bethought of as an initial filtering, before the effects of noise are considered. This way of thinking willbe readdressed in Section 3.1, when we discuss the method of generating realistic, noisy simulations.2.4.2 HYPR denoisingThe noise-resolution trade-off presented by filtering motivates the development of more sophisti-cated operators that can still lower noise-induced variance in the image, while preserving resolution.The bias induced by a filter means quantitative approaches to analyzing PET images will yield inac-curate results, albeit more precise than non-denoised data. In kinetic modeling, for instance, afterfiltering the neurophysiological parameters estimated often become negatively biased in regions ofhigh concentration, since the time-course of the activity has been combined with the lower activitymeasurements of adjacent regions. Furthermore, high activity regions usually are most importantwhen considering quantification, since the radiotracer is evidently targeting these regions.When we have a dynamic dataset, I(x, y, z, t), we can take advantage of this extra dimension byemploying HighlY constrained backPRojection (HYPR) [42], a denoising operator that uses a highersignal-to-noise ratio (SNR) composite image, C(x, y, z), to guide the denoising of the individualframes. HYPR was originally developed for use in time-resolved MRI as a reconstruction methodfrom k-space projections [43]. In such dynamic datasets, the composite image contains valuablehigh frequency features that are present in the individual frames, though they are mixed withthe high frequency contributions from noise. Thus the aim of HYPR is to borrow the denoisingabilities of filtering, but to preserve the high frequency features—using the composite—that aretraditionally removed in a filtering operation. In other words, the composite serves as a guide tothe denoising process.The HYPR operator applied to frame t` works as followsH{I(x, y, z, t`)} = C(x, y, z)I(x, y, z, t`)⊗ F (x, y, z)C(x, y, z)⊗ F (x, y, z) (2.30)Again, the filter can have an arbitrary functional form, but has been shown to be most effectivein terms of minimizing bias when a Gaussian is used [44]. The composite image is a weighted sumof the framesC(x, y, z) =N∑l=1w`I(x, y, z, t`) (2.31)182.4. DenoisingThe weights w` need not be normalized in this case, because unfiltered and filtered versions ofthe composite appear in the numerator and denominator, respectively, of equation (2.30). Notethat one may elect to redefine the weights for each frame, meaning the composite is different foreach frame. For example, a set of weights resembling a sliding window may be applied such thatonly nearby frames to the frame being denoised are highly-weighted [44].The operation can be broken down as follows: First, the raw image is filtered, which we knowwill denoise the image, but also introduce some bias. Next, the high SNR, unfiltered compositeis divided by the filtered composite, effectively extracting the high frequency features from thecomposite. This ratio is multiplied by the filtered frame to “re-inject” the high frequency featuresinto the filtered frame.The HYPR denoised image will be slightly biased towards frames with higher weights. Thereforethe sliding window technique described above sounds attractive to minimize this bias. However,only weighting nearby frames means some of the counts from the dynamic dataset are ignored, thuscreating a lower SNR composite. As a result, general practice is to include all decay-uncorrectedframes3 to maximize the denoising for a small trade-off in bias. Although this argument may soundidentical to that of the noise-bias tradeoff in the filtering operation, the biases we are consideringwith the HYPR operator are much less significant than those induced by filtering. In other words,the greedy researcher is not punished as severely for their wanting the best denoising possible.Bias can also be reduced by iterating the HYPR operator [45]. By iterating, this means updatingthe composite image to be the HYPR denoised frame of the previous iteration and reapplying theoperator to the original, non-denoised frame. However, as always this reduction of bias is paidfor by an increase in noise (since the new composite, while less biased, has lower SNR than theoriginal composite). Therefore, in practice only a single iteration of HYPR is applied to maximizethe denoising capabilities.A single frame from a dynamic RAC PET image is shown in Figure 2.4 without any denoising,as well as the filtered and HYPR post-processed versions of the same image (both using an isotropic5 mm Gaussian kernel) for a visual comparison of the denoising ability—in the spatial domain–ofeach method. To assess the effect of denoising in the temporal domain, we also present a voxel-levelTAC from the same dynamic dataset in Figure 2.5.2.4.3 Incorporating HYPR within image reconstructionThe HYPR operator of equation (2.30) can be abstracted as follows: We start with a dataset ofdimensionality N . By performing a weighted sum along one of these dimensions, we form a higherSNR composite image of dimensionality N − 1. This composite is then multiplied by a ratio ofthe filtered version of the dataset (decomposed along the same dimension as the summation tomatch the dimensionality of the composite) and a filtered version of the composite. The HYPRoperator is effective at denoising so long as the composite image has sufficiently higher SNR andthe multiplication step introduced minimal bias into the denoised image.In this abstract formulation, the composite image is not limited to a summation in the temporaldimension for a dynamic PET dataset. Wherever an extra dimensions lurks, we have the potential toemploy HYPR denoising. Referring back to the formulation of the OSEM reconstruction algorithm,in reconstructing a 3D image we add another dimension by dividing the LORs into subsets. Foreach subset, we run the traditional MLEM algorithm to arrive at a separate image estimate basedon a certain subset of the LORs. Therefore upon defining L subsets, we have L estimations of the3The decay-uncorrected frames are used to maximize the SNR of the composite, and hence the efficacy of thedenoising operation as well. This comes with the trade-off of biasing low count frames. See Section 4.2.2 for a morein-depth discussion of these biases.192.4. DenoisingFigure 2.4: Comparing a raw RAC image (a) with a filtered (b) and HYPR post-processed (c) versionof the same image. Both operations used an isotropic 5 mm Gaussian kernel. However, the loss of highfrequency features is visually evident in the filtered image as compared to the HYPR post-processed image.Colourbar units are Bq/cm3.Figure 2.5: Comparing a raw voxel-level TAC (in the putamen) from a dynamic RAC dataset with thespatially filtered and HYPR post-processed versions. Both operations used an isotropic 5 mm Gaussiankernel. However, negative bias is evident in the filtered image as compared to the HYPR post-processedimage.3D image during each iteration, forming a “4D” dataset. These images can then be summed toform a higher SNR composite image, allowing the HYPR operator to be applied to denoise eachsubset image.This denoised reconstruction is called HYPR-OSEM [35]. In order to form the initial compositeimage, we require one iteration of OSEM to be performed. Afterwards, a composite image forthe nth iteration, Cn(x, y, z), can be defined by summing over the subset images of the (n − 1)thiterationCn(x, y, z) =L∑`=1In−1,`(x, y, z) (2.32)where In−1,`(x, y, z) is the image estimate after the `th subset of the (n−1)th iteration. The HYPR202.5. Kinetic modellingoperator then becomesH`{In,`(x, y, z)} = Cn(x, y, z)F (x, y, z)⊗ In,`(x, y, z)F (x, y, z)⊗ Cn(x, y, z) (2.33)where F (x, y, z) is a spatial Gaussian kernel. We have a few options as to where we wedge theHYPR operator into the OSEM formulation of equation (2.21), but it has been demonstrated fromsimulations and human data that the most accurate and precise reconstruction is produced byapplying the HYPR operator after each subset update (HYPR-AU-OSEM) [35]λn+1,`j = H` λn,`j∑i∈S`Pij∑i∈S`PijYi∑j∈IiPijλn,`j (2.34)Here the HYPR operator is applied to the output image of the (n + 1)th iteration of thereconstruction. This formulation is used as the sole 3D denoised reconstruction method in thiswork.In applying HYPR-AU-OSEM to dynamic data, one reconstructs each frame separately using acomposite defined only using information from that frame. Further denoising can then be achievedby applying the HYPR operator of equation (2.30) as a post-processing step, where the compositeis defined using information from all denoised frames. This is a pseudo-4D approach to denoising,in the sense that the temporal information is not used to guide the denoising within the reconstruc-tion, but nonetheless is exploited in post-processing after the spatially-denoised reconstruction iscomplete.2.5 Kinetic modellingAfter the reconstruction and denoising of the raw PET signal, the next step for the analysis ofdynamic PET images is most often performed by kinetic modeling. Here, neurophysiological modelsbased on the theory of compartmental modeling are fitted to the denoised TACs. The TACs maybe defined on a single voxel, or a grouping of voxels either comprising a subregion of an anatomicalstructure or the entire structure itself.Anatomical structures are generally not defined by contrast in PET images—instead contrastpoints to functional patterns in the brain. Therefore an MRI scan (which shows anatomy throughcontrast) is often taken of the participant being scanned, registered in PET space, and then seg-mented into anatomical regions.To understand kinetic modeling, we move away from the physics of detecting where a positronemission occurred, and towards the biophysics of why a radiotracer molecule ended up in thatposition in the first place.The brain is an immensely complicated organ, requiring vast simplification to model its intricatenetworks. Nonetheless, we can start with the thought experiment of traveling alongside a singleradiotracer molecule, and speculating where it may travel. Radiotracers are chosen to target aspecific system of the brain, often through binding to receptors within that system. Since this workfocuses exclusively on the dopamineric system, we shall use this system as an example.If we inject a molecule targeting the dopaminergic system into the bloodstream, it will be carriedby the plasma through the circulatory system. From here, the molecule may be exchanged fromthe plasma via capillaries across the blood-brain barrier, into brain tissues.212.5. Kinetic modellingSome radiotracers specifically target function within a region of the neurons known as a synapse[46]. A synapse is the junction between axon terminals of a neuron and another neuron. When anelectromagnetic signal within a neuron reaches an axon terminal, neurotransmitters are required toallow the signal to continue to the post-synaptic neuron. One such neurotransmitter is dopamine(DA), common in neurons projecting to a structure of the brain known as the striatum. Ourradiotracer molecule may be targeting presynaptic structures/processes, i.e. targeting DA synthesisand transport (6-[18F]FDOPA, 6-[18F]FMT, [11C]DTBZ), or postsynaptic, i.e. targeting dopaminereceptors ([11C]raclopride, [18F]fallypride) [47, 48]. These radiotracer molecules either bind to theproteins involved in presynaptic function or to the DA receptors. We will then have more bindingin regions with higher synthesis of DA, higher transport of DA, or higher density of DA receptors.With more molecules binding in these regions, the concentration of positron emissions (in otherwords, activity) will also be higher.2.5.1 Compartmental modelling: plasma input modelsWe can model these tissue regions where binding occurs as biological compartments, holding acertain concentration of radiotracer. In general, we can have an arbitrary number of compartments,exchanging both with the plasma and with each other.In a general PET system made up of N compartments, the concentration of radiotracer inthe plasma, Cp(t), is exchanged with the total tissue concentration, CT (t) with rate constant K1[49]. Historically, tissue concentration CT (t) is measured in Bq per cm3 of tissue, whereas plasmaconcentration Cp(t) is in Bq per mL of blood [46]. Thus K1 has units of mL blood per cm3 tissueper unit time.In addition to radiotracer in tissue compartments, a PET measurement will also include ra-dioactivity coming from tracer present in the blood, i.e. the blood radioactivity concentration,CB(t) [50]. If we define the vector of tissue compartment concentrations CT(t)CT(t) =[CT1(t) CT2(t) · · · CTN (t)](2.35)then we can describe a general compartmental system of differential equations in matrix formC˙T(t) = ACT(t) +[K1e1 0N] [CP (t)CB(t)](2.36)where A is a matrix whose elements are combinations of rate constants (often called the statetransition matrix) and e1 =[1 0 · · · 0]T is a standard length-N unit vector and 0N is length-N vectors of zeros. “T ” here denotes the transpose operator. If VB is the blood volume fractionof radiotracer in the system, then the total tissue concentration is a weighted sum of the tissuecompartment concentrations and the whole blood concentrationCT (t) = (1− VB)1NTCT(t) +[0 VB] [CP (t)CB(t)](2.37)Finally, since the scan begins upon injection of the radiotracer into the body, the compartmentsinitially are devoid of radiotracerCT(0) = 0 (2.38)222.5. Kinetic modellingUnder these initial conditions, it can be shown that the solution to such a system is [50]CT (t) = (1− VB)HTP (t)⊗ CP (t) + VBCB(t) (2.39)where HTP (t) is the impulse responseHTP (t) =N∑i=1φie−θit (2.40)Conceptually, this solution can be understood by first analyzing a simpler system. Consider asystem with a single compartment exchanging radiotracer with plasma. If C1(t) is the concentrationof radiotracer in the target and k2 is the rate constant for the exchange of radiotracer from thecompartment back into the plasma, then this system is described by the following differentialequationdC1(t)dt= K1Cp(t)− k2C1(t) (2.41)where k2 necessarily has units of inverse time. Now consider a single bolus injection of magnitudeB into the tissue compartment at t = 0. We would model this input function as a delta functioncentred at t = 0 [i.e. Cp(t) = Bδ(t)]. If the radiotracer does not recirculate, then the solution tothis simple system isC1(t) = BK1e−k2t (2.42)This is the tissue response to the bolus input function. If we instead present the bolus B1 at atime t1 > 0, then the tissue response becomesC1(t) = θ(t− t1)B1K1e−k2(t−t1) (2.43)where θ(t) is the Heaviside function. Since the differential equation is linear, if we introduce asecond bolus B2 at a time t2 > t1, then the tissue response is the sum of the individual responses,by the Principle of SuperpositionC1(t) = θ(t− t1)B1K1e−k2(t−t1) + θ(t− t2)B2K1e−k2(t−t2) (2.44)So for M bolus injections of magnitudes {B1, B2, · · · , BM} at times {t1, t2, · · · , tM}, the tissueresponse isC1(t) =M∑i=1θ(t− ti)BiK1e−k2(t−ti) (2.45)In the limit M →∞, the input delta functions form a continuous input function, and the sumabove becomes an integral. Since the radiotracer is delivered to the tissue compartment throughthe plasma, this input function is simply Cp(t), which can be measured from arterial blood samples.The tissue response to the plasma input function is thenC1(t) =t∫0Cp(τ)K1e−k2(t−τ)dτ (2.46)232.5. Kinetic modellingThis integral is nothing but the convolution of the plasma input function with the tissue reponseto a single bolus input (i.e. K1e−k2t)C1(t) = Cp(t)⊗K1e−k2t (2.47)This is the solution for a single compartment model. If we generalize to N compartments,insofar as the exchange of radiotracer between the compartments and plasma is linearly related tothe concentration in the compartments, the same principle above applies. That is, if the responseof the N -compartment system to a bolus input function is HTP (t), then the tissue response fromthe compartments alone, CNT (t), is given byCNT (t) = Cp(t)⊗HTP (2.48)Including the contribution from whole blood, we arrive at the solution for the total tissueresponse, CT (t), as seen in equation (2.39). If the exchange between the compartments is alsolinear, and the system is reversible, then the impulse response is a summation of exponentialswhose decay constants, θi, are combinations of the rate constants between the compartmentsHTP (t) =N∑i=1φie−θit (2.49)Since the θi are combinations of rate constants (units of inverse time), and φi has units of K1(mL per cm3 per unit time), VD has units of mL/cm3.2.5.2 Compartmental modeling: reference tissue modelsObtaining arterial plasma samples, required to determine the plasma input function, is an invasiveprocedure for the participant being scanned, and modeling the input plasma function introducesadditional assumptions into the model. Thus it is desirable to determine an input function purelyfrom data in the image itself. If there exists a large enough region of the brain where specificbinding is negligible—i.e. a reference region—then a reference tissue model can be used instead ofa plasma input model [46].As before, the model assumes radiotracer exchanges with tissue in a target region, CT (t) (i.e.outside the reference region), containing N compartments, with a rate constant K1. Radiotraceris also exchanged between the plasma and a reference region with M compartments with a rateconstant K ′1.If V ′B is the blood volume distribution in the reference region, then the total concentration inthe target tissue and reference regions are given by [50][CT (t)CR(t)]=[(1− VB)1T 0T0T (1− V ′B)1T] [CT(t)CR(t)]+[0 VB0 V ′B] [Cp(t)CB(t)](2.50)where 0 and 1 are matrices of zeros and ones, respectively. With the reference region time course,CR(t), taking the place of the plasma input function, our solution is now of the formCT (t) = CR(t)⊗HTR(t) (2.51)242.5. Kinetic modellingIf the whole blood concentration is sufficiently low to be ignored (i.e. VB ≈ 0, V ′B ≈ 0), then[50]HTR(t) = φ0δ(t) +M+N−1∑i=1φie−θit (2.52)As t→∞, HTR → φ0δ(t). This yieldsφ0 = limt→∞CT (t)CR(t)=K1K ′1≡ R1 (2.53)Thus general solution for a reference tissue model with N target compartments and M targetvolumes with negligible whole blood contribution isCT (t) = R1CR(t) + CR(t)⊗(M+N−1∑i=1φie−θit)(2.54)2.5.3 Simplified Reference Tissue ModelA specific example of a reference tissue model that is used extensively in this work is the SimplifiedReference Tissue Model (SRTM) [51]. SRTM assumes a single compartment in both the targetregion and reference region (N = M = 1).The “Simplified” in SRTM comes from reducing a two tissue compartment model—theoreticallyappropriate for target regions containing free (Cf ) and specifically-bound (Cb) tracer—to a singlecompartment. This is often a reasonable approximation, since the two compartments are not alwayskinetically distinguishable in the PET data.Radiotracer is transferred from the plasma to the free compartment at a rate K1, and from thefree compartment back to the plasma at a rate k2. Similarly, radiotracer is exchanged from theplasma to the reference region at a rate K ′1, and in the opposite direction at a rate k′2.As well, radiotracer can also be exchanged between the free compartment and the specifically-bound compartment at a rate k3, and in the opposite direction at a rate k4.SRTM combines these two compartments by making the assumption that the distribution vol-ume of non-specific binding tracer is the same in the target and reference regions, such thatK ′1k′2=K1k2(2.55)Recognizing K1K′1≡ R1 =⇒ k′2 = k2R1 . Once the two target compartments are combined to one,the transfer of radiotracer from the plasma to the compartment is still described by K1, but nowthe exchange from the combined compartment back to the plasma is replaced by an effective rateconstant, k2a. This process of combining the two target compartments to the single compartmentof SRTM is illustrated in Figure 2.6.With the relevant constants defined, from equation (2.54) with N = M = 1, SRTM then takeson the general formCT (t) = R1CR(t) + φ1CR(t)⊗ e−k2at (2.56)252.5. Kinetic modelling(a) (b)Figure 2.6: Visualizing tracer exchange in SRTM. The model combines a free compartment and specifically-bound compartment in (a) to a single compartment in (b).From a fitting perspective, this model has three free parameters that can be fit with nonlin-ear regression methods, or through linear least squares if a set of NBF basis functions, B, withreasonable estimates of k2a are pre-specifiedB = {CR(t)⊗ e−k(1)2a , CR(t)⊗ e−k(2)2a , · · · , CR(t)⊗ e−k(NBF )2a } (2.57)From a physiological perspective, R1 describes the relative rates of radiotracer delivery from theplasma to the target region compared to from the plasma to the reference region. As mentioned,k2a specifies the rate of radiotracer exchange from the target region to the plasma. φ1 has noimmediate physiological meaning, but can be shown to be a combination of rate constants that areinformative. To do so, we can go back to first principles and write down the differential equationsdescribing this single compartment systemdCT (t)dt= K1Cp(t)− k2aCT (t)dCR(t)dt= K ′1Cp(t)− k′2CR(t) (2.58)Eliminating Cp(t) via substitution yieldsdCT (t)dt=K1K ′1(dCR(t)dt+ k′2CR(t))− k2aCT (t) (2.59)Replacing k′2 withk2R1, this yieldsCT (t) = R1CR(t) + k2t∫0CR(u)du− k2at∫0CT (u)du (2.60)262.5. Kinetic modellingSubstituting equation (2.56) for CT (t) in the integrandCT (t) = R1CR(t) + k2t∫0CR(u)du− k2at∫0R1CR(u) + φ1CR(u)⊗ e−k2audu= R1CR(t) + (k2 −R1k2a)t∫0CR(u)du− k2aφ1t∫0CR(u)⊗ e−k2audu (2.61)Now we will consider the limiting condition of t→∞. Taking the limit allows us to employ thefollowing property of the convolution∞∫−∞f(x)⊗ g(x)dx =∞∫−∞f(x)dx∞∫−∞g(x)dx (2.62)Since CR(t) = 0 ∀t < 0, we havelimt→∞CT (t) = R1 limt→∞CR(t) + (k2 −R1k2a)∞∫0CR(t)dt− k2aφ1∞∫0CR(t)⊗ e−k2atdt= R1 limt→∞CR(t) + (k2 −R1k2a)∞∫0CR(t)dt− k2aφ1∞∫0CR(t)dt∞∫0e−k2atdtlimt→∞CT (t) = R1 limt→∞CR(t) + (k2 −R1k2a − φ1)∞∫0CR(t)dt (2.63)However, from equation (2.53), or from taking t→∞ in equation (2.56) we also know thatlimt→∞CT (t) = R1 limt→∞CR(t) (2.64)which implies k2 −R1k2a − φ1 = 0 or∞∫0CR(t)dt = 0. The latter is physiologically impossible, so itmust be thatφ1 = k2 −R1k2a (2.65)Therefore, by fitting R1, φ1 , and k2a in equation (2.56), one can derive the final parameterof interest for the system, k2, using the above result. Alternatively, in this derivation we haverecovered an integral equation for SRTM in equation (2.60) whose coefficients are R1, k2, and k2a.If one inputs the measured CR(t) and CT (t) in the respective integrands, then linear least squaresfitting can be used to directly fit the parameters of interest.Finally, one can also define the tracer non-displaceable binding potential [49]BPND =k2k2a− 1 (2.66)272.5. Kinetic modellingBPND is a PET estimate for the receptor binding of the radiotracer (BPND ∼ Bmax/Kd whereBmax is the maximum free binding site density and Kd is the ligand-target affinity), and henceserves as a proxy for free receptor density.If equation (2.56) is used in fitting, then BPND can be derived fromBPND =φ1 −R1k2ak2a− 1 = R1 − φ1k2a− 1 (2.67)Otherwise if equation (2.60) is used in fitting, then BPND can be calculated using the originaldefinition.While SRTM has proven to be a robust model in most cases [51–53], the assumption that thesystem can be approximated by a single compartment is not valid for all radiotracers. Nevertheless,this assumption for the radiotracer exclusively used in this work, RAC, is valid.2.5.4 Linear parametric neurotransmitter PET modelRAC is a D2 receptor antagonist, and therefore competes for binding with DA. Thus DA releasein the target region manifests itself as a drop in CT (t), as compared to CT (t) in a baseline scenariowhere DA remains at basal levels.Historically, quantifying this drop has been achieved using two scans in a double RAC protocol[54–56]. One scan is performed in the baseline condition, where the subject is not presented withany stimuli to induce DA release. Upon fitting with SRTM, the baseline binding potential BPbaselinecan be calculated at a regional or voxel-level. A second scan is performed with the goal of achievinga DA level some constant value above basal levels. For instance, a subject can be administereda dose of L-DOPA, a DA precursor, and scanned one hour later, when the dopamineric system isassumed to be in equilibrium [55]. SRTM is fit once again, yielding BPstimulus, presumably lessthan BPbaseline since the increased DA concentration will compete with RAC for binding.From these two binding potentials, a PET-estimated DA release can be estimated [49]DAR = 100%×(1− BPstimulusBPbaseline)(2.68)The assumption of equilibrium in the dopamineric system may not be valid, and requiring twoscans is inefficient from logistical and economic standpoints, as well as potentially limiting whenwe consider radiation safety of the participant and confounding effects of scan-to-scan variability.These issues motivate a single scan protocol, where the transient effects of inducing a DA releasecan be quantified instead of relying on the assumption of equilibrium. For a transient DA releasethat evolves proportionally to some function h(t), if we assume the release removes radiotracer ata rate linearly related to h(t)CT (t), we can modify the differential form of SRTM accordinglydCT (t)dt= K1Cp(t)− k2aCT (t)− γh(t)CT (t) (2.69)Since we assume there exist no D2 receptors in the reference region and hence no DA release,we can eliminate Cp(t) using the differential equation for CR(t) as before, yieldingdCT (t)dt= R1dCR(t)dt+ k2CR(t)− k2aCT (t)− γh(t)CT (t) (2.70)282.5. Kinetic modellingIntegrating yields an integral equation for CT (t) that can be fitted using linear least squares,provided a set of basis functions for h(t) are pre-specifiedCT (t) = R1CR(t) + k2t∫0CR(u)du− k2at∫0CT (u)du− γt∫0h(u)CT (u)du (2.71)This is the linear parametric neurotransmitter PET model (lp-ntPET) [17]. h(t) is a unitlessfunction, so γ is a rate constant with units of inverse time. Intuitively, we expect the time course ofDA release to have an initial rise to some peak value, then a decay back to basal values. A gammavariate function can describe this generic shape with great versatility [57] using three parameters:tD (the time at which the DA release begins to rise above basal levels), tP (the time of peak release),and α (the “sharpness” of the response: high α describes an impulse-like release, low describes aspread-out release)h(t) = θ(t− tD)(t− tDtP − tD)αexp[−α(t− tDtP − tD)](2.72)Several possible gamma variate functions are presented in Figure 2.7 to illustrate the variety ofrelease patterns lp-ntPET can model. Corresponding simulated TACs with these release patternsare shown in Figure 2.8(a-c). The TACs for a fixed release pattern but different release magnitudes(i.e. by modulating γ) are shown in Figure 2.8(d).With the measured amount of RAC binding decreasing with increasing DA release, we candefine a time-varying binding potential as [58]BP (t) =k2k2a + γh(t)− 1 (2.73)The corresponding PET-estimated DA release level is then also time-varyingDAR(t) = 100%×(1− BP(t)BP (0))(2.74)These metrics, in theory, allow one to quantify the temporal pattern of DA release at anylocation in the brain. DAR(t) corresponding to the noiseless TACs of Figure 2.8(d) are shown inFigure 2.9.Since a vast literature already exists quantifying static DAR using the double scan protocol,some cross-comparison of the metrics would be useful. As well, to design our simulations we wantto ensure our choices of γ and h(t) in equation (2.71) are physiologically realistic. However, thiswould require reducing the time-varying quantity DAR(t) to a single representative value. Simpleaveraging will not suffice, since during the uptake portion of the scan, the patient is still in baseline,meaning BPND(t < tD) = BPND(0). Averaging only during the release period also has its flaws,since this depends on how long after the intervention the scan lasts.Instead of reducing DAR(t) directly to a static quantity, in this work we suggest a new methodthat converts a mid-scan intervention time course to a hypothetical set of equivalent double RACscans.292.5. Kinetic modelling(a) (b)(c)Figure 2.7: Modulating parameters in the gamma variate function, h(t). In (a), tP = 60 min and α = 1are fixed, while tD is modulated. In (b), tD = 35 min and α = 1 are fixed, while tP is modulated. In (c),tD = 35 min and tP = 40 min are fixed, while α is modulated.302.5. Kinetic modelling(a) (b)(c) (d)Figure 2.8: (a-c) Simulated noiseless lp-ntPET TACs based on a RAC dataset from a healthy control, withh(t) chosen to match the gamma variate release functions of Figure 2.7. The other parameters are R1 = 1.06,k2 = 0.005 s−1, k2a = 9.04× 10−4 s−1, and γ = 5× 10−4 s−1. In (d), the release strength is modulated bychanging γ for a fixed release pattern of tD = 35 min, tP = 40 min, and α = 2.312.5. Kinetic modellingFigure 2.9: Time-varying PET-estimated % DA release, DAR(t), corresponding to the TACs of Fig-ure 2.8(d). R1 = 1.06, k2 = 0.005 s−1, and k2a = 9.04 × 10−4 s−1 are fixed for all TACs used to deriveDAR(t).Consider a voxel with a baseline TAC C0(t). Adding a DA release-induced perturbation intro-duced by the choice of γ and h(t) modifies the TAC to CDA(t; γ, h). The loss in area under thecurve (∆AUC) after the intervention time tD can be computed for a given choice of γ and h(t)∆AUCDA(γ, h) =tend∫tDC0(t)dt−tend∫tDCDA(t; γ, h)dt (2.75)To resort to the single metric of DAR used in double RAC protocols, we imagine a hypotheticalthird TAC C∆BP (t) where a constant BPND decrease from that in C0(t) also results in AUC lossin the post-intervention period∆AUC∆BP =tend∫tDC0(t)dt−tend∫tDC∆BP (t)dt (2.76)By determining the ∆BP that matches these two ∆AUCs, we force the cumulative loss in RACsignal in the post-intervention period to be equivalent, and thus time-varying DAR(t; γ, h) can becondensed to a single metric DAR(γ, h)∆AUCDA(γ, h) = ∆AUC∆BPtend∫tDC0(t)dt−tend∫tDCDA(t; γ, h)dt =tend∫tDC0(t)dt−tend∫tDC∆BP (t)dttend∫tDCDA(t; γ, h)dt =tend∫tDC∆BP (t)dt (2.77)322.5. Kinetic modellingFigure 2.10: The steps involved in reducing the time-varying metric DAR(t) to the static metric DAR.The baseline TAC, C0(t), was simulated with SRTM using R1 = 1.06, k2 = 0.005 s−1, and k2a = 9.04×10−4s−1, corresponding to BP0 = 4.58. The TAC with a time-varying DA release incorporated, CDA(t), wassimulated with lp-ntPET using the same, R1, k2, and k2a, with tD = 35 min, tP = 40 min, α = 2, andγ = 4.1× 10−4. The baseline TAC is modified to generate C∆BP (t) with identical parameters as C0(t), savefor k2a = 9.9 × 10−4 (corresponding to ∆BP = 0.46), at which point ∆AUCDA = ∆AUC∆BP . Finally,DAR is calculated using equation (2.78) using the ∆BP and BP0 above to yield DAR = 10%.Figure 2.11: Lookup curves for converting between γ and DAR for the four subregions of the striatum,based on the simulations generated in this work.Baseline binding potential BP0 is decreased by increasing the k2a from C0(t) until equa-tion (2.77) holds. This allows us to compute ∆BP (γ, h), which in turn defines DAR(γ, h)DAR(γ, h) = 100%× ∆BP (γ, h)BP0(2.78)This method is illustrated graphically in Figure 2.10 for DAR = 10%. As will be elaboratedon in Section 3.1, we construct our simulations based on the kinetic parameters extracted from a332.6. Detecting DA releasereal human baseline RAC dataset. Therefore the noiseless simulated TACs in the subregions of thestriatum (L/R caudate and L/R putamen) will be slightly different. For each TAC, we can add anidentical DA-induced perturbation with lp-ntPET (tD = 35 min, tP = 40 min, α = 2) using a rangeof realistic γ values. After performing the conversion to DAR detailed above, we then generatelookup curves for each subregion. These are presented in Figure 2.11, and show little subregionalvariation, demonstrating the robustness of the conversion within high binding regions.2.6 Detecting DA releaseThe previous section has shown that lp-ntPET is able to fit to TACs where DA release is occurring,whereas SRTM cannot. This works by adding a transient term to the integral form of SRTM inequation (2.60). This additional term adds four new parameters (γ, tD, tP , and α) to the originalthree parameters of SRTM (R1, k2, and k2a).If we are given a TAC and wish to make a binary classification of release versus non-release,then a logical method is to see if the additional transient term in lp-ntPET produces a significantlybetter fit [17]. Statistically, this can be checked using an F-test. The F-statistic compares theweighted residual sum of squares (WRSS) between the two fits, as well as considering the numberof free parameters in each modelFν1,ν2 =(WRSSSRTM−WRSSlp−ntPETPlp−ntPET−PSRTM)(WRSSlp−ntPETNf−Plp−ntPET) = ν2ν1(WRSSSRTMWRSSlp−ntPET− 1)(2.79)where ν1 = Plp−ntPET−PSRTM = 4 and ν2 = Nf −Plp−ntPET = Nf −7 are degress of freedom, Nf isthe number of data points in the TAC (i.e. number of frames) and P is the number of parameters inthe model. We start with the null hypothesis that lp-ntPET fits no better than SRTM. The abovestatistic follows an F-distribution where the cumulative distribution function (CDF) indicates theprobability that the null hypothesis is not true (conversely, p = 1−CDF indicates the probabilitythat the null hypothesis is true). Therefore we can reject the null hypothesis at the standard 95%confidence threshold ifp(Fν1,ν2) = 1− CDF(Fν1,ν2) < 0.05 (2.80)This defines some critical F-statistic value Fc where p(Fc) = 0.05. For a first pass, we can fit bothSRTM and lp-ntPET to each voxel-level TAC and compute the F-statistic Fν1,ν2 . If Fν1,ν2 > Fc, weregister this voxel as one containing DA release. Otherwise, it is classified as a non-release voxel.Repeating for all voxels yields a binary significance map, where voxels classified as release are givena value of 1 and non-release voxels a value of 0.In a noisy, voxel-level baseline TAC, we can expect that the four additional parameters oflp-ntPET may try allow overfitting to noise, since there is no DA release contribution to the signal.Therefore in these cases WRSSlp−ntPET < WRSSSRTM even though there is no release in that voxel.However, we can see from the rightmost formulation of F in equation (2.79) that the four additionalparameters in lp-ntPET are penalized to account for this scenario.Nevertheless, this statistical thresholding is evidently not perfect by design; if all baseline voxelsfollow a true F-distribution, then by thresholding at 95% confidence we still have 1 in 20 baselinevoxels registering as a release voxel (i.e. a false positive). Furthermore, reconstruction algorithmsand post-filtering introduce spatial correlations between voxels. Therefore it is not uncommon to342.6. Detecting DA releaseFigure 2.12: Histogrammed cluster sizes found in the striatum of 50 simulations with no true DA release.As a first pass, false positives can be controlled by thresholding the CDF at 99%. However, the long tail ofthe distributions means the cluster size threshold (CST) that ensures 90% of simulations are completely freeof false positives is higher.find groups of false positive voxels (i.e. false positive clusters) in the significance map. Figure 2.12plots the distribution (blue histogram) and cumulative distribution (red line) of false positive clus-ters sizes from 50 baseline simulations (i.e. no DA release in any of the voxels). Although thevast majority of false positive clusters are a few voxels or smaller in size, there is a long tail in thedistribution with false positive clusters up to hundreds of voxels in size.This phenomenon motivates a further thresholding procedure [59]. From several baseline simula-tions, we can compute a cumulative histogram (i.e. CDF) of all false positive cluster sizes detected.This is shown as the red line of Figure 2.12; for each false positive cluster size on the x-axis, thecorresponding value of the CDF on the y-axis refers to the percentage of all false positive clustersthat are that size or smaller. For instance, the vertical dotted line of Figure 2.12 at a false positivecluster size of 95 voxels intersects the CDF at 0.99, meaning 99% of all false positive clusters are95 voxels in size or smaller.We are now prepared to define a cluster size threshold (CST) that removes a certain percentageof the false positive clusters in the simulations. All clusters smaller in size than this CST arerejected from further consideration as true release clusters. For instance, from our example above,if we choose our CST to be 95 voxels, then statistically we remove 99% of all false positive clusters.However, such a CST does not leave 99% of the simulations completely free of false positiveclusters. For instance, in 100 baseline simulations we may have 2000 false positive clusters thatexist. With a CST of 95 voxels, 99% of these 2000 false positive clusters are removed, leaving 20false positive clusters to be distributed amongst the 100 simulations. As a result, in order to have99% of the simulations completely free of false positive clusters, only one simulation would have tocontain all 20 remaining false positive clusters—an extremely unlikely event.What we want instead is a CST that leaves a certain percentage of simulations completely freeof false positives. In practice, [57, 60] elected to increase the CST until 90% of noisy baselinesimulations were devoid of false positives, which we shall adopt in this work as well. This CST isdenoted by the solid vertical line in Figure 2.12. Note that this value is higher than the thresholdof 95 voxels that removes 99% of all false positive clusters, meaning the CST removes more than99% of all false positive clusters, but only renders 90% of all simulations completely free of falsepositives. Another way to conceptualize this threshold is on a scan-by-scan basis: for any givenbaseline scan, there is a 10% chance that at least one false positive will exist.352.7. A novel Monte Carlo method for detecting DA releaseThis F-test based approach is the current best method for detecting transient DA release, andshall be referred to as the “standard method” for the remainder of this work.2.7 A novel Monte Carlo method for detecting DA release2.7.1 MotivationThe above statistical thresholding methods rely on the WRSS of the SRTM and lp-ntPET fitsfollowing an F-distribution, where the null hypothesis is that both fits represent the data equallywell. If this is indeed true, then in thresholding at 95% confidence we expect 5% of the voxels in abaseline scan to be identified as false positives.In a mid-scan intervention where DA release is present in clusters (i.e. not all voxels have DArelease), then the 5% false positive rate still holds only in non-release voxels. In voxels with DArelease, due to the high noise in the TACs the distribution of F values is not well separated fromthe distribution of F values for non-release voxels [c.f. Figure 2.13(a)]. As a result, enforcing a 95%confidence threshold means that we will also have a number of false negative voxels (i.e. voxels withtrue DA release that have F < Fcrit). Furthermore, some true release voxels may pass the F-testthreshold, but may be members of a cluster that is too small to pass the subsequent CST. Theserejected voxels—either by the initial F-test threshold or subsequent CST—are classified as non-release, or baseline voxels, which causes a deviation from a true F-distribution [c.f. Figure 2.13(b)].The idea of the Monte Carlo (MC) method is to return the distribution of F values for non-release voxels back to the expected theoretical F-distribution through regularization. In essence,some voxels having F < Fcrit (i.e. classified as baseline voxels by the standard F-test basedmethod) are re-classified as release voxels, and some voxels having F > Fcrit (i.e. classified asrelease voxels by the standard F-test method) are reclassified as baseline voxels, so as to maintainthe F-distribution for non-release voxels.2.7.2 FormulationIf X(F ) is a true F-distribution and X̂(F ) is the measured distribution of F values for all baselinevoxels (as classified by the detection algorithm), then such a regularization term could look asfollowsR = 1FˆmaxFˆmax∫0(X̂(F )−X(F ))2X(F )dF (2.81)where Fˆmax is the maximum measured F-value in the set of voxels being investigated. Upondiscretizing the integral for numerical computation, this term is essentially the χ2 statistic incomparing the expected F-distribution and distribution of F values for non-release voxels.The question then becomes how do we decide which F < Fcrit voxels are accepted as releasevoxels, and conversely which F > Fcrit are rejected as baseline voxels? There exists a large numberof sets of possible release ↔ baseline reclassifications that result in the distribution of non-releaseF values well-approximating an F-distribution. However, we want the solution that maximizesthe sensitivity and minimizes the number of false positives. In other words, we want to reclassifybaseline voxels for release voxels where true release actually exists, and reclassify release voxels forbaseline voxels where no release exists.Since we expect DA release to exist in clusters in the brain, we can design a cost function Ethat rewards local parity. Considering the local region Ni as the set of 6 nearest neighbours in 3D362.7. A novel Monte Carlo method for detecting DA release(a) (b)Figure 2.13: In (a), the poor separation between the F-distributions of ground truth baseline voxelsand ground truth release voxels is highlighted—contributing to a high number of false negatives. In (b),the ground truth baseline voxels correctly follow an F-distribution, but the rejected voxels (which, if thedetection algorithm was perfectly accurate would be the same set of ground truth baseline voxels) show anoticeable deviation from a true F-distribution, as a result of the high number of false negatives.with shared faces4 for the ith voxel, its parity can be described byPi = −∑j∈Nisisj (2.82)where si is the classified state of the ith voxelsi ={+1 release voxel−1 non-release voxel (2.83)Inspired by Monte Carlo approaches to the Ising spin model [61, 62], we then design the followingcost functionE = −∑iPi + βFˆmaxFˆmax∫0(X̂(F )−X(F ))2X(F )dF = −P + βR (2.84)where β is a hyperparameter controlling the degree of regularization.Since we wish to minimize E, the negative sign in front of the parity term rewards local ho-mogeneity, whereas the regularizing term ensures that the distribution of F values in non-releasevoxels does not deviate too far from a true F-distribution.To minimize E, we modify the Monte Carlo Metropolis algorithm commonly used in Ising spinmodels [62, 63]. To start, we initialize the state image (i.e. each voxel has s = ±1 representing itscurrent state) as a random pattern of −1 and +15. From here, the algorithm proceeds as follows:4Larger neighbourhoods were investigated, but did not offer any advantages, so we do not present these results inthis work.5The effect of different random initializations in considered in Section 2.7.3372.7. A novel Monte Carlo method for detecting DA release1. Choose a voxel at random and set si → −si.2. Calculate the change this flip induces on the cost function∆E = −2Pi + β∆R (2.85)3. If ∆E < 0, then accept the change. Otherwise:(a) If changing from release → baseline, accept with probability PDF (Fi)(b) If changing from baseline → release, accept with probability 1− PDF (Fi)where PDF (Fi) is the value of the probability density function for the F-distribution for theF value of the ith voxel, Fi. In other words, this is the probability that the null hypothesis istrue given Fi.4. Repeat steps 1—3 for all remaining voxels.Conceptually, one could think of this algorithm as follows. For each ∆E < 0 calculated, either:• Both the parity change and ∆R were favourable (i.e. higher parity, lower deviation of X̂(F )from X(F ))• The change in parity was unfavourable (the local neighbourhood in the state image is lesshomogeneous), but this was offset by a favourable change in R. This could occur near theedge of a true cluster, where homogeneity is no longer preferred.• The change inR was unfavourable (the change made the global X̂(F ) deviate fromX(F )), butthis was offset by a favourable parity change. Such cases will occur to avoid over-regularizingthe deviation of X̂(F ) from X(F ).When ∆E > 0, this often indicates unfavourable changes, but we occasionally allow suchchanges with a frequency related to the F value at the voxel as another means to approaching atrue F distribution for baseline voxels.One can liken the MC method to applying a spatial filter to the binary significance mask,though the kernel would be unique for each voxel, taking into account both the local parity and theglobal adherence to a correct F-distribution for all voxels; simply filtering the binary significancemask would isotropically increase the size of both true clusters and false positive clusters, whereasthe MC method will only probabilistically increase cluster size if the global deviation from a trueF-distribution suggests a subset of false negative voxels exists in that region.This algorithm can be repeated for many iterations, with the initial state image being theresult of the previous iteration. In MC approaches applied to Ising models of ferromagnets, macro-scopically after several iterations the global state (in this case the energy of the system) reachesequilibrium. However, microscopically small local fluctuations in spins states still occur in subse-quent iterations, simulating thermal equilibrium. Similarly, in our MC method we expect the globalstate of our system, as reflected by the value of E, to reach equilibrium after a certain numberof iterations. In subsequent iterations, the algorithm will shuffle through all realistic global statesthrough small local voxel classification changes.As a result, unlike other minimization problems we do not expect the progression of E tobe strictly monotonically decreasing, but rather a decrease to within an equilibrium value, about382.7. A novel Monte Carlo method for detecting DA releasewhich E varies slightly above and below. During these equilibrium iterations, since the algorithmis producing a series of probable global states, the number of times a given voxel has si = +1 canbe converted to a percentage confidence, ci, that DA release exists in that voxel.The corresponding confidence map is then the image of the confidence at each voxel. Voxelswith ci 50 indicate high probability that there is no DA release. Voxels with ∼ 50 are likelyoften being changed, meaning it is uncertain whether or not there is DA release. Finally, voxelswith ci → 100 indicate increasing probability that there is DA release.Therefore, the MC approach provides a confidence map rather than a simple binary sensitivitymap for DA release. One may claim that calculating PDF (Fi) for each voxel would also do thetrick. However, the confidence map from the MC method factors in information from the localneighbourhood, whereas the operation of PDF (Fi) does not.In practice, one would still often require a defining boundary for a cluster of DA release. There-fore the confidence map should be thresholded at some critical confidence value ccrit such thatc >= ccrit are temporarily accepted. Using the same false positive cluster thresholding analysis inthe standard method, accepted voxels then pass through a cluster size threshold to be permanentlyselected. The result in still a confidence map, but voxels with c < ccrit or those a member of acluster below the cluster size threshold are set to zero.Presumably, higher confidence voxels have higher SNR with respect to the drop in RAC signalfrom baseline kinetics—either due to the denoising being more effective in this voxel, or the voxelhaving greater release in the first place.2.7.3 Optimizing the MC methodOne disadvantage of the MC method is that it requires some fine-tuning of the sole hyperparameter,β, to optimize its performance. Too large a β means the regularization term (keeping the rejectedvoxels following a true F-distribution) is more dominant over the parity term. This may result inoverconstraining the minimization process, resulting in less clustering of release voxels and hencelower sensitivity (albeit a lower false positive rate as well). On the other hand, too small a β meansthe parity term now dominates, and over-clustering will occur. This will potentially lead to highersensitivity (albeit at the cost of a higher false positive rate).We therefore can make the educated guess that β should be such that the two terms are moreor less evenly balanced. To do this, we applied the MC algorithm to simulations with DA releaseand kept track of the individual contributions of ∆P and ∆R to each ∆E. Histogramming thedistribution of changes at the end allows us to change β until the two distributions roughly havethe same variance. This occurred at β ∼ 8× 105 (results not shown). In Chapter 4, we test othervalues and evaluate the performance at each β, in case an asymmetric weighting proves beneficial.As well, since the algorithm employed in the MC method is non-deterministic, we want toensure that the metric of interest (detection sensitivity) has low variability after multiple trials ofthe same input (parametric F image). We ran 100 trials of the MC method applied to a singlenoisy realization of the cluster release simulation to test this.39Chapter 3MethodsTo summarize, the purpose of this work is to i) test various denoising protocols on simulatedtransient DA release data where the ground truth is known; this will allow us to objectively deter-mine the optimal denoising protocol when using the best current voxel-level transient DA releasedetection algorithm. We then ii) introduce the new MC detection algorithm, and show that it out-performs the best current algorithm, before iii) testing these methods on real human data wherethe ground truth is no longer known. We begin this section by describing the design and functionof the various simulation types to address the goals listed above.3.1 Generating a simulated phantomTo design the most realistic simulation possible, we base our simulated brain phantom on real RACdataset from a healthy control, scanned in the baseline state on the HRRT after an 8 mCi bolusinjection of RAC. The dynamic images were reconstructed using OP-OSEM [64] with 16 subsetsand 6 iterations. A frame-to-frame realignment in SPM (https://www.fil.ion.ucl.ac.uk/spm/) wasperformed to reduce the effects of subject motion. Moving from the real data to a noisy sinogramof a simulated RAC brain scan is done in five steps: segmentation, regional fitting, PSF filtering,forward-projection, and noise addition.3.1.1 SegmentationThe healthy control who was scanned in a baseline state also underwent an anatomical MRI scanon a Philips Achieva 3T scanner. The MRI image provides anatomical detail through T1 contrast.Running the Freesurfer segmentation algorithm [65–68] on the MRI image segments the imageinto 41 structures of the brain. The MRI image was then aligned to PET image6, defining atransformation from MRI space to PET space. This transformation is then applied to the segmentedMRI image such that the anatomical segments can be used to define the relevant anatomical regionsin the PET image.3.1.2 Regional fittingFor each segment, we averaged the voxel-level TACs to form a regional TAC. Each region was fitwith SRTM, using the cerebellum as the reference region. Using these fits, a noiseless baseline brainphantom was then created by filling the voxels of each segment with the corresponding SRTM fit.From here, the voxel-level TACs can be modified to include DA release as necessary by solving theintegral equation for lp-ntPET [equation (2.71)] iteratively, using the SRTM fit as the initial guess.In these voxels, we keep the release dynamics constant, using tD = 35 min, tP = 40 min, and α = 2in equation (2.72). The release magnitude is changed by tuning γ.6Specifically, to the image formed by summing all dynamic frames of the 4D dataset.403.1. Generating a simulated phantomThe fits are continuous functions, so we discretely sample these functions at times correspondingto the midpoint of each frame duration. The framing for our simulations is 3 × 1 minute framesfollowed by 24× 3 minute frames, for a total scan duration of 75 minutes.The segments only define tissue structure, so to add more detail to our simulated brain scanwe use the attenuation map from the PET scan (generated via a transmission scan using a Cs137rotating rod source) to identify regions of outer skin, bone, air cavities, and background regionin between the tissue structure. Bone and air cavities were assigned TACs of zeros, and skin andbackground region were assigned a TAC equal to the regional fit of the original data.These concentration values in each voxel can be scaled at this point to generate a noiseless brainphantom at any dose. The original phantom is considered a 8 mCi dose, reflecting the true dosegiven to the healthy control3.1.3 PSF filteringThe intrinsic resolution of the HRRT is defined through its PSF, which is simplistically modeledby an isotropic 2.5 mm full width at half maximum (FWHM) Gaussian kernel along the centralaxis. To incorporate the effects of the PSF in our simulation, we filtered our noiseless images withan isotropic 2.5 mm FWHM Gaussian kernel. The averages of all frames for both the non-filteredand filtered noiseless simulations (without any true release added) are shown in Figure 3.1.If the true object is considered to be the noiseless simulated brain phantom, then the filterednoiseless image serves as the ground truth; even with infinite counts the reconstruction couldonly reproduce the blurred noiseless phantom, since we do not model the PSF within any of ourreconstructions.3.1.4 Forward-projectionNow having blurred noiseless images, the noiseless sinograms can be generated by forward-projectingthese images using the system matrix of the HRRT. Now in sinogram space, the effects of attentua-tion and normalization can be incorporated into the simulated data by multiplying the attentuationsinogram from the healthy control’s scan, as well as by the normalization sinogram, reflecting thedetector sensitivities and detector gaps of the HRRT.3.1.5 Noise additionThe dominating source of noise in a real sinogram is from the Poisson process of the decay events.At this point each value in our simulated sinogram indicates a number of counts for each LOR.These values reflect the mean number of counts that would be recorded along each LOR after avery large number of repeated scans of the blurred noiseless phantom.However for any individual single scan, we would expect the number of counts to follow a Poissondistribution and hence have a variance equal to the number of counts as well. Therefore we takethe ith LOR from the noiseless sinogram, having Ni counts, and randomly sample from a Poissondistribution with mean Ni to arrive at N˜i, the counts for the ith LOR of the noisy sinogram.This process can be repeated several times to generate a series of unique noisy realizations ofthe sinogram, each of which are then reconstructed, denoised, and analyzed.413.2. Simulation variablesFigure 3.1: A single coronal slice, focusing on the striatum, from average of all frames of a simulated brainscan (no true PSF release added) before (a) and after (b) applying PSF filtering. The colour bar units arein Bq/cm3.3.2 Simulation variablesTo achieve the aims of this work we need to design simulations that modulate important variableswhich likely have a significant effect on whether or not DA release can be detected on the voxel-level. These include dose given to the participant being scanned, DA release level (i.e. DAR),time of the intervention, size of the region containing release, and location of the release within thebrain. A description of how these parameters are modulated is presented below:1. Dose. The higher the dose, the more counts recorded by the scanner, and hence the lowerthe noise level in the reconstructed images. The simulations of [60] suggest that there is aneffective lower limit of 16 mCi for voxel-level DA release detection. Therefore we generatecopies of each simulation at 16 mCi, 20 mCi, and 24 mCi7, to assess the performance ofeach denoising protocol as a function of varying dose level. We expect, intuitively, that ahigher dose level will lead to improved detection abilities as a result of lower statistical noise.However, we also note that dose levels are not the sole factor in determining statistical noise;patient size, metabolism, and scanner sensitivity also play important factors. Nonetheless,comparing the results of the dose levels tested will give a relative sense of how important doselevel is in determining DA release.2. DA release level. With higher DAR, the drop in the RAC TAC is larger and hence easierto distinguish from a baseline TAC. DAR can widely vary based on the type of intervention.Pharmaceutical interventions in rats have shown up to 10-fold increases in DA concentrationsabove basal levels [69], whereas smoking in humans may only induce a 2-fold increase [70, 71].Note that a large change in DA concentration results in a much smaller percentage changein measured binding potential. From double RAC studies, the correspondence may be ∼1:40 [72], meaning a 10-fold increase in DA concentration corresponds to DAR ∼ 25%, and a2-fold increase to DAR ∼ 5%. Thus to cover a reasonable range of DA release levels, in oursimulations DAR spans 3%-21%.724 mCi is chosen as an upper limit to reflect radiation safety considerations which limit the doses given to researchvolunteers423.3. Types of simulations3. Time of intervention. The simulations of [60] have determined tD = 35 min to be theoptimal time to present an intervention to the participant. This allows full uptake of theradiotracer—effectively allowing the system to come to equilibrium—before perturbing thesystem by means of an intervention. Therefore for all simulations we keep tD constant.4. Size of region containing release. As mentioned before, as physiological structures ap-proach the scanner resolution, they become significantly degraded by PVE. However, evenfor larger structures with sharp boundaries, such as a cluster of DA release surrounded bytissue with no release, PVE will affect the edges by effectively mixing together release voxelsnear the edges with non-release voxels surrounding the cluster. Our simulations model theHRRT, one of the highest resolution dedicated brain scanners. The PSF of the HRRT ismodelled with a 2.5 mm FWHM Gaussian kernel, representing the effective resolution of thescanner. The voxel size of the HRRT is isotropic with a side length of 1.21875 mm. Thereforerelease voxels located within a few voxels of the cluster edge will be noticeably contaminatedby non-release voxels. Previous simulations [57, 60] generate clusters of 16, 32, 64, and 128isotropic 2 mm voxels, corresponding to 71, 142, 283, and 566 HRRT-sized voxels in size8.However, these simulations did not include the PSF of the scanner. Therefore in our clustersimulations, we elect to generate clusters of 300, 400, and 500 voxels, since initial simulationswith smaller clusters and modeling the PSF produced extremely low detection rates.5. Location of release. The striatum has a high density of D2 receptors, and hence is theregion where DA release is most expected to be detected. The striatum can be subdividedinto the L/R caudate and L/R putamen. The caudate, being a smaller region, is more affectedby PVE. Therefore we hypothesize that DA release in this region will be harder to detect.To test this, regions of release in various locations of both the caudate and putamen in thesimulations.3.3 Types of simulationsA baseline simulation (with no true DA release incorporated) is designed to determine the clusterthreshold size as detailed in Section 2.6. Two further simulation types are designed to test thevariables above. As mentioned, for each of the three simulation types, three copies are made atdoses of 16 mCi, 20 mCi, and 24 mCi. For each of these, 50 noisy realizations with independentnoise patterns are generated for increased statistical power during analysis.3.3.1 Baseline simulationThe baseline simulation, as evident from its name, has no voxels with DA release. The other twosimulations both contain DA release, but are used to achieve separate goals.3.3.2 Random release magnitude simulationThis simulation contains DA release in every voxel of the L/R caudate and L putamen. The R puta-men is devoid of release to verify false positive thresholding—through both the initial significancethresholding (either using the Fcrit in the standard method or ccrit in the MC method) and theCST—is performing to the 90% false positive-free rate established from the baseline simulations.8Hereafter, any mention of voxel refers to a HRRT-sized voxel of 1.21875 mm by 1.21875 mm by 1.21875 mm.433.3. Types of simulationsFigure 3.2: Axial, coronal, and sagittal views of the DAR values within the striatum for the random releasemagnitude simulation.Figure 3.3: Distribution of voxel-level DAR before and after PSF filtering. The slight upward trend inthe distribution of DAR before PSF filtering is because the conversion of the uniformly-sampled γ to DAR(using the lookup curves of Figure 2.11) is slightly sublinear at higher γ. The falloff at the end is becauseonly putamen voxels can produce DAR > 24% for γ ∈ [0, 12]×10−4. The narrowing of the DAR distributionafter PSF filtering is from low release voxels having a greater likelihood of being surrounded by higher releasevoxels, and vice versa.The γ in each DA release voxel was sampled randomly in the range γ ∈ (0, 12) × 10−4 s−1.After the PSF filtering, nearby voxels are correlated and hence the release levels still vary spatially,but become locally more homogeneous. We then calculate the ground truth DAR from the newvoxel-level TACs (i.e. after the PSF filtering). The ground truth parametric DAR image is shownin Figure 3.2.The purpose of this simulation is to quantify the detection sensitivity (i.e. the percentage ofvoxels with DA release that are detected by the algorithm) as a function of DAR. Due to thegeometry of the striatum, certain regions are more affected by the PVE than others. By makingthe release pattern spatially invariant, we attempt to ensure that the release levels—on average—are equally affected by the PVE, such that any differences in sensitivity are solely due to DAR.An ideal way to achieve this would be to perform a separate simulation for each DAR. However,due to the large range of potential DAR, this would require an impractical number of simulations,so the compromise was made to incorporate the full range of DAR within a single simulation type.443.3. Types of simulationsA consequence of this compromise is that lower release voxels are more likely to be surroundedby higher release voxels, since sampling from a uniform distribution of release levels means thereexists a larger number of release levels higher than some low value. For instance, if one sampleplaces a value of γ = 3 × 10−4 s−1 in a voxel, there is 75% chance that the neighbouring samplewill have γ > 3 × 10−4 s−1, since all γ ∈ (0, 12) × 10−4 s−1 are equally likely. Therefore, afterthe PSF filtering, such a voxel on average will have its release level “boosted” by being correlatedwith adjacent higher release voxels. Conversely, voxels originally having relatively large γ will onaverage have their release level decreased by adjacent lower release voxels. This in effect narrowsthe distribution of release level in the simulation, as evident from Figure 3.3.After reconstructing a noisy realization of a random release magnitude simulation and applyingfiltering as a denoising operation, further correlation will be introduced between the high and lowrelease voxels, which likely will have the effect of boosting detection sensitivity for low release voxels(as compared to a region of uniform low release), and lowering the sensitivity for high release voxels(as compared to a region of uniform high release).3.3.3 Cluster release simulationThe second DA release simulation places six clusters within the L/R caudate and L/R putamen:one in each side of the caudate, and two in each side of the putamen. Amongst the six clusters, wemodulate cluster size and release level (c.f. Table 3.1. Since the caudate is more affected by thePVE than the putamen (being a smaller, more narrow structure), an observation can be made onthe effect of PVE with regards to sensitivity.Table 3.1: Cluster properties in the cluster release simulation. DAR specified here is the mean DAR,computed after PSF filtering, of all voxels within the cluster.ID Location Size [voxels] Size [cm3] DARC1 L putamen 500 0.91 15%C2 L putamen 300 0.54 7.5%C3 L caudate 400 0.72 15%C4 R caudate 400 0.72 7.5%C5 R putamen 300 0.54 15%C6 R putamen 500 0.91 7.5%Multiple slices of the striatum showcasing the voxel-level DAR values for each cluster, bothbefore and after PSF filtering, is shown in Figure 3.4. It is evident that PSF filtering results in adistribution of voxel-level DAR within each cluster; thus we have taken the cluster-level mean tocharacterize each cluster.As well, PSF filtering results in spillover beyond the original boundaries of the clusters, ef-fectively introducing release in a region where no release exists in reality. Since the PSF is anunavoidable consequence of imaging, we neither punish nor reward the detection algorithms fordetecting DA release in these voxels with spillover release; that is if release is detected in such avoxel, then it does not contribute to the overall sensitivity of the method nor to the false positiverate.453.4. Reconstruction and denoising protocolsFigure 3.4: Axial, coronal, and sagittal views of the striatum for the cluster release simulation, indicating(a) the cluster locations and corresponding cluster IDs; (b) the voxel-level DAR before PSF filtering; and (c)the voxel-level DAR after PSF filtering. A decrease in DAR near the edges of the original cluster boundariesas well as spillover of release outside the original boundary can be seen.3.4 Reconstruction and denoising protocolsFor each noisy sinogram realized for a given simulation type, we reconstruct and denoise the imageswith a variety of protocols to assess the effect of denoising on cluster threshold size, sensitivity, andparametric accuracy.3.4.1 ReconstructionsThe most commonly-employed iterative reconstruction method is OP-OSEM [64]. For our simu-lations, we adopt the protocol of our centre, reconstructing each noisy sinogram with OP-OSEMusing 16 subsets and 6 iterations. We compare this method with the denoised reconstructionHYPR-AU-OSEM. Our centre has determined the optimal number of HYPR-AU-OSEM iterationsusing a 2.5 mm FWHM Gaussian kernel size in the HYPR operator. This optimization has beenperformed in terms of locating the optimal noise-bias tradeoff. Since the denoising ability of thereconstruction is variable on the number of counts in the composite image, the optimal iteration463.4. Reconstruction and denoising protocolswill also be a function of counts. The following empirical formula has been derived for the optimalnumber of iterations, Nit, given the total number of counts, Nc, in a sinogramNit(Nc) = d2.87× 10−7 ·Nc + 5.35e (3.1)As an example, in our 16 mCi baseline simulations we have Nit = 6 for the lowest count frames,and Nit = 18 for the highest count frames.3.4.2 Post-processingPrevious work looking into optimizing transient neurotransmitter release protocols found significantimprovements in sensitivity when using HYPR post-processing [57, 60]. They elected to use a 3 x3 x 3 boxcar kernel (2 mm voxels) within the HYPR operator.Based on the results of [44], we employ HYPR post-processing using Gaussian kernels. A widerkernel will better denoise the image, since the filtered target image will be less noisy after filteringwith a wider kernel. As well, the difference between the composite and filtered composite will begreater for high frequency features, allowing one to better recover such features. However, too widea kernel may oversmooth small features in the filtered target. Therefore we test two kernel widthsof 2.5 mm and 5 mm FWHM for each reconstruction type (recall our voxels are isotropic with1.21875 mm side length).HYPR post-processing alone is limited by the noise in the composite image. Therefore we tryan additional denoising approach where the images are filtered before applying the HYPR operator,such that the composite itself is better denoised. We filter with 2 mm and 5 mm Gaussian kernels tocompare narrow versus wide kernels. Since a wider filtering kernel will better denoise the images butwill induce a larger bias, we hypothesize that the sensitivity will be higher for the wider kernel, butthe parametric accuracy will be quite poor. These filtered images subsequently post-processed withHYPR can be considered similar in their approach to HYPR-AU-OSEM + HYPR post-processing,since HYPR-AU-OSEM will denoise the raw images within the reconstruction—thus producing aless noisy composite—but without inducing much bias.In total then we have eight separate reconstruction plus denoising protocols to test on eachsimulation type. They are summarized in Table 3.2, along with the abbreviations we will use forthese protocols for the remainder of this work.3.4.3 Comparison to previous noise-modeling methodsPrevious work [57, 60] incorporated noise into their simulations by adding Gaussian noise to noise-less simulated TACs (i.e. without simulating the effects of the PSF, normalization, attenuation,and reconstruction). Their noisy simulations were then denoised using HYPR post-processing witha 6 mm by 6 mm by 6 mm boxcar kernel. Noise for each time point is sampled from a Gaussiandistribution having a standard deviationσ = µeλt√(CT (t)× e−λt)/∆t (3.2)where λ is the decay constant of the tracer and ∆t is the frame duration. We also add noise to ournoiseless datasets in this manner to both evaluate consistency with published results and the effectof our more accurate data modeling on detection sensitivity outcomes.In their simulations, µ = 25 (min·Bq/cm3)1/2 was used. However, since our voxel sizes differwe set µ = 70 (min·Bq/cm3)1/2 to match the standard deviations observed in a 20 mCi RAC scanon the HRRT with 1.22 mm isotropic voxels.473.5. Analysis workflowTable 3.2: Reconstruction and denoising protocols tested in this work. “NF” = no filtering before applyingHYPR. Initial kernel is either the filtering kernel width for OSEM reconstructions, or the kernel width forthe HYPR operator within HYPR-AU-OSEM.Abbreviation ReconstructionInitial HYPRkernel kernelOSEM-NF-2.5mmOSEM-2.5 mmOSEM-NF-5mm 5 mmOSEM-2mm-2.5mm2 mm2.5mmOSEM-2mm-5mm 5 mmOSEM-5mm-2.5mm5 mm2.5 mmOSEM-5mm-5mm 5 mmHYPR-2.5mm-2.5mmHYPR-AU-OSEM 2.5 mm2.5 mmHYPR-2.5mm-5mm 5 mmFrom here, we apply HYPR post-processing with a 5 mm Gaussian kernel, to keep consistentwith the kernels used to post-process our reconstructed simulations. These simulations generatedwith Gaussian noise are then analyzed identically to the reconstructed simulations, and the resultsare compared.3.5 Analysis workflowThe analysis of each simulation type will differ slightly, though at the core essentially the DA releasedetection algorithm remains the same. The theory of DA detection is outlined in Section 2.6. Thistheory forms the basis of the algorithm. Given a denoised image, we collect the voxel-level TACsonly in the striatal region (defined by the segmented MRI image). The algorithm runs in two steps:a fitting step and a thresholding step.3.5.1 FittingFor each voxel-level TAC in the striatum:1. Fit SRTM and lp-ntPET to the TAC and compute the WRSS.2. Given WRSSSRTM and WRSSlp−ntPET, calculate Fν1,ν2 using equation (2.79).SRTM, in the integral form of equation (2.60), is fit using weighted least squares. For lp-ntPET,the transient term of equation (2.71) needs to be pre-computed using a set of basis functions for h(t).These basis functions should span the realistic regions of the parametric space of θ′ = {tD, tP , α}.Given an intervention time of 35 minutes, in a real scan we would constrain tD ∈ [30, 50] min in1.5 minute increments (i.e. half a frame duration) to account for any anticipatory or delayed-onsetrelease. Since we must have tP > tD, we constrain tP ∈ [tD + 1.5, tend − 5] minutes, again in 1.5minute increments. Finally, we allow α ∈ [0.25, 1, 4]; a small number of values to keep the numberof basis functions from growing too large and ultimately resulting in lengthy computation times andpotential over-fitting. This range also intentionally does not include ground truth basis functionof θ = {tD = 35 min, tP = 40 min, α = 2}, since the probability of specifying the correct basisfunction a priori is exceedingly low.Permutations of these parameters result in 819 unique basis functions. For each basis function,weighted least squares fitting is applied to determine the best fit for R1, k′2, k2a, and γ. Of these483.5. Analysis workflow819 best fits, the one with the lowest WRSSlp−ntPET is selected to determine the best fitting tD,tP , and α.3.5.2 Weighting schemesIn weighted least squares fitting of both SRTM and lp-ntPET, we attempt to provide the bestlinear unbiased estimator (BLUE) of the parameters in the respective models. It can be shownthat the BLUE is produced when the weighting scheme used is the reciprocal of the variance. Ina voxel-level TAC at a given time point t, the variance cannot be known because it would requiremany identical scans—save for a new random noise pattern. However, we know that the variancein the number of counts during a single frame, N(t), in raw sinogram for the frame correspondingto that time point goes as N(t). After backprojection in the reconstruction, the noise in the imageis no longer purely Poisson. But if we assume that the noise level in a target region voxel roughlyscales as the Poisson noise level in the total frame, then the weights should bew(t) =1N(t)(3.3)However, the model is fitted to a corrected version of the raw reconstructed image. One ofthese corrections is for tracer decay since the beginning of the scan, since radioisotope decay doesnot imply the molecule concentration is decreasing. If the initial activity of the radiotracer with ahalf-life t1/2 at the start of the scan (t = 0) is A0, then the activity at time t during the scan isA(t) = A0e−λt (3.4)where λ = ln 2/t1/2 is the decay constant. During a frame starting at t1 and ending at t2, thenumber of counts measured is found by integrating this decaying activity over the frame durationA(t1 → t2) =t2∫t1A(t)dt = A0t2∫t1e−λtdt =A0λ(e−λt1 − e−λt2)(3.5)If there was no radioactive decay (t1/2 → ∞ =⇒ λ → 0), the counts measured during theframe would simply bet2∫t1A0dt = A0∆t (3.6)where ∆t = t2 − t1. Since we want our counts to be corrected for decay, the decay correction isnothing but the ratio of the integrated non-decaying activity to the integrated decaying activityDC(t1 → t2) = A0∆tA0λ (e−λt1 − e−λt2) =λ∆teλt11− e−λ∆t (3.7)This corrects the measured counts during the frame. Our TACs, however, are concentrationvalues as a function of time, having units of Bq/cm3. Thus the corrected counts need to be dividedby the frame duration, ∆t, to arrive at units of Bq, and multiplied by a global calibration factorthat accounts for scanner sensitivity. The activity within a voxel is then divided by the voxelvolume to finally arrive at units of Bq/cm3 for the voxel-level CT (t).These factors not only apply the measured value in the voxel, but also to the variance as well,meaning we should update our weights with the product of all the above listed factors. However,493.5. Analysis workflowbecause the global calibration factor and voxel size are spatially- and temporally-invariant, theycan be ignored in our weights. Thus the final form of our weights isw1(t) =(∆t)2N(t) [DC(t)]2(3.8)These are the weights proposed originally by [52] in fitting SRTM using the basis functionmethod, and we shall refer to this weighting scheme as “traditional weighting”. The weightsintuitively make sense when only using the information contained within a single frame. However,upon applying HYPR post-processing, we are incorporating information from all frames via thecomposite. Since the same composite is used for all frames, we may na¨ıvely think that all framesshould have similar noise levels, and hence uniform weights should be usedw2(t) = 1 (3.9)This weighing scheme will be referred to as “uniform weighting”. However, the composite iscomposed of a sum of the decay-uncorrected images (in order to maximize SNR—the decay-correctedimages amplify the noise levels in later frames). Hence the decay-uncorrected, HYPR post-processedimages will have roughly the same noise level across time. Therefore uniform weighting is mostappropriate in the decay-uncorrected domain. Since the TACs we are fitting to are from the decay-corrected images, to achieve this we simply undo the decay correctionw3(t) =1DC(t)(3.10)This weighting scheme will be referred to as “uncorrected weighting”. A plot comparing eachweighting scheme (normalized to one) for the baseline simulation is shown in Figure 3.5. Weseparately apply these three weighting schemes to each simulation. The optimal weighting schemecan then be determined with respect to detection sensitivity, cluster size threshold, and parametricaccuracy.Figure 3.5: The three weighting schemes tested in this work, normalized to one. The jump seen intraditional weighting after the third frame is from the change in frame duration from one to three minutes.503.5. Analysis workflow3.5.3 ThresholdingAfter performing voxelwise fitting within the striatum, we can compute Fν1,ν2 at each voxel toarrive at a parametric F image. In the standard method, the first threshold keeps only voxels withFν1,ν2 ≥ Fcrit. For our 27 frame protocol, we have ν1 = 4 and ν2 = 20, yielding Fcrit = 2.866. Voxelswith F values higher than this threshold violate the null hypothesis that SRTM and lp-ntPET fitequally well with 95% confidence. By placing a 1 in each of these voxels and zero everywhere else,we form an initial binary significance mask.When the MC method is being used, it is applied as described in Section 2.7.2 with the para-metric F image as the input, and the parametric c image as the output. From here, thresholdingat c ≥ ccrit is performed to compute the initial binary significance mask.Next, for either method we remove any significant voxels whose best fit has γ < 0, since this isunphysical when DA release is expected. γ < 0 implies a positive perturbation to the TAC duringthe post-intervention period. Since we only expect negative perturbations, the only way for this fitto occur is if lp-ntPET is fitting to noise in the TACs.Next, clusters in the binary significance mask are located and their sizes are computed. Clusterswith sizes below the CST are then rejected. The CST is determined from the baseline simulations,and is explained in more detail below. All remaining clusters are considered clusters of significantrelease, and thus comprise the final binary significance mask.3.5.4 Determining the cluster size thresholdWhen either the standard or MC method is applied to the baseline simulations themselves, a clustersize threshold of 1 is used at the last thresholding step to keep all significant voxels (evidently falsepositives, since this is a baseline scan) having γ > 0. After applying this to 50 baseline simulations,we compute the CDF of the distribution of false positive cluster sizes. The cluster size that removes99% of the clusters from all 50 baseline simulations is used as an initial cluster size threshold. Fromhere, the cluster size threshold is increased until 90% of all baseline simulations (i.e. 45 out of 50)are completely free of false positive clusters. This number forms the final CST when analyzing theother simulation types.3.5.5 Assessing detection sensitivityIn non-baseline simulations, the detection sensitivity can be computed from the binary significancemasks and ground truth knowledge. As mentioned, detection sensitivity is a measure of howfrequently true release is detected in the data (i.e. the true positive rate). We can define sensitivityin a number of ways.The first way is to look at every voxel with true release, and see what percentage of the timethis voxel is classified as a release voxel by the algorithm. This is referred to as voxelwise sensitivity[60]. For instance, if the voxel is classified as a release voxel 44 out of the 50 noisy simulations,then we assign a voxelwise sensitivity of 88% to that voxel. This metric is most useful whendetermining detection sensitivity as a function of release level, looking at the spatial patterns ofdetection sensitivity within a cluster, and in the overall detection sensitivity in the random releasemagnitude and cluster release simulations (i.e. by averaging the voxelwise sensitivity over a truecluster).Alternatively, from the point of view of a researcher applying this method, we may not be toodisappointed if the algorithm only detects a certain portion of the voxels within a true releasecluster. Therefore we can define binary sensitivity [60] as what percentage of the time any portionof the cluster is detected, and binaryX sensitivity as what percentage of the time X% of the voxels513.5. Analysis workflowwithin a true cluster are detected. Since 50% is a reasonable number to assess a pass/fail rate forthe algorithm in a given true cluster, we use binary50 sensitivity as a metric in this work.3.5.6 Parametric accuracyThe most obvious benefit of a voxel-level approach for detecting transient DA release is to allowthe detection algorithm to identify the specific regions of the brain having any DA release, ratherthan pre-specifing an ROI and searching for DA release on a regional level. In a way, the voxel-level approach defines ROIs (via the binary significance mask) which can be subsequently used forregional parametric analysis.In a study employing this algorithm for analysis, the ROIs returned can shed light on neuro-physiological processes. Different subregions of the striatum are involved in various neurologicalprocesses, such as cognition, motor function, controlling reward systems, and executive function[73–79]. Therefore the location and size of the clusters detected serve as a qualitative metric ofhow active these processes are in response to the intervention.To move into a quantitative regime, we need to study the voxel-level as well as the regional (i.e.cluster-level) parametric accuracy. Knowing where DA release occurs in the brain is important,but quantitatively comparing the magnitude and/or the temporal pattern of the release would beuseful in a study where the effect of the intervention is compared across multiple cohorts, or wherethe effect of different interventions is the purpose of the study. For instance, if a gambling taskserves as the intervention, what is the release magnitude (and where, spatially, does this releasetend to occur) in gamblers versus non-gamblers?The parametric accuracy can be visually assessed by comparing the detected parameter dis-tributions in a detected cluster with the true distributions. For a quantitative comparison, themean absolute error (MAE) can be used. MAE averages over the absolute difference between theparameter estimate at the jth voxel, θˆj , and its ground truth value θjMAE(θ) =1NvNv∑j=1|θˆj − θj | (3.11)where Nv is the number of true release voxels detected. MAE thus gives an indication of the averagevoxel-level deviation of the parameter in question from its ground truth value.In our opinion, with regards to studying system responses to a stimulus, the most valuableinformation would be how long after the stimulus does release commence (or before if there areanticipatory effects), when does this response reach its peak magnitude, and, more generally, whatis the time-varying magnitude of the response? The first two questions can be answered by theestimated tD and tP , respectively, and thus we quantify their average voxel-level accuracy throughMAE. However, because we do not expect changes in release timing in small regional neighbour-hoods (such as the cluster sizes simulated in this work), we also assess the cluster-level meanestimates of tD and tP , since these estimates will likely be closer to the ground truth than thevoxel-level estimates and consequently may be more informative within the analysis of a study.To estimate the release magnitude of the response as a function of time, we can analyze theparametric accuracy of BP (t) and DAR(t). This will be elaborated on in the next section.Due to the nature of the fitting portion of the analysis method, the voxel-level estimates of tDand tP are constrained to the pre-defined values within the basis functions, h(t). This approachleads to two issues for tD: the ground truth tD = 35 mins does not fall on a frame midpoint nor anavailable tD for fitting, pre-specificed in the lp-ntPET basis functions. The closest frame midpointafter the true release time (i.e. data point on a TAC where the negative perturbation due to DA523.5. Analysis workflowrelease has begun) is 37.5 minutes, and the closest pre-specified tD for fitting are tD = 34.5 ortD = 36 minutes. Therefore in the refitting procedure it may be expected for the refitted tD tobe positively biased. Though this trouble could be rectified by specificing the true tD to land ona frame midpoint and pre-specified tD for fitting, in a realistic scan the probability of having theonset of DA release land exactly on a frame midpoint is extremely unlikely. Similarly for the groundtruth tP = 40 minutes, the closest frame midpoint is 40.5 minutes, and the closest pre-specified tPfor fitting is also 40.5 minutes, making this issue less pervasive for tP . Though this effect couldbe reduced by decreasing frame duration, the noise in the voxel-level TACs will increase, likelylowering detection sensitivity.3.5.7 Comparing DA release curvesWhile the above methods compare a single parameter at a time, the metric DAR(t) of equa-tion (2.74) is a combination of all parameters within the lp-ntPET transient term. Thereforecomparing these DA release curves to their ground truth curves is an effective way to assess theaccuracy of all parameters simultaneously, as well as capturing the interplay between parameters.For instance, since the true α = 2 is not specified in the basis functions, by definition an incorrectα will be selected (again because the likelihood of guessing the ground truth parameters in a scana priori is very unlikely), and the remaining parameters tD, tP , and γ will have to compensate insome way to provide an adequate fit.Since the ground truth image is the noiseless image filtered by the PSF, as mentioned in Sec-tion 3.3 the ground truth γ required for the calculation of ground truth DAR and DAR(t) aredetermined by applying the algorithm to this ground truth image (i.e. refitting to the noiseless,PSF-filtered images with lp-ntPET). The PSF filtering modifies the shape of voxel-level TACswith DA release by replacing these noiseless unfiltered TACs by a weighted average of the neigh-bouring TACs—the neighbourhood defined by the kernel size of the filter. Voxels near the edgesof a DA release cluster are therefore more affected than voxels near the centre of the cluster, sincethe neighbourhood of an edge voxel contains a large proportion of non-release voxels, whereas theneighbourhood for a central voxel is almost exclusively release voxels.The overall effect for a cluster of release is that after PSF filtering, the refitted voxel-level γare negatively biased. Based on the reasoning above, edge voxels in a cluster have larger negativebiases than central voxels. Therefore, after PSF filtering, the initial γ placed at every voxel in thecluster is transformed to a distribution of γ within the cluster having a mean less than the initialvalue. As mentioned before, the PSF filtering also causes “spillover”, meaning non-release voxelsnear the edge of a DA release cluster in the unfiltered noiseless image become release voxels (albeitwith low γ) after filtering. We do not assess the parametric accuracy in these voxels.3.5.8 Parametric thresholding: A more representative picture of releasedynamicsAs a part of our analysis of voxel-level timing parameter accuracy, we will compute the MAE intD and tP to determine, on average, how well the voxel-level timing parameters deviate from theground truth values. Voxels that have estimated timing parameters that exceed these MAE valuesare therefore less informative. In a real scan, however, we will not have the ground truth valuesused in the computation of MAE. Nevertheless, instead we can compute the average deviationof the estimated parameter of the ith detected voxel from the estimated parameters of its local533.5. Analysis workflowneighbourhood within the clusterAD(θˆi) =1||Ni||0∑j∈Ni|θˆi − θˆj | (3.12)where Ni is the set of up to the 124 nearest neighbours in 3D for the ith detected voxel. We canthen express θˆi as the (unknown) ground truth value, θi, plus some error iAD(θˆi) =1||Ni||0∑j∈Ni|θi + i − θˆj | (3.13)Using the triangle inequality, we haveAD(θˆi) ≤ 1||Ni||0∑j∈Ni[|θi − θˆj |+ |i|](3.14)The first term in the summation is the MAE within Ni, so we can sayAD(θˆi) ≤MAE + |i| (3.15)Thus in the limit that θˆi → θi =⇒ i → 0, the absolute deviation is bounded by MAE9.Estimates with more error have absolute deviation bounded by a number greater than MAE. If weassume the MAE computed on Ni is reasonably close to the MAE from our simulations, then wecan exclude voxels—from subsequent parametric analyses—having values of AD(θˆi) greater than acertain predetermined value. These excluded voxels still are informative in terms of spatial extent ofthe cluster, but are not sufficiently parametrically accurate to warrant their inclusion of computingthe meaningful cluster-level metrics of BP (t) and DAR(t) for use in a study.The parametric thresholding approach (PTA) proceeds as follows. Compute AD(tˆD) andAD(tˆP ) at each detected voxel. Based on the MAE results from our simulations (as will be dis-cussed in Section 4.5), we employ voxelwise upper limits of AD(tD) ≤ 4 min and AD(tP ) ≤ 6 min.Voxels obeying both of these limits are kept for computing cluster-level BP (t) and DAR(t).We note that in the case where true gradients in tD and tP exist, this method may rejectvoxels with fairly accurate timing parameters. However, this should be a rare occurrence, since itis unlikely that gradients on the level of our sliding window neighbourhood of 125 voxels (0.23 cm3)will be as large as the above thresholds of AD(tD) = 4 min and AD(tP ) = 6 min.Another instance where PTA may prove important is in separating nearby clusters that havedistinct release timing behaviour. As mentioned, PSF filtering spatially correlates nearby voxels,having the effect of spreading signal of DA release into voxels that actually do not have release.Further, reconstruction and post-processing will introduce additional spatial correlation. Thereforetwo clusters that, for instance, have tD separated by 10 mins, but are spatially located nearby oneanother, may be merged together and detected as a single cluster. PTA—in theory—will be able toseparate these clusters, since voxels in the region between the clusters will deviate significantly onaverage from the voxels in their local neighbourhood. This should not be of use in our simulations,since we simulate all clusters with identical timing parameters, but could be of use in real datasetswhere clusters have unique release timing behaviours.9In fact, it can be seen that the equality holds in this case.543.6. Application to real human data3.6 Application to real human dataA healthy control was given 20 mCi bolus dose of RAC at the beginning of a 75 minute scan. Thiswas a pilot scan as a part of a planned gambling study, using the mid-scan intervention protocol.The participant was presented with a gambling task at 36 minutes into the scan, which lastedapproximately 15 minutes. A detailed description of the full study and its aims is outside thescope of this thesis work; the pilot data were used exclusively to test our approach to DA releasedetection.The images were reconstructed with both OP-OSEM and HYPR-AU-OSEM. The data wereframed as follows: 1 minute x 5 frames, 2 minutes x 5 frames, and 3 minutes x 20 frames. The twobest denoising protocols (as determined from our simulations) were applied, and the resulting datawere analyzed using both the standard and MC methods. Since the ground truth is not known, wecannot assess the detection sensitivity, false positive rate, nor the parametric accuracy. However,we can compare the various reconstructions and analysis methods and provide possible explanationsfor any significant differences, as well as relate these differences to our simulated results.55Chapter 4Results and discussionWith many variables being tweaked—simulation design, reconstruction type, post-processing, andanalysis method—a lot of simulated data was analyzed and thus a road map is necessary to guidethe reader through the results.Since the both the standard and MC methods require a CST, we begin the results by presentingthe CSTs determined for each denoising protocol, and plotting them as a function of dose. Theseresults are presented for all three different weighting schemes.With a CST determined, the standard method can be applied to the random release magnitudesimulations to assess the sensitivity as a function of release level. As well, since we simulate notrue release in the R putamen of the random release magnitude simulation, we also check that theCST is performing adequately. Again, this analysis is repeated at the three dose levels, each usingthe three separate weighting schemes. From these data, we can eliminate some denoising protocolsthat perform unsatisfactorily, as well as decide on the optimal weighting scheme moving forwardfor the last simulation.This leaves us with the cluster release simulation to analyze. Here we compute the voxelwise andbinary sensitivity as well as the parametric accuracy for each cluster, for all remaining denoisingprotocols, at the three dose levels. The parametric accuracy analysis is subdivided into comparingthe cluster-level mean timing parameters (tD and tP ), the MAE of the timing parameters computedover each cluster, and visual comparisons of estimated BP (t) and DAR(t). We then employ theparametric thresholding approach to recalculate BP (t) and DAR(t) with the subset of survivingvoxels from the thresholding, for all remaining denoising protocols, at a single dose level.From here we briefly summarize the procedure of optimizing the hyperparameter, β, in the MCmethod, and arrive at two distinct versions of the MC method. These methods are applied to therandom release magnitude and cluster release simulations and compared to the results from thestandard method, at a single dose level.Finally, the standard and MC methods are applied to the real human RAC dataset from thegambling study pilot scan. Comparisons are made between methods, and any similarities or differ-ences are checked for consistency with the results from our simulated data.4.1 Baseline simulation4.1.1 Cluster size thresholdsTables 4.1-4.3 show the CST for each denoising protocol that result in 45 out of the 50 noisyrealizations being completely free of false positives, for all three weighting schemes, for doses of 16mCi, 20 mCi, and 24 mCi respectively.The effect of weighting schemeFor all protocols, using uncorrected weights produces the smallest CST. The traditional weightsproduce, on average, 54% larger CSTs than the uncorrected weights at 16 mCi, 49% larger at 20564.1. Baseline simulationTable 4.1: CST for a 16 mCi bolus dose. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrected[voxels] [cm3] [voxels] [cm3] [voxels] [cm3]OSEM-NF-2.5mm 42 0.076 54 0.098 23 0.042OSEM-2mm-2.5mm 59 0.107 97 0.176 39 0.071OSEM-5mm-2.5mm 239 0.433 302 0.547 153 0.277HYPR-2.5mm-2.5mm 63 0.114 94 0.170 42 0.076OSEM-NF-5mm 178 0.322 270 0.489 119 0.215OSEM-2mm-5mm 200 0.362 294 0.532 141 0.255OSEM-5mm-5mm 370 0.670 436 0.789 214 0.387HYPR-2.5mm-5mm 141 0.255 218 0.395 108 0.196Table 4.2: CST for a 20 mCi bolus dose. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrected[voxels] [cm3] [voxels] [cm3] [voxels] [cm3]OSEM-NF-2.5mm 38 0.069 41 0.074 20 0.036OSEM-2mm-2.5mm 54 0.098 74 0.134 39 0.071OSEM-5mm-2.5mm 185 0.335 245 0.444 125 0.226HYPR-2.5mm-2.5mm 58 0.105 76 0.138 42 0.076OSEM-NF-5mm 142 0.257 198 0.358 107 0.194OSEM-2mm-5mm 167 0.302 228 0.413 118 0.214OSEM-5mm-5mm 309 0.559 473 0.856 187 0.339HYPR-2.5mm-5mm 132 0.239 167 0.302 97 0.176Table 4.3: CST for a 24 mCi bolus dose. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrected[voxels] [cm3] [voxels] [cm3] [voxels] [cm3]OSEM-NF-2.5mm 34 0.062 43 0.078 25 0.045OSEM-2mm-2.5mm 58 0.105 67 0.121 37 0.067OSEM-5mm-2.5mm 229 0.415 269 0.487 131 0.237HYPR-2.5mm-2.5mm 56 0.101 78 0.141 42 0.076OSEM-NF-5mm 153 0.277 220 0.398 98 0.177OSEM-2mm-5mm 202 0.366 238 0.431 123 0.223OSEM-5mm-5mm 344 0.623 430 0.778 235 0.425HYPR-2.5mm-5mm 134 0.243 188 0.340 98 0.177mCi, and 51% larger at 24 mCi. The uniform weights perform much more poorly, with increases inCSTs over uncorrected weights of 118%, 97%, and 92% at 16 mCi, 20 mCi, and 24 mCi, respectively.Since lp-ntPET is constrained to only fit to temporal patterns characteristic of release for t ≥ 30min, any overfitting that results in a false positive will occur during this time frame. Therefore,it is not surprising that the weighting scheme that penalizes most heavily in this time frame—574.1. Baseline simulationuncorrected weighting (c.f. Figure 3.5)—results in the lowest CST; overfitting is discouraged bynot weighting the WRSS of the fit as much in later times as compared to earlier times.The effect of HYPR post-processing kernel sizePredictably, using a 2.5 mm kernel in HYPR post-processing also produces lower CST than a 5 mmkernel; wider kernels introduce more spatial correlation between voxels, meaning clusters containingTACs with temporal patterns characteristic of false positives are extended in size. This trend is quitesimilar for all dose levels and relatively stable across weighting schemes; one can expect the CST toincrease by ∼ 370% when changing from OSEM-NF-2.5mm to OSEM-NF-5mm, by ∼ 230% whenchanging from OSEM-2mm-2.5mm to OSEM-2mm-5mm, by ∼ 60% when changing from OSEM-5mm-2.5mm to OSEM-5mm-5mm, and by ∼ 130% when changing from HYPR-2.5mm-2.5mm toHYPR-2.5mm-5mm. Though this increase is smallest for OSEM with 5 mm filtering before HYPRpost-processing, the absolute size of the cluster sizes thresholds—no matter what kernel size is usedin HYPR post-processing—is the highest amongst the protocols.The effect of reconstruction method and initial filteringIn terms of reconstruction methods, if 2.5 mm HYPR post-processing is used then OSEM withno filtering always produces the smallest CST, followed by OSEM with a 2 mm initial filtering,then HYPR-AU-OSEM, and finally OSEM with a 5 mm initial filtering. However, if 5 mm HYPRpost-processing is used then HYPR-AU-OSEM always offers the smallest CST, followed by theOSEM reconstructions ordered from smallest to largest kernel sizes in the initial filtering.These two protocols can be expected to have the smallest CST: OSEM without any initialfiltering introduces less spatial correlation than OSEM with any amount of initial filtering, andHYPR-AU-OSEM performs filtering within the reconstruction but preserves high frequency fea-tures through the HYPR operator, meaning cluster boundaries—including those of false positiveclusters—are not extended like in standard filtering.Additionally, the fact that the increase in CST for HYPR-AU-OSEM reconstructions movingfrom 2.5 mm to 5 mm HYPR post-processing is not as drastic as in OSEM reconstructions suggeststhat HYPR-AU-OSEM is less sensitive to the kernel size used in HYPR post-processing.The effect of dose levelAn interesting note is that the CST does not monotonically decrease with increasing dose for allprotocols, as one may expect. This effect is especially present when uniform weighting is usedand/or when a larger kernel size is used for any type of post-processing. This indicates that evenafter 50 independent noisy realizations for each dose level, some dose levels had more cluster sizeslocated in the long tail of the distribution of all theoretically possible cluster sizes than others.Therefore, when the false positive thresholding was performed to remove 90% of all false positiveclusters, the resulting CST may be underestimated for the dose levels that have less populatedtails.Notably, it seems that the 50 noisy realizations for 20 mCi seem to contain fewer false positivecluster sizes in the tail; for most protocols and weighting schemes we see a drop in size from 16mCi to 20 mCi, but then a slight increase from 20 mCi to 24 mCi. This effect is highlighted inFigure 4.1, where distributions of false positive cluster sizes for OSEM-5mm-5mm (uncorrectedweighting) at 20 mCi and 24 mCi are compared, showing the tail of the 20 mCi distribution to beless populated than the 24 mCi distribution.584.1. Baseline simulationFigure 4.1: Comparing false positive cluster size distributions in OSEM-5mm-5mm: 20 mCi vs 24 mCi.Solid lines indicate the CST for the respective dose levels.Overall, if the noise pattern is kept fixed and we increase the dose level (i.e. decrease the noiselevel), then the CST should decrease. However, in a real study one will also have the effect ofdifferent weight and metabolism of the participants, as well as different scanner sensitivity if theHRRT is not used. For a fixed dose level, increasing body mass will result in lower concentrationvalues and hence lower SNR. This variability, coupled with the difficulty of properly sampling thetheoretical (unknown) distribution of possible CSTs with only 50 realizations, can be dealt with byerring on the side of caution and selecting the CST for the analysis of the other simulations (and inreal human data) to be the largest cluster size found across all doses. These values are presentedin Table 4.4.Table 4.4: Cluster size thresholds used for analyzing random release magnitude and cluster release simula-tions, as well as real human data.Weighting traditional uniform uncorrected[voxels] [cm3] [voxels] [cm3] [voxels] [cm3]OSEM-NF-2.5mm 42 0.076 54 0.098 25 0.045OSEM-2mm-2.5mm 59 0.107 97 0.176 39 0.071OSEM-5mm-2.5mm 239 0.433 302 0.547 153 0.277HYPR-2.5mm-2.5mm 63 0.114 94 0.170 42 0.076OSEM-NF-5mm 178 0.322 270 0.489 119 0.215OSEM-2mm-5mm 202 0.366 294 0.532 141 0.255OSEM-5mm-5mm 370 0.670 473 0.856 235 0.425HYPR-2.5mm-5mm 141 0.255 218 0.395 108 0.1964.1.2 Comparing results using previous noise-modeling methodsThese values imply that using the best weighting scheme, uncorrected weighting, we can theoreti-cally detect10 structures as small as 25 voxels (0.045 cm3) if using 2.5 mm HYPR post-processing(OSEM-NF-2.5mm), or 108 voxels (0.196 cm3) if using 5 mm HYPR post-processing (HYPR-2.5mm-5mm).10The CST defines the smallest cluster size that, if detected, can be assumed to be a true positive with 90%confidence. Of course, true release clusters of this size are not guaranteed to be detected, since, as we shall show,detection sensitivity various as a function of release level.594.2. Random release magnitude simulationPrevious work, using a 6 mm by 6 mm by 6 mm boxcar kernel (2 mm voxels) for HYPR post-processing found that with a 90 minute 20 mCi bolus scan, one can detect structures as small as124 HRRT-equivalent voxels (0.224 cm3) [60]. However, these simulations only added Gaussiannoise scaled to levels encountered in typical OSEM-reconstructed HRRT image to the noiselessTACs, rather than adding Poisson noise in sinogram space followed by reconstruction. The closestdenoising protocol in our work would be then be OSEM-NF-5mm, having a CST ranging from 119to 270 voxels depending on the weighting scheme used.Additionally, we reproduce their noise addition methods to our noiseless, PSF-filtered baselinesimulations, and subsequently denoise with HYPR using a 5 mm Gaussian kernel. In these simu-lations, find a CST of 103 voxels (traditional weighting), 138 voxels (uniform weighting), and 86voxels (uncorrected weighting) after 50 noisy realizations.Since this is lower than our more realistic simulation methods, this suggests that the effect ofnormalization, attenuation, and reconstruction increases the CST. The cause is likely from thereconstruction introducing spatial correlation both in signal and noise along LORs, whereas whenGaussian noise is simply added, the noise sample at every voxel is independent of the others, andno spatial correlation is introduced between voxels until HYPR post-processing is applied.Another explanation is that the method of determining noise level (i.e. setting the parameterµ) may not be accurate. In this method, a voxelwise coefficient of variation (COV) was taken in thestriatum over the last 15 mins of a bolus plus constant infusion (B/I) healthy RAC scan. B/I wasused such that the COV is not biased by changes in the true concentration levels. The voxelwiseCOVs were then histogrammed, and a Gaussian was fit to determine the mean COV, from whichµ can be compute. This method may not be robust, since the voxelwise standard deviations aretaken over a relatively short time frame. As well, built in is the assumption that the temporaland spatial noise characteristics are Gaussian, when the latter may be better approximated by alognormal distribution [80, 81]. This may potentially lead to this method underestimating the noiselevel in the TACs, making false positive removal an easier task.4.2 Random release magnitude simulation4.2.1 Sensitivity as a function of release levelThe main goal of the random release magnitude simulation is to quantify sensitivity as a functionof DA release level for each protocol. Figure 4.2 plots the voxelwise sensitivity as a function ofDAR for all protocols (and all weighting schemes) at a dose level of 16 mCi.The global sensitivity, defined as the average sensitivity over all release levels, is presented inTables 4.5-4.7 for all weighting schemes at dose levels of 16 mCi, 20 mCi, and 24 mCi, respectively.Initial observations from these data show a sublinear trend in general; the increases in sensitivitywith increasing DAR are steady at low to moderate DAR, but eventually become less significant asthe sensitivity approaches 100% at higherDAR. DAR greater than those tested in these simulationswill likely consistently produce sensitivities of ∼ 100%. However, such high release levels are likelynot to be expected in response to neurophysiological tasks or stimuli, such as a gambling task orcigarette smoking [57, 69–71].The data also make it clear that the lack of a 5 mm kernel in the post-processing step—either asa post-filter or in the HYPR operator—is detrimental to voxelwise sensitivity. This means OSEM-NF-2.5mm, OSEM-2mm-2.5mm, and HYPR-2.5mm-2.5mm can immediately be disregarded in thesearch for an optimal denoising protocol.For the remaining protocols, OSEM-5mm-5mm produces the highest sensitivity by a significantmargin. OSEM-NF-5mm, OSEM-2mm-5mm, HYPR-2.5mm-5mm and OSEM-5mm-2.5mm tend to604.2. Random release magnitude simulationTable 4.5: Global sensitivity in the 16 mCi random release magnitude simulations across all release levelsfor all protocols and weighting schemes. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrectedOSEM-NF-2.5mm 12.0± 0.5% 10.9± 0.6% 7.3± 0.3%OSEM-2mm-2.5mm 24.5± 0.7% 21.6± 0.8% 17.8± 0.6%OSEM-5mm-2.5mm 63.8±1.0% 60.0±1.0% 57.6±1.1%HYPR-2.5mm-2.5mm 26.2± 0.8% 23.4± 0.9% 20.3± 0.8%OSEM-NF-5mm 55.6± 1.0% 51.6± 1.1% 49.1± 1.1%OSEM-2mm-5mm 61.3± 1.0% 57.1± 1.1% 54.8± 1.1%OSEM-5mm-5mm 76.5±0.8% 73.6±0.9% 71.7±1.0%HYPR-2.5mm-5mm 58.5± 1.1% 54.0± 1.2% 52.1± 1.2%Table 4.6: Global sensitivity in the 20 mCi random release magnitude simulations across all release levelsfor all protocols and weighting schemes. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrectedOSEM-NF-2.5mm 14.9± 0.5% 13.6± 0.5% 10.0± 0.4%OSEM-2mm-2.5mm 29.3± 0.7% 26.1± 0.7% 22.7± 0.6%OSEM-5mm-2.5mm 72.1±0.8% 68.7±0.9% 66.7±0.9%HYPR-2.5mm-2.5mm 32.1± 0.8% 28.7± 0.8% 25.6± 0.7%OSEM-NF-5mm 64.6± 0.8% 60.5± 0.9% 58.4± 0.9%OSEM-2mm-5mm 69.7± 0.8% 66.1± 0.9% 64.0± 0.9%OSEM-5mm-5mm 82.7±0.7% 80.4±0.8% 78.9±0.8%HYPR-2.5mm-5mm 67.3± 0.9% 63.6± 1.0% 61.0± 1.0%Table 4.7: Global sensitivity in the 24 mCi random release magnitude simulations across all release levelsfor all protocols and weighting schemes. The results are subdivided by weighting scheme and HYPR post-processing kernel size; boldface indicates the best performing denoising protocol within each subdivision.Weighting traditional uniform uncorrectedOSEM-NF-2.5mm 17.0± 0.5% 15.5± 0.5% 11.9± 0.4%OSEM-2mm-2.5mm 32.3± 0.6% 29.1± 0.6% 25.8± 0.6%OSEM-5mm-2.5mm 75.5±0.7% 72.2±0.8% 70.7±0.8%HYPR-2.5mm-2.5mm 35.0± 0.8% 32.4± 0.7% 28.1± 0.7%OSEM-NF-5mm 68.5± 0.8% 64.5± 0.9% 62.8± 0.9%OSEM-2mm-5mm 73.2± 0.8% 69.8± 0.9% 68.3± 0.8%OSEM-5mm-5mm 84.9±0.6% 82.7± 0.7% 81.6± 0.7%HYPR-2.5mm-5mm 71.1± 0.8% 67.8± 0.9% 65.3± 0.9%614.2. Random release magnitude simulation(a) traditional weighting (b) uniform weighting(c) uncorrected weightingFigure 4.2: Voxelwise sensitivity as a function of DAR for the three weighting schemes investigated (16mCi dose level).perform similarly across all release levels, producing better-than-chance voxelwise sensitivity (i.e.> 50%) above ∼ 12% DAR.In the next sections, we will provide a more detailed comparison of the methods across thevarious simulation and analysis variables: weighting scheme, HYPR post-processing kernel size,reconstruction type and post-filter size, and dose level.To avoid confusing the reader with cluttered plots, we draw from the observation that OSEM-5mm-5mm evidently produces the highest voxelwise sensitivity and the remaining protocols using a5 mm kernel in some denoising operation all perform similarly; in what follows we thus only presentcomparisons of OSEM-5mm-5mm and HYPR-2.5mm-5mm in subsequent plots (the conclusions forOSEM-NF-5mm, OSEM-2mm-5mm, and OSEM-5mm-2.5mm are sufficiently similar to HYPR-2.5mm-5mm).624.2. Random release magnitude simulationFigure 4.3: Comparing the effect of weighting scheme used on voxelwise sensitivity as a function of DAR.The effect of weighting schemeIn Figure 4.3 we plot sensitivity as a function of release level for OSEM-5mm-5mm and HYPR-2.5mm-5mm for each weighting scheme (16 mCi dose level only; the comparisons from other doselevels are comparable).From these simulations, the best performance in terms of voxelwise sensitivity comes fromusing the traditional weights, followed by uniform weights and then uncorrected weights. Sincethe decay correction is larger at later frames, the fitting algorithm “cares” less about these points.Therefore when computing an F-value at a voxel, perturbations in the post-intervention periodindicating potential DA release are not given as much weight in the weighted residuals as comparedto uniform or traditional weighting, meaning a larger perturbation (i.e. higher release level) isneeded in the uncorrected weighting scheme to produce a significant F-value.As evident from the previous section, this fact is useful for rejecting false positives from post-intervention perturbations due to noise—resulting in a lower CST—but this comes with the dis-advantage of being less sensitive to post-intervention perturbations due to a true drop in the RACsignal.Conversely, in uniform or traditional weighting, we are more sensitive to post-interventionperturbations, so we will have more false positive detection and hence a higher CST, but alsohigher sensitivity—at least if the region detected is larger than the CST.This latter point needs to be emphasized: in this particular simulation, entire regions of thestriatum contain DA release of varying release level, meaning the true clusters are all over 1000voxels in size—much larger than the CST for any protocol. Therefore we do not need to worry asmuch about the CST potentially limiting sensitivity, as most of the clusters detected are also muchlarger than the cluster size threshold for the adequately performing methods. However in a casewhere smaller true clusters exist, such as the cluster release simulations, this may not be the caseas we shall investigate in Section 4.3.Combining these results with the CST data from the previous section, we can conclude thatuniform weighting neither produces the optimal sensitivity nor CST. Therefore, for all futureanalysis we will only analyze using traditional weighting (for optimal sensitivity) and uncorrectedweighting (for optimal CST).634.2. Random release magnitude simulationThe effect of HYPR post-processing kernel sizeAs mentioned earlier, these results indicate that all protocols with 2.5 mm HYPR post-processingproduce unacceptably low sensitivities across all release levels, with the exception of OSEM-5mm-2.5mm. In other words, a 5 mm filtering operation—either initially or within the HYPR operator—needs to be performed in order to produce acceptable sensitivity.Spatial filtering introduces spatial correlation in neighbouring voxels within true release clusters,such that the true signal containing the DA-induced perturbation adds coherently and the noiseadds decoherently, thus increasing SNR. With narrower spatial filtering kernels, the data remainsufficiently noisy such that the ability to discern a DA-induced perturbation is difficult, and as aresult only small fractions of true release clusters are reliably detected, as evident from the lowvoxelwise sensitivity when using 2.5 mm.The effect of reconstruction method and initial filteringHolding the weighting scheme and dose level fixed, we see that OSEM with post-filtering, OSEMwith a 2 mm post-filter, and HYPR-AU-OSEM all produce similar sensitivity across all releaselevels. OSEM with a 5 mm filter produces significantly higher sensitivity than the others across allrelease levels, but separation in sensitivity narrows as one moves to higher release levels.This can be explained by the wider spatial kernel correlating true low release voxels with anynearby higher release voxels, thus increasing their probability of being detected. Of course in a truerelease cluster that only contains low release voxels (such as the low release clusters of the clusterrelease simulation), then this effect should disappear.Keeping the CST vs sensitivity trade-off in mind, OSEM-5mm-5mm fitted with traditionalweights would be the ideal protocol if one is not concerned about detecting smaller regions of DArelease (i.e. smaller than the OSEM-5mm-5mm and OSEM-5mm-2.5mm CSTs): it has a ∼ 16%greater global sensitivity than the next best protocol.However, if one values detecting smaller regions of release, then HYPR-2.5mm-5mm with uncor-rected weighting would be the ideal protocol: the 23% reduction in global sensitivity—as comparedto OSEM-5mm-5mm with traditional weighting—is offset by a 71% reduction in CST. Keep inmind that if one also values parametric accuracy, then they will have to wait until Section 4.2.2 toobtain the complete picture.The effect of dose levelIn Figure 4.4 we compare sensitivity as a function of dose level for OSEM-5mm-5mm and HYPR-2.5mm-5mm (uncorrected weighting only; the comparisons for traditional weighting are compara-ble). Predictably, an increased dose results in increased sensitivity. Averaging across release levels,moving from 16 mCi to 20 mCi the sensitivity is increased by 9% for OSEM-5mm-5mm and 17%for HYPR-2.5mm-5mm. Moving from 20 mCi to 24 mCi, the improvements are 3% for OSEM-5mm-5mm and 6% for HYPR-2.5mm-5mm. This suggests that we receive diminishing rewards forincreasing the dose level much more above 20 mCi, especially when weighing these benefits againstthe increased radiation dose delivered to the subject.4.2.2 Parametric accuracyEffects of post-processing and incorrect basis functionsThere are a number of parameters that are determined in the lp-ntPET fitting process. In mostscans, the parameters of interest will be the ones that describe the dynamics of the DA release:644.2. Random release magnitude simulationFigure 4.4: Comparing the effect of dose level on sensitivity as a function of DAR.tD and tP . The parameters α and γ are less important on their own, but become more impor-tant as they enter the formula for time-varying binding potential or DAR(t) alongside tD and tP(equations 2.73 and 2.74 respectively).We know that the PSF filtering changes the concentration values in the images, so we haveaccounted for this by defining our ground truth DAR(t) as the results from refitting the noiseless,PSF-filtered images while specifying the correct basis function (i.e. tD = 35 min, tP = 40 min,and α = 2). When noise is added to the simulations, after reconstruction we expect the voxelwiseestimated parameters to be somewhat biased, but the average over a relatively large uniform regionin the parametric images should be quite accurate. Thus we apply our denoising protocols, which inprinciple will improve the voxelwise estimates, but likely will introduce some bias into the parameterestimates.A factor that will also introduce bias comes from the finite number of basis functions specifiedin the lp-ntPET fitting procedure. In a real scan, we will not know a priori the dynamics of thetrue DA release. By specifying only 819 basis functions, in fact it is exceedingly unlikely thatwe will specify the exact parameters that represent the ground truth dynamics. Therefore in ourbasis functions we have intentionally not supplied the correct values of tD, tP , and α. The closestestimates of these parameters from the supplied basis functions are tˆD = 34.5 min, tˆP = 40.5, andαˆ = 1.To isolate the significance of these two effects, we applied the fitting procedure to the noiselesssimulations in two ways. First, the noiseless ground truth images are post-processed and fit withthe correct basis function provided, so as to only assess the effect of post-processing on parametricaccuracy. Second, we fit the noiseless ground truth images with the standard set of 819 basisfunctions, to assess the effect of not having the correct basis function available.In the random release magnitude simulation, since we have a distribution of DAR over the entirestriatal region except for the R putamen, it is impossible to compare all voxel-level DAR(t) curves.While all their timing parameters are constant, the magnitude of the release level is different atevery voxel. Therefore we choose to perform this analysis separately on two narrow representativeintervals of DAR: 7% < DAR < 8% for low release, and 15% < DAR < 16% for high release.Figure 4.5 shows the mean bias in the DAR(t) curves computed from the ground truth imagesand the ground truth images after various post-processing [i.e. post-processed DAR(t) subtractedby ground truth DAR(t)], when the correct basis function is supplied.Since, in this simulation, high release voxels tend to be surrounded by lower release voxels,654.2. Random release magnitude simulation(a) low release (b) high releaseFigure 4.5: The bias in DAR(t) introduced by post-processing, in the noiseless random release magnitudesimulation.(a) low release (b) high releaseFigure 4.6: Demonstrating bias introduced into the TACs by post-processing in the noiseless randomrelease magnitude simulation.664.2. Random release magnitude simulation(a) low release (b) high releaseFigure 4.7: Demonstrating bias introduced by not having the correct basis function for fitting, in thenoiseless random release magnitude simulation.we expect the high release DAR(t) to be negatively biased by post-processing due to the spatialcorrelation introduced by the filters, which is what we observe. However for low release DAR(t),though we would expect a positive bias, we actually only observe this for when a 5 mm post-filterfollowed by 2.5 mm HYPR post-processing is applied.5 mm HYPR post-processing alone is negatively biased, which suggests that HYPR post-processing has the effect of negatively biasing any release level, whereas spatial post-filtering willnegatively bias higher release DAR(t) and positively bias lower release DAR(t). As a result, whenspatial post-filtering plus HYPR post-processing with a relatively wide kernel (5 mm) is applied, thetwo effects appear to cancel out in low release DAR(t). When a narrower HYPR post-processingkernel is used instead (2.5 mm), the positive bias from spatial filtering dominates and hence weobserve a net positive bias in low release DAR(t).In high release DAR(t), since both HYPR post-processing and spatial post-filtering both intro-duce negative biases, we always observe a net negative bias. However, when the narrower HYPRpost-processing kernel is used, the negative bias is not as severe.The corresponding mean TACs for the high release and low release voxels are shown in Fig-ure 4.6. Especially in the high release TACs, we can see how HYPR post-processing slightly distortsthe temporal pattern such that the drop in the TAC from DA release is shallowed out, whereasspatial post-filtering has the effect of a relatively constant negative bias across all time points.In the HYPR operator, each individual filtered frame is multiplied by the ratio of the unfilteredand filtered composite image. The composite image is a sum of the uncorrected frames, and hencecontains more signal from the higher count frames at the beginning of the scan. During theseframes there is no release occurring, and so the contrast between regions where release is occurringand not occurring is not as strong in the individual frames where the release is peaking. Therefore,upon multiplying these voxels by the ratio of the unfiltered and filtered composite, these regionsbecome slightly negatively biased. Furthermore, because this ratio extracts the high frequencyfeatures from the composite, any bias from the spatial filters first applied to the frames in theHYPR operator is mitigated, such that the net effect is always a negative bias in DAR(t).674.2. Random release magnitude simulationThese biases can be expected for a given post-processing protocol, but only if the correct basisfunction is selected by the fitting algorithm. If we instead fit the noiseless data without any post-processing, but with the standard set of 819 basis functions (which does not include the correctbasis function), then from Figure 4.7 we also observe a subtle bias in the resulting fitted DAR(t).Here we observe tD positively biased and tP negatively biased for both high and low release. Tomaintain roughly the same area under the curve, γ is then positively biased to increase the peakmagnitude. The biases in the timing parameters likely result from not having the ground truthα = 2 and instead fitting α = 1; in the gamma variate function the term (tP − tD)/α—controllingthe decay of the curve when the exponential portion dominates—will increase from the groundtruth when fitting α = 1, so the fitting algorithm compensates by decreasing the difference tP − tD.When fitting to the noisy, post-processed datasets, we can expect a mixture of these parametricbiases to occur. These results are presented in the next few sections.Timing parametersWith the effect of post-processing on parametric accuracy noiseless data in mind, let us now proceedto assess how parametric accuracy is affected in the noisy simulations. In the noiseless case, theestimated timing parameters had zero variance. Evidently, this will not be the case when noiseis present in the voxel-level TACs, and so a closer look at the individual timing parameters iswarranted before we study the DAR(t) curves.The mean fitted tD and tP as a function of DAR are shown for the five best denoising protocolsin the 16 mCi simulations, using traditional and uncorrected weighting, in Figure 4.8.In general, the mean of the fitted parameters are within 30 seconds for all protocols and weight-ing schemes, except for estimates at the lowest DAR when using traditional weighting, where theestimates are within one minute.At moderate to high DAR, the mean estimates are slightly negatively biased for tP and slightlypositively biased for tD. This reflects our earlier observation that the fit is likely compensating(here on average) for the unavailability of α = 2 by instead choosing αˆ = 1 with a smaller tP − tDto best match the release dynamics with the given basis functions.Between the different denoising protocols, there is little variation in the estimates, thoughOSEM-5mm-5mm tends to have the lowest estimates at both tD and tP . These differences are lessprevalent when using uncorrected weighting versus traditional weighting.These results suggest that the mean voxel-level timing parameter estimates (over the entiredetected cluster) are very accurate—provided the conditions of the random release magnitudesimulation are valid: large extended regions of DA release exist with spatially-invariant timing, butspatially-variant release magnitude.To study the voxel-level accuracy, we compute MAE(tD) and MAE(tP ) as a function of releaselevel for the five best denoising protocols, again comparing traditional and uncorrected weighting(16 mCi dose level), and plot them in Figure 4.9.As one may expect, as DAR increases we observe a decrease in MAE(tD) and MAE(tP ). How-ever, MAE(tD) is always less than MAE(tP ). This is likely because the onset of the release is moreidentifiable; in the TAC this is event is characterized by a kink due to the sudden increase in thegamma variate function describing the release, whereas the peak release level is a more smoothevent in the TAC.Nevertheless, depending on the denoising method used, on average we can localize voxel-level tDto within ∼ 3 min (a single frame duration) and tP to about ∼ 5 min for the highest release levels.For the lowest release levels, these bounds increase to tD within ∼ 4 min and tP within ∼ 6 min.Since different framing protocols were not investigated in this work, it is unknown whether using684.2. Random release magnitude simulation(a) traditional weighting (b) uncorrected weightingFigure 4.8: Comparing the mean fitted timing parameters between the best-performing denoising protocol,for traditional and uncorrected weighting schemes. Both plots use the same colour scale. Solid lines indicatethe ground truth tD = 35 min and tP = 40 min.(a) traditional weighting (b) uncorrected weightingFigure 4.9: Comparing the mean absolute error in tD and tP between the various denoising protocols, fortraditional and uncorrected weighting schemes.shorter frames would increase this precision. We can speculate though that a potential increasein precision from more time points available about the true tD and tP may be offset by decreasedprecision from increased temporal noise resulting from less counts per frame.OSEM-5mm-5mm consistently yields lower MAE in the timing parameters than the other de-noising protocols, though at high DAR the differences in MAE(tP ) vanish. The wide spatial kernelsof OSEM-5mm-5mm homogenize the data, such that there is less variance in the timing parametersestimates. For instance, at DAR = 10%, the standard deviation of voxel-level estimates when usingtraditional weighting is 3.6 min (3.8 min) for tD (tP ), whereas the other protocols average to 4.1694.2. Random release magnitude simulationmin (4.7 min). Because the mean estimates are fairly accurate, the lower variance subsequentlyresults in lower MAE as well.Another caveat with our analysis here is that our simulations have identical timing parametersat each voxel within a cluster. If a gradient in timing parameters were to exist within a cluster,then it is likely that the MAE of the timing parameters would increase, since the post-processingfilters would correlate neighbouring voxels that in reality have different release dynamics. However,with each voxel on the order of a cubic millimetre, we do not expect local variations in timing tobe a realistic scenario.When comparing traditional versus uncorrected weighting, the former produces slightly higherMAE(tD), and slightly lower MAE(tP ), though these differences are small; the more importantfactor here is the denoising protocol used.As was seen before, as dose level is increased the sensitivity improves—especially at lower releaselevels. Although more voxels are accepted, these new voxels may have poorer parametric accuracy,effectively canceling or overpowering the potential increasing in timing accuracy from less noisyTACs. This is observed in Figure 4.10, where we compare MAE in the timing parameters whenusing OSEM-5mm-5mm or HYPR-2.5mm-5mm, at all three dose level tested: MAE is actuallyhigher at 24 mCi as compared to 20 mCi at some release levels. The cluster-level mean estimatesdo not change appreciably with dose level, and so we do not present those results in this work.Figure 4.10: Comparing the effect of dose level on the mean absolute error in the timing parameters as afunction of DAR.BP (t) and DAR(t)With an understanding on how all factors contribute to timing parameter accuracy, we can nowobtain a more concrete picture of the bias in the estimated temporal dynamics by analyzing theBP (t) and DAR(t) curves.Just as in investigating the effect of post-processing on the noiseless data, we chose to comparethe BP (t) and DAR(t) curves on the same two narrow intervals of release level: 7% < DAR < 8%and 15% < DAR < 16%.Since outlier fits are to be expected at the voxel-level, we choose to compute the median valuesof BP (t) and DAR(t) at every time point (90 second temporal resolution). To form an idea of thevoxel-level variance in these estimates, we plot the median ± one median absolute deviation (MAD).The results for BP (t) are presented in Figure 4.11 for uncorrected weighting (16 mCi simula-tion); the effect of dose and weighting scheme are negligible when it comes to these estimates, thuswe do not present them in this work. The corresponding DAR(t) are in Figure 4.12.704.2. Random release magnitude simulation(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.11: Fitted BP (t) as compared to the ground truth curves for the five best denoising protocols.BP (t) is presented at both low release (7% < DAR < 8%) and high release (15% < DAR < 16%). Colouredlines indicate the median value at each time point (90 second temporal resolution) and shaded regionsindicate ± one median absolute deviation. Black lines and shading indicate the ground truth. A spread inthe ground truth exists since the subregions of the striatum have different ground truth BP (t).(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.12: Fitted DAR(t) as compared to the ground truth curves for the five best denoising protocols.DAR(t) is presented at both low release (7% < DAR < 8%) and high release (15% < DAR < 16%).Coloured lines indicate the median value at each time point (90 second temporal resolution) and shadedregions indicate ± one median absolute deviation. Black lines and shading indicate the ground truth. Aspread in the ground truth exists because the range of DAR used is not infinitesimally small.From visual inspection, OSEM-5mm-5mm provides the highest precision (i.e. lowest MAD) atlow release, whereas OSEM-NF-5mm provides the poorest precision. At higher release, however,OSEM-5mm-5mm provides similar precision and the worst parametric accuracy. This is especially714.2. Random release magnitude simulationclear for BP (t), but less so for DAR(t); both the pre-intervention and post-intervention valuesof BP (t) are severely negatively biased, so they partially cancel out when computing DAR(t)[equation 2.74)].The severe negative bias in BP (t)—especially during the tail end of the post-interventionperiod—stems from the wide spatial kernels of OSEM-5mm-5mm negatively biasing the TACs,yielding lower fitted binding potential values. The remaining protocols all produce slight positivebiases. Since all protocols use HYPR post-processing, it is likely that these biases are from theHYPR operator biasing the shape of the TAC, since the overall bias in the magnitude of the TACsis much smaller than that induced by a spatial filtering alone.Comparing the estimates of DAR(t) at high and low release levels, the low release DAR(t)estimates are consistently overestimated and conversely the high release DAR(t) estimates areconsistently underestimated. These observations agree with our analysis of post-processing inthe noiseless case for high release: all denoising methods produced a negatively biased DAR(t)estimate. However for low release, we discovered only OSEM-5mm-2.5mm produced a positive biasin DAR(t); HYPR post-processing tends to negatively bias DAR(t) at all release levels, but spatialfilters will positively bias DAR(t) at low release, as these voxels are surrounded by higher releasevoxels on average.Thus there is likely a selection effect being observed here, where low release voxels are moredifficult to detect (i.e. have lower detection sensitivity). The set of low release voxels will be skewedby those that—due to noisy perturbations or other biases in their TACs—appear to have higherrelease than their ground truth values.As detection sensitivity increases, however, this effect should be less prominent—exactly whatwe observe. The highest sensitivity method, OSEM-5mm-5mm, provides the least positively biasedDAR(t) estimate at low release, whereas the lowest sensitivity method, OSEM-NF-5mm, providesthe most positively biased estimate.As a function of dose level, the median estimates of BP (t) and DAR(t) do not change ap-preciably, and so we do not present these results here. However, a small reduction in MAD isobserved—indicating that the voxel-level estimates are slightly more consistent from the less noisyvoxel-level TACs at higher dose.4.2.3 Validity of the cluster size thresholdWe intentionally did not simulate release in the R putamen of the random release magnitudesimulation so as to check whether the CST is functioning correctly. In Table 4.8, we tabulatethe percentage of random release magnitude simulations that are completely free of false positiveclusters.All protocols at all dose levels remove all false positive clusters from at least 90% of the simu-lations. In the baseline simulations used to determine the CST, all four subregions of the striatumdo not contain release, whereas in the random release magnitude simulation only the R putamendoes not contain release. Thus, statistically it is less likely for a false positive to exist, and so it isnot surprising that the CST is functioning at its intended false positive-free rate or better.724.2. Random release magnitude simulationTable 4.8: Assessing the functionality of the CST by computing the percentage of random release magnitudesimulations that are completely free of false positives in the R putamen of the random release magnitudesimulation. If functioning correctly, the CST should statistically remove all false positive clusters from 90%of simulations.False positive-free rate [%]Traditional weighting Uncorrected weightingMethod 16mCi 20mCi 24mCi 16mCi 20mCi 24mCiOSEM-NF-5mm 98 94 100 96 90 98OSEM-2mm-5mm 96 94 100 96 92 98OSEM-5mm-5mm 92 98 100 94 92 98HYPR-2.5mm-5mm 96 92 100 98 94 100OSEM-5mm-5mm 96 96 100 96 92 984.2.4 Comparing results using previous noise-modeling methodsThe Gaussian noise addition method of previous work [57, 60] was applied to the random re-lease magnitude simulation (16 mCi), followed by HYPR post-processing with a 5 mm kernel. InFigure 4.13, we plot the voxelwise sensitivity as a function of DAR for this method (using both tra-ditional and uncorrected weighting). The curves from HYPR-2.5mm-5mm and OSEM-5mm-5mmare also shown for reference.(a) traditional weighting (b) uncorrected weightingFigure 4.13: Comparing the voxelwise sensitivity as a function of DAR between the Gaussian noise additionmethod and the more realistic noise addition method used in this work.When Gaussian noise addition is used, the resulting voxelwise sensitivity is much higher thaneven the best performing method presented in this work thus far (OSEM-5mm-5mm).Our noise modeling methods add Poisson noise in sinogram space while also accounting fordetector normalization and attenuation, and therefore are a better approximation of the imagingprocess. The sum of these effects evidently negatively affect detection sensitivity. This may bethe result of the noise becoming spatially correlated after reconstruction in our simulation meth-ods, whereas the noise patterns are spatially-invariant when simply adding Gaussian noise to the734.2. Random release magnitude simulationFigure 4.14: Comparing the mean timing parameter estimates from the Gaussian noise simulations withthose from OSEM-5mm-5mm and HYPR-2.5mm-5mm. Ground truth values are indicated by solid blacklines.Figure 4.15: Comparing MAE(tD) and MAE(tP ) from the Gaussian noise simulations with those fromOSEM-5mm-5mm and HYPR-2.5mm-5mm.noiseless image. If a certain noise pattern makes it more difficult to detect the drop in RAC signalfrom DA release, then this pattern may be spread locally during reconstruction. The region withthis pattern is then preserved by HYPR post-processing, ultimately resulting in a region of falsenegatives.In contrast, it is very unlikely that such an extended spatial noise pattern will result fromsuccessive random samples from a Gaussian noise distribution, and thus noise patterns that nega-tively affect sensitivity are more randomly distributed amongst the voxels. Non-extended temporalpatterns are more difficult to preserve through HYPR post-processing (since the HYPR kernel sizeis several voxels in size), and thus on average such noise patterns will be neither extended norpreserved, ultimately resulting in higher sensitivity.In terms of parametric accuracy, the mean timing parameter estimates when using Gaussiannoise addition are generally similar to those from the simulations with our noise modeling methods,especially as one moves to higher release levels (c.f. Figure 4.14). This tells us that, on the cluster-level, the effect of which noise modeling method is used is less important.However, when considering voxel-level timing accuracy, MAE(tD) is significantly lower on av-744.2. Random release magnitude simulation(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.16: Fitted BP (t) as compared to the ground truth curves for the Gaussian noise addition method.The curves from OSEM-5mm-5mm and HYPR-2.5mm-5mm are also shown for reference. BP (t) is presentedat both low release (7% < DAR < 8%) and high release (15% < DAR < 16%). Coloured lines indicate themedian value at each time point (90 second temporal resolution) and shaded regions indicate ± one medianabsolute deviation. Black lines and shading indicate the ground truth. A spread in the ground truth existssince the subregions of the striatum have different ground truth BP (t).(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.17: Fitted DAR(t) as compared to the ground truth curves for the Gaussian noise additionmethod. The curves from OSEM-5mm-5mm and HYPR-2.5mm-5mm are also shown for reference. DAR(t)is presented at both low release (7% < DAR < 8%) and high release (15% < DAR < 16%). Coloured linesindicate the median value at each time point (90 second temporal resolution) and shaded regions indicate ±one median absolute deviation. Black lines and shading indicate the ground truth. A spread in the groundtruth exists because the range of DAR used is not infinitesimally small.erage when Gaussian noise addition is used (c.f. Figure 4.15). As well, MAE(tP ) is significantlyreduced at the lowest DAR, but this reduction is smaller with increasing DAR, until no reductionis observed at the highest DAR.These reductions may occur because the effects of normalization, attenuation, and reconstruc-tion add bias from the ground truth; Gaussian noise will only add variance about the ground truthbut preserve the mean. Therefore even though the noise levels may be matched, the error aboutthe ground truth concentration will be higher in our simulation method and thus presents theopportunity for greater voxel-level errors in timing parameters on average.This hypothesis is further supported when observing the median BP (t) and DAR(t) in Fig-ure 4.16 and Figure 4.17, respectively, where the median estimates resulting from the Gaussiannoise addition method are similar to those of HYPR-2.5mm-5mm, albeit with lower MAD. Theone exception is the lower release estimate of DAR(t), where the Gaussian noise addition methodproduces a less positively biased curve than HYPR-2.5mm-5mm. However, again this is likely aselection effect with HYPR-2.5mm-5mm, where voxels having greater apparent peak DAR(t) are754.2. Random release magnitude simulationpreferentially selected. This selection effect presumably disappears in the Gaussian noise additionmethod, since the sensitivity is much higher.Alternatively, there is also the explanation that this method of noise addition underestimatesthe true noise levels in the image, making voxelwise detection of DA release an easier task.4.2.5 ConclusionsThe random release magnitude simulation, while not perfect, provides an estimation of the effectof denoising on sensitivity as a function of release level. Since the regions of release are large, theeffect of CST is likely small or even negligible in determining sensitivity. Therefore the comparisonof sensitivity at various release levels quoted in this section are most valid when large clusters ofrelease at these release levels are present.Under these conditions, the results from these simulations reveal that OSEM-5mm-5mm is theideal denoising method to achieve optimal sensitivity. However, as release level becomes higher,the differences in sensitivity between OSEM-5mm-5mm and the other methods are not as large. Atmuch higher release levels (i.e. beyond those tested in this work to simulate a realistic task-inducedrelease), the difference can be expected to vanish altogether, such that any denoising methodpresented here will suffice.The effect of weighting scheme used on sensitivity is small, though traditional weighting pro-duces slightly higher sensitivity than uniform or uncorrected weighting. From Section 4.1.1, uniformweighting was shown to require a much larger CST than the other schemes. However, because CSTis not an issue when cluster size is large, we can conclude in this scenario that the choice of weightsis not as important as the choice of denoising method.Increasing dose in these simulations show a modest increase in sensitivity, though the benefitgained from higher doses is smaller as one continues to increase dose. As well, the ethical andeconomical concerns will likely outweigh slightly higher sensitivities achieved through increasingdose beyond about 20 mCi.In terms of parametric accuracy, provided the large region of release has similar release dynam-ics, OSEM-5mm-5mm estimates tD and tP slightly better than the other methods (to ∼ 3 min and∼ 5 min respectively, averaged across release levels). When considering both the magnitude anddynamics of the release together using BP (t) and DAR(t), OSEM-5mm-5mm produces the highestprecision but poorest accuracy at high release levels. Better accuracy is instead supplied by a de-noising method that uses less filtering and/or smaller HYPR kernel sizes, such as OSEM-NF-5mm,HYPR-2.5mm-5mm, or OSEM-5mm-2.5mm.These conclusions would likely be different if the cluster size was approaching the size of theCST for some methods. With smaller clusters, there is the potential that the CST will becomeimportant, meaning some initially-detected voxels after significance thresholding will be rejectedafter cluster size thresholding. As well, the fact that smaller clusters are more effected by spill infrom non-release regions, the parametric accuracy will be more severely affected by larger filteringoperations.This leads us to the next section, where we analyze the cluster release simulations containingmuch smaller clusters of release.764.3. Cluster release simulation4.3 Cluster release simulation4.3.1 SensitivityAs a reminder, the cluster release simulation is composed of six clusters throughout the stria-tum, modulating large/small cluster size and high/low release. The location of the cluster withinthe striatum is also important, as the narrower and smaller caudate is more susceptible to PVE.Table 4.9 characterizes the clusters, which we label C1-C6.Table 4.9: The six clusters included in the cluster release simulation.Cluster Location Size [voxels] Size [cm3] % DARC1 L putamen 500 0.91 15%C2 L putamen 300 0.54 7.5%C3 L caudate 400 0.72 15%C4 R caudate 400 0.72 7.5%C5 R putamen 300 0.54 15%C6 R putamen 500 0.91 7.5%The voxelwise sensitivity is tabulated in Tables 4.10 and 4.11 for all denoising methods at16 mCi, using traditional and uncorrected weighting respectively. The same is done for binary50sensitivity in Tables 4.12 and 4.13, as well as for binary sensitivity in Tables 4.14 and 4.15.Effect of denoising methodOSEM-5mm-5mm generally provides the highest voxelwise and binary50 sensitivity for any dose andweighting scheme, except for some smaller, lower release clusters at lower dose, as will be discussedlater. This means that when any portion of a true release cluster is detected, OSEM-5mm-5mm isthe superior denoising method to use.However, this ignores the fact that sometimes OSEM-5mm-5mm may not detect any release atall. If we pose the question in the opposite direction, asking how often any portion of a true releasecluster is detected (i.e. binary sensitivity), then OSEM-5mm-5mm is always outperformed by othermethods for clusters in the caudate, where PVE is more severe and hence detected portions of theclusters may be smaller than the CST. As well, OSEM-5mm-5mm performs similarly well to theother methods in clusters of high release, where detection is much easier.Therefore, if one wishes to reap the benefits of improved voxelwise sensitivity by employingOSEM-5mm-5mm, they must note that this runs the risk of missing certain clusters completely.Conversely, if one uses a method like HYPR-2.5mm-5mm, they will detect portions of these clusterswith higher probability than OSEM-5mm-5mm (i.e. higher binary sensitivity), but the portion ofthe cluster detected will be smaller on average (i.e. lower voxelwise sensitivity).Effect of cluster characteristicsFor any given denoising method, the general trend observed is that the release level is the mostimportant factor in determining sensitivity. Clusters with DAR = 15% (C1, C3, and C5) all havehigher sensitivity than clusters with DAR = 7.5%, regardless of their size or location within thestriatum.The next most important effect appears to be the size and shape of the anatomical region wherethe cluster is located. The release TACs in the more PVE-affected caudate are contaminated to agreater degree by neighbouring non-release voxels, ultimately making a larger portion of the cluster774.3. Cluster release simulationTable 4.10: Voxelwise sensitivity in each cluster for the 16 mCi simulation, using traditional weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 65.2% 19.7% 49.6% 10.9% 67.7% 15.0%OSEM-2mm-5mm 70.7% 22.9% 54.1% 11.2% 72.4% 18.4%OSEM-5mm-5mm 84.7% 37.2% 55.9% 9.5% 71.7% 27.2%HYPR-2.5mm-5mm 65.6% 21.2% 51.0% 16.4% 68.7% 14.1%OSEM-5mm-2.5mm 73.0% 24.3% 55.2% 10.9% 73.6% 19.0%Table 4.11: Voxelwise sensitivity in each cluster for the 16 mCi simulation, using uncorrected weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 60.3% 18.9% 43.3% 9.9% 63.9% 13.9%OSEM-2mm-5mm 66.5% 22.6% 48.3% 10.4% 68.7% 16.0%OSEM-5mm-5mm 81.8% 34.0% 58.7% 12.6% 82.5% 24.2%HYPR-2.5mm-5mm 60.1% 20.0% 44.6% 10.5% 64.5% 12.5%OSEM-5mm-2.5mm 69.2% 24.8% 50.2% 11.1% 71.2% 17.5%Table 4.12: Binary50 sensitivity in each cluster for the 16 mCi simulation, using traditional weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 90% 18% 52% 4% 94% 6%OSEM-2mm-5mm 96% 28% 70% 8% 98% 6%OSEM-5mm-5mm 100% 44% 78% 12% 82% 26%HYPR-2.5mm-5mm 92% 18% 52% 4% 90% 0%OSEM-5mm-2.5mm 96% 28% 80% 10% 96% 8%Table 4.13: Binary50 sensitivity in each cluster for the 16 mCi simulation, using uncorrected weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 72% 6% 34% 2% 92% 2OSEM-2mm-5mm 90% 18% 56% 2% 94% 2%OSEM-5mm-5mm 100% 42% 76% 8% 100% 20%HYPR-2.5mm-5mm 74% 10% 42% 2% 84% 0%OSEM-5mm-2.5mm 94% 24% 64% 4% 94% 4%784.3. Cluster release simulationTable 4.14: Binary sensitivity in each cluster for the 16 mCi simulation, using traditional weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 100% 54% 96% 28 % 100% 50%OSEM-2mm-5mm 100% 56% 98% 26% 100% 58%OSEM-5mm-5mm 100% 66% 82% 16% 82% 62%HYPR-2.5mm-5mm 100% 58% 98% 48% 100% 52%OSEM-5mm-2.5mm 100% 56% 96% 24% 98% 56%Table 4.15: Binary sensitivity in each cluster for the 16 mCi simulation, using uncorrected weighting.Boldface indicates the highest sensitivity.Method C1 C2 C3 C4 C5 C6OSEM-NF-5mm 100% 52% 96% 32% 100% 54%OSEM-2mm-5mm 100% 54% 98% 30% 100% 54%OSEM-5mm-5mm 100% 64% 96% 26% 100% 56%HYPR-2.5mm-5mm 100% 54% 98% 38% 100% 52%OSEM-5mm-2.5mm 100% 58% 98% 30% 100% 56%more difficult to detect. This can be seen by noting that C5 (R putamen, 300 voxels, 15% release)has a higher sensitivity than C3 (L caudate, 400 voxels, 15% release)—even though they haveequivalent release levels and C3 is larger in size than C5. Along the same vein, C2 (L putamen,300 voxels, 7.5% release) has higher sensitivity than C4 (R caudate, 400 voxels, 7.5% release).Finally, cluster size seems to have a negligible effect unless the cluster size is comparable tothe CST of the denoising protocol used. For instance, when using traditional weights with HYPR-2.5mm-5mm (141 voxel CST), C1 (L caudate, 500 voxels, 15% release) and C5 (R putamen, 300voxels, 15% release) have almost identical sensitivity despite C5 being much smaller in size. How-ever, if OSEM-5mm-5mm is used (370 voxel CST), then the sensitivity is lower for C5 since thecluster will only pass through the CST if most of the true voxels plus some spillover voxels fromthe PSF filtering and/or post-processing are detected (to bring the detected cluster size above 370voxels).Effect of weighting schemeTraditional weighting generally provides a slight increase in voxelwise sensitivity as compared touncorrected weighting, as is shown for OSEM-5mm-5mm and HYPR-2.5mm-5mm in Figure 4.18(16 mCi dose level only; the comparisons from other dose levels are comparable).Notable exceptions are in C3, C4, and C5 for OSEM-5mm-5mm. With traditional weightingincreasing the CST, these clusters are in the caudate or small in size, and thus are more likely tobe rejected by the CST.Similar trends are observed when considering binary50 and binary sensitivity. Here we are moreclearly seeing the effect of CST; especially in C5, using uncorrected instead of traditional weightingwith OSEM-5mm-5mm allows one to more often recover a portion of the cluster, since the CST issmaller.794.3. Cluster release simulation(a) voxelwise sensitivity (b) binary50 sensitivity(c) binary sensitivityFigure 4.18: The effect of weighting scheme used on voxelwise, binary50, and binary sensitivityEffect of dose levelIn most cases, increasing the dose results in higher voxelwise sensitivity, as shown in comparingthe results from OSEM-5mm-5mm and HYPR-2.5mm-5mm across the three doses investigated inFigure 4.19 (uncorrected weighting only; the comparisons for traditional weighting are comparable).Exceptions tend to occur in C4 (R caudate, 400 voxels, 7.5% release) when moving from 16 mCito 20 mCi, where the sensitivity is very low to begin with.This counterintuitive result is likely a statistical aberration; since we only simulated 50 noisyrealizations, outliers can have a significant impact on averaging statistics like voxelwise sensitivity—either a few realizations had much greater sensitivity than average in the 16 mCi simulation, or afew realizations had much smaller sensitivity than average in the 20 mCi simulation.Further support for this claim can be found in the comparisons of binary sensitivity for C4:the 20 mCi simulations have lower binary sensitivity than the 16 mCi simulations, meaning fewer804.3. Cluster release simulation(a) voxelwise sensitivity (b) binary50 sensitivity(c) binary sensitivityFigure 4.19: The effect of dose level on voxelwise, binary50, and binary sensitivity.realizations contribute detected voxels to the voxelwise sensitivity metric.Voxelwise and binary sensitivity between 16 mCi and 24 mCi is almost always higher—mostnotably so in low release clusters. However the gains in sensitivity, on the order of ∼ 20% atmost, may not be worth the 50% increase in dose—a common point harped on in this work whenconsidering the effect of dose.814.3. Cluster release simulation4.3.2 Parametric accuracyEffects of post-processing and incorrect basis functionsLike in the random release magnitude discussion, we also examine the effect of post-processing onthe noiseless cluster release data. Here there are some important differences. Unlike the randomrelease magnitude data, here the clusters have a more homogeneous distribution of DAR values intheir voxels. Therefore the issue of low release voxels being more likely to be surrounded by higherrelease voxels (and vice versa) should not be as pronounced, meaning the filtering kernels of thepost-processing operations should introduce less intra-cluster bias.Additionally, the clusters are surrounded by non-release voxels from the striatum. In the randomrelease magnitude simulations, only the edge voxels of the striatum would be contaminated byadjacent non-release voxels after post-processing. However, these adjacent voxels were not fromthe striatum, and so their TACs have much lower concentration values in the post-interventionregime, effectively corrupting the release TACs less severely. Consequently, in the cluster releasesimulations post-processing should introduce more intra-striatal bias.Thus we have two competing factors: a voxel with release surrounded by both release andnon-release voxels in its local neighbourhood should see less bias introduced from the neighbouringrelease voxels (intra-cluster bias) and more bias introduced from the neighbouring non-release voxels(intra-striatal bias).In Figure 4.20 we present the bias in the estimated DAR(t) for each cluster after post-processingthe noiseless data. In all clusters, the overall effect is that DAR(t) becomes increasingly negativelybiased as the release peaks; this is when the difference between release TACs and non-release TACsare most different.In higher release clusters (C1, C3, and C5), the peak negative bias is more pronounced. This isbecause the drop in the RAC signal in these voxel-level TACs is larger, and therefore when mixedwith adjacent non-release voxels the resulting drop has a greater absolute difference from its groundtruth than that from a smaller ground truth drop.The size of the cluster has a much smaller effect on the amount of bias, with larger clustersbeing slightly more negatively biased. This may run counter to initial intuition: with a largercluster, a smaller fraction of the voxels are heavily affected by PVE, and so the mean bias shouldbe expected to be lower. However, because the simulation was designed to have specific mean DARvalues after the PSF filtering, this required a higher initial release level simulated in smaller clustersbefore the PSF filtering. As a result, after filtering a greater gradient exists (highest release at thecentre, lower release towards the edges, c.f. Figure 3.4). Therefore, for an edge voxel in a smallercluster, the signal is being contaminated by both adjacent non-release voxels (negative bias) andhigher release voxels towards the centre of the cluster (positive bias). Though the same scenariooccurs in larger clusters, because the gradient is not as large, the positive bias from higher releasevoxels towards the centre is smaller—ultimately generating a slightly greater negative bias for largerclusters.In terms of the post-processing methods themselves, since the clusters are always surroundedby non-release voxels, a larger post-filter applied results in a greater negative bias. In terms ofthe kernel size in the HYPR operator, a smaller kernel size produces less negative bias—consistentwith the results from the random release magnitude simulation. This effect is especially noticeablein the caudate—a smaller structure. This supports the idea that the HYPR operator can betterpreserve contrast in small structures with a smaller kernel, though the denoising efficacy will bereduced.Switching to the effect of not having the correct basis function supplied, we again fit the noiseless824.3. Cluster release simulation(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.20: The bias introduced to the estimated DAR(t) of post-processed noiseless data, for each clusterof the cluster release simulation.834.3. Cluster release simulation(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.21: The bias introduced to the estimated DAR(t) of noiseless data when the correct basis functionis not supplied, for each cluster of the cluster release simulation.844.3. Cluster release simulationdata with the standard set of basis functions and plot the estimated DAR(t) alongside the groundtruth curve, for each cluster, in Figure 4.21.For clusters C1 and C3 (high release clusters with size > 400 voxels), the biases are consistentwith those observed in the random release magnitude simulations: the estimated tD is positivelybiased and the estimated tP is negatively biased, in an attempt to conserve (tP −tD)/α when α = 1is selected instead of the ground truth (but unavailable) α = 2.However, for all remaining clusters (lower release and/or smaller in size) there is a tendency fortD to be unbiased or negatively biased, and tP to unbiased or positively biased. This is becausein these clusters α = 4 is most commonly being selected, such that the numerator of (tP − tD)/αneeds to be increased to preserve the area under the curve of DAR(t). This can be achieved by anegative bias in tD and positive bias in tP .Timing parametersRegional means of tD and tP over the detected clusters are plotted in Figure 4.22 for each cluster,using both traditional and uncorrected weighting (16 mCi simulation).All reconstructions perform similarly using either weighting scheme. On average, tD is well-estimated to within a minute of the true tD = 35 min. For tP , the means are well-estimated for thehigher release clusters, but are positively biased for lower release clusters. This is consistent withwhat was observed when plotting tP as a function of release level in the random release magnitudesimulations.As well, neither cluster size nor location within the striatum appear to have a significant effecton the means; rather release level is the limiting factor.Figure 4.23 examines how these mean estimates change with dose for OSEM-5mm-5mm andHYPR-2.5mm-5mm (uncorrected weighting; the comparisons for traditional weighting are similar).In higher release clusters (C1, C3, C5), the results are fairly stable across dose. For lower releaseclusters (C2, C4, C6), the results are less consistent. C2 shows the mean tP increasing moving from20 mCi to 24 mCi, whereas C4 shows the mean decreasing (only for OSEM-5mm-5mm), and C6 isfairly consistent across dose.Though we may anticipate better mean estimates with increasing dose, this variability in lowerrelease clusters can be expected since the sensitivities are low in these clusters, and so a few outliervoxels in terms of estimated timing parameters can cause counterintuitive changes to the meanwith increasing dose.Nevertheless, even in the higher release clusters there is no real gain in the mean estimatedtiming parameters. Thus we turn to MAE in the timing parameters to see how well on average themethods perform in each cluster on the voxel-level.Figure 4.24 plots MAE(tD) and MAE(tP ) for each method in each cluster, using either tradi-tional or uncorrected weighting (16 mCi simulation). Overall, in most clusters OSEM-5mm-5mmproduces the lowest MAE for both timing parameters, with HYPR-2.5mm-5mm and OSEM-5mm-2.5mm performing similarly. In high release clusters, on average OSEM-5mm-5mm (HYPR-2.5mm-5mm) has MAE(tD) = 2.96 min (3.23 min) and MAE(tP ) = 5.06 min (5.13 min). In lower releaseclusters in the putamen (i.e. C2 and C6), these values climb to MAE(tD) = 3.55 min (3.78 min)and MAE(tP ) = 5.50 min (5.68 min). The worst performance occurs in the sole low release clusterin the caudate, C4, which has MAE(tD) = 4.40 min (4.29 min) and MAE(tP ) = 6.05 min (6.01min).Compared to the results from the random release magnitude simulations, these values are con-sistent for HYPR-2.5mm-5mm when interpolating at 7.5% and 15% DAR in the MAE as a functionof release level plots. However at OSEM-5mm-5mm, the MAE values here are slightly higher. In all854.3. Cluster release simulation(a) traditional weighting (b) uncorrected weightingFigure 4.22: Regional means of tD (solid bars) and tP (lightly shaded bars) for each cluster, using traditionalor uncorrected weighting (16 mCi simulation). Ground truth values are indicated by dashed lines.Figure 4.23: The effect of dose level on regional mean estimates of tD (solid bars) and tP (lightly shadedbars) for OSEM-5mm-5mm (yellow) and HYPR-2.5mm-5mm (purple), for each cluster. Ground truth valuesare indicated with dashed lines.likelihood, since the timing parameters are identical in each release voxel, the large filter of OSEM-5mm-5mm is beneficial for MAE in the large clusters of the random release magnitude simulationby spatially correlating all the nearby voxels, but is disadvantageous in smaller clusters when weconsider the greater degree of intra-striatal contamination from nearby non-release voxels. Con-versely, the resolution-preserving properties of HYPR-2.5mm-5mm allow for voxel-level estimatesin smaller clusters to become less biased after denoising, and hence MAE is less affected by the sizeof the cluster in the simulation.Also like in the random release magnitude simulations, traditional weighting produces slightlylower MAE than uncorrected weighting, though this difference is small. The effect of dose on MAEis similar to the effect of dose on mean timing parameter estimates; in high release clusters MAEis relatively stable, and in low release clusters MAE does not consistently improve (Figure 4.25).864.3. Cluster release simulation(a) traditional weighting (b) uncorrected weightingFigure 4.24: MAE(tD) (solid bars) and MAE(tP ) (slightly shaded bars) for each cluster, using traditionalor uncorrected weighting (16 mCi simulation).Figure 4.25: The effect of dose level on MAE(tD) (shaded bars) and MAE(tP ) (unshaded bars) forOSEM-5mm-5mm (yellow) and HYPR-2.5mm-5mm (purple), for each cluster.BP (t) and DAR(t)BP (t) estimates for all clusters are plotted for the five best denoising protocols in Figure 4.26. Thecorresponding DAR(t) estimates are presented in Figure 4.27. Like before in the random releasemagnitude discussion on BP (t) and DAR(t) estimates, we only show the results for uncorrectedweighting, as the effect of weighting scheme used has little effect on these estimates.Consistent with the parametric accuracy analysis of the random release magnitude simulation,OSEM-5mm-5mm provides the highest precision in BP (t) and DAR(t) in all clusters, though withsevere negative biases and hence poor parametric accuracy. Overall, HYPR-2.5mm-5mm providesthe best accuracy in both BP (t) and DAR(t)—especially in the low release clusters. This isa testament to the ability of the denoised reconstruction to provide less noisy voxel-level TACswhile simultaneously preserving functional structures (i.e. clusters of DA release) on the order of874.3. Cluster release simulation(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.26: Comparing BP (t) estimates for all clusters. Coloured lines indicate the cluster-level medianestimate. Shaded bands indicate ± one median absolute deviation. The ground truth is presented in black.There is a spread in the ground truth since the PSF filtering induces a gradient in the release level of eachcluster.884.3. Cluster release simulation(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.27: Comparing DAR(t) estimates for all clusters. Coloured lines indicate the cluster-level medianestimate. Shaded bands indicate ± one median absolute deviation. The ground truth is presented in black.There is a spread in the ground truth since the PSF filtering induces a gradient in the release level of eachcluster.894.3. Cluster release simulationhundreds of voxels. For DAR(t) in high release clusters, all protocols except OSEM-5mm-5mmperform similary in terms of accuracy, suggesting the effect of denoising is less important at highrelease.As was seen in the random release magnitude simulation and consistent with our analysis ofthe post-processed noiseless data, the peak DAR(t) is negatively biased for high release clusters asa result of post-processing. In the noiseless case, peak DAR(t) was also negatively biased for lowrelease clusters. However, in the noisy case we now observe peak DAR(t) being unbiased or slightlypositively biased (except for OSEM-5mm-5mm, which is always negatively biased as a consequenceof its wide spatial filters).This again may be explained by a selection effect: in the cluster release simulation, a gradient ofDAR exists from the PSF filtering, so the higher DAR are easier to detect and may be preferentiallydetected on average, positively biasing the median DAR(t) estimates. As a result, this positivebias either cancels or overpowers the negative bias from post-processing.Another feature of these estimates is an extension of the tail of DAR(t), when the groundtruth DA release has essentially decayed back to baseline levels. This is especially prominent inlow release clusters and in clusters within the caudate (C2, C3, C4, and C6). Since these twocluster categorizations have lower sensitivity, this again suggests a selection effect: voxels withnoisy perturbations or biases that tend to extend the apparent release level will have a greaterdifference in WRSS between SRTM and lp-ntPET fits, resulting in a higher F value and hence abetter chance of being detected.When considering cluster size, the smaller clusters tend to have higher overall accuracy inDAR(t) (compare the smaller C2 with the larger C6, or the smaller C5 with the larger C1). Thismay be because the larger clusters are more elliptical in shape, since they have to fit within theconfines of the “nutshell”-shaped putamen. A more elliptical shape means more voxels are locatedcloser to the edge of the cluster. This introduces more intra-striatal contamination to a greaterfraction of the voxels in these clusters, ultimately resulting in lower accuracy. Thus the shape ofthe ground truth cluster may also factor into the resulting parametric accuracy.While extra-striatal contamination negatively biases the entire TAC to a greater degree thanintra-striatal contamination, the extra-striatal contamination roughly equally targets all time points,whereas the intra-striatal contamination greatly affects the post-intervention period. As a result,the fitted release level is more negatively biased by intra-striatal contamination than extra-striatalcontamination (cf. Figure 4.28).An example of this effect in action can be seen by comparing the peak DAR(t) estimates inC3 (L caudate, 15% DAR, 400 voxels) and C1 (L putamen, 15% DAR, 500 voxels). Since weknow the set of striatal voxels close to the edge of the caudate have already been negatively biasedby adjacent extra-striatal voxels (and thus their TACs are closer to non-release voxels outside thestriatum), and there is significant overlap between this set of voxels and the set of non-releasevoxels surrounding C3, the intra-striatal contamination introduced by post-processing in effect isa mixture of pure extra-striatal contamination and pure intra-striatal contamination, meaning theoverall effect is less bias introduced into DAR(t) by post-processing in C3 than in C1. This maybe the primary effect resulting in a less biased peak DAR(t) in C3, although we acknowledge theaforementioned selection effect in the lower sensitivity caudate may also play a role.In terms of weighting scheme, again there is negligible difference in the precision and accuracy inthe estimated BP (t) and DAR(t), so we can conclude that the effect of weighting is only meaningfulwhen considering CST and sensitivity.Finally, as was the case in the random release magnitude simulation, when dose level is increasedthere is little to no gain in accuracy in estimating BP (t) and DAR(t), and so we do not show theseresults.904.3. Cluster release simulationFigure 4.28: Demonstrating the difference between intra-striatal bias introduced to a release voxel by anon-release striatal voxel, and PVE bias introduced by a non-release extra-striatal voxel.4.3.3 Comparing results using previous noise-modeling methodsFigures 4.29-4.31 compare the voxelwise, binary50, and binary sensitivity, respectively, resultingfrom the simulations using the Gaussian noise addition methods of previous work [57, 60], and fromOSEM-5mm-5mm and HYPR-2.5mm-5mm in this work.(a) traditional weighting (b) uncorrected weightingFigure 4.29: Comparing voxelwise sensitivity between simulations using Gaussian noise addition methodsand simulations using the more realistic noise addition methods of this work.914.3. Cluster release simulation(a) traditional weighting (b) uncorrected weightingFigure 4.30: Comparing binary50 sensitivity between simulations using Gaussian noise addition methodsand simulations using the more realistic noise addition methods of this work.(a) traditional weighting (b) uncorrected weightingFigure 4.31: Comparing binary sensitivity between simulations using Gaussian noise addition methods andsimulations using the more realistic noise addition methods of this work.All types of sensitivities are higher in simulations using Gaussian noise addition. The gains areparticularly significant in low release clusters: at least a portion of C2, C4, and C6 are detectedupwards of 90% of the time, and voxelwise sensitivities are on the order of ∼ 50%. In contrast, thehighest sensitivity denoising protocol when using the more realistic noise addition methods of thiswork (OSEM-5mm-5mm) have maximum binary sensitivities of ∼ 60% and maximum voxelwisesensitivities of ∼ 40%.The cluster-level mean estimates of the timing parameters as well as MAE in the timing pa-rameters are compared between methods as well, in Figure 4.32 and Figure 4.33, respectively.924.3. Cluster release simulationFigure 4.32: Comparing the cluster-level mean estimates of timing parameters in the cluster release sim-ulation using Gaussian noise addition, and in the same simulation using the more realistic noise-modelingmethods. Ground truth values are indicated by dashed lines.Figure 4.33: Comparing MAE(tD) and MAE(tP ) in the cluster release simulation using Gaussian noiseaddition, and in the same simulation using the more realistic noise-modeling methods.When Gaussian noise addition is used, the mean timing parameter estimates are closer to what isobserved in the noiseless case: tD is slightly positively biased, and tP is slightly negatively biased.These results are consistent across all clusters. While tD estimates are also slightly positivelybiased for OSEM-5mm-5mm and HYPR-2.5mm-5mm, tP estimates are generally positively biasedin clusters with lower sensitivity (C2, C3, C4, and C6). This lends credence to the hypothesis thatTACs biased to appear to have extended DAR(t) (and hence increased tP − tD) are more easilydetected, and thus generate a selection effect when sensitivity is lower. This effect is not presentin the Gaussian noise simulations, since the sensitivity is significantly higher.MAE in the timing parameters is also consistently lower when Gaussian noise addition is used,indicating that the voxel-level estimates also are more accurate and have lower variance. This is alsoreflected in the estimated BP (t) and DAR(t), where the curves better resemble the shape of theground truth (albeit still with the usual biases in peak DAR(t) induced by HYPR post-processing)and have lower MAD. Since these comparisons are similar to those in the random release magnitude934.3. Cluster release simulationsimulation, we do not include these plots in this work for brevity.Again, these results suggest that simply adding Gaussian noise to noiseless data and not in-cluding the effects of attenuation, normalization, and reconstruction can overestimate the efficacyof the detection algorithm, or that the method of adding Gaussian noise underestimates the truenoise levels, such that the detection ability is improved.4.3.4 ConclusionsThe cluster release simulation serves as a more physiologically realistic simulation (when consideringa neurophysiological task) as compared to the random release magnitude simulation, with smallerclusters scattered throughout the striatum instead of entire subregions of the striatum being filledwith release.Whereas OSEM-5mm-5mm always provides the highest sensitivity in the random release mag-nitude simulations, here the difference in sensitivity shrinks between OSEM-5mm-5mm and theother viable methods, or in some cases OSEM-5mm-5mm no longer becomes the best method.In smaller, lower release clusters, when using a method with a higher CST such as OSEM-5mm-5mm, the portion of the true cluster detected may fall below the CST, and hence will berejected in some realizations. Therefore when we compare binary sensitivity—the percentage ofnoisy realizations where any portion of the true cluster passes through the CST—OSEM-5mm-5mm performs similarly or worse than the other methods. However, when a portion is detected,the percentage of true release voxels detected tends to be higher for OSEM-5mm-5mm, leading toa higher voxelwise sensitivity.So which metric should we use to decide which method to use? Considering a realistic applica-tion of this protocol, only one scan is performed per patient, so binary sensitivity will tell us theprobability that any release will be detected for the given cluster type in this single scan scenario.Since binary sensitivity is roughly the same for all methods, OSEM-5mm-5mm may still be the bestmethod to use since it tends to detect larger portions of true clusters when they pass through theCST. However, with OSEM-5mm-5mm one runs the risk of missing small clusters (especially onessmaller than those simulated in this work) in PVE-affected regions as a result of the high CST.Another line of thinking would be that since all methods at least detect a portion of trueclusters with the same probability, one may then defer to quantitative accuracy to break the tie.In this simulation, having smaller clusters surrounded by non-release voxels makes the effect ofintra-striatal contamination more prevalent after post-processing, especially for methods with widefiltering kernels such as OSEM-5mm-5mm.Indeed we find that the accuracy in estimating the entire release dynamics is poorest for OSEM-5mm-5mm, despite having the best precision. The best parametric accuracy comes from the methodwithout any spatial post-filter, OSEM-NF-5mm, or the resolution preserving HYPR-2.5mm-5mm.Considering both of these methods share similar sensitivity, but the precision is much improved forHYPR-2.5mm-5mm, if one values parametric accuracy as well as sensitivity then HYPR-2.5mm-5mm would be the optimal choice.In terms of weighting schemes, again we find a small increase in sensitivity by using traditionalversus uncorrected weighting. However, again if we are using denoising protocols with CSTs com-parable to the size of true clusters, then because traditional weighting produces a higher CST werun the risk of missing these smaller clusters. Additionally, weighting scheme has very little effecton parametric accuracy.When considering dose, in almost all cases we see gains in sensitivity after increasing dose.However these sensitivity gains do not come with significant improvements in parametric accuracy.Therefore we reiterate the conclusions from the random release magnitude simulation that a dose944.4. Monte Carlo methodof 20 mCi balances high sensitivity with ethical and economic considerations.These two simulations have identified OSEM-5mm-5mm and HYPR-2.5mm-5mm as the twomost optimal denoising protocols for detecting transient DA release, when using F-test and clustersize thresholding methods. Now we consider modifications to the theory of the detection method:using the MC method in an attempt to improve the sensitivity beyond the results shown here,without compromising parametric accuracy.4.4 Monte Carlo methodIn the previous section, we have found that OSEM-5mm-5mm provides high voxelwise sensitivitywith the trade-off of poor parametric accuracy, whereas HYPR-2.5mm-5mm has lower voxelwisesensitivity but much higher parametric accuracy. As well, the CST, which defines the smallestpossible region of release that can be detected, is 162% (118%) larger for OSEM-5mm-5mm whenusing traditional (uncorrected) weighting.In an ideal scenario, we would like to draw from the high sensitivity of OSEM-5mm-5mm whilealso drawing from the high parametric accuracy of HYPR-2.5mm-5mm. A na¨ıve first approachwould be to use the clusters detected from OSEM-5mm-5mm as ROIs for the HYPR-2.5mm-5mmimages. However, this approach suffers from two disadvantages: i) two separate reconstructionsneed to be performed; and ii) we would still be constrained by the high CST of OSEM-5mm-5mm.The alternative approach is to apply the MC method to the HYPR-2.5mm-5mm images. Sincethe MC method does not modify the images, by definition the voxel-level parametric accuracy ispreserved (though the additionally-detected voxels may introduce a non-negligible different to thecluster-level mean of the voxel-level parameter estimates). Theoretically, the MC method shouldextend the size of the clusters to recover false negatives and correctly reclassify them as releasevoxels, thus increasing sensitivity. Simultaneously, the regularization of maintaining a proper F-distribution should also reduce the CST for a given sensitivity.Thus in this section we present the results from applying the MC algorithm to the HYPR-2.5mm-5mm images for all three simulation types, at 16 mCi. We elect to use uncorrected weightingso as to start with the smallest possible CST, even though the overall sensitivity is slightly lower: theaverage sensitivity increase when using traditional as opposed to uncorrected weighting is 12%11,but this comes with a 31% increase in CST.We shall compare the MC method to the standard F-test plus cluster size thresholding methodapplied to OSEM-5mm-5mm (SM-OSEM) and to HYPR-2.5mm-5mm (SM-HYPR-AU).As a reminder, the MC method replaces the F-test step in the original method with a newstatistic c derived from the MC algorithm. Each voxel in the subregions of the striatum areclassified as either a release on non-release voxel, comprising a classification state. The goal of theMC algorithm is to shuffle through the likely classification states of each subregion probabilistically,determined by an objective function that considers both the local parity of states and the globaladherence of the non-release voxels to an F-distribution. After many iterations of the MC algorithm,c at a given voxel represents what percentage of the time that voxel was classified as a release voxel.If the assumptions of the objective function hold true, then c can be interpreted as a percentageconfidence that release exists in a given voxel. Another way to view c is as a more sophisticated F-statistic, which takes the local spatial information into account while maintaining the assumptionsof the null hypothesis of the F-test.11This change is computed from the average voxelwise sensitivity over all clusters in the cluster release simulation.954.4. Monte Carlo method4.4.1 Choosing the optimal hyperparameter β and thresholding value ccritWe begin by presenting the results from determining the optimal values for the hyperparameter βand thresholding value ccrit. Heuristically, β controls the amount of regularization in the algorithm.Too low a β means the parity term dominates in the objective function, effectively “overclustering”and increasing the false positive rate. Too high a β will minimize the false positive rate, but willalso reduce the sensitivity.Therefore, as mentioned in Section 2.7.3, our first estimate for β is the value that roughlyweights each term equally in the objective function. After histogramming the ∆Pi and ∆R fromeach step of each iteration when applying the algorithm to the cluster release simulation, we findβ ∼ 8× 105 to produce roughly equal contribution.To consider other nearby values, we test β ∈ {8× 104, 4× 105, 8× 105, 4× 106, 8× 106} for allrealizations of the baseline and cluster release simulations, at 16 mCi. We determine the best valuefor the hyperparameter by analyzing curves of sensitivity versus CST, created by thresholding theoutputted c images from each β at ccrit ∈ [65%, 90%]. The sensitivity value is the average voxelwisesensitivity over all clusters. We choose to use sensitivity from the cluster release simulation insteadof the random release magnitude simulation as we believe it to be a more physiologically realisticsimulation.Figure 4.34 shows the curves of sensitivity versus CST for different values of β. Moving fromleft to right along a curve corresponds to decreasing ccrit from 90 to 65 (i.e. allowing less confidentvalues to be accepted, both increasing CST and sensitivity).Unlike the F-test where theory exists to convert the F-statistic to a probability, the choiceof ccrit has no equivalent conversion. However, we can compare the results of the original F-testmethod in terms of equivalent CST or equivalent sensitivity to determine the optimal {β, ccrit}pair.From the curves of Figure 4.34, it is clear that all values of β perform similarly, but that thestarting points (ccrit = 90%) and ending points (ccrit = 65%) are different. This indicates that theparameters β and ccrit are coupled; each β can achieve a similar sensitivity at a given CST, butthe ccrit at which point this occurs will be different for each β.For instance, if we choose a relatively low β, the algorithm will not regularize clustering asmuch. If we instead chose a higher value of β, the voxelwise values of c will be lower from greaterregularization. Therefore to achieve identical performance in terms of sensitivity and CST, thelower choice of β will require a higher ccrit: the lower regularization raises voxelwise c, and henceto avoid a large CST a greater thresholding value, ccrit, is required.This is why the lowest value tested, β = 8×104, starts and ends the farthest right of Figure 4.34:lower regularization means that the true positive rate—as well as the false positive rate—willincrease, so for any given ccrit the CST will be higher for the lowest value of β.These results are encouraging, since the choice of β seems relatively unimportant and hence thealgorithm does not require extensive fine-tuning: any β within the range tested in this work maybe selected, and a corresponding ccrit may be chosen to achieve a desired sensitivity or CST.Two optimizations of the MC methodThe optimal location on the plot is the top left corner (low CST, high sensitivity). Comparedto SM-OSEM and SM-HYPR-AU, for any β chosen the MC method has a higher sensitivity atequivalent CST to the standard methods, or lower CST at equivalent sensitivity.These results motivate designing two optimizations of the MC method. The first maximizes sen-sitivity without exceeding the CST of the best performing standard method in terms of CST—SM-964.4. Monte Carlo methodFigure 4.34: Curves of sensitivity versus CST for each tested β of the MC method, used to determine anoptimal implementation. Moving from left to right along a curve corresponds to decreasing ccrit from 90 to65. The optimal location is the top left corner of the plot. The sensitivity and CST pairs for SM-OSEM andSM-HYPR-AU are provided for reference.HYPR-AU with a CST of 108 voxels. We shall call this method MC-CST. The second optimizationmaximizes sensitivity without exceeding the CST of the best performing standard method in termsof sensitivity—SM-OSEM with a CST of 235 voxels. This method will be called MC-SENS.Starting with MC-CST, of the {β, ccrit} pairs tested, {β = 8 × 105, ccrit = 70%} producesthe highest voxelwise sensitivity with a CST that does not exceed 108 voxels (46.9% voxelwisesensitivity with a CST of 104 voxels). Although in Figure 4.34 it may appear that at CST ∼ 108voxels β = 4× 105 or β = 8× 104 performs slightly better in terms of sensitivity than β = 8× 105.However, these points on the curve correspond to interpolations between one tested ccrit to thenext, and so we cannot be sure that these values of sensitivity are indeed the value we would getby testing and intermediate ccrit. Besides, the various values of β tested perform very similarly atthis CST, and thus the choice of β here makes little difference.For MC-SENS, the highest voxelwise sensitivity produced by the MC method is naturally fromthe least regularized value tested, β = 8×104, with the lowest value of ccrit, yielding 56.0%. Despitethis, the CST is still lower than that of the best performing standard method in terms of sensitivity,SM-OSEM, at a value of 173 voxels. Thus these MC parameter values are used for MC-SENS.Lower values than ccrit = 65% (not shown in Figure 4.34) cause the CST to increase in largesteps for small changes in ccrit, since c ∼ 50% represents a voxel whose classification is uncertainand thus constantly being changed. This is why we selected a maximum of ccrit for our finaloptimization procedure.As well, lower values of β effectively do not regularize since the parity term of equation (2.84)dominates. Though within the MC Metropolis framework a weaker form of F-distribution regular-ization occurs when classification states having ∆E < 0 are probabilistically accepted based on theF value in that voxel, in effect much lower values of β or even β = 0 cause the algorithm to performmore like a filter on the binary sensitivity map by encouraging clustering everywhere, instead ofpreferentially clustering around areas with F values that would lead to a better global adherenceto an F-distribution.974.4. Monte Carlo methodDemonstrating cost function equilibrium and reproducibility of the algorithmFor both MC-CST and MC-SENS, the cost function tends to reach a constant minimum value bythe third iteration (Figure 4.35), before oscillating in “equilibrium” about this constant value. Wethus use 100 subsequent iterations to calculate parametric c images. During these oscillations inthe cost function, the algorithm is presumably shuffling through likely states that satisfy both aglobal F-distribution and local parity. Voxels with high probability of DA release will remain inthe state s = 1 for almost all iterations, whereas in voxels where this is less certain the state willbe switched much more frequently.Figure 4.35: Demonstrating how the cost function of the MC algorithm quickly reaches equilibrium aftera few iterations, before oscillating about a constant minimum value.In terms of reproducability, after running 100 trials of MC-CST on a single noisy realizationof the cluster release simulation, we found a voxelwise sensitivity of 57.2 ± 0.3%, demonstratinglittle variance in the measured outcome of the algorithm albeit being non-deterministic—especiallycompared to the range of sensitivities across all noisy realizations (46.0± 10.9%).Optimizing on other scannersIf these methods are to be applied to a different scanner with roughly equivalent scanner sensitivity(i.e. percentage of real annihilation events detected), the curves of Figure 4.34 will likely give areasonable indication of the performance of the MC method. Higher scanner sensitivity may resultin less noisy TACs, and thus likely produce higher average voxelwise sensitivity at a given CST.Since the HRRT (the scanner used for the simulations of this work) is currently the highestresolution PET human brain scanner, all other scanners will have lower resolution and thus likelylarger voxel sizes. As a rough approximation, the CSTs of this work can be scaled by the ratios ofthe scanner resolutions (i.e. for a scanner with isotropic 5 mm resolution, the CSTs of that scannercould be estimated by multiplying the CSTs of this work by (2.5/5)3 = 0.125). However, as wehave seen, a multitude of factors contributed to CST and voxelwise sensitivity, and so it wouldbe recommended to perform a similar process of simulating baseline and DA release scans on thedifferent scanner—with any preferred reconstruction method—to arrive at the best estimates of theCST required and the corresponding performance in terms of voxelwise sensitivity.984.4. Monte Carlo method4.4.2 Baseline simulationAs mentioned, MC-CST has a CST roughly equal to that of SM-HYPR-AU (104 and 108 voxelsrespectively) with an average voxelwise sensitivity approximately equal to SM-OSEM. Thus we canconsider these methods approximately equal to each other in terms of performance—if the goal isto have the highest possible true positive rate without worrying about CST. However, if one isconcerned about detecting small structures as well, then MC-CST is a more optimal overall, beingable to detect structures 2.3 times smaller by virtue of its much lower CST.MC-SENS, on the other hand, offers only a 26% reduction in CST over SM-OSEM, but hashigher average voxelwise sensitivity than both SM-OSEM and MC-CST.The reductions in CST in both MC methods stem from the enforcement of the baseline classifiedvoxels obeying a true F-distribution. In the baseline simulations, since all voxels are baseline, theF-distribution of all voxels generally well approximate a true F-distribution, meaning we expect asmall subset of voxels within each baseline simulation to have F > Fcrit. In the standard method,these voxels are always initially accepted as release voxels. However, the MC method is built onthis assumption that a set of baseline voxels will obey an F-distribution, and so these F > Fcritvoxels can still be rejected. As a result, overall the CST is reduced for equivalent sensitivity.The MC method does not classify all voxels in the baseline simulations as baseline voxels(though there are some noisy realizations where this does occur), since this is only one of manyglobal classification states which satisfy a true F-distribution; in theory one can proportionallyreclassify a subset of voxels from the span of all possible F values, such that the set of classifiedbaseline voxels remains obeying a true F-distribution.As well, there are anomalous noisy realizations in the baseline simulations where the F valueshave a large number of F > Fcrit. In the standard method, these realizations have the highest falsepositive cluster sizes since all F > Fcrit pass through significance thresholding. In the MC method,since the tail-heavy set of F values in these realizations deviates from an F-distribution, a subsetof these F > Fcrit voxels are classified as release voxels such that only an appropriate number ofF > Fcrit remain classified as baseline voxels. Therefore, these realizations also have the highestfalse positive cluster sizes, albeit still being smaller than the those in the standard method.4.4.3 Random release magnitude simulationSensitivity as a function of release levelApplying MC-CST and MC-SENS to the random release magnitude simulation using the sameset of ccrit, we find that both methods offer significant improvements over SM-HYPR-AU (cf.Figure 4.36). Since all three of these methods use the same initial input data (i.e. parametric Fimage derived from the HYPR-AU-OSEM images), the this demonstrates the power of the MCmethod over the standard method.MC-CST improves sensitivity by 77% at the lowest release levels and by 17% at high releaselevels over SM-HYPR-AU, whereas the improvements for MC-SENS are 103% and 21%, respec-tively.However, when comparing to SM-OSEM—the best performing standard method in terms ofsensitivity—MC-CST still performs similarly and MC-SENS performs slightly better over all releaselevels (a 15% improvement at the lowest release levels and a 1% improvement at the highest releaselevels). The improvements at higher release are not as drastic, since the sensitivity at these levelsis already fairly high.These results demonstrate that the MC methods offer similar performances in term of sensitivityto SM-OSEM in the random release magnitude simulation, though at a much smaller CST and using994.4. Monte Carlo methodFigure 4.36: Comparing sensitivity as a function of DAR in the random release magnitude simulation,between the standard and MC methods.Figure 4.37: Axial, coronal, and sagittal views of the voxelwise sensitivity in the random release magnitudesimulation, for the standard and MC methods.the more parametrically accurate images of HYPR-AU-OSEM. The much smaller CST of the MCmethods is not as important in these simulations, since the region of release is extended throughoutthe striatum and thus is much larger in size than the CST.The improvements in voxelwise sensitivity can be visualized in Figure 4.37 for representativeslices of the striatum. We observe that the improvements in sensitivity are preferentially targetedtowards voxels near the edges and/or the caudate, whose TACs are more susceptible to PVE andhence whose F values more commonly fall below Fcrit set in the standard method.Validity of the CSTAs mentioned in Section 4.2.3, both SM-OSEM and SM-HYPR-AU achieve a 94% and 98% falsepositive-free rate, respectively, in the R putamen of the random release magnitude simulation,which contained no release to verify the CST is indeed working to its intended 90% false positive-free rate. The CSTs determined for MC-CST and MC-SENS also are functioning correctly, both1004.4. Monte Carlo methodachieving a 96% false positive-free rate over the 50 simulations tested.Timing parameter accuracyWhile the MC undoubtedly provides a boost in sensitivity, a valid question to ask is what effectdo these additionally-detected voxels have on parametric accuracy. After all, the additional voxelsdetected were rejected by the SM-HYPR-AU, suggesting that the SNR of these TACs is likelypoorer and thus potentially could have more inaccurate parameter estimates in their lp-ntPET fits.Note, however, that because the MC methods do not modify the original HYPR-AU-OSEM images,the voxel-level parametric accuracy cannot be changed from SM-HYPR-AU. What we want to ask,instead, is what impact do the additionally-detected voxels in the MC method have on the meanvoxel-level parametric accuracy within detected clusters.Starting with the mean estimates of the timing parameters, we compare the results of theMC-CST and MC-SENS to the standard methods in Figure 4.38.For tD, the both MC-CST and MC-SENS yield slightly higher estimates on average than thestandard methods—further away from the ground truth tD = 35 min. However, for tP the MCmethod also estimates a slightly higher value on average, but this is closer to the ground truthtP = 40 min.Since the MC methods draw from the same set of fitted TACs as the original method, this meansthat the additionally-detected voxels are higher on average than the mean tD and tP estimates fromthe standard methods.Glancing at the distributions of tD and tP for all four methods in Figure 4.39, we see thisis indeed the case. For tD, the MC methods proportionally detect more voxels with estimatedtD ≥ 39 min than the original method and less voxels with tD ≤ 37.5 min. For tP , the MC methodsproportionally detects more voxels with estimated tP ≤ 40.5, and more voxels with tP ≥ 42 min.These additionally-detected, positively biased voxels may be from poorer fits to noisier data, suchthat the F values within these voxels are lower and hence were rejected by SM-HYPR-AU.The distributions of tD also demonstrate that although the mean estimated tD is close to theground truth tD = 35, the distributions are peaked at tD = 37.5 min, but a number of voxels withtD = 30 min skews the mean back towards tD = 35 min. This suggests the fitting algorithms maybe overfitting to noise just before the true intervention time in a number of voxels. Conversely, fortP the distribution is well-centered at the ground truth tP = 40 min for all methods.These different distributions also predictably have an effect on the MAE of the timing param-eters. Figure 4.40 compares the MAE for each timing parameter as a function of release levelbetween the MC and original methods.In a similar result, MAE(tD) is slightly higher in the MC method than SM-HYPR-AU, andMAE(tP ) is slightly lower, though these differences are quite small: for MC-CST, MAE(tD) is 11seconds higher on average and MAE(tP ) is five seconds lower on average; for MC-SENS these are14 seconds higher and seven seconds lower, respectively.The values for SM-OSEM are shown for reference; from Section 4.2.2 we have already discussedhow the voxel-level TACs are less noisy in the OSEM-5mm-5mm images, thus making the accuracyof the timing parameters higher. As well, the larger spatial kernel homogenizes the TACs locally,such that the voxel-level variance in estimated timing parameters is low, resulting in lower MAE.From this, we can conclude that the voxel-level timing parameter estimates from additionaldetected voxels from the MC method, as well as the voxels now rejected by the MC method,are slightly less biased for tP , but slightly more biased for tD, when compared to the standardmethod using the same input images (SM-HYPR-AU). When the cluster-level mean of the voxel-level estimates is taken, correspondingly the estimated tP is slightly closer and the estimated tD is1014.4. Monte Carlo methodFigure 4.38: Comparing the mean estimated tD (upper plot) and tP (lower plot) as a function of DAR,between the standard and MC methods. Black lines indicate ground truth values.Figure 4.39: Comparing distributions of estimated tD and tP between the standard and MC methods.Figure 4.40: Comparing the MAE of the estimated timing parameters in the random release magnitudesimulation, between the MC and original methods.1024.4. Monte Carlo methodslightly farther from the ground truth. That said, the differences are so small that we can effectivelysay that there is nothing significantly gained or lost in terms of timing parameter accuracy by usingthe MC method.BP (t) and DAR(t)While the timing parameter accuracy is only slightly affected by using the MC method, the nextquestion to ponder is whether the combination of estimating both the timing and magnitude of therelease—captured by BP (t) and DAR(t)—is affected by using the MC method.The median BP (t) for low release (7% < DAR < 8%) and high release (15% < DAR < 16%) inthe MC methods are presented side-by-side with those from the standard methods in Figure 4.41.The corresponding DAR(t) are shown in Figure 4.42.At both release levels, BP (t) is similarly estimated during peak release in SM-HYPR-AU andthe MC methods. Overall, the curves appear to be slightly shifted downwards. This suggests thatrelease voxels with lower pre-intervention binding (those near the edges of the striatum affected byextra-striatal contamination; the extra-striatal voxels have lower BPND) are being detected morefrequently, which is what was observed in the visualization of voxelwise sensitivity in Figure 4.37.For DAR(t), the higher sensitivity of the MC methods means the selection effect at low re-lease (where voxels having positive DAR(t) biases are preferentially detected) appears to be lessprominent, since DAR(t) is better estimated than SM-HYPR-AU. At higher release, the estimatedDAR(t) are fairly similar, likely because the boost in sensitivity is not as drastic at this releaselevel. However, the MAD in DAR(t) estimates from the MC methods are noticeable higher forthe MC methods. This indicates that the additionally-detected voxels by the MC methods havea greater variance, which can be expected since these voxels failed to pass through significancethresholding in SM-HYPR-AU.(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.41: Comparing estimated BP (t) for low release (7% < DAR < 8%) and high release (15% <DAR < 16%) between the standard and MC methods. Solid lines indicate the median estimate at each timepoint. Shaded bands indicate ± one median absolute deviation.1034.4. Monte Carlo method(a) 7% < DAR < 8% (b) 15% < DAR < 16%Figure 4.42: Comparing estimated DAR(t) for low release (7% < DAR < 8%) and high release (15% <DAR < 16%) between the standard and MC methods. Solid lines indicate the median estimate at each timepoint. Shaded bands indicate ± one median absolute deviation.4.4.4 Cluster release simulationSensitivityFrom the above section, we found that the MC methods greatly improve voxelwise sensitivity—especially in low release regions—while maintaining relatively low CST. We then expect the clusterswith lower release (C2, C4, C6) to show greater relative increases in sensitivity. SM-OSEM also hasshown dramatic increases in low release sensitivity in the random release magnitude simulation, butthese improvements were mitigated in the cluster release simulation where the release was confinedto a smaller region, and hence increased the likelihood of rejection by the relatively large CST ofSM-OSEM.Table 4.16 compares the voxelwise sensitivity in each cluster between the standard and MCmethods. Comparing the MC methods to the standard method using the same input images (SM-HYPR-AU), MC-CST offers an average improvement of 43% in low release clusters (C2, C4, andC6) and a 16% improvement in high release clusters (C1, C3, and C5). For MC-SENS, theseimprovements 149% and 38%, respectively. These latter improvements surpass SM-OSEM in allclusters as well, with average increases in voxelwise sensitivity of 57% in low release clusters and5% in high release clusters.Table 4.16: Comparing voxelwise sensitivity in each cluster between the standard and MC methods.Boldface indicates the best performing method for each cluster. See Table 4.9 for the characteristics (location,size, release level) of each cluster.Method C1 C2 C3 C4 C5 C6SM-HYPR-AU 60.1% 20.0% 44.6% 10.5% 64.5% 12.5%MC-CST 76.7% 39.7% 52.0% 15.9% 72.3% 22.3%SM-OSEM 81.8% 34.0% 58.7% 12.6% 82.5% 24.2%MC-SENS 84.4% 48.8% 64.0% 23.1% 84.0% 35.1%1044.4. Monte Carlo methodTable 4.17: Comparing binary50 sensitivity in each cluster between the standard and MC methods. Bold-face indicates the best performing method for each cluster. See Table 4.9 for the characteristics (location,size, release level) of each cluster.Method C1 C2 C3 C4 C5 C6SM-HYPR-AU 74% 10% 42% 2% 84% 0%MC-CST 98% 44% 62% 6% 82% 10%SM-OSEM 100% 42% 76% 8% 100% 20%MC-SENS 98% 64% 90% 14% 98% 30%Table 4.18: Comparing binary sensitivity in each cluster between the standard and MC methods. Boldfaceindicates the best performing method for each cluster. See Table 4.9 for the characteristics (location, size,release level) of each cluster.Method C1 C2 C3 C4 C5 C6SM-HYPR-AU 100% 54% 98% 38% 100% 52%MC-CST 100% 78% 100% 48% 100% 64%SM-OSEM 100% 64% 96% 26% 100% 56%MC-SENS 100% 84% 98% 50% 100% 80%Table 4.17 and Table 4.18 compare binary50 and binary sensitivity, respectively, between thestandard and MC methods. MC-SENS provides the highest binary50 and binary sensitivity forall methods in all clusters (except for a couple cases where MC-SENS has a value of 98% and theleading method has a value of 100%).Since in a real application of this method, there is only a single scan performed per participant,the most important metric indicating the efficacy of the method is binary sensitivity: what is theprobability that any portion of a true cluster will be detected. High release clusters predictablyhave binary sensitivity approaching 100%, meaning that at least a portion of the cluster will almostalways be detected. From here, voxelwise sensitivity acts as a supplementary metric, informing uswhat fraction of the extent of the true release cluster will be detected. Since MC-SENS offers thehighest voxelwise sensitivity in all clusters, it is the optimal method for detecting higher releaseclusters.However, for low release clusters—where binary sensitivity is much less than 100%—the largeimprovements of MC-SENS are even more valuable: in low release clusters of the putamen, binarysensitivity when using MC-SENS is raised to an average of 82% from 53% (60%) when using SM-HYPR-AU (SM-OSEM). In the caudate, the binary sensitivity is raised to 50% from 38% (26%).Thus the probability of completely missing a low release cluster in any given scan is markedlydecreased.This is where the MC method is most useful. In SM-HYPR-AU, some voxels in the lowerrelease clusters may have F values sufficiently high to pass the significance threshold, and manyof the remaining voxels will have F values slightly below the threshold. Collectively then, thesignificant voxels will be rejected by the CST. Over many noisy realizations, this leads to lowbinary sensitivity. In other words, a higher probability that this specific cluster will not be detectedin any given scan.However, once this spatial distribution of F values is passed to the MC method, rejecting allvoxels of the cluster will distort the F-distribution, and hence the MC method immediately correctsthis by allowing some voxels with F values just below the threshold and clustered about voxels withF values above the threshold to be accepted. The end result is a larger cluster of accepted voxels,1054.4. Monte Carlo methodFigure 4.43: Axial, coronal, and sagittal views of the voxelwise sensitivity in the cluster release simulation,for the standard and MC methods.allowing the cluster to pass through the CST and ultimately boost binary sensitivity.These improvements can be visualized in Figure 4.43, where the voxelwise sensitivity in repre-sentative slices of the striatum is provided. As expected, most of the improvements are near theedges of clusters and in the caudate clusters, where PVE are most prominent.One can immediately notice, however, that SM-OSEM and MC methods occasionally mergetogether clusters within the putamen; for SM-OSEM this is likely the wide spatial kernel correlatingthe baseline voxels separating the two clusters, allowing for these voxels to be incorrectly assignedas release voxels on occasion. For the MC methods, the processed HYPR-AU-OSEM images alsolikely introduce some correlations between these voxels but to a lesser degree, meaning the F valuescommonly fall below the threshold, but are reclassified by the MC algorithm as release voxels onoccasion.Timing parameter accuracyThe mean voxel-level timing parameter estimates for each cluster are compared between the stan-dard and MC methods in Figure 4.44. The results for tD largely mirror those in the random releasemagnitude simulation analysis: the estimates in the MC methods are slightly higher, but not ap-preciably so. This comes from the MC method more frequently selecting positively biased tD, asseen when comparing the distribution of accepted tD in each method (Figure 4.45).The difference now is that the mean estimates of tP in the MC method are also slightly higherin most clusters, and most noticeably so in C4 (R caudate, 400 voxels, 7.5% release), where thedifference between the two methods is almost a minute.Since in C4 the sensitivity is boosted by almost 200%, this suggests that the large number ofadditional voxels added by the MC method have tP positively biased. This is confirmed by againcomparing the distributions of estimated tP accepted in each method (Figure 4.45).Compared to the distributions from the random release magnitude simulation, in tD againa fraction of the detected voxels have fits before the intervention at the lowest available fittingparameter time of tD = 35 min. The MC methods tend to have less fits before the true interventiontime, and more fits after the true intervention time, resulting in the overall positive bias. For tD,the distribution is again centred around the correct tP = 40 min, but the variance is much higher.These observations suggest that the smaller region of release combined with the larger amount of1064.4. Monte Carlo methodFigure 4.45: Comparing distributions of estimated tD and tP in the cluster release simulation, betweenthe standard and MC methods.Figure 4.44: Comparing the cluster-level mean estimates of tD and tP in each cluster between the standardand MC methods. Dashed lines indicate the ground truth values.1074.4. Monte Carlo methodFigure 4.46: Comparing the MAE of the estimated timing parameters in the cluster release simulation,between the MC and original methods.PVE from intra- and extra-striatal contamination result in poorer voxel-level timing parameterestimates.Since the MC method will predominately add voxels more affected by PVE, this leads to agreater bias in the cluster-level mean estimates of timing parameters—most appreciably in lowrelease clusters where the sensitivity improvements are most drastic.The MAE in the timing parameters is compared in Figure 4.46 between the standard and MCmethods. MAE(tD) can be seen to increase when using the MC method, but as in the case ofthe random release magnitude simulation, MAE(tP ) decreases. While the mean tP is higher inthe MC method than the standard methods, in the high release clusters (where the sensitivity ishighest and thus the voxels here contribute most to MAE) the mean is closer to the ground truthtP = 40 min, leading to the lower MAE.Recall, however, that in the case of noiseless post-processed data without providing the correctlp-ntPET basis function, the algorithm tends to negatively bias tP . Therefore this suggests thatthe additional voxels added by the MC method—while true release voxels—have TACs modified bynoise and/or bias in the reconstruction that results in fits with positive bias. Nonetheless, at leastin high release clusters these competing biases cancel each other out, leading to lower MAE in tP .BP (t) and DAR(t)Estimated BP (t) are compared between the standard and MC methods in Figure 4.47. The corre-sponding DAR(t) are presented in Figure 4.48.For both BP (t) and DAR(t), in the high release clusters (C1, C3, and C5), there is littledifference in the estimates. Compared to SM-OSEM, these BP (t) estimates still remain relativelyunbiased, and in DAR(t) the peak magnitudes are slightly less negatively biased.In low release clusters (C2, C4, and C6), BP (t) estimates in the MC method are relativelyunchanged from SM-HYPR-AU. DAR(t) estimates, however, are not as positively biased, againsuggesting that the improved sensitivity reduces the selection effect of voxels with fits having higherpeak DAR(t).Note that in the random release magnitude simulation we found that overall the BP (t) curveswere shifted downwards, which suggested that the edge voxels of the striatum (having lower pre-intervention binding as well as lower release levels, both from PSF filtering) made up the bulk of theadditionally detected voxels in the MC methods. Further evidence for this came from the voxelwise1084.4. Monte Carlo method(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.47: Comparing estimated BP (t) in the cluster release simulation between the standard and MCmethods. Solid lines indicate the median estimate. Shaded bars indicate ± one median absolute deviation.The ground truth is presented in black.1094.4. Monte Carlo method(a) C1 (b) C2(c) C3 (d) C4(e) C5 (f) C6Figure 4.48: Comparing estimated DAR(t) in the cluster release simulation between the standard and MCmethods. Solid lines indicate the median estimate. Shaded bars indicate ± one median absolute deviation.The ground truth is presented in black.1104.5. Parametric thresholdingsensitivity images, where large improvements in voxelwise sensitivity near the edge of the striatumwere visualized. In the cluster relase simulation, since the release is confined to clusters locatedmore centrally in the striatum, intra-striatal contamination is more prevalent than extra-stiatalcontamination. Thus the lower binding edge voxels of the striatum mostly do not contain release.As a result, the additionally detected voxels in the MC method tend to be near the edges of thecluster, which only have lower ground truth release levels—instead of lower pre-intervention bindingas well—from PSF filtering. Ultimately, this means BP (t) estimates in a scenario presented in thecluster release simulation are already fairly accurate with SM-HYPR-AU, and the additional voxelsdetected by the MC methods do not modify the estimates significantly—for better or for worse.4.5 Parametric thresholdingThe idea of the parametric thresholding approach (PTA) is use only a subset of detected clustervoxels—those with suspected higher parametric accuracy—to provide a more accurate cluster-levelrepresentation of BP (t) and DAR(t).From Section 3.5.8, PTA computes the absolute deviation (AD) of the estimated tD and tP ineach detected voxel from the corresponding estimates in its up to 124 nearest detected neighboursin 3D. From here, the parametric AD(tD) and AD(tP ) images are thresholded by some value, andonly voxels who pass both thresholds are used to calculate the median BP (t) and DAR(t) estimatesin the cluster.A proof [equation (3.15)] demonstrated that AD is bounded by the local MAE in the timingparameter, plus the error in the voxel-level estimate. Thus parametrically “average” voxels havingerrors equal to the MAE are bounded by twice MAE.From the cluster release simulation analysis of the timing parameters, we found that in theleast accurate cluster (C4: having both low release and located in caudate where PVE is moreprominent) had MAE(tD) ∼ 4 min and MAE(tP ) ∼ 6 min, whereas the better performing clustersin terms of parametric accuracy had MAE(tD) ∼ 3.5 min and MAE(tP ) ∼ 5 min. Therefore if weset AD(tD) ≤ 4 min and AD(tP ) ≤ 6 as our parametric thresholds, for the least parametrically-accurate clusters, in the worst case scenario where the equality of equation (3.15) holds, only avoxel with zero error will be accepted. However, in more parametrically-accurate clusters, andcases where the inequality holds, then voxels with some—but still relatively low—parametric errorwill pass through the thresholds.Using these thresholds, we first perform a sanity check by observing the tD and tP distributionsfor the voxels surviving both thresholds, presented in Figure 4.49.Compared to the distributions without PTA (Figure 4.45), for both tD and tP the estimates arebetter centered about the ground truth, and the variance in the voxel-level estimates is lower—asshould be expected from our MAE-based thresholding. Particularly in the distribution of estimatedtD, the outlier fits with tD = 30 are much reduced after PTA.Next, we recompute estimated BP (t) and DAR(t) with the subset of voxels passing through thethresholds for each cluster of the cluster release simulation, using both the standard and MC meth-ods. The new BP (t) estimates are presented in Figure 4.50 (high release clusters) and Figure 4.51(low release clusters). The corresponding DAR(t) are presented in Figure 4.52 and Figure 4.53,respectively.For all clusters, the shape of both BP (t) and DAR(t) estimates better conform to the groundtruth after PTA. This is a result of lower MAE in the timing parameters, such that the time ofrelease initiation, peak, and decay back to baseline is more representative of the ground truth.On the whole, this results in more accurate magnitude estimates as well when we compare1114.5. Parametric thresholding(a) no PTA (b) with PTAFigure 4.49: Comparing distributions of estimated tD and tP in the cluster release simulation, betweenthe standard and MC methods, before (a) and after (b) parametric thresholding has been applied.across all time points. In some cases, such as the estimates of BP (t) and DAR(t) in C4, the peakmagnitude has greater bias after PTA, though the bias is reduced before the intervention and oncethe release begins to decay to baseline levels. This may be the result of yet another selection effect,where the voxels with higher release levels also have better timing parameter estimates (since theDA-induced perturbation makes up more of the signal). In all clusters, a distribution of releaseexists from the initial PSF filtering, so even in low release clusters there exists a small subset ofrelatively high release voxels (c.f. Figure 3.4). In PTA the voxels with lower parametric accuracyare not included in the computation of BP (t) and DAR(t), so the resulting estimates would bebiased towards these voxels with higher release. Since C4 has low sensitivity, low release, andextensive PVE, it is prone to having many of its originally-detected voxels not passing through theparametric thresholds. Thus if the above reasoning is correct, the surviving voxels will skewed tothe highest release voxels within the cluster.These improvements occurs for all detection methods, though in some instances for SM-OSEM,such as the estimates of DAR(t) in C3 and C6, even after PTA the tail end of the release remainsseverely overestimated. This may be because the images are of poor parametric accuracy to beginwith, so PTA can only compute the new estimate from a more parametrically accurate—albeit stillseverely biased—subset of voxels.These results demonstrate that it is useful to separate the analysis of where within the brain theclusters exists, from what the release dynamics are within the cluster; the final binary significancemask from one of the MC methods provides the former, and the cluster-level median BP (t) andDAR(t) estimates after PTA provide the latter.1124.5. Parametric thresholding(a) C1: no PTA (b) C1: with PTA(c) C3: no PTA (d) C3: with PTA(e) C5: no PTA (f) C5: with PTAFigure 4.50: Comparing estimated BP (t) in the cluster release simulation between the standard and MCmethods, with and without using the parametric thresholding approach, for high release clusters. Solid linesindicate the median estimate. Shaded bars indicate ± one median absolute deviation. The ground truth ispresented in black.1134.5. Parametric thresholding(a) C2: no PTA (b) C2: with PTA(c) C4: no PTA (d) C4: with PTA(e) C6: no PTA (f) C6: with PTAFigure 4.51: Comparing estimated BP (t) in the cluster release simulation between the standard and MCmethods, with and without using the parametric thresholding approach, for low release clusters. Solid linesindicate the median estimate. Shaded bars indicate ± one median absolute deviation. The ground truth ispresented in black.1144.5. Parametric thresholding(a) C1: no PTA (b) C1: with PTA(c) C3: no PTA (d) C3: with PTA(e) C5: no PTA (f) C5: with PTAFigure 4.52: Comparing estimated DAR(t) in the cluster release simulation between the standard and MCmethods, with and without using the parametric thresholding approach, for high release clusters. Solid linesindicate the median estimate. Shaded bars indicate ± one median absolute deviation. The ground truth ispresented in black.1154.5. Parametric thresholding(a) C2: no PTA (b) C2: with PTA(c) C4: no PTA (d) C4: with PTA(e) C6: no PTA (f) C6: with PTAFigure 4.53: Comparing estimated DAR(t) in the cluster release simulation between the standard and MCmethods, with and without using the parametric thresholding approach, for low release clusters. Solid linesindicate the median estimate. Shaded bars indicate ± one median absolute deviation. The ground truth ispresented in black.1164.6. Application to a human pilot study4.6 Application to a human pilot studyThis work so far has shown that a MC-based method applied to the parametrically accurate recon-struction HYPR-AU-OSEM optimizes the sensitivity versus CST trade-off when detecting clustersof DA release in response to a task or stimulus. The MC methods will provide a binary significancemask that indicates the location of such a cluster within the brain, and after a parametric threshold-ing of the voxel-level parameter estimates within the cluster, estimates of the time-varying bindingpotential, BP (t), and DA release level, DAR(t) can be provided for each cluster that reasonablyapproximate the ground truth.We are now able to apply this optimized detection and analysis suite to a real human study,where the ground truth is not known. As a reminder, we have scanned a single participant usinga 20 mCi RAC injection for 75 mins in the HRRT. At t = 36 mins, the participant was subjectedto a gambling task that lasted ∼ 15 mins.We apply both optimizations of the MC method (MC-CST and MC-SENS) to the dataset, fol-lowed by PTA to arrive at BP (t) and DAR(t) estimates for each detected cluster. For comparison,the two representative standard methods (SM-HYPR-AU and SM-OSEM) are also applied to thedataset. In Figure 4.54, we present the collective results from these analyses.MC-CST is able to detect seven parametrically-unique clusters, whereas MC-SENS merges twoof these clusters together for a total of six detected clusters, though the regional extent of therelease detected by MC-SENS largely overlaps with that of MC-CST. The separation of these twoclusters in MC-CST occurred during PTA; the final binary significance mask had the two clustersmerged together in both methods. As well, PTA was responsible for separating clusters 3 and 5 inMC-CST and MC-SENS (Figure 4.55); these two clusters clearly have unique release dynamics, butwere located in close proximity to each other, thus becoming merged during the initial detectionprocess.SM-OSEM detects five parametrically-unqiue clusters. Relating these to the seven of MC-CST,clusters 1 and 7 are merged together and cluster 5 does not pass through the parametric threshold.Again, the regional extent of release detected by SM-OSEM is very similar to that of the MCmethods, as was observed in the voxelwise sensitivity images from our simulations (Figure 4.43).SM-HYPR-AU only detects four of the seven clusters from MC-CST. This is again consistentwith our simulations: the missing clusters are likely a reflection of the lower sensitivity of thismethod. However, since the CST thresholding is only designed to produce a false postive-free setof clusters 90% of the time, so there is a chance that some of the clusters detected here are falsepositives. Nevertheless, in our simulations the methods commonly had the same noisy realizationswhere false positive clusters exceeded the size of the CST, which would suggest that any falsepositive in this real human scan would be common amongst the four methods tested.In terms of release dynamics, the estimated DAR(t) have similar peak magnitudes to oursimulations (40%-60%), as well as similar timing characteristics: most clusters have release that isshort-lived, quickly reaching its peak magnitude before decaying back to baseline before the end ofthe scan. One difference from the simulations is the existence of clusters with initiation of releasewell after the intervention: clusters 5 and 6 both initiate release at tD ∼ 50 mins and peak ∼ 5mins later. This may suggest a secondary release associated with the completion of the task, whichalso occurred at t ∼ 50 minFor clusters 1-4, the release dynamics are similar, with release initiating just before the inter-vention and peaking at or just after the time of intervention. This may indicate an anticipatoryeffect. In cluster 7, release starts around the time of intervention and peaks ∼ 10 mins later.Comparing the dynamics across the different methods, we see a trend similar to our simulations,where the DAR(t) estimates for SM-OSEM have lower peak magnitudes, and an extended tail.1174.6. Application to a human pilot studyFigure 4.54: A summary of the results from applying all methods with PTA to a human pilot scanwith a gambling task presented at 36 mins into a 75 min RAC scan. Identified clusters are colour-coded,whereas voxels identified as non-release are coloured in dark gray. Voxels that were detected but failedparametric thresholding are coloured in light gray. DAR(t) for each method is presented for each cluster,when applicable. An insert of a representative coronal slice showcasing the cluster in included in each plot.1184.6. Application to a human pilot studyFigure 4.55: Demonstrating the effect of PTA on real human data: clusters 3 and 5 are parametricallyunique, but were initially merged together in the initial detection process, only to be separated by PTA.Shown here are the DAR(t) estimates from MC-CST, before and after PTA.DAR(t) estimates between SM-HYPR-AU and the MC methods are quite similar.The results from this single human scan demonstrate that the MC method, as well as PTA,produce results largely consistent with that of our simulations. This supports not only the validityof the novel methods of this work, but also our simulation methods: both in a physiological sensewhen designing our noiseless datasets, and in the simulation of the imaging process.119Chapter 5ConclusionsThe purpose of this thesis has been twofold:1. We first presented a comprehensive review on the effect of post-processing on the detectionability of transient DA release in PET RAC imaging. In the standard method of analysis,voxelwise fits using SRTM (unable to model a transient response to a task or stimulus) andlp-ntPET (able to model such a response) are statistically compared using an F-test. Falsepositives are controlled by first thresholding at an Fcrit corresponding to p < 0.05, followedby a regional cluster size threshold designed to ensure 90% of baseline scans are completelyfalse positive-free. Detection sensitivity (i.e. the true positive rate) was shown to improvewith more extensive denoising, but at the cost of greater CST; uncovering a fundamentaltrade-off in the standard method.2. This motivated a novel Monte Carlo-based method that replaces the F statistic with a con-fidence metric, c, derived from several iterations of a Metropolis algorithm. From here,significance thresholding at c > ccrit followed by cluster size thresholding proceeds as before.By designing the cost function within the Metropolis algorithm to correct observed flaws inthe standard method, the MC method is able circumvent the above trade-off and improvedetection sensitivity while not increasing CST.In what follows, we break down these points in greater detail.5.1 Improved noise-modeling over previous workThese results were obtained through realistic simulations based on real human RAC datasets. Usinga RAC scan of a healthy control in the baseline state, we divided the brain into anatomical ROIsusing a segmented MRI image, and performed regional fits using SRTM. These fits then served asthe noiseless TACs for our simulated baseline image. Release was incorporated into any given TACusing lp-ntPET.Previous work [57, 60] has modeled the imaging process by simply adding Gaussian noise, scaledto level observed in real human datasets. We instead elect for a more realistic simulation of theimaging process, first incorporating the effect of the PSF of the HRRT scanner by filtering ournoiseless images. From here, the images are forward-projection to sinogram space, where Poissonnoise is added and the effects of attenuation and detector normalization are accounted for, usingthe attenuation sinogram of the healthy control and the normalization sinogram of the scanner.The noisy sinograms could then be reconstructed and post-processed as necessary.Compared to the Gaussian noise addition method, we found that our more realistic simula-tions tended to achieve lower sensitivity and higher CST. This may be because the noise associ-ated with each LOR generates spatially-correlated noise patterns, whereas adding noise throughspatially-independent Gaussian samples does not—effectively enlarging false positive regions andcomparatively mitigating the effect of denoising when trying to detect true release. Alternatively,1205.2. Variables in simulation design, image processing, and analysisthe method of scaling Gaussian noise to scanner levels may not be robust, potentially leading tolower noise levels in the TACs and thus making false positive removal and voxelwise DA releasedetection easier.Using our more realistic noise-addition methods, we generated 50 realizations of three uniquesimulation types: baseline, random release magnitude, and cluster release. Baseline had no releaseincorporated in order to determine the CST. random release magnitude had spatially-invariantrelease levels introduced into all subregions of the striatum, save for the R putamen. The clusterrelease simulation had six clusters incorporated throughout the striatum, with sizes of either 300,400, or 500 voxels, and mean release levels of either DAR = 7.5% or DAR = 15%.5.2 Variables in simulation design, image processing, andanalysisEach simulation was generated at three dose levels to investigate the effect of dose level. The simula-tions were reconstructed using either OSEM or the denoised reconstruction HYPR-AU-OSEM (2.5mm Gaussian kernel). The OSEM images were either not post-filtered, or post-filtered with a 2 mmor 5 mm Gaussian kernel. In all cases (the HYPR-AU-OSEM images, the unfiltered OSEM images,and the filtered OSEM images) either 2.5 mm or 5 mm HYPR post-processing was subsequentlyapplied.When fitting the voxel-level TACs of the reconstructed and post-processed images, using SRTMand lp-ntPET, three weighting schemes were investigated: a traditional weighting that accountsfor decay correction and assumes the noise level scales as the Poisson noise level associated withthe number of counts within a frame, a uniform weighting scheme (i.e. ordinary least squaresfitting), and an uncorrected weighting where the decay correction is undone, such that ordinaryleast squares is performed on the decay-uncorrected TACs.5.3 Baseline simulation resultsThe baseline simulations yielded the CST for each unique triad of dose level, reconstruction andpost-processing protocol, and weighting scheme, by determining the cluster size that removed allfalse positives from 90% of all simulations.We found that CST is dramatically increased when using a 5 mm kernel width instead of2.5 mm in the HYPR post-processing operator. As well, wider kernels used for post-filtering theOSEM images also increases the CST. Comparing reconstruction type, we found that HYPR-AU-OSEM provides a smaller CST than OSEM, since the denoised reconstruction removes noise whilepreserving high frequency features.Dose level often reduces CST, but this is not consistent for all triads, suggesting that more than50 realizations may be needed to generate a more statistically robust estimate of the false positivecluster size distribution. Thus we took the maximum CST computed across all dose levels as aconservative measure. In terms of weighting scheme, uniform weighting consistently produced thehighest CST whereas uncorrected weighting yields the lowest.5.4 Random release magnitude simulation resultsThe random release magnitude simulation allows us to characterize the detection sensitivity as afunction of release level, without the CST becoming important since the region of release is an1215.5. Cluster release simulation resultsorder of magnitude larger than most CSTs.We found that a 5 mm kernel width is required (either in the spatial post-filtering operation, orHYPR post-processing, or both) to have useful detection sensitivity (i.e. sensitivity above 50% formost release levels). Predictably, as release level increases, detection sensitivity increases as well.When OSEM is post-filtered with a 5 mm Gaussian kernel and then 5 mm HYPR post-processing(OSEM-5mm-5mm), the detection sensitivity is noticeably improved over all other reconstructionand denoising protocols. HYPR-AU-OSEM images with 5 mm HYPR post-processing (HYPR-2.5mm-5mm) performed similarly to OSEM images with no post-filtering and 5 mm HYPR post-processing (OSEM-NF-5mm) and with 2 mm post-filtering plus 5 mm HYPR post-processing(OSEM-2mm-5mm).As dose level increased, modest increases in sensitivity across all release levels were observedmoving from 16 mCi to 20 mCi, but the improvement from 20 mCi to 24 mCi was not as significant.Thus we concluded that using 20 mCi is the optimal dose to use when also weighing economicaland ethical considerations.In terms of weighting schemes, traditional weighting produced the highest sensitivity, followedby uniform weighting and then uncorrected weighting. Because uniform weighting produced neitherthe highest sensitivity nor lowest CST, we conclude that it is not an optimal weighting scheme touse.When analyzing parametric accuracy, provided relatively large regions of release exist withinthe striatum, we found that the mean of the voxel-level timing parameter estimates (tD: whenthe release initiates, and tP : when the release level peaks) are accurate to within a minute or lessacross all release levels, whereas the voxel-level estimates have mean errors of MAE(tD) ∼ 3 minand MAE(tP ) ∼ 5.In general, OSEM-5mm-5mm produces the lowest MAE, and for all reconstruction and denoisingprotocols MAE decreases with increasing release level. However, when analyzing the time-varyingbinding potentials, BP (t), and DA release levels, DAR(t), we found that OSEM-5mm-5mm severelyunderestimates BP (t) and has slightly more negative biases in DAR(t) than the other methods.Since these metrics offer a more complete picture of the evolution of the DA release in response tothe task or intervention, we conclude that OSEM-5mm-5mm is the least useful in terms of para-metric accuracy. In constrast, the resolution-preserving qualities of HYPR-2.5mm-5mm providethe highest parametric accuracy when studying BP (t) and DAR(t).5.5 Cluster release simulation resultsThe cluster release simulation renders the CST more important, since the cluster size is now com-parable to the CST for some methods.We found that OSEM-5mm-5mm still has higher voxelwise sensitivity, but now performs simi-larly or worse than the other methods in binary sensitivity: the percentage of time any portion ofa true release cluster is detected. This occurs because i) in high release clusters, all methods areable to detect at least some fraction of true release voxels essentially 100% of the time; and ii) inlow release clusters, it is likely that only a fraction of the true release voxels are detected, so thedetected portion of the cluster now commonly falls below the CST of OSEM-5mm-5mm.For a given release level, however, the size of the cluster did not appreciably change its sensitivity.Rather, the location of the cluster had a larger impact; clusters in the caudate consistently had lowersensitivity than those with identical release levels and similar cluster sizes in the putamen. Thisstems from the voxel-level TACs—especially those near the edges of the cluster—having greaterPVE as a consequence of the smaller size and elongated shape of the caudate.1225.6. Take-aways from analyzing the simulations using the standard methodIn terms of parametric accuracy, the results largely agreed with those from the random releasemagnitude simulation. We also discovered a potential selection effect in lower sensitivity clusters,where voxels having fits with higher release magnitudes (either due to noise, bias, or simply havinghigher release than most other voxels within the cluster) are easier to detect, hence positivelybiasing their estimates of DAR(t) in low sensitivity clusters.5.6 Take-aways from analyzing the simulations using thestandard methodCollectively, the simulations point to an intrinsic trade-off when using the standard method: one hasto pay for increased voxelwise sensitivity with increased CST. OSEM-5mm-5mm with traditionalweighting produces the highest voxelwise sensitivity (an average of 76.5% across all release levelsin the random release magnitude simulation) at a CST of 370 voxels. On the other hand, HYPR-2.5mm-5mm with uncorrected weighting produces the lowest CST of 108 voxels, but with an averagevoxelwise sensitivity of 52.1% in the random release magnitude simulation.As well, while voxelwise sensitivity is useful for determining the full regional extent of theDA release, binary sensitivity can be considered a more important metric in the sense that a realapplication of this method only is comprised of a single scan per patient: for harder-to-detect lowrelease clusters, the higher voxelwise sensitivity of OSEM-5mm-5mm indicates a larger poriton ofthe true cluster is detected—when this portion passes through the CST—but its inferior binarysensitivity means that it has a greater likelihood of completely missing the cluster in a single scan.Furthermore, the smallest cluster we have simulated is 300 voxels (0.64 cm3), but one in theorycould imagine smaller clusters existing in response to a task or stimulus. Thus we left this reviewof the standard method seeking an improved method that is able to combine the high sensitivitycharacteristics of OSEM-5mm-5mm with the low CST and superior parametric accuracy of HYPR-2.5mm-5mm.5.7 Monte Carlo methodSince the MC method only modifies the significance thresholding step of the detection algorithm anddoes not modify the denoised images, we decided to apply the MC method to HYPR-2.5mm-5mmto maximize parametric accuracy. As well, to minimize the CST we used uncorrected weighting,which reduces the CST from traditional weighting by 23%, at the cost of only a 11% reduction inaverage voxelwise sensitivity.The goal of the MC method is compete with OSEM-5mm-5mm in terms of sensitivity, with mucha much smaller CST. Thus we optimized the MC method (through choices of the sole hyperparam-eter in the cost function, β, and the significance thresholding value, ccrit) in two ways: MC-CSTmaximizes sensitivity without increasing the sensitivity of the best-performing standard method interms of CST (HYPR-2.5mm-5mm with uncorrected weighting, hereafter SM-HYPR-AU), and MC-SENS maximizes sensitivity without exceeding the CST of the best-performing standard method interms of voxelwise sensitivity in the cluster release simulation (OSEM-5mm-5mm with uncorrectedweighting, hereafter SM-OSEM).In the random release magnitude simulation, we found that MC-CST performs similarly to SM-OSEM across all release levels, albeit with a 56% reduced CST and superior parametric accuracyfrom the HYPR-AU-OSEM images. MC-SENS improves sensitivity over SM-OSEM by up to 15%,with a 26% reduction in CST. Since the MC methods start from the same images as SM-HYPR-AU,it is also worth mentioning that MC-CST had 77% higher voxelwise sensitivity than SM-HYPR-AU1235.8. Parametric thresholdingat low release and 17% at high release, and MC-SENS improved improved voxelwise sensitivity by103% at low release and 21% at high release.The results from the cluster release simulation were largely similar. In low release clusters,MC-CST (MC-SENS) improved voxelwise sensitivity over SM-HYPR-AU by 43% (149%), and inhigh release clusters the improvements were 16% (38%). MC-CST was slightly outperformed bySM-OSEM in high release clusters, but in low release clusters MC-CST performed similarly orslightly better. MC-SENS offered improvements in both low release clusters (57% on average) andhigh release clusters (5% on average).The improvements when using the MC methods were generally higher in this simulation, sincethe true cluster sizes are on the order of CST, so initially-detected clusters are more frequentlyrejected by the CST in the standard method. However, the MC methods improve sensitivity fromSM-HYPR-AU while either maintaining CST, as in MC-CST, or maintaining relatively low CSTas compared to SM-OSEM, as in MC-SENS, so this rejection issue is less prevalent. This is bestencapsulated with the metric of binary sensitivity: in low release clusters of the putamen, SM-OSEM had an average binary sensitivity of 60%, which is improved to 71% by MC-CST and to82% by MC-SENS. In the sole low release cluster in the caudate, extensive PVE results in SM-OSEM having a binary sensitivity of 26%, whereas MC-CST improves this to 48%, and MC-SENSto 50%. Therefore, the MC methods have a higher likelihood of detecting at least a portion (i.e.higher binary sensitivity) of the difficult-to-detect lower release clusters, and when they do theyalso detect a similar or greater extent of the true release voxels within those clusters (i.e. similaror higher voxelwise sensitivity).The MC methods, by definition, also preserve voxel-level parametric accuracy since they do notmodify the original images. Thus by applying the MC methods to HYPR-AU-OSEM images—themost parametrically accurate images—we have designed a method with both superior detectionsensitivity and parametric accuracy than SM-OSEM, the current best-performing method in termsof detection sensitivity.However, since the set of voxels detected within a true release cluster is larger in the MC methodsthan SM-HYPR-AU, the cluster-level or average voxel-level parametric accuracy will differ; we haveshown that the additionally detected voxels by the MC methods lead to slightly higher MAE(tD)and slightly lower MAE(tP ). Additionally, in capturing the full release dynamics and magnitudethrough BP (t) and DAR(t), visual inspection of the estimated curves show little difference from theSM-HYPR-AU estimates, and large improvements over the less parametrically accurate SM-OSEMestimates.5.8 Parametric thresholdingBP (t) and DAR(t) estimates can be further improved upon by introducing a voxelwise parametricthresholding that only includes a subset of the more parametrically accurate voxels (in terms oftheir timing parameter estimates) from the set of all detected voxels within a cluster. Applyingthis thresholding approach, as expected the estimates show better agreement to the ground truthvalues across time, but most notably during the initial onset of release, and in the decay back tobaseline after peak release.5.9 Application to a human pilot studyThe theory and optimization of the MC method and PTA was tested on simulated data for themajority of this thesis, but we also applied the final optimized methods to a single human pilot1245.10. Final thoughtsscan with a gambling task as the mid-scan intervention.Since the ground truth is evidently not known, the conclusions drawn are necessarily speculativein nature. Nevertheless, consistent with our simulations we find that MC-CST, MC-SENS, andSM-OSEM detect a similar regional extent of DA release, though SM-OSEM has one of the sevendetected clusters fail to pass through parametric thresholding, and SM-HYPR-AU fails to detectthree of the seven detected clusters from the other methods. Also consistent with our simulatedresults, SM-OSEM DAR(t) estimates have lower peak magnitude and an extended tail as comparedto the other methods.Together, these results suggest that the conclusions from our simulations also apply to a realdataset: the MC methods offer similar or better performance to SM-OSEM, but with superiorparametric accuracy.5.10 Final thoughtsIn the difficult goal of detecting voxelwise DA release in the brain in response to a mid-scan taskor stimulus, we have first shown that extensive denoising is required to bring the voxel-level TACsto a sufficiently high SNR for distinctions between release TACs and non-release TACs to be made(through an F-test).The purpose of a voxelwise approach is to detect subtle neurophysiological changes on thesmallest possible volume—the voxel—yet this extensive denoising has the effect of increasing CST,thus reintroducing a regional constraint on detectability. Furthermore, the most effective denoisingprotocol with respect to this application also provides the poorest parametric accuracy in themetrics of interest in a DA release study—BP (t) and DAR(t).The MC method developed in the work, in conjunction with a parametric thresholding ap-proach, can then be applied to a denoising protocol with lower detection sensitivity but the highestparametric accuracy. When this is done, the detection sensitivity can match or exceed the lev-els of the best-performing standard method, while maintaining relatively low CST and reasonablyaccurate approximations of the ground truth BP (t) and DAR(t).The MC method can be likened to applying a spatial filter to the binary mask indicatingthe detected release voxels of the standard method, though the kernel would be unique for eachvoxel, taking into account both the local parity and global adherence to a correct F-distributionfor all baseline voxels; simply filtering the binary significance mask would also improve sensitivitybut simultaneously increase CST by isotropically increase the size of both true clusters and falsepositive clusters, whereas the MC method will only probabilistically increase cluster size if theglobal deviation from a true F-distribution suggests a subset of false negative voxels exists in thatregion.Through simulation and a pilot scan, we have shown that the MC method can be easily im-plemented into the standard analysis pipeline for voxelwise detection of transient DA release, andachieves the desired goal of being able to detect small regions of release—through its relativelysmall CST—with high probability (i.e. high binary and voxelwise sensitivity) and high paramet-ric accuracy. Based on these results, the method can now reliably be applied to full-scale humanstudies.125Bibliography[1] Bevington, C. W. J., Klyuzhin, I. S., Aceves, L., Doudet, D., Cheng, J.-C. K., and Sossi, V.Denoising and DA release: effect of denoising on the ability to identify voxel-level neurophys-iological response. In 2018 IEEE Medical Imaging Conference. Conference Record, 2018.[2] Bevington, C. W. J., Klyuzhin, I. S., Cheng, J.-C. K., and Sossi, V. A Monte Carlo approachfor improving transient dopamine release detection sensitivity. Journal of Cerebral Blood Flow& Metabolism, 0000. [in preparation].[3] Hooker, J. M. and Carson, R. E. Human positron emission tomography neuroimaging. AnnualReview of Biomedical Engineering, 21(1):551–581, 2019.[4] Barrio, J. R., Huang, S. C., Melega, W. P., Yu, D. C., Hoffman, J. M., Schneider, J. S.,Satyamurthy, N., Mazziotta, J. C., and Phelps, M. E. 6-[18F]Fluoro-L-DOPA probes dopamineturnover rates in central dopaminergic structures. Journal of Neuroscience Research, 27(4):487–493, 1990.[5] Logan, J., Dewey, S. L., Wolf, A. P., Fowler, J. S., Brodie, J. D., Angrist, B., Volkow, N. D.,and Gatley, S. J. Effects of endogenous dopamine on measures of [18F]N-methylspiroperidolbinding in the basal ganglia: Comparison of simulations and experimental results from PETstudies in baboons. Synapse, 9(3):195–207, 1991.[6] Gifford, A. N., Gatley, S. J., and Ashby Jr., C. R. Endogenously released dopamine inhibits thebinding of dopaminergic PET and SPECT ligands in superfused rat striatal slices. Synapse,22(3):232–238, 1996.[7] Barrett, S. P., Boileau, I., Okker, J., Pihl, R. O., and Dagher, A. The hedonic response tocigarette smoking is proportional to dopamine release in the human striatum as measured bypositron emission tomography and [11C]raclopride. Synapse, 54(2):65–71, 2004.[8] Trevor Young, L., Wong, D. F., Goldman, S., Minkin, E., Chen, C., Matsumura, K., Schef-fel, U., and Wagner Jr., H. N. Effects of endogenous dopamine on kinetics of [3H]N-methylspiperone and [3H]raclopride binding in the rat brain. Synapse, 9(3):188–194, 1991.[9] Cumming, P., Wong, D. F., Dannals, R. F., Gillings, N., Hilton, J., Scheffel, U., and Gjedde,A. The competition between endogenous dopamine and radioligands for specific binding todopamine receptors. Annals of the New York Academy of Sciences, 965(1):440–450, 2002.[10] Ruhe´, H., Mason, N., and Schene, A. Mood is indirectly related to serotonin, norepinephrineand dopamine levels in humans: A meta-analysis of monoamine depletion studies. MolecularPsychiatry, 12:331–359, 2007.[11] Meeusen, R. and Roelands, B. Fatigue: Is it all neurochemistry? European Journal of SportScience, 18(1):37–46, 2018.126Bibliography[12] Fetissov, S., Meguid, M., Chen, C., and Miyata, G. Synchronized release of dopamine andserotonin in the medial and lateral hypothalamus of rats. Neuroscience, 101(3):657–663, 2000.[13] Jarcho, J. M., Feier, N. A., Labus, J. S., Naliboff, B., Smith, S. R., Hong, J.-Y., Colloca, L.,Tillisch, K., Mandelkern, M. A., Mayer, E. A., and London, E. D. Placebo analgesia: Self-report measures and preliminary evidence of cortical dopamine release associated with placeboresponse. NeuroImage: Clinical, 10:107–114, 2016.[14] M Sullivan, J., Jin Kim, S., P Cosgrove, K., and D Morris, E. Limitations of SRTM, Logangraphical method, and equilibrium analysis for measuring transient dopamine release with[(11)C]raclopride PET. American Journal of Nuclear Medicine and Molecular Imaging, 3:247–260, 2013.[15] Alpert, N. M., Badgaiyan, R. D., Livni, E., and Fischman, A. J. A novel method for noninvasivedetection of neuromodulatory changes in specific neurotransmitter systems. NeuroImage, 19(3):1049–1060, 2003.[16] D. Morris, E., Yoder, K., Wang, C., D. Normandin, M., Zheng, Q.-H., Mock, B., F. Muzic Jr,R., and C. Froehlich, J. ntPET: A new application of PET imaging for characterizing thekinetics of endogenous neurotransmitter release. Molecular Imaging, 4:473–489, 2005.[17] Normandin, M. D., Schiffer, W. K., and Morris, E. D. A linear model for estimation ofneurotransmitter response profiles from dynamic PET data. NeuroImage, 59(3):2689–2699,2012.[18] Phelps, M. E. PET: Molecular Imaging and Its Biological Applications. Springer, 1 edition,2004.[19] Wienhard, K., Schmand, M., Casey, M., Baker, K., Bao, J., Eriksson, L., Jones, W., Knoess,C., Lenox, M., Lercher, M., Luk, P., Michel, C., Reed, J., Richerzhagen, N., Treffert, J.,Vollmar, S., Young, J., Heiss, W.-D., and Nutt, R. The ECAT HRRT: Performance and firstclinical application of the new high resolution research tomograph. IEEE Transactions onNuclear Science, 49:104–110, 2002.[20] Bauer, C. E., Brefczynski-Lewis, J., Marano, G., Mandich, M.-B., Stolin, A., Martone, P.,Lewis, J. W., Jaliparthi, G., Raylman, R. R., and Majewski, S. Concept of an upright wearablepositron emission tomography imager in humans. Brain and Behavior, 6(9):e00530, 2016.[21] Ahmed, A. M., Tashima, H., and Yamaya, T. Investigation of spatial resolution improve-ment by use of a mouth-insert detector in the helmet PET scanner. Radiological Physics andTechnology, 11(1):7–12, 2018.[22] Noble, R. Ambulatory micro dose positron emission tomography: Wearable PET scanner forneurological imaging [epub]. Journal of Nuclear Medicine Technology, 2019.[23] Mullani, N. A., Markham, J., and Ter-Pogossian, M. M. Feasibility of time-of-flight recon-struction in positron emission tomography. Journal of Nuclear Medicine, 21(11):1095–1097,1980.[24] Grant, A. M., Deller, T. W., Khalighi, M. M., Maramraju, S. H., Delso, G., and Levin, C. S.NEMA NU 2-2012 performance studies for the SiPM-based ToF-PET component of the GESIGNA PET/MR system. Medical Physics, 43(5):2334–2343, 2016.127Bibliography[25] Wagatsuma, K., Miwa, K., Sakata, M., Oda, K., Ono, H., Kameyama, M., Toyohara, J.,and Ishii, K. Comparison between new-generation SiPM-based and conventional PMT-basedTOF-PET/CT. Physica Medica, 42:203–210, 2017.[26] Melcher, C. L. Scintillation crystals for PET. Journal of Nuclear Medicine, 41(6):1051–1055,2000.[27] Peng, H. and Levin, C. S. Recent developments in PET instrumentation. Current Pharma-ceutical Biotechnology, 11(6):555–571, 2010.[28] Meikle, S. R. and Badawi, R. D. Quantitative Techniques in PET, pages 93–126. SpringerLondon, London, 2005.[29] Politte, D. G. and Snyder, D. L. Corrections for accidental coincidences and attenuation inmaximum-likelihood image reconstruction for positron-emission tomography. IEEE Transac-tions on Medical Imaging, 10(1):82–89, 1991.[30] Chesler, D. A. Three-dimensional activity distribution from multiple positron scintigraphs.Journal of Nuclear Medicine, 12:347–348, 1971.[31] Shepp, L. A. and Vardi, Y. Maximum likelihood reconstruction for emission tomography.IEEE Transactions on Medical Imaging, 1(2):113–122, 1982.[32] Lange, K. and Carson, R. EM reconstruction algorithms for emission and transmission tomog-raphy. Journal of Computer Assisted Tomography, 8:306–316, 1984.[33] Tong, S., Alessio, A. M., and Kinahan, P. E. Image reconstruction for PET/CT scanners: pastachievements and future challenges. Imaging in Medicine, 2(5):529–545, 2010.[34] Hudson, H. M. and Larkin, R. S. Accelerated image reconstruction using ordered subsets ofprojection data. IEEE Transactions on Medical Imaging, 13(4):601–609, 1994.[35] Cheng, J.-C. K., Matthews, J., Sossi, V., Anton-Rodriguez, J., Salomon, A., and Boellaard, R.Incorporating HYPR de-noising within iterative PET reconstruction (HYPR-OSEM). Physicsin Medicine & Biology, 62(16):6666–6687, 2017.[36] Rahmim, A., Qi, J., and Sossi, V. Resolution modeling in PET imaging: Theory, practice,benefits, and pitfalls. Medical Physics, 40(6):064301, 2013.[37] Rahmim, A., Cheng, J.-C., Blinder, S., Camborde, M.-L., and Sossi, V. Statistical dynamicimage reconstruction in state-of-the-art high-resolution PET. Physics in Medicine & Biology,50(20):4887–4912, 2005.[38] Levin, C. S., Tai, Y. ., Hoffman, E. J., Dahlbom, M., and Farquhar, T. H. Removal of theeffect of compton scattering in 3-D whole body positron emission tomography by Monte Carlo.In 1995 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, volume 2,pages 1050–1054, 1995.[39] Watson, C. C., Newport, D., and Casey, M. E. A Single Scatter Simulation Technique forScatter Correction in 3D PET, pages 255–268. Springer Netherlands, Dordrecht, 1996.[40] Watson, C. New, faster, image-based scatter correction for 3D PET. IEEE Transactions onNuclear Science, 47:1587–1594, 2000.128Bibliography[41] Cherry, S. R. and Huang, S.-C. Effects of scatter on model parameter estimates in 3d petstudies of the human brain. IEEE Transactions on Nuclear Science, 42(4):1174–1179, 1995.[42] Christian, B. T., Vandehey, N. T., Floberg, J. M., and Mistretta, C. A. Dynamic PETdenoising with HYPR processing. Journal of Nuclear Medicine, 51(7):1147–1154, 2010.[43] Mistretta, C. A., Wieben, O., Velikina, J., Block, W., Perry, J., Wu, Y., Johnson, K., andWu, Y. Highly constrained backprojection for time-resolved MRI. Magnetic Resonance inMedicine, 55(1):30–40, 2006.[44] Floberg, J. M., Mistretta, C. A., Weichert, J. P., Hall, L. T., Holden, J. E., and Christian,B. T. Improved kinetic analysis of dynamic PET data with optimized HYPR-LR. MedicalPhysics, 39(6):3319–3331, 2012.[45] Cheng, J. K., Matthews, J., Boellaard, R., and Sossi, V. A MR guided de-noising for PETusing IHYPR-LR. In 2017 IEEE Nuclear Science Symposium and Medical Imaging Conference(NSS/MIC), pages 1–3, 2017.[46] Carson, R. E. Tracer Kinetic Modeling in PET, pages 127–159. Springer London, London,2005.[47] Elsinga, P. H., Hatano, K., and Ishiwata, K. PET tracers for imaging of the dopaminergicsystem. Current Medicinal Chemistry, 13(18):2139–2153, 2006.[48] Kanthan, M., Cumming, P., Hooker, J. M., and Vasdev, N. Classics in neuroimaging: Imagingthe dopaminergic pathway with PET. ACS Chemical Neuroscience, 8(9):1817–1819, 2017.[49] Innis, R. B., Cunningham, V. J., Delforge, J., Fujita, M., Gjedde, A., Gunn, R. N., Holden,J., Houle, S., Huang, S.-C., Ichise, M., Iida, H., Ito, H., Kimura, Y., Koeppe, R. A., Knudsen,G. M., Knuuti, J., Lammertsma, A. A., Laruelle, M., Logan, J., Maguire, R. P., Mintun,M. A., Morris, E. D., Parsey, R., Price, J. C., Slifstein, M., Sossi, V., Suhara, T., Votaw, J. R.,Wong, D. F., and Carson, R. E. Consensus nomenclature for in vivo imaging of reversiblybinding radioligands. Journal of Cerebral Blood Flow & Metabolism, 27(9):1533–1539, 2007.[50] Gunn, R. N., Gunn, S. R., and Cunningham, V. J. Positron emission tomography compart-mental models. Journal of Cerebral Blood Flow & Metabolism, 21(6):635–652, 2001.[51] Lammertsma, A. A. and Hume, S. P. Simplified reference tissue model for PET receptorstudies. NeuroImage, 4(3):153–158, 1996.[52] Gunn, R. N., Lammertsma, A. A., Hume, S. P., and Cunningham, V. J. Parametric imagingof ligand-receptor binding in PET using a simplified reference region model. NeuroImage, 6(4):279–287, 1997.[53] Van Laeken, N., Taylor, O., Polis, I., Neyt, S., Kersemans, K., Dobbeleir, A., Saunders, J.,Goethals, I., Peremans, K., and De Vos, F. In vivo evaluation of blood based and referencetissue based PET quantifications of [11C]DASB in the canine brain. PLOS ONE, 11(2):1–12,2016.[54] Floel, A., Garraux, G., Xu, B., Breitenstein, C., Knecht, S., Herscovitch, P., and Cohen, L.Levodopa increases memory encoding and dopamine release in the striatum in the elderly.Neurobiology of Aging, 29(2):267–279, 2008.129Bibliography[55] Sossi, V., Dinelle, K., Schulzer, M., Mak, E., Doudet, D. J., and de la Fuente-Ferna´ndez,R. Levodopa and pramipexole effects on presynaptic dopamine PET markers and estimateddopamine release. European Journal of Nuclear Medicine and Molecular Imaging, 37(12):2364–2370, 2010.[56] Black, K., Piccirillo, M., Koller, J., Hseih, T., Wang, L., and Mintun, M. Levodopa effects on[11C]raclopride binding in the resting human brain. F1000Research, 4(23), 2015.[57] Kim, S. J., Sullivan, J. M., Wang, S., Cosgrove, K. P., and Morris, E. D. Voxelwise lp-ntPETfor detecting localized, transient dopamine release of unknown timing: Sensitivity analysis andapplication to cigarette smoking in the PET scanner. Human Brain Mapping, 35(9):4876–4891,2014.[58] Me´rida, I., Olivier, F., Nzali, G. S., Hammers, A., Redoute´, J., Reilhac, A., Costes, N., andIrace, Z. Kinetic modelling for endogenous neurotransmitter discharge characterization usingPET imaging: optimization of lp-ntPET. Presented at NRM18, London, 2018.[59] Holmes, A. P., Blair, R. C., Watson, J. D. G., and Ford, I. Nonparametric analysis of statisticimages from functional mapping experiments. Journal of Cerebral Blood Flow & Metabolism,16(1):7–22, 1996.[60] Wang, S., Kim, S., Cosgrove, K. P., and Morris, E. D. A framework for designing dynamiclp-ntPET studies to maximize the sensitivity to transient neurotransmitter responses to drugs:Application to dopamine and smoking. NeuroImage, 146:701–714, 2017.[61] Lenz, W. Beitra¨ge zum versta¨ndnis der magnetischen eigenschaften in festen ko¨rpern. (Ger-man) [Contributions to the understanding of magnetic properties in solid bodies]. PhysikalischeZeitschrift, 21:613–615, 1920.[62] Newman, M. E. J. and Barkema, G. T. Ising Model and the Metropolis Algorithm, pages 43–84.Oxford University Press, Oxford, 1999.[63] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. Equationof state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087–1092, 1953.[64] Comtat, C., Bataille, F., Michel, C., Jones, J., Sibomana, M., Janeiro, L., and Tre´bossen,R. OSEM-3D reconstruction strategies for the ECAT HRRT. In 2004 IEEE Nuclear ScienceSymposium and Medical Imaging Conference (NSS/MIC), volume 6, pages 3492–3496, 2004.[65] Fischl, B., Salat, D. H., Busa, E., Albert, M., Dieterich, M., Haselgrove, C., van der Kouwe,A., Killiany, R., Kennedy, D., Klaveness, S., Montillo, A., Makris, N., Rosen, B., and Dale,A. M. Whole brain segmentation: Automated labeling of neuroanatomical structures in thehuman brain. Neuron, 33(3):341–355, 2002.[66] Fischl, B., Salat, D. H., van der Kouwe, A. J., Makris, N., Se´gonne, F., Quinn, B. T., andDale, A. M. Sequence-independent segmentation of magnetic resonance images. NeuroImage,23:S69–S84, 2004.[67] Han, X. and Fischl, B. Atlas renormalization for improved brain MR image segmentationacross scanner platforms. IEEE Transactions on Medical Imaging, 26(4):479–486, 2007.130Bibliography[68] Se´gonne, F., Dale, A., Busa, E., Glessner, M., Salat, D., Hahn, H., and Fischl, B. A hybridapproach to the skull stripping problem in MRI. NeuroImage, 22(3):1060–1075, 2004.[69] Carboni, E., Spielewoy, C., Vacca, C., Nosten-Bertrand, M., Giros, B., and Di Chiara, G.Cocaine and amphetamine increase extracellular dopamine in the nucleus accumbens of micelacking the dopamine transporter gene. Journal of Neuroscience, 21(9):RC141, 2001.[70] Pontieri, F., Tanda, G., Orzi, F., and Chiara, G. Effects of nicotine on the nucleus accumbensand similarity to those of addictive drugs. Nature, 382:255–257, 1996.[71] Domino, E. F. and Tsukada, H. Nicotine sensitization of monkey striatal dopamine release.European Journal of Pharmacology, 607(1):91–95, 2009.[72] Schiffer, W. K., Alexoff, D. L., Shea, C., Logan, J., and Dewey, S. L. Development of asimultaneous PET/microdialysis method to identify the optimal dose of 11C-raclopride forsmall animal imaging. Journal of Neuroscience Methods, 144(1):25–34, 2005.[73] Steeves, T. D. L., Miyasaki, J., Zurowski, M., Lang, A. E., Pellecchia, G., Van Eimeren, T.,Rusjan, P., Houle, S., and Strafella, A. P. Increased striatal dopamine release in Parkinsonianpatients with pathological gambling: a [11C] raclopride PET study. Brain, 132(5):1376–1385,2009.[74] de la Fuente-Fernndez, R., Phillips, A. G., Zamburlini, M., Sossi, V., Calne, D. B., Ruth, T. J.,and Stoessl, A. J. Dopamine release in human ventral striatum and expectation of reward.Behavioural Brain Research, 136(2):359–363, 2002.[75] Pappata, S., Dehaene, S., Poline, J., Gregoire, M., Jobert, A., Delforge, J., Frouin, V., Bot-tlaender, M., Dolle, F., Giamberardino, L. D., and Syrota, A. In vivo detection of striataldopamine release during reward: A PET study with [11C]Raclopride and a single dynamicscan approach. NeuroImage, 16(4):1015–1027, 2002.[76] Haruno, M. and Kawato, M. Different neural correlates of reward expectation and reward ex-pectation error in the putamen and caudate nucleus during stimulus-action-reward associationlearning. Journal of Neurophysiology, 95(2):948–959, 2006.[77] Ba¨ckman, L., Nyberg, L., Soveri, A., Johansson, J., Andersson, M., Dahlin, E., Neely, A. S.,Virta, J., Laine, M., and Rinne, J. O. Effects of working-memory training on striatal dopaminerelease. Science, 333(6043):718–718, 2011.[78] Choi, E. Y., Yeo, B. T. T., and Buckner, R. L. The organization of the human striatumestimated by intrinsic functional connectivity. Journal of Neurophysiology, 108(8):2242–2263,2012.[79] Pauli, W. M., O’Reilly, R. C., Yarkoni, T., and Wager, T. D. Regional specialization withinthe human striatum for diverse psychological functions. Proceedings of the National Academyof Sciences, 113(7):1907–1912, 2016.[80] Barrett, H. H., Wilson, D. W., and Tsui, B. M. W. Noise properties of the EM algorithm. I.Theory. Physics in Medicine & Biology, 39(5):833–846, 1994.[81] Teymurazyan, A., Riauka, T., Jans, H.-S., and Robinson, D. Properties of noise in positronemission tomography images reconstructed with filtered-backprojection and row-action maxi-mum likelihood algorithm. Journal of Digital Imaging, 26(3):447–456, 2013.131
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Designing an optimized protocol for detecting transient...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Designing an optimized protocol for detecting transient dopamine release Bevington, Connor William James 2019
pdf
Page Metadata
Item Metadata
Title | Designing an optimized protocol for detecting transient dopamine release |
Creator |
Bevington, Connor William James |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | A single positron emission tomography (PET) scan can be used to detect voxel-level transient dopamine release in response to mid-scan task or stimulus. A standard approach uses F-test (significance) thresholding followed by a cluster size threshold to limit false positives. The F statistic is computed from two neurophysiological kinetic models, one able to handle dopamine-induced perturbations to the time-varying radiotracer concentration, and the other unable to do so. Through simulation we first demonstrate that extensive denoising of the dynamic PET images is required for this method to have high detection sensitivity, though this often leads to a large cluster size threshold—limiting the detection of smaller cluster sizes—and poorer parametric accuracy. Our simulations also show that voxels near the peripheries of clusters are often rejected—becoming false negatives and ultimately distorting the F-distribution of rejected voxels. We then suggest a novel Monte Carlo method that incorporates these two observations into a cost function, allowing erroneously-rejected voxels to be accepted under specified criteria. In simulations, the proposed method boosts sensitivity by up to 77% while preserving cluster size threshold, or up to 180% when optimizing for detection sensitivity. A further parametric-based voxelwise thresholding is then suggested to better estimate the dopamine release dynamics in detected clusters. Finally, we apply the Monte Carlo method coupled with the parametric thresholding approach to a pilot scan from a human gambling study, where additional parametrically-unique clusters are detected as compared to the current best method. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-08-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NoDerivatives 4.0 International |
DOI | 10.14288/1.0380541 |
URI | http://hdl.handle.net/2429/71385 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_september_bevington_connor.pdf [ 15.89MB ]
- Metadata
- JSON: 24-1.0380541.json
- JSON-LD: 24-1.0380541-ld.json
- RDF/XML (Pretty): 24-1.0380541-rdf.xml
- RDF/JSON: 24-1.0380541-rdf.json
- Turtle: 24-1.0380541-turtle.txt
- N-Triples: 24-1.0380541-rdf-ntriples.txt
- Original Record: 24-1.0380541-source.json
- Full Text
- 24-1.0380541-fulltext.txt
- Citation
- 24-1.0380541.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0380541/manifest