Proton range uncertainty caused by synthetic computed tomography generated with deep learning from pelvic magnetic resonance imaging

Abstract

Background

In proton therapy, it is disputed whether synthetic computed tomography (sCT), derived from magnetic resonance imaging (MRI), permits accurate dose calculations. On the one hand, an MRI-only workflow could eliminate errors caused by, e.g., MRI-CT registration. On the other hand, the extra error would be induced due to an sCT generation model. This work investigated the systematic and random model error induced by sCT generation of a widely discussed deep learning model, pix2pix.

Material and Methods

An open-source image dataset of 19 patients with cancer in the pelvis was employed and split into 10, 5, and 4 for training, testing, and validation of the model, respectively. Proton pencil beams (200 MeV) were simulated on the real CT and generated sCT using the tool for particle simulation (TOPAS). Monte Carlo (MC) dropout was used for error estimation (50 random sCT samples). Systematic and random model errors were investigated for sCT generation and dose calculation on sCT.

Results

For sCT generation, random model error near the edge of the body (∼200 HU) was higher than that within the body (∼100 HU near the bone edge and <10 HU in soft tissue). The mean absolute error (MAE) was 49 ± 5, 191 ± 23, and 503 ± 70 HU for the whole body, bone, and air in the patient, respectively. Random model errors of the proton range were small (<0.2 mm) for all spots and evenly distributed throughout the proton fields. Systematic errors of the proton range were −1.0(±2.2) mm and 0.4(±0.9)%, respectively, and were unevenly distributed within the proton fields. For 4.5% of the spots, large errors (>5 mm) were found, which may relate to MRI-CT mismatch due to, e.g., registration, MRI distortion anatomical changes, etc.

Conclusion

The sCT model was shown to be robust, i.e., had a low random model error. However, further investigation to reduce and even predict and manage systematic error is still needed for future MRI-only proton therapy.

Keywords:

• MRI-only
• proton range uncertainty
• synthetic CT
• deep learning
• pelvis

Introduction

Conventional proton therapy uses computed tomography (CT) as the standard imaging modality for dose calculation [1]. Recently, it has been proposed to use magnetic resonance imaging (MRI) as an alternative to CT with the potential to improve radiotherapy by exploiting the advantages of MRI, e.g., the unparalleled high soft-tissue contrast and the absence of ionising radiation exposure [2]. However, MRI does not provide information on the proton-stopping power required for dose calculation [1]. A strategy to obtain this information is to, first, generate a so-called synthetic CT (sCT) from the MRI and, second, to transfer the sCT into a stopping power ratio (SPR) map.

Besides conventional analytical tools for sCT generation, i.e., bulk density [3–9] and atlas-based [5,6,10–14] methods, deep learning (DL) based methods have shown great potential in terms of, e.g., better adaptation to tissue heterogeneity and lower computational resource requirements [1]. Dose calculation on sCT has been widely investigated for photon radiation, while proton dose has been much less discussed. According to a recent review in 2021 [2], only a comparably small number (9 out of 57) of state-of-the-art studies investigated proton dose on generated sCT [15–23].

Proton and other ion therapies benefit from and are characterised by their dosimetric advantages: a low entrance dose and a distinct dose maximum, called Bragg peak, near the end of their finite beam range. Besides, the linear energy transfer (LET) of a proton beam, which relates to the biological effectiveness [24], also steeply increases near the Bragg peak at the end of the range. The proton range in patients is uncertain due to, e.g., CT number conversion, and intra- and inter-fractional anatomical changes. A safety margin around the target volume has to be applied in proton therapy to consider the uncertainties of the position of Bragg peaks that are placed at the edge of the target volume. The size of the safety margin depends on proton range uncertainty, i.e., large proton range uncertainty leads to the extension of the margin, and thereby, substantially increases the high dose volume in normal tissue. Calculating proton dose on sCT, which is not directly measured but derived from MRI using models, would induce extra proton range uncertainties, as the proton range is sensitive to CT number errors. It is important to know and ideally reduce the induced uncertainty caused by sCT generation. Until recently, only a few studies investigated proton range on generated sCT [2] and these were mostly related to brain and head and neck patients [15,16,18,23].

This study focuses on quantifying model-induced proton range errors resulting from sCT generation. For sCT generation, a prototype model of a widely investigated [25–36] conditional generative adversarial network (cGAN), pix2pix [37] was applied. The open-access dataset of Gold Atlas [38] was employed.

Material and methods

Patient data

