PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Med Biol. Author manuscript; available in PMC 2010 May 11.
Published in final edited form as:
PMCID: PMC2867044
NIHMSID: NIHMS197571

Incorporating system latency associated with real-time target tracking radiotherapy in the dose prediction step

Abstract

System latency introduces geometric errors in the course of real-time target tracking radiotherapy. This effect can be minimized, for example by the use of predictive filters, but cannot be completely avoided. In this work, we present a convolution technique that can incorporate the effect as part of the treatment planning process. The method can be applied independently or in conjunction with the predictive filters to compensate for residual latency effects. The implementation was performed on TrackBeam (Initia Ltd, Israel), a prototype real-time target tracking system assembled and evaluated at our Cancer Institute. For the experimental system settings examined, a Gaussian distribution attributable to the TrackBeam latency was derived with σ = 3.7 mm. The TrackBeam latency, expressed as an average response time, was deduced to be 172 ms. Phantom investigations were further performed to verify the convolution technique. In addition, patient studies involving 4DCT volumes of previously treated lung cancer patients were performed to incorporate the latency effect in the dose prediction step. This also enabled us to effectively quantify the dosimetric and radiobiological impact of the TrackBeam and other higher latency effects on the clinical outcome of a real-time target tracking delivery.

1. Introduction

A concern that frequently emerges in today's modern radiotherapy practice is ‘if the target is where we think it is supposed to be’ (Webb 2006). It is important that this concern is properly addressed in order to minimize setup errors, before the beam is turned on. Other errors to be minimized are the intra-fraction errors, which are induced in the course of a fraction delivery and the inter-fraction errors, which are induced in the course of radiation treatment.

The increased interest in modalities that synchronize the irradiating beam with the mobile target in real time, also known as real-time target tracking radiotherapy, is in part a direct consequence of this concern. Real-time target tracking radiotherapy, henceforth referred to as four-dimensional tracking radiotherapy (4DTRT), is arguably the most efficient validated form of respiratory motion management technique recommended for compensation of organ motion in radiotherapy (Potters et al 2004, Keall et al 2006). A common feature of 4DTRT is an online motion detection and target localization subsystem with an integrated feedback communication loop that relays the detected target trajectory to the beam controller which then initiates the synchronization process. Only a few clinical applications of 4DTRT have been reported. The robotic radiosurgery system was the first clinical implementation of 4DTRT (Adler et al 1999, Ozhasoglu et al 2000, Schweikard et al 2000, Murphy 2002, 2004), while MLC-based 4DTRT using a linear accelerator is still in the early stages of development.

In an ideal implementation, the synchronization of the irradiating beam with the mobile target following target localization is instantaneous. This is not the case in clinical implementation. Rather a finite time exist between the moment when a new target location is detected and the execution of a synchronization event to compensate for the detected target motion. This is called system latency. Latency can be attributed to the time taken by each of the sub-processes involved in tracking, namely motion detection, the calculation of new leaf positions and the time required by the MLC leaves to reach the physical positions sent by the MLC controller (Sharp et al 2004, Sawant et al 2009). Ultimately, system latency leads to geometric error and uncertainties in the delivered dose. There are two ways to properly account for the system latency effect, namely minimize them or quantify and incorporate the effects as part of the treatment planning process.

In the first approach, system latency can be minimized but not avoided leading to residual system latency. Minimization algorithms rely on the ability to predict future target locations and incorporate this information into the current detected state. In other words, the goal of prediction is to be able to make adjustments to the treatment based on measurements at some time t1 that can give a measure of the target position at some future time t2 > t1 with the collimation adjusted for t2 not t1 (Webb 2006). In the second approach, the geometric errors resulting from system latency are modeled as a probability density function (PDF) and via convolution can be incorporated in the dose prediction step.

The first approach, which uses predictive filters, is a widely accepted method. Linear filtering, Kalman filtering (Kalman 1960), neural networks and most recently local regression are among the most commonly investigated predictive methods (Shirato et al 2000, Murphy 2002, Sharp et al 2004, Ruan et al 2007, Ren et al 2007).

