Calculation of the Beam-Modulation Effect of the Lung in Carbon Ion and Proton Therapy With Deterministic Pencil Beam Algorithms

Ion beams passing through lung tissue show more pronounced energy straggling than expected for solid materials of the same thickness. Energy straggling in active scanning particle therapy can enlarge the pencil beam Bragg peaks in-depth as well as displace them, deteriorating the dose coverage of a target within the lung. While this is not yet considered in any known treatment planning system, we implement a mathematical model to be used for treatment planning, using TRiP98, which relies on a deterministic pencil beam algorithm. Through a randomization process based on a continuous Poisson probability distribution, the HU values of lung voxels are replaced with a modified value in successive iterations. The beam-modulation effect of the lung can thus be taken into account in treatment planning by recalculating the dose n times for n randomized CTs using the raster scan file of a plan that was optimized on the nonmodulated CT. The evaluation follows by averaging the resulting n dose distributions and comparing to the corresponding nonmodulated dose distribution, attending at dosimetric indices and dose-volume histograms. In this work, the functionality of these routines was tested for proton and carbon ion plans for two selected lung cancer patient cases with deep-seated tumors, showing that, with existing standard tools, it is possible to calculate the beam-modulation effect of the lung in TRiP98 in a transparent way. Variable model parameters, such as modulation power, voxel size and density voxel selection range, were evaluated. Furthermore, a systematic study for spherical geometries in a lung tissue CT cube is presented to investigate general trends.


INTRODUCTION
With the dose deposited locally in the energy-dependent narrow Bragg peak (BP) after which little-to-no-dose is given, particle therapy offers potentially more conformal treatments, as compared to photon irradiation, with improved normal tissue sparing. In particular with the pencil beam scanning method, where the lateral beam displacement is done by magnets and the range of each BP is changed by modulating the beam energy, a more conformal dose delivery, with fewer field directions, can be achieved compared to photons [1,2].
Due to the porous nature of the lung, ion beams passing through a lung (or lung substitute material) show a more pronounced energy straggling than for solid materials [3,4]. A substantial broadening of the BP shape has been found when irradiating lung tissue [5,6]. So far, current treatment planning systems (TPS) assume no additional BP broadening due to tissue heterogeneities. Ordinary planning CT images cannot capture this beam-modulation effect of the lung since the CT voxel size spatially undersamples the microscopic lung alveoli. However, the beam-modulation effect of the lung could be of clinical relevance for lung tumor treatments, which motivates this study. If this broadening effect is not taken into consideration, it can potentially lead to an undesired inhomogeneity in the dose distribution, underdosage of the planning target volume (PTV) and overdosage of the healthy surrounding normal tissue, especially in the distal region close to the PTV [3,7].
Through the use of a 3D printed binary voxel model (the voxels in the geometry are either composed by air or by the printing plastic material) the BP degradation of lung irradiation for protons has been shown to cause for initial beam energies of 140 and 200 MeV an increase of up to 60% for the 80-20% distal beam fall-off as compared to the unperturbed reference curve [5]. Using a similar binary voxel model, a mathematical method for estimating the broadening effect of a porous material has been developed [4,6,8]. The degree of this broadening depends on the average size of the microscopic structures as well as on the material density and thickness. For materials with a fine (< 500 μm) and homogeneous microscopic structure, the BP enlargement has been shown to be well-described by a modulation of the unperturbed reference Bragg curve for water with a normal distribution and its sigma as a parameter for the strength of the modulation effect.
Treatment planning studies investigating the effects on the dose uncertainties in treatment plans of lung cancer patients are limited to a few published works. Flatten et al. [9] investigated the effects of the BP degradation in proton beams by using CT phantoms with spherical targets of variable sizes placed at different depths in the lung. The lung modulation has been implemented in the TOPAS Monte Carlo tool [6] using sets of modulated DICOM images. It was shown that the underdosage of the target volume, due to the BP degradation, increases with an increasing depth of the tumor in lung and a decreasing tumor volume. In a further study by the same group [10] the effects were investigated on proton treatment plans of lung cancer patients. It was found that the maximum underdosage of the target volume due to the BP degradation was in the order of 5% while the average underdosage was roughly 2%. In another study [11] the effects on clinical treatment plans in proton beams were investigated using an analytical approach: Instead of using sets of modulated DICOM images and Monte Carlo tools, the pristine proton BPs were convolved using the MatRad toolkit [12] with several superimposed normal distributions of sigma values defined by certain modulation power values of the lung tissue. The underdosage of the target volumes was comparable to the results from Baumann [6]. However, to the best of our knowledge, there are no studies investigating the modulation effects for treatment planning with carbon ion beams. For carbon ion beams a more pronounced modulation effect is expected, because of the much sharper carbon BPs compared to protons. The lack of a carbon ion TPS for carbon ion beams that can calculate the lung modulation effect was a main motivation for this work.
To take the lung modulation effect into account directly, the CT resolution is crucial [5]. The resolution of a clinical CT image used for treatment planning is typically around 1 mm but the biological porous structure of the lung leading to the modulation is in the μm range. To work with CTs with μm resolution, assuming such images can be obtained, would lead to an impractically large data array and time increase in the handling of the DICOM data. The method presented in this work bypasses this problem.
The aim of the present work is to show that, with existing deterministic pencil beam dose calculation algorithms like the one in TRIP98 [13][14][15], it is possible to calculate the beammodulation effect of the lung tissue. Here the term, "deterministic" applies to dose calculation algorithms with analytical functions that describe a 3D pencil beam dose distribution and deliver reproducible results without random seeds. TRiP98 is a treatment planning system that was developed for research and patient treatment at the Helmholtzzentrum für Schwerionenforschung, Darmstadt, Germany (GSI). It was used to treat 440 patients during the GSI pilot project. The basic algorithms and the pencil beam model of TRIP98 were implemented into the carbon ion planning module of the Siemens Syngo PT TPS [16] that is currently used at the ion beam facilities in Heidelberg (HIT), Marburg (MIT), Shanghai (SPHIC) and Pavia (CNAO). Although we are using a deterministic pencil beam dose calculation, the results of the procedure cannot be considered entirely deterministic because the sampling of the HU values is done with random functions.
The beam-modulation effect calculation is done by implementing a mathematical model based on a continuous Poisson probability distribution within the open-source PyTRiP tool [17]. Through a randomization process, the HU value of each selected lung voxel is changed into a modified value. The lung beam-modulation effect on patient treatment plans can then be estimated by recalculating the dose n times for n randomized CTs, using raster scan files for plans optimized on nonmodulated CTs and afterward averaging the n dose distributions. Thus the result can be compared to the corresponding nonmodulated dose distribution. Using spherical PTV geometries and two selected patient cases, we evaluated a number of different model and planning parameters to show their respective relevance in terms of PTV coverage deterioration as a result of the lung modulation effect for carbon ion as well as for proton plans.
In former studies [18,19], facility-specific baseline data, including depth-dose distributions and lateral profiles, were generated with SHIELD-HIT12A [20] for both carbon ions and protons.
PyTRiP [17,21] is an open-source python package for, among other purposes, facilitating the work with TRiP98 and with TRiP native files. PyTRiP can import, export, convert and modify such files and execute TRiP98 locally or remotely. PyTRiP enables scripting large parameter studies of treatment plans and more advanced and automatized (programmable) manipulation than what a commercial TPS usually allows. PyTRiP handles the CT and planning contour files, as well as dose cubes from calculated plans and permits extensive manipulation on these objects, making it possible to randomize n CTs, recompute a plan on each of them and average all the calculated dose cubes.
The presented concept follows Baumann et al. [6] but is done in a fully automated and integrated environment and works for a standard pencil beam algorithm.