This study used MRI-CT pairs of nineteen patients with cancer in the pelvis acquired by three different Swedish departments. They were used as provided and described by the open-access dataset Gold Atlas, which was initiated for the purpose of establishing segmentation and sCT generation methods [38]. In brief, the number of patients contributed by department 1, 2 and 3 are eight, seven, and four, respectively. The MRI images of departments 1, 2 and 3 were obtained with different machines, namely, 3 T GE Discovery, 1.5 T Siemens, and 3 T GE Signa, respectively [38]. Correspondingly, CT scanners were: Siemens Somatom Definition AS+, Toshiba Aquilion, and Siemens Emotion 6, respectively [38]. Furthermore, the authors deformably registered the CT to the T2-weighted MRI images resulting in a set of two-dimensional (2D) MRI-CT slice pairs for each patient. The deformed and resampled CT and the T2-weighted MRI images both had a voxel size of 0.875 × 0.875 × 2.5 mm³ [38].

In this work, it was assumed that the CT in the dataset was the planning CT on which the treatment plan dose was supposed to be calculated. Thus, the CT and the dose calculated on the CT were regarded as ground truth. For modelling, the nineteen patients were randomly split into training (ten patients, 880 pairs of 2D MRI-CT slices in total), test (five patients, 412 pairs of 2D MRI-CT slices), and validation (four patients, 237 pairs of 2D MRI-CT slices).

Preprocessing

First, all images were converted to have a common image size of 500 × 376 pixels (pixel size unchanged) by adding/removing air-filled areas at the edges of all paired CT and MRI images. The image intensities of MRI images from the three sites were linearly transferred so that the differences between histograms of MRI intensities of the three MR scanners were minimised. Based on the experience of [25], where the same dataset and a similar DL model were used, in this work, the following preprocessing steps were applied to the MRI data: 1) N4ITK bias field correction [39] with parameters of convergence threshold 0.004, number of histogram bins 200, Wiener filter noise 0.01, downsampling 3. 2) A curvature anisotropic diffusion filter was applied with the settings of conductance 1, scaling update interval 1, time step 0.03, iterations 10. The MRI intensities were truncated from 0 to 2165 (98% maximum), i.e., voxels with initial intensities higher than 2165 were set to 2165. For CT images, the intensities were truncated from −600 HU to 1400 HU, as CT numbers below −600 were considered air (and set to −600), and CT numbers above 1400 were identified (and set to 1400) only in a limited number of cases [25]. The percentage of voxels with CT numbers higher than 1400 HU and between −600 HU and −1000 HU in the body of the initial CT before truncation were both less than 0.1%. Afterwards, both MR and CT images were normalised to the range of [-1] for DL modelling.

Model architecture and hyperparameter settings

For the sCT prediction, a cGAN model extension, pix2pix [37], was employed. The original authors observed that this model allows training using a comparably small number of image pairs, e.g., 400 [37]. Different models based on pix2pix have been discussed in some detail elsewhere [25–35] and it was not the intention of this work to make a comparison between those approaches of applying pix2pix. Thus, the prototype architecture of pix2pix was applied. The input and output of the model were MRI and its corresponding CT, both were 3-channel 2D images. The three channels included a central 2D image slice with the two adjacent slices above and below it. Previous work on hyperparameter search showed that different settings of learning rate, number of epochs, and the trade-off λ had no significant improvements over default settings [25]. Thus, confirmed with patients in the validation set, the corresponding default hyperparameters, 2 × 10⁻⁴, 200, 100, respectively, were used.

For training, paired image sections with a size of 256 × 256 pixels (i.e., the default input size for pix2pix) were randomly cropped from the MRI-CT image pairs and then entered into the model. The pix2pix default dropout rate of 50% was applied, i.e., 50% of the nodes were randomly and temporarily excluded during training to avoid overfitting [40]. For model testing, four partially overlapping 256 × 256-pixel image parts, each aligned with a corner of the MRI, were cut from the input MRI and then fed into the model. The resulting four model outputs were combined to produce a full sCT prediction of this MRI. For the pixels where partial images overlapped, averages were taken from the image intensities.

Particle transport simulation for proton range

The tool for particle simulation (TOPAS) [41] was used to simulate the proton irradiation on the real CT and sCT for proton range comparison. A grid of 11 × 11 proton pencil beams (called spots in the following) was simulated. Proton spots were evenly distributed within a square area of 10 × 10 cm² and entered the patient under 90° from the right side. The energy for all spots was set to 200 MeV (proton range in water of about 26 cm), which is usually among the highest energies for clinical treatments. The relative energy spread was set to 0.7 and the lateral spot width was set to 0.65 cm. For each spot, 1 × 10⁵ protons were simulated and a 3D dose distribution was recorded using the dose to water scorer and a dose grid corresponding to the CT grid. The proton range was calculated by the distance between spot entrance into the patient and the 80% distal dose fall-off of the laterally integrated depth dose profile.

