Stima della correlazione interno/esterno nella radioterapia ad alta precisione per la compensazione dei movimenti respiratori

The goal of this work is to compare different methods to estimate the internal/external correlation during the radiation treatment of extra-cranial tumors. The problem is generated by the fact that during radio-therapy treatment of tumors located in the abdominal and thoracic regions, organ motion (i.e. the movement of the organs caused by physiological processes) negatively affects the treatment efficacy. There are 2 different kinds of organ motion: 1. Inter-fractional organ motion: movements that can be only detected between one fraction and the following ones. This type of organ motion can be generated by phenomena like patient weight loss, tumor shrinkage or growth, modifications due to bladder or intestine physiological processes. 2. Intra-fractional organ motion: it is generated by processes that are detectable during a single therapy fraction as well, i.e. heart or aorta pulsation or breathing movement. In our research we have focused our attention on intra-fractional organ motion, in particular breathing motion, and we have briefly summarized some of the existing techniques to quantify and compensate tumor movements: Breath Hold, Respiratory Gating or Tumor Tracking. In particular, tumor tracking techniques try to track the tumor position in real-time and dynamically re-position the radiation beam. In order to do that, it’s necessary to follow the tumor position in real time and develop a predictive model that can anticipate tumor motion; this is necessary to take into account the latencies in beam positioning. An example of such approach is the SynchronyTM Respiratory Tracking System, a sub-system of the treatment robotic device Cyberknife (Accuray, Inc., Sunnyvale CA) used for the treatment of extra-cranial tumors, which move due to respiration. The RTS system combines the information coming from an infra-red positioning system, that monitors patient abdominal movements, with the information coming from a couple of X-ray imagers, that provide accurate information on the internal tumor position; this is possible thanks to gold markers which are implanted near the lesion, and which are clearly visible in the X-ray images. The infrared emitters are placed on the chest and abdomen of the patient, and the infrared tracking system records the emitters’ movements at very high frequency. Even if the system provides a real-time information about the breathing movement of the patient, this could be not enough to determine the exact position of the lesion, because the external markers alone may reflect improperly the internal displacement generated by breathing motion. Consequently, X-ray images are acquired at the beginning of each fraction; these images show both the internal and the external markers and are therefore used to compute the existing correlation between internal and external movement. Due to the fact that breathing movements can severely change during the fraction, additional X-ray images are acquired during treatment and are used to update the model parameters. The RTS system builds the correlation using a linear model, if a simple correspondence exists between internal and external markers, or a quadratic model otherwise. Following the RTS system analysis, we have considered different methods for the computation of the internal/external correlation: linear and quadratic filters, Kalman filters, ARMA models, state models and neural networks. We have then decided to implement the state models and the neural networks for our application. The state model infers tumor position by using information about the current coordinates of the external surrogates and a delayed observation where τ is the time delay. In this way it is possible to detect the first order breathing signal dynamics. Once analyzed the equations that control the model, we have implemented different models and analyzed their performance varying the model parameters, such as the input signal provided to the model, the number of samples to be used for training, the control frequency of the output, the error threshold beyond which the correspondence is considered invalid and the time delay τ. Starting from a preliminary evaluation, we chose to implement the following 3 models (I=input, O=output): 1. Quadratic model with I: first principal component of each marker coordinate, O: the modeler (x, y, z) coordinates; 2. Quadratic model with linear control and I: first principal component of marker 1 coordinates, O: the modeler (x, y, z) coordinates; 3. Quadratic model with I: first principal component of marker 1 coordinates, O: the modeler (x, y, z) coordinates; The 3 models have been implemented by using quadratic equations, except for model number 2, where a linear model has been added as well; in this case, the linear or quadratic model are used according to specific positions of the external surrogate. The other parameters were set as follows, for all models: - Initial training samples number: 500; - Intra-fraction training samples number: 500; - Error threshold (above which a re-training is required): 3 mm; - Output control frequency: each 500 samples (i.e. every 20 seconds) The procedure has been similar for the neural networks as well: we have tested different types of networks, changing the input, the number of samples to be used for training, the network output check frequency, the error threshold, the number of samples used for intra-fraction training, the transfer function, the number of hidden layers and the number of neurons per layer. The networks that provided the best results and that have been chosen for our analysis are the following: 1. Network with 3 hidden layers and 1 output layer, 2 neurons for each hidden layer (1 for the output layer), sigmoidal transfer function for the hidden layers and linear between the last layer and the output, quasi-Newton learning algorithm, input composed by all marker coordinates (the input matrix dimension is 9xN, where N represents the fraction duration in number of samples acquired at 26 Hz) and output consisting of the 3 coordinates of the lesion (dimension 3xN); the network was provided with 1000 samples for initial training, 200 samples for the re-training of the existing network (if the error is between 3 and 5 mm) and 1000 samples for the training of a new network (if the error exceeds the 5mm threshold); 2. Similar as above, but in this case the network was provided with 500 samples instead of 1000 for initial training and 500 samples for the training of the new network (in case the error exceeds the 5 mm threshold); 3. In this case we have implemented 3 different networks, each one with 3 hidden layers and 1 output layer, 2 neurons for each hidden layer (1 for the output layer), sigmoidal transfer function for hidden layers and linear between the last layer and the output, quasi-Newton learning algorithm. The input for each network consisted of one marker coordinates (that is, the input matrix dimension is 3xN) whereas the output consisted of the 3 coordinates of the lesion (3xN dimension); the final output is obtained by computing the average of the 3 estimations calculated by the 3 different networks. The network was provided with 1000 samples for initial training, 200 samples for the re-training of the existing network (if the error is between 3 and 5 mm) and 1000 samples for the training of a new network (if the error exceeds the 5mm threshold); 4. Similar as above, but in this case the network was provided with 500 samples instead of 1000 for initial training and 500 samples for the training of the new network (in case the error exceeds the 5mm threshold); The 3 models and the 4 networks described above have been then used to estimate the correspondence in 20 fractions coming from 16 patients of the Georgetown University Medical Center (Washington, DC) database. Among those, there are 10 cases indicated as Control (C), which are characterized by a low standard deviation breathing signal and 10 cases indicated as worst (W) where the standard deviation is significantly higher. The results can be summarized as follows (where the models and the networks are referred using the numbers listed above):  The results obtained, in terms of mean square error, are worse for patients in the W class with all the analyzed methods;  Among the implemented state models, model number 1 has provided the most accurate output estimation in all the control cases; however, it features consistent errors in almost all the patients belonging to the W class. For these patients, the best results were given by model number 2, which is the only one that uses a mixed quadratic and a linear method;  As far as neural networks are concerned, network number 2 has computed the best estimation for most of C patients. For the W class patients the situation is different, and there is no real prevalence of one of the investigated networks. What can be seen in those cases is that the 2 networks provided with 500 samples (for initial training and for the re-training during treatment) have smaller errors if compared to the corresponding networks provided with 1000 samples; in addition to this, they have also used a smaller portion of the signal for the intra-fraction training;  By comparing the models and the networks that have provided the best results, it can be noticed that neural networks obtained better performance if compared to state models, especially for the W class patients, where in some cases the mean square error was one order magnitude smaller. By observing the breathing signal standard deviation it can be inferred that:  If the standard deviation is smaller than 3 mm the use of network 2 or model 1 is more or less equivalent;  For those patients where the standard deviation is in the range of 3 – 5 mm, networks 2 and 4 and model 2 have provided the best results;  If the standard deviation is greater than 5 mm, then the best results have been obtained with networks 2 and 4 in more or less all cases. Possible future developments may include:  The implementation of a linear control for state model number 1;  The simultaneous prediction with multiple models/networks with on-line computation of prediction errors, in order to implement error-based selection criteria;  To select the architecture parameters for ANN models by implementing optimization criteria for the initialization conditions.

