2.1 Target setup

A spherical target with a diameter of 84 mm was placed as the clinical target volume (CTV) in a numerical water phantom at a center depth of 130 mm (Fig. 1a). In this study, a simple water-equivalent system was adopted for the primary purpose of comparison with simulations and measurements, assuming clinical application to bone and soft tissue or liver tumors with densities similar to water in and around tumors. The CTV motion directions were assumed to be in the lateral and/or longitudinal (proximal, beam upstream) to investigate the motion effect. If the CTV motion occurred simultaneously in both directions, an elliptical CTV motion was performed with a phase difference of π/4. The amount of the CTV motion within the planning target volume (PTV) was set to a combination of 0–5, 7, and 20 mm physical length in the lateral direction and 0–5-mm physical length [equivalent to water-equivalent path length (WEPL)] in the proximal direction. Table 1 lists the experimental conditions for the CTV setting and irradiation.

Fig. 1

a Schematic diagram of the numerical phantom for the spherical target and the measured plane. The target moves in the lateral and/or proximal directions. b Schematic and c photograph of the measurement setup. The 2D array was positioned to measure the dose through the target center. White arrows in (a) and (b) correspond to the arrows in Fig. 3a and b

Table 1 Experimental conditions for the respiratory motion of the CTV2.2 Treatment planning

The XiO-N treatment planning system (ELEKTA, Stockholm, Kingdom of Sweden and Mitsubishi Electric, Tokyo, Japan) [7] was used to derive plan parameters for both experiments and simulations. The treatment plans in XiO-N were created for PTVs inclusive of internal margins (IM) [7], although they were assumed to be in a static condition. A 3 mm setup margin (SM) [7] was applied in all directions and adjustments were made to beam range and leaf margins to ensure sufficient dose coverage within the PTV in the static condition. For the respiratory-gated irradiation, the IM for a gating window of 30% lateral motion was provided in each direction. In this paper, the CTV motion amount within the gating window is referred to as “CTV motion amount during irradiation”.

2.3 Measurement setup

Irradiation was performed under static and moving conditions described in Sect. 2.1, based on the planned data described in Sect. 2.2. The treatment plans with the prescribed dose of 5 Gy (relative biological effectiveness, RBE) to the PTV were used. The respiratory cycle was set to 3.5 s, offset from the 3 s of the synchrotron operating cycle. The physical dose distributions were measured at least twice for each condition using a cross-calibrated 2D array [OCTAVIUS detector 729XDR (T10031), PTW, Freiburg, Germany] with carbon-ion beams. The plane physical dose distributions were measured by placing the 2D array at a 30-degree angle (Fig. 1b and c). A respiratory motion phantom [KMF-GU-001, Accelerator Engineering Corporation (AEC), Chiba, Japan] was used to simulate the respiratory motion. The lateral CTV motion was performed by shifting the moving stage of the phantom. The longitudinal CTV motion was simulated using a pair of acrylic wedges moving in the lateral direction to cause beam range variation. A laser-type synchronizer, AZ-733VI (Anzai Medical, Tokyo, Japan), was used for respiratory-gated irradiation. The synchronizer was placed to enable the laser to detect the surface of the acrylic block and ensure that the laser was projected in the lateral direction.

2.4 Simulation of 3D physical and clinical dose distributions for layer-stacking irradiation with respiratory motion

In-house C++ software was developed to simulate the 3D physical and clinical dose distributions including the target motion. The simulation utilized the broad-beam method [20, 21] and incorporated simplified wobbling effects and beam scattering, as described later in this section to closely replicate the conditions in XiO-N. The plan parameters, such as MLC, range shifter, range compensator, weighting factor for each layer, tissue–phantom ratio (TPR), and linear-quadratic (LQ) model parameters α and β for RBE calculations [22], were obtained from the corresponding plan made by XiO-N for the conditions in Table 1.

The dose distributions were simulated on a 1-mm grid within a 3D water volume. The parameters of the set values are summarized in Table 2. The temporal relationships between several parameters are shown in Fig. 2. The 3D physical dose distribution \(D\left(\textbf}\right)\) was derived as follows:

$$D\left(\textbf} \right) = \sum_\sum__d\left(}+\Delta \textbf}}_\right),$$

(1)