Evaluation

In this work, model-induced errors [42–46], i.e., systematic and random model errors, were investigated. The systematic error was estimated by the difference between the estimation of the DL prediction and the reference while the random model error was measured by the standard deviation (STD) of the prediction samples generated from the same input. For this purpose, a Monte Carlo (MC) dropout strategy [47] was employed to generate sCT prediction samples: by applying MC dropout on a trained model, a sub-model with a predefined rate of randomly excluded nodes can be obtained, allowing for an individual complete sCT prediction of an input MRI. MC dropout can be repeated n times to obtain n individual sub-models and sCT predictions for the same MRI. In the following, these sCT predicted by individual sub-models for the same MRI are called sub-sCT. In this work, the dropout rate was set to 50% and for each input MRI, 50 sub-sCT were predicted.

For the sCT generation, the mean of the CT numbers of all sub-sCT in a voxel represents the estimation of the sCT voxel prediction for a given MRI. Likewise, the random model error of an sCT voxel prediction is measured by the STD of the CT numbers of the sub-sCT in this voxel. The systematic error was measured by the mean absolute error (MAE) and mean error (ME) of the CT numbers between the real CT and the estimation of sCT.

For the proton range, the errors associated with the random model error and the systematic error of the prediction were determined. For a proton spot, the proton ranges were simulated on all 50 sub-sCTs and the obtained STD in the range was considered as the random model error of the range for that spot. The systematic error of the range prediction for a proton spot was measured by the difference between its range on the real CT and on the estimation of the sCT, i.e., the mean of the sub-sCT. The Spearman correlation between the systematic proton range error of a spot and the mean along its beam path of 1) the sCT systematic error, 2) the square of the sCT systematic error and 3) the sCT random model error was calculated. Random model error of the proton range was not considered in the correlation analysis as this error was observed to be very small for all proton spots (c.f. Section ‘Model-induced proton range error’').

Results

Model-induced sCT errors

For sCT prediction (example cf. Figure 1), both systematic and random model errors were not uniformly distributed over a slice and were found to be higher at the edges of different (anatomical) objects, especially at the body contour. The random model error at the edge of the patient was much higher (∼200 HU) than in the other considered areas inside the body (mostly less than 10 HU) and also higher than at the bone edges (∼100 HU).

Figure 1. Example of MRI (a), CT (B), and sCT (C) for a corresponding slice, along with associated random model error (i.e., the standard deviation of the 50 sub-sCTs; D) and the systematic error of sCT (i.e., the difference between the CT and the sCT; E).

The histograms of the systematic and the random model error (displayed for error values larger than 10 HU) for different areas are shown in Figure 2. The four areas body, bone, air, and edge correspond to the entire region within the body contour, regions with CT numbers of the real CT greater than 250 HU, regions within the body contour with CT numbers smaller than −200 HU, and a band at the body contour with a width of 7 voxels, respectively. The systematic error of body, bone, air, and soft tissue, i.e., all voxels within the body contour except those from the air and bone areas, are listed in Table 1. The systematic error within the body contour in areas with air on the real CT was systematically biased, i.e., the model had difficulties predicting air cavities that are found on the CT based on the given MRI.

Figure 2. Histogram of the (A) random model and (B) systematic error for sCT in the test dataset Separated according to area defined by the CT numbers of the real CT. Relative histograms are shown as the volumes of the considered areas differ substantially. For random model error, data with values < 10 HU are omitted to highlight the behaviour at large errors.

Table 1. The mean absolute error (MAE) and mean error (ME) between the sCT and CT Evaluated separately for the training and test dataset in different areas.

	MAE (HU)				ME (HU)
Dataset	Body	Bone	Air	Soft tissue	Body	Bone	Air	Soft tissue
Training	37 ± 3	138 ± 15	421 ± 88	26 ± 2	1 ± 7	28 ± 27	−409 ± 87	1 ± 5
Test	49 ± 5	191 ± 23	503 ± 70	33 ± 5	12 ± 11	102 ± 41	−494 ± 76	3 ± 10

Model-induced proton range errors

The random model error of the predicted proton range was low (< 0.2 mm) for all spots. It was distributed evenly across the 10 × 10 cm² spot grid (Figure 3A) and comparable for all patients.

Figure 3. Beam’s eye view (BEV) of the (A) random model error and (B) systematic error of proton range of 121 individually simulated lateral 200 MeV proton beams on a 10 × 10 cm² beam grid oriented in anterior-posterior (AP) and cranial-caudal direction (CC) for one patient in the test dataset.

On the other hand, the systematic proton range errors were larger and were not evenly distributed in space (Figure 3B) and varied between patients. Absolute and relative values of systematic proton range errors (± STD) were estimated to be −1.0 (± 2.2) mm and −0.4 (± 0.9)%, respectively, for the test dataset and 0.1 (± 1.6) mm and 0.1 (± 0.6)% for the training dataset, respectively. The median absolute proton range error for the test and training datasets were 1.5 and 0.7 mm, respectively. The histograms of absolute and relative values of systematic proton range errors for all spots are shown in Figure 4A and B, respectively.

Figure 4. Histogram of (A) absolute and (B) relative values of systematic proton range errors of all beams simulated on the training and test dataset. The number of proton beams is normalised.

Figure 5 shows laterally integrated proton depth dose profiles as well as CT and sCT numbers along the beam path for four typical situations that resemble one ideal sCT prediction (I) as well as three predictions (II–IV) resulting in large systematic range errors for the considered spots. The corresponding CT, MRI, and sCT slices for these four proton spots are presented in Figure 6. It can be observed that:

Figure 5. Comparison of the distal part of the laterally integrated dose profiles (A) and corresponding CT number profiles along the beam path (B) of four exemplary proton beams (I-IV) simulated on corresponding CT (red) and sCT (blue) slices.

Figure 6. Comparison of CT (a), MRI (B) and corresponding sCT (C) for the same four exemplary patient slices (I-IV) as in Figure 5 with the simulated beam paths indicated as red lines.

For case I, the depth dose distribution remained unaltered as the small differences between CT and sCT numbers averaged out along the beam path. For case II, air cavities in the CT substantially degraded the Bragg peak but were not predicted on the sCT resulting in a shorter range and steeper distal dose fall-off. For case III, the overall agreement between CT and sCT numbers was good, except for the maximum CT numbers of the sCT in the high-density bone areas, which were predicted to be lower than those in the CT resulting in a less degraded shape of the Bragg peak curve. For case IV, the shape of the sCT was consistent with the MRI but differed, particularly, in the outer regions of the image from the CT. There, the sCT appeared to be stretched, resulting in an effectively shorter proton range in the patient compared to the CT.

Correlation of errors

Spearman correlation between the relative values of systematic proton range error and the mean of the systematic sCT error along the beam path was observed to be −0.78. Low correlations were observed between the systematic proton range errors and the mean of the squared systematic sCT errors as well as the mean of the random sCT errors in the beam path, which are −0.04 and −0.03, respectively.

Discussion

In this work, a cGAN DL model, which showed promising performance for sCT prediction in previous studies [25–35], was trained and tested using the open-source dataset Gold Atlas [38], which contains CT-MRI pairs for 19 patients with cancer in the pelvis. Afterwards, the errors induced by the model were analysed. For the systematic error of sCT generation, the model performance was consistent with that in the previous work by Fetty et al. [25] using the same DL model and dataset. The random model error of sCT prediction was high at the edge of a patient, moderate at the bone edges, and much lower inside the body.

The systematic range errors of 200 MeV proton spots resulting only from the sCT prediction were distributed for the test dataset on the absolute and relative scales as −1.0 (± 2.2) mm and −0.4% (± 0.9)%, respectively. The random model error of the proton range, on the other hand, was found to be much smaller and below 0.2 mm for all beams. The median absolute proton range error was 1.5 mm, while other studies reported a median absolute error in proton range estimation with sCT of 1.7 mm for the liver [22] and 0.5 mm for the brain [48]. Although different models for different anatomical regions and datasets are not directly comparable [1], the errors observed in this work are in the same order of magnitude as those reported in the literature.

For 4.5% and 1.1% of the spots in the test and training dataset, respectively, large systematic proton range errors (more than 5 mm) were observed. These large errors were usually related to, e.g., objects such as low-density air cavities or high-density bone on CT that were not predicted 'correctly’ (cf., Figures 5 and 6, cases II and III), or general mismatch between CT and MRI (case IV). This is consistent with the systematic error maps between CT and sCT, where large errors in CT numbers were found near the boundaries of objects, e.g., between bone, tissue, and air. In accordance with this, random model error maps showed high uncertainty values at the edge of the body contour, i.e., the robustness of the predicted CT numbers is low at the edge of the body. We believe that this is also related to the mismatch in spatial extent, especially near the outer contour, between the CT and MRI in the training dataset. Previous work reported [22] that the sCT predicted by their cycle GAN model matched the associated MRI better than the corresponding real CT, which is consistent with our findings. Future applications of MRI-only proton therapy may benefit if the prediction algorithm generates an sCT that matches the input MRI (e.g., from daily imaging), rather than forcing an sCT geometry that matches that of the original CT, because organ locations and anatomy may differ in between MRI and CT imaging [22]. However, on the one hand, unwanted sCT errors could occur due to possible geometrical distortions or artefacts during MRI reconstruction. On the other hand, anatomical mismatch affects both model training (due to spoiled optimisation of the loss function) and testing (due to imperfect 'ground truth’). This mismatch issue might primarily be mitigated by improving data quality. First, image data pre-processing including accurate spatial distortion correction algorithms for MRI image reconstruction and CT-MRI registration could be improved. Research [49] reported that improved registration had a positive impact on sCT generation model training. Second, the similarity between the underlying MRI and CT images could be increased by minimising the effects of, e.g., setup changes (patient positioning) and anatomical changes (e.g., different gas filling in organs) between image acquisitions. Furthermore, improving the ground truth by using specific MRI sequences that provide better bone contrast is an option. A recent work using a Zero Echo Time (ZTE) MRI sequence, which provides better information on cortical bone and faster scanning times (65 s), achieved MAE (±STD) values within external contour, soft-tissue, and bone of 36 (±3), 7 (±1) and 98 (±14) HU for patients with cancer in the pelvis [50].

