algorithms

In this section, the different steps of the computation carried out by the ConformalFLASH algorithm are described.

ConformalFLASH requires inserting several patient specific accessories in the FLASH snout:

• ConformalFLASH energy modulator (CEM) which is a patient specific, non uniform, ridge filter. The CEM recreates a spread out Bragg peak (SOBP) from a single energy layer. The principles are described in this publication. In addition, the CEM acts as a range compensator.
• Range shifter is composed of several aluminum slabs. The number of slabs will be adapted to the maximum range required to reach the distal surface of the PTV. There are several thicknesses of slabs: 16mm, 12mm, 8mm, 4mm. The range shifter will therefore have a total thickness that is a multiple of 16mm, possibly with an additional thinner slab.
• Aperture block to improve the lateral penumbra. The range shifter causes lateral scattering which will result in a large pencil bealm scanning (PBS) spots at isocentre. In order to improve the lateral penumbra, an aperture block is added.

Figure 1: The FLASH snout. The proton beams travels from the left to the right of the figure. The ConformalFLASH energy modulator (CEM) is attached to the green plate. The range shifter is composed of several slabs (orange) inserted in the snout. Empty range shifter slabs are blue. The brass aperture block is at the right of the figure.

Figure 2: Cross section through the high resolution CT scan with all the patient specific accessories virtually inserted in the CT. The dose distribution is overlaid as a color map. The base of the CEM (blue rectangle) acts as a residual range shifter. The shape of the base (red line) of the CEM acts as a range modulator that defines the distal surface of the dose distribution. The height of the spikes (green line) defines the proximal surface of the dose distribution. The shape of the spikes controls the flatness of the dose distribution in the SOBP. The range shifter is made of several slabs represented by the gray rectangles. The border of the aperture is represented by the white bars near the water phantom. The water phantom is represented in gray on the right of the figure.

The optimization of the treatment plan requires the optimization of several correlated degrees of freedom:

• The scanning trajectory of the PBS spots, i.e. the order in which the spots are delivered
• The weight of each PBS spots
• The shape of the CEM

Ideally, one would want to compute the dose map using a Monte Carlo dose engine at each iteration of the optimizer. However, due to the large number of degrees of freedom and the computation speed of a Monte Carlo dose engine, this would lead to unreasonable computation time. Therefore, several simplifications are required to reduce the computation time to a reasonable value:

Firstly, the scanning trajectory is not optimized at each optimizer iteration. The scanning trajectory is defined once, before the start of the optimization and it is no longer changed during the optimization. This section describes the trajectory optimization algorithm.

Secondly, the spot weight optimization and the CEM shape optimization are carried out in two sequential steps. In a first step, the spot weights of a reference IMPT plan (the 'mother plan') with multiple energy layers, a range shifter, and an aperture are optimized. This assumes that the CEM will be able to perfectly reproduce the dose map of this mother IMPT plan. The mother IMPT plan weight optimization is described in the "Optimize spot weights" section. The standard optimization algorithm of MIROPT is used but the objective function has been modified to include dose rate objectives which are described in the "Dose rate in objective function" section. At the end of the optimization, it is assumed that the dose map of the mother IMPT plan is the best achievable dose map.

The mother IMPT plan optimization assumes that the protons have an energy ranging from the maximum energy available in the treatment machine to a lower energy equal to the range difference between the distal and the proximal surface of the PTV. In order to bring the distal Bragg peak to the distal surface of the PTV, a range shifter is virtually inserted in the CT scan. In this way, the mother IMPT plan will be optimized using PBS spots with a spot lateral profile described by a Gaussian with a standard deviation $\sigma$ that is representative of the spot $\sigma$ obtained with the FLASH snout. Because of the lateral scattering in the range shifter, the spot $\sigma$ at isocentre will be large. In order to keep a reasonable lateral penumbra, an aperture block is also included (the aperture is included in the virtual CT scan as described in section "Insert aperture in CT scan" ). Furhtermore, to help the optimizer deal with steep water equivalent thickness (WET) gradients, a range compensator may be optionally added in the CT scan to conform the dose map to the distal surface of the PTV as described in this section.

Finally, the CEM shape is optimized in order to reproduce the dose map of the mother IMPT plan when using a single energy layer plan. The different methods used for this shape optimization are described in the section "Optimize CEM shape".

In the following paragraph, the sequential steps of the ConformalFLASH algorithm are described.

Compute range shifter thickness

For each beam, the water equivalent thickness (WET) from the skin surface to each voxel in the PTV is computed using the REGGUI function WEPL_computation. The HU to stopping power conversion uses the scanner parameters defined in scannerDirectory directory of the MCsquare library. The maximum WET in the PTV defines the maximum range $R_{max}$ required for the proton beam.

Using the MachineType parameter (read from the BDL), the ConformalFLASH library determines the maximum proton beam energy available in the treatment room and uses the energy2range function of REGGUI to determine the maximum range of the beam $R_{beam}$ that the machine can deliver.

The WET of the range shifter $R_{WET}$ is the difference between the two ranges: $R_{WET}=R_{beam}-R_{max}$. From the material properties defined in the MCsquare library, it is possible to convert the range shifter WET into a physical thickness of the range shifter.

The range shifter is composed of an integral number of aluminum slabs 16mm thick. Once the $R_{WET}$ is known, it is possible to divide the range shifting between an integral number of range shifter slabs. The residual range shifting will be provided by one additional thin range shifter slab and additional material being printed at the base of the CEM, as can be seen in Figure 22.

It is possible that the residual range shift would require a CEM base that is thinner than the minimum thickness that is printable (defined by the parameter MinThickness). If this is the case, then the last 16mm range shifter slab is replaced by a thinner range shifter slab and the CEM base thickness is increased accordingly.

Range compensator for mother IMPT plan

