Predicting the evolution of cervical cancer during the treatment using ODE models and validation by mean of CBCT images

In the fight against cancer, radiotherapy represents a highly effective localized treatment, in which image-guidance (Image Guided RadioTherapy, IGRT) plays a fundamental role. Images featuring diagnostic quality are acquired daily, in order to check the occurrence of anatomical changes, to adjust patient setup and to evaluate the need for treatment re-planning. In addition, functional imaging to visualize tumor characteristics (e.g., glucose consumption, hypoxia, repopulation) is applied, not only for the initial tumor staging, but also at given time intervals during the fractionated course of treatment. Despite advances in recognizing both the physics and radiation biological mechanisms underpinning radiotherapy treatment, early prediction of tumor response and toxicity effects remains difficult to achieve, particularly at the level of the individual patient. Over the last 10 years, attempts have been made to apply models to predict tumor response and toxicity occurrence on a patient specific level However the translation of computational approaches into effective clinical practice has been hindered by lack of convincing validation. As far as tumor response and therapy outcome prediction are concerned, many efforts have been tackled to reinforce the clinical potentials of in silico models towards their translation into radiotherapy clinical practice, for the optimization and personalization of the treatment. Biological modeling of the tumor response however is still playing only a very limited role in commercial treatment planning system and is restricted to the use of classical tumor control probability models. These models make their predictions using population-based input parameters. Due to the dynamic character of tumor properties, this is not sufficient to predict the development of tumors in the individual case. Tumor heterogeneity caused by the different radio-sensitivities of cell sub-populations, as well as by varying environmental factors such as oxygenation, prevents researchers to get reliable predictions for the individual patient. This situation can only be improved if individual input parameters, as well as mechanistic models on the spatio-temporal development of the tumor, are systematically considered. To achieve this, specific in-silico models have been developed, which simulate the spatio-temporal growth of tumors, based on basic radiobiological principles. Single-cell-based models using the principles of cellular automata have been developed and applied to study a variety of radio-biological problems. Although these models may well be used to simulate the mechanisms underlying radiation response, they are limited with respect to clinical application as the morphological and functional information from medical images are difficult to include. This led to the development of voxel-based tumor-models, which preserve some of the properties of single-cell-based models, but which have the advantage that they can process information from medical images. Combining mechanisms on a cellular level and treating them with statistical sampling method within a voxel-based geometrical patient model provides the first step towards multi-scale models. Although some approaches in developing such models have already been made in the past, they have not been translated into clinical application yet. In this thesis, we aimed at investigating and retrospectively validating computer based models for tumor response to radiation therapy applied to the cervical cancer, a malignant neoplasm located in the cervix. Imaging data of 9 patients were provided by IEO (European Institute of Oncology, Milan, Italy). For each patient, the dataset was constituted by one planning CT study (Computed Tomography) along with about 10 CBCTs (Cone Beam CT) studies, acquired in correspondence of the irradiation treatment fractions. The CBCTs served as control point for training and prediction evaluation. In the patient dataset, two different groups, who differ for cancer type and therapies, were identified. First group included 5 patients affected by squamous cell carcinoma who were administered with radiation only. Second group included 4 patients affected by adenocarcinoma (AC) who were administered with radiation and chemotherapy. Two different modeling approaches were studied: 1) cell-based discrete model, implemented through Cellular Automata; 2) tissue-level continuous model, implemented through ODE. Due to high computational burden, the cell-based discrete model underwent just preliminary simulated tests and was disregarded from extensive validation focusing the main work on the tissue-level continuous model. The tissue-level continuous model used Gompertzian or Logistic functions for the tumor growth while three different formulations were considered to represent the relation with the dose: a) classical linear/quadratic model; b) exponential; c) linear. The model training was implemented using an Evolutionary computation approach, it was repeated for different formulation of the tumor growth and therapy response terms, and performed across the two different patient groups. The tumor shape was simplified as a spheroid evolving throughout time. Volume and area of the tumor mass, measured in the CBCT datasets (expert manual segmentation and automatic 3D surface reconstruction) and predicted by the model, were chosen as training cost functions and validation quantities. Tumor volume regression rate has been suggested in literature as an important predictor of local control and long-term survival of the patient. The residual errors of predicted volume and area for the first groups of patients were (absolute RMS) 5.3 cm3 (relative error: 27%) and 10.4 cm2 (relative error: 24%), respectively. Interestingly, the Gompertzian model combined with the exponential one provided best predictive results in both tumor volume and area. Smaller values were obtained in the second group, namely, 4.6 cm3 (relative error: 21%) for the volume and 7.1 cm2 (relative error: 16%) for the area. A sensitivity analysis was performed by studying the stability of the tumor mass volume-area prediction to the perturbation of specific model parameters. The natural tumor growth rate, represented by the first coefficient of each model, resulted generally the most critical, but in some cases after a while the weight of the therapy effects resulted to be the major cause of the prediction variation. In order to test extrapolation properties of the models, leave-one-out validation was considered. It consisted in removing one patient at a time from its group, and in training the model across the remaining patient datasets. As expected, the prediction error on the one-out patient was generally larger than on the rest of the group. Interestingly, for some patients, the prediction ability was comparable or even better than that of the remaining patient in the group. For example in the first group, considering the area results, the mean was about 25% while it reached 40% on the test, and the variance among the results was less than 0.05 on the first while above 0.25 on the second one. As expected, the quality of the prediction on the untrained patient was much poorer since the model fitted a small amount of data. As a first consideration, the obtained prediction results were promising as they were in the range of the data uncertainty, qualitatively about 10% and 20% for area and volume, respectively. Noteworthy, such results were basically coherent with recent domain literature (Huang Z. et al. 2010) where relative prediction errors of about 25% were reported. As far as model training and validation is concerned, results highlighted that volume and area prediction criteria conflicted in between for the second group. In detail, while the combination of the Logistic curve and the linear/quadratic model turned out to be the most performing on the area, the volume was better fitted using the Gomperzian model along with the linear model of the radiation response. The different behavior of the two indexes can be motivated by considering the effect of the chemotherapy on the overall treatment. The drugs are indeed delivered by mean of blood vessels that are more densely localized on the outer surface of the tumor and this peculiarity could cause an increase in the relevance of the events occurring to the most superficial layer of cells. While present results are not decisive, this finding reinforces the concept that the representation ability of the model is fundamental as well as the cost function to minimize for model parameter assessment (training). This could play a significant role in elaborating more complex predicting approaches. A number of issues in this work must be highlighted. First, we represented tumor growth and therapy response, using tissue-level models implemented through continuous equations (ODE), thus limiting the overall process representation to just macro-scale mechanisms. Micro-scale processes at cellular level as for example cell cycle, cell-to-cell interaction, and tumor cell invasion weremodeled using cellular automata but the obtained models were not validated. Second, in the mathematical representation of tumor evolution we neglected the spatial dimension and only scalar validation criteria (area and volume) were adopted, approximating the tumor shape as a spheroid. This allowed to avoid registration across CBCT studies but loosing tumor progression in a common coordinate system. Third, the patient cohort was limited in size and this prevented us from getting any statistical evidence about the value of the prediction results and model performance. Fourth, the uncertainty on the data used for validation was just qualitatively evaluated. This aspect is however critical for any consideration about the model training feasibility and performance. Fifth, in the model training, the nominal isodose, administered to the patient, was used. However, the actual dose can be different due to tumor shape modification, soft tissue deformation, physical interaction with the surrounding organ and gravity effects. In order to overcome this drawback, actual dose computation should be carried out exploiting deformable registration between the planning CT volume, in which the nominal does is calculated, and the actual CBCT volume taken at treatment time. Nonetheless, despite such shortcomings, the presented work showed how also not complex tissue scale models can be effective in predicting the interaction between tumor progression and the effects of radiation therapy. In conclusion, while being a preliminary study, the results are meaningful and promising for extending the methodology to include more sophisticated representation of the biological processes and therapy interaction, and making the approach suited for more extensive validation and future clinical translation.

