Optimizing robot motion for robotic ultrasound-guided radiation therapy

An important aspect of robotic radiation therapy is active compensation of target motion. Recently, ultrasound has been proposed to obtain real-time volumetric images of abdominal organ motion. One approach to realize flexible probe placement throughout the treatment fraction is based on a robotic arm holding the ultrasound probe. However, the probe and the robot holding it may obstruct some of the beams with a potentially adverse effect on the plan quality. This can be mitigated by using a kinematically redundant robot, which allows maintaining a steady pose of the ultrasound probe while moving its elbow in order to minimize beam blocking. Ultimately, the motion of both the beam source carrying and the ultrasound probe holding robot contributes to the overall treatment time, i.e. beam delivery and robot motion. We propose an approach to optimize the motion and coordination of both robots based on a generalized traveling salesman problem. Furthermore, we study an application of the model to a prostate treatment scenario. Because the underlying optimization problem is hard, we compare results from a state-of-the-art heuristic solver and an approximation scheme with low computational effort. Our results show that integration of the robot holding the ultrasound probe is feasible with acceptable overhead in overall treatment time. For clinically realistic velocities of the robots, the overhead is less than 4% which is a small cost for the added benefit of continuous, volumetric, and non-ionizing tracking of organ motion over periodic x-ray-based tracking.

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, Camps et al 2018. Ipsen et al (2016) demonstrated 4D tracking for MLC-based motion compensation by acquiring US volumes with 0.3 × 0.8 × 1.6 mm 3 resolution and sizes from 17 × 14 × 45 mm 3 at 250 Hz volume rate to 68 × 105 × 170 mm 3 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 . 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 . 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(Bell et al , Ş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(Gerlach et al , 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.

Material and methods
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 (LINACrobot) 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.

Treatment plan optimization
In robotic radiation therapy, beams are delivered from a finite set of points roughly distributed over a hemisphere around the patient. The points are often called nodes. At each node, beams with different orientation and shape can be delivered. Different collimators shaping the aperture are in use with the CyberKnife, including cylinder collimators, an IRIS collimator (Echner et al 2009), and a multi-leaf collimator (MLC) (Fürweger et al 2016). However, for our purpose the shape of the beams is not important, as the blocking of complete or partial beams mainly depends on the overall beam direction and the configuration of the US-robot. For simplicity we consider circular beams of different diameter and orientation, i.e. starting at one node and passing through the PTV. Furthermore, given the infinite number of possible beams we adopt the typical approach to generate a set of candidate beams heuristically. Given the set of candidate beams the beam activation times or beam weights can be obtained by inverse planning (Schlaefer and Schweikard 2008).
In our scenario, some of the beams might be blocked by the US-robot, i.e. they would pass through either the robot or the ultrasound probe. For each node and each US-robot configuration it is possible to compute the projection of robot and transducer into a plane orthogonal to a line connecting the node and the centroid of the PTV. Using a distance transform on the projection and accounting for uncertainties we can establish whether a beam overlaps with the projection and therefore is blocked (Gerlach et al 2017a). Figure 1(b) summarizes this approach. Considering the set of suitable configurations of the US-robot, we discard and re-sample beams blocked in all of these configurations. Hence, we eventually arrive at the desired number of candidate beams which are feasible in at least one configuration of the US-robot.
Our actual inverse planning uses a step-wise multi-criteria approach based on linear programming (Schlaefer and Schweikard 2008). Note that there are two key advantages of this approach in the context of a low level analysis of treatment planning. First, the fixed bounds on all upper dose bounds prevent any side effects, e.g. small violations of the dose constraints for OAR. Second, having access to the actual objective value from the optimization underlying the inverse planning we can compare different algorithms without confounding by effects due to discretization and interpolation.

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 −1 for the first five joints and is 175 deg s −1 for the last. For the LBR iiwa, the first three joints are limited to 98-100 deg s −1 , the fourth to 130 deg s −1 , the fifth to 140 deg s −1 , and the last two joints to 180 deg s −1 . 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.

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 (1) to each edge (v i , v j ) ∈ 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 In our scenario, each treatment beam b i is related to a vertex such that 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 where θ i g 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 θ = θ + n · 2π, n ∈ Z, and are wrapped to the half-open interval [−π, π). The denominator is the velocity of the joint, which we assume to be a fraction λ LA ∈ (0, 1] of its maximum velocity s LA g . 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.

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 We create a new model where for each beam b j all pairs (b j , c k ), c k ∈ 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) tō 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.
Thus, we have N groups with at least one and at most M vertices each which form the overall set of vertices with edges 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 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 · M N possible solutions.

TSP GTSP
Structures of the graph models. In the TSP model (left), each vertex corresponds to one beam b i which has to be delivered by the linear accelerator. Instead, the GTSP (right) has one vertex for each pair (b i , c k ) of beam and configuration of the US-robot, in which this beam is not blocked. Vertices representing the same beam are not connected, but form a group (dashed ellipses).
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-ofthe-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 compariso n we also solve the TSP optimizing just the LINAC-robot motion without adding the US-robot.

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.

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.

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 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 A . Likewise, all vertices containing the last beam in S 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 .

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.

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(Gerlach et al , 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 3 to 115.4 cm 3 with mean 84.5 cm 3 and standard deviation 13.8 cm 3 . 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).

Robot configurations
To investigate the effect of robot speed on treatment time, we assume 1%, 3%, 5%, and 10% for λ LA and 1%, 4%, 7%, and 10% for λ US . In the actual CyberKnife system, the linear accelerator is moved with approximately 3%-4% of the maximum speed in Cartesian space.

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 −1 . 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.

Results
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 Figure 3. Illustration of method B for four beams and three configurations of the US-robot. Given the optimal beam sequence for the linear accelerator, i.e. b i is the ith beam in the sequence, a directed graph is constructed. The shortest path from source to sink gives the best sequence of US-robot configurations under the fixed beam order. The costs of the edges are equal to those of the GTSP model. Infeasible combinations of beam and configurations are shown grayed out. Therefore, at least one configuration change (from c 1 to c 2 ) is required to deliver all beams in this example.  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. 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 λ LA = 3% and 3.8% for λ LA = 10%. For the slow linac with λ 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.
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 USrobot 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 LINACrobot 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.
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 3 ) and the three patients with the largest PTVs (93.2, 109.7, and 115.4 cm 3 ). 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.
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.
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 4 11.7 ± 5.8 12.8 ± 3.5 3.9 ± 1.6 6.9 ± 1.9 5 77.3 ± 9.3 67.7 ± 4.2 6.9 ± 3.9 8.7 ± 2.9 6 7.8 ± 3.6 7.6 ± 2.2 7 76.6 ± 8.9 66.7 ± 3.6 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. 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.

Discussion
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 andGibbs 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 −1 or 80 mm s −1 , 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 −1 which corresponds to 1%, 3%, 5%, and 10% of the maximum velocity in Euclidean space, which is 2 m s −1 .
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 Table 2. Mean µ and standard deviation σ 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. 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 λ 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 . 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.

Conclusion
Robotic ultrasound guidance in radiation therapy is a promising approach to realize non-invasive and nonionizing 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.