Optimizing robot motion for robotic ultrasound-guided radiation therapy

Stereotatic body radiation therapy (SBRT) is a non-invasive approach to cancer treatment. It uses dose distributions which are designed to cover the target structures (planning target volume, PTV) very well while at the same time sparing other critical structures (organs at risk, OARs). In the case of prostate cancer, the intention is to reduce the dose in bladder and rectum, for example.

The flexibility of robotic radiation therapy to deliver beams from a wide range of positions around the patient can help shaping very conformal dose distributions with steep gradients towards critical structures (Hossain et al 2010, Schlaefer et al 2013). Another important and interesting aspect of robotic radiation therapy is active compensation of target motion. By detecting the motion of the target and moving the treatment beams accordingly, the adverse effect of organ motion can be mitigated and hence the safety margins accounting for organ motion can be reduced (Beltran et al 2008, Ricotti et al 2018). So far, approaches based on x-ray imaging (Krauss et al 2012, Nguyen et al 2017) and surrogate markers (Dürichen et al 2014, Antoni et al 2016) have dominated. Typically, the target region is defined by fiducial markers implanted before treatment. For some targets in the thorax, direct localization is feasible (Jung et al 2015, Yang et al 2017). However, given its ionizing nature, x-ray imaging is typically not used for continuous motion tracking. Moreover, this approach does not extend to targets in the abdomen, which are often not clearly visible in x-ray images.

In contrast, volumetric ultrasound (US) presents a non-ionizing image modality with high spatial and temporal resolution. Recently, a number of approaches to use ultrasound for motion compensated treatments have been proposed (Schwaab et al 2014, Western et al 2015, Schlosser and Hristov 2016, Camps et al 2018). Ipsen et al (2016) demonstrated 4D tracking for MLC-based motion compensation by acquiring US volumes with $0.3 \times 0.8 \times 1.6$ mm³ resolution and sizes from $17 \times 14 \times 45$ mm³ at 250 Hz volume rate to $68 \times 105 \times 170$ mm³ at 17.8 Hz. Note that the dimension of the volumes will cover typical target motion, e.g. of the prostate. Moreover, even if motion exceeding this range occured, the tracking would detect this and similar strategies as in current motion compensated treatments could be employed, e.g. emergency stop and re-setup. Tracking algorithms for both 2D and 3D ultrasound images have been evaluated and showed that this modality can be suitable for tracking of targets like the liver or the prostate (O'Shea et al 2014, De Luca et al 2018).

While it has been demonstrated that 4D ultrasound can be a viable alternative to obtain real-time images of the abdominal anatomy, a remaining challenge is its integration with the treatment systems. Particularly, the ultrasound probe needs to be positioned at the patient such that high quality imaging of the target region is feasible throughout the treatment. This can be realized by using a robotic arm holding the ulrasound probe (Priester et al 2013, Bell et al 2014, Şen et al 2017). Different approaches have been considered to integrate ultrasound imaging with radiation therapy, where beams should not pass through the ultrasound probe or the robot and therefore some beams may be effectively blocked. Elekta proposed a setup where the transducer is placed at the perineum (Abramowitz et al 2012, Li et al 2015). However, this is primarily limited to imaging the prostate. Also the use of specialized robot systems has been studied (Schlosser et al 2016). Nevertheless, there will always be a residual blocking of beam directions by the probe itself (Bazalova-Carter et al 2015), because ultrasound imaging requires direct contact of the probe with the patient close to the target structure. For this reason, a special radiolucent ultrasound system has been proposed which also aims to reduce distortion of computed tomography (CT) images while the probe is present (Schlosser and Hristov 2016). The latter aspect is important because the pressure of the probe induces deformations and a planning CT without the probe does not reflect the actual geometry during US imaging. Therefore, workflows for reproducible probe placement have been developed (Bell et al 2014, Şen et al 2017).

It has also been demonstrated that the freedom in beam delivery of the robotic CyberKnife (Accuray Inc., United States) combined with an optimized configuration of a robot carrying the ultrasound probe will largely compensate for possible adverse effects to the treatment plan quality (Gerlach et al 2017a, 2017b). Considering the expected beam blocking during treatment planning it is possible to determine a set of treatment beams and a configuration for the robot carrying the ultrasound probe such that the change in plan quality compared to a setup without ultrasound is very small. Recently, also automatic approaches to determine a robot setup with minimum impact on the plan quality have been proposed (Schlüter et al 2019). Particularly, employing a kinematically redundant robot like the light-weight robot iiwa (KUKA Roboter GmbH, Germany) to carry the ultrasound probe can further reduce the beam blocking. Its seven degrees of freedom (DoF) allow to continuously change the orientation of the robot's elbow while maintaining the same pose of the end-effector. Thereby, the robot is able to maintain ultrasound transducer contact and provide steady imaging throughout the beam delivery, even when it changes its configuration. This is interesting, as beams blocked by the elbow in one configuration may not be blocked in another configuration and changing configurations would not affect the ultrasound imaging. It has been shown that this extra flexibility can indeed improve the dosimetrical plan quality further (Gerlach et al 2017a, Schlüter et al 2019) and it has been demonstrated how the calibrations between the different systems and coordinate frames can be realized (Gerlach et al 2017a). For the inverse kinematic, this elbow rotation due to the additional degree of freedom can be parametrized by a continuous angle we call LIFT (Kuhlemann et al 2016).

However, while the dose distribution is the key determinant for plan quality, delivery time contributes to treatment quality, workflow, and patient experience. The potential need to wait for the ultrasound robot to change its configuration before continuing beam delivery from a certain direction will take time and prolong the treatment. The problem to optimize the sequence of beam source motions during beam delivery can be modeled as a traveling salesman problem (TSP) (Kearney et al 2017). While in general the TSP is known to be hard, we first demonstrate that the particular instances resulting from the CyberKnife setup can be efficiently solved with state-of-the-art methods. Yet, when including the configuration changes of the ultrasound robot, the problem becomes more difficult. Second, we derive a generalized TSP (GTSP) model for the two-robot setup. We compare two approaches to solve this problem. Finally, we provide an analysis of the effect of the optimization for different assumptions regarding the robots' velocity settings.

In this section, we first describe our general treatment planning setup. We then introduce the robots considered in our US-guided robotic radiation-therapy scenario, namely the robot carrying the linear accelerator (LINAC-robot) and the robot carrying the ultrasound probe (US-robot) as sketched in figure 1(a). Furthermore, we develop models to coordinate and optimize robot motion and beam delivery. We describe two approaches to solve the underlying optimization problems and present a case study to evaluate our methods.

Zoom In Zoom Out Reset image size

2.1. Treatment plan optimization

2.2. Robot parameters

Any motion of the robots carrying the linear accelerator and the ultrasound probe is limited by their respective joint speeds. We consider kinematic parameters of a large articulated arm (KR 300 R2500 ultra, KUKA Roboter GmbH, Germany) for the LINAC-robot. This robot represents a standard serial kinematic chain with six DoF. For the US-robot we consider a light-weight device with seven DoFs (LBR iiwa, KUKA Roboter GmbH, Germany). The additional joint results in a kinematic redundancy allowing for rotation of the elbow while maintaining the end-effector pose. While our methods work for any ultrasound-probe shape, we consider an actual 4D ultrasound probe (GE 3V, General Electric, USA), which is a 2D matrix array probe, for our experiments. The setup with both robots is illustrated in figure 1(a). The maximum joint velocities of the KR range from 101 to 122 deg s⁻¹ for the first five joints and is 175 deg s⁻¹ for the last. For the LBR iiwa, the first three joints are limited to 98–100 deg s⁻¹, the fourth to 130 deg s⁻¹, the fifth to 140 deg s⁻¹, and the last two joints to 180 deg s⁻¹. However, considering the payload and safety concerns in a clinical setting, much lower velocities will be used in practice (McGuinness et al 2015). We study velocities ranging from 1% to 10% of the maximum values.

Following the derivation of an analytic inverse kinematic by (Kuhlemann et al 2016), different configurations of the US-robot can be parametrized by a continuous angle which we call LIFT. By varying this LIFT angle during the treatment, it is possible to reduce the number of blocked beams, i.e. beams blocked in one elbow configurations may be unblocked in another elbow configuration.

2.2.1. Beam delivery without ultrasound robot

After solving the inverse planning problem, we obtain a set B of N treatment beams b_i, $i=1,...,N$ with non-zero weight which have to be delivered by the linear accelerator. Hence, the LINAC-robot has to move to all nodes where treatment beams start. Multiple beams may start at each node and small changes in the beam orientation also require some robot motion. The order in which the CyberKnife delivers the beams determines the time required for robot motion. For example, it is typically less efficient to visit a different node after each delivered beam compared to first delivering all beams from one node before moving to the next node. To minimize the overall treatment time, we are therefore interested in a sequence of beams which minimizes the total robot motion time.

The problem of finding such a sequence of beams can be modeled as a TSP. Consider a complete graph $G=(V,E)$ with a set of vertices $V$ which are connected by a set of edges $E = \{(v_i,v_j):v_i,v_j \in V\}$. We assign a non-negative weight

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d(v_i,v_j)=d_{ij} \geqslant 0 \nonumber \end{align} \tag{ 1 }$$

to each edge $(v_i,v_j) \in E$, which represents the distance or costs for traveling between the two vertices. We use symmetric weights, i.e. $d_{ij} = d_{ji}$ for all $i,j$, and further set d_ii  =  0. A path in such a graph is a sequence of distinct vertices. If only the first and last vertex in a sequence are identical, it is called a cycle. The TSP now asks for an optimal tour, which is a cycle containing each vertex of the graph. Optimality is defined in terms of minimizing the sum of costs d_ij of the involved edges in the tour, i.e. the edges which directly connect the sequent vertices. Given N vertices, there are in total not $N!$ but $(N-1)!/2$ different tours, because a tour does neither have a defined start (factor of N) nor a defined direction in the case of symmetric costs (factor of 2). For example, the tours described by $(v_1,v_2,v_3,v_1)$, $(v_2,v_3,v_1,v_2)$, and $(v_1,v_3,v_2,v_1)$ are all equivalent.

In our scenario, each treatment beam b_i is related to a vertex such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle V_{\rm TSP} = \left\{b_i: i \in \{1,...,N\} \right\} \nonumber \end{align} \tag{ 2 }$$

is the set of vertices to be visited and the cost for traveling between two vertices is the time needed for changing the end-effector pose of the LINAC-robot accordingly. We estimate this time as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e:tla} d_{ij} = t_{\rm LA}(b_i,b_j) = \max_{g \in \{1,...,6\}} \frac{|\theta_g^i - \theta_g^{\,j}|}{\lambda_{\rm LA} s_g^{\rm LA}}, \nonumber \end{align} \tag{ 3 }$$

