Adaptive planning strategy for high dose rate prostate brachytherapy—a simulation study on needle positioning errors

Throughout Task 1 and Task 5, there is the requirement to select automatically the best dose plan among several others. Therefore the quality of a dose plan must be quantified. For that, the energy E proposed by Borot de Battisti et al (2015) defined as:

$$\begin{eqnarray}&&E=\min \left(A,B,C,D\right)\end{eqnarray} \tag{ 5 }$$

with:

$$\begin{eqnarray}&&A={{D}_{95\%~\text{PTV}}}-D_{95\%~\text{PTV}}^{\text{min}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&B=D_{10\%~\text{Ur}}^{\text{max}}-{{D}_{10\%~\text{Ur}}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&C=D_{1~\text{cc}~\text{Rec}}^{\text{max}}-{{D}_{1~\text{cc}~\text{Rec}}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&D=D_{1~\text{cc}~\text{Bl}}^{\text{max}}-{{D}_{1~\text{cc}~\text{Bl}}}\end{eqnarray} \tag{ 9 }$$

is used. The parameters A, B, C and D represent the differences between the dose coverage parameters and the clinical constraints (set as input) of the PTV, the urethra, the rectum and the bladder, respectively. For example, if A is positive, the dose plan reached the clinical prescribed dose of the PTV. Conversely, if A is negative, the clinical prescribed dose of the PTV was not achieved. More precisely, the higher the value of A, the higher the dose coverage of the PTV becomes. Consequently, the minimum value over A, B, C and D in equation (5) corresponds to the 'worst' error in dose deposition between the different regions of interest (PTVs and OARs). E represents therefore the quality of the dose plan: the larger E is, the better the dose plan in terms of quality (if E  >  0 all clinical constraints are fulfilled, conversely if E  <  0 they are not).

Let $\boldsymbol{r}_{{\text{rot}}}=\left({{x}_{\text{rot}}},{{y}_{\text{rot}}},{{z}_{\text{rot}}}\right)$ be the position of the center of rotation of the setup and $\left({{\theta}_{i}},{{\phi}_{i}}\right)$ be the angle of the ith needle track in the spherical coordinate system. In common practice of HDR brachytherapy, the distance $\Delta$ between the dwell positions along the needle is constant. Due to the finite size of the needle, the number ${{N}_{\text{source}}}$ of possible dwell positions of the source along the needle is limited. The dwell positions $\boldsymbol{r}_{k}^{i}\left(\boldsymbol{r}_{{\text{rot}}},{{\theta}_{i}},{{\phi}_{i}}\right),~k\in \left[1,{{N}_{\text{source}}}\right]$ of the ith needle insertion can then be expressed in the Cartesian coordinates as follows:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\boldsymbol{r}_{k}^{i}\left(\boldsymbol{r}_{{\text{rot}}},{{\theta}_{i}},{{\phi}_{i}}\right)=\left[\begin{array}{*{35}{l}} {{x}_{\text{rot}}}+k\Delta \sin \left({{\theta}_{i}}\right)\cos \left({{\phi}_{i}}\right) \\ {{y}_{\text{rot}}}+k\Delta \sin \left({{\theta}_{i}}\right)\sin \left({{\phi}_{i}}\right) \\ {{z}_{\text{rot}}}+k\Delta \cos \left({{\theta}_{i}}\right) \end{array}\right].\end{array}\end{eqnarray} \tag{ 10 }$$

To calculate the initial dose plan for a given number ${{N}_{\text{needle}}}$ of divergent needle insertions, the following parameters must be optimized.

• The position of the center of rotation of the setup $\boldsymbol{r}_{{\text{rot}}}$.
• The angle of the ith needle track $\left({{\theta}_{i}},{{\phi}_{i}}\right)$ ($i\in \left[1,{{N}_{\text{needle}}}\right]$).
• The dwell time $t_{k}^{i}$ of the source position k of the ith needle track ($i\in \left[1,{{N}_{\text{needle}}}\right]$ and $k\in \left[1,{{N}_{\text{source}}}\right]$).

Borot de Battisti et al (2015) described a way to determine those parameters, but other methods can also be applied.

At Task 2, the needle that is about to be inserted is selected by the user. In our simulations, we try to mimic a strategic intent of the user to place the most sensitive needle track first, such that potential errors can still be compensated by the re-optimization of the subsequent needles. To determine the most sensitive needle track, a criterion that quantifies the impact of an inaccurate needle insertion on the dose plan quality is introduced. More precisely, the needle track chosen is the one that may have the most negative impact on the dose plan quality if a deviation occurs during insertion of the needle: inserting the needle into this track leaves the best possibility to compensate for needle placement error by re-optimizing the dose plan for the remaining needle tracks.

Practically, the most sensitive needle track can be selected using a stochastic approach based on simulations as follows.

For all remaining needle tracks i ($i\in \left[1;{{N}_{\text{tracks}}}\right]$):

• (1)  
  Starting from the current dose plan, a random change of angulation of needle track i is simulated. A random angulation change modeled by a Gaussian distribution with a standard deviation of 0.015 rad was found to be a good representation of a typical insertion error (0.015 rad corresponds to a deviation of 3 mm at a distance of 200 mm from the center of rotation).
• (2)  
  Without any re-optimization, the value of E corresponding to the resulting dose plan is calculated.
• (3)  
  Steps (1) and (2) are repeated ${{N}_{\text{angle}}}$ times (with ${{N}_{\text{angle}}}$ an integer big enough to be statistically acceptable). After investigation, a typical value of ${{N}_{\text{angle}}}=100$ was found to be a good compromise between speed and accuracy.
• (4)  
  The 5th percentiles (noted ${{E}_{5\%}}(i)$) of the values of E obtained in step (2) are calculated: this quantifies how potentially negative the impact of an inaccurate needle insertion in needle track i is on the dose plan quality.

The selected most sensitive needle track is therefore ${{i}_{\text{selected}}}$ such that:

$$\begin{eqnarray}&&{{i}_{\text{best}}}=\underset{i\in \left[1;{{N}_{\text{tracks}}}\right]}{\mathop{\text{argmin}}}\left[{{E}_{5\%}}(i)\right].\end{eqnarray} \tag{ 11 }$$

At Task 5, the recently inserted needle underwent a potential deviation as compared to the targeted dose plan. A re-optimization of the remaining parameters (angles of the needle tracks and dwell times) is thus mandatory to compensate for this deviation and to achieve the clinical prescription. This section presents the method of re-optimization of these remaining parameters.

Supposing that the recently inserted needle is associated to the ${{N}_{\text{ins}}}$th needle insertion (${{N}_{\text{ins}}}\leqslant {{N}_{\text{needle}}}$), the following parameters need to be re-optimized:

• The angle of the remaining needle tracks ${{\left({{\theta}_{i}},{{\phi}_{i}}\right)}_{{{N}_{\text{ins}}}<i\leqslant {{N}_{\text{needle}}}}}$.
• The dwell times $\left(t_{k}^{i}\right)_{{{N}_{\text{ins}}}\leqslant i\leqslant {{N}_{\text{needle}}}}^{1\leqslant k\leqslant {{N}_{\text{source}}}}$ of the needle which has just been inserted ($i={{N}_{\text{ins}}}$) and the dwell times of the remaining needle tracks (${{N}_{\text{ins}}}<i\leqslant {{N}_{\text{needle}}}$).

Let $p=({{\theta}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}},{{\phi}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}},t_{1\cdots {{N}_{\text{source}}}}^{{{N}_{\text{ins}}}},\cdots,t_{1\cdots {{N}_{\text{source}}}}^{{{N}_{\text{needle}}}})$ be the vector containing the parameters of the setup to re-optimize and $\Omega$ its corresponding set of feasible solutions. According to the guidelines of AAPM Task Group No. 43 (Nath et al 1995, Rivard et al 2004), the dose $D(\boldsymbol{r},p)$ received at $\boldsymbol{r}$ can be expressed as the sum of the contribution of all source positions:

$$\begin{eqnarray}&&D\left(\boldsymbol{r},p\right)={{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)+\underset{i={{N}_{\text{ins}}}}{\overset{{{N}_{\text{needle}}}}{\mathop \sum}}\underset{k=1}{\overset{{{N}_{\text{source}}}}{\mathop \sum}}\text{d}_{k}^{i}\left(\boldsymbol{r},\boldsymbol{r}_{{\text{rot}}},{{\theta}_{i}},{{\phi}_{i}}\right)t_{k}^{i}\end{eqnarray} \tag{ 12 }$$

where ${{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)$ is the already delivered dose to the patient.

A way to determine the optimal parameters ${{p}_{\text{optimal}}}\in \Omega$ is to approach the dose to a given value ${{D}_{\text{opt}}}\left(\boldsymbol{r}\right)$ for all points $\boldsymbol{r}$. The determination of ${{p}_{\text{optimal}}}$ is therefore a multi-objective optimization and can be determined by solving the following equations:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{p}_{\text{optimal}}}=\underset{p\in \Omega}{\mathop{\text{argmin}}}\left[C\left(\,p\right)\right]\\ \text{with} ~C\left(\,p\right)=\mathop{\int}^{}\mathop{\int}^{}{{\mathop{\int}^{}}_{\boldsymbol{r}\in \mathbb{R}}}\omega \left(\boldsymbol{r}\right){{\left[D\left(\boldsymbol{r},p\right)-{{D}_{\text{opt}}}\left(\boldsymbol{r}\right)\right]}^{2}}\text{d}\,\boldsymbol{r}\end{array}\end{eqnarray} \tag{ 13 }$$

where C(p) is the cost function and $\omega \left(\boldsymbol{r}\right)$ is a spatial weight function (the latter will be detailed in section 2.4.2).

The strategy for the resolution benefits from the linear impact of the dwell times on the deposited dose. If $\left({{\theta}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}},{{\phi}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}}\right)$ are fixed, the determination of $t_{k}^{i}$ reverts to solving a set of linear equations and it is thus feasible to find a direct and efficient solution. However, solving equation (13) for the variables of the needle positioning $\left({{\theta}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}},{{\phi}_{{{N}_{\text{ins}}}+1\cdots {{N}_{\text{needle}}}}}\right)$ is highly non-convex and consequently difficult to overcome. The workflow of the dose plan re-optimization is depicted in figure 2. The proposed method relies on the determination of the optimal angulations of the remaining needle tracks which are deduced using heuristic searches (Task 5a) while the optimal dwell times are determined by the resolution of linear equations (Task 5b).