Lo scopo di questo lavoro è stato quello di confrontare diversi metodi di stima della correlazione interno/esterno nel trattamento di forme tumorali extra-craniche. La problematica ha origine dal fatto che, nel trattamento radioterapico di tumori della zona addominale e toracica, l’organ motion, ovvero il movimento cui gli organi sono soggetti in conseguenza a processi fisiologici, inficia in modo sostanziale la riuscita del trattamento. L’organ motion può essere di 2 tipi: 1. Organ motion inter-frazione, ovvero movimenti che sono apprezzabili solo tra una seduta di trattamento e l’altra (cioè in un tempo sull’ordine dei giorni) e non all’interno della stessa seduta. Esso può essere causato da fenomeni come la perdita di peso del paziente, accrescimenti o riduzioni della massa tumorale, modificazioni dovute a processi fisiologici riguardanti organi quali vescica o intestino. 2. Organ motion intra-frazione, ovvero dovuto a processi i cui effetti sono apprezzabili anche all’interno dell’arco temporale di una singola seduta di trattamento, come la pulsazione del cuore e dell’aorta o il movimento respiratorio. Nel nostro studio abbiamo analizzato con maggior attenzione le problematiche derivanti dall’organ motion intra-frazione ed abbiamo esposto brevemente alcune tra le tecniche a disposizione per quantificare e compensare il movimento della massa tumorale, ovvero le tecniche basate sul Breath Hold, il Respiratory Gating o il Tumor Tracking. In particolare, il tumor tracking è una tecnica che prevede di tener traccia della posizione del tumore in real time e di riposizionare il fascio radiante dinamicamente in modo da seguire costantemente la lesione. Per far ciò, è necessario identificare la posizione del tumore in tempo reale e sviluppare modelli predittivi, in grado di anticipare il movimento del tumore, per tener conto delle latenze nel posizionamento del fascio. Un esempio di sistema per il tracking del respiro è rappresentato dal sistema SynchronyTM Respiratory Tracking System, un sottosistema del dispositivo robotico di trattamento Cyberknife (Accuray, Inc., Sunnyvale CA) per il trattamento di tumori extracranici che si muovono in conseguenza della respirazione. Il sistema RTS combina l’informazione proveniente da un sistema di sensori di posizione ad infrarossi, che monitora il movimento dell’addome del paziente, con quella fornita da coppie di immagini a raggi X, che forniscono invece un’informazione accurata sulla posizione interna del tumore, grazie a marker radiopachi impiantati nella lesione e visibili nelle immagini a raggi X. Gli emettitori ad infrarossi vengono posizionati sul petto e sull’addome del paziente ed il sistema di tracking ad infrarossi registra il movimento degli emettitori a frequenze molto elevate. Anche se il sistema di tracking fornisce un’informazione real time riguardo il movimento respiratorio del paziente, ciò potrebbe non essere sufficiente a determinare la posizione precisa della lesione, in quanto i marker esterni da soli non possono riflettere adeguatamente lo spostamento interno causato dal movimento respiratorio. Di conseguenza, vengono acquisite una serie di immagini a raggi X all’inizio del trattamento, che mostrano i marker interni ed esterni simultaneamente, e viene stabilita la correlazione esistente tra il movimento esterno e quello interno. Dato che il movimento respiratorio può variare anche considerevolmente nel corso della seduta, vengono acquisite ulteriori immagini durante il trattamento e fornite al modello di correlazione affinché venga aggiornata la corrispondenza interno/esterno. Il sistema RTS costruisce la correlazione utilizzando un modello lineare, se vi è una corrispondenza semplice tra i marker interni ed esterni, o un modello quadratico altrimenti. Una volta analizzato il funzionamento del sistema RTS, abbiamo preso in considerazione diversi metodi utilizzati in bibliografia per calcolare la correlazione interno/esterno, ovvero, filtri lineari e quadratici, filtri di Kalman, modelli ARMA, modelli di stato e reti neurali e, tra questi, abbiamo deciso di sviluppare gli ultimi due. Il modello di stato inferisce la posizione del tumore a partire non solo dall’informazione corrente della posizione dei surrogati esterni, ma anche da osservazioni ritardate di un tempo τ degli stessi. In questo modo è possibile cogliere la dinamica di primo ordine del segnale respiratorio. Una volta analizzate le equazioni che governano il modello, abbiamo implementato diversi modelli ed analizzato le loro performance al variare di alcuni parametri, ovvero, l’ingresso fornito al modello, il numero di campioni da utilizzare per l’addestramento, la frequenza di controllo dell’uscita, la soglia dell’errore oltre la quale considerare la corrispondenza non più valida ed il numero di campioni utilizzati per il conseguente riaddestramento, il ritardo τ. In base alle valutazioni effettuate abbiamo deciso di implementare i seguenti 3 modelli (dove I indica l’ingresso e O l’uscita): 1. modello quadratico con I: prima componente principale di movimento di ogni marker e O: tutte e 3 le coordinate del modeler; 2. modello quadratico con controllo lineare con I: prima componente principale calcolata sulle 3 coordinate del marker 1, O: tutte e 3 le coordinate del modeler; 3. modello quadratico con I: prima componente principale calcolata sulle 3 coordinate del marker 1, O: tutte e 3 le coordinate del modeler; Tutti i modelli implementati sono di tipo quadratico, ad eccezione del secondo modello in cui abbiamo provato ad aggiungere anche un modello lineare in modo che venga utilizzato l’uno o l’altro in base al valore assunto dal segnale esterno. Gli altri parametri sono stati mantenuti costanti per tutti e 3 i modelli e sono i seguenti:  Numero di campioni per l’addestramento iniziale: 500;  Numero di campioni per il riaddestramento durante la seduta: 500;  Soglia dell’errore oltre la quale riaddestrare: 3 mm; l’errore è stato valutato in termini di errore quadratico medio;  Frequenza di controllo dell’errore commesso: ogni 500 campioni (quindi ogni 20 s circa). Per quanto riguarda le reti neurali, il procedimento è stato simile: abbiamo testato diversi tipi di reti, variando l’ingresso, il numero di campioni da utilizzare per l’addestramento, la frequenza di controllo dell’uscita della rete, la soglia dell’errore oltre la quale considerare la corrispondenza non più valida ed il numero di campioni utilizzati per il conseguente riaddestramento, l’algoritmo di addestramento, la funzione di trasferimento, il numero di strati nascosti ed il numero di neuroni per strato. Le reti che hanno fornito i migliori risultati e che sono state quindi scelte per il nostro studio sono le seguenti: 1. rete a 3 strati nascosti più uno d’uscita, 2 neuroni per ogni strato nascosto (1 per quello d’uscita), funzione di trasferimento sigmoidale negli strati nascosti e lineare tra l’ultimo strato e l’uscita e algoritmo di apprendimento quasi-Newtoniano. L’ingresso è costituito da tutte le coordinate dei marker (quindi la matrice d’ingresso risulta avere dimensioni 9xN, dove N rappresenta la durata della seduta espressa in n° di campioni acquisiti alla frequenza di 26 HZ) e l’uscita dalle 3 coordinate della lesione (dimensioni 3 x N). Sono stati utilizzati 1000 campioni per l’addestramento iniziale, 200 campioni per il riaddestramento della rete esistente (nel caso l’errore commesso sia compreso tra 3 e 5 mm) e 1000 campioni per l’addestramento della nuova rete (nel caso l’errore sia superiore ai 5 mm); 2. Come la precedente, ma in questo caso sono stati utilizzati 500 campioni anziché 1000 per l’addestramento iniziale e 500 campioni per l’addestramento della nuova rete (nel caso l’errore sia superiore ai 5 mm); 3. In questo caso sono state implementate 3 diverse reti, ognuna a 3 strati nascosti più uno d’uscita, 2 neuroni per ogni strato nascosto (1 per quello d’uscita), funzione di trasferimento sigmoidale negli strati nascosti e lineare tra l’ultimo strato e l’uscita e algoritmo di apprendimento quasi-Newtoniano. L’ingresso di ogni rete è costituito dalle coordinate di 1 marker (quindi la matrice d’ingresso risulta avere dimensioni 3xN) e l’uscita dalle 3 coordinate della lesione (dimensioni 3 x N); l’uscita finale è quindi costituita dalla media delle 3 stime effettuate separatamente dalle 3 reti. Sono stati utilizzati 1000 campioni per l’addestramento iniziale, 200 campioni per il riaddestramento della rete esistente (nel caso l’errore commesso sia compreso tra 3 e 5 mm) e 1000 campioni per l’addestramento della nuova rete (nel caso l’errore sia superiore ai 5 mm); 4. Uguale alla precedente, ma sono stati utilizzati 500 campioni anziché 1000 per l’addestramento iniziale e 500 campioni per l’addestramento della nuova rete (nel caso l’errore sia superiore ai 5 mm). I 3 tipi di modelli e le 4 reti sono quindi stati utilizzati per stimare la corrispondenza in 20 frazioni provenienti da 16 pazienti del Database del Georgetown University Medical Center (Washington, DC). Tra questi vi sono 10 casi indicati come Control (C) caratterizzati da movimenti respiratori con una bassa deviazione standard e 10 casi indicati come Worst (W) in cui la deviazione standard è molto maggiore. I risultati si possono riassumere nei seguenti punti (in cui i modelli e le reti sono indicati secondo i numeri dell’indice sopra riportato):  I risultati ottenuti, in termini di errore quadratico medio, sono peggiori per i pazienti della classe W con tutti i metodi indagati;  Tra i modelli di stato indagati, il numero 1 è risultato il più accurato nella stima dell’uscita in tutti i casi di controllo, probabilmente grazie al fatto che riceve in ingresso le prime componenti principali di movimento di ogni marker; esso commette però errori molto consistenti in quasi tutti i pazienti della classe W. Per questi è invece il modello 2 a fornire la stima migliore, ovvero l’unico che utilizza sia un modello quadratico che uno lineare, in base al valore del segnale esterno.  Tra le reti neurali, la numero 2 ha effettuato la stima migliore nei pazienti di controllo in quasi tutti i casi, mentre nei pazienti della classe W non è stata riscontrata una prevalenza netta di una delle 4 reti indagate. Ciò che si è osservato è che le due reti a cui sono stati forniti 500 campioni per l’addestramento iniziale ed il riaddestramento durante la seduta hanno commesso errori inferiori rispetto alle reti corrispondenti a cui ne sono stati invece forniti 1000, utilizzando una porzione totale di segnale minore per addestrare.  Confrontando tra loro i modelli e le reti che hanno fornito i migliori risultati, si è osservato che queste ultime hanno riportato performance migliori rispetto ai modelli, soprattutto nei pazienti della classe W, per i quali, in alcuni casi, esse hanno commesso un errore quadratico medio inferiore di un ordine di grandezza. Analizzando i risultati ottenuti in rapporto alla deviazione standard del segnale respiratorio si può concludere che:  se la deviazione standard è inferiore a 3 mm circa l’utilizzo della rete 2 o del modello 1 risulta equivalente;  per una deviazione standard compresa tra 3 e 5 mm circa le reti 2 e 4 ed il modello 2 hanno fornito i risultati migliori;  se la deviazione standard supera i 5 mm i risultati migliori sono stati ottenuti con l’utilizzo delle reti 2 e 4 in quasi tutti i casi. Eventuali sviluppi futuri del lavoro potrebbero riguardare:  l’implementazione del controllo lineare anche per il modello di stato numero 1;  il calcolo della predizione con diversi modelli con valutazione in linea degli errori, implementando opportuni criteri di scelta basati sull’errore;  selezionare i parametri che definiscono l’architettura della rete utilizzando criteri di ottimizzazione delle condizioni iniziali.