Cornel Zachiu

Baudouin Denis De Senneville

Chrit Moonen

Mario Ries

A framework for the correction of slow physiological drifts during MR-guided HIFU therapies: Proof of concept

Keywords: motion correction, MR-guidance, non-invasiveness, HIFU, slow physiological drifts

A spatial coherence of all three is particularly important to assure a confluent ablation volume and to prevent remaining islets of viable malignant tissue.

Conclusions:

This study proposes a motion correction strategy for displacements resulting from slowly varying physiological motion that might occur during an MR-guided HIFU intervention. We have shown that such drifts can lead to a misalignment between interventional planning, energy delivery and therapeutic validation. The presented volunteer study and in-vivo experiment demonstrates both the relevance of the problem for HIFU therapies and the compatibility of the proposed motion compensation framework with the work-flow of a HIFU intervention under clinical conditions.

Purpose: While respiratory motion compensation for MR-guided HIFU interventions has been extensively studied, the influence of slow physiological motion due to, for example, peristaltic activity, has so far been largely neglected. During lengthy interventions, the magnitude of the latter can exceed acceptable therapeutic margins. The goal of the present study is to exploit the episodic work-flow of these therapies to implement a motion correction strategy for slowly varying drifts of the target area and organs at risk over the entire duration of the intervention.

Methods:

The therapeutic work-flow of an MR-guided HIFU intervention is in practice often episodic: Bursts of energy delivery are interleaved with periods of inactivity, allowing the effects of the beam on healthy tissues to recede and/or during which the plan of the intervention is re-optimized. These periods usually last for at least several minutes. It is at this time scale that organ drifts due to slow physiological motion become significant. In order to capture these drifts we propose the integration of 3D MR scans in the therapy work-flow during the inactivity intervals. Displacements were estimated using an optical flow algorithm applied on the 3D acquired images. A preliminary study was conducted on 10 healthy volunteers. For each volunteer 3D MR images of the abdomen were acquired at regular intervals of 10min over a total duration of 80min. Motion analysis was restricted to the liver and kidneys. For validating the compatibility of the proposed motion correction strategy with the work-flow of an MR-guided HIFU therapy, an in-vivo experiment on a porcine liver was conducted. A volumetric HIFU ablation was completed over a time span of 2h. A 3D image was acquired before the first sonication, as well as after each sonication.

Results: Following the volunteer study, drifts larger than 8mm for the liver and 5mm for the kidneys prove that slow physiological motion can exceed acceptable therapeutic margins. In the animal experiment motion tracking revealed an initial shift of up to 4mm during the first 10min and a subsequent continuous shift of ∼2mm/h until the end of the intervention. This leads to a continuously increasing mismatch of the initial shot planning, the thermal dose measurements and the true underlying anatomy. The estimated displacements allowed correcting the planned sonication cell cluster positions to the true target position, as well as the thermal dose estimates during the entire intervention and to correct the NPV-measurement.

I. INTRODUCTION

Magnetic resonance guided high intensity focused ultrasound (MRg-HIFU) is an attractive solution for treating tumors inside the human body in a non-invasive manner [START_REF] Cline | MR-guided facused ultrasound surgery[END_REF][START_REF] Hynynen | A clinical, noninvasive, MR imaging-monitored ultrasound surgery method[END_REF] . Such interventions have the potential to target and ablate non-invasively primary and metastatic cancer deep in the human body in a single session [START_REF] Illing | The safety and feasibility of extracorporeal high-intensity focused ultrasound (hifu) for the treatment of liver and kidney tumours in a western population[END_REF] . Due to the associated reduced burden for the patient, treatment related patient morbidity and complication rate, HIFU represents a potential alternative to patient groups which are currently not eligible for invasive or minimally invasive therapy [START_REF] Tempany | Focused ultrasound surgery in oncology: overview and principles[END_REF][START_REF] Foley | Image-guided focused ultrasound: state of the technology and the challanges that lie ahead[END_REF] . MRg-HIFU has been successfully applied for a number conditions such as symptomatic uterine fibroids, prostate cancer, liver tumors, thyroid nodules, palliative treatment for bone metastases, breast fibroadenoma and functional neurological disorders, including essential tremor, neuropathic pain and tremor-dominant Parkinson's disease [START_REF] Foley | Image-guided focused ultrasound: state of the technology and the challanges that lie ahead[END_REF] . However, MRg-HIFU in clinical applications becomes challenging when the tumor is located in the abdomen or the lower thorax due to continuous physiological motion [START_REF] Ries | Realtime 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF][START_REF] Crijns | Proof of concept of MRI-guided tracked radiation delivery: tracking one-dimensional motion[END_REF][START_REF] Diez | Analysis and evaluation of periodic physiological organ motion in radiotherapy treatments[END_REF][START_REF] Goitein | Organ and tumor motion: An overview[END_REF][START_REF] Vijayan | Validation of a non-rigid registration method for motion compensation in 4D ultrasound of the liver[END_REF] .

Uncertainties regarding the exact tumor position might lead to undesired thermal damage to healthy tissue or the delivery of a thermal dose to the tumor that is below the lethal threshold [START_REF] Ries | Realtime 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF] . In addition to shifts caused by respiration, tumors in the abdomen and the lower thorax will also suffer displacements due to digestive and metabolic processes, abdominal and thoracic muscle relaxation and ultimately with the spontaneous motions that the patient might exhibit. Important is to realize that these types of motion do, in general, not occur on the same time-scale:

• Respiratory induced motion of free-breathing patients has a frequency of 0.2-0.4 Hz and displays a rapid in and exhalation of a duration of 1-1.5 s followed by a rest of the diaphragm of 2-3 times of this duration. The majority of the existing scientific work on motion correction schemes for external beam-interventions focuses on the correction of this type of motion 11,12 .

• Although the time-scale of digestive, metabolic or muscle relaxation motion events depends on the event source (e.g. development of digestive gases or the passage of digestive products in the gastrointestinal tract, bladder filling with urine, changes in patient respiration), the resulting shifts in the position of abdominal organs has in most studies been identified to become significant on a time-scale of several minutes [13][14][15] .

For this reason this type of motion will be referred to in this paper as slow physiological motion or drifts. Digestive activity can clinically be moderated by several measures: Peristaltic bowel motion and the development of digestive gases can be reduced by adjusting the diet of the patient prior to the intervention 16 , or by the administration of butylscopolamine and/or glucagon which act as digestive motility inhibiting agents 17,18 [REF for glucagon]. Furthermore, shifts in the lower abdomen due to bladder filling can be reduced by the use of Foley catheters 19,20 .

• Finally, spontaneous motion is a difficult problem -it occurs infrequently and on a very short time-scale. It is particularly problematic for long interventions, which require the patient to remain in an uncomfortable position. In the past, this problem has been addressed by either using restraints such as molds or casts 21 or by sedating the patient 22,23 .

Although a motion compensation technique is required in order to reduce tumor displacement effects on the outcome of the therapy, it should imply minimal discomfort caused to the patient. For respiratory and slow physiological motion this can be achieved through the means of an image guided intervention [START_REF] Ries | Realtime 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF][START_REF] Crijns | Proof of concept of MRI-guided tracked radiation delivery: tracking one-dimensional motion[END_REF][24][25][26] . Beneficial for the development of a motion compensation strategy for HIFU therapies is the fact that the therapeutic work-flow is in practice often episodic: Bursts of energy delivery are interleaved with periods of inactivity. More specifically, interventions usually alternate between short bursts of 30-60s of ultrasonic energy and replanning/cool-down phases of several minutes [START_REF] Foley | Image-guided focused ultrasound: state of the technology and the challanges that lie ahead[END_REF]27 . During energy delivery respiration is the main source of abdominal organ motion, thus the focus will be on compensating for this type of displacements. Given the larger time scale at which organs shift due to slow physiological drifts, this type of motion can be compensated during the replanning/cool-down phases.

While substantial work has already been conducted for the compensation of respiratory motion, this study aims to exploit the episodic work-flow of an MRg-HIFU intervention to propose a motion correction scheme for slow physiological motion. The study is divided in three parts: First, we assess long term shifts of the liver and kidney position in ten healthy volunteers over a duration of 80min. This is achieved by the periodical examination of the area or the upper abdomen with fast 3D MRI and an on-the-fly 3D co-registration of the image volumes, which are both compatible with the therapeutic work-flow of an MRg-HIFU therapy. Subsequently, we assess a lower boundary of the accuracy of the proposed co-registration methods for both organs based on a comparison with manually tracked anatomical landmarks at selected areas. Finally, we demonstrate the clinical benefit of the proposed approach by correcting for the slowly varying organ drifts during an in-vivo high intensity focused ultrasound ablation on a porcine liver.

A similar study was recently conducted by Schmitt et al. [REF] in which temperature and thermal dose measurements were corrected with respect to susceptibility artifacts arising from bowel motion during MRg-HIFU interventions in the prostate. While we propose a correction strategy based on displacements provided by a motion estimation algorithm, their study focuses on developing and implementing a filter that excludes temperature information in regions that manifest a higher temperature fluctuation than the heated area. The approach was validated in five out of six clinical cases showing good potential for further use not only in the prostate but also for targets situated higher in the abdomen.

II. MATERIALS AND METHODS

II.A. General strategy

II.A.1. Proposed MR imaging framework

Figure 1 details the HIFU ablation protocol employed in our study. Our framework consists of performing a series of n ablative sessions, each one having an approximate duration of 5 -10 minutes (which is divided between the sonication and a cool-down interval). Each sonication is usually preceded by a planning/re-planning phase, which consists in establishing/refining the position and size of the area that needs to be ablated by that respective sonication. This step is typically performed on a separate high-resolution T1/T2 weighted image acquired at the beginning of the therapy. Note that the MR acquisition sequence for the planning image might differ from the one of the subsequent 3D T1-weighted abdominal scans (see Figure 1). During the actual sonications, MR thermometry is performed in order to monitor the temperature of the sonicated tissue and detect when the lethal thermal dose has been delivered. In this study we were interested in tracking the liver and kidneys over the entire duration of the intervention. At this timescale the organs may exhibit significant displacements in all three directions. Thus, organ tracking should be performed in 3D. With current technology this becomes, however, a challenging task during sonications due to the real-time feedback loop required for temperature monitoring [START_REF] Ries | Realtime 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF] . Such a time constraint does currently not allow the acquisition of 3D images with reasonable quality coupled with a 3D image based tracking procedure. For this reason imaging and tracking during sonications has to be performed in 2D. During this time interval the focus is mainly on compensating for displacements caused by respiratory motion. In order to extend liver and kidney tracking to the duration of the entire therapy, we propose interlacing 3D T1-weighted abdominal scans in between sonications. These scans are evenly distributed over the time span of the therapy with a time gap T a in between them. Registration of the acquired 3D volumes to a reference position allows compensation of the 3D displacements exhibited during sonication by the targeted organs due to slow physiological motion. The acquisition and registration of the 3D volumes in between sonications is timewise feasible. This is due to the fact that each sonication is followed by a time interval within which the tissue situated in the near and far field of the HIFU transducer is allowed to cool down. The displacements estimated during the 3D registration process are then used to locate the ablated tissue on a 3D reference volume.