Zoom In Zoom Out Reset image size

2.4.1. Re-optimization of the position of remaining needle positions (Task 5a).

In this section, the heuristic approach to re-optimize the remaining needle positions of the dose plan is described. We seek a uniform distribution in space of needle tracks that will cover the PTV without crossing the urethra.

For that, the PTV and the urethra are projected from the center of rotation on a transverse plane behind the PTV (see figure 3). The projection of the needle tracks on the same plane can be represented by ${{N}_{\text{needle}}}$ points. To determine a uniform distribution of points, an extension of the procedure described by Borot de Battisti et al (2015) is applied: k-means clustering with the additional constraints of immobility of the cluster centers corresponding to the tracks where the needle was already inserted. Figure 3 shows a typical example of re-optimization of the needle track positions using this method.

Zoom In Zoom Out Reset image size

Consequently, with this method, we can re-optimize the positions of the needle tracks in order to obtain a distribution that is as uniform as possible and the errors in angulation during insertion of the needle can be thus compensated.

2.4.2. Re-optimization of the remaining dwell times (Task 5b).

A planning study was performed to test our re-optimization procedure. The approach consists in comparing the dose distribution obtained in the scenario with and without feedback of the needle position. For that, complete single fraction HDR brachytherapy procedures with and without the proposed feedback strategy were simulated. In the situation without feedback, no re-optimization of the dose plan was done after the insertion of the needle (i.e. Task 5 was skipped).

