Physical and biological beam modeling for carbon beam scanning at Osaka Heavy Ion Therapy Center

Abstract We have developed physical and biological beam modeling for carbon scanning therapy at the Osaka Heavy Ion Therapy Center (Osaka HIMAK). Carbon beam scanning irradiation is based on continuous carbon beam scanning, which adopts hybrid energy changes using both accelerator energy changes and binary range shifters in the nozzles. The physical dose calculation is based on a triple Gaussian pencil‐beam algorithm, and we thus developed a beam modeling method using dose measurements and Monte Carlo simulation for the triple Gaussian. We exploited a biological model based on a conventional linear‐quadratic (LQ) model and the photon equivalent dose, without considering the dose dependency of the relative biological effectiveness (RBE), to fully comply with the carbon passive dose distribution using a ridge filter. We extended a passive ridge‐filter design method, in which carbon and helium LQ parameters are applied to carbon and fragment isotopes, respectively, to carbon scanning treatment. We then obtained radiation quality data, such as the linear energy transfer (LET) and LQ parameters, by Monte Carlo simulation. The physical dose was verified to agree with measurements to within ±2% for various patterns of volume irradiation. Furthermore, the RBE in the middle of a spread‐out Bragg peak (SOBP) reproduced that from passive dose distribution results to within ±1.5%. The developed carbon beam modeling and dose calculation program was successfully applied in clinical use at Osaka HIMAK.