where \(i\) was the layer number and \(j\) was the elapsed time element within a single layer \(i\). \(d\left(\textbf}\right)\) represents the depth-dose distribution for a single layer, obtained from the TPR, the source to axis distance (SAD), and the depth component at position \(\boldsymbol\). Due to respiratory motion, the coordinate \(\boldsymbol\) receives the dose at the coordinate \(\boldsymbol+\Delta \boldsymbol_\). The weighting factor \(_\) was calculated as follows:

$$_ = \frac_}_}}\times \frac_}_}},$$

(2)

where \(_\) and \(\dot\) were prescribed dose and average dose rate throughout the synchrotron cycle, respectively. The beam was extracted at the extraction time \(_\) out of the synchrotron cycle \(_\). The denominator in Eq. (2) represents the total time required to irradiate all layers. \(\Delta _\) represents the beam-on time period within the time interval \(\Delta t\) at elapsed time \(_\) for layer \(i\). \(\Delta t\) is the calculation unit time interval and is expressed as \(\Delta t=T/_\) using the respiratory cycle \(T\) and the number of time divisions \(_\). In this simulation, we set \(_\) = 20. If the beam on/off occurred within the time interval \(\Delta t\), the net beam-on time period was derived as \(\Delta _\). After the completion of irradiation in each layer, the beam was not extracted for 1 s while the range shifter switching. For gated irradiation, the beam was not irradiated during the time period when the CTV was positioned outside the gating window. The weighting factor \(_\) for beam-on time period \(\Delta _\) was accumulated and when it reached a value equal to the single-layer weighting factor \(_\), the irradiation of the layer \(i\) was completed. The relationship between \(_\) and \(_\) was given by the following equation:

Table 2 Simulation conditions for the respiratory motion of the targetFig. 2

Timing relationship between respiratory motion and irradiation devices

Assuming a horizontal beam (gantry angle = 90°) and defining the 3D coordinates as \(}=\left(x,y,z\right)=\) (posterior-anterior, lateral (inferior-superior), depth (left–right)) direction, the coordinate \(+\Delta\boldsymbol}_\) and its respective elements were determined as follows:

$$+\Delta\boldsymbol}_=_, _,_\right)},$$

(4)

$$_=_+_+B\left(x+_, _\right).$$

(7)

Range shifter thickness was calculated as \(_=_+\left(i-1\right)\times 2.5\) mm-WEPL and inserted after the end of irradiation in each layer to adjust the beam range. \(_\) is the initial thickness at \(i\) = 1. \(B\left(x, \text\right)\) is the thickness of the range compensator (bolus) in the coordinates of the beam vertical plane. \(_\) and \(_\) are the lateral and longitudinal amounts of the respiratory motion relative to the elapsed time, expressed as follows:

$$_=\frac_}\left\\left(\frac_}+_\right)\right\},$$

(8)

$$_=\frac_}\left\\left(\frac_}+\frac+_\right)\right\},$$

(9)

where \(_\) and \(_\) are the amplitudes from the maximum exhalation to the maximum inhalation. When the CTV moved in two directions simultaneously, the phase difference in the proximal direction relative to the lateral direction was set to \(\pi /4\). \(_\) is the initial phase of the respiratory motion with respect to the synchrotron cycle, and seven patterns were verified. The initial phase changes caused a shift in the timing of beam irradiation. The elapsed time \(_\) represents the dose-calculation time point set at the center of the time interval \(\Delta t\)

$$_ = __}+j\Delta t,$$

(10)

where \(_\) is the final time element at the layer number \(i\). We defined \(_} = -\Delta t\), because the starting time point was set as \(_} = 0\).

The RBE and clinical dose \(C\left(\boldsymbol\right)\) were calculated as follows, based on the formulated clinical RBE [4] and the range-modulated beam method [22]:

$$RBE}= \frac_\left(\boldsymbol\right)}_\left(\boldsymbol\right)}^-4_\left(\boldsymbol\right)\text0.1}-_\left(\boldsymbol\right)},$$

(11)

$$C\left(\boldsymbol\right) = RBE\left(\boldsymbol\right)D\left(\boldsymbol\right).$$

(12)

The clinical factor 1.46 employed in the Gunma University Heavy Ion Medical Center (GHMC) was used. \(_\left(\boldsymbol\right)\) and \(_\left(\boldsymbol\right)\) are the coefficients of the LQ model for the mixed beam obtained as dose averages including the motion as follows:

$$_\left(\boldsymbol\right)= \sum_\sum_f\left(}+\Delta }}_\right)\alpha \left(}+\Delta }}_\right),$$

(13)

$$_\left(\boldsymbol\right)}= \sum_\sum_f\left(+\Delta \boldsymbol}_\right)\sqrt+\Delta \boldsymbol}_\right)}},$$

(14)

where \(\alpha (z)\) and \(\beta (z)\) are coefficients of the LQ model obtained from XiO-N corresponding to the depth-dose distribution \(d(z)\) in a single layer. \(f\left(+\Delta \boldsymbol}_\right)\) is given by the following equation:

$$f\left(+\Delta \boldsymbol}_\right)= \frac_d\left(\boldsymbol}_\right)}\right)}.$$

(15)

Parallel beams were used to simplify the simulations due to the fact that the SAD was approximately 7 m, which was considered sufficiently long. An approximated factor that accounts for the slight increase in the lateral dose distribution near the wobbling radius was employed and applied at all depths. The factor was defined as the ratio of the planned value of the lateral dose distribution by XiO-N to the simulated one before the wobbling correction at the isocenter depth. To simply account for multiple-beam scattering, the lateral dose distributions were convoluted using a Gaussian distribution with a standard deviation (SD) of 3.8 mm, which was determined to best fit the planned and simulated physical dose distributions on the isocenter plane. This convolution was carried out for the simulated physical and clinical dose distributions. For static conditions with combinations of the lateral and longitudinal IMs shown in Table 1, average values of gamma passing rates (3%/3 mm) and SDs for the physical and clinical dose distributions between XiO-N and the simulation were 99.94 ± 0.04% and 99.91 ± 0.02%, respectively. The simulation reproduced the dose distribution in XiO-N in static conditions.

Gating was employed at 30% amplitude when the target moved in a single direction. When the CTV moved in two directions simultaneously, the beam was gated with respect to the lateral motion. The gating delay time between the actual motion and the gating signal was set to 0.06 s, which was used in clinical practice.

2.5 Comparison of measured and simulated physical dose distributions

For comparison with measurements, the physical dose distributions were simulated with the setting of the prescribed dose of 5 Gy (RBE). Two methods were used to evaluate the agreement between the measured and simulated physical dose distributions.

First, a gamma analysis [23] was conducted. The criterion was set to 3%/3 mm, considering the sensitive volume of the detectors and experimental setup errors. The γ value with the best agreement at each chamber position was obtained between the average of the measured physical dose distributions and the simulated ones for the various initial phases.

Second, dose coverage inside the CTV was compared between measurement and simulation. The physical doses were not constant in the SOBP for the uniformly prescribed clinical dose. Therefore, the dose variation from the static condition was examined to evaluate the dose coverage for the CTV motion. The simulated physical dose distributions corresponding to the area measured by the 2D array were extracted. The averaged physical doses on the incident surface of the sensitive volume in each ionization chamber were obtained. Finally, the delta 2D conformity index (ΔCI2D), defined as the ratio of the number of ionization chambers within a ± 5% dose difference in the CTV area, was estimated.

2.6 Clinical dose evaluation

The clinical dose distributions were simulated at 4 and 10 Gy (RBE). The respiratory cycle was examined within a range of 2–10 s to investigate the effect on dose uniformity. Additionally, to examine more motion conditions, we defined the variable CTV (CTVV) created by inversely subtracting IMs from the fixed PTV [6], with 20-mm lateral and 5-mm proximal margins to the original CTV presented in Table 1. As a result, the CTVV was adjusted according to changes in the CTVV motion amounts. The CTVV offered the advantage of avoiding the necessity to create individual treatment plans for various respiratory motion, because the single fixed PTV was used.

The acceptable target motion amounts during irradiation were derived from the evaluation of the clinical dose uniformity. The conformity index (CI; defined as (V95 − V105)/CTV), was employed, because the CTV dose coverage was important when considering motion. The differences in CI values between the motion condition and the static condition, denoted as ΔCI (= CImove − CIstatic), were obtained to remove the effect of differences in the CI values in the static conditions. When the average value for all initial phases was within ± 5%, it was deemed clinically acceptable. The worst averaged ΔCI value for various respiratory cycles was defined as ΔCIworst. The maximum ΔCIworst value within ± 5% for motion conditions was regarded as ‘acceptable CTV motion amount during irradiation’. In the case of gating, this refers to the target motion amount within gating window.