{'fraction_non_alphanumeric': 0.022767303803248744, 'fraction_numerical': 0.03998956226759737, 'mean_word_length': 4.631888317413666, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of a spring-dashpot system to clinical lung tumor motion data

A spring-dashpot system based on the Voigt model was developed to model the correlation between abdominal respiratory motion and tumor motion during lung radiotherapy. The model was applied to clinical data comprising 52 treatment beams from 10 patients, treated on the Mitsubishi Real-Time Radiation Therapy system, Sapporo, Japan. In Stage 1, model parameters were optimized for individual patients and beams to determine reference values and to investigate how well the model can describe the data. In Stage 2, for each patient the optimal parameters determined for a single beam were applied to data from other beams to investigate whether a beam-specific set of model parameters is sufficient to model tumor motion over a course of treatment. In Stage 1 the baseline root mean square (RMS) residual error for all individually-optimized beam data was 0.90 plus or minus 0.40 mm. In Stage 2, patient-specific model parameters based on a single beam were found to model the tumor position closely, even for irregular beam data, with a mean increase with respect to Stage 1 values in RMS error of 0.37 mm. On average the obtained model output for the tumor position was 95% of the time within an absolute bound of 2.0 mm and 2.6 mm in Stage 1 and 2, respectively. The model was capable of dealing with baseline, amplitude and frequency variations of the input data, as well as phase shifts between the input tumor and output abdominal signals. These results indicate that it may be feasible to collect patient-specific model parameters during or prior to the first treatment, and then retain these for the rest of the treatment period. The model has potential for clinical application during radiotherapy treatment of lung tumors.


Introduction 51
Lung cancer was the most commonly diagnosed cancer worldwide in 2008, with an 52 incidence rate of over 1.5 million new cases, or 12 % of the total new cancer diagnoses 1 . It 53 was also the most common cause of death from cancer, accounting for 17% of all cancer 54 deaths. Radiation therapy is a common modality for lung cancer treatment, however 55 treatment efficacy is limited by the motion of the lungs during respiration 2 , which is primarily 56 driven by diaphragm motion, and to a lesser extent by chest motion. The magnitude of lung 57 tumor motion depends on patient-specific breathing and tumor characteristics, and is usually 58 most pronounced along the superior-inferior (SI) axis, compared to the anterior-posterior 59 (AP) and lateral directions 3 . One study, which included 39 patients treated with the Real-60 Time Radiation Therapy system in Sapporo, Japan, found a median tumor movement of 1.1 61 mm, 2.3 mm and 5.4 mm in lateral, AP and SI direction, respectively 4 . 62 Typically, for lung tumors the clinical target volume (CTV) is enlarged to the internal target 63 volume (ITV) with the intent of ensuring sufficient tumor coverage in the presence of motion 5 . 64 This strategy can result in excessive irradiation of surrounding healthy tissue, or marginal 65 miss of the tumor 6 . The high doses required for tumor control are close to or above the 66 tolerance level for healthy tissues, resulting in increased side effects, or requiring a reduction 67 in dose to the tumor, decreasing tumor control probability (TCP) 7,8 . In this way, 68 compensating for the presence of respiratory motion of lung tumors can improve the 69 therapeutic ratio and thus survival rates 9 . 70 One means to compensate for tumor motion is by means of real-time tumor tracking and 71 motion compensation 8, 10-14 . Motion tracking is complicated by variations in baseline, 72 frequency and oscillatory amplitude and form, despite the overall superficial regularity of 73 respiratory motion. These effects can differ widely between patients but can also vary over 74 the course of treatment of a single patient 7 . Methods to suppress such variability include 75 abdominal compression to reduce tumor mobility, breath control techniques (active or 76 surrogate motion. All of these techniques have some limitations, such as the need for active 78 patient cooperation, consistent ability to maintain total lung capacity 20 , or a lengthened 79 treatment time for gated therapy due to the required beam-off periods. approaches, for which the internal behavior of the physical system is not or only partially 95 known, respectively. These models consider the input and output data, but not the exact 96 physical relationship between them. Various approaches of this form include linear 97 correspondence models 9, 28 , composites of baseline drift, frequency variation, fundamental 98 pattern change and random observation noise 7 , characterization of the motion with a 99 piecewise linear model of defined stages of the breathing cycle 29 , state-based probabilistic 100 models 30 , least-squares parameter estimation and systems identification 3 or adaptive neural 101 networks 31 . These models generally require a large amount of sample data which 102 encompasses the range of possible relationships between the input and output states. Many 103 marker locations, motion dimensions and type of breathing pattern 3, 4, 28 . One of the more 107 promising models is the model by Cervino et al. 26 using an internal surrogate signal of the 108 diaphragm; however the practical benefits of using an external signal warrant further 109 investigation into external abdominal surrogates.  well-modeled by a single Voigt unit acting in the SI direction, with the tumor attached at one 144 end and the abdominal system providing input at the other. It is this pseudo-one dimensional 145 class of tumor motion which we consider in this paper. In the following condensed derivation, 146 all variables and parameters are non-dimensional. 147 The variable x is the tumor displacement from a resting equilibrium position in the SI 148 direction. The Hookean spring has associated parameter ω , while the dashpot, (1) 159 The model will be shown to hold for all patients in this class of tumor motion, with values of 160 the parameter triplet ( ) derived from a "training" subset of patient data, and 161 successfully applied to model the remaining data. 162 163

Numerical approach 164
Equation (1) was rewritten as two coupled first-order ordinary differential equations, which 165 were solved for each set of training data using the MATLAB™ differential equation solver 166 ode45, which is recommended for non-stiff problems 37 . The resulting output data was  In Stage 2 ( 'non-optimized') we take the optimized parameters from a single beam -e.g. 182 Day 1, Beam 1-and use these to model the tumor motion via (1) for that patient's 183 respiratory data obtained from other beams on the same day and consecutive days. The 184 comparison between the predicted tumor position, x, and the clinical data, X, is presented 185 and discussed in sections 3 and 4. 186 187

Clinical data 188
The clinical data consists of 3D internal fiducial-based motion obtained using radiopaque 189 fiducial markers implanted and visualized in real-time using stereoscopic diagnostic x-ray 190 fluoroscopy, collected using a Mitsubishi Real-Time Radiation Therapy system, at the 191 Nippon Telegraph and Telephone Corporation Hospital, Sapporo, Japan 27 . This was 192 with lung cancer at various sites, with a total of 171 treatment beams for up to ten 196 consecutive treatment fractions (day 1-10) and four beam configurations per fraction (beam 197 . The patients were not a random sample of the general lung cancer population, but 198 were selected on the basis of an estimated internal marker motion of more than 10 mm peak 199 to peak (SI, lateral or AP direction). The data show phase shifts of the internal-external data 200 (one signal lagging the other) 32, 33 in the SI direction that are mostly between 100 and 200 201

ms. 202
The data were not acquired in a common coordinate system with absolute spatial co-203 ordinates and therefore they could not be normalized to a common reference point. Spatial 204 baselines for the tumor and abdominal markers were defined by normalizing on the mean of 205 each beam data. 206 For baseline values in Stage 1 the model was initially applied to data from 52 treatment 207 beams from 10 patients: beam 1 data for all days for all patients (as this type of data was 208 available for all patients), as well as all beams for patients 1-3. More detailed analysis was 209 then carried out on representative individual patients: one whose model performance was 210 close to the mean of all patients (Patient 7); one with relatively poor model performance 211 (Patient 8); and one with excellent model performance (Patient 10). 212 In the second stage, for each patient the optimal parameters determined for an arbitrarily 213 chosen single beam were applied to data from other beams and treatment days to 214 investigate whether a beam-specific set of model parameters is adequate to model tumor 215 motion over a course of treatment. 216 (2) Patient 5 had eight fractions, then a further four two months later. The analysis was kept consistent. The data for patient 11 were unreliable and not included. Additionally, we considered the residual vectors for each beam, which appeared visually to 250 have a normal-like distribution, although a more rigorous testing of these vectors revealed 251 that only four beams were sufficiently normal to satisfy the Jarque-Bera 38 test. The bottom 252 row of Figure 4 shows the histograms of three beams. Note that the best fitting probability 253 distribution is scaled to the area of the histograms. The full width at half maximum was 1.6 254 mm, 3.9 mm and 1.2 mm for patients 7, 8 and 10, (mean, poor, excellent), respectively. For clinical applications, a major concern is the proportion of time, which we term the 262 residency, for which the model output is within a small distance of the actual tumor position, 263 i.e. within the treatment margin. The residency of the 52 beams was calculated for bounds 264 on the absolute difference in positions of 1, 2, 3, 4, 5 mm. The residencies from Patients 7, 8 265 and 10, and all 52 beams are shown in Table 1. As can be seen from the table, to be 266 resident 95% of the time a bound of 2 mm is required on average. In the second stage the 267 parameter results obtained from optimizing the Day 1, Beam 1 data for each patient were 268 applied to other beams and days for that patient in lieu of optimizing the data from each 269 When the robustness of the model was examined with respect to the linear relationship 272 between the tumor position and model prediction, the mean correlation coefficient of the 52 273 beams was 0.96, the same value as when each beam was optimized independently. 274 There was a drop in residency in this second stage as can be seen from Table 3. This time a 275 band width of 2.6 mm was required on average, an increase of 0.60 mm, for 95% residency.    One beam was omitted for patient 10 because it contained a substantial amount of missing data, and obviously spurious measurements.
To test the sensitivity to beam and day parameter optimization, non-optimized results were 304 obtained for Patients 7, 8, and 10 using the parameters obtained from optimizing the Day 3, 305 Beam 2 data, a choice made arbitrarily. The change in the average RMS error and the band 306 width required for 95% residency was small as can be seen from Table 5 below, and there 307 was no change in the mean correlation results. This would indicate that different parameter 308 sets obtained from different beams for the same patient give adequate results. 309 310

Irregular beam data 314
The model was capable of dealing with data containing irregularities such as missing data 315 covering several seconds (poor tracking of the abdominal motion), spikes in the amplitude of 316 tumor motion and baseline drift. Figure 5 shows an example of the model applied to data 317 containing some initial baseline drift of the breathing signal, a period of missing data, and an 318 unusual rapid variation in amplitude, or "spike". These data problems are important 319 considerations for modeling lung tumor motion to ensure that errors are not avoidably high. The RMS error for this beam data is 1.09 mm and the correlation coefficient is 0.98. The 339 band width for 95% residency is 2.1 mm. These results indicate that the model deals well 340 with irregularities in the data. Figure 6 below shows the histogram of residuals and the 341 regression plot for the data presented in Figure 5. The tail of negative residuals in this 342 histogram is long relative to those shown in Figure 4. This is not due to poor systematic 343 performance but rather due to the performance in the short recovery period immediately 344 following the spike in Figure 5, lasting less than one second. 345 One of the useful features of the spring-dashpot system is that it managed to successfully 388 model baseline drifts and irregular tumor motion. It was also capable of quickly restoring 389 accurate model output when data was missing as shown in Figure 5. The determination of 390 optimal modeling parameters did not include a systematic search of the global parameter 391 space; a downhill search algorithm was used for efficiency. It remains to be investigated if a 392 stochastic search technique, e.g. simulated annealing or approximate Bayesian 393 computation, is required to remove dependence on the initial parameter estimates, but this 394 will be at the expense of higher computational costs. For consistency the initial starting 395 condition for the search was kept constant but it was found that the nominal values of the 396 optimized parameters diverged if different starting conditions were used (data not shown). 397 However, it was found that the resulting residuals were almost identical, which indicates that 398 the solution space is relatively flat. Therefore it was not considered critical but remains an 399 area of further investigation. 400 For clinical implementation, several extensions to the current modeling process are needed. 401 One is the extension to a full three dimensional system which also considers the minor axes 402 of tumor motion (AP and lateral) for the sake of completeness for the proportion of patients 403 where this motion is non-negligible. Also, ideally the system should be predictive, as there is 404 to the determined change in position. Ultimately, any algorithm has limitations and is 406 dependent on the quality of the input data and therefore the output of the model should be 407 tested and verified with independent means in a clinical setting. 408 409

Conclusion 410
A semi-physical spring-dashpot model to correlate breathing to tumor motion in the superior-411 inferior direction has been presented and applied to clinical tracking data. Optimized model 412 parameters were found to be robust and transferrable to different beams on the same day 413 and consecutive days. Day-to-day variations in the agreement between the model output 414 and the measured data were small, indicating that the model parameters may be determined 415 prior to or on the first day of treatment and then used throughout the course of treatment.