equipped with a scanning system into clinical use. The Kanagawa facility has also started carbon therapy based on scanning irradiation.
In addition, the basic design and construction of the Osaka Heavy Ion Therapy Center (Osaka HIMAK) started in 2014, and clinical commissioning, including beam modeling for treatment planning, started in 2018.
Carbon scanning treatment planning requires precise beam modeling with absolute dose calculation to determine the number of carbon particles irradiated to every spot. Moreover, a carbon beam has high biological effect around the Bragg peak, so cell killing of tumors and normal tissues must be controlled by assuming an appropriate biological model. In starting carbon scanning treatment at NIRS, a biological model based on a microdosimetric kinetic model (MKM) was proposed. 1,2 In addition to the MKM, an updated clinical dose system adopting a carbon beam as standard radiation and considering the dose dependency of the relative biological effectiveness (RBE) was also introduced. 3 In contrast with the above trend, this study is rather conservative with respect to the biological model. We instead use a well-proven biological model established by passive carbon irradiation, which is based on a conventional linear-quadratic (LQ) model and the theory of mixed radiation fields with different radiation quality. 4,5 This approach fully complies with the carbon dose distribution accumulated so far for passive irradiation, on which a dose fractionation and escalation study was based. 6 The biological and clinical doses for normal tissue vary significantly by introducing the RBE dose dependency, however, and this means that the normal tissue complication data accumulated in past carbon broad-beam therapy can no longer be used. In starting carbon scanning therapy at Osaka HIMAK, we thus consider it too early to introduce the RBE dose dependency, when a comprehensive conclusion as to the dose level for normal tissue complication has not yet been achieved. Therefore, we still use the photon equivalent dose to describe the biological and clinical doses, and we fix the RBE at 10% HSG survival without dose dependency to evaluate the clinical dose for normal tissue in a compatible way.
The aim of this study is to develop physical and biological dose modeling for carbon scanning treatment planning. The treatment planning software is VQA Plan (Hitachi, Ltd.). To apply this software to clinical use, we registered appropriate carbon beam scanning data reflecting the characteristics of the irradiation nozzle, and radiation quality data such as the linear energy transfer (LET) needed to calculate biological effects. Then, the biological dose optimization and dose calculation function were implemented in VQA Plan.  The facility has three treatment rooms, of which two have two fixed irradiation ports with vertical and horizontal incidence, while one has two fixed ports with horizontal and oblique (45°) incidence. A carbon beam accelerated to the range of 100 to 430 MeV/n by a synchrotron accelerator is transported to the irradiation nozzles and scanned in the x and y directions by scanning magnets. A vacuum chamber extends upstream of the beam monitors, and its distance from the isocenter is about 1.4 m. Of the three installed beam monitors, two are ionization chambers to measure the irradiated particle number of the carbon beam, and one is a multi-wire proportional counter to detect the beam position in the x and y directions. The carbon scanning beams are prepared at a 3-mm range interval, and the beam energy is changed by adjusting both the accelerator energy and the binary range shifters in the nozzles. The pristine carbon Bragg peaks are broadened by a ripple filter created by a 3D printing technique. The accelerator provides 12 energies from 100 to 430 MeV/n, and the range shifters provide beams of intermediate energy. The irradiation scheme uses continuous beam scanning, in which the carbon beam is continuously extracted from the synchrotron and irradiated spot by spot without turning the beam off and on. The beam intensity can be adjusted to a constant ranging from 1 to 10 MU/s by tuning the extraction RF power of the synchrotron to specify the intermediate dose between spots when the continuous beam is moved from one spot to the next. Note that MU stands for "monitor unit," which corresponds to the particle number of the carbon beam and is defined later in the text.
The dose calculation algorithm for carbon treatment planning is based on a triple Gaussian model 7 : where d i denotes the physical dose at calculation point (or voxel) i, and d ij is the dose contribution from beam j to point i, and w j is the , d ij (2) , and d ij (3) are calculated by multiplying the integral depth dose (IDD) by the lateral Gaussian distribution: Here, IDD n (z) denotes the IDD of the n-th component, and G n (z) is the lateral two-dimensional Gaussian distribution at depth z: where σ 1,x (z), σ 1,y (z), σ n (z) (n = 2,3) are the beam sizes at depth z, (x i , y i ) is the lateral position of point i and, (x j (z), y j (z)) is the center position of beam j at depth z. The first component in Eqs. (1), (2), and (3) corresponds to the incident 12 (4), IDD n (z) (n = 1,2,3) in the depth direction and the beam sizes σ 1,x (z), σ 1,y (z), σ n (z) (n = 2,3) in the lateral direction are parameters to be determined in the physical beam-modeling process.
Here, we introduce MU units to describe the irradiated particle number of the carbon beam, as follows. In the ionization dose monitor, electrodes collect ion pairs generated by the carbon beam's passage through the air, and then electric circuits convert the accumulated charge into a pulse signal. We introduce the gain G of the dose monitor, ie the number of ion pairs generated by one incident carbon particle, as where 1/dE/dx| air denotes the stopping power of 12 C in air, g is the dose monitor's air gap, d is the density of air, and W denotes the air's W-value to generate one ion pair (35.1 eV). Figure 2a shows the design value of the gain G for the case of a 1-cm electrode gap g in the monitor. We chose an ionization amount of 25 nC as 1 MU so that about 100 MU would correspond to irradiating a physical dose of 1 Gy in a 1-liter volume, as described in 8. In this definition, the calibration factor K(E) of the dose monitor, 9 ie the number of carbon particles per MU, is calculated as where e denotes the elementary charge (1.602 × 10 -19 [C]). Figure 2b shows the calculated design value of K(E). In the figure, 1 MU corre-

