Leveraging respiratory organ motion for non-invasive tumor treatment devices: a feasibility study

In noninvasive abdominal tumor treatment, research has focused on minimizing organ motion either by gating, breath holding or tracking of the target. The paradigm shift proposed in this study takes advantage of the respiratory organ motion to passively scan the tumor. In the proposed self-scanning method, the focal point of the HIFU device is held fixed for a given time, while it passively scans the tumor due to breathing motion. The aim of this paper is to present a treatment planning method for such a system and show by simulation its feasibility. The presented planning method minimizes treatment time and ensures complete tumor ablation under free-breathing. We simulated our method on realistic motion patterns from a patient specific statistical respiratory model. With our method, we achieved a shorter treatment time than with the gold-standard motion-compensation approach. The main advantage of the proposed method is that electrically steering of the focal spot is no longer needed. As a consequence, it is much easier to find an optimal solution for both avoiding near field heating and covering the whole tumor. However, the reduced complexity on the beam forming comes at the price of an increased complexity on the planning side as well as a reduced efficiency in the energy distribution. Although we simulate the approach on HIFU, the idea of self-scanning passes over to other tumor treatment modalities such as proton therapy or classical radiation therapy.


Introduction
High intensity focused ultrasound (HIFU) is a well-known non-invasive thermal ablation modality for tumor treatment which is widely accepted for about a decade (Schwenke et al 2015, ter Haar et al 1989, Cline et al 1992. For image guidance during HIFU sonication, ultrasound (Kennedy 2005) as well as magnetic resonance (MR) can be used (Hynynen 2010. MR not only provides images of the tumor and its motion, but is also used for temperature mapping (Hynynen 2010). With MR thermometry, the focal spot position and the heating procedure are monitored. Further, the thermal dose is calculated to determine the tissue damage (Hynynen 2010, Wijlemans et al 2012. MR-guided HIFU has been successfully applied for tumor ablation in immobile organs, such as uterus, prostate, breast and brain (Hynynen 2010, Ellis et al 2013. However, HIFU treatment of abdominal organs, such as kidney and liver, remains challenging due to respiratory motion and organ drift (Von Siebenthal et al 2007).
Today, there are two established principles to address the problems of liver motion in HIFU sonication: the respiratory-gated and the continuous motion-compensation sonication strategies (De Senneville et al 2012). In the former ones, an almost stationary part of an adult breathing cycle is exploited. Within a temporal window called gating-window of 1-2 s, where the liver remains approximately still, the acoustic energy is periodically deposited . The disadvantages of the respiratory-gated method is the prolonged treatment time and the negligence of the remaining motion during the gating-window.
Continuous motion-compensation sonication, in contrast, is based on the idea of continuously readjusting the focal point position to the current tumor position in order to prevent both undesired tissue damage and energy spread (Ries et al 2010, Arnold et al 2011, Holbrook et al 2014. The advantage of this method compared to the respiratory-gated approach is the near 100% duty cycle, which results in a shorter treatment time. However, in order to steer the focal spot of the HIFU device, a phased-array-transducer with hundreds or even thousands of elements is required, each having its own amplifier and delay time. This complicates the electrical design of the device and is expensive. Moreover, the lateral steering causes an intensity decay at the focal spot . A further challenge of the motion-compensation sonication is the prediction of the liver motion to actively steer the HIFU beam. The idea of direct motion tracking is to detect the tumor in real-time and adjust the focal spot accordingly. Due to processing latency, the tumor position has to be anticipated to provide direct motion tracking (Ries et al 2010, Preiswerk et al 2014. Preiswerk and Cattin (2015) used a non-linear Gaussian process regression method to predict organ motion based on a model-topology independent external respiratory signal.
As solid tumors are much larger than the focal spot of the HIFU device, ablating the whole tumor requires several treatment spots. The development of a sonication plan that reduces the treatment time and generates a uniform lesion has been addressed before in immobile tissue (Salomir et al 2000, Mougenot et al 2004, Hui et al 2009, Zhou 2013. The basic idea is to sonicate points in a regular grid, such as a spiral pattern or raster scanning. The ultrasound waves emitted by the HIFU device cause a pressure field which induces a temperature rise at the focal spot. The temperature can be described by Pennes bio-heat equation (Cline et al 1992, Wissler 1998. The most accepted model to determine how tissue is affected by the induced temperature rise is the thermal dose (Sapareto and Dewey 1984). The thermal dose equation estimates the equivalent minutes at a temperature of 43 °C. Another study (Rosenberg et al 2013) shows that a temperature of 54 °C can be sufficient to induce necrosis.
In this paper, we present a planning method for the proposed self-scanning approach and show its feasibility by simulation. The novel method takes advantage of the perpetual respiratory motion to passively scan the tumor. In other words, we are placing the static focal point of the HIFU into the tumor. The motion caused by breathing moves the tumor through this focal point. Based on the motion model of (Preiswerk and Cattin 2015), we anticipate at which time point tumor tissue is located under the focal spot and thus modulate the HIFU intensity based on this information. Once the tumor has been ablated along the self-scanned trajectory, the focal spot is relocated to a different but static position within the body. The main advantage of this principle is that no electronic steering of the focal spot is needed. Actively compensating for the motion can collide with both avoiding near field heating and covering the sonication volume. As a consequence, it is much easier to find an optimal solution for beam forming with the proposed self-scanning method. Moreover, the amplitude of the respiratory motion is comparable or larger than the steering range of a standard phased-array-transducer. Thus, the spatial extension of the sonicated pattern in the tumors referential is equal or larger with self-scanning techniques than with electronic steering of the beam. However, the reduced complexity on the beam forming comes at the cost of an increased complexity during the planning stage and a reduced efficiency in the energy distribution. With the proposed method, we combine the advantages of the respiratory-gated and motion-compensation method: a fixed focus can be used and a high duty cycle is achieved. Moreover, since with the proposed selfscanning approach no lateral steering of the focal spot is required, position-dependent decay of the focal spot intensity during lateral steering is avoided.
In the following, we present a method to determine optimal sonication plans for such a system ensuring complete tumor ablation. In section 2, we introduce a mathematical formulation of our model and show how to solve it. The model minimizes treatment time and handles complex non-rigid motion and non-time dependent diffusion. In section 3, we simulated our method on a patient specific statistical respiratory motion model. We compared the results to a continuous motion-compensation approach. The different approaches were evaluated in terms of efficiency and treatment time. The efficiency coefficient is a measure to determine how emitted energy can be exploited. With simulation, we showed that our method provides feasible treatment plans for the proposed self-scanning approach, which yields a faster treatment time than the continuous motion-compensation approach.

Method
In the following section, we show how a treatment that takes advantage of the respiratory motion can be planned. The task is to find appropriate tumor points which are sonicated by the HIFU device. The focal spot will stay fixed during a given time to achieve a precalculated temperature rise. During this phase, different liver tissue will pass through the focal spot. The points and the corresponding temperature rises have to be chosen such that the whole tumor is ablated. To avoid overtreatment, the energy has to be distributed mainly on the tumor and as few as possible energy should be spread on healthy tissue. Under these conditions, the treatment time is minimized. The treatment time consists of beam-and changing time. The beam time is the overall time where the HIFU device is focused on one point. The HIFU system used with the self-scanning method is able to electrically steer the focal spot in depth along the acoustic axis. For the other two directions, mechanical displacement is used. Hence, the focal spot can be changed fast along the acoustic axis and slower in the other two directions. The time needed to change in the two slower directions is called changing time.
The sonication plans are calculated before treatment. As described in the next section, we assume a non-rigid, but repetitive motion and breathing pattern. We therefore calculate a lookup table that tells us at every time point which position inside the target has to be sonicated.

Assumptions
In this first proof of principle study, we assume that there is no heat diffusion over time. This means that a liver point will only change its heat value if it is inside the focal spot of the HIFU device for a certain time during treatment. In other words we assume linear dose accumulation. We will also not consider heat changes due to blood perfusion and metabolic processes. Note that volumetric ablation is less sensitive to diffusion and perfusional sink effects as the thermal build up is rather smooth compared to sharp spatial profiles generated by a fixed focus in static tissue (Köhler et al 2009). Further, we assume that the liver motion stays periodic over time. In other words, we assumed that no organ drift occurs and that the breathing pattern of the patient stays constant during treatment. This means that the breathing pattern repeats itself after each breathing cycle. This corresponds to the assumption of a regular breathing pattern of the patient. Note however, that our model can handle complex non-rigid breathing motion.
In order to determine which tissue has been ablated, the most accepted model is the thermal dose equation (Sapareto and Dewey 1984). However, Rosenberg et al (2013) showed that tissue is ablated when heated above a temperature of 54 °C, which we will use in this study.
The above assumptions are simplifications of the real problem in order to show the feasibility of the self-scanning approach. In our ongoing research, we will include diffusion over time as well as perfusion. Moreover, we plan to switch to the thermal dose model.

Mathematical model
For the mathematical formulation, we consider three sets of points in different coordinate systems. The first set consists of points that might be heated up during the treatment. These points may be inside or outside the tumor and are in the liver coordinate system. We will call this set liver points P l . The second set includes the target points, which form a subset of the liver points, ⊆ P P t l . The target points include tumor points as well as some surrounding tissue points. To be sure that all of the tumor cells will be ablated, some surrounding tissue of the tumor has to be treated as well (safety margin). The goal of the HIFU treatment is to heat the target points and ablate them. The last set consists of points in the world coordinate system which can be focused by the HIFU device. We will refer to these as world points P w . The main difference between the two sets P l and P w is that the liver points will move in the world coordinate system during breathing, whereas the points from the second set are still.
On the basis of speed and direction of the liver movement, we will set up a cover matrix M similar as in Trofimov et al (2005). This matrix indicates which liver points will be heated during free-breathing when sonicating a certain point in the world coordinate system. Each matrix column j represents the heat profile of all the liver points after a breathing cycle when a world point p j is sonicated. These heat profiles can be thought of as cylinders. The shape of a cylinder is determined by the liver motion at the specific location (length and direction) and the focal spot size (diameter). A cylinder represents the tissue which is heated during one breathing cycle when a certain world point is sonicated. In figure 1, a schematic 2D setup is N Möri et al Phys. Med. Biol. 61 (2016) 4247 shown. In this case, a column of the covermatrix is drawn as a rectangle that represents the treated liver points. For simplicity, the motion of the liver shown in figure 1 is linear. In fact, the motion is nonlinear and the rectangles would have a more complex form. Note however, that our model handles complex non-rigid breathing motion. The task in finding a feasible treatment plan can be formulated as choosing appropriate cylinders which cover the whole tumor. Mathematically, this corresponds in finding a vector ∈ R x m , which has positive values at the positions where the corresponding cylinders are part of the treatment plan. Note, however, that all the cylinders represent a heat profile, and thus there have to be some overlaps of the cylinders to make sure that all of the target points will be ablated. Moreover, the cylinders are banana-shaped due to the ellipsoidal breathing motion.
Suppose we have n liver points where I(i, j) is the total temperature rise at liver point l i during one breathing cycle when focusing world point p j . Note that for ∈ + R x m , the values of Mx represent how much the temperature rises at each liver point, when p j is sonicated by the HIFU device with intensity x j , for j = 1, ..., m. We assume that tissue is ablated when heated above a temperature of 54 °C (Rosenberg et al 2013). The goal is therefore to provide a minimal temperature rise of α = 18 °C to each ∈ l P i t and simultaneously minimize the temperature rise for points outside the target. The target vector ∈ R T n is defined as: To ablate the target, we want to find ∈ + R x m that satisfies ⩾ Mx T. (2) We call the solution x to the above equation sonication or treatment plan. We assume that the temperature rise induced by the HIFU device can be approximated with a Gaussian blurred ellipsoid with deviations σ z in parallel direction, and σ r perpendicular to the acoustic axis. The formula of I(i, j) used in the definition of the covermatrix M in equation (1) is as follows: Hence, the maximal temperature rise that can be delivered by focusing a certain point p j with intensity x j = 1 yields the lethal temperature rise of 18 °C.

Optimization
The beam time is related to the number of points that have to be sonicated to ablate the entire target volume, i.e. the non-zero entries of x. Hence, we want to minimize equation (2) over the number of non-zero entries of x, which is the 0 -norm. However, the minimization problem with the 0 -norm is an NP-hard problem (Natarajan 1995) and we thus minimize equation (2) over the 1 -norm of x, which is a good approximation.
Note that when minimizing ∥ ∥ x 1 , overtreatment outside as well as inside the target is minimized. Here ⩾ x 0 means that ⩾ x 0 k ∀ k, and the 1 -norm is defined as The above equation can be stated as a linear programming problem in the following way: The advantage of this approach is that the solution will fulfill the inequality constraints exactly. The ablation of the whole target volume is thus guaranteed. However, the disadvantage of the linear programming algorithm is its computational complexity, which is ( ) O m 3.5 (Robere 2012), where m is the number of variables (in our case, ≈ m 63 000). To reduce the complexity, the minimization problem in equation (4) can also be stated in an unconstrained version, which yields a convex problem:  with λ > 0. The above formulation is a Lagrange Multiplier formulation of the compressed sensing problem, which can be minimized with the weighted spectral projected gradient algorithm (wspgl1) proposed by Mansour (2012). The wspgl1 algorithm has a complexity of ( ) O m (Birgin et al 2014). However, the numerical solution of the 1 -formulation in equation (6) will not fulfill the equality constraint exactly. This leads to a treatment plan that does not completely ablate the target volume. To get a feasible sonication plan which ablates the whole target volume, the solution x to equation (6) needs to be scaled, which leads to a non-optimal sonication plan.
For the scaling, we define the index sets Thus, the minimal temperature rise inside the target volume is α and the scaled solution x thus ablates the whole target volume. Note that it is possible that Hence, the scaling not only ensures complete target volume ablation, but it also enhances overtreatment.
In practice, the treatment plans have to be adjusted during treatment due to unpredictable organ drift and changing breathing cycle of the patient during treatment. With these changes in the respiratory motion, the calculated treatment plans have to be adjusted to guarantee complete target volume ablation. Therefore, the shorter computation time of the compressed sensing approach in equation (6) would be preferred.

Reducing overtreatment
There are two different kinds of overtreatment that occur: treating healthy tissue outside of the tumor and heating target volume tissue above the lethal temperature rise α. Both types of overtreatment should be avoided, not only to minimize treatment time, but also to avoid unwanted tissue damage. When the induced temperature rise inside the target is too high, the tissue can start to boil, which is unwanted. To avoid treating tissue outside the target, the HIFU device is turned off when healthy tissue is exposed. To account for this, the covermatrix is slightly changed to with the same notation as in equation (3). The set S j denotes the time where most of the points inside the focal spot are target tissue. In the specific case, if the target volume is a sphere, S j can be defined as , where m(t) denotes the target center coordinates at time t, and R is the radius around the target center in which the HIFU device is allowed to treat tissue.
To reduce overtreatment inside the target, first a treatment plan x has to be calculated by equations (5) or (6). Then, the breathing cycle is split into p time steps (2016) 4247 A new covermatrix is calculated which indicates the temperature rise induced by the HIFU device during this shorter time periods. When thinking of cylinders, the time splitting means that we divide the height of each cylinder into p smaller ones. With the splitting, another minimization problem has to be solved to determine the temperature rises at the different time steps. Let B = {i: x i > 0} denote the indices of the points that are sonicated in the treatment The minimization problem to determine the temperature rise at the different time steps can be formulated as x 0 s 1 s Note that if we would not calculate a solution x first, the minimization problem in equation (9) had ⋅ p n variables. With the knowledge of a treatment plan, this number can be reduced to ⋅| | p B.

Time measurement
To estimate the beam time, we assume that we need 5 breathing cycles to induce the lethal temperature rise α. In other words, a temperature rise of 3.6 °C can be achieved in one breathing cycle. Hence, for the treatment plan x, the beam time can be calculated as The treatment time additionally depends on the time needed to change the HIFU device position. The time needed to change with the mechanical displacement in the two slower directions is denoted by c, the electric steering along the acoustic axis can be neglected. If the fast direction is along the x-axis, the time function for a change between two points ( ) = P x y z , , The total changing time is ⋅ N c, if we have to perform N mechanical displacements. Therefore, we check how many different x-coordinates occur among the focus points of the sonication plan in order to calculate the changing time.

Efficiency coefficient
In order to compare our methods, we calculate the conversion efficiency η. The efficiency coefficient measures how much energy is wasted by overtreatment. Overtreatment occurs when tissue outside the target is treated or when target tissue is heated above α. The formula for η is as follows: where E in is the energy that is put inside the patient by the HIFU device, and E out is the energy that can be used to ablate the target. Note that E out is the minimal energy needed to ablate the target, which is ∥ ∥ T 1 . The efficiency of the optimal solution is η = 1 opt .

Continuous motion-compensation method
To evaluate the different methods, we compare the performance of the results to a continuous motion-compensation approach. In the motion-compensation approach, a covermatrix which respects the focal spot steering is computed. As a consequence, the focal spot range is smaller and the temperature rise is higher as compared to the self-scanning approach. Moreover, we modeled the intensity decay at the focal spot during steering with a Gaussian decay in the perpendicular direction to the acoustic axis of the HIFU device . Along the acoustic axis, the intensity decay can be neglected. We assumed that the highest intensity can be achieved at the target center, and that the deviation σ I of the Gaussian is 12 mm. In the following, the continuous motion-compensation approach is referred to as ContMotion.

Results
As mentioned in section 2.3, we consider two different optimization problems, a compressed sensing and a linear programming approach. With a different choice of the point sets, the compressed sensing method is split into two different methods, OT and WL, which will be explained in the next section 3.1. In sections 3.3-3.9, the results of our simulations are explained. We evaluated the treatment time and the energy distribution efficiency. Further, we looked at sagittal and axial slices through the liver and which temperature rise is induced into the target points. We tested our methods on realistic liver breathing motion data using a patient specific statistical respiratory motion model (Preiswerk and Cattin 2015). We set the target to be a sphere inside the liver such that all the energy delivered by the HIFU device will be emitted into the liver. The target with a radius of 15 mm is located at the liver dome, where the motion is non-rigid.
If nothing else is mentioned, the results in the following are obtained by minimizing with covermatrix M as defined in equation (1) with I(i, j) defined in equation (3).

Choice of point sets
We made the following choice of the three point sets: P w consists of target points, since we do not want to sonicate healthy tissue. For the liver points we had two different approaches. The first one was to consider liver points which fulfill ⊊ P P t l . This leads to a minimization not only inside, but also outside the target. Therefore, as few as possible energy is wasted on healthy tissue. We refer to this method as whole liver method (WL). The second approach was to consider solely target points as liver points, i.e. = P P l t . In other words, we exclusively optimize inside the target. The advantage of this approach is that the temper ature profile inside the target is more homogeneous. We call this method only target method (OT). The two different methods WL and OT are solved with wspgl1. For the linear programming method (LinProg), the choice of the liver point set does not matter. Both suggested choices of liver points achieved similar results. A summary of the methods can be seen in table 1.

Focal spot size
The beam time additionally depends on the focal spot size of the HIFU device, σ z and σ r . In figure 1, the two parameters are shown. Notice that in Figure 1, the focal spot is drawn as a circle, but in fact has an ellipsoidal shape and therefore σ z and σ r are not equal. The focal spot size along the acoustic axis is denoted by σ z , σ r denotes the size perpendicular to the acoustic axis. The smaller the focal spot, the longer it takes to ablate the target. On the other hand, for smaller focal spot size the method is more flexible: When treating a point near the target border, less healthy tissue will be treated with a narrower focal spot. We tried different values of σ z and σ r to see the influence on the treatment time and the efficiency. The ratio between the two deviations σ z and σ r is constant and given by the HIFU device (Cline et al 1994). Here, we assume that the ratio σ σ : z r is equal to 5.

Treatment time
The beam-, changing-and total treatment time of the different methods are shown in figure 2. The duration in figure 2(a) is calculated from the sonication plan x as described in equation (10). It indicates how long the beam time lasts, measured in breathing cycles. It can be seen that the number of breathing cycles used decreases with increasing focal spot size. The reason is that more tissue can be treated if the focal spot size increases. Therefore, less points have to be sonicated to cover the whole target. This yields a shorter beam time. Figure 2(b) shows how many time the HIFU device has to change its position. The number of changes was calculated by equation (11) with c = 1. Notice that when setting c = 1 in equation (11), the changing time and the number of changes are in fact the same. Recall that the continuous motion-compensation approach does not need any change in HIFU transducer position as it can freely steer the beam in all directions.
To estimate the total treatment time, we assume that it takes one breathing cycle to change the HIFU device position along the mechanical axes. The total treatment time for the continuous motion-compensation approach consist only of the beam time. The results can be seen in figure 2(c). WL has a shorter treatment time than OT. However, LinProg achieves the shortest treatment time. The total treatment time of ContMotion is in between of OT and WL.
In figure 2(d), we see the relative treatment time compared to the ContMotion approach. The relative duration shows how the different methods perform compared to the ContMotion solution for different focal spot sizes. We can observe that the relative treatment duration of OT is increasing over the different focal spot sizes, whereas the relative treatment time of LinProg and WL stays almost constant. WL has a relative treatment time around 0.8, LinProg achieves values of about 0.59-0.64. OT has a treatment time of 1.2-1.5 times the sonication time of ContMotion.

Efficiency
In figure 3, we see how the efficiency coefficient given by equation (12) evolves with increasing focal spot size. We can see that the efficiency is almost constant for the different focal spot sizes. The efficiency is the highest for the ContMotion approach. Since ContMotion follows the liver motion, overtreatment can be reduced. However, also with this approach not all overtreatment can be avoided. The reason is that when ablating the border points of the target, also some healthy tissue will be heated. OT performs better than WL, and the LinProg approach has the highest efficiency coefficient of our proposed methods.

Overtreatment
The distribution of the two different kinds of overtreatment is shown in figure 4. The difference between the two methods WL and OT is the distribution of the overtreatment. WL  (11) with c = 1, (c) total treatment time assuming that one change of the HIFU device takes one breathing cycle, (d) relative treatment time compared to the continuous motioncompensation method (ContMotion). induces a higher temperature rise inside the target and treats healthy tissue less. The reason for this behavior are the different strategies. The prior target of OT is to avoid temperature rises above α inside the target. The treatment of healthy tissue is not minimized. Due to the bigger focal spot of the HIFU device with increasing σ z and σ r , more healthy tissue is affected when focusing a point at the target border. This results in an increasing treatment outside the target. LinProg combines the advantages of the two. The treatment of tissue outside the target is almost the same as with WL, and at the same time the overtreatment inside the target is below the level of OT. However, ContMotion achieves the least overtreatment inside and outside the target.

Temperature rise of target points
The dose volume plots for a focal spot size of σ = 5 z mm and σ = 1 r mm can be seen in figure 5. An ideal dose volume plot would have a maximal and minimal temperature rise of 18 °C inside the whole target. We see that OT and LinProg achieve similar dose volume plots. WL induces a temperature rise of above 60 °C, which can cause boiling. However, the ContMotion approach induces the smallest temperature rise inside the target while ablating the whole target. This result corresponds to figure 4(b), where the overtreatment inside the target is shown.

Sagittal and axial slices through the liver
To observe the temperature profile after the HIFU treatment, sagittal and axial slices of the different methods are shown in figures 6 and 7, respectively. In these plots, filled circles represent target points, whereas unfilled circles are healthy tissue. The focal spot deviations are σ = 5 z mm and σ = 1 r mm. In the sagittal slices, we can see that the methods heat up points that lay above the target. This is due to the main motion direction, which is along the z-axis. Note that there is also motion in the xy-plane due to the non-rigid motion, however the main motion is along the z-axis. The LinProg approach induces a lower temperature rise into the points above the target as compared to OT. As a results, less volume will be ablated above the target with the LinProg approach. WL on the other hand induces the highest temperature rise inside the target, this can be seen in figure 8. However, less healthy tissue points are ablated as compared to OT, which can be seen in table 2. The ContMotion approach has the least treatment outside the target, due to the motion-compensation. However, at the right and left hand side of the target more healthy points are treated than below and above the target. This is due to the focal spot shape: more overtreatment occurs along the acoustic axis.
In the axial slices, the heat distribution inside the target can be seen. The ContMotion approach achieves a temperature rise of 18-25 °C inside the target, except for two hotspots, where the temperature is about 35 °C. WL induces the highest temperature rise inside the tumor, We see that all the methods induce the lethal temperature rise to 100% of the target tissue. The optimal solution induces to 100% of the target tissue the lethal temperature rise of 18 °C, and to 0% a temperature rise above 18 °C. and produces two hot spots with a temperature rise of above 60 °C, which can be seen in figure 8. This is problematic, since these temperature rises can cause boiling. OT and LinProg induce a higher temperature rise at the tumor bottom compared to the remainder of the target.
If we look at table 2, we can see how much healthy tissue is ablated in percentage of the target volume. We observe that with the ContMotion approach, healthy tissue of the size of 0.1% of the target volume is ablated. In our proposed methods, more healthy tissue is heated above 18 °C. However, treatment of tissue outside the target can be reduced by turning off the beam when the focal spot position is outside the target.

Minimizing with M f
If minimizing equation (6) with covermatrix M f as defined in equation (8), treatment of tissue outside the target with OT can be reduced. The results can be seen in figure 9. To find appropriate values of R in the definition of S j , we found out that R = 3 mm gives the best results for the focal spot size of σ = 5 z mm, σ = 1 r mm. When choosing R = 0 mm, which means that the beam is immediately turned off when focusing healthy tissue, the method does not improve. Ablation outside and inside the tumor increases (see table 2 and figure 9(e)). The problem when choosing R = 0 mm is that border points of the target have to be heated from points that are inside the target. As a result, these points are heated more. In other words, the energy can not be equally distributed, which results in high temperature rises outside as well as inside the target. However, when choosing R = 3 mm, i.e. allowing the focal beam to heat healthy tissue with a distance of at most 3 mm from the target, better results can be achieved. The energy can be distributed more evenly, which results in lower temperature rises. Hence, ablation of tissue outside the tumor can be reduced (table 2), while the dose volume plot stays almost the same as with OT (figure 9(e)). However, ablation of healthy tissue is still higher compared to the other methods.

Minimizing with M s
The dose volume plot of WL can be further improved by the minimization problem stated in equation (9). We divided the breathing cycle into T = 10 different time steps. In figure 10(c), it can be seen that the dose volume plot of WL can be improved, the maximal temperature rise inside the target is much less compared to WL. However, it is still higher than with the OT and LinProg approach. The ablated healthy tissue is much less compared to WL, which can be seen in table 2. However, the computational effort is larger with this method than computing WL. Further, the treatment time is elongated with this method: the duration is about 1.7 times longer than with WL. This is due to the higher intensities used in the treatment plan.

Discussion
We conclude that LinProg combines the advantages of WL and OT. It minimizes overtreatment inside the target, while avoiding damage on healthy tissue. However, in practice the sonication plan has to be adjusted during treatment. The breathing pattern of the patient can change during treatment which results in a different liver motion. To account for this, the treatment plans have to be adapted regularly. Additionally, the thermo-acoustic parameters of the tissue can not be predicted with a high accuracy (Zhang et al 2009). This requires a fast computation of the treatment plan. The treatment plans computed with the compressed sensing approach (OT & WL) can be calculated in linear time, whereas LinProg has a complexity of ( ) O m 3.5 . However, currently no real-time computation is possible. The presented methods are calculated before the treatment. It will be subject of future work to improve the algorithms towards on-line computation of the treatment plans. As soon as our treatment planning method can be computed in real-time, we will be able to consider non-repetitive breathing patterns of the patient. The goal is then to use a pre-calculated sonication plan. As soon as we observe a change in the motion pattern, the sonication plan will be adapted to the new conditions. However, at the moment this is not yet possible and is subject of future work. The disadvantage of WL is the high temperature rise inside the target which can cause tissue boiling. This problem can be solved by choosing the covermatrix M s and the minimization problem as stated in equation (9). However, the complexity of this problem is higher than the compressed sensing approach. The disadvantage of OT on the other hand is treatment outside the target. By changing the covermatrix M to M f in the compressed sensing formulation, treatment outside the target can be reduced, while preserving linear complexity. The dose volume plots on the other hand stay almost the same.
When comparing our method to ContMotion, we observe that the efficiency is decreased. However, we can achieve shorter treatment plans with the self-scanning approach using the LinProg method.
The idea of self-scanning shows potential to be adapted in other tumor treatment modalities such as proton therapy or classical radiation therapy. The difference between the models is the dose accumulation. Therefore, in our proposed model the covermatrix M has to be adjusted to the particular therapy.
In our simulation, we did not consider diffusion over time and therefore our results are an approximation of the reality. But since the motion-compensation approach was simulated with the same model, we are able to compare the methods under the same conditions. The main thing that will change when including diffusion over time is that after a certain amount of time, the already heated tissue will cool down. As a result, these points have to be heated again, which increases the treatment time. At the border of the target, this effect helps to avoid ablation of healthy tissue. Inside the target on the other side it is probable that the treatment time increases. However, when sonicating the points in an optimal order, especially when sonicating neighboring points, the increase of treatment time can be reduced. Moreover, since in the presented results both the self-scanning as well as the continuous motion-compensating approach were simulated without diffusion, we expect a comparable change in the treatment time. However, it remains to be shown that the self-scanning approach has a shorter treatment time under real conditions.
Another simplification in our model is the assumption that tissue which is heated above 54 °C undergoes necrosis. A more accurate measure on the thermal impact of tissue is the thermal dose equation. According to the thermal dose equation, it might not be enough in some cases to heat tissue above a temperature of 54 °C. However, if tissue is heated for enough time to 54 °C, by the thermal dose equation, necrosis can be induced. Including the thermal dose model in a further study may thus change our estimated treatment time and it remains to be shown that the impact on the self-scanning and the continuous motion-compensating approach is in the same range.
Moreover, there are other limitations, such as the reduced ultrasound window due to ribs and skin heating (Hynynen 2010, Wijlemans et al 2012, Ellis et al 2013. It has to be shown in further studies, how these problems can be solved in detail when using the proposed selfscanning approach. However, the self-scanning approach simplifies beam forming, as compared to the continuous motion compensating approach, we do no longer have to deal with the often contradicting goals of compensation of the motion and avoiding near field heating while at the same time covering the sonication volume.

Conclusion
In this paper, we proposed a novel self-scanning approach for HIFU in liver and showed its feasibility by simulation. We indicated an algorithm to calculate a treatment plan for such a setup. Our method achieved a shorter treatment time than the motion-compensation approach.
However, we simulated our method on a simplified model. It remains to be shown in further studies how the self-scanning method performs under a more realistic model. Moreover, the energy distribution efficiency is reduced. Despite the decreased efficiency, the expected shorter treatment time and the reduced complexity on the beam forming possible with the proposed self-scanning method render the idea of self-scanning attractive.