As the spot $\sigma$ is large, the optimizer cannot compensate steep water equivalent thickness (WET) gradients on the distal surface of the PTV by asjusting the weight of adjacent spots in successive energy layers. In order to provide good conformality of the dose map to the distal PTV surface, a range compensator is added in the CT scan used to compute the mother IMPT plan. The range compensator material is the same as the CEM. The maximum WET ( $WET_{max}$ ) is the water equivalent thickness between the skin and the most distal point of the PTV. The WET of the range compensator at position (x,y) ( $RC_{WET}(x,y)$ ) in the IEC gantry coordinate system is equal to:

$$RC_{WET}(x,y) = WET_{max} - WET(x,y)$$

where $WET(x,y)$ is the WET from skin to the distal surface of the PTV at position (x,y) in IEC gantry coordinate system.

To compute the mother IMPT plan, the range compensator is inserted in the CT scan between the aperture and the range shifter (see Figure 2b ). When computing the mono-layer plan with CEM, the range compensator will be added directly into the CEM.

Figure 2b: Cross section through the CT scan with the virtual patient specific accessory inserted to compute the mother IMPT plan. The range shifter is made of a single slab represented by the gray rectangles. The border of the aperture is represented by the white bars near the water phantom. The range compensator is inserted upstream to the aperture. The water phantom is represented in gray on the right of the figure.

Compute optimum spot Spacing

When the parameter SpotSpacing = 0, ConformalFLASH computes the optimum spot spacing. The optimum spot spacing to achieve uniform dose coverage at isocentre depends on the spot $\sigma$ at isocentre.

ConformalFLASH computes the dose distribution in the patient CT scan (with the range shifter) for one single PBS spot on the optical axis and with the maximum energy of the machine. The intermediate results of the computation are saved in the folder /Beamlets/SingleSpot.

The depth of the Bragg peak maximum is identified in the IDD (integrated depth dose) curve computed by integrating, at each depth z, 3D dose distribution over the plane orthogonal to the proton beam axis. A 2D Gaussian is fitted to the 2D-dose profile in the orthogonal slice at the depth of the Bragg peak maximum. The result of the fit gives and estimate of the spot $\sigma$ at isocentre.

Let’s consider a one-dimensional Gaussian. We want to have the largest possible distance between the Gaussian spots, so as to have the largest possible field. In addition, we want a dose as uniform as possible across the field. Let’s define the maximum dose deviation $\Delta D$ (in %) within this area:

$$\Delta D = 200.\frac{max(D(x))-min(D(x))}{max(D(x))+min(D(x))}$$

where D(x) is the dose at position x.

The optimum distance giving a minimum ΔD and the largest field size is determined for different values of spot $\sigma$. Note that we are interested in the spot spacing at isocenter where the spot is mostly circular due to the scattering in the range shifter and the patient. Therefore, we do not need to consider different $\sigma$ in X and Y directions. The simulation results are shown in Figure 3. A linear regression is fitted to the simulated points in order to obtain a linear equation relating the optimum spot spacing with the value of the spot $\sigma$.

Figure 3 : Optimum distance between PBS spots as a function of the spot $\sigma$. The blue circles are the simulation results. The dotted red line is the linear regression.

Insert aperture in CT scan

The contour of the aperture opening is defined by the parallel projection of the PTV in the isocenter plane. The PTV is first dilated in 3D by the value defined in the yaml variable TargetMargin before the projection. The coordinates of the apexes of the polygon defining the aperture border are saved in the RT plan in the tag (300A,03A2) IonBlockSequence. The coordinates correspond to the polygon projected in the isocentre plane.

Then, the polygonal aperture is projected into the aperture plane taking into account the anamorphic magnification due to the fact that the source to axis distance for the X-scanning magnet and Y-scanning magnets are different.

The aperture block is placed as close as possible to the body contour RT struct while making sure that the brass block does not intersect the body contour. If the parameter AirGap is defined, the defined air gap is added between the skin and the aperture. The distance of the aperture plane from isocentre is saved in the treatment plan in the DICOM tag (300A,00F7) IsocentreToBlockTrayDistance. The aperture block is added into the planning CT scan (using the spatial resolution of the planning CT scan) by setting the HU of the relevant voxels to the HU of brass. If the original CT scan is to small to fit the aperture, the CT is expanded before adding the aperture.

The resulting CT scan in saved in the output folder 'output_path'/Outputs/ct_IMPT/.

Figure 4: Three orthogonal slices through a (virtual) CT scan of a water phantom (gray) and a zoomed axial view. The aperture block (white) is digitally inserted as close as possible to the body contour (blue line). The range compensator and the range shifter are digitally inserted in the CT scan by setting the HU of the voxels to the material of the CEM. For the mother IMPT plan, the range compensator is placed against the aperture while for the monolayer plan the range compensator is merged with the CEM. For the computation of the mother IMPT plan, the range shifter is one single block with the required WET (as there is no CEM to provide additional range shifting). This is the CT scan used to compute the mother IMPT plan. The target RT struct is shown in yellow. The dose map of the mother IMPT plan is overlaid on the CT scan as a color map.

Optimize mother IMPT with range shifter and range compensator

At this point, the aperture, range compensator and range shifter are inserted in the planning CT scan plan (see Figure 4). In this chapter, the spot weights of the mother IMPT plan with multiple energy layers are optimized.

Define spot trajectory

MIROPT place the PBS spots on a period lattice in the beam's eye view (BEV) plane. The lattice can be square or hexagonal (defined by the parameter GridLayout). The method to compute the lattice period (distance between the spots) is describe in the "Compute optimum spot Spacing" section. It is specified by the parameter SpotSpacing.

To define the zone in the BEV plane where to place the spots, the PTV is first dilated in 3D using a convolution by a kernel with a size specified by the parameter TargetMargin. Then the dilated PTV is (parallel) projected into the BEV plane. The external contour of the projection defines the zone where to place the spot lattice. The goal of expanding the PTV is to take into account the penumbra and the patient position errors. For example, in the Figure 5, you can notice that some of the spot center (crosses) are outside of the aperture contour (red line). This is due to this PTV dilation. This indicates that some PBS spots will be partially delivered into the aperture block and only a part of the spot will be delivered into the field. This give more freedom to the optimizer to obtain a flat dose profile near the field edge.