II.B. Estimation of 3D organ deformation

The 3D organ deformations were estimated after each sonication using an optical flow algorithm applied on the T1-weighted volumetric scans. The deformation estimation is based on the optical flow algorithm proposed by Horn&Schunck 28 . In its original form the algorithm is unable to estimate displacements that are larger than the size of one voxel. In order to overcome this limitation, a multi-resolution approach with an iterative refinement of the optical flow has been employed 29 . A multi-threaded C++ implementation of this algorithm has been provided by Team MR-HIFU, UMC Utrecht, The Netherlands. The algorithm requires as input two images: a reference image and an image to register. For the volunteer study and also for the animal experiment, the first image in their corresponding series of 3D T1-weighted abdominal scans was established as reference. Additionally, the optical flow algorithm requires fixating a parameter α which reflects the elasticity of the tracked organs. All 3D images were registered using α = 0.05, value which was established empirically.

Given that the optical flow algorithm is gray-level intensity based 28 , it is important to evaluate its robustness to the noise that might superimpose on the images during the acquisition process. For this purpose, registration has been performed between two un-altered 3D T1-weighted abdominal scans and the resulting motion vector field was established as gold standard. Afterwards, noise was added to the image to register, gradually reducing its signal-to-noise ratio (SNR). Registration was then performed for each noise level and a mean square error(MSE) between the current estimated motion vector field and the gold standard was computed. Such an approach will aid in detecting the lowest value of the SNR for which the estimated motion is still reliable.

II.C. Experimental setup

All MR imaging was performed on a 1.5T Philips Achieva scanner (Philips Healthcare, Best, The Netherlands).

II.C.1. Healthy volunteer study

The slowly varying drifts exhibited by the human liver and kidneys have been assessed for ten healthy volunteers. All volunteers went inside the MR scanner in a head firstprone position and were instructed not to move over the duration of the image acquisition process. A series of nine 3D T1-weighted scans were performed on the abdomen of each volunteer. Timewise the scans were evenly distributed over a duration of approximately T T herapy = 80 minutes with an average time gap of T a = 10 minutes in between the scans.

The MR acquisition sequence was a fat suppressed spoiled gradient recalled echo employing the following parameters: TE = 2 ms, TR = 4.3 ms, image matrix 192 × 192 × 75, 10 • flip angle, bandwidth in phase-encoding direction 373.4 Hz, an isotropic voxel size of 2 × 2 × 2 mm 3 , resulting an acquisition time of ∼ 90 s. Fat suppression was achieved through the spectral presaturation with inversion recovery technique. For the remainder of this paper the images acquired using the MR sequence with the aforementioned parameters will be referred to as 3D anatomical images. An exception to this nomenclature will be made for the 3D anatomical images that play the role of reference in the registration process. These images will be refered to as 3D anatomical anchors.

II.C.2. Animal experiment

HIFU ablation has been performed in-vivo on a porcine liver using a Sonalleve (Philips Healthcare, Helsinki, Finland) ablation system. The pig was anesthetized by an initial intramuscular injection containing ketamine (13 mg/kg), midazolam (0.7 mg/kg), atropine (0.02 mg/kg) and meloxicam (0.4 mg/kg). Subsequently, anesthesia was maintained by intravenous administration of sufentanil (11.3 µg/kg/h), midazolam (1 mg/kg/h) and nimbex (0.09 mg/kg/h). A total number of n = 10 sonications were carried out over an approximate duration of T T herapy = 2 hours. After each sonication a 3D anatomical image has been acquired. An additional 3D anatomical image was acquired before the first sonication, which will play the role of reference during the registration process. Closely after the acquisition of the last 3D anatomical image, a gadolinium based contrast agent has been administered to the animal. Subsequently a contrast enhanced T1-weighted image was acquired. In this manner, a visual feedback of the non-perfused tissue resulting from the 10 sonications is obtained. After the completion of the experiment, the animal was euthanized by intravenous injection of pentobarbital.

The animal experiment was conducted in agreement with the European law on animal experimentation and in compliance with the institution's rules for animal care and use.

II.D. Analysis and validation of the 3D estimated motion

II.D.1. Motion analysis and validation on healthy volunteers

In this study we were interested exclusively in the motion exhibited by the liver and kidneys. In order to restrict motion analysis to the mentioned organs of interest, a semiautomatic segmentation procedure was performed prior to the analysis. An initial segmentation was obtained by the means of active contours which was subsequently refined through a manual correction. The software tool used for segmentation was ITK-Snap v 2.4.0.

The validity of the estimated motion was tested by manually tracking blood vessel bifurcations in the liver and kidneys. For each of the volunteers three vessel bifurcations were tracked: one in the upper part of the liver (lobes 7 and 8), one in the lower part of the liver (lobes 5 and 6) and one in the kidneys. The manually tracked landmarks were then compared to the displacements estimated by the optical flow algorithm.

The manually measured displacements were also used to quantify the positional errors caused by the slowly varying drifts in the absence and respectively in the presence of a motion compensation strategy. This was achieved by comparing the displacements obtained through manual tracking (which correspond to the positional errors we would have in the absence of a motion compensation scheme) with the difference between these displacements and the displacements estimated by the optical flow algorithm (which corresponds to the residual positional errors when a motion compensation scheme is integrated in the therapy workflow). Note that for computing the positional errors, the magnitude of the 3D displacement vectors was used instead of separately evaluating the errors on each of the three directions of motion.

II.D.2. Motion analysis and validation for a real HIFU therapy on an animal

As a result of slow physiological motion, the porcine liver might suffer displacements during the HIFU ablation session. By using the acquired 3D anatomical images, the position of the sonicated areas was retrospectively tracked over the time span of the session. This allowed identifying on the planning image, acquired at the beginning of the HIFU intervention, the actual anatomy that was ablated, which, due to the organ drifts, might differ from the one originally planned. The estimated motion was also used to bring the temperature maps acquired during the 10 sonications to a common reference space determined by the planning image. In this reference space, the delivered thermal dose 30 was re-calculated. This resulted in a correction such that the newly calculated thermal dose will match the anatomy captured by the planning image. Additionally, the position of the non-perfused volume (NPV) resulting from the intervention was also brought in the reference space of the planning image.

Similar to the healthy volunteer study, motion analysis for the animal experiment has been restricted to the sonicated organ. For this purpose a segmentation procedure as the one described in section II.D.1 has been conducted.

III. RESULTS

III.A. Measured displacement for liver and kidneys caused by slow physiological drifts

III.B. Robustness to noise of the optical flow algorithm

Figure 3 displays the relationship between the SNR of the image to register and the MSE of the corresponding estimated vector field with respect to the gold standard. It can be observed that for SNR values less than 5 errors start increasing exponentially. Thus, for 285 images with an SNR value bellow this threshold, the estimated motion vectors are to be considered unreliable.

III.C. Quantitative analysis of the motion amplitude

For each volunteer the output of the registration process was a set of eight motion vector fields, each corresponding to the displacements estimated at the acquisition time of their 290 associated 3D anatomical image. The motion fields from the individual volunteers were The values on the abscissa correspond to the time delay in minutes since the acquisition of the reference image. The values on the ordinate indicate the estimated signed displacements in millimeters. Each boxplot is constructed in the following fashion: the lower limit of the box corresponds to the first quartile, the upper limit of the box is the third quartile, the median is represented by the red horizontal line in the middle of the box, the red cross indicates the mean of the set, the lower whisker corresponds to the 5 th percentile and the upper whisker to the 95 th percentile. For both the organs of interest, the overall median of the displacement in the CC direction has a tendency to increase over time reaching ∼ 3 mm.

The maximum liver displacement in the CC direction is ∼ 8 mm, while for the kidneys this measures ∼ 5 mm. For the LR and AP directions, the median of the displacements does not manifest any tendency. For the liver the maximum displacement in the LR direction is ∼ 3 mm, while for the kidneys it reaches at a value of ∼ 4 mm. In the AP direction the maximum displacement measures ∼ 2.5 mm for the liver and ∼ 4 mm for the kidneys. Figures 4(g) and 4(h) depict the time evolution of the median magnitude of the motion vectors estimated in the liver and respectively kidneys of each volunteer. These graphs reconfirm the tendency of the displacements to increase over time, however, this might happen in non-monotone manner, with the organs randomly moving closer or further from their reference position, and with rather large inter-individual variations.

III.D. Validation of the 3D estimated motion

Validation of the estimated motion was made using the approach described in section II.D.1. Figure 5 displays, for all of the volunteers, the correspondence between the manually determined landmark displacements and the landmark displacements estimated by the optical flow algorithm. The assessment was made in the cranio-caudal (blue circles), left-right (red squares) and anterior-posterior (black crosses) directions respectively. In the case of a perfect correspondence between the manually tracked landmarks and the estimated landmark positions, the points from the graph would lie on the first bisector (depicted by the purple line in figure 5). The precision of the manually tracked displacements is however limited by the discrete nature of the acquired images. Therefore tracking errors of up to one voxel (2 mm) are to be expected. Nevertheless, a paired T-test applied on the displacements obtained through the two techniques indicates that there is no significant difference between the values at 5 % significance level.

Figure 6 displays positional errors in millimeters that would affect the intervention in the selected landmarks, in the presence and respectively in the absence of a motion correction scheme. The errors were computed as described in section II.D.1. The boxplots are built in the same manner as the ones illustrated in Figure 4. More than 75% of the errors become less than 2 mm (one voxel) when motion compensation is present. The median of the errors also drops from ∼ 2.8 mm to ∼ 1.3 mm. A magnification of the planning image with focus on the originally planned cell cluster circles. A positional discrepancy between the cell cluster and the NPV can be observed.