where $\theta_g^i$ is the angle of the gth joint, which is calculated from inverse kinematics for realizing beam b_i. The differences between the angles are calculated by using the identity $\theta=\theta+n \cdot 2\pi$, $n\in \mathbb{Z}$, and are wrapped to the half-open interval $[-\pi,\pi)$. The denominator is the velocity of the joint, which we assume to be a fraction $\lambda_{\rm LA} \in (0,1]$ of its maximum velocity $s_g^{\rm LA}$. Note that the total treatment time would include beam delivery. However, as the beam weights result from inverse planning they are fixed. Also, the time needed to change the effective diameter of the IRIS collimator is neglected as it is not related to the robot motion and could be realized concurrently.

The solution of a TSP is an optimal tour, i.e. the first and the last visited vertex are equal. In terms of our radiotherapy scenario, this means that the LINAC-robot eventually moves back to its initial pose although the related beam has already been delivered. This is neither necessary nor intuitve. Thus, we search for a path containing all vertices exactly once rather than for a cycle. Such a path is called a Hamiltonian path resulting in a shortest Hamiltonian path problem (HPP). Yet, we can transform the HPP into an equivalent TSP by adding an artificial vertex connected to all other vertices with zero costs. In the resulting optimal tour, the successor and predecessor of this artificial node are the first and last node in the optimal Hamiltonian path, respectively. For example, an optimal tour $(v_1,v_2,...,v_5,v^a,v_6,...,v_N,v_1)$ with actual vertices $v_i$, $i=1,...,N$, and artificial vertex $v^a$ is converted to the optimal Hamiltonian path $(v_6,...,v_N,v_1,v_2,...,v_5)$ only containing the N actual vertices. The size of the solution space is $N!/2$ for an HPP with N vertices and thus a factor of N larger than for a TSP, because it has a defined start. The direction is still not relevant due to the symmetric costs. For simplicity, we will subsequently describe the TSP models, but all results are computed for the associated HPP.

2.2.2. Beam delivery including the robotic ultrasound

The kinematic redundancy of the US-robot allows delivering beams by moving the elbow such that the beams are no longer blocked. Such a motion will take time, i.e. the problem to optimize the treatment delivery needs to include the motion of both the LINAC-robot and the US-robot. Note that both robots may move concurrently. We consider M configurations c_k, $k=1,...,M$, of the US-robot such that each beam is feasible, i.e. not blocked, in at least one configuration. We map each beam b_j to its set F_j of feasible configurations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F(b_j) = F_j = \left\{c_k : b_j ~{\rm feasible~in}~ c_k ~{\rm and}~ k\in \{1,...,M\} \right\}. \nonumber \end{align} \tag{ 4 }$$

We create a new model where for each beam b_j all pairs $(b_j, c_k)$, $c_k \in F_j$, are added as vertices. Therefore, the number of vertices increases from N to up to NM. Assuming that both robots do not change their configurations during beam delivery we model the cost for robot motion by extending equation (3) to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{e:tlaus} \bar{d}_{ijkl} = t \left( (b_i,c_k),(b_j,c_l) \right) = \max \left\{t_{\rm LA}(b_i,b_j) , t_{\rm US}(c_k,c_l) \right\}, \nonumber \end{align} \tag{ 5 }$$

