Empirical quenching correction in radiochromic silicone-based three-dimensional dosimetry of spot-scanning proton therapy.

BACKGROUND AND
PURPOSE: Three-dimensional dosimetry of proton therapy (PT) with chemical dosimeters is challenged by signal quenching, which is a lower dose-response in regions with high ionization density due to high linear-energy-transfer (LET) and dose-rate. This study aimed to assess the viability of an empirical correction model for 3D radiochromic silicone-based dosimeters irradiated with spot-scanning PT, by parametrizing its LET and dose-rate dependency.
MATERIALS AND METHODS: Ten cylindrical radiochromic dosimeters (Ø50 and Ø75 mm) were produced in-house, and irradiated with different spot-scanning proton beam configurations and machine-set dose rates ranging from 56 to 145 Gy/min. Beams with incident energies of 75, 95 and 120 MeV, a spread-out Bragg peak and a plan optimized to an irregular target volume were included. Five of the dosimeters, irradiated with 120 MeV beams, were used to estimate the quenching correction factors. Monte Carlo simulations were used to obtain dose and dose-averaged-LET (LETd) maps. Additionally, a local dose-rate map was estimated, using the simulated dose maps and the machine-set dose-rate information retrieved from the irradiation log-files. Finally, the correction factor was estimated as a function of LETd and local dose-rate and tested on the different fields.
RESULTS: Gamma-pass-rates of the corrected measurements were >94% using a 3%-3 mm gamma analysis and >88% using 2%-2 mm, with a dose deviation of <5.6 ± 1.8%. Larger dosimeters showed a 20% systematic increase in dose-response, but the same quenching in signal when compared to the smaller dosimeters.
CONCLUSION: The quenching correction model was valid for different dosimeter sizes to obtain relative dosimetric maps of complex dose distributions in PT.

Introduction

Three-dimensional (3D) dosimetry can be a valuable tool for experimental validation of radiotherapy (RT) delivery due to its capability of providing dose measurements of high spatial resolution [1], [2], [3], [4]. State of the art proton therapy (PT) is based on spot-scanning, where a thin pencil beam of protons covers the target by scanning it in varying energy layers [5]. However, the clinical use of 3D dosimetry in PT is challenged by steep dose gradients and regions of high linear-energy-transfer (LET) [6], [7], [8]. Dose-rate dependency has also been observed for a variety of systems, but there is so far limited data for the dose-rate ranges pertinent to spot-scanning PT [1], [9], [10], [11].

In this study, we investigated a radiochromic silicone-based dosimeter containing leucomalachite green (LMG) as dye material and chloroform as an initiator. The chemical process that gives dose-response involves three main reactions: activation of free radicals via irradiation; activation of LMG by free radicals; and conversion of the colourless LMG to malachite green (MG), which is light-absorbing. The read-out of the integrated dose can be done by measuring the local absorption coefficient with an optical CT-scanner [12]. The dosimeter formulation chosen was optimized in a previous study from our group [9], and has been characterized in terms of dose and energy dependence [8], [9], [13], temporal stability [14] and physicochemical and mechanical properties [7], [15], [16].

Similar to most solid-state and chemical dosimeters, radiochromic dosimeters present a linear-energy transfer (LET) dependent response, with significant quenching at the distal end of the Bragg peak [17], [18], [19]. Quenching is the non-linear dosimeter response for low proton energies that can be understood on the microscopic level by track-structure theory [6], [20], [21]. According to this theory, the spatial dose distribution around the primary proton track is dictated by the range of secondary electrons (-rays). The -rays generated by high energy protons (lower LET) can travel far and thus have a high possibility of initiating a chemical process that leads to a measurable signal. The -rays from low energy protons (higher LET) deposit their energy very close to the primary proton track, and their energy is more likely to be lost in processes that do not generate signal.

Additionally, dose-rate dependency has been reported for radiochromic dosimeters for proton beams [9] but has so far not been investigated in the clinical range of dose-rates used for spot-scanning proton therapy. The decrease in signal for higher dose rates is likely related to an increase in competing radiation-induced chemical reactions in the dosimeter [1] which do not convert LMG to MG.

Track- as well as dose-averaged LET (LETd) have been successfully used in the calibration of a variety of dosimetric systems, such as polyacrylamide gels [21] and plastic scintillators [19]. Its use in radiochromic silicone-based dosimeters is supported by previous studies [8], [22], which found a linear response for proton beams up to 50 Gy for each measured LETd interval.

