A novel adaptive needle insertion sequencing for robotic, single needle MR-guided high-dose-rate prostate brachytherapy

High-dose-rate (HDR) prostate brachytherapy has gained increasing interest as an advanced treatment for patients with localized prostate cancer (Morton and Hoskin 2013, Hauswald et al 2016). It consists of inserting catheters into or close to the tumor and irradiating the localized cancer by stepping a radioactive source (e.g. Ir-192) through the catheters at various dwell positions for certain times according to a calculated dose plan (Venselaar et al 2013). In the current practice of HDR brachytherapy for localized prostate cancer, catheters are inserted manually under transrectal ultrasound (TRUS) guidance in a parallel configuration with the support of a template brachytherapy grid (Hoskin et al 2013). After completion of the implant procedure, needle reconstruction, dose planning and needle position verification are assessed using TRUS, CT or MRI. This procedure is not optimal for several reasons: firstly, US-guided manual insertion of needles may lead to sub-optimal needle configurations due to (1) anatomical modifications such as prostate deformations and swelling, (2) the shadowing effect behind the implant needles which decreases the image quality. Secondly, the procedure can be relatively time consuming (∼few hours) if two imaging modalities are used (TRUS for insertion of the needle and CT or MRI for the needle reconstruction).

Imaging, pathology and dose delivery studies have shown the superior value of multi-parametric MR-imaging for localization of prostate tumor (Atalar and Ménard et al 2005, Groenendaal et al 2010). Consequently, a complete MR-guided focal HDR prostate brachytherapy procedure is currently being investigated: to support needle insertion into the prostate with the patient in the MR bore, MR-compatible robotic methods have been developed recently in several institutes (Muntener et al 2006, Fischer et al 2007, 2008) and DiMaio et al (2007). In particular, a single needle MR-compatible robotic device is being developed at the University Medical Center Utrecht (UMCU) (van den Bosch et al 2010). This robotic device allows the automatic transperineal insertion of the needle into the prostate under MR-guidance. With this setup, the needle is inserted without a template, under different angles in a divergent pattern from a single rotation point situated at the perineum. Moreover, the dose is delivered needle-by-needle. A fast optimization planning method was recently proposed for this setup (Borot de Battisti et al 2015).

However, during the intervention, suboptimal dose delivery can occur due to two unpredictable events as follows:

• (i)  
  Sub-optimal needle positioning can happen due to uncertainties of the robotic needle placement or bending of the needle: Strassmann et al (2011) showed, on a prostate model, that the average needle positioning accuracy was $2.7\pm 0.7$ mm and $1.8\pm 0.6$ mm in the two following scenarios: (1) when the needle is manually inserted by the doctor with support of a template and (2) when using a robot-assisted method, respectively.
• (ii)  
  Intra-operative internal organ motion such as swelling, displacement (Stone et al 2002) or rotation (Lagerburg et al 2005) of the prostate (related to the trauma caused by the needle insertion) or intra-procedural changes in rectum or bladder filling may also induce sub-optimal target coverage or overdose of organ at risk.

In practice, sub-optimal dose delivery cannot be measured directly: the dose deposition relies on model-based dose calculation methods (Beaulieu et al 2012) which requires, as input, tumor, Organs at risk (OARs), needle and source localization. The ultimate goal is therefore to develop a fully automatic control system, where the dose plan and the needle insertion sequence are re-optimized during the intervention according to the two aforementioned perturbations. This fully automatic control system would require (1) to localize accurately the tumor and OARs which can be provided by the MRI (e.g. in current practice of HDR prostate brachytherapy at the UMCU, after insertion of catheters under TRUS guidance, the anatomy of the patient is imaged and delineated using MRI for dose planning), (2) to localize accurately the needle with MRI or other MR compatible localization methods such as fiber Bragg gratings-based technology (Borot de Battisti et al 2016b), (3) to calculate and update the dose plan and (4) the needle insertion sequence during the procedure