III.E. In-vivo experiment on a porcine liver

By using the displacements estimated between the last 3D anatomical image and the 3D anatomical anchor, the contrast enhanced image was registered to match the planning image anatomy. The motion corrected NPV was then extracted and superimposed on the planning image as a high intensity region as depicted by figure 9(c). Note that due to observer subjectivity during the manual NPV delineation, the non-corrected and the corrected NPV masks might differ. Without motion compensation, the distance between the centers of mass of the position corrected cell cluster and the NPV amounts to ∼ 3.5 mm. By correcting the NPV position the distance is reduced to ∼ 1 mm.

IV. DISCUSSION

The goal of the presented study was to exploit the episodic work-flow of an MR-guided HIFU therapy to implement a motion correction strategy for the displacement of abdominal organs under the effect of slow physiological drifts. The strategy involves tracking both the target area and organs at risk (OAR) over the entire duration of the intervention.

In the first part of the study we investigated the slowly varying drifts exhibited by the human liver and kidneys during a time period of 80 minutes in 10 healthy volunteers. The obtained displacements of up to 8 mm for the liver and 5 mm for the kidneys are in good correspon- dence with previous reportings 13,26 . Please note that respiratory motion was compensated by a navigator gated acquisition scheme for all acquired images. Therefore the obtained shifts are also likely to include shifts of the diaphragm with respect to the navigator.

The findings that result from the volunteer study have shown that large inter-individual variations have to be expected. For the liver of volunteer #9 the median of the displacements remained close to 1 mm over the entire 80 minute duration of the study, while for volunteer #3 the median was already close to 4 mm halfway through the experiment (see Figure 4(h)). The measured displacements are also expected to vary with patient position: prone or supine. Thus, this reconfirms that a continuous supervision and correction of the displacements is clearly beneficial for the therapeutic success. For example for the noninvasive treatment of stage 1 tumors in the liver and kidneys the reported shifts represent ∼ 20-50% of the typical tumor dimensions. If the therapy is conducted without a motion compensation scheme, undesired shifts of this magnitude need to be included in the surgical margin and to a substantially longer intervention 31,32 .

One of the challenges that might occur while evaluating the performance of motion estimation algorithms dedicated to non-invasive in-vivo beam interventions in general is the lack of a reliable gold standard [REF]. For minimally invasive beam therapies this problem was addressed by tracking radio-emitters implanted in/or in the proximity of the targeted region [REF]. Such an approach is, however, not feasible for the applications discussed in this study since it would imply a certain degree of invasiveness. In order to validate the motion esti-mated by the employed tracking procedure we have manually tracked a group of landmarks and established the resulting displacements as a gold standard. More than 75% of the errors between the estimated displacements and the gold standard remained under 2 mm which is actually the image resolution and thus the precision limit for manual tracking. The source of these errors is the combined effect of misregistration by the motion estimation algorithm and the subjective nature of the landmark manual tracking procedure. It is also expected that as long as the assumptions made by the employed optical flow algorithm are not violated, including voxel gray-level intensity conservation from the reference to the moving image and the continuous nature of the displacements, the same results can be extrapolated to the organs as a whole. Obtaining a gold standard by the proposed approach is however limited by the presence and/or visibility of specific anatomical landmarks (e.g. vessel bifurcations) in the images considered for manual tracking. If this condition is not met, alternative methods need to be considered. Nevertheless, our findings are in good correspondence with previous reportings [34][35][36] .

The best way to improve the accuracy and precision of the estimated motion is to increase the resolution of the 3D anatomical images used for registration. This however becomes challenging due to the time constraints imposed by the therapy work-flow. If image resolution is increased while keeping the acquisition time at acceptable limits, the SNR of the image might drop to values that deteriorate the estimated motion to the point that the improvement in image resolution becomes practically useless. Accuracy of the estimated motion is also limited by the image distortions intrinsic to the employed MR acquisition sequence. The distortions mainly originate from the inhomogeneities of the MR scanner's static magnetic field and from chemical shifts (ex. water-fat shift) that might occur at organ and body borders. For example, Dragonu et al 37 have reported spatial distortions due to B 0 inhomogeneities larger than 1 cm in the liver. The study has been performed, however, for rapid gradient recalled MRI sequences using an echo-planar readout, which is not the case in the current work (see section II.C.1). Nevertheless, the results in their study indicate that B 0 inhomogeneities of ∼ 1.3 ppm are to be expected in the human abdomen. This would cause the sequence used in our study to generate spatial distortions of ∼ 0.3 mm. Such a value is well beneath the image resolution (which is 2 mm), thus its effects on the estimated motion were deemed to be negligible. Nevertheless, spatial distortions due to B 0 inhomogeneities must be addressed for MR images used for motion estimation, either by optimizing the MR acquisition sequence or by choosing an appropriate correction scheme [38][39][40][41] . In order to eliminate the effects of the water-fat shift, the MR-sequence used in our study has been optimized to be fat suppressing. This will not only remove fat shifts (and thus image artifacts) at the organ and image boundaries but it will also increase the contrast-to-noise ratio of the entire image, which is beneficial for the accuracy of the estimated motion. An additional factor that affects the quality of the estimated motion is the SNR of the images themselves, fact which is to be expected given that the optical flow algorithm employed for registration is gray-level intensity based. According to Figure 3 as the SNR drops bellow 5, the motion estimation errors start increasing exponentially. Thus, it can be considered that beyond this threshold the estimated motion vector field becomes unreliable. This however does not affect in any way the success rate of our method since the 3D anatomical images used for registration have a SNR value of ∼ 44, therefore in the right extrema of the curve depicted by Figure 3.

The Horn&Schunck motion estimation algorithm employed in our study was chosen mainly due to its previous validation for the in-vivo correction of respiratory motion during HIFU ablations in the liver [START_REF] Ries | Realtime 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]25 . There are, however, a few aspects concerning the algorithm that could be further improved. For example, since it relies on the assumption that voxels conserve their gray-level intensity on the followed trajectory, the Horn&Schunck algorithm is sensitive to gray-level intensity variations not attributed to motion. Such a situation might occur when, for instance, anatomical structures appear and/or disappear from the field-ofview. In addition, the algorithm imposes that the estimated displacements are spatially continuous. Since abdominal organs are elastic and incompressible, this is a justified assumption. However, an exception from this rule arises at the interface between the organs and the abdominal and/or the thoracic wall. In these regions shearing effects will occur, thus violating the continuity assumption. In order to overcome these limitations, versions of the optical flow algorithm that reduce the impact of such effects could be, for example, employed. Such algorithms have already been applied successfully in the video domain [REF],

however, to our knowledge, the validity of such approaches for in-vivo medical imaging ap-

plications has yet to be proven.

The focus of this study was on the analysis of the displacements exhibited by the human liver and kidneys due to slow physiological drifts, over the entire duration of a HIFU therapy session. Thus, the MR sequence used to acquire the 3D anatomical images, upon which the tracking procedure is based, was optimized to cover the entire upper abdomen. During a real HIFU intervention only the displacements experienced by the tumor itself are of interest. For this purpose the MR sequence can be easily adjusted to the target area and the surrounding organs at risk (OAR) reducing the required acquisition time from ∼ 90 s to ∼ 30 s.

The HIFU experiments performed in-vivo on a porcine liver have shown that the method proposed in our study can be applied for a real HIFU ablation therapy. In the investigated case, the animal was not fed 24 hours prior to the intervention. During the intervention, the animal was under full anesthesia and thus the thoracic and abdominal muscles were under temporary paralysis. Similar to the volunteer experiments, all respiratory induced displacements were compensated by a gating of both all MR-acquisitions and the therapeutic energy delivery. As a consequence, the observed displacements can be largely attributed to shifts due to digestive gas development in the stomach and the intestines. As the timecourse displays (see Figure 7), the majority of the displacement occurred during the first 10 minutes of the intervention. This demonstrates that even under ideal conditions, the occurrence of shifts of 3.5 -4 mm can occur and are problematic for therapy. The interventional planning was carried out on a T2-weighted high resolution image obtained at the beginning of the intervention. The subsequent energy deliveries were placed on the target delineation. Therefore, the shift after 10 minutes would have the following consequences for the subsequent therapy :

• Parts of the tumor would have not been ablated, or the surgical margin would have not been sufficient.

• Shifts in-between shots potentially lead to small areas of viable tissue in between the ablation cells. Without motion compensation these have to be accounted for by a large overlap of the sonication cells, which significantly decreases the achievable ablation speed and thus the overall duration of the intervention.

It is important to note that during the in-vivo experiment no motion estimation or correction has been performed. The estimation of the displacements was made in post-processing using the 3D anatomical images acquired during therapy. The resulting motion vectors were used to retrospectively update the position and configuration of the sonication cell cluster such that it will match the actual ablated anatomy. The small residual shift of ∼ 0.3 mm between the updated cell cluster and the motion corrected thermal dose and of ∼ 1 mm between the updated cell cluster and the motion corrected NPV is a good indicator that the proposed motion estimation procedure can be successfully applied in-vivo with acceptable error margins.

The current multithreaded C++ implementation of the 3D registration algorithm used in our study has an execution time of approximately 30 s. Although motion compensation was performed in post-processing, the time interval between energy deliveries, which is usually 6-8 minutes, allows an integration. Furthermore, acceleration using an implementation via the compute unified device architecture (CUDA) platform on a graphical processing unit (GPU) 25 allows to bring the motion estimation down to ∼ 10 seconds, which renders the approach feasible with the overall work-flow.

The motion compensation framework proposed in this study is, however, not limited to HIFU To address this problem, a 3D anatomical image could be acquired, for example, after each fraction and registered to a 3D anatomical anchor acquired before delivering the first fraction. The resulting motion vector fields could then be integrated into the therapy planning algorithm and the plan could be re-optimized prior to the next fraction. The proposed framework can also be adapted to operate during the delivery of a single fraction. In order to reduce the damage caused to the OAR's, a single fraction is usually delivered over several changes in the gantry position. A 3D anatomical image could be acquired while changing gantry positions and registered to a 3D anatomical anchor acquired prior to starting dose delivery. The estimated displacements could be used, for example, to re-optimize the plan on-the-fly. However, such an approach is timewise challenging since the gantry positional changes usually take place over time intervals of only several seconds.