Mathematical Implementation of CT Randomization
For each investigated planning case, first a treatment plan on the unmodified planning CT is calculated, producing a raster scan file for each case. Afterward, a dose recalculation, using these raster scan files, is done n times following the mathematical method described below. This is based on the "binary voxel model" for porous targets, where the target is divided into high and low density voxels, which each traversing particle will hit with probabilities of p and (1-p), respectively, [4,6]. The probability for hitting k voxels with high density is thus given by a binomial distribution, which, for a large number of voxels, can be approximated by a normal distribution. While this normal distribution can be subsequently used to modulate the individual PBs, here we will modify the CT voxels of a certain volume of interest (VOI), in this case the ipsilateral lung VOI lung .
First the tissue density in each voxel of the CT image is calculated using the standard HU conversion table of TRiP98. Then, the voxels of VOI lung (within a user-defined density range ρ range from ρ min to ρ max ) are selected and defined as the voxels to be modulated. The modulation power P mod as defined prior is a good estimation for the strength of the lung modulation effect [4,6,[8][9][10] and given as a second user-defined input parameter: with σ being the Gaussian sigma from the normal distribution and t the mean water-equivalent path length. The larger the value of P mod , the broader the BP and the wider its distal fall-off. Both parameters ρ range and P mod are discussed in more detail below and are investigated within this work.
The water-equivalent voxel length t ρ H2O · d is then calculated for each individual voxel, where d is the voxel size and ρ H 2 O is the water-equivalent density of the given voxel. For the case that the beam is parallel to one of the three CT axes, d is simply assigned as either d x , d y , or d z , being the voxel length in the x-, y-and z-directions, respectively. If the beam direction is oblique, then the intersection lengths through each voxel can differ much from each other. We found that for our model, d can be well described by: where n (n x , n y , n z ) is the normalized direction vector parallel to the beam. This formula approximates the mean-quadratic length divided by the mean length of the intersection pieces, (averaged over many shifted beams and voxels) and was tested with a simple ray-tracing program.
In each of the n runs the following procedure is performed for all selected lung voxels. The water-equivalent path length through the voxel is randomized corresponding to P mod . The first approach would be a sampling by a Gaussian distribution with a width σ (P mod · d · ρ H2O ) 1/2 around the water equivalent density ρ H2O of the individual voxel. However, such a Gaussian density distribution, for voxel sizes in the millimeter range, has a relevant contribution in the negative range (for details, see Ref. 6), which must be corrected by an adapted distribution with a singular weight at ρ H 2 O 0. This is in accordance to the above mentioned binary voxel model and the corresponding Poisson distribution P(n), which describes the discrete probability distribution for hitting n non-void voxels when traversing a row of the binary voxels. However, since a discrete Poisson distribution cannot be applied here, a continuous Poisson distribution for the water-equivalent thickness t′ was chosen instead: The delta-Dirac function δ(x) adds the aforesaid singular weight for t′ 0 ( Figure 1). The gamma function replaces the factorial function normally occurring in a non-continuous Poisson distribution for non-integers in order to generalize the formula into the continuous Poisson distribution. The distribution parameters, δ, λ and A 0 , depend on the nonmodulated values t and the modulation power P mod . The dependencies were calculated for 30 different values for P mod , from 0.1 to 0.8 mm, and for 40 different t values, from 0.1 to 2.0, and interpolation between the matrix elements in these precalculated matrices were used to obtain other values. The mathematical principle behind these calculations are given in the subsection below.
For the density sampling of the individual voxels with an initial density ρ H2O and (mean) size d the following density distribution is then applied: At the end HU(ρ′) for all selected voxels are found using the previously loaded lookup-table and a final randomized (henceforth called "modified") CT, usable in TRiP98 for dose computation, is obtained. By recalculating the dose n times on n different modified CTs and averaging these n dose cubes, one obtains a dose distribution that takes the beam-modulation effect of the lung into account. This can then be compared to the original plan dose distribution. In comparison to the previous work [6], where P(t′) was individually optimized as a discrete function for each density value and each P mod , the method presented here is much faster and can be easily reproduced with the provided tables [22]. By simultaneously calculating on 24 CPUs, for n 100, the whole recalculation procedure lasted 15-20 min when using the TRiP dose calculation algorithm "all-point" and 50-100 min when using the more complex "multiple scattering" algorithm [23].

Calculation of the Poisson Distribution Parameters
The mathematical principle for the calculation of the tables of the Poisson distribution P(t′) parameters δ(t, P mod ), λ(t, P mod ) and A 0 (t, P mod ) is that these parameters were adapted such that when folding P(t′) N times by itself, the resulting distribution best fits the Gaussian distribution for the total path length N · t and P mod : N is chosen as a value large enough (typically 10), that the contribution of the negative part of the Gaussian distribution (right term in Eq. 5) becomes negligible. The idea behind this concept is that P(t′) substitutes the Gaussian distribution, avoiding non-physical negative contributions for the density distribution in the single voxels. The tables for δ(t, P mod ), λ(t, P mod ) and A 0 (t, P mod ) and the scripts for their calculation can be downloaded from Mendeley [22].

Benchmark With an Analytical Convolution
Our PyTRiP routine was tested against our well-established analytical convolution method [6] (which has been confirmed experimentally [4,8]) for 130 and 230 MeV/u carbon ion beams. For this test we used n 60, P mod 0.500 mm and ρ range from 0.1 to 0.5 g/cm 3 to calculate the modified doses. A lung tissue cube with constant HU of −741 HU [24] and a length of 15 cm in the beam direction was placed within a water phantom, starting at 3 cm water phantom depth z. For each depth value z of the modulated curve the pristine depth-dose curve (in water equivalent scale) was individually convoluted by a Gaussian filter, with a σ-value given as (t · P mod ) 1/2 , where t is the water equivalent thickness of the lung tissue in the beam path up to the depth z.
The comparisons of the Bragg curves can bee seen in Figure 2. The agreement between the analytical curves and the curves from our PyTRiP routine is clearly seen.

Lung Voxel Selection Range
When selecting the VOI lung voxels to be randomized, instead of using all voxels in the lung, we implemented a voxel selection by a user-defined density range ρ range . By using all lung voxels, one would end up modulating a case-specific variable number of nonporous voxels such as massive vessels.
We adopted values of ρ range from 0.1 to 0.5 g/cm 3 using experimental data from an inflatable pig lung set-up [4,6,8]  It is assumed that voxels within this range are lung voxels as well as any edge voxel that consists half of air and half of tissue. Both of these can lead to a modulation effect. In addition to this preselected range, we evaluated the effect on the patient plans when varying ρ range .

Study of Spherical Planning Geometries
For a systematic analysis of the modulation effect as a function of the isocenter depth d ic (placed at the PTV center) and the target size, plans for simulated spherical geometries were calculated for carbon ions. The spheres, composed of water, were placed within a cubic CT phantom, containing only voxels with a constant single Hounsfield unit value of −741, representing the lung environment [24]. First, a sphere of a constant 17 mm radius was placed at variable d ic values of 40, 70, 100, 125, 150, 175, 200, 225, 250, and 260 mm. Second, spheres with variable radii, from 8 to 15 mm in steps of 1 mm, from 17 to 25 mm in steps of 2 mm, and of 28 mm, were placed at a constant d ic 100 mm.
For the sphere planning study, n 100 and P mod 0.400 mm were used and all lung voxels were taken into account for the modulation. All plans were single field plans with the beam parallel to the z-axis. A 2 × 2 mm 2 lateral grid, a 3 mm depth step size and a lateral contour extension of 0.7·FWHM iso were used. All plans where optimized for the physical dose with the "ap"-algorithm. At d ic 40 mm, a 4 mm bolus was required.

Patient Selection and Planning
Two lung cancer patient cases, henceforth called P1 and P2, with different PTV volumes V PTV and locations were selected. They were originally treated with photons at the university hospital in Marburg. For P1 the PTV is deep-seated and located close to the spine, whereas the PTV for P2 is centrally located. For both cases, the tumors are not close to the thorax wall, as some distance in the lung is required for the degradation effect. The PTVs were sufficiently small as to not drastically decrease the lung volume to be modulated. The isocenter slice of both patient geometries, P1 and P2, is given in Figure 3, showing for each patient case the PTV position and an exemplary dose distribution calculated on the non-modulated CTs.
Multiple plans for each of the two patients were created for protons as well as for carbon ions, focusing on the latter. For each nonmodulated case (planned on the nonmodulated CTs), planning parameters such as raster scan grid size, energy step size and optimization margins yielding the best PTV coverage with respect to dose conformity and target homogeneity were found by iterations. For the nonmodulated plans the planning objective was to deliver at least 95% of the prescribed dose to 98% of the PTV [25]. All plans were 1-or 2-field plans and optimized for physical (absorbed) dose. For the carbon ion cases, a 3 mm thick ripple filter was applied to broaden the BP, as is routine in the clinic [26].
Coplanar fields were used for the study. In order to have data with beams traversing different distances through the lung tissue s lung , variable field angles were chosen. For the carbon ions, a horizontal field direction at 0°as well as −45°, +45°and +90°to horizontal were selected. For carbon ions dual-field plans with 0°a nd +45°field angles and with 0°and +90°field angles were  calculated as well. For protons 0°, −45°, and +45°to horizontal was used for P1 and 0°for P2. The 0°and +45°options match the fixed beamlines at the Marburger Ionenstrahl-Therapiezentrum and Shanghai Proton and Heavy Ion Center. Some of these field directions might not be clinically relevant, since sensitive structures and long distances through the lung are typically avoided, but were included here for investigative purposes. Table 1 shows the s lung for the investigated field angles and the V PTV values. The s lung values were estimated as the distance starting from the surface of the lung to the isocenter, placed in the middle of the PTV, along the central path of the beam. The V PTV values were given by TRiP.
For protons, the "multiple scattering" algorithm was used for optimization and dose calculation, which considers the broadening of the beam as a function of patient depth [15]. For the carbon ion plans the much faster "all-points" algorithm, which is considering only all neighboring raster beam spots that may contribute to a given voxel [15], was used as this is normally sufficient for carbon ion plans [27]. Although comparisons between carbon ion plans with the two algorithms showed a slight difference in the plans (data not shown), the relative coverage difference between a modulated and nonmodulated dose distribution when comparing the algorithms was, even in the worst case, less than 0.5 percentual point (pp).
In addition to the main study, we calculated for the patient case P1 a range of carbon ion plans with lateral contour extensions from 0.6·FWHM iso to 1.2·FWHM iso in steps of 0.1, and for each of these with beam spot spacing of 2 and 3 mm and depth steps of 2 and 3 mm. A coplanar horizontal field was used for all plans. Dose distributions, calculated on the original CT, were compared to modulated doses with P mod 0.400 mm and ρ range 0.1-0.5 g/cm 3 . The goal was to estimate the influence of the planning parameters on the relative difference in PTV coverage between the nonmodulated and modulated dose distributions.

Evaluation of the Parameters of the Lung Modulation Model
For the patient cases, the model parameters n, d, P mod and ρ range were evaluated.
The parameter n was evaluated for values of 24, 48, 72, 100, 120, 150, 200, and 400. For P mod , we chose 0.300, 0.400, 0.500, and 0.600 mm, as our former works show that realistic values of P mod lie below 0.600 mm. For the treatment planning cases with oblique fields for P1, plans with three different values of mean voxel length d were calculated: One following Eq. 2 yielding d 0.8439 mm, one using the unchanged length of the voxel d d x 0.9473 mm, and the last value was obtained using the maximum diagonal length through the voxel calculated with trigonometry as d 2 √ d x 1.339 mm. The parameter ρ range was evaluated by changing the maximum and minimum density value. We evaluated plans with ρ range of 0.2-0.4, 0.2-0.5, 0.1-0.3, 0.1-0.4, 0.1-0.5, 0.0-0.5, 0.0-0.8, and 0.1-0.8 g/cm 3 and additionally ran calculations modulating all lung voxels.
After the parameter study was done, selected values of these parameters were then used throughout the paper: n 100, ρ range 0.1-0.5 g/cm 3 , P mod 0.400 mm and d corresponding to Eq. 2.

Plan Evaluation
PTV coverage, dose homogeneity and planning conformity were used for dosimetric plan evaluations. In each patient case, the doses to lung, spinal cord, trachea and heart were compared between nonmodulated and modulated dose distributions. V 95% , being the volume receiving at least 95% of the prescribed dose, was used for evaluating the PTV coverage. The homogeneity index (HI) and the conformity index (CI) were defined as: where V 95%pi is the volume within the patient receiving 95% of the prescribed dose.

RESULTS
In general, the coverage and dose homogeneity of the PTVs were found to be worse for the modulated doses compared to the dose distributions calculated on the original CTs. This leads to a PTV dose overestimation by the TPS. The change in dose homogeneity was mainly due to a change in the minimum dose D 98% . As expected, the relative difference between the modulated and nonmodulated dose distributions increases for larger values of P mod . From n 72 to 400 the dose distribution for the modified plans only changes with maximum 0.1 pp in V 95% and D 98% , while D 2% remains constant. We adopted the value n 100 throughout the rest of the work. Figure 4 summarizes the results for modulated and nonmodulated dose distributions calculated for the spherical targets with 17 mm radius in a lung tissue CT cube. In Figures 4A,B), the DVHs for the nonmodulated and modulated dose distributions are shown for 6 different d ic values. In Figure 4C, V 95% as a function of d ic is plotted for both the nonmodulated and modulated dose distributions and the differences between the corresponding pairs of values of the PTV coverage are shown in Figure 4D. Figure 4D shows a clear linear dependency of the relative difference in coverage as a function of d ic .