The simulations were done on eight patients, with various numbers of needle insertions (${{N}_{\text{needle}}}$ varying from 4 to 12). The anatomy of the patients was obtained using the delineations of the prostate tumors and the OARs considered (urethra, bladder, rectum and rest of the tissues) on 1 mm³resolution MR images by an experienced oncologist.

For all patients, the clinical constraints were: $D_{95\%~\text{PTV}}^{\text{min}}=19$ Gy, $D_{10\%~\text{Ur}}^{\text{max}}=21$ Gy, $D_{1~\text{cc}~\text{Rec}}^{\text{max}}=12$ Gy and $D_{1~\text{cc}~\text{Bl}}^{\text{max}}=12$ Gy. Those constraints are in line with the study of Hoskin et al (2014) and Prada et al (2012) who showed that HDR brachytherapy as monotherapy is feasible with acceptable levels of acute complications by delivering a single fraction dose of 19 Gy to the target. All the MR images used for the simulations of this study are intra-operative images of patients who underwent HDR prostate brachytherapy as monotherapy in a single fraction at the UMCU. These were patients with localized prostate cancer, a prostate-specific antigen level lower than 10 ng ml⁻¹ and a Gleason score of 7 or less. The PTV volumes ranged from 20.7 cm³ to 31.3 cm³ with a median of 24.75 cm³. The acquisition of patient MR images used in this study was approved by the Institutional Review Board.