In addition, the ConformalFLASH library allows the lattice to be positioned with respect to the PTV in two different ways, depending on the value of the parameter SpotAtIso:

• The lattice is shifted so that one spot is placed at isocentre
• Or the lattice is shifted so that the first spot is placed at the edge of the (dilated) PTV projection

Figure 5: MIROPT places the spots on a periodic lattice in the beam's eye view plane. The lines define the iso-WET contour from the skin to the distal surface of PTV. The red line defines the aperture contour in the isocentre plane. The center of the PBS spots is indicated by the crosses. This is a square lattice positioned with spot at isocentre.

There is one constraint on the position of the PBS spots between the different energy layers for the mother IMPT plan. The spots must be aligned (along the proton beam axis) between the different layers. Indeed, at the next steps, the multi energy layer plan will be converted into a single energy layer plan and the CEM will recreate the Bragg peaks of the spots aligned along the proton beam axis.

The (water equivalent) distance between energy layer is defined by the parameter LayerSpacing.

Figure 6: The spots of the mother IMPT plan in the different energy layers must be aligned along the proton beam axis. The mother IMPT plan is computed with a range shifter and an aperture.

At this point, the spot position is defined and the scanning order of the spots in BEV is defined. When computing the dose rate for the mother IMPT plan, it is assumed that all the PBS spots (in different energy layer) at the same position in the orthogonal plane are delivered simultaneously, i.e. it is assumed that all energy layers are delivered simultaneously and that only (X,Y) scanning takes time. This simulates what effectively happens during the delivery of the single energy layer plan with the CEM.

Beam data library

To compute the PBS spot properties in the nozzle, MCsquare uses a beam data library (BDL) that must be commissioned for each treatment machine. The beam data library of MCsquare is defined in a text file file stored in the folder plugins/openMCsquare/lib/BDL. The BDL file for ConformalFLASH has some specific features.

Firstly, the spot shape must be described accurately at the base of the CEM. Therefore, the spot $\sigma$ and the beam divergence must be measured at the base of the FLASH snout for all snout extensions. As the FLASH snout can be retracted inside the body of the universal nozzle, the PBS spot must be characterized even when the base of the snout is fully retracted inside the nozzle. The parameter Nozzle exit to Isocenter distance in the BDL must be set to the value corresponding to the distance between the isocentre and the base of the snout with the snout holder fully retracted.

Secondly, the PBS spot is elliptical and the main axis of the spot is not aligned with the IEC gantry coordinate axes. MCsquare was modified to describe ellipses that are rotated with respect to the IEC gantry axes. The Figure 10 shows an example of a PBS spot at the base of the CEF.

Figure 10: Dose map in a plane transversal to the proton beam axis overlaid onto the slice of the virtual CT scan at the base of the CEM. One PBS spot is represented in the dose map. The spot is elliptical and the main axis is tilted with respect to the IEC gantry coordinate system.

The spot $\sigma$ along the two axes of the ellipses are measured at several distances from isocentre, as shown in the Figure 11. In addition, the rotation angle between the axes of the IEC gantry coordinate system and the ellipse axes is determined. MCsquare defines a spot as the sum of two 2D Gaussians as described in the calibration procedure of MCsquare. A new parameter is added in the BDL: SpotTilt which represents the rotation angle (in radian) between the main axis of the BDL and the axis of the IEC gantry coordinate system.

Figure 11: Spot $\sigma$ along the two ellipse axes as a function of the distance between the isocentre and the measurement plane. The circles are the measurements. The lines are the fitted Courant Snyder equations.

Secondly, when the proton energy is changed using the degrader and the ESS, the spot properties at isocentre are a function proton energy. In the case of ConformalFLASH, one single energy is used: the maximum energy of the treatment machine. The spot properties must therefore be characterised only once, for this single energy. However, as is described in the section "Optimize spot weights", ConformalFLASH computes the mother IMPT plan which requires several energies to be defined in the BDL. The spot properties of these 'virtual' energy layers will be the same as the spot properties of the maximum energy because they will actually be delivered using a CEM with the maximum energy. Therefore, in the BDL, the spot properties of the maximum energy layer are copied, identically, to all the other beam energies.

Thirdly, when computing the beamlets for the mother IMPT plan (see section "Optimize spot weights"), a range shifter is added in the CT scan (see Figure 2b). In ConformalFLASH, the range shifter is a patient specific accessory. Note that BDL contains the field RS_WET in the section "Range Shifter parameters" where the default thickness of a range shifter is defined. This field is mandatory in the BDL but it will not be used in our application. Indeed, the RS_WET of the BDL is used by default by MCsquare when no range shifter WET is explicitly defined in the plan (e.g. when the plan refers to a range shifter permanently fixed to the nozzle drawer). In our case, ConformalFLASH add the range shifter in the CT scan and removes it from the plan before calling MCsquare. The field RS_WET of the BDL will be ignored by MCsquare in ConformalFLASH.

Compute dose influence matrix

To optimizes the spot weights of the mother IMPT plan, MCsquare computes the 3D dose distribution of each beamlet (for each energy layer) with a weight $\omega_j=1$ where j is the index of the beamlet. The index j run for all (x,y) spot positions in the orthogonal plane and all energy layers E. To compute the spot dose distribution, the following geometry is used in MCsquare (see Figure 2b):

• a range shifter is added in the CT scan. The range shifter thickness is already known. The isocentre to range shifter distance can be computed as the isocentre to aperture block distance is known and the aperture thickness is also known. The scatter occurring inside the range shifter is computed by MCsquare so that the correct PBS spot size is obtained at isocentre.
• a range compensator is added to the CT scan, upstream to the aperture block.
• an aperture block is virtually added in the CT scan. The aperture block is modeled by voxels with the HU of brass inside the CT scan.
• The CT scan resolution is the original CT scan resolution; typically 1x 1 x 2mm.