Spherical Geometries
In Figure 5A, the V 95% values for plans with d ic 100 mm and variable PTV radii are shown. The relative difference in V 95% between the nonmodulated and modulated dose distributions is given in Figure 5B. In general the larger the PTV volume the better the coverage for both the nonmodulated and modulated dose distributions with an approximated exponential decrease of the relative coverage difference with PTV radius. In both cases of Figures 4C, 5A, some outliers can be seen in the individual V 95% -curves. Importantly, they occur for both the nonmodulated and modulated case, leaving their relative difference consistent. This is to be expected as the "modulated plans" are dose recomputations of the nonmodulated ones. Similar curve shapes were seen when using other coverage indices, as for example V 90% , V 80% or the mean dose (data not shown). One deviation from the smooth tendency that can be easily explained is the very superficial case at d ic 40 mm. Here, there is very little lung tissue in front of the PTV, resulting in a marginal difference between the nonmodulated and modulated dose distributions. Additionally, due to the low beam energies required at such low depth, a 4 mm bolus was applied, which might lead to an inferior coverage independently of the lung modulation effect.

Patient Plans
Dosimetric indexes V 95% , D 98% , D 2% and dose homogeneity HI of the carbon ion patient plans are presented in Table 2 (for P1) and Table 3 (for P2) for the nonmodulated cases and for the recalculated modulated dose distributions. The planning cases for P2 for ± 45°field angles have been omitted due to not fulfilling the planning objectives. In Table 2, dosimetric indices are additionally given for the three different values of  Frontiers in Physics | www.frontiersin.org November 2020 | Volume 8 | Article 568176 d used for the oblique beam cases. Dosimetric indexes for the proton patient plans are given in Table 4 for both patient cases. DVHs for PTV, trachea and lung in the nonmodulated and modulated dose distributions in selected one-beam plans for P1 (horizontal and oblique beams) and P2 (horizontal beam) are shown for carbon ions and protons in Figures 6, 7, respectively. No difference was found in the DVHs with and without lung modulation for the organs not displayed in these graphs. Zoomed plots for the DVHs of the PTV help to see the relative differences between nonmodulated and modulated dose distributions. Figures 6, 7, the relative differences between nonmodulated and modulated dose distributions are larger for carbon ions than for protons, which can be attributed to the generally sharper carbon ion BPs. Actually, for all cases of the modulated proton plans for P1, the drop in V 95% never compromised the planning objective of giving at least 95% prescribed dose to 98% of the PTV volume. Our conclusion is that the lung modulation effect is not very severe for proton treatment planning. This is in agreement with other findings [10,11]. We focus on the carbon ions plans for further analysis due to the larger modulation effect observed for this particle type relative to protons.   Of all the displayed DVHs, only those corresponding to the trachea for the case P2 exhibit a slightly changed maximum dose when considering the lung modulation effect (D 2% changes from 45.5% of the prescribed dose for the non-modulated plan to 38.9% for the case of P mod 0.600 mm). In the low dose region, the DVH for the trachea is slightly larger for the modified doses compared to the nonmodulated one. For example V 20% 7.3% of the prescribed dose for the nonmodulated dose compared to 10.4% of the prescribed dose for the dose distribution calculated with P mod 0.600 mm. None of these findings bear however a clinical relevance.