Error prediction indicating the confidence of model prediction would help in the future application of DL-based MRI-only proton therapy. A recent study of [42] reported that the square of the prediction error, including systematic error and random error, can be predicted using DL methods. The Spearman correlation between the systematic proton range error and the systematic sCT error was high. However, the correlation between the systematic proton range error and the currently predictable square of the prediction error was low. Future studies concerning the prediction of proton range error are needed.

This work has several limitations. The underlying CT-MRI dataset [38] was taken directly from the literature. MRI sequence selection and preprocessing, e.g., MRI reconstruction and deformable registration remained unchanged. Specifically, the registration strategy may not guarantee an ideal match of all objects on CT and MRI. The CT images were deformably registered using parameter files that were experimentally optimised through visual inspection. Image similarity was measured by mutual information and B-splines were used as a transformation model. The bladder contour delineated by five human experts was used for guidance [38]. Higher sCT errors were found at the edge of the body where a clear mismatch can be observed (cf. Figure 1, Figure 6 case I and IV). To our knowledge, there is no solution available for 'ideal’ pixel-to-pixel registration that allows 'ideal’ match for all objects. Further investigation of registration strategies, e.g., registration based on air-soft tissue-bone instead of clinical delineation, is desired for future development of sCT generation in proton therapy. Besides, more details on, e.g., CT scan voltage and time gap between the acquisition of MRI and CT are missing. These factors might limit the performance of the model. However, as previously suggested by [1], we also believe that keeping the dataset unchanged is useful for direct comparison between our work and potential future work. We avoided applying any extension of the pix2pix model that has been proposed by others [25–36] as most of these extensions were based on a specific dataset and it was not the intention to compare or validate them in this work. We believe that the most important cause of error was the mismatch between MRI and CT, which can hardly be resolved by changes to the model. In addition, the clinical impact of uncertainties, e.g., in patient setup, image acquisition, and CT-SPR conversion, was not quantitatively considered. Besides, the energy of all proton spots was set to 200 MeV, i.e., the spots did not resemble a clinical treatment field and, accordingly, no clinical dose distributions such as dose volume histograms were considered. Our motivation was to focus generally on range errors in sCT of proton spots that could be placed at the distal edge of the target volume in pelvic cases. Instead of sCT, in the future, SPR maps could be directly predicted using MRI and dual energy CT [51] to further reduce range uncertainty in proton beam simulations [52] and, eventually, the necessary dose margins around the target volume.

Conclusion

A cGAN model using an open-source pelvis MRI-CT dataset (Gold Atlas) was trained for uncertainty analysis. Systematic and random model errors of both the generated sCT and proton ranges on the sCT were estimated. The systematic errors of both sCT generation and proton range were comparable to similar work reported by others. The random model error of proton ranges was in general small and lower than 0.2 mm. The random model error of the sCT near the body contour (∼200 HU) was higher than that near a bone-soft tissue border (∼100 HU) and particularly higher than elsewhere in the body (<10 HU). Eliminating the mismatch between MR-CT pairs, especially at the edge of the body and for air-filled cavities, will be critical to improve proton range accuracy and enable future MRI-only proton therapy.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from [38]. Restrictions apply to the availability of these data. Data are available from https://zenodo.org/record/583096 with the permission of the authors.

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.