2.A.2 | IDD data
We used both dose measurement results and Monte Carlo simulation to model IDD n (z) (n = 1,2,3) data. The Geant4.9.3 10 platform provided the Monte Carlo calculation code, and the irradiation apparatus shown in Figure 1 was implemented accordingly in the simulation. The IDD components of the incident 12 C particles, IDD 12C (z), and of the fragment particles, IDD frag (z), were then determined through both the dose measurements and the Monte Carlo simulation. The calculated results from the simulation were corrected to reproduce the measured IDD data from a large-area ionization chamber, the Stingray (IBA Dosimetry), whose detection diameter is 12 cm. Because the total cross section in the Geant4 Monte Carlo code was assumed to differ from the true value, the calculated IDD did not agree with the Stingray measurement results, so we developed correction methods in terms of the total cross section.
First, we express IDD 12C , ie the dose contribution from 12 C, by multiplying the 12 C fluence distribution Φ(z) at depth z and the convolution of the stopping power S(r) in terms of the range straggling σ strag : where R 0 denotes the initial residual range, S(r) is the stopping power of 12 C for residual range r, and λ is a quantity related to the total cross section. The 12 C particle fluence Φ(z; λ) is then reduced, according to the depth z, by nuclear reaction: where Φ 0 denotes the initial particle fluence, N A is the Avogadro number, is the density of the material, M is the molecular weight, and σ is the total cross section. In Eq. (8), the quantity λ is defined meaning that λ is inversely proportional to σ. Note that, in the convolution calculation in Eq. (7), the stopping power S(r) is tabulated from an ICRU73 data table, and the straggling parameter σ strag is kept constant independently of the depth. By tuning λ and σ strag , we can use Eq. (7) to well reproduce the relative distribution of Geant4 IDD 12C (z) scored results. Correction for IDD 12C (z) only pertains to the fluence term Φ(z; λ), and thus, the correction formula for IDD 12C (z) is derived as where λ and λ' are adjustment parameters for IDD 12C (z) before and after correction respectively. The fragment component IDD frag (z) is | 79 assumed to scale with the total cross section; that is, we assume a linear relationship for IDD frag (z) with the total cross section, giving the following correction formula: Figure 3 shows the results of IDD correction for a 430-MeV/n carbon beam obtained by the above methods, together with corrected Monte Carlo results reproducing the measured result within AE2% in the proximal and distal regions.
After adjusting the relative IDD shape to the dose measurements, the IDD data are converted to absolute dose values reflecting the absolute dose measurements at a shallow depth (2 cm), as described in 11. The dose area product (DAP) at a 2-cm depth was measured by rectangular uniform grid irradiation using a monoenergy beam 12 . The DAP measurement was performed with an advanced Markus chamber, and the absolute dose was measured at the center of the field where lateral dose distribution is flat. The irradiation number of carbon particles is uniform, with a value of Q [MU] per spot, and the spot spacing is 3 mm in both the x and y directions. Then, the DAP per MU is defined as where D meas denotes the measured dose in units of Gy at the center of the field. We introduce the total IDD as Next, using the corrected relative IDD 12C (z) and IDD frag (z) in Eqs.
The dimensions of the DAP and IDD in Eqs. (12) and (14) where ϵ denotes the detection efficiency of the Stingray chamber as explained later, in Section 2.A.4.

2.A.3 | Beam size of 12 C
The first component in Eqs. (2) and (3) corresponds to the incident 12 C particles. Its beam sizes σ 1,x (z), σ 1,y (z) are calculated by considering both the beam's optical parameters and multiple Coulomb scattering. The beam transport equations including multiple scattering according to Fermi-Eyges theory are.
where σ 11 (z), σ 12 (z), σ 22 (z) are the phase-space parameters at depth z, θ 2 is the increment of the angular divergence when passing through distance z, and z = 0 denotes the initial point of transportation. The beam size is calculated by the square root of σ 11 (z) as Equations (17), (18), (19), and (20) are expressions in the x direction, and the same equations also hold in the y direction. The formula for the angular divergence increment was improved to a suitable form for treatment planning by Kanemastsu 13 : (a) Dose monitor gain G with a gap size of 1 cm, and (b) calibration factor K(E) (number of carbon particles per MU).
where R and R' are the initial and residual ranges, respectively, Z is the atomic number, A is the mass number of the incident particles, and C is an adjustment parameter for the scattering power.