V. CONCLUSIONS

So far, studies concerned with motion estimation during MRg-HIFU interventions in the abdomen were focused on tracking the targeted area only over the time course of the therapeutic energy delivery. At this level, the interest was compensating for motion mainly caused by respiration. In the presented study, we demonstrated that motion compensation only during energy delivery is conceptually not sufficient : HIFU therapies rely on a non-invasive target identification and delineation, the subsequent interventional planning of the ablation areas, the therapeutic energy delivery and the validation of the therapeutic endpoint. One of the strengths of MR-guided interventions is that except the initial target identification and delineation, all subsequent steps can be carried out on-the-fly with the patient in place on the interventional table. The advantage of the proposed approach is that MR-guided interventions allow to adjust planning and therapy directly to the current anatomical position for the particular position of the patient on the interventional table.

While static parts of the human anatomy such as the cranium or the limbs are relatively easy to immobilize (stereotactic frames, limb molds), the thoracic and the abdominal area are subject to both respiratory motion and slow physiological drifts. As we have demonstrated, this can lead to a misalignment between interventional planning, energy delivery and therapeutic validation. Since this can potentially lead to both undesired tissue damage of organs at risk and viable residual tumor tissue, a suitable motion correction strategy should encompass all of the essential steps of the intervention: Planning, monitoring of the therapeutic energy delivery and the validation of the therapeutic endpoint. The presented method is a first step towards such a motion correction strategy for interventions in the abdomen.

Fig. 1

 1 Fig. 1 HIFU ablation protocol employed in our study. The 3D T1-weighted abdominal scans performed at regular time intervals allow the estimation of the 3D displacements exhibited by liver and kidneys due to slow physiological motion.

Figure 2

 2 Figure 2 illustrates the evolution of the anatomy changes induced by slow physiological motion in one of the volunteers abdomen. The displayed figures correspond to a coronal (figure 2(a)) and respectively a sagittal (figure 2(b)) slice extracted from a set of images obtained by subtracting from the 3D anatomical anchor subsequent 3D anatomical images.The leftmost images in the two figures depict the anatomy present in the two slices as captured by the 3D anatomical anchor. In this particular data set, considering just the liver and kidneys, the superior part of the liver experiences an obvious displacement in the craniocaudal direction. The peak of this displacement takes place at the acquisition time of the last 3D anatomical image. Given the image's voxel size of 2 × 2 × 2 mm 3 , this displacement is measured to be approximately 13.5 mm. In comparison, from a visual point of view, the displacements in the left-right and anterior-posterior directions were deemed to be negligible.

Fig. 2

 2 Fig. 2 Set of images that result from the subtraction between the 3D anatomical anchor and the 3D anatomical image acquired t minutes after the acquisition of the anatomical anchor in one of the volunteers. The figures correspond to a (a) coronal (b) sagittal slice from the resulting 3D difference images. The leftmost slices give a visual depiction of the anatomy present in the difference images.

Fig. 3

 3 Fig.3The relationship between the SNR of the image to register and the MSE of the corresponding estimated vector field with respect to the gold standard.

Figure 7

 7 Figure 7 illustrates the slow physiological drifts experienced by the porcine liver over the duration of the HIFU intervention. The displayed graphs were constructed in a similar fashion as the ones presented in Figure 4. As opposed to the observations made during the volunteer study, there is no clear tendency of the displacements to increase over time in either of the directions. Nevertheless, displacements of up to 4 mm can be observed in the cranio-caudal direction.

Fig. 4

 4 Fig. 4 Evolution of the liver and kidneys displacement caused by slow physiological motion over a duration of 80 minutes in ten healthy volunteers. (a) -(f) Each boxplot is a representation of a joint set containing the displacements estimated on all volunteers after k × 10 minutes since the acquisition of the 3D anatomical anchor, with k = 1, . . . , 8. (g), (h) The median of the magnitudes of the displacement vectors estimated for the liver and respectively the kidneys of each individual volunteer. The values on the abscissa correspond to the time in minutes since the acquisition of the 3D anatomical anchor. The values on the ordinate are the signed displacements in millimeters in the cranio-caudal, left-right and anterior-posterior directions.

Fig. 5

 5 Fig. 5 Validation of the optical flow algorithm. The graph illustrates the correspondence between the displacements obtained by manually tracking a set of vessel bifurcations in the liver and kidneys (on the abscissa) and the displacements estimated by the optical flow algorithm (on the ordinate). The displacements were determined for all volunteers in three directions: cranio-caudal (blue circles), left-right (red squares) and anterior-posterior (black crosses). The purple line corresponds to the ideal location of the plotted points.

Fig. 6

 6 Fig. 6 Positional errors that would affect the intervention in the selected landmarks, in the presence and respectively in the absence of a motion correction scheme. The values on the ordinate correspond to the positional errors in millimeters computed according to II.D.1.

Fig. 7

 7 Fig. 7 Evolution of the liver displacements caused by slow physiological motion during the in-vivo HIFU pig experiment. The values on the abscissa correspond to the time delay in minutes since the acquisition of the 3D anatomical anchor. The values on the ordinate are the signed displacements in millimeters in the (a) cranio-caudal, (b) left-right and (c) anterior-posterior directions.