The aim of this study was therefore to investigate the viability of a calibration procedure for this 3D dosimetry system for use in spot-scanning PT. A calibration model accounting for LETd and local dose-rate was established by relating the measured signal for single-spot proton beams to Monte-Carlo-simulated dose maps. The model validity was tested on both monoenergetic and clinically relevant fields.

Materials and methods

Dosimeter fabrication

Radiochromic silicone-based dosimeters were produced following the procedure established by Høye et al. [8]. The chemical composition in weight percentage (wt%) was 5.1 wt% curing agent and 93.2 wt% silicone elastomer (both contained in the Sylgard 184 Silicone elastomer kit, Dow Corning), 0.26 wt% leucomalachite green and 1.5 wt% chloroform (Sigma-Aldrich). After thorough mixing, the compound was placed in a low vacuum to remove bubbles, poured into two different cylindrical moulds (small: Ø 50 mm, 100 mm height and large: Ø 74 mm, 75 mm height) and left to cure completely shielded from light at room temperature for 40 h. In total, ten dosimeters were produced; two using the Ø 74 mm moulds, and eight using the Ø 50 mm, out of which five (used in the correction model) had to be cast on a separate occasion at the same site; the procedure above was strictly followed, to minimize variations between the results.

Pre-and post-irradiation read-out procedures and data reconstruction

Pre- and post-irradiation scanning of the dosimeters were performed with the commercially available Vista 16 optical CT-scanner (Modus Medical, London, Canada) using 1000 projections over a 360-degree rotation. The dosimeters were placed in a tank containing a water and glycerol mixture fine-tuned to match the dosimeter's refractive index by placing a checkerboard pattern between the dosimeter and light source and changing the glycerol concentration until minimal distortions could be observed. Additionally, blue food colouring was added to the index-matching liquid to match the dosimeter colour, thereby increasing the dynamical range of the absorption measurement, which is limited by the camera [12], [23].

Optical scanning was performed two hours pre- and post-irradiation, and data reconstruction of the optical CT slices was performed using the ordered subsets convex algorithm with regularization via total variation (OSC-TV) algorithm from the Vista 3-D Reconstruction software (Modus Medical, London, Canada), with 0.5 mm3 voxel size.

Irradiation procedure

The dosimeters were irradiated with proton beams generated by an isochronous cyclotron (Varian ProBeam Proton Therapy System) at the Danish Centre for Particle Therapy in Aarhus, Denmark. Mono-energetic (single spot) plans were delivered in beam quality assurance (QA) mode, where the energy, dose, and machine-set dose-rate were input directly, without the need of a beam-delivery plan. The five dosimeters used to parametrize the LETd and dose-rate correction model were irradiated with an incident energy of 120 MeV, delivering 18 Gy in water at the Bragg peak and with machine-set dose-rates of 56, 79, 101, 124 and 146 Gy/min at 2 cm depth in water. For the dosimeters used to test the correction, three incident energies were selected: 75, 95 and 120 MeV delivering 10, 17 and 18 Gy in water at the Bragg peak, with a machine-set dose-rate of 124 Gy/min at 2 cm depth in water. Because two batches were used, we used the dosimeters irradiated with the 120 MeV beam, with dose-rate of 124 Gy/min beam as a control, to ensure that the dosimeters had the same response. The dosimeters were Ø 50 mm, and build-up of 1–3 cm solid water slabs were used.

The plans with multiple spots were prepared in the clinical treatment planning system (Eclipse Varian Medical Systems, v. 13.7.16) using a virtual water phantom. Two structures were delineated and used as planning target volumes: a rectangular cuboid (3 × 3 × 2 cm3) and an irregular shape with a volume of 2.5 cm3. The plans were optimized with upper and lower dose constraints of 16.6 and 16.1 Gy and delivered an average of 16.5 Gy to the target volume. The plan with the cuboid target consisted of an SOBP composed of six pencil beam layers with energies between 72 and 88 MeV (total of 150 spots), and the irregular shape’s plan was composed of six pencil beam layers with energies between 72 and 91 MeV (285 spots), which will be referred to as an optimized irregular plan (OIP). Both plans were delivered to the flat surface of the Ø74 mm dosimeters, with no solid water build-up slabs. The machine-set dose-rate was defined automatically for each energy layer, so it was retrieved from irradiation log files and observed to vary between 101 and 145 Gy/min at 2 cm depth in water.

Monte Carlo simulations and data analysis