To make a step towards this ultimate goal, a pipeline (see figure 1(a)) was proposed by our group in a previous study (Borot de Battisti et al 2016a): this pipeline consists of calculating an initial dose plan and re-optimizing the dose plan during the intervention with feedback on sub-optimal needle positioning. This study showed that updating the dose plan, during the interventional procedure, may improve the delivered dose and allow to reach the clinical constraints. However, the needle insertion sequence is left to the doctor and may thus be potentially sub-optimal. A tool to assist the doctor in determining the optimal needle insertion sequence would therefore be of great interest. In order to develop a fully automatic control system with feedback on sub-optimal needle positioning and increase the chance of reaching the clinical constraints, the fast, automatic and adaptive determination of the optimal needle insertion sequence is therefore mandatory. In the scenario of TRUS-guided needle insertion with parallel needle configuration, the GEC/ESTRO working group (Hoskin et al 2013) recommends to start to implant with the anterior catheters. That way, issues related to the ultrasound shadowing effect behind the implant needles are reduced. Moreover, interference with the pubic arch can be checked at an early stage of the intervention so the setup can be adjusted to overcome this. However, those issues are not applicable for needle-by-needle delivery or MR-guided needle insertion.

Zoom In Zoom Out Reset image size

In this manuscript, a new needle insertion sequencing strategy for MR-guided focal HDR brachytherapy involving needle-by-needle dose delivery is proposed (figure 1(b)). The needle is inserted along the track which has the largest impact on the dose coverage if a sub-optimal needle positioning occurs. That way, potential sub-optimal needle positioning is compensated by the re-optimization of the subsequent needle insertions. The impact on the dose coverage (or also called the 'sensitivity' later on in this manuscript) of each possible needle track is predicted by a stochastic method based on needle insertion simulations. The proposed needle sequencing strategy was assessed by simulating MR-guided HDR prostate brachytherapy on 11 patients with different sequences of needle insertions.

The main contribution of this study is therefore threefold:

• (i)  
  A new adaptive needle insertion sequencing is introduced, based on the determination of the most sensitive needle track.
• (ii)  
  A new stochastic criterion to determine the sensitivity of each needle track is proposed, based on needle insertion simulations.
• (iii)  
  The performance of the proposed needle insertion sequencing is assessed by simulating MR-guided HDR prostate brachytherapy using MR-data of 11 patients diagnosed with prostate cancer.

In this manuscript, a new adaptive and automatic needle insertion sequencing strategy is introduced. The method consists of inserting the needle into the most sensitive needle track: that way, sub-optimal needle positioning may still be compensated by the re-optimization of the subsequent needle insertions. For this purpose, a criterion designed to determine automatically the sensitivity of the needle tracks is presented in section 2.3. The adaptive needle sequencing is then assessed by simulating complete brachytherapy procedures on 11 patients with varying number of needle insertions (see section 2.4).

2.1. Clinical constraints

Focal prostate HDR brachytherapy aims to deliver an optimal dose distribution with a high irradiation dose to the tumor and the lowest possible dose to surrounding OARs. At the UMCU, a dose plan is considered clinically acceptable when:

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} {{D}_{95 \%\text{PTV}}}>D_{95 \%\text{PTV}}^{\text{min}} \\ {{D}_{10 \%\text{Ur}}}<D_{10 \% Ur}^{\text{max}} \\ {{D}_{1\text{cc}~\text{Rec}}}<D_{1\text{cc}~\text{Rec}}^{\text{max}} \\ {{D}_{1\text{cc}~\text{Bl}}}<D_{1\text{cc}~\text{Bl}}^{\text{max}} \end{array}\right. \end{eqnarray} \tag{ 1 }$$