As can be seen comparing Tables 2-4 and by comparing
The dose to the lung displays minimal changes between modulated and nonmodulated dose distributions. Due to the generally smaller beam spot sizes of carbon ion beams, the ipsilateral lung could be better spared using carbon ions with approximately half the integral dose given to the lung compared to the proton plans. As an example, for P1, the V 5% with a horizontal proton beam is 18.8% of the prescription dose compared to 8.9% delivered with carbon ions.
Due to the tumor location of the selected patients, far from critical organs (the spinal cord for example), we will focus on assessment of the changes in target coverage as presented in Tables 2-4.
In Figure 8 is shown the relative difference between the V 95% of the nonmodulated doses and the modulated doses calculated for the values of P mod 0.300 mm (blue plots) and P mod 0.600 mm (green plots) as a function of the distance traversed through the lung s lung . The general trend of the sphere study, that a larger relative difference is seen for smaller PTV sizes, is present. However, that the relative difference is larger for a larger value of s lung is only partially confirmed. For P1, when comparing the single field plans at 0°and +45°, the relative differences between the nonmodulated and modulated dose distributions are larger for the +45°case. This is expected since at this angle the beam traverses a larger amount of lung tissue. However, the beam at 0°y ields a larger relative difference than for −45°by a large margin, even for the lowest P mod value, although s lung is slightly lower for the −45°setup. For the largest s lung value at +90°, a much smaller relative difference between nonmodulated and modulated doses are observed than for the 0°and −45°angles.
In the clinic, single field plans would rarely be used. Therefore, we also included selected few dual-field plans in the study. The dual-field plans are providing better and more conform target coverage, also for the nonmodulated doses and it can be seen that the drop in target coverage caused by the beam-modulation effect is slightly smaller for the dual-field plans relative to the nonmodulated ones.