2.A.4 | Determination of fragment components
The second and third components' parameters, ie IDD 2 (z), IDD 3 (z), (2) and (4)  one is field size factor (FSF) measurement, as described by Zhu 11 using various sizes of square fields, and the other is frame pattern irradiation, as described by Inaniwa 15 . In both cases, the absolute dose at the center of the field is measured with a pinpoint chamber.
We chose the frame pattern irradiation method, because the primary component with a large dose contribution is absent in the detection area, so we expected to be able to model the second and third components more accurately. We exploited the seven irradiation frame patterns listed from A to G in Table 1. Figure 6 shows the spot con-  (10) and (11) are λ = 170 mm (before correction) and λ'=150 mm (after correction).
fitting of IDD 2 (z), IDD 3 (z), σ 2 (z), σ 3 (z) was performed using Eq. (1) by setting i as the origin, ie the center of the field, and summing the dose contributions from beam j according to the frame pattern con- Table 1.
where R denotes the radius of the Stingray chamber's detection area We again fitted parameters IDD 2 (z), IDD 3 (z), σ 2 (z) with the condition of σ 3 = 25 mm. Figure 8 shows the fitted results for σ 2 (z) with various energies (Table 2) where z p represents the carbon Bragg-peak depth, and the parameter values are C 1 = 5 [mm], C 2 = 30 [mm], C 3 = 0.006. Next, Figure 9 shows the fitted results for σ 2 (z; L) when inserting range shifters with shifters of water equivalent thickness L, is calculated by  Figure 10 shows the fitted results for IDD 3 (z). As seen here, IDD 3 (z) increases up to the Bragg-peak depth and then decreases at depths exceeding the Bragg peak. This tendency holds regardless of the incident energy. Therefore, we approximate IDD 3 (z) in the same exponential function form as where z p denotes the Bragg-peak depth, and the parameters are where IDD 3 (z; L) denotes the third component with range-shifter of water equivalent thickness L.

2.A.5 | Determination of correction factor (IDNF)
Here, we define the integral dose normalization factor (IDNF), which is a correction factor for each IDD to compensate for the difference from the measured absolute dose. Our IDNF factor works like the depth-dose normalization factor (DDNT) in 11. We determine the IDNF by comparing the absolute calculated dose value with absolute dose measurement results for several volume irradiations. We introduce this factor to correct the ambiguity of the absolute dose for the DAP at a 2-cm depth.
where i denotes the calculation point or voxel, S i is the cell survival, where α ij and β ij are the LQ parameter contributions from beam j to point i. The physical dose contribution d ij has a Gaussian distribution, as in Eq. (2), in the lateral direction, while the LQ parameters α ij and β ij are assumed to have no lateral distribution and are only functions of the LET value at the specified depth. Next, we introduce a biological effect at point i as Then, the number of particles of beam j, w j , is optimized to make the biological effect e i match a goal value E i . An evaluation function is calculated by summing the squares of the residuals from the goal value over all points: The clinical RBE is thus defined as the ratio of the clinical and physical doses: The above description is the established theory for passive irradiation since the start of carbon therapy in Japan. Recently, this model has been improved to use two kinds of LQ parameters, one for carbon and the other for helium. 16 where d ij (c) and d ij (frag) are the physical dose contributions of the carbon isotopes and fragment isotopes respectively. Figure 11 shows the carbon and helium LQ parameters α (C) , β (C) , α (He) , β (He) used in this study, with the carbon LQ parameters α (C) and β (C) in red and the helium LQ parameters α (He) and β (He) in blue. These LQ parameters are functions of the LET, so we prepared dose-averaged LET (LETd) data for each scanning beam from the Geant4 calculation results.