To simulate the positioning error during the insertion of the needle, we modeled it by a random angulation error described by a Gaussian distribution with a standard deviation of 0.015 rad. Moreover, we supposed that no error occurred in the measurement of the needle position. In that way, only the contribution of the proposed feedback strategy on the needle positioning errors is assessed.

Source positions were chosen according to the common procedure at the UMCU: for each needle insertion, the active source center positions were separated by a step-size of $\Delta =2.5$ mm and situated inside the PTV with margin of 3 mm.

To simulate the irradiation and calculate the dose delivered to the patient, the dose rate was calculated using the point source approximation model because of its minimum time of computation, with a small adaptation as follows to avoid over-optimization of the dose close to the source:

$$\begin{eqnarray}&&{{d}_{i}}\left(\boldsymbol{r}\right)={{S}_{K}}\Lambda{{g}_{P}}\left[{{R}_{i}}\left(\boldsymbol{r}\right)\right]{{\Phi}_{\text{an}}}\left[{{R}_{i}}\left(\boldsymbol{r}\right)\right]\frac{{{R}_{0}}^{2}}{{{R}_{i}}{{\left(\boldsymbol{r}\right)}^{2}}+\exp \left[-{{R}_{i}}{{\left(\boldsymbol{r}\right)}^{2}}\right]}\end{eqnarray} \tag{ 14 }$$