The Effect of the Voxel Length d
As expected, and as can additionally be seen in Table 2, the relative difference between the plans with different values of d are significantly larger for the +45°cases where the particles traverse FIGURE 6 | DVHs of relevant VOIs for 12 C plans on the non-modified CTs compared to plans on n 100 plans summed for n 100 modified CTs with different P mod values. Three treatment planning cases are shown with P1 planned with a horizontal field as well as with one 45°to horizontal and with patient case P2 with a horizontal field. Plans were optimized and calculated for physical dose with a 2 × 2 mm 2 raster scan size and a 3 mm depth step size.
FIGURE 7 | DVHs of selected VOIs for proton plans on the non-modified CTs compared to plans on n 100 summed modified CTs with different P mod values. Three treatment planning cases are shown for P1 planned with a horizontal field as well as with one oblique field at +45°as well as for P2 with a horizontal field. Plans were optimized and calculated for physical dose with a 3 × 3 mm 2 raster scan size and a depth step size of 3 mm P1 and 2 mm for P2.
Frontiers in Physics | www.frontiersin.org November 2020 | Volume 8 | Article 568176 roughly twice the amount of lung tissue compared to the −45°c ases. Also, it can be seen that the larger the value of P mod , the more critical a variation in d becomes. However, even for +45°the largest decrease in V 95% seen between d 0.8439 mm and d 1.339 mm is 2.9 pp. For −45°, the largest decrease is only 0.4 pp.