Monte Carlo (MC) simulations in TOPAS (version 3.2) [24], [25] were used to estimate dose and LETd distributions based on the treatment plans. The modules used were: g4em-standard_opt4, g4h-elastic_HP, g4h-phy_QGSP_BIC_HP, g4decay, g4ion-binarycascade and g4stopping. The beam was simulated using a phase-space model validated against ionization chamber arrays and scintillator measurements of spot sizes at different depths, parametrizing the beam energy, energy spread and the positional and angular spread [26].

It was not possible to obtain the exact atomic composition of the silicone matrix and curing agent for the dosimeter material. An estimate was therefore made based on manufacturer's information and expected chain length of polydimethylsiloxane (n = 362) [7], [13]. We assumed a negligible variation of the mass density, so the value of 1.05 g/cm3 was used for all simulations.

Dose and LETd were scored to the medium, i.e. the dosimeter material, in 0.5 mm3 voxels. To estimate local dose-rate and LETd of the plans with multiple spots, the contribution of each spot was scored individually. As mentioned in Section 2.3, the spot layers had a defined machine-set dose-rate which was retrieved from the irradiation log files. Machine set dose-rates are measured at 2 cm depth in water, so had to we corrected it to depth in the material, which has a stopping power ratio of 0.97 [27]. We wrote a shell script that initiated TOPAS for each spot, and after a run was completed, initiated a python script to estimate the dose-rate map by normalizing the integrated depth-dose to the machine-set dose-rate for that spot.

For each mono-energetic beam, 107 incident protons were scored, while 109 protons were scored for the SOBP and OIP. All primaries and secondaries particles were scored. To estimate absolute dose, a conversion table between number of Monitor units (MU) from the treatment planning system to number of protons in TOPAS was used [26]. The table was created from absolute dose measurements in our department beam line. Data analysis was implemented in Python, using Scipy, Numba, and PyElastix libraries [28], [29].

Experimentally determined quenching correction factors

For the measurements performed to determine the quenching correction factors, we defined the following quantities: 1) change in the absorption coefficient between pre-and post-irradiation optical CT scans ([Δα] = cm−1), 2) dose absorbed in the dosimeter ([Dmc] = Gy), which was estimated with Monte Carlo simulations, 3) local dose-rate ([ = Gy/min) estimated with Monte Carlo simulations and the machine-set dose-rate values, and 4) the measured dose to the dosimeter (Dα, [Dα] = Gy), found by multiplying α by a constant optical density-to-dose normalization factor (F, [F] = cm/Gy).

The factor F was estimated with the dosimeter irradiated with the 120 MeV beam and the lowest possible machine set dose-rate (56 Gy/min), at 1 cm depth, corresponding to a LETd of 1.03 keV μm−1. However, the preliminary analysis indicated that the Ø75 mm dosimeters had a systematically 20% higher F factor than the Ø50 mm dosimeters. Therefore, the SOBP and OIP measurements will only be reported normalized to the maximum dose.

In the absence of quenching, Dα would change linearly with Dmc for the entire proton range at the dose levels applied in this study. However, with increasing LET and dose-rate, there is a significant reduction in the ratio between measured dose (Dα) and absorbed dose (Dmc). This was modelled with a quenching correction factor (QCF) of the form:where is the local dose-rate for each spot estimated with Monte Carlo simulations. We found ξ by separating the data into high and low LETd groups, applying a 2D non-linear least squares optimization to the two groups, and finding the intersection between the curves, which yielded a value of 2.45 keV μm−1. The fit parameters were then re-estimated by a 2D non-linear least-squares optimization with a fixed boundary at ξ.

The optical measurement and the calculated dose distributions were spatially aligned with an intensity-based minimization. To minimize possible parametrization uncertainties of the beam penumbra, only voxels in a radius of 1 mm from the central axis were included in the estimation of Eq. (1). Following Christensen et al. [30], the LETd values included in the model were scored from the entrance to the 80% distal dose point.

Results

The dosimeter sensitivity estimated with the 120 MeV beams increased slightly for low LETd until 2.45 keV μm−1, decaying exponentially for higher LETd values until the distal-80% depth (Fig. 1). The fit had a mean square error of 0.03, and R2 = 0.92, with the fitting parameters shown in Table 1.

Fig. 1

Measured (circle markers) and Monte Carlo simulated (solid black line) depth-dose curves for the dosimeters irradiated with a 120 MeV proton beam and machine-set dose-rates of 56–145 Gy/min are shown in a), with the LETd curve shown in red. The data points represent the mean value of voxels in a 1 mm radius from the central beam axis, with the colours showing the different machine-set dose-rates of the five measurements. In b), the measured quenching factor is shown alongside the model fit (Eq. (2)) with the same colours as a). Note that, while the machine-set dose-rate was fixed for each measurement, the local dose-rate used in the model changed with depth. Therefore, the fit to the data used the varying dose-rate, instead of the machine-set dose-rate indicated by the colour. Figure c) shows the 2D quenching curve as a function of LETd and local dose-rate as a heat map with contour lines. The curve follows Eq. (2) with η = 1, and for higher LETd, η is set to −1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