The intermediate results are saved in the folder /Beamlets. MIROPT creates the dose influence matrix $P_{i,j}$ where i is the index of the voxel in the CT scan where the dose is delivered and j is the index of the beamlet delivering the dose. MIROPT then defines a weight vector $\omega_j$ where j is one component of the vector representing the weight of the j-th beamlet. The total dose $D_{i}$ delivered in the i-th voxel of the CT scan by all the beamlets is:

$$D_i = P_{i,j} * \omega_j$$

where * is the matrix multiplication.

Optimize spot weights

MIROPT varies the value of the elements of the $\omega_{j}$ vector until the minimum of the objective function $GOF(\omega_{j})$ is reached. The objective function is a sum of several terms:

$$GOF(\omega_j) = \sum_k GOF_k(\omega_j , S_k)$$

where $S_k$ is the structure in which $GOF_k$ is computed.

The terms $GOF_k$ of the objective function can represent objectives on dose, as in the standard MIROPT, but, in ConformalFLASH, they can also represent objectives on the dose rate:

• minDR : The percentile dose rate is computed for each voxel in the structure $S_k$. The average ( $DR$ ) of this dose rate is the computed by including all the voxels that receive a dose higher than parameter Dref. A penalty proportional to $(DR - DRref)^2$ is added to the objective function when $DR \lt DRref$.
• minDRm : The percentile dose rate is computed for each voxel in the structure $S_k$. The median ( $DRm$ ) of this dose rate is the computed by including all the voxels that receive a dose higher than parameter Dref. A penalty proportional to $(DRm - DRref)^2$ is added to the objective function when $DR \lt DRref$.
• minDADRm : The dose average dose rate is computed for each voxel in the structure $S_k$. The median ( $DADRm$ ) of this dose rate is the computed by including all the voxels that receive a dose higher than parameter Dref. A penalty proportional to $(DADRm - DRref)^2$ is added to the objective function when $DADRm \lt DRref$.

where DRref is a parameter defined in the YAML file.

MIROPT varies the weight $\omega_{j}$ of the PBS spots until a minimum of the objective function $GOF(\omega_{j})$ is reached.

At the end of the optimization of the mother IMPT plan, if the parameter ComputeIMPTdose = true, MIROPT saves the dose map of the mother IMPT plan in the file MCsquare_Dose_IMPT.dcm in the /Output folder.

In the following sections, a definition of the different dose rates is given.

Dose rate in objective function

The Figure 12 shows the dose delivery vs time simulated in one voxel of the structure. Every time a PBS spot overlaps with the voxel, there is a burst of dose deposited in the voxel. This is represented by the increase of the instantaneous dose rate observed at some specific times in the graph. When the PBS spot is moved far enough from the voxel, no further dose is deposited to that voxel. The blue curve in the Figure 12 is the cumulative dose, i.e. the integral of the red curve.

FIGURE 12: Red (right scale): Peak dose rate (Gy/s) as a function of time in one voxel of the CT scan. Blue (left scale): cumulated dose as a function of time at the same point. The blue curve is the integral of the red curve.

Percentile dose rate

In each voxel of the structure, a X-percentile dose rate can be defined. The X-percentile dose rate is defined as the ratio of the X-percentile dose divided by the time required to deliver that percentile of the dose. For example, one could define the 98-percentile dose rate or the 95-percentile dose rate. Or even the 100-percentile dose rate. The percentile level is defined by the parameter Vref in the YAML.

$$DR_{98} = \frac{D_{98}}{T_{98}}$$

Dose averaged dose rate

The dose average dose rate was defined in the paper by van de Water et al. [van de Water, S., Safai, S., Schippers, J. M., Weber, D. C. & Lomax, A. J. Towards FLASH proton therapy: the impact of treatment planning and machine characteristics on achievable dose rates. Acta Oncol. (Madr). 58, 1463–1469 (2019)].

Optimize CEM shape

At this point, the weight $\omega_j$ of the PBS spots (of all energy layers) has been optimized in order to minimize the $GOF(\omega_{j})$ objective function. The next step will be to replace the multiple energy layers by a single energy layer and a CEM shape that will give the same SOBP than the mother IMPT plan. The principle is shown in the Figure 14.

Figure 14: Principle of the CEM. Instead of delivering sequentially energy layers, one single energy layer (at maximum energy) is delivered. The protons of the PBS spot will go through different thickness of material in the CEM (represented by the colored blocks) and therefore undergo different range shifts. They will generate Bragg peaks at different depths in the target (represented by the different colored spots, matching the color of the blocks in the CEM). The cross section area of a block determines the number of protons at a defined energy (i.e. with a defined range in water) and therefore the weight of the Bragg peak at that range.

Firstly, let's rename the spot weight $\omega_j = \omega_{s,E}$ to indicate that the index j runs over the spots position s (= (x,y) in a plane orthogonal to the proton beam axis) and over the energy layers E. As the spots were aligned along the proton beam axis, there will be series of spots at the same s position but with different energies E. In the mother IMPT plan, the dose deposited in voxel (x,y,z) by the spread out Bragg peak (SOBP) generated by the spots at position s and with different energy E is:

$$D_s(x,y,z) = \sum_E \omega_{s,E} . BP_{s,E}(x,y,z)$$

where $BP_{s,E}(x,y,z)$ is the dose deposited at voxel (x,y,z) by the spot at position s with energy E. This is equal to the dose influence matrix $P_{i,j}$.

In the single energy layer plan, there is one single spot at position s with weight $\omega_s$ such that:

$$\omega_s = \sum_E \omega_{s,E} . \frac{\alpha_E}{\alpha_{max}}$$

During the optimization of the mother IMPT plan, each layer has a different energy and therefore, the ionization chamber (IC) will measure protons of different energies for each layer. The BDL shows that the calibration between MU (spot weights) and dose at isocentre is a function of proton energy (see Figure 15). However, in the case of the monolayer plan, the IC will measure protons of maximum energy only. Therefore, the weight $\omega_{s,E}$ must be rescaled to protons of maximum energy. In the equation, $\alpha_E$ is the number of protons per MU at energy of the layer E and $\alpha_{max}$ is the number of proton per MU at the maximum energy of the monolayer plan. Note that in the BDL described i nthe section above, the ratio $\frac{\alpha_E}{\alpha_{max}}$ is equal to 1 at all energies.