where ${{D}_{95 \% ~\text{PTV}}}$, ${{D}_{10 \% ~\text{Ur}}}$, ${{D}_{1\text{cc}~\text{Rec}}}$ and ${{D}_{1\text{cc}~\text{Bl}}}$ are the doses received by $95 \%$ of the PTV, by $10 \%$ of the urethra and by 1cc of rectum and bladder, respectively. $D_{95 \% ~\text{PTV}}^{\text{min}}$, $D_{10 \% ~\text{Ur}}^{\text{max}}$, $D_{1\text{cc}~\text{Rec}}^{\text{max}}$ and $D_{1\text{cc}~\text{Bl}}^{\text{max}}$ correspond to the minimal accepted value of ${{D}_{95 \% ~\text{PTV}}}$ and the maximal accepted value of ${{D}_{10 \% ~\text{Ur}}}$, ${{D}_{1\text{cc}~\text{Rec}}}$ and ${{D}_{1\text{cc}~\text{Bl}}}$ in order to obtain a clinically acceptable and optimal dose plan. Their values are equal to 19Gy, 21Gy, 12Gy and 12Gy, respectively. 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 19Gy to the target.

2.2. Assessment of the quality of the dose plan

To determine automatically the sensitivity of each needle track, the quality of a dose plan must be assessed. An evaluation parameter is therefore mandatory to compare the quality of different dose plans. In the scope of this study, we used the evaluation parameter E introduced by Borot de Battisti et al (2015) called the 'energy parameter' and is such that the greater E is, the better quality the dose plan becomes. E is expressed as the minimum of the four parameters, A, B, C and D where:

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} A={{D}_{95 \%\text{PTV}}}-D_{95 \%\text{PTV}}^{\text{min}} \\ B=D_{10 \%\text{Ur}}^{\text{max}}-{{D}_{10 \% ~\text{Ur}}} \\ C=D_{1\text{cc}~\text{Rec}}^{\text{max}}-{{D}_{1\text{cc}~\text{Rec}}} \\ D=D_{1\text{cc}~\text{Bl}}^{\text{max}}-{{D}_{1\text{cc}~\text{Bl}}} \end{array}\right. \end{eqnarray} \tag{ 2 }$$

The parameters A, B, C and D represent the differences between the dose coverage parameters and the clinical constraints of the planning target volume (PTV), the urethra, the rectum and the bladder, respectively. E represents therefore the smallest difference in Gray between each dose parameter and the clinical constraints. Theoretically, E can take any value between $-\infty$ and $\infty$. Moreover, if E  >  0, all dose parameters achieved the clinical constraints. Conversely, if $E\leqslant 0$, at least one dose parameter did not achieve the clinical constraints.

It is important to note that the use of this particular evaluation parameter is not mandatory to perform the needle insertion sequencing proposed in this manuscript: the needle insertion sequencing strategy is also compatible with the use of any other evaluation parameter expression which assesses the dose plan quality.

2.3. New criterion to determine the sensitivity of each needle track

In this section, a criterion is proposed to predict automatically the sensitivity of each possible needle insertion track. The needle track sensitivity is determined using a stochastic method based on needle insertion simulations. For a given needle track i ($i\in \left[1,{{N}_{\text{tracks}}}\right]$, N_tracks is the number of possible needle insertion tracks), the algorithm assessing the needle track sensitivity is described as follows:

• Step 1: Starting from the current dose plan, N_angle needle insertions are simulated by randomly modifying the needle track angulation (the dwell times and positions along the needle remain unchanged). N_angle is an integer large enough to be statistically acceptable. The random angulation is modeled by a Gaussian distribution with a standard deviation of 0.015 rad: this corresponds to a typical insertion error of 3 mm at a distance of 200 mm from the center of rotation. After investigation, a typical value of ${{N}_{\text{angle}}}=100$ was found to be a good compromise between speed and statistical significance.
• Step 2: The N_angle newly created dose plans are then evaluated by calculating the energy parameter E (defined in section 2.2) for each dose plan.
• Step 3: The 5th percentiles of the determined E-values (noted ${{E}_{5 \% }}(i)$) are calculated. ${{E}_{5 \% }}(i)$ corresponds to the worst impact on the dose that a sub-optimal needle positioning in needle track i may induce. ${{E}_{5 \% }}(i)$ is therefore related to the sensitivity of needle track i: the higher ${{E}_{5 \% }}(i)$ is, the less sensitive the needle track becomes. The 5th percentiles are preferred rather than the minimum in order to avoid outliers.