Figure 8 (

 8 Figure 8(c) displays the location of lethal thermal dose after registering the thermometry images. The distance between the centers of mass of the position corrected sonication cell cluster and the region determined to have received a lethal thermal dose is reduced from ∼ 2 mm without motion correction to ∼ 0.3 mm with motion correction.

Figure 9 (

 9 Figure9(a) displays a slice from the contrast enhanced image acquired at the HIFU therapy endpoint. The NPV resulting from the 10 sonications is enclosed in a red rectangle and magnified for better visibility. A mask covering the NPV was manually extracted from this image and overlaid on the planning scan as a high intensity region as depicted in figure9(b). For positional reference, the corrected sonication locations were also overlaid as red

Fig. 8 (

 8 Fig. 8 (a) Planned cluster of sonication cells (yellow circles) together with the motion corrected location of the sonications (red circles) overlaid on the corresponding anatomy from the planning image. (b), (c) The anatomy that was determined to have received a lethal thermal dose (depicted by the high intensity region): (b) before registering the thermometry images (c) after registering the thermometry images.

Fig. 9 (

 9 Fig. 9 (a) Slice of the contrast enhanced scan acquired at the end of the HIFU intervention upon which the resulting NPV is visible. The NPV resulting from the sonications is enclosed in a red rectangle and magnified for better visibility. (b), (c) Corrected sonication locations overlaid as red circles on the corresponding anatomy from the planning image. The high intensity region in the images corresponds to the NPV: (b) without motion correction (c) with motion correction.

  interventions. With the development of the first MR linear accelerator (MR-LinAc) prototypes[REF], our motion correction strategy could potentially be integrated in the work-flow of an MR-guided radiotherapy (MRg-RT). Similar to HIFU interventions, energy deliveryduring external beam radiotherapy is usually episodic, with the total prescribed radiation dose being divided into multiple fractions. The time interval in between fractions, which can sometimes extend over several days, allows the healthy tissues to partially recover from the effects of the absorbed radiation. These time intervals can be seen as the equivalent of the cool-down periods during HIFU therapies. Naturally, there is a high probability that due to the discussed physiological drifts or due to patient positional errors, the target location will differ from one fraction to the next. If these displacements are not taken into consideration, they might lead to the under-treatment of the tumor or to additional damage to the OAR's.

ACKNOWLEDGMENTS

This work was supported in part by the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and in part by the Dutch Technology Foundation (STW)

(project OnTrack #12813).


Maxence Borot De Battisti

M Maenhout

Baudouin Denis De Senneville

Gilion Hautvast

D Binnekamp

Jan J W Lagendijk

Marco Van Vulpen

M A Moerland

An automated optimization tool for high-dose-rate (HDR) prostate brachytherapy with divergent needle pattern

Introduction

Focal-HDR has gained an increasing interest as an alternative to whole-gland therapy in patients with localized prostate cancer. Focal therapy aims at reducing treatment related side effects and toxicity [START_REF] Pieters | Comparison of three radiotherapy modalities on biochemical control and overall survival for the treatment of prostate cancer: A systematic review[END_REF]. The success of focal therapy will depend on the ability of tumour localization and dose delivery. Imaging, pathology and dose delivery studies have shown the value of multi-parametric MR-imaging for tumour localization [START_REF] Groenendaal | Validation of functional imaging with pathology for tumor delineation in the prostate[END_REF]. Furthermore, [START_REF] Polders | U A Establishing implantation uncertainties for focal brachytherapy with I-125 seeds for the treatment of localized prostate cancer[END_REF] showed an adequate dose coverage is obtained if a 5mm-margin is applied to the MR based tumour delineations. Therefore, a MRI-guided focal-HDR procedure is under investigation. In the daily practice at the UMCU, needles are inserted under ultrasound guidance while needle reconstruction, dose planning and needle positions verification are based on MRI. For optimal MR guidance during therapy, needle insertion should also be under MRI guidance, which is not currently possible due to access restrictions in a closed MR-bore system. To enable MR-guided needle insertion in the space restricted MR environment, MR-compatible robotic devices have been developed at several institutes: [START_REF] Fischer | Robotic Assistant for Transperineal Prostate Interventions in 3T Closed MRI[END_REF][START_REF] Fischer | MRI-Compatible Pneumatic Robot for Transperineal Prostate Needle Placement[END_REF] and [START_REF] Dimaio | Robot-assisted needle placement in open MRI: system architecture, integration and validation[END_REF] recently designed robotic assistants for transperineal prostate needle placement. At the UMCU, a robotic device that can automatically insert needle into the patients prostate under MRI guidance is currently under development [START_REF] Van Den Bosch M R Moman | MRI-guided robotic system for transperineal prostate interventions: proof of principle[END_REF]. This robot is placed between the legs of the patients inside the MR bore. A tapping mechanism is used for needle insertion to restrict prostate movement and tissue deformation [START_REF] Lagerburg | A new robotic needle insertion method to minimise attendant prostate motion[END_REF]. Furthermore, needles are inserted under different angles in a divergent way, from a single rotation point. This rotation point is placed just beneath the perineal skin to have access to the whole gland. With this method, it is warranted to deliver the irradiation dose, needle per needle, guided by an adaptive planning system that takes anatomy changes and needle deviations into account. The ultimate goal is therefore to develop a procedure where the parameters of the dose plan are re-optimized after each needle insertion according to the perturbations of the set-up (change of anatomy, needle bending,...). For this purpose, a fast, accurate and stable optimization algorithm is important for inverse planning, enabling the implementation of dose-adaptive focal-HDR brachytherapy in the future. Currently, there is no commercial treatment planning system available which allows divergent needle insertion and adaptive planning. The development of such a treatment planning method is of great importance to find the most optimal implant in regard to planning target volume (PTV) coverage and sparing of organs at risk.

For the planning of brachytherapy intervention, various optimization methods have been developed. The most common algorithms deal with the dwell-times distribution within already implanted catheters such as geometrical optimization (GO) [START_REF] Visser | Optimization of interstitial volume implants[END_REF] or inverse planning simulated annealing (IPSA) (Lessard andPouliot 2008, Alterovitz et al 2006). Some algorithms have been developed to determine the distribution of catheters within the prostate such as the 2D Centroidal Voronoi Tesselations (CVT) algorithm [START_REF] Poulin | Adaptation of the CVT algorithm for catheter optimization in high dose rate brachytherapy[END_REF]. This CVT algorithm optimizer has been described only for parallel needle configurations. [START_REF] Siauw | NPIP: A skew line needle configuration optimization system for HDR brachytherapy[END_REF] described a needle planning by integer program (NPIP) algorithm to generate needle configurations that avoid critical structures near the penile bulb and other healthy structures, and avoid needle collisions inside the body. It was not applied for a divergent needle pattern with a single rotation point. More recently, hybrid algorithms which optimize the catheter positions and the dwell times have also been under investigation. [START_REF] Holm | Heuristics for integrated optimization of catheter positioning and dwell time distribution in prostate HDR brachytherapy[END_REF] has developed a heuristics method for catheter positioning and dwell time distribution for parallel needle pattern with a runtime of 1h without constraints regarding the nonperforation of urethra by the needles. [START_REF] Gorissen | Mixed integer programming improves comprehensibility and plan quality in inverse optimization of prostate HDR brachytherapy[END_REF] has developed a hybrid optimizer by using mixed integer programming with a runtime of several minutes.

All mentioned optimizers have certain drawbacks: First, most of them use parallel needle insertions and are not applicable for divergent needle insertion with a single point of insertion. Therefore, they cannot be used to determine the needle angulations or the position of the center of insertion. Second, they usually require numerous iterations and may produce sub-optimal dose results due to the trapping in local minima regions of the cost function landscape. Third, most optimizers require the manual determination of several input parameters as minima or maxima for PTV and organs at risk. In clinical practice however, dose coverage to the PTV and dose to the organs at risk are used to evaluate if the plan is clinically acceptable and therefore it is desirable to perform an optimization based on these parameters. Fourth, most optimizers usually require weight penalties (or importance coefficients) as input (see [START_REF] Alterovitz | Optimization of HDR brachytherapy dose distributions using linear programming with penalty costs[END_REF][START_REF] Hsu | Comparison of inverse planning simulated annealing and geometrical optimization for prostate high-dose-rate brachytherapy[END_REF]. These values depend on the patient anatomy. Therefore, weight penalties are not intuitive for a clinician and often need individual adjustment to obtain an acceptable plan. It is therefore necessary to generate a dose plan which will minimize the dose deposition error resulting from the weight penalty dependencies. Finally, in order to be eligible for intra-operative use, the total calculation time to obtain a clinically acceptable plan must be less than several minutes. This is usually not the case for most optimizers.

The aim of this article is to describe the development of a fully automatic inverse dose planning optimization tool for MRI-guided focal-HDR brachytherapy on prostate with a divergent needle pattern. Its goal is to determine the optimal center of rotation, needle angulations, source positions and dwell times within a reasonable time needed for intra-operative use. The optimizer will be tested in a planning study by assessing the dose volume parameters.

Methods

Specification of the optimizer

The coordinate system for the anatomy, the needle and the source positions is defined as follows. Let r = (x, y, z) be the spatial coordinate where x, y and z correspond to the positions in the left-right, the anteroposterior and the inferior-superior direction, respectively. For this divergent needle technique with the needle tracks coming from a single entry point, let r rot = (x rot , y rot , z rot ) be the position of the center of rotation of the set-up and (θ i , φ i ) the angle of the i th needle insertion in the spherical coordinate system. In the common practice of HDR brachytherapy, the distance ∆ between the dwell positions along the needle is constant. Due to the finite size of the needle, the number N source of possible dwell positions of the source along the needle is limited. The dwell positions r i k (r rot , θ i , φ i ), k ∈ [1, N source ] of the i th needle insertion can then be expressed in Cartesian coordinates as follows:

r i k (r rot , θ i , φ i ) =    x rot + k∆ sin(θ i ) cos(φ i ) y rot + k∆ sin(θ i ) sin(φ i ) z rot + k∆ cos(θ i )    (1)
To develop a fully automatic optimizer for a given number of divergent needle insertions (referred to as N needle needles in the scope of this study), the following parameters need to be optimized: (1) the position of the center of rotation r rot = (x rot , y rot , z rot ); (2) the angles of the needle tracks (θ

1•••N needle , φ 1•••N needle ) in
the spherical coordinate system; and (3) the dwell times of the sources t i k at the source position k

(k ∈ [1, N source ]) of the i th needle insertion (i ∈ [1, N needle ]). Let p = (r rot , θ 1•••N needle , φ 1•••N needle , t 1 1•••Nsource , • • • , t N needle 1•••Nsource
) be the vector containing the parameters of the set-up to be optimized and Ω its corresponding set of feasible solutions.

Finally, for a given dose plan, let D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl be the dose received by 95% of the PTV, by 10% of the urethra and by 1cm 3 of rectum and bladder respectively. The optimizer must fulfill the following objectives: 1) highest coverage of the PTV (i.e. a large D 95% P T V ) and 2) dose on defined organs at risks (OAR) as low as possible (i.e. a small D 10% U r , D 1cc Rec and D 1cc Bl ).

The goal of the proposed optimizer is to determine the optimal set of parameters

p optimal = (r rot , θ 1•••N needle , φ 1•••N needle , t 1 1•••Nsource , • • • , t N needle 1•••Nsource )
optimal to obtain the desired coverage in the PTV without exceeding the constraints of the organs at risk. Therefore, the following constraints function as input of the optimizer:

D min 95% P T V , D max 10% U r , D max 1cc Rec
and D max 1cc Bl . They correspond to the minimum value of D 95% P T V and the maximum value of D 10% U r , D 1cc Rec and D 1cc Bl respectively in order to obtain a clinically acceptable and optimal dose plan. Concretely, a dose plan is clinically acceptable when:

D 95% P T V > D min 95% P T V ,D 10% U r < D max 10% U r , D 1cc Rec < D max 1cc Rec and D 1cc Bl < D max 1cc Bl .
As output, the optimizer should give the optimal parameters of the set-up p optimal which corresponds to the maximum D 95% P T V in combination with the minimization of D 10% U r , D 1cc Rec and D 1cc Bl .

Dose computation and proposed optimization workflow

According to the guideline of the AAPM Task group No. 43 [START_REF] Nath | Dosimetry of interstitial brachytherapy sources: Recommendations of the AAPM Radiation Therapy Committee Task Group[END_REF][START_REF] Rivard | Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations[END_REF], the dose D(r, p) received at r is expressed as the sum of the contribution of all source positions:

D(r, p) = N needle i=1 Nsource k=1 d i k (r, r rot , θ i , φ i )t i k (2)
where d i k (r, r rot , θ i , φ i ) is the dose-rate of source position k of the i th needle insertion and, depending on the model of the dose distribution chosen (point source, line source,...), is usually a complex non-linear function of r, r rot , θ i and φ i .

A way to obtain the desired dose plan is to approach the dose to a given value D opt (r) for all points r by solving the following equations:

p optimal = argmin p∈Ω [C(p)] with C(p) = r∈R ω(r)[D(r, p) -D opt (r)] 2 dr (3)
where C(p) is called the cost function, p optimal are the optimal parameters defined previously and ω(r) is the weight coefficient at the point r which will be detailed in section 2.5. The strategy of the proposed optimizer benefits from the linear impact of the dwell times on the deposited dose. If the center of rotation and the angles of the needle tracks are fixed, determining the value of t i k goes back to solving a set of linear equations and it is thus feasible to find a direct and efficient solution. However, solving Eq. 5 for the variables of the needle positioning (r rot , θ 1•••N needle and φ 1•••N needle ) is highly non-convex and consequently difficult to overcome. Therefore, the proposed method relies on the exact determination of the optimal dwell times using the resolution of linear equations while the remaining variables (position of center of rotation and angulations of needles) are deduced using heuristic or exhaustive searches. Since the area of insertion of the needle is reduced to several square centimeters in order of magnitude, an exhaustive enumeration of r rot is employed. However, regarding the determination of the needle orientations (θ i , φ i ) , an exhaustive optimization is highly time-consuming in practice. Therefore, the idea of evaluating a certain number of heuristics N heur for each center of rotation selected was employed to accomplish this task. The complete workflow of the optimizer is depicted in Figure 1. describes the proposed heuristics for the determination of needle orientations. First, the PTV and the urethra are projected from the center of rotation onto a transverse plane (cf. schematic on the left). A uniform distribution of needles is then found by applying k-means clustering on P P T V until no clusters are found in the urethra zone (cf. schematic on the right).

Task 1: Exhaustive enumeration of the rotation center of the set-up

The first task in the workflow is to select the position of the rotation center of the set-up (i.e. r rot ).

The robotic device currently in development at the UMCU for HDR brachytherapy with divergent needle pattern is made such that the center of rotation should be placed at the perineum (to avoid multiple insertion points and to have full access to the prostate). Therefore, r rot is supposed fixed in the inferior-superior direction (i.e. z-axis) such that z rot is determined by the perineum. Furthermore, the ranges of robot movement are 2cm along the x and y-axis due to the restricted space between the legs. Thus, the values of x rot and y rot are initially chosen such that the point P is in line with the center of the largest PTV contour in the transversal plane. An exhaustive search for the optimal center of rotation was performed by evaluating 9 possible candidates on a 1cm grid (see Figure 2(a)).

Task 2: Heuristic approach of the needle orientations

At this point, the center of rotation is fixed. The next step is to construct several heuristics for the needle orientations. The heuristic approach proposed in this manuscript consists in finding a uniform needle distribution in space with the additional condition that none of the needle tracks perforates the urethra.

For this purpose, the method of [START_REF] Poulin | Adaptation of the CVT algorithm for catheter optimization in high dose rate brachytherapy[END_REF] for parallel needle set-up was adapted as follows: the PTV and urethra volume are projected from the center of rotation on a transverse plane behind the PTV (see Figure 2(a)). The latter projections are noted P P T V and P U rethra respectively. The aim is to find a uniform distribution of N needle points in P P T V without any of those points inside P U rethra . First, N needle initial generators are positioned randomly in P P T V . Then, k-means clustering is applied on the indices of the surface P P T V with the additionnal condition that no cluster centers must stand in P U rethra . The final cluster centers reflect the positions of the needle (see Figure 2(b)).

This method of heuristic selection relies on a random initial needle set up, therefore several heuristics were evaluated. The number of heuristics evaluated N heur in section 2.7.2.

Task 3: Linear optimization of dwell times

While entering this step, the center of rotation and the needle orientations are fixed. The determination of the dwell times can now be expressed using a linear optimization problem as follows.

In the following, N tot source is the total number of active source positions and t m and

d m (r) (m ∈ [1, N tot source ]
) are the dwell time and dose rate at each active source position respectively.

The final objective is a high PTV coverage while the OAR are spared as much possible. A way to reach this goal is to approach the dose to a certain value D P T V in the PTV and to 0 in the organs at risk. Considering a number N OAR organs at risk, the weight coefficients are supposed constant for the PTV and OAR l , l ∈ [1, N OAR ] and are noted ω P T V and ω OAR l respectively in the following. In a discrete space, the cost functions of the PTV and the OAR l can be defined by referring to the definition of C(p) in Eq. 3, as follows:

C P T V (p) = ω P T V V P T V r∈P T V [D(r, p) -D P T V ] 2 (4) ∀l ∈ [1, N OAR ], C OAR l (p) = ω OAR l V OAR l r∈OAR l D(r, p) 2 (5)
These cost functions were divided by the volume of the organs to avoid volume dependency in the optimization. Therefore the basic cost function could be expressed as:

C(p) = C P T V (p) + N OAR l=1 C OAR l (p) (6)
The optimal values of ω OAR l•••N OAR , ω P T V and D P T V are obviously dependent of the anatomy of the patient. Therefore, in the following sections, the algorithm to minimize C(p) for a given value of ω OAR l•••N OAR , ω P T V and D P T V is presented (section 2.5.1). In section 2.5.2, the method to determine the optimal values of ω OAR l•••N OAR , ω P T V and D P T V according to the patients anatomy is described.

2.5.1. Solution using matrix inversion. In this section, a direct and fast method to determine the dwell times of the source without unrealistic negative values is described for fixed values of ω OAR l•••N OAR , ω P T V and D P T V . This problem of finding the optimal source dwell times for fixed needle positions is analogous to finding the optimal beam intensities in external beam radiotherapy: this reverts to finding a solution of the inverse problem which does not yield unphysical negative values. Therefore, a similar algorithm described by [START_REF] Goldman | Feasibility of a fast inverse dose optimization algorithm for IMRT via matrix inversion without negative beamlet intensities[END_REF]2009) in the case of Intensity Modulated Radiation Therapy (IMRT) known as Fast Inverse Dose Optimization (FIDO) is used. Goldman described a method to obtain a direct solution of the inverse problem that avoids negative beamlet weights. It involves reformulating the organs at risk cost functions C OAR l (p): in Eq. 5, D(r, p) 2 is replaced by

N tot source m=1 [d m (r)t m ] 2
. By modifying the organs at risk cost functions, C OAR l (p) will not be null through destructive interference effects between dwell times and most of the unphysical negative solutions are therefore excluded. Furthermore, the optimization problem is reduced to a set of linear equations.

With this method, the optimal dwell times are obtained by the matrix inversion:

T = α -1 β (7) with T = (t m ) m∈[1,N tot source ]
is the vector of N tot source elements containing the dwell times at all active source positions. β is also a vector of N tot source elements and α is a N tot source ×N tot source matrix defined respectively as follows:

∀m ∈ [1, N tot source ], β m = ω P T V D P T V V P T V r∈P T V d m (r) (8) ∀m, n ∈ [1, N tot source ], α mn = ω P T V V P T V r∈P T V d m (r)d n (r) + N OAR l=1 ω OAR l V OAR l r∈OAR l d m (r)d n (r)δ mn (9)
where δ mn is the Kronecker delta function.

2.5.2. Exhaustive search of the weight coefficients (ω P T V , ω OAR 1•••N OAR ) and the dose approached for the PTV (D P T V ). As the matrix inversion in Eq. 7 is little time consuming to execute, it is possible to obtain the overall results by an exhaustive search of ω OAR l•••N OAR , ω P T V and D P T V . Consequently, no individual adjustment of these parameters is needed to obtain an acceptable plan and the dwell times found will be optimal for any anatomy. The exhaustive optimization of ω OAR l•••N OAR , ω P T V and D P T V will be investigated in section 2.7.1.

Task 4: Selection criterion

The exhaustive enumeration of Task 1 (center of rotation of set-up), the several heuristics generated in Task 2 (needle orientations) and the exhaustive search ω OAR l•••N OAR , ω P T V and D P T V in Task 3 (dwell times) will gives several dose plans. The objective is now to select the best plan. The following criterion is proposed: the best plan has the highest value of energy E (in Gy) defined as:

E = min(A, B, C, D) (10)
with the relative parameters:

A = D 95% P T V -D min 95% P T V B = D max 10% U r -D 10% U r C = D max 1cc Rec -D 1cc Rec D = D max 1cc Bl -D 1cc Bl
The parameters A, B, C and D represent the relative difference between the dose coverage parameters and the clinical constraints (set as input) of the PTV, the urethra, the rectum and the bladder, respectively. The higher the value of A, the better the dose coverage of the PTV. Moreover, the clinical constraint at the PTV is achieved if, and only if A > 0. Consequently, the minimum value over A, B, C and D in Eq. 10 corresponds to the maximal dose covering error with respect to the clinical constraints within the region that receives the "worst" dose deposition. E therefore represents the quality of the dose plan: the greater E is, the better the dose plan becomes (the dose plan is clinically acceptable if, and only if E > 0).

In conclusion, by maximizing E, the plan which has all parameters optimized (maximum D 95% P T V in combination with minimization of D 10% U r , D 1cc Rec and D 1cc Bl ) is selected.

Experimental validation

The different steps in the optimizer workflow were assessed before the complete optimizer workflow was validated. In the first experiment, the dwell time optimization (see section 2.5) was analyzed by assessing the automatic search of ω OAR l•••N OAR , ω P T V and D P T V . In the second experiment, regarding the determination of the needle angulations (see section 2.4), the required number N heur of heuristic to evaluate was analyzed. The last experiment consisted in testing the complete optimizer workflow in a planning study by assessing the dose volume parameters.

The experiments were performed in a retrospective evaluation using clinical data from 10 patients. The delineations of the prostate tumors and the OAR considered (urethra, bladder, rectum and rest of the tissues) were made on a 1mm 3 resolution MRI image by an experienced oncologist. The PTV volumes ranged from 8.5cm 3 to 23.3cm 3 with a median of 16.1cm 3 . For all experiments, the clinical values as input were: D min 95% P T V = 19Gy, D max 10% U r = 21Gy, D max 1cc Rec = 12Gy and D max 1cc Bl = 12Gy. Those are the clinical constraint values for single fraction HDR brachytherapy as monotherapy used at the UMCU. For all experiments, varying numbers of needle insertions (from 2 to 14) were tested.

Regarding the source position the common procedure at the UMCU was adopted: for each needle insertion, the active source center positions were separated by a step-size of ∆ = 2.5mm and situated inside the PTV with an extra margin of 3mm.

For dose calculation, the dose rate was calculated using the point source approximation model due to the minimum time of computation, with a small adaptation as follows to avoid over-optimization of the dose close to the source:

d i k (r, r rot , θ i , φ i ) = S K Λg P [R i k (r, r rot , θ i , φ i )]Φ an [R i k (r, r rot , θ i , φ i )] R 0 2 R i k (r, r rot , θ i , φ i ) 2 + exp[-R i k (r, r rot , θ i , φ i ) 2 ] (11)
where S K is the air-kerma strength, Λ the dose-rate constant in water, Φ an (R) the one-dimensional anisotropy factor, R 0 denotes the reference distance which is specified to be 10mm, g P (R) corresponds to the radial dose function in the case of point source approximation model and R i k (r, r rot , θ i , φ i ) is the distance (in millimeters) between the source position r i k (r rot , θ i , φ i ) and

r (R i k (r, r rot , θ i , φ i ) = ||r i k (r rot , θ i , φ i ) -r|| 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 k (r, r rot , θ i , φ i ) 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) (Λ = 1.108cGy.h -1 .U -1 ) and an arbitrary source strength S K = 40.80mGy.h -1.m 2 was chosen. The multiplication of the radial dose function g P (R) and the anisotropy factor Φ an (R) were approximated by a 2 nd order polynomial fit (g

P (R) • Φ an (R) = a 0 + a 1 R + a 2 R 2 ).
The coefficients for the fit were a 0 = 1.11, a 1 = -3.30 • 10 -3 and a 2 = 3.12 • 10 -6 , where R is in millimeters.

As well as an evaluation of dose to target and OAR, dose homogeneity and conformity were also investigated. The parameters homogeneity (HI 150% and HI 200% ) and conformal index (COIN ) are defined respectively as:

HI i% = V 100% P T V -V i% P T V V 100% P T V with i ∈ {150, 200} (12)
COIN = (V 100% P T V ) 2 V P T V V 100% body (13)
where V i% P T V (V 100% body respectively) is the volume inside the PTV (the total volume respectively) that receive i% (100% respectively) of the prescribed dose i.e. D min 95% P T V ). The HI i% measures the volumes fraction that receives between 100% and i% (i ∈ {150, 200}) of the prescribed dose and the COIN compares the reference dose coverage in the PTV with the total volume in the reference isodose volume.

2.7.1. Assessment of the exhaustive search of ω P T V , ω OAR 1•••N OAR and D P T V ) for the dwell times optimization (Task 3). In this experiment, the automatic search of ω OAR l•••N OAR , ω P T V and D P T V is assessed. The center of rotation was fixed in line with the center of the largest PTV contour in the transversal plane and the mid-plane of the tumor. A heuristic for the needle orientations was determined as described in section 2.4. The dwell times for all source positions were calculated using Eq. 7 for all possible combinations of weights coefficients of PTV, urethra, bladder, rectum and other tissues (respectively ω P T V , ω U r , ω Bl , ω Rec and ω tissue ) and D P T V such that:

D P T V ∈ {20, 21, • • • , 79, 80} (14) log 10 (ω P T V ) ∈ {0, 0.001, • • • , 9.999, 10} (15)
(log 10 (ω U r ), log 10 (ω Bl ), log 10 (ω Rec ), log 10 (ω tissue )) ∈ {0, 1, • • • , 9, 10} 4 (16)

with the additional condition that:

log 10 (ω P T V ) + log 10 (ω U r ) + log 10 (ω Bl ) + log 10 (ω Rec ) + log 10 (ω tissue ) = 10 ( 17)

to avoid redundancy of results. From all resulting dose plans, the best one was selected with the selection criterion defined in section 2.6

In the following experiments (described in section 2.7.2 and 2.7.3), to obtain the optimal dwell time, the optimal values of ω P T V , ω U r , ω Bl , ω Rec and ω tissue and D P T V were determined by applying the same exhaustive search. However, in order to decrease the time of calculation, the exhaustive search range of ω P T V was reduced to:

log 10 (ω P T V ) ∈ {0, 1, • • • , 9, 10} (18)
2.7.2. Analysis of the required number of heuristic for the needle angulations to be evaluated in Task 2. In this experiment, the number of heuristic evaluated N heur in Task 2 was assessed. The complete optimizer workflow was performed several times without an exhaustive search of the center of rotation (which was fixed in line with the center of the largest PTV contour in the transversal plane and the mid-plane of the tumor) for N heur varying from 1 to 80. The energy, E, of the final dose plan was determined as a function of N heur . For the following experiment described in section 2.7.3, N heur was fixed to 10 to limit the time of calculation. 2.7.3. Analysis of the number of inserted needles. In this experiment, the optimizer was tested in a planning study by assessing the dose volume parameters. p optimal was determined to obtain the desired coverage without exceeding the constraints of the organs at risk for 2, 4, 6, 8, 10, 12 and 14 needle insertions by using the earlier mentioned constraints as inputs.

Results

In this section, the results of the experiment described in section 2.7 are presented below. Firstly, the results of the exhaustive search of ω P T V , ω U r , ω Bl , ω Rec and ω tissue and D P T V to obtain optimal dwell times are described in section 3.1. Regarding the optimization of the needle angulations, the results of the analysis of the required number N heur of heuristic evaluations are described in section 3.2. Finally, the results of the planning study for one typical patient and afterwards, for all patients are shown in section 3.3. 3.1. Assessment of the exhaustive search of ω P T V , ω OAR 1•••N OAR and D P T V ) for the dwell times optimization (Task 3)