FIGURE 15: Proton charge (nCb) in a PBS spot for one MU vs the proton energy. The relation is computed using the function MU_to_NumProtons.m. The calibration curve is defined in the BDL of MCquare for each treatment machine.

First approximation to area $A_E$

To recreate an SOBP from the mono-energy layer plan, the CEM contains different height $h_E$ of material covering different fraction $A_E$ of the surface of the PBS spot. Let's define the range $R_E$ in water of the layer E and the range $R_{max}$ in water of the proton of maximum energy. The physical height $h_E$ of material to add to the CEM to shift the proton range from maximum energy down to energy E will be:

$$h_E = \frac{R_{max} - R_E}{SPR}$$

where $SPR$ is the relative stopping power of the material used to manufacture the CEM.

Let's suppose that the CEM is divided in multiple cells and let's assume that the cells are roughly of the size of a PBS spot. As there can be some overlap between PBS spots at the base of the CEM, the cells can be slightly smaller than a PBS spot. The cells form a lattice with the same geometry than the PBS spot lattice. The lattice period of the cells is the same as the lattice period of the PBS spots, projected in the plane of the base of the CEM. Inside the cell at position s, the material of height $h_E$ will cover an area $A_E$. In a first approximation, to recreate the weight $\omega_{s,E}$, the area $A_E$ shall respect the following equation:

$$\omega_{s,E} . \frac{\alpha_E}{\alpha_{max}} = \frac{\int \int_{A_E}F(x,y).dx.dy}{\int \int_{A_{base}}F(x,y).dx.dy}$$

where F(x,y) is the fluence (protons/mm2) at position (x,y) in a plane normal to the proton beam axis and $A_{base}$ is the area of the cell.

The parameter PreOptimization defines the function to use to describe the fluence F(x,y). If the PBS spots are sufficiently close at the base of the CEM, one can assume a uniform fluence, i.e. $F(x,y)=const$. In that case PreOptimization = uniform.

If the PBS spots are sufficiently far form each other at the base of the CEM, than the fluence can be described by a Gaussian function:

$$F(x,y) = \frac{A}{2 . \pi . \sigma_x . \sigma_y . \sqrt{1-r^2}} . e^{-\frac{(X/\sigma_x)^2 - 2.r.(X/\sigma_x).(Y/\sigma_y) + (Y/\sigma_y)^2}{2.(1-r^2)} }$$

In that case PreOptimization = Gaussian. $\sigma_x$ and $\sigma_y$ are the sigma of the lateral spread of the spot at the base of the CEM.

Second approximation to area $A_E$ including lateral scattering

In the first approximation, the lateral scattering occurring inside the material of the CEM was neglected. This is not an acceptable simplification because the protons may be scattered to a zone with a thicker or thinner $h_E$ and therefore their range will be changed: they will contribute to a different Bragg peak of the SOBP and this will change the relative weight of the Bragg peaks. This effect must be taken into account. For example, compare the fluence map for protons at 216MeV at the base and at the top of the CEM in the Figure 16. At the base of the CEM, the structure of the CEM is clearly visible in the fluence map. At the top of the CEM, a "blur" in the fluence map is visible. This represents the protons that were scattered out of the areas of different heights with an energy 216MeV. This will increase the weight of the 216MeV Bragg peak and lead to a less flat SOBP.

One option to model the lateral scattering is to use a Monte Carlo simulation. However, iteratively optimizing the CEM using a Monte Carlo simulation at each iteration would be time consuming. In order to reduce the computation time, an analytical simplification is done. The lateral scattering is modeled using the Moliere theory.

FIGURE 16: Color map of the proton fluence for different proton energy. Fluence at the base of the CEM. In this example, there is one single PBS spot delivered through the CEM. The Gaussian distribution of the fluence of the PBS spot is visible in the maps. If there were no scattering and the protons were moving in straight lines, this map would represent the proton fluence (for different energies) at the exit of the CEM.

Assuming that the protons travel in a straight line (i.e. without scattering), the fluence of proton leaving (with an energy E) a step of height $h_E$ the CEM is the sum of the fluence of all the PBS spots times a mask defining the CEM shape:

$$F^b_E(x,y) = M_{E}(x,y) . \sum_s F_{s,E}(x,y)$$

with

$$F_{s,E}(x,y) = \frac{\omega_s}{2 . \pi . \sigma_x . \sigma_y . \sqrt{1-r^2}} . e^{-\frac{((X-s_x)/\sigma_x)^2 - 2.r.((X-s_x)/\sigma_x).((Y-s_y)/\sigma_y) + ((Y-s_y)/\sigma_y)^2}{2.(1-r^2)} }$$

where the s-th PBS spot is centered at position $ (s_x,s_y)$ and has a weight $\omega_s$. The spot $\sigma_x$ and $\sigma_y$ are those at the base of the CEM. $M_E(x,y)$ is the binary mask defining the (x,y) position of the step with height $h_E$.

Let's now consider the lateral scattering inside the step $h_E$. For each fluence map $F^b_{E}(x,y)$, the increase of the spot $\sigma$ induced by the thickness $h_E$ of material can be predicted by the Moliere theory. The researchers can use user defined functions to define the lateral scattering taking place in the CEM. This is described in more details in the section User defined function for lateral scattering.

Downstream to the CEM, at the entrance of the range shifter, the fluence $F^{CEM}_{E}(x,y)$ of the proton with energy E at position (x,y) is:

$$F^{CEM}_E(x,y) = F^b_E(x,y) * G_{\sigma_M(h_E)}(x,y)$$

where (x,y) are expressed in the IEC gantry coordinate system and G() is a 2D-Gaussian function. The symbol '*' represents a convolution and:

$$\sigma_M(h_E) = \sqrt{ (\theta_M(h_E) . h_E)^2 + (\theta_M(h_E) . d )^2 + (\theta_B . h_E)^2 + (\theta_B . d )^2 }$$