2.B.2 | Beam data for biological dose
The physical dose contributions of the carbon and fragment iso- , d ij (frag) in Eqs. (34) and (35), are related to the triple Gaussian components in the following way: where IDD frag (z) denotes the sum of IDD 2 (z) and IDD 3 (z), G 2 (z) is the lateral Gaussian distribution of the second component as in Eq. (4), and R(z) is the ratio of the IDD frag values including and excluding carbon isotopes. Here, we introduce R(z) because the second component of the triple Gaussian model is supposed to be a mixture of carbon isotopes such as 11 C and other fragment isotopes from Z = 1 to Z = 5. This ratio R(z) is calculated by Geant4 as where the IDD frag (z; scored except Z = 6) are scored without Z = 6 isotopes, and the IDD frag (z; scored except 12 C) are scored without incident 12 C particles. The numerator IDD frag (z; scored except Z = 6) thus excludes the effect of carbon fragment isotopes such as 11 C, while the denominator IDD frag (z; scored except 12 C) includes such effects. Therefore, R(z)×IDD frag (z) corresponds to the fragment contribution from Z = 1 to Z = 5, while (1-R(z))×IDD frag (z) corresponds to the component for carbon fragment isotopes such as 11 C. Table 4 summarizes the relationship between the triple Gaussian model and the applied LQ parameters. Figure 12 shows the Geant4 calculation results for the ratio of IDD frag without carbon isotopes (red line) to Next, the LETd data for obtaining the LQ parameters is also prepared by Geant4 according to where k denotes the k-th event in Geant4, S(E) is the stopping power, and d k and E k are the dose and energy deposit respectively.
We use MSTAR 17 Figure 11 via the dose-averaged LET in Eq. (39). In addition, the LQ parameter α for carbon α (C) in Eq. (34) is directly calculated by Geant4 to consider a broad LET spectrum around the Bragg-peak region as explained in Inaniwa 3 : where event k is scored involving only Z = 6 particles. The doseaveraged α (C) for carbon is directly scored and prepared as a function of depth z, as shown in Figure 14.

2.C | Cell experiments
To verify the calculation accuracy of the physical and biological doses described above, we performed cell survival measurements under a biologically optimized spread-out Bragg peak (SOBP) condition. The methods of these cell experiments were as follows. in a humidified atmosphere of 5% CO 2 , at less than 90% confluence.

| RESULTS
The appendix summarizes the beam data, such as the IDD, LET, and beam size, described in Section 2. The dose calculation accuracy was verified by comparing the results of dose measurement and past publication results described in 16.

3.A | Physical dose
The calculated absolute physical dose was compared to absolute dose measurements after tuning the IDNF parameters. An SOBP was created by applying treatment planning software with the registered beam data, and dose measurements were performed with a

| 87
Markus ionization chamber at the center of the SOBP. Figure 15 shows the results of IDNF correction and the final dose calculation accuracy. In the figure, the horizontal axis represents the range in water, while the vertical axis represents (a) the IDNF value and (b) the difference between the absolute calculated and measured doses.
The final dose calculation accuracy was verified under the irradiation conditions listed in Table 3, and Figure 15b shows the field size dependency as a function of depth. The IDNF correction is within AE1.5%, and the final dose calculation accuracy is within AE2% in the range of 4 to 30 cm. We will report verification results in case of heterogeneity in a later submission. 19 3.B | Biological dose  We next performed the HSG cell irradiation experiment and measured the cell survival in the SOBP region, as shown in Fig. 17.
Measurements were taken at three depth points-distal, middle, and proximal-in the 6-cm SOBP flat region with a range of 15 cm (Fig. 16a), and the observed survival fraction was uniform, as shown in Fig. 17a. The dose dependency of the survival fraction was also measured at the middle of the 6-cm SOBP, as shown in Fig. 17b. Instead, we could directly score α (He) in the Geant4 simulation as a function of depth like α (C) , but we consider this discrepancy tolerable. This is because the fragment α (He) does not affect the RBE in the SOBP, and overestimation in the tail region is on the side of safety for normal, nontumorous tissue.
In this study, we chose a biological model in which the biolog- where D oar is the dose constraint. After dose calculation under the 10% survival condition, this scaling factor is again applied to calculate the clinical dose, physical dose, and the number of particles at each spot, w j : where w j (S = 10%) and d i (S = 10%) are the number of particles of beam j and the physical dose at point i, respectively, under the con-  of β-term. The same scheme of multifield irradiation is possible in scanning irradiation because the biological model is the same.

D A T A A V A I L A B I L I T Y S T A T E M E N T
Research data are not shared, because the data of this study are obtained from openly available Monte-Carlo calculation programs and commercially available dosimetric equipments and we think sufficient reproducibility is ensured.

APPEN DIX
This appendix gives the registered beam data used in this study. Figures A1, A2, and A3 show the IDD data, LET data and LQ parameters, and beam size data respectively.