Analysis with tacmagic

Eric Brown

2020-12-08

Background

Positron emission tomography (PET) is a research and clinical imaging modality that uses radioactive tracers that bind to target molecules of interest. A PET scanner identifies the tracer location by virtue of the tracer’s radioactive decay, providing information to determine the location of the target in the body. As the spatial resolution of PET is relatively poor, analysis is frequently combined with higher resolution imaging such as magnetic resonance imaging (MRI), which can be spatially co-registered to the PET image. Subsequently, radiotracer activity (over time) can be identified by spatial region of interest (ROI).

An image analysis pipeline is required to extract regional time activity curves (TACs) from a dynamic PET image. There are various pipelines available including widely-used commercial solutions (e.g. PMOD) and newer open-source options (e.g. magia1). Pipelines generally implement the following steps:

Various pipelines save TAC, volume and related data in various formats. This package enables the loading and analysis of TAC and ROI volume data from image analysis pipelines for further analysis in R.

Vignette data

The sample data for this vignette uses an anonymized scans of a participant with Alzheimer’s dementia, data from http://www.gaain.org which was generously made available for unrestricted use.2 The radiotracer used is Pittsburgh Compound B (PIB) which binds to beta-amyloid, a protein found in high concentration in the brains of individuals with Alzheimer’s dementia.

There are two approaches to using the tacmagic package to analyze PET time-activity curve data: either by loading participant data individually and using the various functions to analyze it, or via the batch functions to list and analyze data from multiple participants. Here, we illustrate the main features of tacmagic, by walking through the analysis of a single participant. We provide an explanation of the batch mode at the end.

Time-activity curve operations

Data loading

Time-activity curve (TAC) data is loaded via load_tac(), which is a wrapper for format-specific functions. To specify which file format the TAC data is stored as, use the format parameter. Supported formats can be viewed in help(load_tac).

The minimal amount of information required is the TAC data for one or more ROI, including the start and stop times of each frame, the time units and the activity units. This information may be in 1 or more files depending on the format and software that created it.

For example, PMOD’s .tac files contain all of the information, but the TAC .voistat files do not contain start and stop times, but this information could be specified using a .acqtimes file. Support is also available for DFT format, which contains both TAC and volume data.

We processed the PIB PET and T1 MRI data with the PMOD PNEURO software suite to produce a .tac file with TACs for all ROIs in the Hammer’s atlas. The .tac file can be loaded with load_tac():

A TAC object is a data frame with extra attributes including time and activity units. A summary can be printed with the generic print() function.

PMOD’s suite also produces .voistat and .acqtimes formats, than can be used to produce the same data if you do not have .tac files:

We also used Turku’s magia pipeline to process the same data. It can be loaded similarly, though with units explicitly entered because the information is not encoded in the .mat file:

Manually-created TAC objects

For other data sources, tacmagic TAC objects can be created from data.frame objects with as.tac(). The time and activity units must be specified as arguments if not already set as attributes in the data.frame. The columns of the data.frame are the regional TACs, with the column names the names of the ROIs.

Radioactivity unit conversion

Most modern PET tools use kBq/cc as the standard activity units required by the software (e.g. TPCCLIB, PMOD). Most often data will be in this format and not require conversion. However, conversion of TAC objects to these units, or to other radioactivity units for that matter, is possible with change_units(). This is a generic function that works on tac objects as well as numeric objects. For numeric objects, both “to” and “from” units need to be specified. The function works regardless of whether the units are per volume (i.e. kBq is treated the same was as kBq/cc or kBq/mL).

With tac objects, as the activity units are stored in the object, they should not be provided to change_units():

ROI merging

Often it is desirable to combine TAC ROIs into larger ROIs. For example, if the PET analysis pipeline created TACs for each atlas ROI, your analysis may call for merging these atomic ROIs into larger regions, such as merging all of the atlas ROIs that make up the frontal lobe into a single frontal lobe ROI.

If this is done, the means should be weighted for the relative volumes of the atomic ROIs. If volume information is available, tac_roi() provides this functionality.

In PMOD’s software, volume information is available in .voistat files. Units do not matter because it is the relative volume information that is needed.

In addition to TAC and volume information, we must specify which atomic ROIs make up the merged ROI. This is done by providing a named list, where the names are the merged ROIs and the list items are themselves lists of the atomic ROIs that make up each merged ROI. For the Hammer’s atlas, and as an example, typical data is provided in roi_ham_stand(), roi_ham_full(), or roi_ham_pib().

Plotting

Basic TAC plotting can be done by calling plot, which accepts two TAC objects, e.g. from 2 participants or group means. The ROIs to plot are specified as a vector of ROI names as they appear in the TAC object. As the TAC object contains time unit information, the plot can convert to desired units, which can be specified with the time argument.

Model calculation

Standardized uptake value (SUV)

As the activity in the TAC is impacted by the dose of the radiotracer administered and the participant’s body weight, a value adjusted for these factors is sometimes used, the SUV:

\[SUV = \frac{Ct}{\frac{Dose}{Weight}}\]

Where activity is measured in kBq/mL (kBq/cc), the dose is in MBq, and the weight is in kg, the radioactivity units cancel and the SUV units are g/mL.

With the tac_suv() function, a tac object can be converted to SUV values with units g/mL, as demonstrated below (note the weight and dose are fabricated here for the demonstration).

More often, an everage value over a certain time period, or a maximum value may be desired. This can be calculated with the suv() function.

SUV ratio (SUVR)

The standardized uptake value ratio (\(SUVR\)) is a simple quantification of PET activity that is commonly used from many tracers including PIB. It is the ratio of the tracer activity over a specified time period (\(Ct\)) in a target ROI to a reference region. Using a ratio allows factors that are normally required to calculate an \(SUV\) to cancel out, namely tracer dose and patient body weight, and therefore \(SUVR\) can be calculated from TAC data alone, i.e. without the need to specify tracer dose or body weight:

\[SUVR = \frac{SUV_{TARGET}}{SUV_{REF}} = \frac{Ct_{TARGET}}{Ct_{REF}}\]

In the literature, SUVR is variably described and calculated using the mean of activity for the frames of the specified time period, or the area under the curve. For PIB, the mean/summed activity has been used, and the time windows have varied from starting at 40-50 minutes and ending at 60-90 minutes.3

The suvr() function calculates SUVR for all regions in a TAC file based on the provided time information (as a vector of frame start times) and the specified reference region (a string). If the frames used are of different durations, the weighted mean is used.

An alternative method, using the area under the curve with the mid-frame times as the x-axis is available with suvr_auc() and should provide very similar results.

DVR

The Distribution Volume Ratio (DVR) is a method of quantifying tracer uptake that is used as an alternative to the SUVR in PIB studies, for example. Like SUVR, it can be calculated from TAC data without the need for arterial blood sampling, by making use of a reference region. In this case, it is called the non-invasive Logan plot method. It is calculated with a graphical analysis technique described by Logan et al.4

In addition to the TAC data, depending on the tracer, a value for k2’ may need to be specified. For PIB, this has limited effect on the result, but can be specified, and a value of 0.2 has been recommended.5

The non-invasive Logan plot works by finding the slope of the line of the following equation after time \(t*\) where linearity has been reached:

\[\frac{\int_0^{T}C_{roi}(t)dt}{C_{roi}(t)} = DVR[\frac{\int_0^{T}C_{cer}(t)dt + C_{cer}(t) / k2`}{C_{roi}(T)}] + int \]

Find t*

The time, \(t*\) (t_star), after which the relationship is linear can be found by testing the point after which the error is below a certain threshold (default is 10%). If t_star=0, then tacmagic can find the suitable value.

To visually confirm that the model behaved as expected with linearity, there is a plotting function:

The right plot shows the Logan model, with the vertical line representing the identified \(t*\), and the linear model fitted to the points after that time. In this case, the line after \(t*\) can be seen to fit well. The slope of that line is the DVR.

Similarly, DVR can be calculated for all ROIs, either by setting t_star manually or to 0 as before. If 0, a different value will be identified for each ROI.

For this data, the DVR calculation has been shown to produce equivalent results as an existing tool.6

A wrapper function dvr() is available to conveniently calculate DVR for a target ROI or all ROIs, and currently defaults to using the Logan reference method:

Batch analysis

In most cases, a project will involve the analysis of multiple participants. The above workflow can be used to test and visualize an analysis, but a batch workflow will likely be preferred to analyze multiple participants.

All analyses can be run using 2 steps: a batch data loading step and a batch analysis step.

Batch loading

Data loading is done by batch_load(). See help(batch_load) for the required arguments.

The first argument is a vector of participant IDs that corresponds to file names, e.g.:

participants <- c("participant01", "participant02") if the files are located e.g. /mypath/participant01.tac and /mypath/participant01_TAC.voistat. In this case, the function call might look like:

my_data <- batch_load(participants, dir="/mypath/", tac_format="PMOD", roi_m=T, vol_file_suffix="_TAC.voistat", vol_format="voistat", ROI_def=roi_ham_stand(), merge=F)

The above would load the appropriate TAC and voistat files, perform the ROI merging specified by ROI_def, because roi_m = TRUE, and would return a list where each element represents a participants, e.g. the first participant would be my_data$participant1.

To calculate SUV in batch, the participants dose and weight must be specified when loading with batch_load(), as it is then added to the respective tac objects.

Batch analysis

Once the TAC data is loaded, all analyses can be run using batch_tm(). The output from batch_load() is the first argument for batch_tm(). The models implemented in tacmagic can be specified using the models argument, e.g. models = c("SUVR", "Logan") to calculate both SUVR and Logan DVR. The relevant model parameters will also need to be specified, so see help(batch_tm) for all possible arguments.

Batch example

For the purpose of the vignette, the list of participants will be a list of the full TAC filenames (hence tac_file_suffix=""). In real-world data, the participants parameter can be a list of participant IDs that correspond to the actual filenames, i.e. the filename is made up of dir + participant + tac_file_suffix.

We will also choose not to use the roi_m option in batch_load(), which could be used to combine ROIs as outlined above.

Cut-off calculations

In the analysis of PIB/amyloid PET data, often researchers want to dichotomize patients into PIB+ vs. PIB-, i.e. to identify those with significant AD-related amyloid pathology (PIB+).

There are a number of approaches to this depending on the available data. We have implemented a method described by Aizenstein et al.7 which uses a group of participants with normal cognition to establish a cutoff value above which participants are unlikely to have minimal amyloid pathology.

The method identifies a group of participants out of the normal cognition group with higher-PIB outliers removed. An outlier is a participant with any ROI with a DVR higher than the upper inner fence, from a set of ROIs known to be associated with amyloid deposition. Such participants are removed from the group, and this process is done iteratively until no more outliers are removed. Then, cutoff values are determined from this new group for each ROI, set again as the upper inner fence. Then these cutoff values are applied to all participants, and a participant is deemed PIB+ if they have at least 1 ROI above its cutoff.

To demonstrate, a fake dataset of DVR values for 50 fake participants was generated and is available as fake_DVR. This would be equivalent to using batch_tm() on a group of participants with the "Logan" model specified.

fake_DVR[1:5,]
#>    ROI1_DVR ROI2_DVR  ROI3_DVR  ROI4_DVR
#> p1 2.235058 3.026140 1.1420583 4.0557974
#> p2 1.677543 2.834710 1.2529424 3.3368371
#> p3 1.466090 2.856528 0.2837647 2.5383574
#> p4 1.046542 1.818056 1.2176133 2.5827267
#> p5 2.934105 1.971775 1.7611040 0.8371429

To calculate the cutoff values using this iterative method, cutoff_aiz() takes 2 arguments: the DVR data, and the names of the variables of the ROI DVRs to use (and there must be at least 2 for this method).

cutoffs <- cutoff_aiz(fake_DVR, c("ROI1_DVR", "ROI2_DVR", "ROI3_DVR", "ROI4_DVR"))
#> Iteration: 1 Removed: 10
#> Iteration: 2 Removed: 1
#> Iteration: 3 Removed: 0

cutoffs
#> ROI1_DVR ROI2_DVR ROI3_DVR ROI4_DVR 
#> 1.370899 1.332725 1.330504 1.273470

The final step is to apply the cutoffs to the full set of participants. We will use the same sample data:

positivity_table <- pos_anyroi(fake_DVR, cutoffs)

positivity_table
#>    p1    p2    p3    p4    p5    p6    p7    p8    p9   p10   p11   p12   p13 
#>  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE 
#>   p14   p15   p16   p17   p18   p19   p20   p21   p22   p23   p24   p25   p26 
#> FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
#>   p27   p28   p29   p30   p31   p32   p33   p34   p35   p36   p37   p38   p39 
#> FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE 
#>   p40   p41   p42   p43   p44   p45   p46   p47   p48   p49   p50 
#> FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

The algorithm identified 11 PIB+ participants. In the generation of the sample data, the DVRs from the first 10 participants were drawn from a normal distribution with mean 1.9, sd 0.6 and for the latter 40 participants, from mean 1.3, sd 0.3; thus this pattern is in line with what we would expect: all 10 of the first participants are PIB+, and just 1 of the latter 40 was (by chance).

  1. Tomi Karjalainen, Severi Santavirta, Tatu Kantonen, Jouni Tuisku, Lauri Tuominen, Jussi Hirvonen, Jarmo Hietala, Juha Rinne, Lauri Nummenmaa. Magia: Robust automated modeling and image processing toolbox for PET neuroinformatics. bioRxiv 604835; https://doi.org/10.1101/604835↩

  2. Klunk, William E., Robert A. Koeppe, Julie C. Price, Tammie L. Benzinger, Michael D. Devous, William J. Jagust, Keith A. Johnson, et al. “The Centiloid Project: Standardizing Quantitative Amyloid Plaque Estimation by PET.” Alzheimer’s & Dementia 11, no. 1 (January 2015): 1-15.e4. https://doi.org/10.1016/j.jalz.2014.07.003.↩

  3. Lopresti, B. J., W. E. Klunk, C. A. Mathis, J. A. Hoge, S. K. Ziolko, X. Lu, C. C. Meltzer, et al. “Simplified Quantification of Pittsburgh Compound B Amyloid Imaging PET Studies: A Comparative Analysis.” J Nucl Med 46 (2005): 1959–72.↩

  4. Logan, J., Fowler, J. S., Volkow, N. D., Wang, G.-J., Ding, Y.-S., & Alexoff, D. L. (1996). Distribution Volume Ratios without Blood Sampling from Graphical Analysis of PET Data. Journal of Cerebral Blood Flow & Metabolism, 16(5), 834-840. https://doi.org/10.1097/00004647-199609000-00008↩

  5. http://www.turkupetcentre.net/petanalysis/analysis_11c-pib.html↩

  6. https://gitlab.utu.fi/vesoik/tpcclib↩

  7. Aizenstein HJ, Nebes RD, Saxton JA, et al. 2008. Frequent amyloid deposition without significant cognitive impairment among the elderly. Arch Neurol 65: 1509-1517.↩