Fitting parameters for the quenching correction curve (Eq. (2)).

μ	(a±σ)×10^-1	(b±σ)×10^-1	(c±σ)×10^-3	(d±σ)
1	0.32 ± 0.12	−4.22 ± 1.16	6.63 ± 2.11	−0.81 ± 0.11
−1	1.06 ± 0.05	−0.12 ± 0.24	7.51 ± 1.92	−1.13 ± 0.07

The corrected measurements of the monoenergetic beams had gamma pass-rates of >96% using 3%-3 mm thresholds, and >90% in 2%-2 mm gamma analysis, using a lower threshold of 10% maximum simulated dose (Fig. 2, Table 2). The mean deviations in a voxel-by-voxel comparison between the corrected measurements and Monte Carlo simulations were 0.9 ± 1.6% for the 120 MeV beam, 1.7 ± 2.4% for the 95 MeV beam and 5.6 ± 1.8% for the 75 MeV beam (percentages are relative to the maximum Monte Carlo dose).

Fig. 2

Measured, calibrated and Monte Carlo simulated depth-dose curves for the monoenergetic beams (upper panel). The Monte Carlo dose is shown in black, measured signal normalized to dose in yellow (dotted line) and corrected signal (square markers). The lower panel shows the error (in Gy) between Monte Carlo dose and corrected signal. Lines corresponding to ±0.5 Gy (dotted-dash lines) and ±1 Gy (dotted lines) were included. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2

Summary of the dosimeters; simulated, measured, and corrected dose maps. The Machine-set dose-rate for each energy was estimated from irradiation log files and corresponded to the dose-rate at 2 cm depth in water. The dose reported indicated whether absolute or relative dose maps were obtained from the given dosimeter. The dosimeters were the symbol ‘–’ appears were used to parametrize the correction model.

Dosimeter Ø [mm].	Energy [MeV]	Machine-set dose-rate [Gy/min]	Dose reported	Mean deviation ± σ [% max Dmc]	γ pass 3%-3 mm	γ pass 2%-2 mm
50	120	56	–	–	–	–
50	120	78	–	–	–	–
50	120	101	–	–	–	–
50	120	124	–	–	–	–
50	120	145	–	–	–	–
50	75	124	Absolute	5.6 ± 1.8	97%	92%
50	95	124	Absolute	1.7 ± 2.4	95%	91%
50	120	124	Absolute	1.9 ± 1.6	98%	94%
75	SOBP	101–145	Relative	3.7 ± 4.3	96%	90%
75	OIP	101–145	Relative	4.4 ± 2.6	94%	89%

After quenching correction, the SOBP showed a 96% pass-rate with a 3%-3 mm gamma analysis, and 90% with a stricter 2%-2 mm constraint, and a mean dose error of 3.7 ± 4.3% (% of maximum Monte Carlo dose) in a voxel-by-voxel comparison (Fig. 3, Table 2). The OIP field showed a 94% and 89% pass-rates for a 3%-3 mm and 2%-2 mm gamma analysis, and a mean voxel-by-voxel dose error of 4.4 ± 2.6% (Fig. 4, Table 2).

Fig. 3

Cross-sections of the simulated and measured dose distributions of the SOBP field, including the lateral (panels a–d) and top-to-bottom views (panels e–h). Plots a) and e) represent the measurement (normalized to dose), b) and f) the Monte Carlo simulated dose and c) and g) the measurement after the quenching correction. In d) and f) the voxel-by-voxel deviation between corrected measurement and Monte Carlo dose can be seen.

Fig. 4

Cross-sections of the simulated and measured dose distributions of the OIP field, including the lateral (panels a–d) and top-to-bottom views (panels e–h). Plots a) and e) represent the measurement (normalized to dose), b) and f) the Monte Carlo simulated dose and c) and g) the measurement after the quenching correction. In d) and f) the voxel-by-voxel deviation between corrected measurement and Monte Carlo dose can be seen. The optical artefacts in figures a) and c) at the distal end of the field are due to the proximity to the bottom of the dosimeter, which creates distortions on the optical CT-scanner.

Discussion