with $\theta_M(h_E)$ is the characteristic multiple angle scattering angle at the top of the step of height $h_E$ as predicted by the Moliere theory, d is the distance between the top of the CEM and the base of the range shifter. $\theta_B$ is the intrinsic divergence of the PBS beamlet at the base of the CEM.

We now need to take into account the fact that some protons will be scattered outside of the step $h_E$:

• The proton scattered out of step $h_E$ into air are assumed to continue their journey with within the same energy E
• The proton scattered into a step with larger thickness $h_E'$ with $E' < E$ will be allocated to the fluence map $F^{CEM}_{E'}(x,y)$

The scattered protons are allocated to the fluence map with the correct energy by using the binary masks $M_E(x,y)$ defining the shape of the differents CEM steps in the plane normal to the proton beam axis.

FIGURE 17: This CEM is made of 3 paralellipipeds with different heights. The CEM is represented in two orthogonal projection planes. In the horizontal plane, the masks $M_E(x,y)$ define the (x,y) position of steps with height $h_E$ . In the vertical plane, the the lateral scattering of the proton through the step $h_E$ is represented. The protons remaining in the step and the proton scattered into air are assumed to exit the CEM with energy $E$ . The protons scattered into a taller step $h_E'$ (with $E' < E$ ) are assumed to exit the CEM with at the lower energy $E'$.

Taking into account the protons scattered into higher steps, the equation for the fluence map at the entrance of the range shifter becomes:

$$F^{CEM}_E(x,y) = F^b_E(x,y) * G_{\sigma_M(h_E)}(x,y) . \sum_{\epsilon=1}^E M_{\epsilon}(x,y) + \sum_{\epsilon=E+1}^{E_{max}} F^b_{E}(x,y) * G_{\sigma_M(h_\epsilon)}(x,y) . M_{\epsilon}(x,y)$$

Further lateral scattering will occur in the range shifter, and the Moliere theory is used to compute the standard deviation $\sigma_M(RS)$ of the Gaussian convolution kernel describing the scattering in the range shifter.

FIGURE 18: Color map of the proton fluence for different proton energies. Fluence at the exit of the range shifter. Because of the lateral scattering inside the CEM, thefluence map are 'blured'.

Further scattering will take place inside the water tank. This scattering is modeled in the same way using the Moliere theory and a convolution of the fluence maps with a Gaussian kernel. The resulting spot $\sigma$ is computed in a plane corresponding to an hypothetical water phantom located at the position of the 'body' RT-struct. This results in the 'blur' visible in the Figure 19. The fluence maps for the different energies in the water tank are shown in the Figure 19.

FIGURE 19: Color map of the proton fluence for different proton energies after going through the CEM, the range shifter and water. There is one single PBS spot. Fluence inside the water tank.

This Figure 19 represents the modeled fluence map $F^m_E(x,y)$ (in a plane orthogonal to the proton beam axis) for protons of energy E in a water tank when the the PBS spot at position s has been through the CEM, a range shifter and water.

Finally, the effect of the range straggling in the CEM, range shifter and water is modelled analytically by a convolution:

$$F^{strg}_E(x,y) = \sum_{\epsilon = E}^{E_{max}} F^m_E(x,y) . e^{\frac{-(R(E)-R(\epsilon))^2}{2.\sigma(\epsilon)^2} }$$

where $R(E)$ is the range in Lexan of proton of energy E. $\sigma(E)$ is a semi-analytical expression of the $\sigma$ of the Gaussian kernel for the range straggling convolution as a function of the proton energy. The method used to establish the semi-analytical model of range straggling is described in the next sub-section.

The fluence $F^{strg}_E(x,y)$ modelled for the geometry CEM + range shifter + water is compared to the reference fluence $F^r_E(x,y)$ that one expects for the PBS spots of energy E :

$$F^r_E(x,y) = (\sum_s \omega_{s,E} . F_{s,E}(x,y) ) * G_{\sigma_M(water)}(x,y)$$

where '*' is a convolution. $G_{\sigma_M(water)}(x,y)$ is a Gaussian convolution kernel representing the scattering in water. Moliere theory is used to compute the standard deviation $\sigma_M(water)$ of the convolution kernel describing the scattering in the water phantom.

The objective function is then defined as the difference between the reference fluence and the modeled fluence:

$$GOF_s(m) = \sum_{x,y} \sum_E (F^r_{E}(x,y) - F^m_{E}(x,y))^2$$

The objective function depends on the modeled fluence map $F^m_{E}(x,y)$ which depends on the shape of the CEM. Starting from the first approximation of the CEM, the shape of the CEM is iteratively modified until a minimum is found for $GOF_s(m)$. The iterative optimization is carried out for each cell separately. All cells are optimized sequentially in order to reduce the number of degrees of freedom for the optimizer.

FIGURE 20: Elevation maps of the CEM in one cell. This example is a fractal CEM. The different square spikes are visible in the elevation map. The range compensator is added to the CEM and appears as a smooth gradient accross the elevation map.

FIGURE 21: Difference between the modeled fluence $F^m_{E}(x,y)$ and the reference fluence $F^r_{E}(x,y)$ in one cell for the different proton energies. The maps are represented in a plane normal to the proton axis. In the last row, the integral of the fluence map is represented as a function of the energy layer. The blue line is the reference fluence and the red line is the fluence after CEM. The last plot is the histogram of the number of towers in the cell as a function of the tower height.

Every time the shape of the CEM is optimised for one beam, a temporary file Plan_ridge_tmpX.mat (where x is the beam number) is saved in the Output path folder. When the CEM of all beams have been computed, the temporary file planRidge.mat is saved in the same folder.

At this point, the shape of the CEM is optimized so that the proton fluence (for different proton energies) in a (virutal) water phantom downstream to the range shifter and the CEM is as close as possible to the proton fluence expected for the mother IMPT plan (with range shifter). The lateral scattering in the CEM, range shifter and water tank is taken into account analytically using the Moliere theory.

User defined function for lateral scattering