where $t_{\rm US}(c_k,c_l)$ is calculated similarly to $t_{\rm LA}(b_i,b_j)$ but with a 7th joint as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle t_{\rm US}(c_k,c_l) = \max_{g \in \{1,...,7\}} \frac{|\theta_g^k - \theta_g^l|}{\lambda_{\rm US} s_g^{\rm US}}. \nonumber \end{align} \tag{ 6 }$$

Solving this problem as an ordinary TSP would result in a sequence including all pairs $(b_j,c_k)$, i.e. each feasible US-robot configuration would be considered for each beam. To identify a sequence containing each beam just once we need to solve the generalized TSP, also known as set TSP or group TSP (Noon and Bean 1993). Here, the vertices form exclusive groups and a tour has to contain exactly one vertex of each group. The standard TSP can be seen as a special instance of the generalized TSP, in which the size of each group is one. See figure 2 for a comparison of the two models. In our model, there is one group G_j per beam b_j and it contains one element per configuration in which b_j is feasible, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_j = \left\{(b_j,c_k): c_k \in F_j \right\}. \nonumber \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Thus, we have N groups with at least one and at most M vertices each which form the overall set of vertices

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle V_{\rm GTSP} = \bigcup_{j \in \{1,...,N\}} G_j \nonumber \end{align} \tag{ 8 }$$

with edges

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\rm GTSP} = \left\{\left( (b_i,c_k),(b_j,c_l) \right) : (b_i,c_k) \in G_i, (b_j,c_l) \in G_j , i \neq j \right\} \nonumber \end{align} \tag{ 9 }$$

of our GTSP model. The size of the solution space depends on the sizes of the individual groups. The $(N-1)!/2$ possible tours in terms of visiting each group in a sequence $(G^1,G^2,...,G^N,G^1)$ have to be multiplied by the product $\prod_j |G_j|$ of the number of vertices of each group. In the worst case with all groups containing M vertices, we end up with $(N-1)!/2 \cdot M^N$ possible solutions.

While the GTSP provides a suitable model for our problem, finding its solution is typically hard. There exist algorithms to transform any GTSP to a symmetric standard TSP (Noon and Bean 1993) of similar size. However, this introduces very large costs on some of the edges and results in a TSP with structures which hamper state-of-the-art solvers. Some branch & cut algorithms for solving GTSPs exist (Fischetti et al 1997), but they are limited to to very small problem instances. Hence, several heuristic solvers have been introduced, which do not guarantee optimality but are able to find good approximations in reasonable time. The GLKH solver (Helsgaun 2013) based on the Lin–Kernighan–Helsgaun heuristic (Lin and Kernighan 1973, Helsgaun 2009) has been shown to work well on GTPS problems (Helsgaun 2015). Other heuristic approaches include methods based on genetic algorithms (Gutin and Karapetyan 2010), ant-colony search (Jun-man and Yi 2012), or particle-swarm optimization (Shi et al 2007). Subsequently we present three methods to solve the GTPS including the US-robot. For comparison we also solve the TSP optimizing just the LINAC-robot motion without adding the US-robot.

2.2.3. Method A: LINAC-robot motion only

As a baseline we solve the TSP resulting from transforming the HPP to optimize the path of the LINAC-robot using Concorde (Applegate et al 2003). This solver is based on a mixed-integer programming (MIP) approach and employs efficient branch & cut strategies for TSP. Hence, the solver is complete and optimal and known to be very efficient, even in case of relatively large TSP instances. Solving the TSP results in a sequence S_A of beams which encodes the motion of the LINAC-robot.

2.2.4. Method B: decoupled, no-optimization

The solution to the TSP from method A can be extended in a GTPS solution in a straightforward fashion. We keep the same sequence S_A and start with an arbitrary US-robot configuration such that the first beam is feasible. Subsequently, the beams in the sequence are delivered until a beam is not feasible given the current US-robot configuration. In this case, the closest US-robot configuration is selected and the cost for the motion of both robots is considered.

As a post-processing step, we shift the resulting sequence S_B such that the last and first vertex are connected by the edge with highest weight. This weight does not contribute to the final costs as we are interested in a Hamiltonian path and not a tour.

2.2.5. Method C: decoupled, optimization

Although method B will always result in a feasible sequence of robot motion it may incur unnecessary cost as the configuration of the US-robot is considered at each beam delivery separately. This problem is tackled by our method C. In this approach, we again keep the beam sequence $S_{\rm A}$ from the problem without the US-robot fix. However, we now optimize the sequence of configurations. We do this by setting up a directed graph with artificial source and sink vertex as shown in figure 3. The source is connected by zero-weighted edges to all vertices containing the first beam in $S_{\rm A}$. Likewise, all vertices containing the last beam in $S_{\rm A}$ are connected to the sink with zero costs. We now search for the shortest path from source to sink, which can be done by using Dijkstra's algorithm (Dijkstra 1959). This path is our resulting sequence S_B of beams and configurations. We apply the same post-processing as for method A, in order to remove the edge with largest cost from the costs of the sequence S_C.

Zoom In Zoom Out Reset image size

2.2.6. Method D: combined, optimization

To obtain an estimate to what extent optimizing the full GTSP can improve the treatment time, we solve the problem using GLKH. We use its standard parameters and apply one run with 1000 iterations. For comparison, we also evaluate the influence of more iterations and restarts. As for all other methods, the highest edge of the tour is removed afterwards.

2.3. Case study

To illustrate the methods, we consider a prostate cancer scenario for which a previous plan quality analysis has shown that the effect of adding the US-robot can be mitigated by using multiple configurations (Gerlach et al 2017a, 2017b). Using the same constraints for inverse planning but different parameters for the robot configurations we now study to what extent the proposed methods can also mitigate the US-robot's effect with respect to treatment time. We apply the methods to 10 patient cases previously treated with the CyberKnife and having PTV sizes ranging from 52.7 cm³ to 115.4 cm³ with mean 84.5 cm³ and standard deviation 13.8 cm³.

2.3.1. Robot configurations

For each patient, a promising ultrasound probe pose was determined based on the results of the previous planning studies. Moreover, three different robot base positions around the patients were selected and are shown in figure 4(a). One robot base is located between the patients' legs and allows to vary the LIFT angle collision-freely within the interval $[-69^{\circ},69^{\circ}]$. For the two positions left of the patient, the intervals are $[-90^{\circ},90^{\circ}]$ and $[-115^{\circ},115^{\circ}]$. During planning we consider using 3, 5, and 7 LIFT angles, which we select equidistantly in these intervals. For the base position with the widest range, the seven LIFT angles in use are illustrated in figure 4(b).

Zoom In Zoom Out Reset image size

To investigate the effect of robot speed on treatment time, we assume 1%, 3%, 5%, and 10% for $\lambda_{\rm LA}$ and 1%, 4%, 7%, and 10% for $\lambda_{\rm US}$. In the actual CyberKnife system, the linear accelerator is moved with approximately 3%–4% of the maximum speed in Cartesian space.