where S_K is the air-kerma strength, $\Lambda$ the dose-rate constant in water, ${{\Phi}_{\text{an}}}(R)$ the one-dimensional anisotropy factor, R₀ the reference distance which is specified to be 10 mm, g_P(R) to the radial dose function in the case of the point source approximation model and ${{R}_{i}}\left(\boldsymbol{r}\right)$ the distance (in millimeters) between the ith source position at coordinate $\boldsymbol{r}_{{i}}$ and $\boldsymbol{r}$ (${{R}_{i}}\left(\boldsymbol{r}\right)=\,\Vert \boldsymbol{r}_{{i}}-\boldsymbol{r}{{\Vert}_{2}}$). With this adaptation of the point source model, the dose has an upper limit value close to the source, therefore reducing the numeric instabilities for ${{R}_{i}}\left(\boldsymbol{r}\right)$ approaching 0. TG43 constants, anisotropy factor and radial dose function for the microSelectron-HDR (Elekta/Nucletron, Veenendaal, The Netherlands) 192-Iridium source were taken from a study of Daskalov et al (1998) ($\Lambda=1.108$ cGy.h⁻¹.U⁻¹) and an arbitrary source strength S_K  =  40.80 mGy.h⁻¹.m² was chosen. The multiplication of the radial dose function g_P(R) and the anisotropy factor ${{\Phi}_{\text{an}}}(R)$ was approximated by a 2nd order polynomial fit (${{g}_{P}}(R)\cdot {{\Phi}_{\text{an}}}(R)={{a}_{0}}+{{a}_{1}}R+{{a}_{2}}{{R}^{2}}$). The coefficients for the fit were a₀  =  1.11, ${{a}_{1}}=-3.30\cdot {{10}^{-3}}$ and ${{a}_{2}}=3.12\cdot {{10}^{-6}}$, where R is in millimeters.

Due to the randomness of the needle deviations during insertions, both simulations with and without feedback were repeated 100 times to verify the robustness of both procedures.

The dose parameters obtained at the end of the procedure were compared to the aforementioned clinical constraints to verify if the dose distributions were clinically acceptable.

As a typical example of Task 2, patient 1 of the planning study with six needle insertions (${{N}_{\text{needle}}}=6$) is detailed in figure 4. We first focus our interest on the specific instant of the interventional procedure for which the second insertion of the needle must be simulated (the needle track corresponding to the first needle insertion is noted as $\text{needle}~\text{track}~1$). The most sensitive needle track of the current dose plan must be selected. For that, the method described in section 2.3 is applied in order to predict the impact of an inaccurate needle insertion for each individual needle track. The value of E corresponding to each individual needle track ($\text{needle}~\text{tracks}~2$ to 6) is presented in figure 4(a).

Zoom In Zoom Out Reset image size

The distribution of energy values is different from one needle track to another. For this case, the median(minimum,maximum) of energy corresponding to $\text{needle}~\text{tracks}~2$, 3, 4 and 5 are 0.9(−0.4, 1.2) Gy, 0.5(−1.3, 1.2) Gy, 0.7(−0.9, 1.2) Gy, 1.0(0.0, 1.2) Gy and 1.1(0.9, 1.2) Gy respectively. The 5th percentile of energy was 0.1 Gy,−0.8 Gy,−0.5 Gy, 0.7 Gy and 1.0 Gy for $\text{needle}~\text{tracks}~2$, 3, 4 and 5 respectively. The energy value in case of insertions without error was 1.2 Gy. The most sensitive needle track is therefore $\text{needle}~\text{track}~3$ (conversely, $\text{needle}~\text{track}~6$ is the least sensitive needle track). Consequently, the needle track selected for insertion is $\text{needle}~\text{track}~3$ in that case.

For all tested patients, the percentages of clinically acceptable plans for the procedures with and without feedback are presented in table 1. When the value is not presented (for example patient 5 for ${{N}_{\text{needle}}}=4$), it corresponds to the case where the initial dose plan (without errors in needle placement) obtained at Task 1 did not achieve the clinical constraints: in those cases there was no point in performing the experiment since the dose plan quality will only decrease due to error in needle positioning at each insertion. In all other cases, the dose plan without error in needle placement achieved the clinical constraints.

The results presented in table 1 show that in 36 of the 37 cases tested ($97\%$), the number of clinically acceptable plans obtained at the end of the procedure was larger with feedback as compared to the scenario without feedback. The only exception was for patient 5 with ${{N}_{\text{needle}}}=6$: the values in that case were quite similar ($7\%$ and $9\%$ with and without re-optimization respectively). Furthermore, the percentage of clinically acceptable dose plans tended to increase with the number of needle insertions employed in both situations with and without re-optimization.

