Individual dosimetry system based on PHITS
Figure 1 shows the flowchart of our developed system. It can be divided into three processes: [1] conversion from PET-CT images to PHITS input files, [2] calculation of absorbed doses using PHITS, and [3] estimation of EQDX(α/β) based on the PHITS results coupled with the MK model. EQDX(α/β) as well as total dose and deposition energy in each voxel are converted in DICOM RT-DOSE format. Thus, they can be imported to commercial DICOM software for further analysis. Details of each process are described below.
Conversion from PET-CT images to PHITS input files
Firstly, the patient-specific voxel phantom in the PHITS input format is created from his/her CT image using the Image2PHITS module included in the DICOM2PHITS package [19]. Then, we adopted the correlation between CT numbers (Hounsfield Unit) and tissue parameters proposed by Schneider et al. [24] in this conversion, though users can define their own formula to represent the correlation in our system. The tallies for scoring the absorbed doses in Gy and deposition energies in MeV are also generated during this process. The resolutions of the created voxel phantom and mesh tallies are the same as the CT image.
A new module named PET2PHITS was developed in this study to create the maps of the cumulative activities as well as biological decay constants of the radionuclides based on the PET images. There are two types of patient-specific dosimetry systems; one is to create time-dependent activity maps and execute the particle transport simulations for each time step, and the other is to create a cumulative activity map and execute a single particle transport simulation. Using the former method, dynamical dose evaluation is possible by fitting the calculated doses for each time step. However, it is very time consuming because the Monte Carlo simulation needs to be continued until sufficiently small statistical uncertainties of the calculated doses in each voxel and time step are obtained to achieve the meaningful fitting. We therefore adopted the latter method; our system determines the cumulative activities and the biological decay constants of the radionuclides by fitting the dynamic PET images. Then, the dose rates are estimated under the assumption that they are proportional to the sum of the physical and biological decay constants of nearby voxels. The detail procedures for determining the cumulative activities and the biological decay constants are shown in Appendix A.
Calculation of absorbed doses using PHITS
Using the input files created from CT and PET images, PHITS simulation is performed to calculate the absorbed doses in the patient. In this study, PHITS version 3.20 was employed, and the EGS5 mode [25] was used for the photon, electron, and positron transport. The fluxes of the source particles including the contributions from daughter nuclides are determined from the RI source generation function in PHITS, based on ICRP Publication 107 [26]. The absorbed doses due to the ionisation induced by α and β± particles (referred to α and β doses, respectively) were separately calculated in the simulation. Note that the kerma approximation was not adopted, and thus, the photon doses were categorised as their secondary particle doses, i.e. β dose.
Before performing the particle transport simulation inside the patient body, another PHITS simulation must be performed to calculate the dose probability densities (PD) of lineal energy, d(y), in water for α and β doses, which are to be provided to the MK model for the RBE estimation. This simulation is required once for each radionuclide because it is not specific at each patient. The microdosimetric function of PHITS [27] is utilized for this calculation because the site size of y needed to be evaluated for the MK model is too small (less than 1 µm) to be handled with the condensed history method employed in EGS5. Note that the microdosimetric function was developed by fitting the results of track-structure simulation. Thus, it can analytically determine the PD of y down to the nanometer scales, considering the dispersion of deposition energies from the production of δ-rays. Figure 2 shows examples of the calculated PD of y for α and β doses of 211At.
Estimation of EQDX(α/β)
EQDX(α/β) is defined as the total absorbed dose delivered by the reference treatment plan (fraction size X) leading to the same biological effect as a test treatment plan [8]. Assuming that the biological effectiveness is proportional to the cell surviving fraction following a linear-quadratic (LQ) relationship, EQDX(α/β) for a test treatment with the surviving fraction S can be calculated by
where α and β are the LQ parameters for the reference treatment. Based on the MK model with the extensions of the saturation correction due to the overkill effect [28] and the dose rate effect [29], the cell surviving fraction in any radiation field with an absorbed dose D can be estimated by
where α0 is the linear coefficient of the surviving fraction with the limit of LET → 0, G is the correction factor due to the dose rate effect, and is the saturation-corrected dose-mean specific energy, deduced by
where y* is the saturation-corrected lineal energy, rd is the radius of a subcellular structure referred to as domain, y0 is a so-called saturation parameter that indicates the lineal energy above which the saturation correction due to the overkill effect becomes very important, and d(y) is the dose probability density in domain. d(y) in each voxel can be determined from its α and β doses, Dα and Dβ, respectively, as written by
where dα(y) and dβ(y) are their dose PD for each radionuclide precalculated by PHITS using the microdosimetric function. Assuming that the dose rates of TRT are expressed as a mono-exponential function with a decay constant of λphy + λbio, where λphy and λbio are the physical and biological decay constants, respectively, the value of G can be calculated using [13]
where μ is the recovery rate constant. The parameters α, β, μ, α0, rd, and y0 depend on the cell line. Among them, α0, rd, and y0 are specific to the MK model, and their determination requires the experimental data of cell surviving fractions for various ion irradiations, which are generally not available. Thus, we fixed rd and y0 to 0.282 μm and 93.4 keV/μm, respectively, which were evaluated from the surviving fractions of the HSG cell irradiated by various radiations including He ions [30, 31], and calculated α0 from α and for the reference radiation. Then, the user input parameters to our dosimetry system are α, β, and μ, which can be obtained from the measured surviving fractions of the reference radiation, as well as the fraction size X. Referring to our previous works [23, 30], we set α = 0.251 Gy-1, β = 0.0615 Gy-2, μ = 1.5 h-1 and X = 2 Gy in the test simulations performed in this study. Consequently, EQDX(α/β) calculated in this study can be expressed as EQD2(4.08).
EQDX(α/β) in a certain voxel can be simply calculated from Eq. 1 by substituting the surviving fraction in the voxel obtained from Eq. 2. In contrast, special care should be taken when EQDX(α/β) in a certain volume of interest (VOI) consisting of multiple voxels such as tumor and normal tissue is calculated because of the non-linear relationship between the EQDX(α/β) and the surviving fraction. In such cases, the mean surviving fraction in VOI, SVOI, is given by
where Si, Di and mi is the surviving fraction, dose, and mass, respectively, of voxel i made up of VOI. EQDX(α/β) in VOI can be obtained from Eq. 1 by supplying SVOI to S in similar to the concept of the equivalent uniform dose (EUD) [32]. DMH in VOI is the key quantity in this evaluation, which can be also calculated from our dosimetry system.
Dynamic PET-CT acquisition and analysis
This study was approved by the institutional review board, and written informed consents were obtained from all participants. The performance of the system was examined using the dynamic PET-CT data of four healthy volunteers after injecting 18F-labeled NKO-035 with 221.6 ± 3.8 MBq, which is a specific substrate of L-type amino acid transporter-1 (LAT1). The dynamic PET data were acquired in nine frames (total scan duration: 90 min) with low-dose CT scan. All images were depicted by OSIRIX (Newton Graphics, Inc. Sapporo, Japan). Details of the data acquisition procedures were described in Appendix B. In the dose estimation, NKO-035 was assumed to be labeled with not only 18F but also 211At, 131I, and 177Lu with the same distribution in the body. Volume of interest were placed in major organs on dynamic PET images using PMOD software (PMOD Technologies Ltd., Zurich, Switzerland) with reference to CT images. The residence times in major organs and tissues were estimated for each patient based on their dynamic PET data using the method described in Appendix A. Supplying those data into OLINDA 2.0, the organ doses were calculated and compared with the corresponding data obtained from our dosimetry system by Bland-Altman analysis.