2.3.2. Treatment planning parameters

In our study, we are mainly interested in the effect different robot motion patterns have on treatment time. We adopt the same bounds as in the previous study (Gerlach et al 2017a), constraining the dose in PTV, rectum, and bladder to 4050 cGy, 3600 cGy, and 3600 cGy, respectively. SHELL structures are used to maintain the same level of conformality with patient-depending upper bounds from 2700 to 2975 cGy and from 3525 to 3660 cGy, respectively. The PTV has a target dose of 3625 cGy. For dose delivery, we assume a rate of 800 MU min⁻¹. The upper MU limit is set to 300 MU per beam and 40 000 MU in total. For each combination of patient and robot setting, 50 random sets of 6000 candidate beams are generated and used for treatment planning. In total, the delivery of 4500 different sets of candidate beams is analyzed for 12 different combinations of robot speeds.

The objective for the optimization is minimization of underdosage. This is formulated as a linear function by summing the difference to target dose over all voxels with a current dose below the target dose. Note that solving the optimization problem will result in many of the candidate beams having zero weight, i.e. they are not used for treatment. Only the beams having a positive weight are subsequently considered for optimizing the robot paths.

Using multiple US-robot configurations during the treatment allows for more beam directions, because each beam only has to be feasible in at least one of the configurations. Table 1 shows that when 3 LIFT angles are considered, on average 4.1% of the sampled candidate beams are only feasible in one of the three configurations with an additional 16.1% being feasible in two configurations. The remaining 79.8% are feasible in all three configurations. Similarly, 77.3% and 76.6% are feasible in all configurations when 5 and 7 LIFT angles are used. If we only consider the treatment beams instead, the fraction of beams which are feasible in all configurations drops by about 10 percentage points to 67.7% and 66.7% for 5 and 7 LIFTs, respectively. Thus, beams which are not always feasible are disproportionately selected by the optimizer. The same holds for the case of 3 LIFT angles. If only the three patients with the smallest prostates are considered, the average percentage of treatment beams visible in all configurations are 69.5%, 66.7%, and 66.1% for 3, 5, and 7 LIFT angles, respectively. For the three patients with the largest prostates, these values are 71.6%, 69.5%, and 68.2%, respectively.

Table 1. Distribution of the number of US-robot configurations in which the beams are feasible for different numbers of LIFT angles. Columns labeled with treatment only consider beams having a non-zero weight after dose optimization, columns labeled with candidates consider all 6000 candidate beams.

	Feasible beams in % (mean and std)
	3 LIFTs		5 LIFTs		7 LIFTs
# Feas	Candidates	Treatment	Candidates	Treatment	Candidates	Treatment
1	$4.1 \pm 1.3$	$9.1 \pm 3.7$	$1.0 \pm 0.6$	$2.0 \pm 1.1$	$0.6 \pm 0.3$	$1.2 \pm 0.9$
2	$16.1 \pm 7.1$	$20.6 \pm 4.4$	$3.5 \pm 1.1$	$7.5 \pm 2.7$	$1.4 \pm 0.8$	$3.1 \pm 2.0$
3	$79.8 \pm 7.0$	$70.2 \pm 3.3$	$6.5 \pm 3.0$	$10.0 \pm 2.3$	$2.8 \pm 0.8$	$5.8 \pm 1.6$
4			$11.7 \pm 5.8$	$12.8 \pm 3.5$	$3.9 \pm 1.6$	$6.9 \pm 1.9$
5			$77.3 \pm 9.3$	$67.7 \pm 4.2$	$6.9 \pm 3.9$	$8.7 \pm 2.9$
6					$7.8 \pm 3.6$	$7.6 \pm 2.2$
7					$76.6 \pm 8.9$	$66.7 \pm 3.6$

Figure 5 summarizes the results of the robot path optimization without consideration of the US-robot (method A). The median fraction of treatment time, i.e. time for beam delivery and the motion of linac, is 11.6% for $\lambda_{\rm LA}=3$% and 3.8% for $\lambda_{\rm LA}=10$%. For the slow linac with $\lambda_{\rm LA}=1$% it increases to 28.3%. For comparison, results for a random trajectory are shown. Our simulated velocities are compared in figure 6 with velocities of the manipulator which were logged during two real CyberKnife treatments. Only traveling distances above 100 mm are considered in the histograms for both our simulation and the real data. This excludes situations in which not the maximum velocity but rather the feasible acceleration dominates the traveling time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering the US-robot, generally its motion is related to an increase in the overall treatment time, i.e. time for beam delivery and robot motion. The results are presented in figure 7 for different combinations of robot speeds and different numbers of LIFT angles in use. The 100% reference is the overall time by method A, i.e. representing motion of the LINAC-robot only. The results for methods B, C, and D move towards 100% if the US-robot is considerably faster than the LINAC-robot, i.e. the overhead due to the US-robot vanishes. For the other extreme, i.e. when the LINAC-robot moves about ten times faster than the US-robot, the median increase in treatment time is approximately 20% and 10% for method B and D, respectively. For the realistic case of LINAC-robot and US-robot operating at 3% and 4% of their maximum speeds, respectively, the median increase is about 4%. More detailed results are shown in table 2, which compares only the times needed for moving the robots. The time for coordinated linac and US-robot motion by methods B, C, and D is expressed relative to the reference motion by method A. Method B is outperformed by methods C and D. Whether method C or D produces better results depends on the robot speed settings. The differences are mostly significant in terms of a Wilcoxon rank sum test. However, the variation of the results for method D is typically smaller than for method C.

Zoom In Zoom Out Reset image size

Table 2. Mean $\mu$ and standard deviation $\sigma$ of the time for robot motion by method B, C, and D (linac and US-robot motion) relative to method A (linac-only motion) for different combinations of LIFT angles and robot speeds. A Wilcoxon rank sum test is used to check the null hypothesis that two methods' results have equal median and the p-value is marked bold if it cannot be rejected on a 5% level.

#LIFTs	$\lambda_{\rm LA}$	$\lambda_{\rm US}$	Method B	Method C	Method D	p B,C	p B,D	p C,D
3	1%	1%	$\mu=1.60, \ \sigma=0.25$	$\mu=1.42, \ \sigma=0.20$	$\mu=1.42, \ \sigma=0.10$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.13, \ \sigma=0.06$	$\mu=1.08, \ \sigma=0.04$	$\mu=1.12, \ \sigma=0.07$	${0.00}$	${\bf 0.21}$	${0.00}$
		10%	$\mu=1.03, \ \sigma=0.02$	$\mu=1.01, \ \sigma=0.01$	$\mu=1.01, \ \sigma=0.01$	${0.00}$	${0.00}$	${\bf 0.50}$
	3%	1%	$\mu=2.79, \ \sigma=0.78$	$\mu=2.28, \ \sigma=0.62$	$\mu=1.81, \ \sigma=0.20$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.44, \ \sigma=0.19$	$\mu=1.31, \ \sigma=0.15$	$\mu=1.37, \ \sigma=0.09$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.16, \ \sigma=0.07$	$\mu=1.10, \ \sigma=0.05$	$\mu=1.18, \ \sigma=0.08$	${0.00}$	${0.00}$	${0.00}$
	5%	1%	$\mu=4.00, \ \sigma=1.31$	$\mu=3.14, \ \sigma=1.04$	$\mu=2.20, \ \sigma=0.30$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.74, \ \sigma=0.32$	$\mu=1.52, \ \sigma=0.25$	$\mu=1.47, \ \sigma=0.12$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.28, \ \sigma=0.12$	$\mu=1.19, \ \sigma=0.09$	$\mu=1.31, \ \sigma=0.09$	${0.00}$	${0.00}$	${0.00}$
	10%	1%	$\mu=7.02, \ \sigma=2.67$	$\mu=5.30, \ \sigma=2.10$	$\mu=3.16, \ \sigma=0.56$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=2.49, \ \sigma=0.65$	$\mu=2.07, \ \sigma=0.52$	$\mu=1.71, \ \sigma=0.18$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.59, \ \sigma=0.25$	$\mu=1.42, \ \sigma=0.20$	$\mu=1.43, \ \sigma=0.10$	${0.00}$	${0.00}$	${0.00}$
5	1%	1%	$\mu=1.57, \ \sigma=0.27$	$\mu=1.40, \ \sigma=0.20$	$\mu=1.44, \ \sigma=0.12$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.11, \ \sigma=0.06$	$\mu=1.07, \ \sigma=0.04$	$\mu=1.11, \ \sigma=0.06$	${0.00}$	${\bf 0.75}$	${0.00}$
		10%	$\mu=1.03, \ \sigma=0.02$	$\mu=1.01, \ \sigma=0.01$	$\mu=1.01, \ \sigma=0.01$	${0.00}$	${0.00}$	${0.00}$
	3%	1%	$\mu=2.72, \ \sigma=0.83$	$\mu=2.23, \ \sigma=0.63$	$\mu=1.89, \ \sigma=0.21$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.42, \ \sigma=0.20$	$\mu=1.29, \ \sigma=0.15$	$\mu=1.38, \ \sigma=0.11$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.15, \ \sigma=0.07$	$\mu=1.09, \ \sigma=0.05$	$\mu=1.16, \ \sigma=0.08$	${0.00}$	${0.00}$	${0.00}$
	5%	1%	$\mu=3.87, \ \sigma=1.41$	$\mu=3.07, \ \sigma=1.07$	$\mu=2.27, \ \sigma=0.32$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.71, \ \sigma=0.33$	$\mu=1.50, \ \sigma=0.25$	$\mu=1.53, \ \sigma=0.13$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.27, \ \sigma=0.13$	$\mu=1.18, \ \sigma=0.09$	$\mu=1.28, \ \sigma=0.09$	${0.00}$	${0.00}$	${0.00}$
	10%	1%	$\mu=6.79, \ \sigma=2.78$	$\mu=5.16, \ \sigma=2.16$	$\mu=3.20, \ \sigma=0.62$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=2.43, \ \sigma=0.68$	$\mu=2.02, \ \sigma=0.53$	$\mu=1.80, \ \sigma=0.19$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.57, \ \sigma=0.27$	$\mu=1.40, \ \sigma=0.20$	$\mu=1.46, \ \sigma=0.12$	${0.00}$	${0.00}$	${0.00}$
7	1%	1%	$\mu=1.51, \ \sigma=0.22$	$\mu=1.32, \ \sigma=0.18$	$\mu=1.45, \ \sigma=0.11$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.10, \ \sigma=0.05$	$\mu=1.05, \ \sigma=0.03$	$\mu=1.12, \ \sigma=0.06$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.03, \ \sigma=0.01$	$\mu=1.01, \ \sigma=0.01$	$\mu=1.01, \ \sigma=0.01$	${0.00}$	${0.00}$	${0.00}$
	3%	1%	$\mu=2.56, \ \sigma=0.67$	$\mu=2.00, \ \sigma=0.56$	$\mu=1.99, \ \sigma=0.21$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.38, \ \sigma=0.17$	$\mu=1.23, \ \sigma=0.13$	$\mu=1.37, \ \sigma=0.10$	${0.00}$	${0.03}$	${0.00}$
		10%	$\mu=1.13, \ \sigma=0.06$	$\mu=1.07, \ \sigma=0.04$	$\mu=1.16, \ \sigma=0.07$	${0.00}$	${0.00}$	${0.00}$
	5%	1%	$\mu=3.64, \ \sigma=1.14$	$\mu=2.67, \ \sigma=0.95$	$\mu=2.44, \ \sigma=0.32$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=1.64, \ \sigma=0.28$	$\mu=1.40, \ \sigma=0.22$	$\mu=1.54, \ \sigma=0.12$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.24, \ \sigma=0.11$	$\mu=1.14, \ \sigma=0.08$	$\mu=1.28, \ \sigma=0.08$	${0.00}$	${0.00}$	${0.00}$
	10%	1%	$\mu=6.31, \ \sigma=2.31$	$\mu=4.39, \ \sigma=1.92$	$\mu=3.50, \ \sigma=0.58$	${0.00}$	${0.00}$	${0.00}$
		4%	$\mu=2.30, \ \sigma=0.57$	$\mu=1.83, \ \sigma=0.47$	$\mu=1.88, \ \sigma=0.19$	${0.00}$	${0.00}$	${0.00}$
		10%	$\mu=1.51, \ \sigma=0.22$	$\mu=1.32, \ \sigma=0.18$	$\mu=1.47, \ \sigma=0.11$	${0.00}$	${0.00}$	${0.00}$

Generally, the size of the PTV can have an impact on the treatment plan quality. To analyze the influence on our study, we grouped the three patients with the smallest PTVs (52.7, 67.1, and 70.3 cm³) and the three patients with the largest PTVs (93.2, 109.7, and 115.4 cm³). We use a Wilcoxon rank sum test to test the null hypothesis that the results for both groups are from distributions with equal median. As shown in table 3, significant differences are observed for most results of method D. For 31 of the 36 presented parameter combinations this hypothesis can be rejected on a 5% significance level.

Table 3. Time for robot motion of method D (LINAC-robot and US-robot) relative to method A (only LINAC-robot) for the three patients with the smallest prostate volumes and the three patients with the largest prostate volumes. Additionally, the p-values for the null hypothesis that both results are samples from distributions with equal median (Wilcoxon rank sum test) are given for each combination of parameters. Values above 5% are marked bold.

#LIFTs	$\lambda_{\rm LA}$	$\lambda_{\rm US}$	small (3)	p	large (3)
3	1%	1%	$\mu=1.394, \ \sigma=0.111$	${0.03}$	$\mu=1.411, \ \sigma=0.095$
		4%	$\mu=1.084, \ \sigma=0.050$	${0.00}$	$\mu=1.138, \ \sigma=0.070$
		10%	$\mu=1.010, \ \sigma=0.006$	${0.00}$	$\mu=1.014, \ \sigma=0.010$
	3%	1%	$\mu=1.768, \ \sigma=0.220$	${\bf 0.09}$	$\mu=1.807, \ \sigma=0.191$
		4%	$\mu=1.344, \ \sigma=0.103$	${0.00}$	$\mu=1.370, \ \sigma=0.088$
		10%	$\mu=1.136, \ \sigma=0.078$	${0.00}$	$\mu=1.192, \ \sigma=0.086$
	5%	1%	$\mu=2.159, \ \sigma=0.327$	${\bf 0.88}$	$\mu=2.184, \ \sigma=0.294$
		4%	$\mu=1.445, \ \sigma=0.132$	${0.01}$	$\mu=1.470, \ \sigma=0.113$
		10%	$\mu=1.287, \ \sigma=0.092$	${0.00}$	$\mu=1.314, \ \sigma=0.084$
	10%	1%	$\mu=3.096, \ \sigma=0.599$	${\bf 0.49}$	$\mu=3.143, \ \sigma=0.544$
		4%	$\mu=1.684, \ \sigma=0.191$	${\bf 0.53}$	$\mu=1.700, \ \sigma=0.176$
		10%	$\mu=1.410, \ \sigma=0.111$	${\bf 0.10}$	$\mu=1.423, \ \sigma=0.102$
5	1%	1%	$\mu=1.398, \ \sigma=0.115$	${0.00}$	$\mu=1.448, \ \sigma=0.105$
		4%	$\mu=1.084, \ \sigma=0.051$	${0.00}$	$\mu=1.125, \ \sigma=0.061$
		10%	$\mu=1.007, \ \sigma=0.006$	${0.00}$	$\mu=1.012, \ \sigma=0.008$
	3%	1%	$\mu=1.829, \ \sigma=0.207$	${0.00}$	$\mu=1.901, \ \sigma=0.198$
		4%	$\mu=1.340, \ \sigma=0.100$	${0.00}$	$\mu=1.380, \ \sigma=0.106$
		10%	$\mu=1.123, \ \sigma=0.066$	${0.00}$	$\mu=1.167, \ \sigma=0.074$
	5%	1%	$\mu=2.190, \ \sigma=0.293$	${0.00}$	$\mu=2.259, \ \sigma=0.300$
		4%	$\mu=1.480, \ \sigma=0.125$	${0.00}$	$\mu=1.535, \ \sigma=0.123$
		10%	$\mu=1.248, \ \sigma=0.089$	${0.00}$	$\mu=1.293, \ \sigma=0.090$
	10%	1%	$\mu=3.030, \ \sigma=0.571$	${0.00}$	$\mu=3.170, \ \sigma=0.579$
		4%	$\mu=1.737, \ \sigma=0.178$	${0.00}$	$\mu=1.803, \ \sigma=0.171$
		10%	$\mu=1.418, \ \sigma=0.116$	${0.00}$	$\mu=1.465, \ \sigma=0.112$
7	1%	1%	$\mu=1.427, \ \sigma=0.109$	${0.00}$	$\mu=1.457, \ \sigma=0.107$
		4%	$\mu=1.097, \ \sigma=0.057$	${0.00}$	$\mu=1.126, \ \sigma=0.059$
		10%	$\mu=1.008, \ \sigma=0.007$	${0.00}$	$\mu=1.011, \ \sigma=0.009$
	3%	1%	$\mu=1.928, \ \sigma=0.209$	${0.00}$	$\mu=2.000, \ \sigma=0.209$
		4%	$\mu=1.339, \ \sigma=0.093$	${0.00}$	$\mu=1.378, \ \sigma=0.104$
		10%	$\mu=1.140, \ \sigma=0.066$	${0.00}$	$\mu=1.167, \ \sigma=0.067$
	5%	1%	$\mu=2.348, \ \sigma=0.334$	${0.00}$	$\mu=2.447, \ \sigma=0.311$
		4%	$\mu=1.506, \ \sigma=0.122$	${0.00}$	$\mu=1.546, \ \sigma=0.126$
		10%	$\mu=1.254, \ \sigma=0.083$	${0.00}$	$\mu=1.281, \ \sigma=0.085$
	10%	1%	$\mu=3.383, \ \sigma=0.577$	${0.03}$	$\mu=3.481, \ \sigma=0.576$
		4%	$\mu=1.834, \ \sigma=0.195$	${0.00}$	$\mu=1.895, \ \sigma=0.181$
		10%	$\mu=1.436, \ \sigma=0.112$	${0.00}$	$\mu=1.471, \ \sigma=0.115$

Another potentially important parameter is the selected position of the US-robot's base. The individual results for the three positions in use are presented in table 4 together with pairwise p-values for the same statistical test as described before. The null hypothesis is rejected in 31 cases for position 0 and 1, in all 36 cases for position 1 and 2, and in 31 cases for position 0 and 2. The locations of the positions are illustrated in figure 4(a) and position 1 is the one between the legs of the patient. But note that also the feasible ranges of the LIFT angle are different for each position.

Table 4. Results of method D relative to method A (no US-robot) for each robot base position separately. Additionally, the p-values for the null hypothesis that both results are samples from distributions with equal median (Wilcoxon rank sum test) are given for each combination of parameters. Values above 5% are marked bold.

#LIFTs	$\lambda_{\rm LA}$	$\lambda_{\rm US}$	robot 0	p 0,1	robot 1	p 1,2	robot 2	p 2,0
3	1%	1%	$\mu=1.39, \ \sigma=0.10$	${\bf 0.65}$	$\mu=1.39, \ \sigma=0.11$	${0.00}$	$\mu=1.47, \ \sigma=0.06$	${0.000}$
		4%	$\mu=1.10, \ \sigma=0.05$	${0.00}$	$\mu=1.13, \ \sigma=0.08$	${0.04}$	$\mu=1.13, \ \sigma=0.07$	${0.000}$
		10%	$\mu=1.01, \ \sigma=0.01$	${0.00}$	$\mu=1.02, \ \sigma=0.01$	${0.00}$	$\mu=1.01, \ \sigma=0.01$	${0.000}$
	3%	1%	$\mu=1.73, \ \sigma=0.19$	${0.00}$	$\mu=1.76, \ \sigma=0.22$	${0.00}$	$\mu=1.93, \ \sigma=0.11$	${0.000}$
		4%	$\mu=1.35, \ \sigma=0.10$	${\bf 0.81}$	$\mu=1.34, \ \sigma=0.10$	${0.00}$	$\mu=1.42, \ \sigma=0.06$	${0.000}$
		10%	$\mu=1.16, \ \sigma=0.07$	${0.00}$	$\mu=1.18, \ \sigma=0.09$	${0.00}$	$\mu=1.20, \ \sigma=0.08$	${0.000}$
	5%	1%	$\mu=2.07, \ \sigma=0.27$	${0.00}$	$\mu=2.13, \ \sigma=0.34$	${0.00}$	$\mu=2.40, \ \sigma=0.15$	${0.000}$
		4%	$\mu=1.45, \ \sigma=0.11$	${\bf 0.95}$	$\mu=1.43, \ \sigma=0.13$	${0.00}$	$\mu=1.53, \ \sigma=0.08$	${0.000}$
		10%	$\mu=1.29, \ \sigma=0.08$	${\bf 0.44}$	$\mu=1.29, \ \sigma=0.10$	${0.00}$	$\mu=1.35, \ \sigma=0.06$	${0.000}$
	10%	1%	$\mu=2.90, \ \sigma=0.49$	${0.00}$	$\mu=3.07, \ \sigma=0.61$	${0.00}$	$\mu=3.52, \ \sigma=0.33$	${0.000}$
		4%	$\mu=1.65, \ \sigma=0.17$	${0.00}$	$\mu=1.66, \ \sigma=0.20$	${0.00}$	$\mu=1.83, \ \sigma=0.10$	${0.000}$
		10%	$\mu=1.41, \ \sigma=0.10$	${\bf 0.26}$	$\mu=1.39, \ \sigma=0.12$	${0.00}$	$\mu=1.49, \ \sigma=0.07$	${0.000}$