Nella lotta contro il cancro, la radioterapia rappresenta un trattamento localizzato altamente efficace, nel quale la guida per immagini (Image Guided RadioTherapy, IGRT) svolge un ruolo fondamentale. Le immagini con qualità diagnostica sono acquisite giornalmente, allo scopo di appurare la presenza di cambiamenti anatomici e valutare la necessità di una riorganizzazione del trattamento del paziente. Inoltre, l’imaging funzionale per la visualizzazione delle caratteristiche tumorali (ad esempio consumo di glucosio, ipossia, ripopolamento) è applicato, non solo per l’iniziale staging tumorale, ma anche a intervalli di tempo definiti durante il trattamento frazionato. Nonostante l’avanzamento nel riconoscimento dei meccanismi fisici e biologici legati alla radioterapia, la predizione tempestiva della risposta del tumore e degli effetti tossici rimane difficoltosa da raggiungere, in particolare a livello del singolo paziente. Negli ultimi 10 anni, si è tentato di applicare modelli alla predizione di questi eventi sui singoli casi di trattamento, tuttavia l’introduzione di approcci computazionali nella pratica clinica è stata impedita dalla mancanza di validazioni convincenti. Per quanto riguarda la predizione dello sviluppo tumorale e degli effetti delle terapie, molti sforzi sono stati compiuti per rinforzare le potenzialità cliniche dei modelli in silico attraverso la loro traslazione nella pratica clinica per l’ottimizzazione e personalizzazione del trattamento. La modellizzazione biologica della risposta tumorale tuttavia ricopre un ruolo molto ridotto nei sistemi commerciali di pianificazione del trattamento, che si limita ai classici modelli di probabilità del controllo tumorale. Questi modelli compiono le loro predizioni attraverso parametri d’ingresso basati su popolazioni. A causa delle dinamiche delle caratteristiche del tumore, questo non è sufficiente a predire lo sviluppo dei tumori sui casi individuali. L’eterogeneità dei tumori, causata dalla differente radio-sensitività delle sottopopolazioni cellulari e dalla variazione di parametri ambientali come l’ossigenazione, impedisce predizioni attendibili paziente specifiche. Questa limitazione può essere superata unicamente considerando sistematicamente sia modelli meccanici dello sviluppo spaziotemporale del tumore che parametri tarati sul singolo individuo. Per raggiungere questo , sono stati sviluppati specifici modelli computazionale, in grado di simulare la crescita spazio temporale della massa tumorale, basati su principi radiobiologici. Modelli di interazioni tra le singole cellule sono stati sviluppati utilizzando i principi dei Cellular Automata e sono stati applicati allo studio di diversi problemi radio-biologici. Sebbene questi modelli possano essere usati per simulare i meccanismi che regolano la risposta alle radiazioni, sono limitati nell’utilizzo clinico dalla difficoltà di includere informazioni morfologiche e funzionali sufficientemente accurate tratte da immagini mediche. Questo porta allo sviluppo di modelli basati su voxel, che conservano alcune delle proprietà dei modelli a singola cellula, ma nei quali si ha il vantaggio di poter utilizzare le informazioni derivanti dall’imaging clinico. Combinando meccanismi a livello cellulare con metodi statistici all’interno di un voxel dalle dimensioni ridotte, si compie un passo verso l’approccio multi-scala. Nonostante alcuni tentativi in questa direzione siano stati approcciati in passato, non sono ancora stati tradotti in pratica clinica. In questa tesi intendiamo investigare e validare retrospettivamente modelli computazionali sulla risposta tumorale alla radioterapia applicati al cancro della cervice uterina . Dallo IEO (European Institute of Oncology, Milano, Italia)sono stati forniti dati di imaging relativi a 9 pazienti . Per ciascun pazente , il dataset era costituito da uno studio CT (Computed Tomography) insieme a circa 10 studi CBCT (Cone Beam CT), utilizzati come punti di controllo per il training e la valutazione della predizione. Nel dataset dei pazienti sono stati identificati due distinti gruppi. Il primo gruppo includeva 5 pazienti affetti da carcinoma squamocellulare (SC) i quali sono stati trattati unicamente per mezzo di radiazioni. Del secondo gruppo facevano parte 4 pazienti colpiti da adenocarcinoma(AC) che hanno subito sia trattamento radioterapico che chemioterapico. Sono stati analizzati due differenti approcci modellistici: 1) modello discreto a livello cellulare implementato attraverso la tecnica ei Cellular Automata; 2) modello continuo a livello tissutale che prevede l’utilizzo di ODE. A causa dell’alto costo computazionale, il modello basato sulle interazioni tra cellule è stato spinto solo fino a test e simulazioni preliminari, mentre la validazione estesa è stata focalizzata unicamente sul secondo tipo di approccio. Il modello su scala tissutale utilizzava funzioni Gomperziane o Logistiche per la crescita tumorale mentre sono state considerate tre possibili formulazioni per rappresentare la relazione con la dose: a) classico modello lineare/quadratico; b)esponenziale; c)lineare. L’addestramento del modello è stato implementato utilizzando tecniche evolutive ed è stato ripetuto per le diverse combinazioni delle formulazioni del termine di crescita tumorale e di risposta alla terapia per entrambi i gruppi. La forma del tumore è stata semplificata considerandolo uno sferoide in evolvuzione nel tempo. Il volume e l’area della massa tumorale, ricavata dal dataset delle immagini CBCT (attraverso segmentazione manuale di un esperto e ricostruzione automatizzata della superficie 3D) e predetta dal modello, sono state utilizzate nella funzione di costo per l’ottimizzazione del modello e per la sua validazione. Il tasso di regressione del volume tumorale è stato presentato in letteratura come un importante indicatore di riuscita del controllo locale e predittore della sopravvivenza a lungo termine del paziente. Gli errori residui relativi a volume ed area predetti per il primo gruppo sono (in valori assoluti RMS) 5.3 cm3 (errore relativo: 27%) and 10.4 cm2 (errore relativo: 24%), rispettivamente. È interessante notare come il modello Gomperziano combinato con quello esponenziale abbia fornito i risultati migliori sia per il volume che per l’area. Valori più bassi sono stati ottenuti per il secondo gruppo, ossia, 4.6 cm3 (errore relativo: 21%) per il volume e 7.1 cm2 (errore relativo: 16%) per l’area. Un’analisi di sensitività è stata svolta studiando la stabilità della predizione di area e volume in risposta alla perturbazione di specifici parametri. Il tasso di crescita naturale del tumore, rappresentato dal primo coefficiente di ciascun modello, è risultato generalmente il più critico, ma in alcuni casi durante il corso del trattamento il peso degli effetti della terapia è risultato essere la maggiore causa di variazione nella predizione. Allo scopo di testare la capacità di estrapolazione del modello è stata utilizzata la tecnica di validazione leave-one-out. Questa consiste nel rimuovere on paziente alla volta dal suo gruppo ed addestrare il modello sui rimanenti. Come ci si aspettava, l’errore di predizione commesso sul paziente lasciato fuori è risultato generalmente maggiore di quello ottenuto mediamente sul resto del gruppo. Interessante notare come per alcuni pazienti rimossi la capacità predittiva sia risultata comparabile, o addirittura migliore, a quella del resto del gruppo. Per esempio nel primo gruppo, considerando i risultati ottenuti per l’area, la media dell’errore è risultata del 25% mentre raggiungeva il 40% per i pazienti di test, e la varianza, molto contenuta tra i primi (0.05) era molto più elevata nei secondi (0.25) Come ci si aspettava, la qualità della predizione sui pazienti non utilizzati per l’addestramento è risultata ridotta dal momento che il modello aveva a disposizione solo pochi esempi, insufficienti per una corretta generalizzazione. Come prima considerazione, i risultati di predizione ottenuti sono stati promettenti poiché non lontani dal range di incertezza dei dati , qualitativamente intorno al 10% e 20% rispettivamente per area e volume. Per quanto riguarda i modelli, i risultati hanno evidenziato una discrepanza tra il criterio di predizione di area e volume nel secondo gruppo. In dettaglio, mentre la combinazione del modello lineare/quadratico e della funzione Logistica si è rivelata la più performante per quanto riguarda l’area, nel caso del volume il fitting migliore è stato ottenuto con la curva Gomperziana combinata con la risposta al trattamento lineare con la dose. Il diverso comportamento può essere spiegato considerando l’effetto ella chemioterapia sul trattamento complessivo. I farmaci sono trasportati per mezzo dei vasi sanguigni maggiormente presenti sulla superficie esterna della neoplasia, questa caratteristica potrebbe portare ad una particolare rilevanza dei fenomeni che si presentano nello strato cellulare più esterno. Sebbene questi risultati siano preliminari, è possibile supporre che oltre alla tipologia di modello anche il metodo di ottimizzazione e validazione giochino un ruolo importante nella realizzazione di approcci predittivi più complessi. Vanno evidenziate alcune criticità di questo lavoro. Per prima cosa, abbiamo rappresentato la crescita tumorale e la risposta alla terapia usando modelli continui a livello tissutale, implementati per mezzo di equazioni differenziali (ODE), limitando la rappresentazione ai meccanismi presenti su scala macroscopica. I processi su micro sala, come il ciclo cellulare e le interazioni tra cellule sono stati studiati attraverso l’approccio dei Cellular Automata ma il modello ottenuto non è stato validato. In secondo luogo, nella rappresentazione matematica dell’evoluzione tumorale abbiamo trascurato la spazialità, ma abbiamo fatto riferimento unicamente a grandezze scalari (area e volume) approssimando la massa tumorale ad uno sferoide. Questo ci ha permesso di evitare la registrazione dei vari studi CBCT tra loro, ma ci ha fatto perdere la progressione morfologica del tumore in un sistema di riferimento cartesiano univoco. Terzo, la dimensione limitata del dataset ci ha impedito di ottenere evidenze statistiche significative riguardo il valore dei risultati di predizione e performance del modello. Quarto, l’incertezza sui dati utilizzati per la validazione è stata valutata solo qualitativamente. Questo aspetto è tuttavia fondamentale per qualunque considerazione riguardo i risultati della predizione, una valutazione oggettiva del rumore sui dati sarebbe indispensabile alla valutazione della qualità dei risultati ottenuti. Quinto, nell’addestramento del modello è stata utilizzata la dose nominale somministrata al paziente, tuttavia la dose reale può differire sensibilmente dalla precedente a causa di modificazioni morfologiche del tumore, deformazioni dei tessuti molli e fisiologici, spostamenti degli organi circostanti e altri fenomeni. Per superare questo problema, bisognerebbe sfruttare una registrazione con deformazione per calcolare la dose effettivamente rilasciata nella regione di interesse e comprendere quale sia il volume realmente trattato. Tuttavia, nonostante questi inconvenienti, il presente lavoro mostra come, anche modelli non complessi e su scala tissutale, possano essere efficaci nel predire l’interazione tra la progressione tumorale e gli effetti della radioterapia. I principali risultati della predizione sono stati essenzialmente coerenti con quanto rintracciabile nella recente letteratura (Huang Z. et al., 2010) dove sono stati riportati errori relativi di circa 25%. In conclusione, pur essendo uno studio preliminare, i risultati ottenuti sono significativi e promettenti per estendere la metodologia, includere rappresentazioni più sofisticate dei processi biologici e delle interazioni dovute al trattamento, e rendere il modello adatto per una validazione più estesa con la conseguente traduzione in pratica clinica.