The default function predicting the spot $\sigma$ at the exit of a defined thickness of CEM is based on the Moliere theroy. However, the researchers can provide their own function to described this scattering process. The parameter Scattering is set to User, then the algorithm will call a user defined function to compute the lateral scattering in the CEM.

Semi-analytical model for range straggling

The range straggling is modeled by applying a convolution kernel to the fluence map (see previous section). The $\sigma(E)$ of the convolution kernel as a function of the proton energy is described by a semi-analytical model (a polynomial):

$$\sigma(E) = \sum_i a_i . E^i$$

In this section, the procedure to determine the value of the polynomial coefficients $a_i$ is described.

The principle to determine the polynomial coefficient is to compare the integrated depth dose (IDD) computed with the analytical model described in the previous section (which depends on the $a_i$ ) with the IDD computed from a Monte Carlo simulation. The comparison is repeated for differnet geometries. The polynomial coefficient are iteratively optimised so that the difference between the semi-analytical IDD and the Monte Carlo IDD is minimum.

The comparison was done for the following geometries :

PTV size (mm)	Range shifter thickness (mm)	Energy layer spacing (mm)
20 * 20 * 50	50	5
20 * 20 * 50	100	5
20 * 20 * 50	150	5
20 * 20 * 20	50	2
20 * 20 * 20	50	5
20 * 20 * 70	50	5

The range shifter was Lexan. The target was a water phantom. The PTV is a parallelepiped. The CEM were pyramids.

The dose map in the water phantom was computed using MCsquare. The integrated depth dose profile $IDD_{MC}(z)$ is computed by integrating, at each depth z, the dose in the whole plane orthogonal to the proton beam.

The analytical dose map is computed by first computing the fluence $F_{E,a_i}^{strg} (x,y)$ at different energies using the analytical method described in the previous section. These fluence map depends on the choice of the polynomial coefficients $a_i$. The fluence for the beamlet of energy E is the average of the fluence map at the correspodning energy. The depth dose profile (per unit of fluence) $Bort(z)$ is described by the Bortfeld model of Breagg peak. The analytical IDD $IDD_{a_i}(z)$ is then the weighted sum of the Bortfeld profiles:

$$IDD_{a_i}(z) = \sum_E Bort(z) . \int \int F^{strg}_{E,a_i}(x,y) .dx.dy$$

where z is the depth paralell to the proton beam axis. The analytical IDD depends on the choice of the polynomial coefficients $a_i$.

The objective function to optimise is:

$$GOF(a_i) = \sum_z (IDD_{a_i}(z) - IDD_{MC}(z))^2$$

The polynomial coefficients $a_i$ are varies until the minimum of $GOF(a_i)$ is reached.

CEM geometries

In the next paragraph, we will describe the geometries that can be used in one cell of the CEM.

Pyramid

The parameter SpikeType is up. In one cell, there is one single "aztec" pyramid made of the stacking of square or hexagonal blocks (depending on the parameter GridLayout) of decreasing radius. There is one pyramid per cell, i.e. one pyramid per PBS spot. The width of each terrace is optimiszd to deliver the correct weight $\omega_{s,E}$ of proton with energy E.

In one PBS spot, the proton fluence is the highest at the center of the PBS spot, i.e. at the position of the highest step of the pyramid. The highest step of the pyramid will give the Bragg peak with the smallest range in the SOBP, which is also a Bragg peak with a small weight in the SOBP. This is unfortunate as this will result in a very small area attributed to the highest step of the pyramid. The tip of the pyramid will therefore be very narrow and mechanically fragile. In addition, a small printing error in the pyramid will significantly increase the weight of the proximal Bragg peak. This pyramid geometry is therefore not very robust to printing errors.

FIGURE 22: There is one Aztec pyramid per cell (i.e. one pyramid per PBS spot). The width of each step of the pyramid is optimized to give the correct weight to each Bragg peak in the SOBP. The highest step (corresponding to the proximal Bragg peak, with small weight) is at the centre of the cell, in the region of highest fluence. Note the large height of the base of the CEM which is required to provide the additional range shifting not possible with an integral number of range shifter slabs.

Chimney

The parameter SpikeType is ellipse. In one cell, there is one single chimney made of the stacking of elliptical blocks. The heigh in the centre of the chimney is smaller than the height on the periphery of the chimney. The width of each terrace is optimized to deliver the correct weight $\omega_{s,E}$ of proton with energy E.

In one PBS spot, the proton fluence is the highest at the centre of the PBS spot, i.e. at the position of the shallowest step of the chimney. The shallowest step of the chimney will give the Bragg peak with the largest range in the SOBP, which is also a Bragg peak with a largest weight in the SOBP. Therefore, the Bragg peak with the largest weight will corresponds to a step located in the region of highest fluence in the PBS spot. This gives a good match between proton fluence and spot weight. This makes the chimney more robust to printing erros than the aztec pyramids.

FIGURE 23: There is one chimney per cell (i.e. one chimney per PBS spot). The width of each step of the chimney is optimized to give the correct weight to each Bragg peak in the SOBP. The highest step (corresponding to the proximal Bragg peak, with small weight) is at the periphery of the cell, in the region of smallest fluence. To compensate the reduced fluence, the width of this step is increased, therefore making it mechanically more robust.

The chimney has a shape matching the PBS spot shape. The diameter of the chimney matches the diameter of the PBS spot. The main axis of the elliptical chimney is aligned with the main axis of the PBS spot. In addition, the central axis of the chimney is aligned with the axis of the beamlet. The beamlets are not paralell to each other. They are diverging from a common source located at the centre of the magnet (in fact, there is one source for the X axis and one source for the Y axis). The axis of the chimney is also tilted: the axes are divergent with the same source as the proton beamlets.

FIGURE 24: Section through the virtual CT of the CEM in a plane orthogonal top the proton beam axis. The dose map of one PBS spot is overlaid as a colour map.

Fractal

The parameter SpikeType is fractal. In one cell, there are several columns. All columns have the same base area (typically 1x1mm or 2x2mm). The number of columns along a side of the cell is defined by the parameter NbColumns. There will be several sets of columns. In one set, all the columns have the same height. The number of columns in one set is optimised to deliver the correct weight $\omega_{s,E}$ of protons with energy E. The columns of a set are randomly placed in the cell to make the design more robust to error in the relative position of the PBS spot with the CEM.