5	1%	1%	$\mu=1.46, \ \sigma=0.10$	${0.00}$	$\mu=1.38, \ \sigma=0.10$	${0.00}$	$\mu=1.49, \ \sigma=0.12$	${0.000}$
		4%	$\mu=1.12, \ \sigma=0.06$	${0.00}$	$\mu=1.08, \ \sigma=0.05$	${0.00}$	$\mu=1.14, \ \sigma=0.06$	${0.000}$
		10%	$\mu=1.01, \ \sigma=0.01$	${0.00}$	$\mu=1.01, \ \sigma=0.01$	${0.00}$	$\mu=1.01, \ \sigma=0.01$	${\bf 0.077}$
	3%	1%	$\mu=1.92, \ \sigma=0.17$	${0.00}$	$\mu=1.80, \ \sigma=0.19$	${0.00}$	$\mu=1.97, \ \sigma=0.23$	${0.000}$
		4%	$\mu=1.38, \ \sigma=0.09$	${0.00}$	$\mu=1.32, \ \sigma=0.09$	${0.00}$	$\mu=1.43, \ \sigma=0.10$	${0.000}$
		10%	$\mu=1.16, \ \sigma=0.06$	${0.00}$	$\mu=1.12, \ \sigma=0.07$	${0.00}$	$\mu=1.19, \ \sigma=0.07$	${0.000}$
	5%	1%	$\mu=2.29, \ \sigma=0.24$	${0.00}$	$\mu=2.11, \ \sigma=0.28$	${0.00}$	$\mu=2.41, \ \sigma=0.34$	${0.000}$
		4%	$\mu=1.55, \ \sigma=0.11$	${0.00}$	$\mu=1.47, \ \sigma=0.12$	${0.00}$	$\mu=1.57, \ \sigma=0.14$	${0.000}$
		10%	$\mu=1.29, \ \sigma=0.07$	${0.00}$	$\mu=1.23, \ \sigma=0.08$	${0.00}$	$\mu=1.33, \ \sigma=0.09$	${0.000}$
	10%	1%	$\mu=3.22, \ \sigma=0.47$	${0.00}$	$\mu=2.86, \ \sigma=0.55$	${0.00}$	$\mu=3.52, \ \sigma=0.63$	${0.000}$
		4%	$\mu=1.83, \ \sigma=0.15$	${0.00}$	$\mu=1.72, \ \sigma=0.17$	${0.00}$	$\mu=1.86, \ \sigma=0.21$	${0.001}$
		10%	$\mu=1.47, \ \sigma=0.10$	${0.00}$	$\mu=1.40, \ \sigma=0.11$	${0.00}$	$\mu=1.51, \ \sigma=0.13$	${0.000}$
7	1%	1%	$\mu=1.42, \ \sigma=0.09$	${0.00}$	$\mu=1.50, \ \sigma=0.08$	${0.00}$	$\mu=1.43, \ \sigma=0.13$	${\bf 0.175}$
		4%	$\mu=1.12, \ \sigma=0.06$	${0.00}$	$\mu=1.14, \ \sigma=0.06$	${0.00}$	$\mu=1.10, \ \sigma=0.06$	${0.000}$
		10%	$\mu=1.01, \ \sigma=0.01$	${0.00}$	$\mu=1.01, \ \sigma=0.01$	${0.00}$	$\mu=1.01, \ \sigma=0.01$	${0.000}$
	3%	1%	$\mu=1.98, \ \sigma=0.16$	${0.00}$	$\mu=2.07, \ \sigma=0.14$	${0.00}$	$\mu=1.93, \ \sigma=0.27$	${0.000}$
		4%	$\mu=1.34, \ \sigma=0.07$	${0.00}$	$\mu=1.41, \ \sigma=0.08$	${0.00}$	$\mu=1.36, \ \sigma=0.11$	${0.001}$
		10%	$\mu=1.16, \ \sigma=0.06$	${0.01}$	$\mu=1.17, \ \sigma=0.06$	${0.00}$	$\mu=1.15, \ \sigma=0.07$	${0.000}$
	5%	1%	$\mu=2.40, \ \sigma=0.24$	${0.00}$	$\mu=2.55, \ \sigma=0.22$	${0.00}$	$\mu=2.37, \ \sigma=0.43$	${0.029}$
		4%	$\mu=1.51, \ \sigma=0.10$	${0.00}$	$\mu=1.58, \ \sigma=0.10$	${0.00}$	$\mu=1.53, \ \sigma=0.15$	${0.047}$
		10%	$\mu=1.26, \ \sigma=0.07$	${0.00}$	$\mu=1.32, \ \sigma=0.08$	${0.00}$	$\mu=1.26, \ \sigma=0.09$	${\bf 0.932}$
	10%	1%	$\mu=3.40, \ \sigma=0.45$	${0.00}$	$\mu=3.71, \ \sigma=0.42$	${0.00}$	$\mu=3.39, \ \sigma=0.75$	${\bf 0.204}$
		4%	$\mu=1.88, \ \sigma=0.14$	${0.00}$	$\mu=1.95, \ \sigma=0.13$	${0.00}$	$\mu=1.81, \ \sigma=0.25$	${0.000}$
		10%	$\mu=1.44, \ \sigma=0.09$	${0.00}$	$\mu=1.51, \ \sigma=0.09$	${0.00}$	$\mu=1.45, \ \sigma=0.14$	${\bf 0.129}$

Method D uses the heuristic GTSP-solver GLKH. As for all iterative solvers, we can expect better solutions if we apply more iterations and restarts with different initial guesses. Figure 8 shows the relative improvements compared to our baseline settings. Note that these results are only based on a subset of our parameter combinations. Improvements up to 35% are observable in the diagrams. However, the GLKH runtime is basically linear in number of iterations and restarts (if not calculated in parallel). Tripling the computation time by applying 3000 iterations led to an improvement of at most 6% in 90% of the cases. The same holds for tripling the time by applying 2 restarts.

Zoom In Zoom Out Reset image size