The second approach, namely modeling the geometric errors associated with system latency as a PDF and incorporating the effect via convolution in the treatment planning process, has been exploited to a lesser extent even though the use of convolution to incorporate geometric errors in treatment planning is not new. For example, in compensating for geometric errors associated with patient daily setup uncertainties (inter-fraction errors) or errors associated with organ motion (intra-fraction errors), the conventional approach has been the use of bigger margins (setup or motion margins). A widely studied alternative is to model the geometric error as a PDF, and then by convolving with the static dose distribution, one can effectively incorporate their impact in the treatment planning process (Lujan et al 1999, 2003, Keall et al 1999, Beckham et al 2002, Mavroidis et al 2002, Chetty et al 2003, 2004, Naqvi and D'Souza 2005, Vedam et al 2005, Li et al 2008).

This work will focus on the convolution approach, which can be used in conjunction with predictive models to incorporate residual system latency. Following convolution, the plan is re-evaluated and the treatment plan parameters are then manually adjusted or automatically optimized, until the geometric error modified dose distributions meet the clinical objectives (Bortfeld et al 2004, Webb 2006).

Accurate determination of the system latency PDF is crucial for the proper implementation of the convolution technique. It is important to be able to verify the system latency effect independently from the manufacturer's specifications, if provided. In this work, the geometric errors associated with system latency were quantified based on phantom measurements performed on an actual real-time tracking radiotherapy system, TrackBeam (Initia Ltd, Israel), a prototype system assembled and evaluated at our cancer Institute. It should be noted that a linear adaptive filter is inbuilt into the TrackBeam prototype system. Therefore, strictly speaking, residual latency effects are involved. The methodology to be presented in determining system latency or residual system latency is the same and we shall simply refer to the residual system latency associated with the TrackBeam system as system latency or latency for short. The convolution approach was then verified by treatment planning and delivery on a dynamic phantom, and then by comparing the measured dose distribution following tracking with the convolved dose distribution incorporating system latency effects.

Finally, the convolution approach was applied to treatment plans involving six (n = 6) 4DCT volumes of lung cancer patients to: (i) estimate the dosimetric and radiobiological impact of the geometric errors associated with system latency for patients treated on the system and (ii) apply the quantified dosimetric or radiobiological impact to adjust the treatment plan parameters and re-optimize the plan in anticipation of the system latency effects.

2. Materials and methods

Let Dideal be the planned dose distribution based on an ideal real-time target tracking delivery. In other words, Dideal is obtained from a treatment plan assuming that the system latency is zero. Also, let the geometric errors that are associated with the system latency be quantifiable as a probability density function (PDF). Then, we can incorporate the system latency effect in the treatment planning process and obtain the real dose distribution (Dtracking) using the following convolution formula:

Dtracking(x,y,z)=x,y,zDideal(x,y,z)PDF(xx,yy,zz)dxdydz.
(1)

In the 4DTRT system used in this study, tracking can be performed in 2D only, namely the plane perpendicular to the beam direction. We shall consider for simplicity 1D tracking along the superior/inferior (S/I) direction where tumor motion is predominant. The 1D version of equation (1) can then be written as

Dtracking(x,y,z)=zDideal(x,y,z)PDF(zz)dz=DidealPDF.
(2)

Modeling system latency as a PDF is therefore crucial in this work.

2.1. Quantifying system latency

2.1.1. Modeling system latency as a PDF

The experimental setup is shown in figure 1. A prototype 4DTRT, TrackBeam (Initia Ltd, Israel), assembled at our cancer center was used for the measurements. The system comprises a dual-layer micro multi-leaf collimator (DmMLC) attached as an add-on to the Varian Clinac 600C linear accelerator (Varian Medical Systems, Palo Alto, CA). Using an electronic portal imaging device (EPID) and integrated image processing software, the system tracks the trajectory of an implanted fiducial Gold marker in real time, and via a feedback communication loop, it triggers the DmMLC controller to re-position the leaves in response to the detected fiducial marker's location. The EPID uses a high resolution CV-M300 industrial monochrome CCD camera (KAI Corporation, Kanagawa, Japan). The setup was such that a fixed source-to-image plane distance of 150 cm was maintained (see figure 1).

Figure 1
Experimental setup showing the TrackBeam dual-layer micro MLC (DmMLC) and EPID integrated to the CLINAC 600 linac. The QUASAR dynamic phantom set on the linac couch is also shown.

The TrackBeam system latency is inherent in any real-time tracking measurement.

By assuming that equation (1) holds for phantom irradiations, the goal was to design experiments, which would measure the dose distributions Dideal and Dtracking and by deconvolution, derive the PDF. We used a circular reference field of 4 cm diameter to irradiate the QUASAR phantom (Modus Medical Devices Inc., Ontario, Canada) in two instances: (i) deliver 200 MU using a static field on the static phantom. (ii) With the same setup and using a real-time dynamic synchronized beam, deliver 200 MU on a dynamic phantom, where the phantom has been programmed to move in a sinusoidal pattern with 10 mm amplitude and 6 s period. We used EBT gafchromic film (ISP, NJ) inserted in a coronal plane in the phantom for the dose measurements. It should be noted that a radio-opaque marker, a ball bearing of 2 mm diameter was securely inserted into the phantom for the real-time tracking application.

The dose distribution measured in (i) is similar to the ideal case where system latency is zero and constitutes the ideal dose distribution (Didealref). On the other hand, situation (ii) represents a typical real-time target tracking radiotherapy, and therefore the system latency effect is embedded in the measured dose distribution. We shall denote the measured dose distribution in this case as (Dtrackingref). The superscript (ref) has been used to emphasize that the measured dose distributions are for reference fields. Direct substitution of the measured dose distributions in equation (2) gives the following:

Dtrackingref=DidealrefPDF.
(3)

The derivation of the PDF from equation (3) using the measured dose distributions is a deconvolution process. Solving the equation in the space domain is unrealistic in most cases. Using the convolution theorem, namely that the Fourier transform (FT) of a convolution is the pointwise product of FT, the equation can be easily solved in the inverse space (or Fourier space) domain. The solution is

PDF=FFT1(FFT(Dtrackingref)FFT(Didealref)).
(4)

This approach is not suitable for clinical application mainly due to two reasons. First we run into problems when we encounter measurements that lead to zero divisors. Secondly, the deconvolution process is quite sensitive to noise in the measured dose data and the accuracy by which the dose distribution is measured. Therefore, perfectly reasonable attempts at deconvolution can lead to physically unrealistic results such as negative PDF values (Press et al 1992). It is for this reason that this technique was abandoned for another approach.

In the second approach, we assumed a Gaussian-like distribution for the system latency PDF and formulated a least-squares optimization problem by defining a theoretically predicted dose distribution (Dtheoretical) and an experimentally determined dose distribution (Dtracking). Our choice of Gaussian function for the system latency was based on preliminary qualitative observations of the geometric error distribution derived from an offline analysis of the acquired EPID tracking video. In a related work, Vedam et al (2005) confirmed the Gaussian nature of the PDF attributable to system latency by performing the Kolmogorov–Smirnov (K-S) test. This is a goodness-of-fit test used to decide if a test sample comes from a population with a specific distribution. The Gaussian form and the least-squares formulation of the problem is mathematically expressed as follows:

PDF(z,m,σ)=e(zm)22σ2ze(zm)22σ2
(5)
Dtheoretical(z,m,σ)=Didealref(z)PDF(z,m,σ)
(6)
S(m,σ)=z(Dtheoretical(z,m,σ)Dtrackingref(z))2.
(7)

Minimizing equation (7) with respect to the parameters m and σ is a classical nonlinear least-squares optimization problem, which can be solved via many methods. We applied the Levengberg–Marquardt algorithm (LMA) (Levenberg 1944, Marquardt 1963) implemented in MATLAB (The Mathworks Inc., Natick, MA) to solve for m and σ. The user supplied the initial guess (m and σ) and the LMA iteratively determines the best parameters that minimized the sum S (m, σ). The LMA algorithm is robust and converges to a solution even when the initial guess is far off. This was confirmed by repeating the parameter determination process using different initial points. Both m and σ are modeled in the LMA algorithm but moving forward, the dose blurring effect is predominantly a property of σ and hence we shall simply quote the σ value when referring to the latency PDF.

2.1.2. Quantifying system latency from the acquired EPID tracking video

Another way to quantify system latency is through offline analysis of the recorded real-time portal images acquired in the course of the tracking experiment. The EPID video was recorded by the TrackBeam system at a frame rate of 25 Hz. By analyzing the video frame by frame, we can estimate the trajectories of the implanted marker (zm(t)) and the geometric center of the reference field (zGC(t)) as a function of elapsed time (t). The geometric error associated with system latency (εSL) can be obtained directly from the measured trajectories via

εSL(t)=zm(t)zGC(t).
(8)

It should be noted that equation (8) requires separate estimations of the marker and the field geometric center trajectories as a function of elapsed time. It is important to note this distinction from other works that have quantified εSL based on a single measurement, for example using fluoroscopic data of the marker trajectory. In this case, a response time (Δ) of the tracking system is then used to derive a similar form of equation (8) as εSL(t) = zm(t)−zm(t−Δ). In our experience and based on the work by Sawant et al (2009) involving measurements performed on the Calypso four-dimensional localization system (Calypso Medical Technologies), the response time of the tracking system may vary as a function of time, and therefore, our approach is the more accurate.

A probability distribution function of geometric errors associated with system latency can be estimated directly from equation (8) if statistically significant measurements are made. Note that in this case, it will then be possible to derive a Gaussian-like PDF from the distribution which should be comparable to that obtained via the convolution method. For a video acquired over 10 breathing cycles, there can be as many as 1250 frames (images) to analyze therefore requiring a computational approach to the problem.

A software code was developed in SIMULINK (The Mathworks Inc., Natick, MA) that took as input the acquired portal image tracking video and a searchable pattern or template (2D region of interest) and generated as output the location of the pattern in the video as a function of elapsed time. The searching strategy was based on the pattern matching algorithm (PMA) through the application of two-dimensional cross correlations (2DXCorr). Let t (x, y) and I (x, y) be respectively the intensity distribution of a searchable template and an image frame from the portal tracking video. The normalized 2D cross correlation at any given location (u, v) was computed following a formula similar to that used by Lewis (1995):

2DXCorr(u,v)=xy(I(x,y)I¯u,v)(t(xu,yv)t¯)xy(I(x,y)I¯u,v)2xy(t(xu,yv)t¯)2.
(9)

In this equation, t is the average intensity of the template and Īu,v is the average intensity of I (x, y) in the region under the template. A cross-correlation coefficient was computed for each pixel of a given video frame, and the pattern location per frame was considered as the location of the maximum coefficient per frame. Using the normalized form of the 2DXCorr for image application involving template tracking ensures that temporal changes such as the brightness of the image and template have a minimal effect on the PMA algorithm.

Finally, we note that a more intuitive albeit less useful characterization of system latency is by the average time delay (ms). For example, for a tracking beam irradiating the QUASAR dynamic phantom, pre-programmed such that the implanted marker performs sinusoidal motion with period (T), the acquired trajectories would be fitted to sine functions of the form A*sin ((2πt/T) + [var phi]). The time delay (Δt) is then estimated from the phase difference (Δ[var phi]) as Δt=Δϕ2πT.

2.2. Phantom investigations

Phantom investigations were performed to verify the validity of the PDF derived for the system latency and the methodology used to incorporate the system latency effect in the dose prediction step. A treatment plan was generated on the QUASAR dynamic phantom consisting of six 3D conformal fields defined by the gantry angles at 45°, 135°, 180°, 225°, 315° and 360°. The treatment plan was based on a simple spherical target of diameter 4 cm. No further margins such as to account for setup or motion errors were used, and therefore, the PTV coincided with the CTV. The plan was delivered on the QUASAR dynamic phantom programmed to perform 1D sinusoidal target motion in the S/I direction with 10 mm amplitude. The dose distribution for the static fields irradiating a static phantom (Dideal) was considered as the dose distribution that would be delivered in a tracking delivery where real-time synchronized beams irradiate the dynamic phantom under ideal conditions; that is zero system latency. We applied equation (2) to incorporate system latency and estimate the dose that would be delivered under normal tracking conditions, that is non-zero system latency effects (Dtracking). We then performed the tracking delivery and compared the measured dose distribution with Dtracking as well as with Dideal. We used EBT gafchromic film to measure a planar dose distribution. The gamma index analysis (Low et al 1998) was used for the dose comparison with a dose difference tolerance of 3% and distance to agreement (DTA) tolerance of 3 mm.

In describing the TrackBeam experiments leading to the PDF derivation, we have assumed perfect detection of the implanted marker throughout the dose delivery period which is crucial for an uninterrupted tracking process. This is possible only if we can guarantee the implanted marker visibility in the acquired MV video, expressed in terms of contrast-to-noise ratio (CNR) beyond a certain threshold. It tends out that the marker visibility depended on several factors including LINAC dose rate (DR), the video acquisition frame rate (FR), material type and size of fiducial and the phantom or patient thickness. In a preliminary and extensive work, we characterized the TrackBeam system such that given a marker size and phantom thickness, the optimal DR and FR combination that will produce the highest marker CNR was determined. All experiments were performed under these optimal settings generating CNR values that exceeded the threshold for tracking by over 50% therefore guaranteeing an uninterrupted tracking process.

2.3. Patient studies

4DCT image sets from six patients previously treated for lung cancer at our department were used for the study. Each patient's 4DCT volume together with the respiratory signal acquired simultaneously in the course of CT scanning was transferred to the Advantage GE workstation (General Electric, WI) where ten consecutive image sets corresponding to the various phases of the respiratory cycle were further generated. A single clinician delineated the target and organs at risk (OAR) on each image set. To accurately estimate the dose distribution deliverable for the TrackBeam real-time target tracking delivery, it was necessary to produce treatment plans for each of the ten phases of the respiratory cycle. The current implementation of the TrackBeam system is based on 3D conformal beam delivery where a beam aperture is synchronized with the implanted marker trajectory from one respiratory phase to another. The conformal beams were planned on a clinical target volume (CTV) with no further setup, motion or latency margins; therefore, the CTV coincided with the PTV. Each phase plan was performed on the ADAC Pinnacle3® treatment planning system (TPS), version 8.0 m (Philips, Fitchburg, WI) with a dose grid setting of 3 mm × 3 mm × 3 mm. It should be noted that the Pinnacle TPS used in the present study, as well as most currently available commercial TPS, do not calculate dose distributions on 4DCT image data sets. Therefore, the TPS was used to calculate the 3D dose distributions for each of the phases and the results were exported to in-house developed software (Roland et al 2009) where a 4D dose distribution was derived from the multiple plans by applying a validated, non-rigid deformable image registration (DIR) algorithm.

The implemented DIR method used is primarily based on Thirion's diffusion model also known as the ‘demons’ algorithm (Thirion 1995, 1998). The demons DIR algorithm was implemented using the National Library of Medicine Insight Toolkit (ITK), an open source cross-platform C++ software toolkit (Ibanez et al 2003). The result of the DIR application was a deformation field between each image set and a reference image set, which was then used to deform the dose distributions from each of the phases back to the reference phase. In this study, the reference image set was chosen to be that of the end-of exhale phase. A weighted sum of the different dose distributions was computed to constitute the 4D dose distribution:

DR,4DTRTideal=i=110WiτiR[Di,3D].
(10)

Here, DR,4DTRTideal is the 4D dose distribution accumulated on the reference phase image set, Wi is the weight associated with the ith phase, and set to 0.1 as the image sets were equally separated in time, and τiR is the transformation resulting from the DIR application and it takes as its argument a dose distribution from the ith phase and deforms it to the reference phase (R). Note that we have used the subscript 4DTRTideal to emphasize that this dose does not take into account the system latency effect and therefore corresponds to the real-time tracking delivery under ideal conditions.

To account for system latency and therefore incorporate this effect in the dose prediction step, we convolved the dose distribution in equation (10) with the derived PDF, which is attributable to system latency. The resulting dose distribution shall be subscripted with 4DTRTcorrectedSL to underscore that this is the dose distribution that has been corrected for system latency. The two distributions derived for all the patient plans were then evaluated considering dosimetric and radiobiological analyses.

2.3.1. Dosimetric analysis

The impact of system latency on the dose delivery was evaluated in terms of its effect on normal tissue sparing and target coverage. For normal tissue sparing, we computed and compared the mean lung dose (MLD) and lung V20 for the 4DTRTideal and 4DTRTcorrectedSL dose distributions. Lung V20 is the volume of lung irradiated by at least 20 Gy. In order to determine the loss of target coverage adequacy due to system latency effects, we computed and compared the minimum dose irradiating 95% of the planning target volume (PTV D95).

It should be noted that the capability to incorporate system latency in the dose prediction step enables us to first generate a dose distribution assuming an ideal situation, and then by comparing the system latency corrected dose distribution with the ideal dose distribution, we can manually adjust the plan parameters until the original treatment plan clinical objectives are met.

2.3.2. Radiobiological analysis

A more sensitive assessment of the effect of system latency on the treatment outcome can be achieved by performing a complete radiobiological analysis, which applies models that take into account both the physical dose distribution and the varying organ radiobiological dose-response characteristics. The models used to describe tumor control or normal tissue complication probability, P(D), were the linear quadratic Poisson and the relative seriality models and have been thoroughly described elsewhere (Mavroidis et al 1997, Lind et al 1999). A short summary and application of these models is presented here:

P(D)=exp{N0e(D/D50)(eγlnln2)}=exp{eeγαndβnd2}.
(11)

In this equation, D is the uniform dose, d = D/n is the dose per fraction and n is the number of fractions applied. D50 is the dose that causes response to 50% of the patients, γ is the steepness of the dose–response curve and α and β are the fractionation parameters of the model and account for the early and late effects expected. Both D50 and γ depend on N0, the initial number of the clonogenic cells for tumors or functional subunits for healthy tissues. These parameters are specific for each organ and injury type and are derived from clinical data.

Note that equation (11) above applies to a uniform dose. For heterogeneous dose distributions (D) often encountered in clinical practice, we can use equation (11) to derive two useful quantities, the tumor control probability, also known as the probability of benefit, PB(D), and the normal tissue complication probability also referred to as the probability of injury, PI(D). This can be achieved by applying the P(D) model to multiple subvolumes (Δj) of the region of interest where each subvolume is irradiated with a uniform dose (Dj). The result is a P(Dj) for each subvolume (Δj) from which PB(D) or PI(D) can be deduced by considering the joint probability.

The radiobiological parameters evaluated for the treatment plans 4DTRTideal and 4DTRTcorrectedSL were the biological effective uniform dose (BEUD) and the complication-free tumor control probability (P+) derivable from PB(D) and PI(D). BEUD is the uniform dose ([D with double overline]) that causes the same tumor control or normal tissue complication probability as the real dose distribution in a complex target or normal tissue case (Lind et al 1999, Mavroidis et al 2000, 2001). The target BEUD ([D with double overline]B) or normal tissue BEUD ([D with double overline]I) was derived from the expressions of PB(D) and PI(D), respectively, by a straightforward application of the definition of BEUD. For example, to calculate the [D with double overline]B, we set up the equation PB([D with double overline]) = PB(D), and then solve for [D with double overline].

The P+ objective on the other hand is a very effective radiobiological index that combines a treatment plan's advantages in terms of tumor control and disadvantages in terms of normal tissue complications (Mavroidis et al 2001). The general expression for P+ (Kallman et al 1992, Mavroidis et al 1997) is given by

P+=PBPI.
(12)

By computing and comparing the [D with double overline] and P+ parameters of the 4DTRTideal and 4DTRTcorrectedSL plans generated on each patient's 4DCT volume, we estimated the effect of system latency on the clinical outcome of the treatment delivery.

3. Results

3.1. Quantifying system latency

Applying the deconvolution technique and LMA algorithm, the standard deviation (σ) of the Gaussian PDF attributable to the system latency (figure 2) for the TrackBeam experimental system was estimated to be 0.372 cm. Furthermore, application of the PMA algorithm on an acquired tracking video of the trajectory of an implanted fiducial marker irradiated by a diamond-shaped field resulted in an estimated average system latency of 172 ms. A 38 × 37 pixel searchable template containing the marker fluence projection was used as the fiducial marker template, while a surrogate template of similar size containing the edge of the diamond field was used to track the field geometric center (figure 3(a)). Finally, a probability distribution function of the geometric error between the marker location and the geometric center location was derived from a histogram of the variation of the marker and irradiating field geometric center locations (figure 3(c)). The histogram was based on a 19.5 s video acquired over 3.25 cycles at a frame rate of 25 Hz, thus generating approximately 488 images (frames), which were analyzed.

Figure 2
System latency modeled as PDF via deconvolution and the LMA algorithm. (a) Profiles for static beams, one involving a static target and the other involving a dynamic target (no tracking) thus producing the motion blurring effect. In (b), tracking reduces ...
Figure 3
Quantifying system latency from the acquired tracking portal video. (a) A frame from the video with marker ROI and an edge ROI used as marker and geometric center templates, respectively. (b) The acquired trajectories as a function of elapsed time. (c) ...

3.2. Phantom investigation

We measured dose distributions corresponding to three delivery scenarios, namely Dideal for the treatment plan (static beams) delivered on a static phantom, Dtracking for the treatment plan (realtime synchronized beams) delivered on a dynamic phantom and Dno-tracking for the treatment plan (static beams) delivered on a dynamic phantom. The delivered dose distributions were measured using EBT Gafchromic film, and a gamma analysis was performed between the distributions with the RIT software.

The percentage of pixels satisfying the set gamma tolerance criteria and therefore with a computed gamma index of ≤ 1 was 89% for the Dideal and Dno–tracking comparison (figure 4(a)) and 93% for the Dideal and Dtracking comparison (figure 4(b)). The 11% failure in the first comparison is attributable to motion blurring effects, while the 7% failure in the second comparison is due to system latency effects. By convolving Dideal with the derived system latency PDF, we obtained a dose distribution that has been corrected for system latency DcorrectedSL. Comparing DcorrectedSL with Dtracking resulted in 98% of pixels satisfying the gamma tolerance criteria (figure 4(c)). It can be observed from the gamma map that the 2% failure occurs mainly at the location of the high-density radio-opaque implanted fiducial marker. The marker produced a very high heterogeneous media, and the convolution algorithm does not perform well for high inhomogeneous media.

Figure 4
Gamma pass/fail map for (a) Dideal versus Dno–tracking—89% pass. (b) Dideal versus Dtracking—93% pass and (c) DcorrectedSL versus D ...

3.3. Patient application

The dosimetric impact of system latency on dose delivery for a simple case considering 1D tracking in the SI direction was estimated for all six patients. Dosimetric parameters computed for dose distributions assuming an ideal tracking situation (4DTRTideal) and the corresponding dose distribution corrected for system latency (4DTRTcorrectedSL) showed very small differences. For example, the average observed variation in lung V20, MLD and PTV D95 was (0.4 ± 0.3)%, (0.5 ± 0.4)% and (1.9 ± 0.6)%, respectively (figure 5). Consistent with this trend were the radiobiologically quantified parameters where the observed variation in the target BEUD, normal tissue BEUD and the complication-free tumor control probability, P+ was (0.5 ± 0.2)%, (0.6 ± 0.5)% and (0.6 ± 0.5)%, respectively (figure 5).

Figure 5
Dosimetric and radiobiological analyses of the effect of system latency on delivered dose distribution for the TrackBeam system.

Overall, target coverage was affected more than normal tissue sparing. It was observed that increasing the total delivered monitor units (MU) by approximately 1 to 3% compensated for the anticipated loss of coverage while keeping the normal tissue parameters and P+ within acceptable limits. It should be noted that even without the compensation, our results show that the system latency effect of the experimental tracking system will have a negligible impact on the delivered dose or on the clinical outcome. This may not be the case for systems with higher latency effects. To illustrate this, we repeated the analysis for two other latencies (σ = 0.777 cm and σ = 1.172 cm). Two patients were considered in this case (P1, P2). The dosimetric and radiobiological impact were considerable in this case, for example in figure 6 a 4D DVH plot shows the dosimetric impact for patient P2 and the probability versus BEUD plot shows the radiobiological impact on the same patient.

Figure 6
Illustrating the dosimetric and radiobiological impact of system latency for patient P2. In the radiobiological analysis, the probability versus BEUD plot is used for plan assessments in the same way as DVHs are used in physical dosimetry.

As one would expect, higher system latencies affected the outcome of the treatment plan more. This is illustrated in figure 7. The ultimate goal of incorporating system latency in the dose prediction step is the proper adjustment of the treatment plan parameters in order to compensate for the degradation. This was illustrated for patients P1 and P2. For simplicity, we aimed at an adjustment that would restore PTV D95 to within 0.5% of the ideal and ensure that lung V20 and MLD would not deteriorate by more than 5% following the adjustments. Note that adjustments can be as simple as renormalizing the dose distribution, that is rescaling the delivered monitor units, or as extended as re-optimizing the entire plan. We used the simple MU rescaling for this illustration and the results are summarized in table 1. In this table, it can be observed that the simple MU adjustment compensated for PTV D95 while keeping the OAR parameters at acceptable limits for the system with σ = 0.777 cm only. This was not the case with the higher system latency values applied, which means that in such cases more extended parameter adjustments such as re-weighting the beams, using non-coplanar fields or re-optimizing the whole plan will be necessary.

Figure 7
Loss of PTV coverage shown for three system latencies, σ = 0.372 cm; σ = 0.777 cm and σ = 1.172 cm.
Table 1
Incorporating system latency in the treatment planning process—a simple illustration for patients P1 and P2. The goal of the adjustment was to increase the delivered MU to offset the decrease in D95 and ensure that lung V20 and MLD do not deteriorate ...

Finally, it is important to note that the radiobiological and dosimetric impact of the latency effect were quantified based on a simple situation involving 1D motion and a population-based target trajectory. In general, the latency PDF can depend on several factors including target motion amplitude, breathing period, the LINAC dose rate, the geometry of the reference beam aperture and the general form of the tumor trajectory. Phantom investigations showed negligible dependence of the PDF on amplitude, breathing period and LINAC dose rate within the ranges [7.5 mm, 30 mm], [3 s, 7 s], and [240 MU min−1, 400 MU min−1] respectively. We also observed negligible dependence on the beam aperture used for the reference field; circular versus diamond shaped fields. The dependence on the form of the target trajectory was not investigated, but for the patients studied, the simple sinusoidal curve was sufficient. Overall, our methodology can be performed for each patient study prior to treatment where each patient's specific tumor trajectory is programmed into the dynamic phantom to acquire the reference field measurements. In this case, the derived PDF is patient specific and we do not have to pay further attention to variation in the various factors mentioned.

4. Discussions

Many studies have demonstrated that prediction filters lead to improvement in the target localization. For example Sharp et al (2004) investigated the localization error for systems using predictive models such as linear prediction, neural networks and Kalman filtering against the errors for systems which used no prediction and found significantly reduced root mean square error in the former case. Their work was based on actual tumor trajectories measured for 14 lung cancer patients. In another study based on diaphragm-motion fluoroscopy recordings from five patients and using the linear adaptive filter as the predictive model, Vedam et al (2005) showed that the dosimetric impact of geometric errors due to motion prediction was smaller than the dosimetric impact of latency, where response times from 0 to 0.6 s were considered. Therefore, the integration of prediction filters in real-time target tracking radiotherapy systems is a widely acceptable approach to minimize the effect of system latency. Current and future efforts seem to focus on improving the implemented algorithms and produce more robust filters.

The goal of our study was to examine an approach where system latency can be modeled as a PDF and then incorporated as part of the treatment planning process via convolution. The technique is particularly useful when applied in concert with predictive filters to incorporate residual system latency effects in the treatment planning process. Accurate modeling of the system latency effect as a PDF was crucial for the proper implementation of the convolution method. A comprehensive method was presented that uses reference field measurements on a dynamic phantom to extract the system latency PDF based on the deconvolution process. This technique was applied to derive the standard deviation of the Gaussian PDF associated with the system latency of the TrackBeam system, which was 0.372 cm. Another approach was also presented to extract the implanted marker and irradiating field geometric center trajectories based on a pattern matching algorithm, applicable in situations where offline analysis of the acquired EPID real-time tracking video was possible. The estimated system latency, expressed as the average time delay, was derived in this analysis to be 172 ms.

The capability to perform real-time tracking delivery provided an opportunity to verify the convolution approach based on phantom measurements. A six-field treatment plan was delivered on a dynamic phantom under normal real-time tracking conditions using the TrackBeam system. If the convolution approach was plausible, the measured dose following tracking should be in good agreement with the convolved dose that has been corrected for system latency. Using the gamma index analysis, we verified that this was indeed the case since 98% of the pixels satisfied the gamma tolerance criteria of 3% dose difference and 3 mm DTA. Comparing the measured dose distribution following tracking with the planned dose distribution assuming ideal system latency of zero, the gamma pass rate was reduced to 93%. The 5% decrease was thus attributable to system latency effects.

Following the phantom verification, patient studies were performed to quantify the dosimetric and radiobiological impact of the effect of system latency if the patients were treated on the TrackBeam system under the experimental settings. For simplicity, we only considered one-dimensional tracking in the direction with predominant lung motion, namely the SI axis. In the dosimetric analysis, the PTV D95 for target coverage adequacy and lung V20 and MLD for normal tissue sparing were computed and compared. On average, a 2% loss of tumor coverage was observed and the target coverage loss was consistent across all the six patients. Normal tissue sparing was on average within 0.5%. However, in this case, the observed variation was an improvement in OAR sparing for some patients and degradation for others. This is important as this means that the system latency blurring effects can improve the normal tissue sparing for some plans, understandably so as the blurring may result in the distribution of high doses to a larger volume thus reducing the effective average dose. Consistent with the dosimetric analyses were the radiobiological assessments in which the BEUD values that were computed for the target and OAR as well as the complication free tumor control probability (P+) were affected on average by less than a percentage point.

Quantifying the impact of system latency on the delivered dose distribution via the convolution approach was useful as this could be incorporated as part of the treatment planning process. The necessary action to compensate for degradation caused by system latency could be as simple as adjusting some treatment plan parameters in anticipation of the effect or in complex cases it could involve the re-optimization of the entire treatment plan (Li and Xing 2000, Bortfeld et al 2004, Webb 2006). For simplicity, we focused on the compensation requiring treatment plan parameter adjustments. For the six patient plans studied, which were based on treatments using the TrackBeam system, it was observed that we can recover the system latency degradation by simply increasing the planned MU by 40–115 MUs. To illustrate this further, we chose two higher levels of system latency characterized by σ = 0.777 cm and σ = 1.172 cm. It was observed that the system latency-induced degradation could be compensated for the former system by simply increasing the MUs but this was not the case in the latter system. This underscored the fact that simple treatment plan parameter adjustments may not be enough to compensate for system latency and some situations may require more complex adjustments including re-optimization.

Finally, it is important to note some of the assumptions that allowed us to apply the convolution approach in incorporating geometric errors in the dose prediction step: (i) shift invariance, namely the assumption that the dose cloud is unaffected by geometrical shifts in the medium being irradiated. Surface curvature and media heterogeneities undermine this assumption. (ii) Infinite number of fractions. (iii) A constant dose delivered per fraction. The limitations of the convolution method resulting from these assumptions have been studied by Craig et al (2003b), (2003a) and Song et al (2004) respectively. Some authors have proposed a fluence convolution method to circumvent the limitation due to the shift-invariance assumption (Beckham et al 2002, Chetty et al 2003). In spite of these limitations, the convolution approach appears in general to be a time-efficient, analytical approach to incorporate the effects of random geometric uncertainties in the treatment planning process.

5. Conclusions

A convolution technique to incorporate geometric errors associated with system latency was presented and successfully implemented in the TrackBeam real-time target tracking system. Accurate modeling of the system latency as a PDF was crucial for the proper application of the technique and this was done via deconvolution involving reference field measurements on a dynamic phantom. The technique allows us to accurately quantify the geometric errors associated with either system or residual latency and therefore to determine their impact on the delivered treatment. Furthermore, the convolution approach can be used in conjunction with the predictive filters to incorporate residual latency effects as part of the treatment planning process.

Acknowledgments

This work was partially funded by Initia Corporation and by National Institutes of Health/National Library of Medicine grant (1 R01 LM009362). The authors would like to thank Yaxi Liu, PhD, and Virginia Goytia, MD, for their valuable contribution.

References

  • Adler JR, Jr, Murphy MJ, Chang SD, Hancock SL. Image-guided robotic radiosurgery. Neurosurgery. 1999;44:1299–306. [PubMed]
  • Beckham WA, Keall PJ, Siebers JV. A fluence-convolution method to calculate radiation therapy dose distributions that incorporate random set-up error. Phys Med Biol. 2002;47:3465–73. [PubMed]
  • Bortfeld T, Jiang SB, Rietzel E. Effects of motion on the total dose distribution. Semin Radiat Oncol. 2004;14:41–51. [PubMed]
  • Chetty IJ, Rosu M, McShan DL, Fraass BA, Balter JM, Ten Haken RK. Accounting for center-of-mass target motion using convolution methods in Monte Carlo-based dose calculations of the lung. Med Phys. 2004;31:925–32. [PubMed]
  • Chetty IJ, Rosu M, Tyagi N, Marsh LH, McShan DL, Balter JM, Fraass BA, Ten Haken RK. A fluence convolution method to account for respiratory motion in three-dimensional dose calculations of the liver: a Monte Carlo study. Med Phys. 2003;30:1776–80. [PubMed]
  • Craig T, Battista J, Van Dyk J. Limitations of a convolution method for modeling geometric uncertainties in radiation therapy: II. The effect of a finite number of fractions. Med Phys. 2003a;30:2012–20. [PubMed]
  • Craig T, Battista J, Van Dyk J. Limitations of a convolution method for modeling geometric uncertainties in radiation therapy: I. The effect of shift invariance. Med Phys. 2003b;30:2001–11. [PubMed]
  • Ibanez L, Schroeder W, Ng L. ITK Software Guide. New York: Kitware Inc.; 2003.
  • Kalman RE. A new approach to linear filtering and prediction problems. Trans ASME D. 1960;82:35–45.
  • Kallman P, Agren A, Brahme A. Tumor and normal tissue responses to fractionated non uniform dose delivery. Int J Radiat Biol. 1992;62:249–62. [PubMed]
  • Keall PJ, Beckham WA, Booth JT, Zavgorodni SF, Oppelaar M. A method to predict the effect of organ motion and set-up variations on treatment plans. Australas Phys Eng Sci Med. 1999;22:48–52. [PubMed]
  • Keall PJ, et al. The management of respiratory motion in radiation oncology report of AAPM Task Group 76. Med Phys. 2006;33:3874–900. [PubMed]
  • Levenberg K. A method for the solution of certain non-linear problems in least squares. Q Appl Math. 1944;2:164–8.
  • Lewis JP. Fast template matching. Vis Interface. 1995;95:120–3.
  • Li G, Citrin D, Camphausen K, Mueller B, Burman C, Mychalczak B, Miller R, Song Y. Advances in 4D medical imaging and 4D radiation therapy. Technol Cancer Res Treat. 2008;7:67–82. [PubMed]
  • Li JG, Xing L. Inverse planning incorporating organ motion. Med Phys. 2000;27:1573–8. [PubMed]
  • Lind BK, Mavroidis P, Hyödynmaa S, Kappas C. Optimization of the dose level for a given treatment plan to maximize the complication free tumor cure. Acta Oncol. 1999;38:787–98. [PubMed]
  • Low DA, Harms WB, Mutic S, Purdy JA. A technique for the quantitative evaluation of dose distributions. Med Phys. 1998;25:656–61. [PubMed]
  • Lujan AE, Balter JM, Ten Haken RK. A method for incorporating organ motion due to breathing into 3D dose calculations in the liver: sensitivity to variations in motion. Med Phys. 2003;10:2643–9. [PubMed]
  • Lujan AE, Larsen EW, Balter JM, Ten Haken RK. A method for incorporating organ motion due to breathing into 3D dose calculations. Med Phys. 1999;26:715–20. [PubMed]
  • Marquardt D. An algorithm for least-squares estimation of nonlinear parameters. SIAM J Appl Math. 1963;11:431–41.
  • Mavroidis P, et al. Comparison of conformal radiation therapy techniques within the dynamic radiotherapy project ‘Dynarad’ Phys Med Biol. 2000;45:2459–81. [PubMed]
  • Mavroidis P, et al. Effects of positioning uncertainty and breathing on dose delivery and radiation pneumonitis prediction in breast cancer. Acta Oncol. 2002;41:471–85. [PubMed]
  • Mavroidis P, Kappas C, Lind BK. A computer program for evaluating the probability of complication-free tumor control incorporated in a commercial treatment planning system. J BUON. 1997;3:257–64.
  • Mavroidis P, Lind BK, Brahme A. Biologically effective uniform dose for specification, report and comparison of dose response relations and treatment plans. Phys Med Biol. 2001;46:2607–30. [PubMed]
  • Murphy MJ. Fiducial-based targeting accuracy for external-beam radiotherapy. Med Phys. 2002;29:334–44. [PubMed]
  • Murphy MJ. Tracking moving organs in real time. Semin Radiat Oncol. 2004;14:91–100. [PubMed]
  • Naqvi SA, D'Souza WD. A stochastic convolution/superposition method with isocenter sampling to evaluate intrafraction motion effects in IMRT. Med Phys. 2005;32:1156–63. [PubMed]
  • Ozhasoglu C, et al. Real-time tracking of the tumor volume in precision radiotherapy and body radiosurgery-a novel approach to compensate for respiratory motion. In: Lemke HU, Vannier MW, Inamura K, Farman AG, Doi K, editors. Proc 14th Int Conf on Computer Assisted Radiology and Surgery (CARS 2000) San Francisco, CA, USA: 2000. pp. 691–6.
  • Potters L, et al. Int J Radiat Oncol Biol Phys. Vol. 60. 2004. American Society for Therapeutic Radiology and Oncology and American College of Radiology Practice Guideline for the performance of stereotactic body radiation therapy; pp. 1026–103. [PubMed]
  • Press WH, et al. Numerical Recipes in C: The Art of Scientific Computing. 2nd. New York: Cambridge University Press; 1992. pp. 542–3.
  • Ren Q, Nishioka S, Shirato H, Berbeco RI. Adaptive prediction of respiratory motion for motion compensation radiotherapy. Phys Med Biol. 2007;52:6651–61. [PubMed]
  • Roland T, Mavroidis P, Gutierrez A, Goytia V, Papanikolaou N. A radiobiological analysis of the effect of 3D versus 4D image-based planning in lung cancer radiotherapy. Phys Med Biol. 2009;54:5509–23. [PubMed]
  • Ruan D, Fessker JA, Balter JM. Real-time prediction of respiratory motion based on local regression methods. Phys Med Biol. 2007;52:7137–52. [PubMed]
  • Sawant A, et al. Towards submillimeter accuracy in the management of intrafaction motion: the integration of real-time internal position monitoring and multileaf collimator target tracking. Int J Radiat Oncol Biol Phys. 2009;74:575–82. [PubMed]
  • Schweikard A, Glosser G, Bodduluri M, Murphy MJ, Adler JR. Robotic motion compensation for respiratory movement during radiosurgery. Comput Aided Surg. 2000;5:263–77. [PubMed]
  • Sharp GC, Jiang SB, Shimizu S, Shirato H. Prediction of respiratory tumour motion for real-time image-guided radiotherapy. Phys Med Biol. 2004;49:425–40. [PubMed]
  • Shirato H, et al. Physical aspects of a real-time tumor-tracking system for gated radiotherapy. Int J Radiat Oncol Biol Phys. 2000;48:1187–95. [PubMed]
  • Song W, Battista J, Van Dyk J. Limitations of a convolution method for modeling geometric uncertainties in radiation therapy: the radiobiological dose-per-fraction effect. Med Phys. 2004;31:3034–45. [PubMed]
  • Thirion JP. Technical Report, Research Report RR-2547 Epidure Project. INRIA Sophia; 1995. Fast non-rigid matching of 3D medical image.
  • Thirion JP. Image matching as a diffusion process: an analogy with Maxwell's demons. Med Image Anal. 1998;2:243–60. [PubMed]
  • Vedam S, Docef A, Fix M, Murphy M, Keall P. Dosimetric impact of geometric errors due to respiratory motion prediction on dynamic multileaf collimator-based four-dimensional radiation delivery. Med Phys. 2005;32:1607–20. [PubMed]
  • Webb S. Motion effects in (intensity modulated) radiation therapy: a review. Phys Med Biol. 2006;51:R403–25. [PubMed]