The dose parameters of two typical examples (patient 1 with ${{N}_{\text{needle}}}=6$ and patient 4 with ${{N}_{\text{needle}}}=4$) are depicted in figures 5(a) and (b). Typical dose distributions and the dose-volume histograms (DVHs) obtained from the simulations are shown in figures 5(c)–(f) for patient 1 with ${{N}_{\text{needle}}}=6$. Both the situations with and without re-optimization are presented.

Zoom In Zoom Out Reset image size

For patient 1 with ${{N}_{\text{needle}}}=6$, the dose parameters (${{D}_{95\%~\text{PTV}}}$, ${{D}_{10\%~\text{Ur}}}$, ${{D}_{1~\text{cc}~\text{Rec}}}$ and ${{D}_{1~\text{cc}~\text{Bl}}}$) of the dose plan without error in needle placement were 20.3 Gy, 12.7 Gy, 10.8 Gy and 4.8 Gy respectively. The distributions (median(minimum,maximum)) of ${{D}_{95\%~\text{PTV}}}$, ${{D}_{10\%~\text{Ur}}}$, ${{D}_{1~\text{cc}~\text{Rec}}}$ and ${{D}_{1~\text{cc}~\text{Bl}}}$ in the case with error in needle placement were 19.4(14.9, 21.3) Gy, 12.7(11.0, 15.7) Gy, 10.9(8.9, 14.1) Gy and 4.8(4.3, 5.4) Gy respectively without re-optimization and 19.9(19.4, 20.3) Gy, 13.2(11.4, 18.1) Gy, 11.0(10.6, 11.6) Gy and 4.8(3.8, 6.4) Gy respectively with re-optimization. Moreover, without re-optimization, $53\%$ of the dose plans with error in needle placement did not achieve the clinical constraints whereas $100\%$ of them did with re-optimization (those values are reported in table 1). The example presented in figures 5(d)–(f) show better PTV coverage for the plans with reoptimization compared to the plans without.

For patient 4 with ${{N}_{\text{needle}}}=4$, the dose parameters of the dose plan without error in needle placement were 20.2 Gy, 19.1 Gy, 10.9 Gy and $10.4''$ Gy respectively. The distributions (median(minimum,maximum)) of ${{D}_{95\%~\text{PTV}}}$, ${{D}_{10\%~\text{Ur}}}$, ${{D}_{1~\text{cc}~\text{Rec}}}$ and ${{D}_{1~\text{cc}~\text{Bl}}}$ in the case with error in needle placement were 19.4(16.9, 20.9) Gy, 19.5(15.2, 26.8) Gy, 10.9(9.4, 12.7) Gy and 10.5(9.7, 11.3) Gy respectively without re-optimization and 19.7(18.4, 20.3) Gy, 19.4(15.8, 21.4) Gy, 11.1(9.9, 12.6) Gy and 10.8(10.4, 11.7) Gy respectively with re-optimization. Without re-optimization, $42\%$ of the dose plans with error in needle placement were clinically acceptable without re-optimization, as compared to $94\%$ for the scenario with re-optimization (values depicted in table 1).

To find $({{t}_{1}}\cdots {{t}_{N_{\text{source}}^{\text{opt}}}})$ which minimizes C(p) with the additional condition that for all $n\in \left[1,N_{\text{source}}^{\text{opt}}\right]$, ${{t}_{n}}\geqslant 0$, we modify the cost function of the OAR as follows:

$$\begin{eqnarray}&&{{C}_{\text{OAR}}}\left(\,p\right)=\frac{{{\omega}_{\text{OAR}}}}{{{V}_{\text{OAR}}}}\underset{\boldsymbol{r}\in \text{OAR}}{\mathop \sum}\left(\underset{n=1}{\overset{N_{\text{source}}^{\text{opt}}}{\mathop \sum}}{{\left[{{d}_{n}}\left(\boldsymbol{r}\right){{t}_{n}}\right]}^{2}}+\underset{n=N_{\text{source}}^{\text{opt}}+1}{\overset{N_{\text{source}}^{\text{tot}}}{\mathop \sum}}{{\left[{{d}_{n}}\left(\boldsymbol{r}\right){{t}_{n}}\right]}^{2}}\right).\end{eqnarray} \tag{ A.5 }$$

By modifying the OAR cost function, ${{C}_{\text{OAR}}}\left(\,p\right)$ will not be null through 'destructive interference' effects between the dwell times and most of the unphysical negative results solutions are therefore excluded.

From now, we note ${{\tilde{\omega}}_{\text{PTV}}}=\frac{{{\omega}_{\text{PTV}}}}{{{V}_{\text{PTV}}}}$, ${{\tilde{\omega}}_{\text{OAR}}}=\frac{{{\omega}_{\text{OAR}}}}{{{V}_{\text{OAR}}}}$ and the already delivered dose ${{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)=\underset{n=N_{\text{source}}^{\text{opt}}+1}{\overset{N_{\text{source}}^{\text{tot}}}{\mathop \sum}}{{d}_{n}}\left(\boldsymbol{r}\right){{t}_{n}}$. We then develop the same demonstration described by Goldman et al (2009). The minimum of C(p) is obtained by the following condition:

$$\begin{eqnarray}&&\forall j\in \left[1,N_{\text{source}}^{\text{opt}}\right],~\frac{\partial C\left(\,p\right)}{\partial {{t}_{j}}}=0.\end{eqnarray} \tag{ A.6 }$$

Using equation (A.4), we obtain $\forall j\in \left[1,N_{\text{source}}^{\text{opt}}\right]$:

$$\begin{eqnarray}&&2{{\tilde{{\omega}}}_{\text{PTV}}}\underset{\boldsymbol{r}\in \text{PTV}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right)\left[\underset{n=1}{\overset{N_{\text{source}}^{\text{opt}}}{\mathop \sum}}{{d}_{n}}\left(\boldsymbol{r}\right){{t}_{n}}+{{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)-{{D}_{\text{PTV}}}\right]+2{{\tilde{{\omega}}}_{\text{OAR}}}\underset{\boldsymbol{r}\in \text{OAR}}{\mathop \sum}{{d}_{j}}{{\left(\boldsymbol{r}\right)}^{2}}{{t}_{j}}=0.\end{eqnarray} \tag{ A.7 }$$

By exchanging the order of summations, equation (A.7) becomes:

$$\begin{eqnarray}\begin{array}{*{20}{l}}{{\tilde{{\omega}}}_{\text{PTV}}}\underset{\boldsymbol{r}\in \text{PTV}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right)&\left[{{D}_{\text{PTV}}}-{{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)\right]\\ &=\underset{n=1}{\overset{N_{\text{source}}^{\text{opt}}}{\mathop \sum}}{{t}_{n}}\left[{{\tilde{{\omega}}}_{\text{OAR}}}\underset{\boldsymbol{r}\in \text{OAR}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right){{d}_{n}}\left(\boldsymbol{r}\right){{\delta}_{nj}}+{{\tilde{{\omega}}}_{\text{PTV}}}\underset{\boldsymbol{r}\in \text{PTV}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right){{d}_{n}}\left(\boldsymbol{r}\right)\right]\end{array}\end{eqnarray} \tag{ A.8 }$$

where ${{\delta}_{nj}}$ is the Kronecker delta function. Therefore, the values of dwell times can be determined by a matrix inversion:

$$\begin{eqnarray}&&T={{\alpha}^{-1}}\beta \end{eqnarray} \tag{ A.9 }$$

where $T={{\left({{t}_{j}}\right)}_{j\in \left[1,N_{\text{source}}^{\text{opt}}\right]}}$ is the vector of dwell times (vector of $N_{\text{source}}^{\text{opt}}$ elements), α is a square symmetrical positive matrix of $\left[N_{\text{source}}^{\text{opt}},N_{\text{source}}^{\text{opt}}\right]$ elements, and β is a vector of $N_{\text{source}}^{\text{opt}}$ elements such that:

$$\begin{eqnarray}\begin{array}{*{20}{l}}\forall \left(\,j,n\right)\in &{{\left[1,N_{\text{source}}^{\text{opt}}\right]}^{2}},\\ &{{\alpha}_{jn}}={{\tilde{\omega}}_{\text{OAR}}}\underset{\boldsymbol{r}\in \text{OAR}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right){{d}_{n}}\left(\boldsymbol{r}\right){{\delta}_{nj}}+{{\tilde{\omega}}_{\text{PTV}}}\underset{\boldsymbol{r}\in \text{PTV}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right){{d}_{n}}\left(\boldsymbol{r}\right)\end{array}\end{eqnarray} \tag{ A.10 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}\forall j\in &\left[1,N_{\text{source}}^{\text{opt}}\right],\\ &{{\beta}_{j}}={{\tilde{\omega}}_{\text{PTV}}}\underset{\boldsymbol{r}\in \text{PTV}}{\mathop \sum}{{d}_{j}}\left(\boldsymbol{r}\right)\left[{{D}_{\text{PTV}}}-{{D}_{\text{deliv}}}\left(\boldsymbol{r}\right)\right].\end{array}\end{eqnarray} \tag{ A.11 }$$

It is interesting to note that the matrix inversion in equation (A.9) is low time consuming. It is therefore possible to find the optimal value of ${{\left({{\omega}_{\text{OA}{{\text{R}}_{l}}}}\right)}_{l\in \left[1,{{N}_{\text{OAR}}}\right]}}$, ${{\omega}_{\text{PTV}}}$ and ${{D}_{\text{PTV}}}$ using an exhaustive search. Therefore, we propose to calculate the dwell times using equation (A.9) for different combinations of ${{\left({{\omega}_{\text{OA}{{\text{R}}_{l}}}}\right)}_{l\in \left[1,{{N}_{\text{OAR}}}\right]}}$, ${{\omega}_{\text{PTV}}}$ and ${{D}_{\text{PTV}}}$ within the same ranges described by Borot de Battisti et al (2015), which are:

$$\begin{eqnarray}&&{{D}_{\text{PTV}}}\in \left\{20,21,\cdots,79,80\right\}\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&{{\log}_{10}}\left({{\omega}_{\text{PTV}}}\right)\in \left\{0,1,\cdots,9,10\right\}\end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray}&&\forall l\in \left[1,{{N}_{\text{OAR}}}\right],~{{\log}_{10}}\left({{\omega}_{\text{OA}{{\text{R}}_{l}}}}\right)\left\{0,1,\cdots,9,10\right\}\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&\underset{l=1}{\overset{{{N}_{\text{OAR}}}}{\mathop \sum}}{{\log}_{10}}\left({{\omega}_{\text{OA}{{\text{R}}_{l}}}}\right)+{{\log}_{10}}\left({{\omega}_{\text{PTV}}}\right)=10.\end{eqnarray} \tag{ A.15 }$$

The exhaustive search of ${{\left({{\omega}_{\text{OA}{{\text{R}}_{l}}}}\right)}_{l\in \left[1,{{N}_{\text{OAR}}}\right]}}$, ${{\omega}_{\text{PTV}}}$ and ${{D}_{\text{PTV}}}$ provides several dose plans. For each obtained dose plan, the metric E, which represents the quality of the dose plan, is calculated (see equation (5)). The optimal dose plan is then automatically determined by selecting the dose plan with the highest value of E.