In this study, we tested the feasibility of constructing a LETd and local dose-rate based correction model for radiochromic silicone-based dosimeters based on mono-energetic proton spots. We demonstrated that such a model can be used to measure relative 3D dose distributions.

Clinical proton beams present energy straggling and a non-negligible amount of dose is deposited by secondary particles [31]. Consequently, the LET in each point is described by a spectrum, instead of a single quantity. Dose- and fluence-averaged LET have been successfully implemented in quenching correction models for a variety of other dosimetric systems: the study by Robertson et al. used the Birks model for scintillator light emission and LET calculations with Monte Carlo to correct 2D scintillator measurements [19], Resch et al. used LETd as a descriptor of beam quality for radiochromic films [32], and Herrmann et al. used track structure theory to create a correction model for alanine dosimeters [33]. The QCF curve of our study was similar to relative effectiveness curves for other dosimeters estimated with track structure theory [20], [21], [34]: a nearly constant dosimeter response for low LET and rapid decrease for higher values. The agreement between MC and corrected measurements was poorer for points after the distal-80%; this could be due to mistakes in the alignment, uncertainties in the dosimeter composition, or in the MC beam model and the QCF parametrization. As mentioned in Section 2.5, LETd and dose-rate values after the distal-80% dose depth were not considered in the parametrization of the QCF curve for a similar reason; the steep dose gradient would make it more susceptible to uncertainties in the MC beam model and positioning errors between simulated and measured dose distributions [26], [35], [36], [37], [38]. In future experiments, positioning uncertainty could be minimized by placing markers in the dosimeters, which appear in both cone-beam and optical CT scans. Uncertainties in the material composition could also be removed by measuring the dosimeter in a mass spectrometer.

The dose-rate dependency was non-negligible in the clinical range, exemplified by the more complex plans included in this study, where there would be an additional 3% deviation in the measured dose if the dose-rate had not been parametrized. Ideally, the signal would be normalized to the dosimeter response to photons instead of the low-LET region of the 120 MeV beam. This was not possible, as medical linear accelerators deliver dose at a much lower dose-rate than protons, and the dosimeter response depends on dose-rate.

The large dosimeters presented a 20% systematic increase in dose-response, although the same quenching in signal was observed when compared to the smaller dosimeters, meaning that the quenching correction model was still valid. Several parameters can affect the sensitivity of radiochromic dosimeters, including chemical composition, light exposure, temperature variations, the timing of read-out and the dosimeter size [39]. The former four should not have a large effect, as the difference in response was observed in dosimeters produced in the same batch, stored, and transported in the same container, and time between irradiation and read-out was kept constant. However, size could be important because of chloroform evaporation; it acts as a sensitizer, so its concentration affects the dosimeter response [13]. Smaller dosimeters have a higher surface to volume ratio, which potentially accelerates the evaporation process from the bulk, as discussed by Wheatley et al. [7]. The size of the dosimeter could also affect the curing time, which our preliminary studies showed can significantly change the response of this dosimeter composition. This topic will be investigated in further detail in future studies.

Our results showed that LETd and local dose-rate are good quenching descriptors, also in the correction of plans with multiple spots. However, the model presented in this study has two main limitations: 1) Due to the dose-rate dependency, several dosimeters are required to parametrize the correction model. To date, all chemical dosimeters present a dose-rate dependency to some degree [1], [9], [10], [11]. However, if the dependency could be reduced to an acceptable level in the clinical dose-rate range for protons, a correction model could potentially be established using a dosimeter irradiated with a single proton beam. 2) A different dose-response was observed depending on the dosimeter size. Similar problems have been reported in studies using cuvettes to estimate dose-response curves for larger dosimeters [1]. For that reason, 3D measurements using chemical dosimeters often do not report absolute dose values.

Although two batches were used in this study, the general applicability of the model for different batches still must be investigated carefully. The dose-response of radiochromic dosimeters depends on different factors which can be controlled with varying degrees of difficulty; such as small variations in the chemical composition, curing time and sample temperature history – both during storage and irradiation [39].

The practicality of 3D dosimetry and substantial workload of the correction model create barriers for routine clinical usage. However, given its 3D readout, high spatial resolution, flexibility, and the possibility of being cast in anthropomorphic shapes, the dosimeter formulation investigated still has potential to become a complementary tool for relative dosimetry in patient- and protocol-specific verification.

In conclusion, we presented an empirical quenching correction method for radiochromic silicone-based 3D dosimeters irradiated with spot-scanning proton beams, showing that verification of complex dose distributions using this system is feasible.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.