Figure 9 compares the number of necessary configuration changes by the US-robot in the solutions of methods C and D. More configuration changes correspond to more robot motion close to the patient. The faster the robot with the linear accelerator moves, the clearer becomes the difference between the method. Method C is only able to optimize the sequence of US-robot configurations but uses a fixed sequence of beams. Method D optimizes both sequences and results in fewer changes of the configuration.

Zoom In Zoom Out Reset image size

We present methods to coordinate and optimize the motion of two robots in a scenario where one robot carries a linear accelerator to deliver beams of radiation and the second robot positions an ultrasound transducer for continuous image guidance. The LINAC-robot motion will always contribute to the overall treatment time if beams are delivered from different directions. Different approaches to optimize the path planning problem have been described (Dieterich and Gibbs 2011, Kearney et al 2018). We confirm that the TSP related to this problem (figure 5) can be very efficiently solved by state-of-the-art solvers. While the robot's travel time would depend on the desired velocity limits, these are typically far below the robot's specification and do not change, i.e. the same TSP solution would be valid even if the robot was moving twice as fast.

The CyberKnife M6 limits the end-effector velocity of its robot in Cartesian space to 60 mm s⁻¹ or 80 mm s⁻¹, depending on the distance to certain obstacles. Instead, we limit the joint velocities in our simulations and therefore obtain differently shaped distributions for the velocities in Cartesian space (figure 6). However, the US-robot does not move its end-effector and therefore we chose to consistently control the joint velocities for both robots. The observed velocities in our simulations are roughly centered around 20, 60, 100, and 200 mm s⁻¹ which corresponds to 1%, 3%, 5%, and 10% of the maximum velocity in Euclidean space, which is 2 m s⁻¹.

Considering an actively moving US-robot in this setup allows for more feasible beam directions. Table 1 showed that inverse planning led to disproportionately many beams which are blocked in at least one of the considered US-robot configurations being used for treatment. However, the main purpose of our work is to study the impact that integrating an US-robot would have on the treatment time. If the US-robot is able to change its configuration in a shorter time than the LINAC-robot needs to move to deliver the next beam, then the total time for robot motion is the same as without the US-robot. However, this depends not only on the distance which each robot has to travel but also on their speeds. The results in figure 7 and table 2 illustrate that a combination of a 1% LINAC-robot speed and 10% US-robot speed would indeed result in a negligible increase in travel time of approximatley 1%. In practice, the speed of the robotic CyberKnife is limited to at most 4% of the maximum robot speed (figures 6(b) and (c)). Considering $\lambda_{\rm LA}=3$% and a similarly conservative motion setting of 4% of maximum speed for the US-robot, we get an increase in travel time in the order of 30%. This can be explained by the fact that the US-robot sometimes has to move its elbow substantially to avoid beam blocking. For all methods and numbers of LIFT angles in figure 7, the additional robot motion corresponds to a median increase of the overall time by less than 5%. However, it is important to note that for x-ray-based image-guided treatments the linac also needs to move out of the x-ray beam path. Thus, our solution from method A had to be corrected by the time needed for these movements as well, if the current clinical setup with x-ray imaging is considered. Additionally, the x-ray imaging is limited to the time when the linac does not block the beam path. Instead, our setup allows for continuous imaging and tracking. This is especially an advantage if spontaneous motions of the organs occur, because these would not be detected until the next x-ray images are acquired and thereby corrupt the actually achieved dose distribution.

Comparing the results for method B to the results for methods C and D (figure 7), it is clear that simply using the same beam sequence as without the US-robot is not optimal. Interestingly, the GTSP instances resulting from an inclusion of the US-robot motion are no longer solvable exactly in a reasonable amount of time. Both methods C and D employ heuristics and generally there is no guarantee to find a global optimum. While method D employs a more sophisticated heuristic included in the GLKH solver, it does not always result in the shortest motion. Particularly when the speed of the LINAC-robot is low, our proposed heuristic C compares favorably. This is remarkable, as method C is both deterministic and substantially faster than method D. The computation time for the baseline GLKH version with 1000 itarations and no restart is already in the order of 1 h compared to few seconds for method C (inluding solving the TSP). However, in a clinical setting the allowed velocity of the US-robot will be limited. For the cases with the US-robot being slower than the LINAC-robot, method D is able to further improve the solution in most situations, except when a large number of LIFT angles is considered.

The quality of the solutions by method D can be further improved in general by increasing the number of iterations and restarts (figure 8), but this comes with very high costs in terms of additional computational time. However, method D often gives solutions with fewer configuration changes for the US-robot, i.e. motion close to the patient, compared to method C (figure 9). Note that if the linac is moved slowly, the US-robot will still change its configuration more than 50 times.

Our analysis also indicates that the US-robot's placement has an impact on the total travel time. We mostly observed significant deviations of the medians for the three robot setups under consideration (table 4). Similarly, the results for patients with a small prostate and patients with a large prostate showed significant deviations (table 3). For the large prostates a higher temporal overhead was observed.

While approaches to minimize the impact of intra-fractional ultrasound imaging have been proposed, they typically come with some limitations, too. The Clarity system (Elekta) is used for setup and motion detection during radiation therapy of the prostate (Abramowitz et al 2012, Western et al 2015, Grimwood et al 2018). The system uses a transperineal approach and is thus limited to a relatively small volume. Another approach is to construct the ultrasound transducer such that beams can pass without substantial attenuation (Schlosser and Hristov 2016). However, the design does not include the robotic arm and the CyberKnife's relative beam motion due to active motion compensation may complicate dose calculations accounting for any residual attenuation. In contrast, we have demonstrated that using a commercial 7-DoF robotic arm and coordinating the motion of both robots allows for practical treatment times.

So far, configuration changes by the US-robot are only considered between delivered beams. Instead, one could also consider to move the US-robot during beam delivery to further decrease the overhead of its motion. However, separate motion and beam delivery might be more realistic in terms of safety, because it allows to check beam-wise for a conflict between the beam and the US-robot configuration. Otherwise, a fine-grained evaluation of the whole trajectory might be required. Further improvements may also be possible by allowing arbitrary LIFT angles, which would potentially reduce the amount of motion until subsequent beams can be delivered.

Robotic ultrasound guidance in radiation therapy is a promising approach to realize non-invasive and non-ionizing real-time tracking of the motion and deformation of the target and surrounding structures. Using a kinematically redundant robot carrying the ultrasound probe, the effect of beam blocking can be largely mitigated. We study methods to also account for the potential prolongation of the treatment delivery due to configuration changes of the ultrasound robot. While we illustrate that the optimization problem to coordinate robotic beam delivery and robotic ultrasound is hard, we also show that practical solutions can be obtained with the proposed methods. We also demonstrate that the effect of robot coordination depends on the speed settings, the ultrasound robot's base position, and the target size. Considering that continuous ultrasound motion monitoring would be particularly helpful to detect spontaneous motion in the abdomen, the relatively small overhead in treatment time of less than 4% would be acceptable.

This work was partially funded by Deutsche Forschungsgemeinschaft (Grant SCHL 1844/3-1).