Selecting Modulation Lung Voxels
Dosimetric indexes for plans with P mod 0.400 mm and variable ρ range values are listed in Table 5 for carbon ion plans and Table 6 for proton plans. Values for the nonmodulated cases are given for comparison. The same trends for V 95% and HI as for ρ range 0.1-0.5 g/cm 3 are seen. For protons, the dependency of the ρ range parameter on the plans is generally much lower when compared to carbon ions, as the relative differences between modulated and nonmodulated dose distributions are generally lower for protons. In the case of the largest variation for the protons plans, there is a decrease in V 95% of 1.4 pp when using the narrowest ρ range of 0.2-0.4 g/cm 3 compared to using all the voxels in the ipsilateral lung. There HI has an increase of approximately 43%. In comparison, for the carbon ions, the extreme case of P1 shows a decrease in V 95% of 4.2 pp when comparing ρ range 0.2-0.4 g/cm 3 with using all the voxels and an 58% HI increase. The less severe case of P2 shows for the same comparison a coverage drop of only 1.8 pp and a small HI increase of 16%.
Both of the extreme cases, ρ range 0.2-0.4 g/cm 3 or no fixed ρ range , are unrealistic choices. For P2, the V 95% values of the experimentally estimated ρ range 0.1-0.5 g/cm 3 lie in between the values of these two extremes for both protons and carbon ions. For protons, this is also the case for P1. For carbon ions for P1, the difference in V 95% from ρ range 0.2-0.4 g/cm 3 to ρ range 0.1-0.5 g/cm 3 is 1.7 pp, while for ρ range 0.1-0.5 g/cm 3 compared to the case with no ρ range limit is 2.5 pp. Figure 9 shows a boxplot with outliers for coplanar horizontal field plans for P1 with variable planning parameters (see the materials and methods section). Nonmodulated and modulated doses, calculated with P mod 0.400 mm, are shown separately. All plans that in the nonmodulated cases had V 95% below 98.0% or a All plans are calculated with horizonzal fields, e.g., with the +0°option with planning parameters given in the material and methods section. All plans are calculated with horizonzal fields, e.g. with the +0°option with planning parameters given in the material and methods section.

FIGURE 8 |
The difference in target coverage between the nonmodulated doses and the doses calculated with P mod 0.300 mm (blue plots) and P mod 0.600 mm (green plots) as a function of the amount of lung traversed by the particle beam (see Table 1). Circles mark data for P1, squares mark data for P2. For the P mod 0.600 mm plots the matching field angles are shown. All plans are for carbon ions.
Frontiers in Physics | www.frontiersin.org November 2020 | Volume 8 | Article 568176 CI above 2.5 were omitted from this study, as they did not fulfill the planning objectives. It is seen how V 95% , for the modulated dose distributions, varies more from a change in planning parameters. We observed that the better the nonmodulated plan was, in terms of PTV coverage and dose homogeneity, the smaller the relative difference between the nonmodulated and the modulated dose distribution.
For a small number of plans with specific parameter settings, the dose distribution was even improved by the lung modulation, in terms of PTV coverage and homogeneity. It was seen that, in some cases, it is possible to first calculate an inferior plan for the nonmodulated case, which due to the broader and slightly displaced BPs caused by the lung-modulation effect can result in a plan fulfilling the planning objectives for the modulated case.
As an example, we show in Figure 10 DVHs for a plan calculated on patient case P2, for a non-modified case (solid lines) compared to the sum of doses calculated on 100 modified CTs (dashed lines). A horizontal field was used, with beam spot spacing of 3 mm, energy step 3 mm, and lateral contour extension 1.2·FWHM iso with values of P mod 0.400 mm and ρ range 0.1-0.5 g/cm 3 for the calculation of the modulation effect. V 95% 97.7% and HI 0.102 were found for the nonmodulated doses and V 95% 97.0% and HI 0.061 were found for modulated doses, respectively. Although the values of the PTV coverage were similar for these two cases and slightly lower still for the modulated case, the HI value was reduced by one third after the modulation, as the much steeper fall-off for that DVH curve in Figure 10 also illustrates.

DISCUSSION
The beam-modulation effect caused by the porous material of the lung results in a lower PTV coverage and dose homogeneity. This effect is clearly more pronounced for carbon ions than for protons, as the sharper the initial BP the stronger the modulation effect. The modulation effect is generally negligible in proton treatments, as compared to the other dose delivery and treatment uncertainties related to protons, which confirms the work of Baumann et al. [10]. Among such uncertainties are particle range uncertainties, patient setup errors, interfractional anatomy changes (e.g., weight loss, tumor shrinkage, etc.) and intra-fractional tumor motion. Additionally, TPS algorithms often have problems handling the air cavities and sharp density gradients present in the lung.
Using spherical geometries, we show that when 1) increasing the distance the beam travels through the lung tissue and/or when 2) reducing the PTV size, the effect of the beammodulation increases. Both tendencies have been shown previously only for protons [9]. The latter effect fits the general knowledge that a good PTV coverage is more difficult to obtain for smaller PTVs when using broader BPs. For the non-modulated plans, the PTV coverage increases with an increasing depth, while for the modulated cases, the larger the depth of the target in the lung the larger the beammodulation effect. This altogether yields the increasing relative difference with depth between the non-modulated and modulated dose distributions.
These trends were only partially confirmed for the patient cases. The systematic study on spherical target geometries represents a simplified scenario, while for patient cases there are many degrees of freedom and many further parameters in FIGURE 9 | A standard boxplot (median, 1. and 3. quartiles and whiskers including the outliers) for V 95% , HI and CI data sets obtained for patient case P1 using a horizontal beam for non-modulated as well as modulated dose distributions for carbon ions. A set of variable planning parameters was used as given in the materials and method section. consideration, such as case-dependent planning parameters and much more complex patient geometries. The angle from which the beam comes might lead to a large change in the relative difference between modulated and nonmodulated dose distributions. Our conclusions are based on our limited number of planning cases, but in this work, a wider general test of the lung modulation model as well as the systematic study with spheres were prioritized. The largest uncertainties in our model (and as an extension in any similar routine dealing with the lung modulation) lie in the parameters P mod and ρ range . As P mod has been investigated in other works, we focused on ρ range , which has not been studied before. A too large ρ range would mean that non-porous biological lung structures, such as massive vessels, are falsely taken into account as porous voxels and the calculated modulation effect would be slightly larger than the actual modulation. On the other hand, by choosing a too small ρ range the modulation effect of the lung would be underestimated, which in particular for heavier ions can be important (as the effect is more pronounced for heavier ions as indicated in this work). We found that the upper limit might be more critical than the lower limit, since a too low upper limit would filter out too much porous lung material from the calculation. As long as the lower limit is kept within a reasonable value, which can easily be determined from CT voxel histograms, this value only marginally influences the modulation results. One should do at least a manual filtering of voxels containing visible big vessels and/or similar biological structures, which should not be modified.
Another, less important, parameter in the model is the voxel length d. For beam angles running parallel to the CT axis d is trivially the voxel length in that direction whereas, for any other angles, it must be calculated. The effect on the lung modulation of the value of d was found to be not so dramatic and only important when the distance through the lung traveled by the particles is large. We recommend using the tested formula, given in this work as Eq. 2.
The specifically chosen example shown in Figure 10 might indicate that, in some rare planning cases, an inferior dose coverage of the PTV can actually be improved by the lung modulation effect, although, in the large majority of cases, not taking the lung modulation effect into consideration in treatment planning is a net disadvantage for the patients.
The mathematical model described in this work could, in theory, be built into other TPSs too. We note, although, that the "future" of implementing the lung modulation into a TPS lies within MC codes when these will be coupled to treatment planning in a more general way. MC codes allow for on-the-fly modulation of the voxels during the actual treatment planning process and not by recomputing the dose on many different modulated CTs, followed by summing up the dose distribution to estimate the modulation effect. However, even with MC available, for the case that a deterministic TPS mode would still be opted for-for instance to calculate biologically effective dose with the Linear Effect Model (LEM) [14,28]-the concept presented in this work is an appropriate option for recalculation of treatment plans and for assessment of the lung modulation effect. For carbon ion therapy, the authors estimate that TPS with deterministic pencil algorithms are assumed to remain still in clinical use in the next years. Additionally, deterministic TPS programs, like TRiP98 and matRad [11], will play an important role for educational purposes.

CONCLUSION
With existing standard tools, it is possible in a transparent and reproducible way, to implement and calculate the beammodulation effect of the lung in a known TPS for particles. Our evaluation of two selected lung cancer patients is in line with previous results. For protons, the lung modulation effect is found to be negligible compared to other dose delivery uncertainties related to particle therapy treatments of the lung. For carbon ions, the effect is significantly more pronounced and cannot be ignored. The decrease in both PTV coverage and dose homogeneneity, caused by the modulation effect, could lead to underdosage within the PTV and to overdosage in the healthy normal tissue surrounding the PTV. Planning parameters and beam angles also have an impact on the relative difference between nonmodulated dose distributions and those taking the beammodulation effect into account, caused by the porous nature of lung tissue.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent from the participants was not required to participate in this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
UW devised the project, the main conceptual ideas and project outline. TR, AS, KZ, and UW further developed the project with AS and KZ supervising the work and providing financial support. TR, LG, and NB developed the software tools for the manuscript, based on a mathematical model developed by UW. TR performed the computations, the data gathering and the data analysis. TR and AS created the images and graphs. KB, VF, and RE provided the patient data. TR wrote the manuscript with input from all authors.