With this method, we can predict the most and least sensitive needle track ${{i}_{\text{least}~\text{sens}}}$ and ${{i}_{\text{most}~\text{sens}}}$ which are expressed as follows:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{i}_{\text{least}~\text{sens}}}=\underset{i\in \left[1;{{N}_{\text{tracks}}}\right]}{{\text{argmax}}}\,\left[{{E}_{5 \% }}(i)\right] \end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{i}_{\text{most}~\text{sens}}}=\underset{i\in \left[1;{{N}_{\text{tracks}}}\right]}{{\text{argmin}}}\,\left[{{E}_{5 \% }}(i)\right] \end{array}\end{eqnarray} \tag{ 4 }$$

The optimal needle sequence corresponds thus to inserting the needle into the needle track ${{i}_{\text{most}~\text{sens}}}$ of the current dose plan. Conversely, the least optimal sequence corresponds to inserting the needle into the needle track ${{i}_{\text{least}~\text{sens}}}$ of the current dose plan.

2.4. Experimental evaluation: simulation of HDR prostate brachytherapy to assess the proposed adaptive needle sequencing

In this section, the experiment to assess the proposed adaptive needle insertion sequencing is described. Focal HDR prostate brachytherapy procedures were simulated on 11 patients with different number of needle insertions (4, 6, 8, 10 and 12). The simulated brachytherapy procedures followed the strategy described in figure 1(b). Sub-optimal needle positioning was simulated for each needle insertion (see section 2.4.3). Two situations were compared:

• (i)  
  When the least sensitive needle track was selected for insertion: that corresponds to the 'least optimal' needle insertion sequence.
• (ii)  
  When the most sensitive needle track was selected for insertion: that corresponds to the 'optimal' needle insertion sequence.

All brachytherapy simulations were repeated 100 times to assess the robustness of the two tested needle sequences. For each scenario, the final dose parameters ${{D}_{95 \% ~\text{PTV}}}$, ${{D}_{10 \% ~\text{Ur}}}$, ${{D}_{1\text{cc}~\text{Rec}}}$ and ${{D}_{1\text{cc}~\text{Bl}}}$ and the final energy parameter E (assessing the quality of the dose distribution) were calculated.

The details of the simulations are presented as follows.

2.4.1. Dose planning algorithm.

The calculation of the initial dose plan and the re-optimization of the dose plan were performed using the algorithm described by Borot de Battisti et al (2015) and Borot de Battisti et al (2016a), respectively: the algorithm involves an exhaustive search of the center of rotation of the setup, a heuristic optimization of the needle track angulations and the determination of the dwell times by solving linear equations.

2.4.2. Anatomy.

Between May 2013 and April 2016, 30 patients were treated by focal MR-guided focal HDR prostate brachytherapy as monotherapy at the UMCU according to the current practice: catheters were first inserted manually under TRUS-guidance with the support of a template grid, and needle reconstruction and dose planning were then assessed using MRI. These patients had localized prostate cancer with a prostate-specific antigen (PSA) level lower than 10 ng ml⁻¹ and a Gleason score of 7 or less. The intra-operative MR-data—taken directly after the needles were in place—of the first 15 consecutive patients were included in the simulation study. The anatomy of the patients was obtained using the delineations of the prostate tumors and the OAR considered (urethra, bladder and rectum) on $1~\text{m}{{\text{m}}^{3}}$ resolution MR images by an experienced oncologist. 4 patients were left out of the study because the initial simulated dose plan did not reach the clinical constraints. Consequently, the brachytherapy procedures were simulated using MR-data of 11 patients. The PTV ranged from 9.5 cm³ to 31.3 cm³ with a median of 22.8 cm³. The acquisition of patient MR images used in this study was approved by the Institutional Review Board (IRB). It is noticeable that Patient 9, 10 and 11 were considered as more 'difficult cases' for dose planning compared to the other patients because the PTV was very close (if not touching) an OAR (see figure 2): patient 9 and 10 had the PTV wound around the urethra and for Patient 11, the PTV was very close to the rectum. Consequently, the arrangement of needle tracks was limited for those patients in order to avoid the needle going through an OAR.

Zoom In Zoom Out Reset image size

