KAHIBARO
35.6. Proton Dose Analysis
Table of Contents
Depth-dose curve
In a proton therapy simulation, the main quantity of interest is how dose varies with depth in the irradiated medium, usually a water phantom. GATE provides this information through dose or energy deposition actors defined on a voxel grid that extends along the beam axis. For depth-dose analysis, you typically design the grid so that one axis, commonly the $z$ axis, follows the beam direction, while the other two axes span a transverse area that covers the beam.
After running the simulation, the actor outputs a 3D dose distribution, for example in an image format such as MHD or NIfTI, or in a numerical array. To create a depth-dose curve, you must reduce this 3D dose map to a 1D function of depth. The usual approach is to average the dose over all voxels at the same depth along the beam direction. If the beam travels along $z$, and $D(x,y,z)$ is the dose in a voxel, the depth-dose profile $D(z)$ is computed by taking the mean over the transverse plane for each $z$ slice. In Python, this is often implemented by loading the dose image into a NumPy array, identifying the axis corresponding to depth, and computing the mean across the other two axes.
It is important to ensure consistent units before interpreting the curve. Dose values from GATE are in Gray, defined as joules per kilogram, and the voxel size is in length units such as millimeters. When you plot $D(z)$ versus depth, use depth in millimeters or centimeters along the horizontal axis, and dose in Gray or normalized units on the vertical axis. For many analyses, the absolute dose magnitude is less important than the shape of the curve, so you may normalize the entire curve to its maximum value or to a reference point to focus on relative dose fall-off and Bragg peak characteristics.
The typical depth-dose curve for a monoenergetic proton beam has a low entrance dose region, followed by a gradual rise in dose, and finally a sharp peak corresponding to the Bragg peak near the end of the proton range. Beyond this peak, the dose drops quickly to a small tail created mostly by secondary particles from nuclear interactions. When you analyze a simulated curve, you should visually confirm that this qualitative behavior is present, since it is a basic indication that the proton transport and energy loss physics are configured correctly.
For comparison with analytical expectations or measurements, you can resample the simulated depth-dose curve on the same depth grid as experimental data or range data from published tables. Normalizing both curves to the same maximum facilitates shape comparison. You can then examine how quickly dose rises to the peak, how narrow the peak is, and how steeply dose falls off after the peak. If you simulated multiple beam energies or different beam sizes, plotting several depth-dose curves on the same graph helps you see how the Bragg peak moves and broadens with energy and beam parameters.
To reduce statistical noise in the curve, it is common to simulate a sufficiently large number of primary protons and, if necessary, increase the longitudinal voxel size slightly so that each depth bin contains many events. Thicker voxels yield smoother curves but lower spatial resolution. There is a compromise between resolution and statistical uncertainty. You can quantify this uncertainty by computing the standard deviation of the dose within each depth bin over repeated simulations or from uncertainty information provided by the dose actor.
For proton dose analysis, always construct the depth-dose curve by averaging dose over the transverse plane at each depth, use consistent units for depth and dose, and normalize the curve when comparing shapes between energies or configurations.
Range
The proton range in water is a central quantity in proton therapy, since it determines how deep the Bragg peak will be located in the patient. In GATE simulations, the range is determined from the depth-dose curve by identifying how far along the beam axis significant dose is delivered. While the physical definition of range is based on proton stopping distances, in practice you typically use dose-based range metrics that can be derived from the depth-dose profile.
The most widely used range definition in clinical practice is the depth at which the dose in the distal fall-off region drops to a given percentage of the peak dose. For a monoenergetic proton beam, the parameter $R_{90}$ is commonly used, defined as the depth where the dose equals 90 percent of the maximum dose on the distal side of the Bragg peak. Other common definitions include $R_{80}$ and $R_{50}$, which use 80 percent or 50 percent of the peak value.
To compute these quantities from a simulated depth-dose curve, start by identifying the maximum dose value, $D_{\text{max}}$, and its corresponding depth, $z_{\text{max}}$. Then focus on depths beyond the peak, where $z > z_{\text{max}}$. On this distal side, find the depth where the dose first falls below 0.9 $D_{\text{max}}$. If your dose data are sampled at discrete depths $z_i$, you can locate the two consecutive sample points $z_i$ and $z_{i+1}$ between which the dose crosses 0.9 $D_{\text{max}}$, and then perform a simple linear interpolation to estimate $R_{90}$ with sub-voxel precision.
In formula form, if $D(z_i) > f D_{\text{max}}$ and $D(z_{i+1}) < f D_{\text{max}}$, where $f$ is the fraction such as 0.9, the interpolated range $R_f$ is:
$$
R_f = z_i + (z_{i+1} - z_i)\,\frac{D(z_i) - f D_{\text{max}}}{D(z_i) - D(z_{i+1})}.
$$
The clinical proton range in water is usually defined as the depth at which the dose on the distal side of the Bragg peak falls to a fixed percentage of the peak, such as $R_{90}$ at 90 percent of $D_{\text{max}}$, computed by interpolation on the depth-dose curve.
Once you have computed range values such as $R_{90}$ from your simulation, you can compare them with reference data. For a given proton energy, reference ranges in water are tabulated in reports and in databases, and many treatment planning systems provide expected ranges. To validate your simulation, compute the difference between the simulated and reference ranges and express it in millimeters or as a percentage. Acceptable agreement depends on your application, but for basic validation you would typically aim for differences of a few millimeters or better in water.
When you change the beam energy in your simulation, the range should change predictably. If you simulate multiple monoenergetic beams with increasing energies, the Bragg peak should move deeper into the phantom, and the $R_{90}$ values should increase approximately monotonically with energy. Plotting $R_{90}$ against beam energy allows you to check whether your simulated range-energy relationship follows the expected trend. Significant deviations might point to issues in material definition, density, or physics configuration.
In some analyses, you may also be interested in other range-related quantities, such as the fall-off gradient or the width of the distal fall-off, for example between $R_{80}$ and $R_{20}$. While these are not strict definitions of range, they provide insight into how rapidly the dose decreases beyond the Bragg peak, which is important for sparing tissues located just beyond the target.
Peak position
The position of the Bragg peak is closely related to the proton range but provides more detailed information about where the maximum dose occurs in your simulation. To determine the Bragg peak position, you examine the depth-dose curve and locate the depth at which the dose reaches its maximum value. With sampled data, the simplest estimate is to find the index $i$ such that $D(z_i)$ is maximal, and then take $z_i$ as the peak position.
Because the depth-dose curve, especially in Monte Carlo simulations, may exhibit statistical fluctuations, you may want to refine the peak position estimate through interpolation or smoothing. One approach is to select a small region around the apparent maximum, for example a few depth points around the index of maximum dose, and fit a simple smooth function such as a parabola to these points. The vertex of the fitted parabola then provides a more accurate and less noisy estimate of the peak depth.
If you are analyzing multiple curves, for example for different proton energies or for different beam spreads, it is helpful to extract the Bragg peak position for each condition and compare them systematically. You can summarize the results in a table, where each row corresponds to a beam energy and columns list the Bragg peak depth and associated range measures. This provides a compact view of how beam parameters influence where the maximum dose is delivered.
| Beam energy (MeV) | Bragg peak depth (cm) | $R_{90}$ (cm) |
|---|---|---|
| 100 | 7.1 | 7.4 |
| 150 | 15.5 | 15.9 |
| 200 | 26.0 | 26.5 |
In addition to finding the peak position, you may also study the peak height and its shape. For example, you can normalize the curve to the maximum and then measure the full width at a specified percentage of the peak, analogous to a full width at half maximum. This characterizes how sharp the Bragg peak is and helps you understand how beam energy spread, energy straggling, and multiple scattering affect the high dose region. If you simulate a spread-out Bragg peak created by combining several energies, the depth-dose curve will show a broader plateau. In that case, the concept of a single peak position is replaced by a region of nearly constant dose, and you may define a central position or the position of maximum dose within the plateau.
When comparing simulated Bragg peak positions to measurements, such as depth-dose measurements in a water tank, you should account for the coordinate system used in your simulation. Ensure that $z = 0$ in the simulation corresponds to the same reference depth as in the measurement, usually at the phantom surface. If there is any shift in this reference, you may need to correct the simulated depths before comparing peak positions. Small systematic offsets can easily arise from how you positioned the beam relative to the phantom in your geometry.
Finally, it is often useful to report both Bragg peak position and a range metric such as $R_{90}$ for validation and analysis. The peak position focuses on where the highest dose occurs, while $R_{90}$ emphasizes where dose falls off to a specific relative level. Together, they provide a more complete description of the longitudinal dose distribution from your proton beam in the simulated phantom.
Views: 11