FIGURE 25: There are several columns per cell (i.e. several column per PBS spot). The number of columns with a given height is optimized to give the correct weight to each Bragg peak in the SOBP. The columns all have the same base area. The number of columns of a specific height depends on the weight of the Bragg peak. The position of the columns in a cell is random.

Export CEM shape

If the parameter makeSTL is set to true, then an STL file is generated with the CEM and is saved in the folder 'output_path'.

Compute high resolution dose maps

At this point, the scanning trajectory, the spot weight (of the mono energy layer plan) and the CEM shape are optimized. It is now possible to run a Monte Carlo dose simulation through the full stack of elements to predict the final dose 3D distribution. The following geometry is used for the Monte Carlo simulation:

• The original CT scan is interpolated at the high resolution defined in the parameter intrpCTpxlSize. The resolution must be high enough to correctly describe the small structures of the CEM in this virtual CT scan. It is not possible to properly pixelize the CEM at the resolution of the initial CT scan because its voxels are too large.
• The aperture block is virtually added to the CT scan by defining the HU of brass in the corrsponding voxels. The distance between the isocentre and the downstream side of the aperture block was already computed.
• The range shifter slabs are modeled in the virtual CT scan by defining the HU of aluminum in the corresponding voxels of the CT scan. The gap between the range shifter slabs is modeled in the CT scan.
• The CEM is modeled by defining the HU of the CEM material in the corresponding voxels of the CT scan. The residual range shifting is done by adjusting the thickness of the base of the CEM. The range compensator (that was used in the mother IMPT plan) is merged with the CEM: the surface of the range compensator defines the base onto which the spike of the CEM are added (see Figure 2).
• The distances between the different elements match the mechanical drawings of the FLASH snout. The snout extension is defined by the position of the downstream side of the aperture block.

All the elements (CEM, range shifter, aperture) are virtually inserted in the CT scan. MCsquare will model the particles trajectory and the dose deposition directly in the high resolution CT scan. No additional range shifter is added in the beam line model in MCsquare.

The resolution of the CT scan (parameter intrpCTpxlSize) is typically of the order of 0.2mm. If the full CT scan was converted at high resolution, this would use a very large amount of computer memory. In order to avoid out-of-memory errors, the initial CT scan is divided in sub-volumes centered on the axis of each PBS beamlet. There is one sub-volume per beamlet. The pixels axis of the sub-volume are aligned with the proton beam axis. The width of the sub-volume is sufficiently large to include the full cross section of the beamlet. The length of the sub-volume starts at the base of the CEM and extends beyond the distal surface of the PTV to include the complete volume in which the beamlet will deliver dose. The dose distribution is computed for each sub-volume (and hence PBS spot) separately.

The each high resolution CT scan is saved in the folder /Outputs/CEM_beam1/CT_with_CEM_spt_X_0_Y_15/. There is one folder per sub-volume. The name of the folder corresponds to the spot position. The axes of the Ct scan are aligned with the axes of the IEC gantry coordinate system.

FIGURE 26: Sub-volume of the high resolution CT scan around one PBS beamlet. The axes of the CT scan of the subvolume are aligned with the axes of the IEC gantry coordinate system. The CEM, the range shifter and the aperture are inserted in the CT scan. The slabs of the range shifter are modeled in the virtual CT scan. The distances between the elements follow the mechanical drawings of the FLASH snout. The shape of the range compensator is visible as a curvature at the base of the spikes. The thick base of the CEM acts as a residual range shifter.

The dose in each sub-volume is computed using the MCsquare Monte Carlo dose engine. If the dose was tallied in the 0.2x0.2x0.2mm voxels, a very large number of protons (parameter protonsHighResDose) would need to be included in the MC simulation in order to reduce the statistical noise. However, there is no need to have a dose map with a resolution of 0.2x0.2x0.2mm. A typical dose map resolution is of the order of the CT scan voxels (2x2x2mm). Therefore, during the Monte Carlo simulation, while the particle trajectory will take place in (a sub-volume of) the high resolution virtual CT scan (at 0.2x0.2x0.2mm resolution), the dose will be tallied using larger voxels (defined by parameter CEFDoseGrid). This allows to reduce the number of protons in the MC simulation (parameter protonsHighResDose) and therefore reduce the computation time.

The high resolution dose map is saved in the file dose_HighRes_spt6_spt_X_0_Y_15.dcm in the folder /Outputs/CEM_beam1/Outputs//. There is one file per sub-volume. The name of the file corresponds to the spot position.

If the parameter SaveDoseBeamlets is set to TRUE, then an additional high resolution dose map of each PBS spot is saved in a DICOM file with the same coordinate system as the low resolution CT scan. This allows the researcher to overlay each beamlet with the original CT scan. In addition, a text file is stored with the spot timing. This allows the researchers to re-compute dose rate using different definitions.

The dose map of the sub-volume is re-sampled at the resolution of the initial CT scan and it is added to the dose from the other sub-volumes. This creates the final dose map for the protons having gone through the CEM, range shifter, aperture and CT scan. If the parameter ComputeCEMdose = true , the resulting dose map is saved in the file Outputs/Dose_withCEF_beam1.dcm in the /Output folder.

Compute dose rate statistics

At this point, the spot trajectory (and hence the timing of spot delivery) and the spot weights are known. The dose influence matrix of each PBS spot is also kwnon. For each voxel in the CT scan, it is now possible to compute the timing of dose delivery. One such plot can be drawn per voxel in the CT scan.

Using the timing of dose delivery, it is possible to compute, in each voxel, different parameters:

• Percentile dose rate : This paragraph definess this dose rate. The results are saved in DICOM files in the output folder.
• Dose averaged dose rate : This paragraph define this dose rate. The results are saved in DICOM files in the output folder.

the ConformalFLASH library saves the dose rate maps in the /Output folder.