To clarify the benefits of the proposed needle sequencing strategy on the overall results, the anatomy of the patient (delineated before the initial dose plan) was supposed to remain constant throughout the procedure in order to exclude possible artifacts and uncertainties of image registration or dose accumulations.

2.4.3. Sub-optimal needle positioning.

Sub-optimal needle positioning was modeled as follows: each needle insertion had a random angulation error described by a Gaussian distribution with a standard deviation of 0.015 rad (i.e. 3 mm at the tip of the needle). The value of the standard deviation was chosen in line with the study of Strassmann et al (2011) who assessed the accuracy of needle positioning on a prostate model: in the worst case of the two tested scenarios (manual and robotic needle insertion), the average needle positioning accuracy was $2.7\pm 0.7$ mm. Additionally, we assumed that, during the brachytherapy intervention, the needle position could be measured without error.

2.4.4. Model of dose calculation.

To simulate the distributed dose to the patient, the point source approximation model was chosen to calculate the dose rate of the source 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}&&\begin{array}{*{35}{l}} {{d}_{i}}\left(\boldsymbol{r}\right)={{S}_{\text{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{array}\end{eqnarray} \tag{ 5 }$$

where ${{d}_{i}}\left(\boldsymbol{r}\right)$ is the dose rate of the ith source position at $\boldsymbol{r}$, 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) 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)=\,\parallel \boldsymbol{r}_{{i}}-\boldsymbol{r}{{\parallel}_{2}}$). With this adaptation of the point source model, the dose has an upper limit value close to the source, therefore it reduces 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~\text{cGy}\centerdot {{\text{h}}^{-1}}\centerdot {{\text{U}}^{-1}}$) and an arbitrary source strength ${{S}_{\text{K}}}=40.80\,\text{mGy}\centerdot \text{h}{{\,}^{-1}}\centerdot {{\text{m}}^{2}}$ 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)\centerdot {{ \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\centerdot {{10}^{-3}}$ and ${{a}_{2}}=3.12\centerdot {{10}^{-6}}$, where R is in millimeters.

The results of the brachytherapy simulations corresponding to Patient 2 with 4 needle insertions are depicted in figure 3. Figures 3(a) and (b) represent the final needle trajectory locations and the order of needle insertions resulting from 1 simulation in the situation of optimal and least optimal needle insertion sequence respectively. Figure 3(c) presents the final dose parameters resulting from the 100 simulations: the final dose parameters [${{D}_{95 \% ~\text{PTV}}}$, ${{D}_{10 \% ~\text{Ur}}}$, ${{D}_{1\text{cc}~\text{Rec}}}$, ${{D}_{1\text{cc}~\text{Bl}}}$] in the situation without sub-optimal needle positioning (calculated from the initial dose plan) were [20.0Gy, 19.8Gy, 8.6Gy, 10.9Gy]. With sub-optimal needle positioning, the median(minimum,maximum) of the final dose parameters was [19.2(16.4, 20.2)Gy, 20.5(19.0, 23.6)Gy, 8.2(6.9,12.0)Gy, 11.2(10.4,12.9)Gy] for the least optimal needle sequencing scenario. Using the proposed optimal needle sequencing, these values were equal to [19.4(16.0, 20.6)Gy, 20.3(18.6, 23.7)Gy, 7.9(6.5, 11.0)Gy, 10.7(9.5, 11.8)Gy].

Zoom In Zoom Out Reset image size

The final energy parameters corresponding to Patient 2 with 4, 6, 8, 10 and 12 needle insertions are presented in figure 4 in the situation of optimal and least optimal needle insertion sequence. For (4, 6, 8, 10, 12] needle insertions, the median(minimum,maximum) of the energy parameter was [0.2(−2.6, 1.2), 0.8(−1.8, 1.6), 1.5(0.5, 2.1), 1.7(−0.1, 2.3), 1.8(0.8, 2.4)) in the situation of least optimal sequencing and [0.4(−3.0, 1.6)Gy, 1.3(−1.7, 2.0)Gy, 1.6(0.2, 2.1)Gy, 1.8(0.8, 2.3)Gy, 2.0(1.0, 2.5)Gy] in the situation of optimal sequencing.

Zoom In Zoom Out Reset image size

Table 1 presents the percentage of final dose distributions which fulfilled all clinical constraints (i.e. $E\geqslant 0$) for the two tested needle sequencing protocols. According to the table, the percentage of clinical acceptable final dose distributions was equal or higher for 50 scenarios tested with the optimal needle sequence compared to the least optimal needle sequence and lower for 5 scenarios (Patient 5 with 4 needle insertions, Patient 8 with 8 needle insertions, Patient 9 with 6 needle insertions, and Patient 11 with 8 and 10 needle insertions). This corresponds to a dose improvement in 91% of the scenarios tested. It is noticeable that the percentage of clinically acceptable dose plans was low for Patient 9, 10 and 11 (considered as the 'difficult cases') compared to other patients in both situations of least optimal and optimal needle sequencing.

In this study, a new automatic and adaptive needle sequencing strategy was proposed for MR-guided HDR prostate brachytherapy. This sequencing is based on the determination of the most sensitive needle track. For that, a new criterion was introduced to predict and quantify the sensitivity of each possible needle insertion track.

The simulation study showed that inserting the needle into the most sensitive needle track will statistically lead to an improved dose coverage compared to inserting into the least sensitive track. Most specifically, the final dose parameters presented in figure 3 illustrate this trend: when the needle insertion sequence was optimal (i.e. the needle is inserted into the most sensitive needle track), the dose tended to increase at the PTV and to decrease at the OARs in comparison to the situation of least optimal needle insertion sequence (i.e. the needle is inserted into the least sensitive needle track). Inserting the needle into the most sensitive needle track allows sub-optimal needle positioning to be compensated by the re-optimization of the subsequent needle insertions. Moreover, for some patients, the difference in the dose coverage was more pronounced for small number of needle insertions (see figure 4). This could be explained by the fact that the dose delivered by each needle decreases with the number of needle insertions: for a small number of needle insertions, the dose delivered by each needle is high and there are few needles to compensate for sub-optimal positioning. Therefore, the robustness of the dose plan to sub-optimal needle positioning decreases when a small number of insertions is used. Consequently, the difference of impact on the distributed dose due to sub-optimal needle positioning between least optimal and optimal needle sequence situation is more important with a small number of needle insertions. This effect is attenuated for a large number of needle insertions. In this study, there was no clear correlation between the percentage of clinically acceptable dose plans and the prostate volume. Although a large prostate (or PTV) may be more difficult to cover compared to a small one, the accessibility of the target by the needle (which may be impaired by the urethra or the pubic bone), or the localization of the PTV with respect to the OARs may also play a crucial role to obtain a clinical acceptable dose plan. Another important point concerns Patient 9, 10 and 11: for those patients, the percentage of clinically acceptable dose plans was lower compared to the other patients (see table 1). For those patients, the PTV was very close to an OAR. Because of this, (1) the arrangement of needle tracks was limited and (2) the dose coverage of the PTV was restricted by the neighboring OAR. Therefore, those patients had less chance to receive a clinically acceptable dose. Finally, the brachytherapy simulations were implemented in MatLab 2015: the computation time increased linearly with the number of needle tracks involved in the sequencing process, and was less than 6 seconds per needle track.

In the presented study, an open number of needle insertions was chosen. In practice, a compromise has to be found: it is clear that the more needle insertions there are, the better the dose distribution will be. However, the influence of adding one or more needle insertions may lead to additional toxicity: Vargas et al (2005) and Boyea et al (2007) showed the urinary toxicity following HDR brachytherapy is significantly increased by using more than 14 needle insertions. For those reasons, we believe that an automatic determination of the number of needle insertions requires additional work and we decided to perform the simulations for different number of needle insertions. A range of 4 to 12 needle insertions was chosen because, in this range, clinically acceptable dose plans in most patient cases were obtained. It is also noticeable that Patient 11 is a difficult case because the percentage clinically acceptable dose plans obtained was less compared to the other patients. The implication for this patient may be to increase the number of needle insertions until the clinical constraints are reached.

The limitations of this study are mostly related to the limitations in the proposed brachytherapy simulation. Firstly, the point source approximation employed in this study (see equation (5)) was chosen for implementation simplicity as well as computational speed considerations. However, the proposed needle sequence determination allows more precise source models such as line source approximation. Moreover, the anatomy was supposed to remain constant throughout the procedure in order to exclude possible artifacts and uncertainties of image registration or dose accumulations. In practice, the intra-operative internal organ motion due to, for example, swelling, displacement (Stone et al 2002) or rotation (Lagerburg et al 2005) of the prostate (related to the trauma caused by the needle insertion) or intra-procedural changes in rectum or bladder filling can cause uncertainties in the delivered dose during the intervention. Those issues can be solved by performing an adaptation of the dose plan and needle insertion sequence with feedback on the anatomy changes (in addition to feedback on sub-optimal needle positioning). The dose plan and needle sequencing adaptation with feedback on anatomy movements is possible while performing brachytherapy treatment under MRI-guidance (but also US-, CT-guidance) for which the anatomy can be captured during the interventional procedure. Finally, the MR data used for the simulation study were taken after insertion of the catheters. Since, one of the factors that affect needle insertion is prostate swelling, an interesting investigation can be to compare the magnitude of dose difference if the simulations are performed on images before versus after implant. This will be investigated in future works.

The error associated with the measurement of the needle position should also lead to uncertainties on the delivered dose. This issue may be tackled by determining with a high accuracy and update rate the needle position in order to: (1) track the needle during insertion and consequently help the steering and the placement of the needle in the correct position, (2) reconstruct the needle after insertion to re-optimize the dose plan with feedback on sub-optimal needle positioning. To insure the fast and accurate measurement of the needle position, we believe fiber Bragg gratings (FBG)-based needle tracking described by Borot de Battisti et al (2016b) may be a great candidate. FBG-based needle tracking involves a stylet that can be inserted into brachytherapy needles to measure quickly and accurately the needle shape. The error propagation due to FBG-based tracking, adaptive dose planning and needle sequencing combination will be studied in future work.

The proposed needle sequencing is compatible with the experimental setup developed at the UMCU for MR-guided HDR prostate brachytherapy:

• in terms of hardware: the robotic device is designed to support needle insertion under MR-guidance
• in terms of software: the computation time of the needle insertion sequencing strategy is eligible for intra-operative use

Furthermore, the proposed needle sequencing may also be used (to some extent) for other applications with different setup as follows:

• using parallel and divergent needle patterns
• with or without brachytherapy template
• involving needle-by-needle delivery or delivery after all needles are in place
• with different imaging modalities such as CT, TRUS or MRI

since the stochastic approach only relies on the knowledge of the PTV and OARs localization.

In the presented study a new automatic and adaptive needle insertion sequencing was proposed for focal MR-guided HDR prostate brachytherapy involving needle-by-needle dose delivery. The approach consists of inserting the needle into the most sensitive needle track. To predict the sensitivity of each possible needle track, a stochastic criterion was proposed based on needle insertion simulations. To assess this needle sequencing strategy, HDR prostate brachytherapy was simulated on 11 patients with varying number of needle insertions (from 4 to 12). An improvement in the distributed dose was observed in 91% of the tested scenarios for which the needle was inserted into the most sensitive compared to the least sensitive needle track. Finally, the computation time for sequencing was less than 6 seconds per needle track (implementation on MatLab 2015). This novel adaptive sequencing tool can therefore assist the doctor during the intervention. Furthermore, it is a step towards the development of a fully automatic control system with feedback on unpredictable events occurring during brachytherapy, such as sub-optimal needle positioning and/or intra-operative internal organ motion.

This study was funded by Philips Medical Systems Nederland B.V. M Borot de Battisti is funded by Philips Medical Systems Nederland B.V. M A Moerland is principal investigator on a contract funded by Philips Medical Systems Nederland B.V. G Hautvast and D Binnekamp are full-time employees of Philips Medical Systems Nederland B.V. The authors also thank the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and the ITEA (project 12026, SoRTS).