In this section, the results of the experiment described in section 2.7.1 are presented in Figure 3 for one patient with 10 needle insertions as a typical example.

The best dose plan (i.e. the one maximizing E) was found with the following parameters: {D P T V , ω P T V , ω U r , ω Bl , ω Rec , ω tissue } = {75, 10 0.587 , 10 2 , 10 4 , 10 2 , 10 1.413 }. The values of D 95% P T V , D 10% U r , D 1cc Rec , D 1cc Bl and E are plotted in Figure 3 as a function of D P T V and log 10 (ω P T V ). For the representation, {ω U r , ω Bl , ω Rec } were set to their optimal values {10 2 , 10 4 , 10 2 } and ω tissue was deducted from Eq. 17.

Figure 3(a), 3(b), 3(c) and 3(d) show an increase of the dose deposition in all the different volumes considered (PTV, urethra, bladder and rectum) for increasing D P T V and log 10 (ω P T V ). In Figure 3(e), the best solutions (where E is maximal) are disposed on a continuous line.

For log 10 (ω P T V ) > 1.42, one or more dwell times found with FIDO had negative values. These nonphysical solutions are due to the resolution of the equation: by modifying the organ at risk cost functions in FIDO, most of the solutions which give negative dwell times are removed but not all of them. In this section, the results of the experiment described in section 2.7.2 are presented in a situation of a large number of needle insertions (N needle = 14): This corresponds to the worst case scenario where the calculations are the most timeconsuming. For each patient, the energy E of the selected dose plan was calculated within the number of heuristics evaluated. The average energy E on all 10 patients was plotted as a function of N heur (Figure 4 and the calculation time of the complete optimization workflow was also represented.

Analysis of the required number of heuristic for the needle angulations to be evaluated in Task 2

Figure 4 shows a fast initial increase of E with the number of evaluated heuristics.

Afterwards, the gradient decreases progressively. For N heur = 10, E achieves 64% of its value for N heur = 80. Furthermore, the number of heuristics evaluated in Task 2 is directly proportional to the calculation time of the complete optimization workflow. The time of optimization for N heur = 10 was approximatly 20 minutes on a PC with a 3.10GHz Intel R CoreTM i5-2400 processor and 8GB RAM using MATLAB R2013a.

Analysis of the number of inserted needles

The dose distribution was calculated on 10 patients (PTV volumes ranged from ranged from 8.5cm 3 to 23.3cm 3 with a median of 16.1cm 3 ) for 2, 4, 6, 8, 10, 12 and 14 needle insertions by using the earlier mentioned constraints as input. For each case, D 95% P T V , D 10% U r , D 1cc Rec , D 1cc Bl , HI 150% , HI 200% and COIN was computed.

The result for one typical patient is shown in Figure 5(a)(b)(c) and (d). The dose constraints set as input are also represented. D 95% P T V increased with the number of needle insertions, and D 1cc Rec and D 1cc Bl decreased. Furthermore, for this patient, D 10% U r did not seem to have a clear trend. All the dose constraints were already reached at four needle insertions. Moreover, HI 200% and COIN increased with the number of needle insertions, but the influence of N needle on HI 150% was relatively small.

Figure 5(c) shows the MRI image of a sagittal plane of the same patient with delineations and figure 5(d) and 5(e) show the dose distribution in this same sagittal plane for 2 and 10 needle insertions respectively. The cold spots in the PTV and the hot Once again, D 95% P T V increased on average with the number of needle insertions, while D 10% U r , D 1cc Rec and D 1cc Bl . The large ranges of values were due to the different anatomies of the 10 patients tested. HI 200% and COIN also increased with the number of needle insertions but HI 150% did not show any clear trend. On average, a clinical acceptable plan is already reached by using four needle insertions.

Discussion

In this manuscript, an automatic inverse dose planning optimization tool for MRIguided focal-HDR brachytherapy on prostate with divergent needle pattern is proposed.

The aim was to determine the optimal parameters of the set-up (point of rotation, needles angles and dwell times) which corresponds to the maximization of D 95% P T V in combination with minimization of D 10% U r , D 1cc Rec and D 1cc Bl . For that, the linear impact of the dwell times on the deposited dose was exploited and the remaining variables were determined by evaluating several heuristics for the needle angulations and by an exhaustive search for the position of the point of rotation.

Unlike most optimizers such as HIPO or IPSA, the proposed method does not require individual adjustments of several input parameters such as minimum dose, maximum dose or weight coefficients for PTV and organs at risk to obtain an acceptable plan [START_REF] Dinkla | Comparison of inverse optimization algorithms for HDR/PDR prostate brachytherapy treatment planning Brachytherapy[END_REF]. Figure 3 illustrates the importance of an exhaustive optimization of the weight coefficients and approached dose at the PTV to determine the optimal dwell times. It shows that the dose deposition is made at the expense of all the other regions. In particular for this patient, the rectum played a very important role in the dose deposition. However, the dose deposition is a little less critical for the urethra and has little impact for the bladder. Thus, it is difficult to predict how the weight coefficients influence the overall dose plan. As an example, increasing ω P T V will certainly increase D 95% P T V , but its influence on D 10% U r , D 1cc Rec or D 1cc Bl is unknown: it could be insignificant or dramatic depending on the anatomy of the patient. However, the gain in energy would not significantly drop by using a larger step in the exhaustive search of log 10 (ω P T V ). Therefore, log 10 (ω P T V ) is varied from 0 to 10 with a step of 1 in order to reduce the time of optimization and still obtain a good dose plan.

An important step in brachytherapy and thus in the optimizer workflow is to find the optimal needle positions. The algorithms developed recently are usually based on iterative methods [START_REF] Holm | Heuristics for integrated optimization of catheter positioning and dwell time distribution in prostate HDR brachytherapy[END_REF][START_REF] Gorissen | Mixed integer programming improves comprehensibility and plan quality in inverse optimization of prostate HDR brachytherapy[END_REF][START_REF] Siauw | NPIP: A skew line needle configuration optimization system for HDR brachytherapy[END_REF] which cause two problems. First, such algorithms strongly depend on the initialization and therefore could produce sub-optimal solutions due to the trapping in local minima regions of the cost function landscape. The second problem is that, due to the high non-convexity of the problem, the optimization may require a long calculation time. To avoid these problems, a different approach was chosen: The angulation of the needle tracks is determined by evaluating several heuristics chosen carefully using the k-means clustering. According to Figure 4, for a given point of rotation, evaluating 10 heuristics already gives a good dose plan compared to 80 heuristics and a significantly better dose plan compared to one heuristic. By evaluating 10 heuristics, the total time of the optimization for a given number of needle insertions is below 20min. It is important to note that the program has been developed on Matlab and has not been optimized for speed. The proposed pipeline and the employed numerical schemes thus show great perspectives for a further reduction of the computation time using Graphical Processing Units (GPU's) in a parallel architecture. This will be studied in future work.

Finally, the planning study of the proposed optimizer workflow shows promising results. The clinical constraints set as input were reached on average with 4 needle insertions which is better compared to the current clinical procedure where 13 to 17 needles are usually inserted (see [START_REF] Hsu | Comparison of inverse planning simulated annealing and geometrical optimization for prostate high-dose-rate brachytherapy[END_REF], Menard et al 2004). More precisely, the increase of D 95% P T V as well as the decrease of D 1cc Rec or D 1cc Bl with the number of needle insertions for the patient shown in Figure 5(a) was expected. D 10% U r did not depict a trend with the number of needle insertions: it already fulfilled by far the constraints for the urethra set as input (D max 10% U r = 21Gy) and therefore it was not the parameter to be optimized in priority. Moreover, HI 200% and COIN , expressing the dose homogeneity and conformity respectively, increased with the number of needle insertions as expected. HI 150% did not show a clear trend. In this study, the PTV volumes had a median of 16.1cm 3 and all the dose constraints were already reached at 4 needle insertions on average. Therefore, those results are in line with the study of [START_REF] Steggerda | Minimizing the number of implantation needles for prostate 125 I brachytherapy: An investigation of possibilities and implications[END_REF] which shows, in the case of Low-Dose-Rate (LDR) brachytherapy, satisfactory dose coverage for, on average, 0.3 needles per cm 3 prostate volume. Furthermore, Vargas et al (2004) and Boyea et al (2007) showed the urinary toxicity following HDR brachytherapy is significantly increased by using more than 14 needle insertions. The cold spot in the PTV for 2 needle insertions as presented on Figure 5(d) was reduced by an increasing number of needle insertions (see Figure 5(e)). This illustrates how the number of needle insertions could improve the dose distribution. Figures 5(g A limitation in using less needle insertions is that the robustness of the dose plan may drop. It must be kept in mind that these are only simulations: in practice, there may be errors in needle positioning or in the position of the rotation point which may lead to modifications in the dose distribution. The impact of an error in position of one needle on the dose plan must be studied in future work. Furthermore, it is important to note that in this study, no modification of the anatomy (for example due to the insertion of the needle) was taken into account. [START_REF] Lagerburg | A new robotic needle insertion method to minimise attendant prostate motion[END_REF] however showed that the prostate motion was significantly less when using a robotic device that taps the needle compared to hand insertion.

A way to further increase the quality of the dose plan would be to couple this method with a gradient-based optimizer on all parameters at the same time (point of rotation, needles angles and dwell times) based on the cost function described in Eq. 3. However, the dose distribution may not be drastically changed because it would only improve locally the parameters. Since on average, a clinically acceptable plan is already reached by using four needle insertions, coupling this method with a gradient-based optimization might not be necessary.

Moreover, the pipeline proposed is compatible with a re-optimization of the dose plan parameters after each insertion of the needle thanks to the proposed experimental set-up: in terms of hardware, the robotic device developed in our institution is such that the needle can be inserted needle under MRI guidance, and, in terms of software, the calculation time of the optimizer could be further reduced to make it eligible for intra-operative use. The development of a procedure which re-optimize of the dose plan parameters after each insertion of the needle will be also studied in future work.

Although the optimizer described here was implemented with the point source approximation for simplicity of calculation, the proposed method allows the use of more precise source models such as the line source approximation.

This optimizer workflow completes all objectives: it is developed for divergent needle patterns with a single rotation point and it optimizes the clinically relevant dose parameters of HDR brachytherapy for prostate cancer, specifically D 95% P T V , D 10% U r D 1cc Rec and D 1cc Bl . Unlike most optimizers such as HIPO [START_REF] Karabis | Optimization of Catheter Position and Dwell Time in Prostate HDR Brachytherapy using HIPO and Linear Programming World Congress on Medical Physics and Biomedical Engineering[END_REF] and IPSA [START_REF] Hsu | Comparison of inverse planning simulated annealing and geometrical optimization for prostate high-dose-rate brachytherapy[END_REF], no manual determination of several input parameters as minima or maxima dose for PTV and organs at risk is necessary. Moreover, most optimizers require weight coefficients to be defined (see [START_REF] Dinkla | Comparison of inverse optimization algorithms for HDR/PDR prostate brachytherapy treatment planning Brachytherapy[END_REF], whereas the proposed optimizer workflow does not since an automatic exhaustive search for optimal coefficient values is performed. The results prove that the optimizer workflow presented is able to obtain a clinically accepted plan with a few needle insertions already, whereas during the current clinical procedure, 13 to 17 needles are usually inserted. However, the impact of practical error in needle positioning must be studied to find the optimal number of needle insertions for real clinical procedures. Finally, the proposed optimizer took less than 20 minutes to compute although it has not been optimized in speed. Therefore, it shows great perspectives for a further reduction of the computation time by parallelizing the calculation in order to be eligible for intra-operative use.

Conclusion

In this paper, a complete inverse dose planning optimization workflow for focal-HDR prostate brachytherapy with a divergent needle pattern was presented. It can determine the optimal center of rotation, needle positioning (thus source positions) and dwell times of the sources in order to deliver the desired dose distribution for a given number of needle insertions. The optimization is made on the dose coverage (meaning the D 95% P T V , D 10% U r D 1cc Rec and D 1cc Bl ) and does not need to set importance factors for the organ doses considered as input. Clinically accepted plans were obtained for on average 4 needle insertions for the 10 patients tested.

Figure 1 .

 1 Figure 1. Complete optimization workflow for the determination for the position of the center of rotation, the needle orientations and the dwell times of the source.

Figure 2 .

 2 Figure 2.Figure 2(a) depicts the schematic representation of the exhaustive search of the position of the center of rotation. The area of insertion is represented in light green and the center of rotation evaluated are represented by the black dots. Figure 2(b) describes the proposed heuristics for the determination of needle orientations. First, the PTV and the urethra are projected from the center of rotation onto a transverse plane (cf. schematic on the left). A uniform distribution of needles is then found by applying k-means clustering on P P T V until no clusters are found in the urethra zone (cf. schematic on the right).

Figure 3 .

 3 Figure 3. Results of the simulation obtained on a typical patient with 10 needle insertions. The values of D 95% P T V , D 10% U r , D 1cc Rec , D 1cc Bl and E (computed from Eq. 10) are presented as function of log 10 (ω P T V ) and D P T V (Figure 3(a), 3(b), 3(c), 3(d) and 3(e) respectively). log 10 (ω tissue ) was deducted from Eq. 17. On figure 3(a), 3(b), 3(c) and 3(d), the horizontal green lines represent the clinical constraints set as input.

Figure 4 .

 4 Figure 4. Energy of the selected dose plan averaged on ten patients as a function of the number of heuristics evaluated. The total computation time of the optimizer (linear with the number of heuristics evaluated) is depicted in red.

Figure 5 .

 5 Figure 5. Results of the planning study for using 2, 4, 6, 8, 10, 12 and 14 needle insertions. Figure 5(a)(b)(c)(d) and (e) present the optimization results for a typical patient case. Figure 5(a) shows the output parameters (D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl ).Figure 5(b) shows the Homogeneity (HI 150% and HI 200% ) and the conformal index (COIN ).Figure 5(c) presents a slice of MRI image in the sagittal plane with the delineations of the volumes of interest. Figures 5(d) and (e) present the corresponding dose distributions for 2 and 10 needle insertions respectively on the same sagittal slice. Figure 5(f) and (g) present the optimization results for 10 different patients. Figure 5(f) shows, the average output parameters on all patients D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl .Figure 5(g) represents the average on all patients of homogeneity and conformal index (HI 150% , HI 200% and COIN respectively). For the two latter graphs, the whole range of values is represented by the solid vertical line. For Figure 5(a)(b)(f) and (g), The red solid lines represent the clinical constraints set as input (D min 95% P T V , D max 10% U r , D max 1cc Rec and D max 1cc Bl respectively)

  Figure 5. Results of the planning study for using 2, 4, 6, 8, 10, 12 and 14 needle insertions. Figure 5(a)(b)(c)(d) and (e) present the optimization results for a typical patient case. Figure 5(a) shows the output parameters (D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl ).Figure 5(b) shows the Homogeneity (HI 150% and HI 200% ) and the conformal index (COIN ).Figure 5(c) presents a slice of MRI image in the sagittal plane with the delineations of the volumes of interest. Figures 5(d) and (e) present the corresponding dose distributions for 2 and 10 needle insertions respectively on the same sagittal slice. Figure 5(f) and (g) present the optimization results for 10 different patients. Figure 5(f) shows, the average output parameters on all patients D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl .Figure 5(g) represents the average on all patients of homogeneity and conformal index (HI 150% , HI 200% and COIN respectively). For the two latter graphs, the whole range of values is represented by the solid vertical line. For Figure 5(a)(b)(f) and (g), The red solid lines represent the clinical constraints set as input (D min 95% P T V , D max 10% U r , D max 1cc Rec and D max 1cc Bl respectively)

  Figure 5. Results of the planning study for using 2, 4, 6, 8, 10, 12 and 14 needle insertions. Figure 5(a)(b)(c)(d) and (e) present the optimization results for a typical patient case. Figure 5(a) shows the output parameters (D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl ).Figure 5(b) shows the Homogeneity (HI 150% and HI 200% ) and the conformal index (COIN ).Figure 5(c) presents a slice of MRI image in the sagittal plane with the delineations of the volumes of interest. Figures 5(d) and (e) present the corresponding dose distributions for 2 and 10 needle insertions respectively on the same sagittal slice. Figure 5(f) and (g) present the optimization results for 10 different patients. Figure 5(f) shows, the average output parameters on all patients D 95% P T V , D 10% U r , D 1cc Rec and D 1cc Bl .Figure 5(g) represents the average on all patients of homogeneity and conformal index (HI 150% , HI 200% and COIN respectively). For the two latter graphs, the whole range of values is represented by the solid vertical line. For Figure 5(a)(b)(f) and (g), The red solid lines represent the clinical constraints set as input (D min 95% P T V , D max 10% U r , D max 1cc Rec and D max 1cc Bl respectively)

  ) and 5(f) confirms the observation shown in Figures 5(a) and 5(b) for a typical patient. D 95% P T V increased with the number of needle insertions, whereas D 10% U r , D 1cc Rec and D 1cc Bl decreased. The large range between the maximum and minimum values is explained by the different anatomy of each patient especially the position of the tumor with respect to the other organs at risk. HI 150% , HI 200% and COIN showed the same trend as in Figures 5(a) and 5(b).

Acknowledgements

This study was funded by Philips Medical Systems Nederland B.V. The authors also thank the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma).