Mathematical models of tumor response in advanced radiotherapy techniques

INTERNATIONAL DOCTORAL SCHOOL OF THE USC Isabel González Crespo PhD Thesis Mathematical models of tumor response in advanced radiotherapy techniques Santiago de Compostela, 2024 Doctoral Programme in Mathematical Modelling and Numerical Simulation in Engineering and Applied Science Ph.D. Thesis MATHEMATICAL MODELS OF TUMOR RESPONSE IN ADVANCED RADIOTHERAPY TECHNIQUES Isabel Gonz´alez Crespo Director: Juan Pardo Montero Director-Tutor: ´ Oscar L´opez Pouso ESCOLA DE DOUTORAMENTO INTERNACIONAL DA UNIVERSIDADE DE SANTIAGO DE COMPOSTELA PROGRAMA DE DOUTORAMENTO EN M´ ETODOS MATEM´ ATICOS E SIMULACI´ ON NUM´ ERICA EN ENXE ˜ NAR´ IA E CIENCIAS APLICADAS SANTIAGO DE COMPOSTELA 2024 Dedications A meus pais e meu irm´an, por crer sempre en min. Acknowledgements En primeiro lugar, quero agradecerlle aos meus directores de tese, Juan e ´ Oscar, a s´ua gu´ıa e consello, o tempo dedicado e a confianza depositada desde que comecei o meu Traballo de Fin de M´aster. A Juan, grazas tam´en por acompa˜narme aos congresos, polos caf´es de “media ma˜n´a” ´a unha da tarde, e pola paciencia no d´ıa a d´ıa. M´agoa que non conseguiras que apunte sempre as cousas na libreta. E a ´ Oscar, grazas por ser un gran docente e por achegarme de novo ´a beleza das formalidades da Matem´atica. Grazas a meus pais, Manolo e Marga, por permitirme chegar ata aqu´ı e acompa˜narme en cada paso; todos os meus logros son tam´en vosos. Ao meu irm´an, Nico, grazas por estar sempre ao meu lado; espero que te sintas tan orgulloso de min como me sinto eu de ti. Grazas aos meus av´os, Faustino, Herminia, Isabel e Manolo, que son un referente de resiliencia e esforzo. En especial ´a abuela Herminia, que se sentou comigo a facer os deberes tantas tardes despois do colexio, grazas por axudarme en todo. ´ A mi˜na familia ao completo, grazas por apoiarme e desexarme sempre o mellor. A Ruben, grazas pola comprensi´on, o apoio e a paciencia inesgotable en cada etapa da tese. O cami˜no ´e m´ais f´acil contigo. Grazas ´as mi˜nas amigas, Ana, Isa e Sara, polas risas e os momentos de desconexi´on; a Susi e Tere, por animarme a seguir desde antes de ser “tres matem´aticas”; e a Ux´ıa por escoitarme, estar presente, e mandarme trevos da sorte. ´ As mi˜nas compa˜neiras e compa˜neiros do Servizo de Radiof´ısica e Protecci´on Radiol´oxica do CHUS, e ao Grupo de F´ısica M´edica e Biomatem´aticas do IDIS, grazas por acollerme durante estes anos. En especial a Araceli, Arnau e Julio, polas conversas e os caf´es. Por ´ultimo, grazas a todas as persoas que formaron parte desta aventura. November 11, 2024 Additional acknowledgements The development of this doctoral thesis has been partially funded by FIDIS (Fundaci´on do Instituto de Investigaci´on Sanitaria de Santiago) through the “Predoutorais FIDIS 2021–2024” call, for which I was awarded a 3-year contract (January 2021–January 2024) within the Medical Physics and Biomathematics Group of IDIS (Instituto de Investigaci´on Sanitaria de Santiago). This work has received financial support from Xunta de Galicia within the project “Optimization and individualization of treatments in advanced radiotherapy techniques” (2022–2026, IN607D2022/02); and the Ministerio de Ciencia e Innovaci´on within the projects “Optimization of the prescription dose and dose/volume limits in the treatment of prostate cancer: an approach from big data” (2023–2025, PID2021), and “Advanced dosimetry for novel radiotherapy approaches in brain tumors (DOSEBRAIN)” (2022–2025, PLEC2022-009476). Isabel Gonz´ alez Crespo Contributions: conceptualization, data curation, formal analysis, software implementation, validation, and writing. Affiliations: 1Group of Medical Physics and Biomathematics (IDIS) 2Department of Applied Mathematics (USC). 3Department of Particle Physics (USC). 4Department of Medical Physics (CHUS). The contents of this article are partially reproduced in Chapter 4. *2024 data was not available at the date of finalization of this document. 2 Resumo O termo cancro comprende m´ais de duascentas enfermidades xen´eticas distintas, caracterizadas pola aparici´on de c´elulas anormais que se multiplican de maneira incontrolada, levando na maior´ıa dos casos ´a formaci´on de tumores s´olidos. Ademais, poden diseminarse nos tecidos circundantes e chegar a zonas distantes do organismo derivando en tumores metast´aticos. Actualmente, a radioterapia ´e a t´ecnica m´ais empregada no tratamento do cancro (en arredor do 60 % dos pacientes segundo a OMS), ben sexa como monoterapia ou en combinaci´on con outras t´ecnicas, tales como a quimioterapia e a cirurx´ıa. Este tratamento consiste en depositar altas doses de radiaci´on sobre o tumor para eliminalo ou minguar a progresi´on da enfermidade. A radiaci´on causa modificaci´ons no ADN celular derivando en danos letais, que causan a morte directa da c´elula, ou subletais, que no caso de non ser reparados poden acumularse coa aplicaci´on de fracci´ons de dose sucesivas levando ´a morte da c´elula. Por outra banda, a radioterapia tam´en afecta a outros axentes que forman parte do tecido tumoral, como o sistema vascular ou as c´elulas do sistema inmune, derivando en mecanismos de morte indirecta. Desde o desenvolvemento dos primeiros tratamentos de radioterapia a principios do s´eculo XX, os numerosos estudos realizados e os sucesivos avances tecnol´oxicos permitiron desenvolver tratamentos m´ais efectivos ao mesmo tempo que reducir os efectos secundarios prexudiciais para os pacientes. Observouse que as c´elulas tumorais son m´ais sensibles ao efecto da radiaci´on que as c´elulas normais, xa que a s´ua taxa de multiplicaci´on acelerada fainas m´ais susceptibles de recibir dano no seu ADN. Ademais, as c´elulas tumorais son menos eficaces ´a hora de reparar o dano provocado pola radiaci´on debido ´as mutaci´ons que sofren. Isto levou ao dese˜no dos tratamentos de radioterapia convencionais, nos cales a dose total pautada (medida en Gy) div´ıdese en fracci´ons de dose administradas en d´ıas consecutivos (habitualmente 30–40 fracci´ons cunha fracci´on diaria de 1,8–2 Gy de luns a venres) ata a finalizar o tratamento. Desta maneira as c´elulas normais poden reparar parte do dano non letal ocasionado entre fracci´ons de dose, mentras que a s´ua acumulaci´on contrib´ue ´a eliminaci´on das c´elulas tumorais. Pouco despois da aparici´on da radioterapia comezaron a desenvolverse modelos matem´aticos co obxectivo de predicir o resultado dos tratamentos, investigar o seu Isabel Gonz´ alez Crespo efecto sobre o tecido tumoral e o tecido san, e asistir na planificaci´on da terapia no ´ambito cl´ınico. O modelo lineal-cadr´atico (Linear-Quadratic, LQ) formulouse m´ais de medio s´eculo atr´as de maneira paralela en varios traballos. Tr´atase dunha ecuaci´on sinxela que consta dun termo lineal e outro cadr´atico, e relaciona a dose de radiaci´on recibida coa fracci´on de supervivencia celular resultante da s´ua administraci´on. Esta f´ormula obt´ıvose de maneira fenomenol´oxica mediante o axuste de datos experimentais in vitro e serve para estimar o efecto directo da radiaci´on sobre as c´elulas. O modelo involucra dous par´ametros, un para cada termo da ecuaci´on, que toman valores espec´ıficos para cada tumor, e que se relacionan coa sensibilidade das c´elulas ao dano letal e subletal, sendo respectivamente α(lineal) e β(cadr´atico). A fracci´on de supervivencia empr´egase para obter outras m´etricas con aplicaci´on cl´ınica, como a probabilidade de control tumoral (Tumor Control Probability, TCP), que permite estimar a probabilidade de que un tratamento de radioterapia elimine completamente o tumor; e a probabilidade de complicaci´on no tecido normal (Normal Tissue Complication Probability, NTCP), que serve para calcular a probabilidade de que a dose administrada aporte toxicidade nos tecidos sans. Un factor que infl´ue na supervivencia celular ´a radioterapia ´e o estado de osixenaci´on das c´elulas. ´ E sabido que os baixos niveis de osixenaci´on (hipoxia) est´an asociados a unha maior resistencia ao dano por radiaci´on. Isto d´ebese tanto a factores biol´oxicos como fisicoqu´ımicos. Por un lado, as c´elulas hip´oxicas empregan certas prote´ınas chamadas factores inducibles por hipoxia (Hypoxia Inducible Factors, HIFs) para promover a s´ua proliferaci´on. As c´elulas tumorais benef´ıcianse deste mecanismo para potenciar a s´ua progresi´on. Por outro lado, a presenza de os´ıxeno ten un papel clave na eliminaci´on das c´elulas tumorais, xa que promove a formaci´on de especies reactivas que danan o ADN celular. Polo tanto, a hipoxia ´e aceptada como un factor que compromete a efectividade dos tratamentos de radioterapia. Para modelar o papel do os´ıxeno no resultado dos tratamentos, introduciuse o modelo LQ coas raz´ons de mellora polo os´ıxeno (Oxygen Enhancement Ratios, OERs), que tende ´a f´ormula cl´asica do modelo LQ en condici´ons totalmente ´oxicas. De acordo con estudos publicados, a introduci´on da radioterapia hipofraccionada (menor n´umero de fracci´ons con maiores doses por fracci´on, arredor de 8–30 Gy) supuxo unha mellora na calidade de vida dos pacientes, aportando maior control tumoral e menos efectos secundarios no tratamento de determinados tipos de cancro. Numerosos estudos levaron a debate a capacidade do modelo LQ para reproducir o efecto dos tratamentos hipofraccionados, apuntando a un fen´omeno de saturaci´on do dano para as altas doses por fracci´on empregadas, froito da reparaci´on do dano subletal. Xurdiron as´ı o modelo lineal-cadr´atico-lineal (Linear-Quadratic-Linear, LQL) e outras formas derivadas do modelo LQ, que incorporan par´ametros adicionais para modular o efecto das altas doses na fracci´on de supervivencia. Nas ´ultimas d´ecadas xurdiron novas t´ecnicas de radioterapia avanzada para as cales o modelo LQ e as s´uas variantes son insuficientes. Os novos tratamentos involucran mecanismos de morte celular indirecta, que non resultan unicamente do dano celular 4 Resumo causado pola radiaci´on. ´ E necesario incorporar no modelo o efecto da terapia sobre outros axentes presentes no microambiente tumoral, como as c´elulas do sistema inmune, e outros factores que alteran o resultado do tratamento, como os cambios na osixenaci´on durante a s´ua administraci´on. Nesta tese pres´entanse novos modelos de resposta tumoral a t´ecnicas avanzadas de radioterapia. Por un lado, para o uso combinado de radioterapia e inmunoterapia, co˜necido como radioinmunoterapia. Por outro lado, para a radioterapia FLASH, que consiste en empregar taxas de dose m´ais altas permitindo a administraci´on de maiores doses de radiaci´on en menos tempo de tratamento. Os modelos presentados constr´uense tomando como referencia os modelos cl´asicos mencionados anteriormente. A inmunoterapia ´e unha t´ecnica que consiste en potenciar o efecto antitumoral do propio sistema inmunitario do paciente interferindo nalgunha das fases do chamado ciclo inmune do cancro. Este ciclo consta de distintas fases mediante as cales a presenza de axentes nocivos no organismo, como as c´elulas tumorais, d´a lugar ´a liberaci´on de sinais espec´ıficos (ant´ıxenos) que activan as c´elulas T do sistema inmune, encargadas de reco˜necer e eliminar o axente contra o cal foron especializadas. O organismo conta con puntos de control inmunitario para evitar unha resposta inmune excesiva, que poder´ıa derivar en enfermidades inflamatorias. As c´elulas do cancro benef´ıcianse destes mecanismos para evitar a morte inmune manifestando inhibidores dos puntos de control. O receptor CTLA-4 ´e un mecanismo de control que interfire na fase de activaci´on das c´elulas T. Para a correcta especializaci´on e activaci´on das c´elulas T, o ligando CD28, presente na s´ua superficie, debe unirse co CD80 ou CD86 na superficie das c´elulas encargadas de transportar e presentar os ant´ıxenos. Pola contra, co fin de evitar unha sobreactivaci´on, o CD28 pode unirse co receptor CTLA-4 dando lugar a unha c´elula T inactiva. Para potenciar a activaci´on, desenvolveuse a inmunoterapia con anticorpos inhibidores do CTLA-4 (anti-CTLA4). Nesta terapia, o anticorpo admin´ıstrase de maneira intravenosa e ´unese co receptor CTLA-4 favorecendo a uni´on dos ligandos que levan ´a activaci´on das c´elulas T. O eixe PD-1/PD-L1 constit´ue outro punto de control e afecta ´a identificaci´on das c´elulas tumorais por parte das c´elulas T xa activadas. O PD-1 ´e unha prote´ına que est´a presente na superficie das c´elulas T, e o seu ligando, o PD-L1, manif´estase noutros axentes do microentorno tumoral. Cando estes dous ligandos se unen, a c´elula T “desact´ıvase”. As c´elulas tumorais son capaces de manifestar PD-L1 na s´ua superficie e escapar da acci´on das c´elulas T. As inmunoterapias inhibidoras do PD-1 e o PD-L1 (anti-PD(L)1) admin´ıstranse tam´en de maneira intravenosa co obxectivo de favorecer a actividade inmune. Os primeiros f´armacos de inmunoterapia anti-CTLA4 e anti-PD1 aprob´aronse para o seu uso cl´ınico nos anos 2011 e 2014, respectivamente. Malia resultar beneficiosos no tratamento do cancro, ofrecen unha mellora limitada en gran parte dos pacientes, xa que a s´ua efectividade depende do bo funcionamento do propio sistema inmune. Por este motivo, ad´oitase usar en combinaci´on con outras t´ecnicas como a radiotera5 Isabel Gonz´ alez Crespo pia, dando lugar ´a radioinmunoterapia. Ademais, o uso de radiaci´on d´a lugar tanto a efectos inmunosupresivos como inmunox´enicos. Por un lado, dana as c´elulas inmunes presentes no tumor, limitando a s´ua acci´on. Por outro lado, incrementa a liberaci´on de ant´ıxenos, promovendo unha maior activaci´on de c´elulas T. Numerosos estudos experimentais apuntan a que a radioinmunoterapia leva a maiores taxas de control tumoral en comparaci´on co uso das d´uas t´ecnicas de maneira independente. Na actualidade, esta t´ecnica combinada contin´ua sendo investigada de cara a dese˜nar tratamentos ´optimos onde se atope a combinaci´on de ambas terapias que aporte maior beneficio aos pacientes. O modelo de radioinmunoterapia proposto nesta tese ´e de tipo mecanicista e compartimental. Describe a din´amica dos principais axentes involucrados na resposta ´a radioinmunoterapia seguindo o ciclo inmune do cancro: c´elulas tumorais non danadas pola radiaci´on ou viables, c´elulas tumorais danadas que ser´an eliminadas, c´elulas T activas e ant´ıxenos (ou c´elulas presentadoras do ant´ıxeno). Ademais incorpora un submodelo de activaci´on no que participan os ant´ıxenos, as c´elulas inmunitarias dispo˜nibles para ser activadas, as c´elulas T activas e as c´elulas T bloqueadas polo efecto do PD-1/PD-L1. Este modelo est´a formado por ecuaci´ons diferenciais ordinarias, con retardos e con impulsos. A´ında que se trata dun modelo compartimental, tense en conta a compo˜nente espacial de maneira indirecta mediante a distinci´on de d´uas localizaci´ons do organismo que levan a inclu´ır retardos temporais na din´amica das variables consideradas. Por un lado, a liberaci´on de ant´ıxenos, a morte inmune, e o efecto da radioterapia e a inmunoterapia con anti-PD(L)1 ocorren no tumor. Por outro lado, a activaci´on das c´elulas T e o efecto da inmunoterapia con anti-CTLA4 d´ase nos ´organos linfoides. Polo tanto, as ecuaci´ons con retardo empr´eganse para ter en conta o desprazamento dos ant´ıxenos desde o tumor ata os ´organos linfoides, as´ı como das c´elulas T activadas facendo o percorrido inverso. En canto ao efecto dos tratamentos, introd´ucense no modelo mediante impulsos nos tempos de administraci´on. Por un lado, a radioterapia mod´elase mediante o c´alculo da fracci´on de supervivencia. Nos instantes de tratamento, a poboaci´on de c´elulas tumorais viables (que aumenta seguindo unha forma lox´ıstica para simular a s´ua multiplicaci´on) vese reducida a unha fracci´on da mesma, que pasa a formar parte do compartimento de c´elulas danadas para ser eliminada progresivamente cunha taxa de morte exponencial. Ao mesmo tempo, consid´erase que unha certa fracci´on das c´elulas T activas no tumor ´e eliminada de maneira instant´anea. No c´alculo da fracci´on de supervivencia investig´aronse tres modelos: o LQ, o LQL e unha variante do LQ. Ademais, estudouse o efecto inmunosupresor da radiaci´on que xorde ao danar a vasculatura tumoral limitando o acceso das c´elulas T ´a zona. Para isto, empregouse un factor modulador no compartimento de c´elulas T activas que chegan ao tumor. Por outro lado, a inmunoterapia mod´elase inclu´ındo d´uas variables m´ais que describen a concentraci´on dos f´armacos anti-CTLA4 e anti-PD(L)1, respectivamente, con taxas de eliminaci´on exponenciais. Estas variables modulan os termos de activaci´on das c´elulas T ou da morte inmune das c´elulas tumorais. 6 Resumo O modelo de resposta presentado validouse co axuste de datos precl´ınicos procedentes de dous traballos experimentais publicados. Neles comparouse a evoluci´on de volumes tumorais en poboaci´ons de ratos para distintos grupos: grupo control en ausencia de tratamento, s´o radioterapia, s´o inmunoterapia con anti-PDL1 ou antiCTLA4, e tratamento combinado de radioinmunoterapia co respectivo anticorpo. Para obter os valores dos par´ametros do modelo que axustan as medias poboacionais dos distintos grupos, empregouse un algoritmo de optimizaci´on (Simulated Annealing). Posteriormente, tomando como referencia os valores dos par´ametros obtidos do axuste dos volumes precl´ınicos, lev´aronse a cabo estudos sobre o calendario de administraci´on da radioimmunoterapia con anti-CTLA4 que maximiza a probabilidade de control tumoral. Nun primeiro estudo fix´aronse os d´ıas de administraci´on da radioterapia e vari´aronse os da inmunoterapia. Noutro segundo estudo, fix´aronse os d´ıas de administraci´on da inmunoterapia e variouse o n´umero de fracci´ons de dose de radioterapia. Os resultados desta investigaci´on, en li˜na cos dos traballos experimentais, apuntan a que administrar a inmunoterapia dous d´ıas despois de comezar coa radioterapia proporciona maiores taxas de control que comezar o mesmo d´ıa ou retrasala catro d´ıas, sendo a ´ultima a que aporta peores resultados. Con base no estudo realizado, hipotet´ızase que os retardos temporais asociados aos desprazamentos de ant´ıxenos e c´elulas T activas pode ter relevancia ´a hora de dese˜nar os tratamentos. Espaciar o inicio da inmunoterapia con respecto ´a radioterapia pode dar lugar a que o n´umero de c´elulas T activas no tumor non se mante˜na e dimin´ua o seu efecto terap´eutico. Por outro lado, obt´ıvose que para o caso de estudo a radioterapia moderadamente hipofraccionada leva a maiores taxas de control tumoral que os tratamentos convencionais ou extremadamente hipofraccionados. A radioterapia FLASH consiste na administraci´on da radiaci´on con altas taxas de dose, superiores aos 40 Gy s−1en comparaci´on cos 0.05-0.40 Gy s−1da terapia convencional. Esta terapia cobrou importancia nos ´ultimos anos xa que os estudos precl´ınicos realizados apuntan a que mant´en a efectividade da modalidade convencional aportando menos toxicidade aos tecidos sans veci˜nos. Estudos experimentais en soluci´ons acuosas e outros compostos que imitan as c´elulas, as´ı como estudos precl´ınicos in vivo, demostraron que a radioterapia FLASH provoca unha diminuci´on de os´ıxeno durante o tratamento nos medios irradiados, ao contrario que a convencional, onde este efecto non ´e observable coa tecnolox´ıa actual. Este fen´omeno pode deberse a que o consumo de os´ıxeno asociado a taxas de doses baixas se compensa coa continua aportaci´on do sistema vascular tumoral, que resulta insuficiente no caso da radioterapia FLASH. Dado que os niveis reducidos de os´ıxeno aportan maior resistencia ao dano por radiaci´on, ´e amplamente aceptado que este efecto causa unha diminuci´on da toxicidade do tratamento no tecido san, as´ı como outros factores de tipo inmune ou fisicoqu´ımico. Pola contra, non est´a claro por que ese efecto protector non se observa no tumor, o cal levar´ıa a unha menor efectividade coa radioterapia FLASH respecto da convencional. Mentres que diversos estudos atrib´uen este fen´omeno a diferenzas entre os tecidos tumoral e san, como na produci´on de radicais libres, outros suxiren que 7 Isabel Gonz´ alez Crespo o efecto protector no tumor existe pero non ten evidencia significativa nas pequenas mostras experimentais que se adoitan empregar. Nesta tese invest´ıgase a aparente isoefectividade de ambas terapias prestando atenci´on ao papel do os´ıxeno como elemento diferenciador. Para elo, inicialmente empregouse un modelo de osixenaci´on obtido da literatura para axustar datos de eliminaci´on de os´ıxeno en distintas soluci´ons e tumores precl´ınicos. Este modelo consta dunha ecuaci´on en derivadas parciais de reacci´on-difusi´on que describe a din´amica do os´ıxeno no tecido tumoral, inclu´ındo o sistema vascular como termo fonte, e as c´elulas e a radiaci´on como termos de consumo. Para ter en contra a heteroxeneidade que amosan os tumores en canto ´as distribuci´ons de os´ıxeno, os modelos de osixenaci´on resolv´eronse sobre un dominio bidimensional que simula o tecido tumoral cunha distribuci´on aleatoria de capilares que simula o seu sistema vascular. Do mesmo traballo obt´ıvose unha expresi´on para calcular a fracci´on de supervivencia asociada ´a radioterapia FLASH baseada no modelo LQ coa modificaci´on dos OERs. Ademais, dese˜nouse un modelo que describe a evoluci´on do volume tumoral despois do tratamento partindo das fracci´ons de supervivencia obtidas. O modelo de resposta presentado consta de d´uas ecuaci´ons diferenciais con impulsos que describen respectivamente a din´amica das c´elulas tumorais viables e das c´elulas danadas que son eliminadas por efecto do tratamento. Este modelo axustouse a datos experimentais de volumes precl´ınicos empregando o algoritmo de optimizaci´on Simulated Annealing. Tomando como referencia os valores ´optimos obtidos para os par´ametros que describen o modelo, xerouse unha poboaci´on cun gran n´umero de tumores simulados, divididos en grupos control, radioterapia FLASH e radioterapia convencional, e estud´aronse as diferenzas en volume e TCP entre os distintos grupos. O estudo suxire, de acordo con experimentos precl´ınicos, que as diferenzas en volume entre os grupos asignados a ambas terapias non son significativas na maior´ıa dos casos, polo que a radioterapia FLASH poder´ıa semellar isoefectiva respecto da convencional. A continuaci´on, estudouse a isoefectividade en termos da TCP. Fix´eronse dous estudos independentes sobre a poboaci´on de tumores simulados para analizar posibles factores que infl´uen no control tumoral. Inicialmente asign´aronse tres cocientes α/β distintos (3, 10 e 20) para estudar como afecta a radiosensibilidade das c´elulas ´a perda de control tumoral, atopando que a maior diferenza entre as radioterapias FLASH e convencional se d´a para tumores con cocientes baixos. Isto pode explicarse ´a vista da ecuaci´on do modelo LQ coa modificaci´on dos OERs, na cal se ve que o par´ametro β (cadr´atico) engade maior sensibilidade aos cambios de osixenaci´on. Por outro lado, tomando como referencia o caso para α/β = 10 Gy, dividiuse a poboaci´on en funci´on da s´ua osixenaci´on mediana en tres grupos: osixenaci´on baixa (menor de 10 mmHg), media (entre 10 e 20 mmHg) ou alta (maior de 20 mmHg). Obt´ıvose que a maior perda de TCP se d´a nos tumores m´ais osixenados. Pola contra, cando se analiza a supervivencia celular das c´elulas para distintos niveis de os´ıxeno, a maior diferenza entre ambas terapias d´ase para osixenaci´ons baixas, arredor dos 2–6 mmHg. Esta aparente contradici´on est´a causada pola heteroxeneidade nas osixe8 Resumo naci´ons tumorais. As zonas de baixa osixenaci´on est´an presentes tam´en en tumores “ben osixenados”e son as que marcan a TCP, ao ser as m´ais resistentes ao tratamento. Dado que as doses necesarias para acadar unha certa TCP coa terapia convencional son menores en tumores ben osixenados, poden resultar insuficientes para eliminar as c´elulas que entran en hipoxia como consecuencia da ca´ıda de os´ıxeno durante a terapia FLASH, mentres que as doses m´ais elevadas asociadas aos tumores ´oxicos seguen sendo efectivas para as c´elulas que sofren dita transici´on. Polo tanto, ´e posible que na incorporaci´on da radioterapia FLASH ao ´ambito cl´ınico se deban aumentar as doses administradas con respecto ´a terapia convencional para acadar a mesma efectividade en determinados casos. En conclusi´on, esta tese presenta dous modelos matem´aticos de resposta tumoral a t´ecnicas de radioterapia avanzada que xurdiron nos ´ultimos anos, a radioinmunoterapia e a radioterapia FLASH. Responde as´ı ´a necesidade de desenvolver novos modelos para estudar o efecto dos tratamentos e explorar os mecanismos asociados aos mesmos que non son efecto directo do dano celular por radiaci´on. Os modelos cl´asicos como o LQ e as s´uas variantes non son suficientes para describir mecanismos de morte indirecta que involucran m´ultiples axentes do microentorno tumoral ou outros efectos biol´oxicos e fisicoqu´ımicos asociados ´a terapia, non sendo ´optimos para a planificaci´on dos tratamentos no ´ambito cl´ınico. Os modelos presentados valid´aronse co seu axuste a datos experimentais e empreg´aronse para estudar os mecanismos subxacentes ´a efectividade dos tratamentos. En canto ´a radioinmunoterapia, os resultados apuntan a que os retardos biol´oxicos asociados ao transporte das c´elulas T desde os ´organos linfoides cara o tumor poden ter un papel relevante na planificaci´on dos tratamentos. Por outro lado, viuse que as curvas de crecemento tumoral poden ser insuficientes para establecer a isoefectividade das radioterapias FLASH e convencional en termos da TCP, de maior interese na pr´actica cl´ınica. As variaci´ons de os´ıxeno asociadas ´a radioterapia FLASH poden derivar nun aumento da fracci´on de supervivencia tumoral, que non leva a diferencias significativas nas curvas de volume (para os tama˜nos mostrais empregados nos estudos precl´ınicos) pero trasl´adase ´a TCP, proporcionando taxas de control m´ais baixas que as da modalidade convencional. 9 Chapter 1 Introduction “No one told me there would be math!” — Bart Simpson, The Simpsons: Mathlete’s Feat. Since the inception of radiation as a cancer therapy, mathematical models have emerged as a valuable tool for designing treatments and predicting their outcomes. This thesis aims to present novel mathematical models that may contribute to a better understanding and planning of advanced radiotherapy treatments that have emerged in recent years, such as radioimmunotherapy and FLASH radiotherapy. This introduction provides a historical overview of the use of radiation in cancer treatment, from its first clinical applications in the early-20th century to the development of more advanced techniques used today. Additionally, this chapter offers insights into classical mathematical models employed since the 1950s, which have assisted oncologists and medical physicists in designing clinical treatments. These models constitute the foundation for developing more complex and suitable mathematical models to describe the effect of novel radiation treatments like radioimmunotherapy and FLASH radiotherapy; these techniques are introduced later in this chapter. In the final section, the motivation and objectives of this thesis are stated. Isabel Gonz´ alez Crespo being αox and βox the αand βparameters under fully aerobic conditions; OERαand OERβthe maximum oxygen enhancement ratios; and km(mmHg) the oxygen partial pressure at which OERs achieve the half-maximum value. The effect of the OERs tends to be zero as the oxygen pressure increases, thus the α(p) and β(p) values given by equations (1.5) and (1.6) approximate to αox and βox, respectively. Analogously, equation (1.4) tends to equation (1.1) under high oxygen levels. Typically, OERαand OERβvalues are within the range 2.50–3 and kmis set to 3.28 mmHg [42]. Figure 1.2 shows SF–dose curves obtained for different oxygen statuses together with the response under aerobic conditions obtained through the LQ model. 0 5 10 15 20 10−6 10−5 10−4 10−3 10−2 10−1 100 Dose (Gy) SF Fully aerobic p=5 mmHg p=10 mmHg p=20 mmHg p=40 mmHg Figure 1.2: Surviving Fraction (SF) versus dose curves obtained from equation (1.4) for different oxygen partial pressures, p, and equation (1.1) under aerobic conditions, with the parameters αox = 0.20 Gy−1,βox = 0.02 Gy−2,OERα= 2.50, OERβ= 2.50 and km= 3.28 mmHg. The vertical axis is on a logarithmic scale. 1.3.2 Linear-Quadratic-Linear (LQL) model The Linear-Quadratic-Linear (LQL) model [32] is an extension of the LQ model (1.1) that includes a term that compensates for the LQ underestimation of sublethal damage repair. The SF obtained with the LQL model is given by: SF(d) = exp−αd −2β γ2(γd −1 + exp (−γd)),(1.7) where γ(Gy−1) is an extra parameter that modules the slope of the curve. The above expression results from considering that the repair of sublethal damage, represented by the quadratic term, depends on the delivered dose. Thus, the maximum 18 Chapter 1. Introduction effect of radiation is given by the LQ model, but when considering sublethal damage repair, equation (1.1) becomes: SF(d) = exp−αd −Gβd2,(1.8) where Gis named the dose protraction factor [32,43], given by: G=2 d2Z∞ 0 ˙ D(t)dt Zt 0 exp (−λ(t−t′)) ˙ D(t′)dt′,(1.9) being λthe repair rate, and ˙ Dthe dose rate. For T, the delivery time of the dose d, and a constant dose rate, it follows that ˙ D=d/T for 0 ≤t≤T. Thus, integrating the above expression gives: G=2 (λT + exp (−λT)−1) (λT)2,(1.10) that considering γd =λT, becomes equation (1.7). As shown in Figure 1.3, when the dose increases, the calculated SF from the LQL model is higher compared to the LQ model. 0 5 10 15 20 10−6 10−5 10−4 10−3 10−2 10−1 100 Dose (Gy) SF LQL LQ Figure 1.3: Surviving Fraction (SF) versus dose curves obtained from the Linear-Quadratic (LQ) and Linear-Quadratic-Linear (LQL) models with parameters α= 0.20 Gy−1and β= 0.02 Gy−2in both models, and γ= 0.05 Gy−1in the LQL model. The vertical axis is on a logarithmic scale. 1.3.3 Biologically Effective Dose (BED) The Biologically Effective Dose (BED) is a concept used to compare the biological effects of different radiation fractionation schemes, taking into account both physical 19 Isabel Gonz´ alez Crespo magnitudes, such as the dose per fraction, d, and the total dose, D; and biological aspects, like the α/β ratio. Based on the LQ model, the Biologically Effective Dose (BED) is calculated as follows [44,45]: BED =D1 + d α/β .(1.11) Tumors with low α/β ratios are more sensitive to fractionation as shown in Figure 1.4. Thus, these tumors benefit less from increasing the number of fractions. 1 2 3 4 5 6 7 8 9 10 20 40 60 80 100 120 140 160 Number of fractions BED α/β = 3 Gy α/β = 10 Gy α/β = 20 Gy Figure 1.4: Biologically Effective Dose (BED) curves for different α/β ratios obtained from equation (1.11) with a total dose D= 20 Gy varying the number of fractions. Using the BED facilitates the search for optimal treatments that give the desired therapeutic effect while minimizing damage to surrounding healthy tissues. However, two different treatments with the same BED can promote different outcomes due to particular aspects of the TME. 1.3.4 Tumor Control Probability (TCP) The Tumor Control Probability (TCP) is a concept used to quantify the likelihood that a RT treatment effectively controls or eliminates a tumor. The most common approach [46] to calculate TCP assumes that a single tumor cell has the potential to proliferate and derive tumor regrowth, known as the clonogenic cell hypothesis [47]. Moreover, this model assumes that cells are eliminated randomly and independently, with each cell having a survival probability SF. From the above considerations, TCP is defined as the probability that all tumor cells are eliminated. Let Nbe the variable describing the number of cells that survive the treatment, and N0be the initial number of cancer cells within a tumor. Considering that N 20 Chapter 1. Introduction follows a Poisson distribution with S=N0SF as the expected number of surviving cells, it follows: P(N) = exp(−S)SN N!.(1.12) Thus, TCP is given by the probability of having no surviving cells, resulting in: TCP =P(0) = exp (−N0SF).(1.13) Particularly, the LQ-Poisson methodology to calculate the TCP based on the LQ model (1.1) gives: TCP(d) = exp −N0exp −αd −βd2.(1.14) 1.4 Advanced techniques of RT The technological advances of the last decades enabled the development of new RT modalities that aim at maximizing the radiation effect on tumor cells while decreasing the toxicity on surrounding organs and tissues. This thesis focuses on two techniques introduced in this section: the combination of RT with ImmunoTherapy (IT), namely RadioImmunoTherapy (RIT), and FLASH-RT. 1.4.1 ImmunoTherapy (IT) In contrast to traditional cancer treatments such as surgery, chemotherapy, and RT, which directly target the tumor, IT aims to stimulate and enhance the natural ability of the immune system to recognize and eliminate cancer cells [48,49]. The immune response comprises a sequence of steps where immune cells identify and eliminate foreign harmful agents. In the context of cancer, this process is referred to as the immune cycle of cancer [50], detailed later in this section. In the course of the immune response, the organism relies on regulatory mechanisms known as immune checkpoints at different stages of the cycle. These checkpoints prevent excessive immune reactions that may otherwise lead to inflammatory diseases. However, certain tumors exploit them to counteract the immune response and evade anti-tumor defenses. The most extensively studied immune checkpoints are the Cytotoxic T-Lymphocyte Antigen 4 (CTLA-4) receptor, found on the surface of immune cells and playing a role in their activation; and the Programmed Death 1 (PD-1)/Programmed Death-Ligand 1 (PD-L1) pathway, which significantly influences the identification of target cells. The history of IT dates back to 1891 when Dr. William Coley, considered the father of IT, observed that some cancer patients experienced tumor regression after bacterial infections. Inspired by these findings, he developed a treatment using bacteria which, although showed some success, had limited acceptance by oncologists due 21 Isabel Gonz´ alez Crespo to inconsistent results. While several IT advancements emerged in the latter half of the 20th century, it was not until the 2010s that Immune Checkpoint Inhibitors (ICIs) demonstrated significant success in treating diverse cancers, particularly when combined with RT [51]. The immunity cycle of cancer The cancer immunity cycle [50] is a recurring process where the immune system identifies and responds to harmful agents within the organism. This complex response involves components of the TME, such as tumor cells, antigens (molecular structures enabling the immune system to discern the presence of cancer cells), and immune cells, like T-cells. Immune cells can directly target and eliminate tumor cells displaying antigens, being a defense mechanism against the development and progression of cancer. Figure 1.5 illustrates the different phases of the anti-tumor immune response. Initially, the presence of tumor cells promotes the release of antigens, which are different from those released by cells infected with other pathogens (step 1). Then, Antigen Presenting Cells (APCs), particularly Dendritic Cells (DCs) capture these antigens and migrate to lymphoid organs (step 2), where T-cells identify them and become active against the specific presented antigen, triggering the immune response (step 3). Subsequently, activated T-cells migrate towards the tumor (step 4) and infiltrate it (step 5). Finally, T-cells recognize the tumor cells and eliminate them (step 6), resulting in the release of new antigens, thereby initiating the immunity cycle again. Immune Checkpoint Inhibitors (ICIs) As previously mentioned, the organism relies on certain mechanisms to regulate the immune response and prevent excessive reactions that could lead to inflammatory diseases. Cancer takes advantage of these mechanisms to evade the anti-tumor immune response. Two examples of scape tactics are the stimulation of CTLA-4 receptor on the surface of T-cells and the PD-1/PD-L1 pathway. The CTLA-4 receptor is involved in step 3 of the immunity cycle, where T-cells are activated against a specific antigen presented on the surface of APCs. As depicted in Figure 1.6, T-cell activation requires at least two signals [52]. Initially, the T-Cell Receptor (TCR) binds to the antigen sustained in a structure known as Major Histocompatibility Complex (MHC) on the surface of APCs, initiating T-cell activation. Subsequently, the co-stimulatory molecule Cluster of Differentiation (CD)28 must interact with one of its ligands, CD80 or CD86, to complete the activation process. During the activation procedure, the CTLA-4 receptor matures on the cell surface and competes with CD28 for binding to its ligands. If CTLA-4 binds to CD28, it inhibits T-cell activation, acting as an ICI to prevent premature T-cell overactivation. 22 Chapter 1. Introduction Figure 1.5: Phases of the immunity cycle of cancer response. 1) Tumor cells release antigens; 2) Dendritic Cells (DCs) take antigens to activation sites (lymphoid nodes); 3) T-cells against tumor cells are activated; 4) T-cells migrate towards the tumor and 5) infiltrate it; 6) T-cells attack tumor cells and kill them. [Created in BioRender.com.] In 2011, ipilimumab became the first approved human antibody against the receptor CTLA-4 [52]. Treatments with anti-CTLA-4 (αCTLA4) involve intravenous administration of the antibody, which binds to the CTLA-4 receptor, preventing its interaction with CD28 and promoting T-cell activation. Figure 1.6: Step 3 of the immunity cycle. T-cell activation requires two signals: the binding of T-Cell Receptor (TCR) and the Major Histocompatibility Complex (MHC) sustaining the antigen, and the union of the Cluster of Differentiation (CD), which may be modulated by the Cytotoxic TLymphocyte Antigen 4 (CTLA-4) receptor and ImmunoTherapy (IT) with anti-CTLA-4 (αCTLA4). [Created in BioRender.com.] 23 Isabel Gonz´ alez Crespo In contrast, the PD-1/PD-L1 pathway is involved in step 6 of the immunity-cycle, as depicted in Figure 1.7.PD-1 is a protein on the surface of T-cells, and its ligand, PD-L1, is present on the surface of DCs and other cells in the TME. When PD-1 binds to PD-L1, the T-cell becomes “deactivated”, serving as another ICI to avoid an excessive immune response. Tumor cells can present PD-L1 on their surface, using it as an escape mechanism [53]. The first anti-PD-1 (αPD1) and anti-PD-L1 (αPDL1) drugs (generally named as anti-Programmed Death-(Ligand) 1 (αPD(L)1)) emerged in the last decade. In 2014, the monoclonal antibodies nivolumab and pembrolizumab were approved to counteract the PD-1/PD-L1 evasion pathway. Similar to αCTLA4 drugs, αPD(L)1 antibodies are administrated intravenously and bind to their respective ligand to prevent T-cell deactivation. Figure 1.7: Step 6 of the immunity cycle. T-cells detect tumor cells regulated by the Programmed Death 1 (PD-1)/Programmed Death-Ligand 1 (PD-L1) pathway and ImmunoTherapy (IT) with antiProgrammed Death-(Ligand) 1 (αPD(L)1). [Created in BioRender.com.] 1.4.2 RadioImmunoTherapy (RIT) RT is known to promote an immune anti-tumor response in the organism, although the mechanisms behind this effect are not fully understood. Radiation acts like an in situ vaccine that increases the recruiting of T-cells in the tumor site, motivated by the release of specific signals (antigens) [54]. It can even influence an immune response far from the irradiated field, contributing to eliminating metastasic tumors [55]. This phenomenon has been known as the abscopal effect since 1953. Recently, IT has arisen as a promising complement to RT, leading to RIT treatments [56,57]. In the context of cancer treatment, IT pursues the goal of stimulating or enhancing the immune response to target and eliminate cancer cells. As mentioned 24 Chapter 1. Introduction in the previous section, the use of IT drugs with ICIs blocks inhibitory signals in the immune system, allowing T-cells to recognize and attack cancer cells more effectively. Particularly, IT with αCTLA4 and αPD(L)1 drugs enhances the natural anti-tumor response carried by the immune system, improving RT outcomes [18,58,59]. However, the success of RIT depends on patient-specific features, like the good behavior of their immune system [58]. Moreover, an excessive immune response can damage normal cells and promote inflammatory diseases, increasing toxicity. For these reasons, RIT continues to be a field of ongoing research. The synergy between RT and IT RT not only induces an anti-tumor immune response, which can be enhanced by RIT, but also leads to some immunosuppressive effects. Radiation damages both tumor cells and T-cells present within the tumor at the irradiation time, as well as the vascular system which serves as the infiltration way of immune cells and, if damaged, can lead to hypoxia and radioresistance (although the effect of vascular damage on tumor control has been debated [11,14,21–23]). Hypofractionated RIT has been demonstrated to improve tumor control rates in preclinical experiments compared to the independent use of RT and IT as monotherapies [58,59]. One of the reasons behind this may be that in conventional RT, when T-cells infiltrate the tumor days after the first RT fraction, the tumor continues to be irradiated, and T-cells are damaged; while in hypofractionated regimens, T-cells that progressively infiltrate the tumor are exposed to fewer radiation fractions. Previous preclinical research has studied the effect of different fractionations on the immune cells, finding that hypofractionated treatments increase T-cell accumulation within the tumor a few days after the first fraction compared to more fractionated procedures [15, 19,58]. However, this phenomenon must be further investigated. Finding the combination of both therapies that maximizes their synergistic effect is not straightforward and it is necessary to find new approaches that help in designing optimal treatments. Modeling RIT requires the necessary description of indirect immune-mediated effects and the two-faced effect of radiation as an immunosuppressive (damages circulating T-cells and tumor vasculature) and immunostimulatory factor (increases the recruiting of T-cells in the tumor). Thus, simple models, such as the LQ model and its derived forms presented in the previous sections, are insufficient to explain the complex interactions of radiation and IT drugs with the immunity cycle of cancer. In recent years, some research studies have presented new mathematical models in RIT, attempting to explore the mechanisms behind its effect and to assist in treatment planning [60]. 25 Isabel Gonz´ alez Crespo 1.4.3 FLASH RadioTherapy (FLASH-RT) FLASH-RT is a novel treatment strategy that consists of delivering radiation doses at UHDRs, exceeding the 40 Gy s−1instead of the 0.05–0.40 Gy s−1used in conventional RT (conv-RT) [61,62]. In 1959, Dewey and Boag [63] reported the FLASH effect for the first time when comparing the radiation sensitivity to conv-RT and FLASH-RT in bacteria, finding that UHDRs led to lower damage. Similar results were observed later in mammalian cells and preclinical trials [64,65]. More recent in vivo studies have shown the potential of FLASH-RT to maintain the effectiveness of conv-RT on tumors while sparing healthy surrounding tissue [66–74]. Considering these findings, FLASH-RT has emerged as a promising RT modality for future clinical practice. The biological mechanisms involved in FLASH-RT are very complex [64] and constitute an active field of research. Many studies have suggested that the protective effect of healthy tissue arises from the Radiolytic Oxygen Depletion (ROD) process caused by FLASH irradiation [62,75–78], which increases cell radioresistance. As mentioned in Section 1.3.1, it is established that radiation can deplete oxygen in irradiated cells and tissues through the radiolysis of water molecules. The use of UHDRs promotes a quick ROD process in FLASH-RT. On the contrary, in conv-RT, the tumor vascular system continuously resupplies oxygen, preventing observable oxygen depletion. As discussed earlier in this chapter, hypoxic cells may be less responsive to RT, a phenomenon described by the OERs, presented in Section 1.3.1. That is mostly accepted as the reason behind the sparing effect of FLASH-RT in healthy tissues [65]. However, other research works propose that ROD alone is insufficient to fully explain the sparing effect. Instead, they attribute it to other factors like immune effects [61,79]. The short exposure to radiation in FLASH-RT may significantly reduce the proportion of circulating immune cells that are irradiated and damaged, leading to a more effective immune system capable of repairing radiation-induced damage to normal tissue and contributing to reducing inflammation [65]. Understanding why the protective effect is not observed in tumors remains an ongoing area of research. Some studies suggest that oxygen depletion has a higher protective effect on healthy tissue as tumors generally present lower oxygen levels, limiting the benefits derived from ROD. Nevertheless, other studies propose that additional physicochemical factors may contribute to this phenomenon [61,65], such as variations in reactive species production between tumor and non-tumor cells [80], the recombination of free radicals [81], or radiation-induced immune effects [79]. In recent years, theoretical work has accompanied experimental studies in an attempt to elucidate the underlying mechanisms of the FLASH effect [81–84]. 26 Chapter 1. Introduction 1.5 Motivation and objectives Mathematical models have long been paramount in analyzing the response of cells and tissues to radiation. In this regard, the LQ model has been widely used to design novel RT fractionations that have improved effectiveness and toxicity. However, in recent years, new RT techniques have emerged, such as RIT and FLASH-RT. These novel treatments involve cell death mechanisms not directly caused by radiation damage. Therefore, the LQ and LQ-derived models cannot fully describe the treatment outcome. These limitations motivate the development of new mathematical models, such as those presented in this thesis. The main objective of this work is to develop novel tumor response models in RIT and FLASH-RT that serve to study the mechanisms behind their anti-tumor effect. The specific objectives are: •To collect experimental data from the literature to validate the proposed models. •To model the effect of biological delays on RIT, associated with the activation of immune cells and their infiltration within tumor tissue. •To propose hypotheses regarding the optimal administration schedules for treatments combining RT and IT drugs. •To model the effect of oxygen depletion during FLASH-RT on the treatment outcome. •To argue about the circumstances under which FLASH-RT and conventional RT may demonstrate iso-effectiveness in clinical tumor control. 1.6 Outline The rest of the document is structured as follows: •Chapter 2presents the mathematical and computational methods used in chapters 3and 4to design new mathematical models in RIT and FLASH-RT. •Chapter 3describes the specific methods and the obtained results in RIT. •Chapter 4describes the specific methods and the obtained results in FLASH-RT. •Chapter 5contains the conclusions of this doctoral thesis. 27 Isabel Gonz´ alez Crespo numerical solution of the oxygenation problem. Typically, tumor oxygenation is modeled by a reaction-diffusion PDE, which takes account of the oxygen supply through the vascular system and metabolic consumption by tumor cells. Furthermore, in the context of FLASH-RT, an additional consumption term associated with ROD is considered. The resulting oxygenation problem is addressed using the Finite Element Method (FEM) after defining the geometry of the problem and establishing initial and boundary conditions. 2.2.1 Mathematical model Oxygen partial pressure, p, in tumor and healthy tissue is typically modeled by the following reaction-diffusion equation: ∂p ∂t (t, x) = DO2∆p(t, x)−g(p(t, x)),(2.8) where (t, x) are the temporal-spatial coordinates, ∆ the Laplacian operator with respect to x,DO2the oxygen diffusion coefficient and ga consumption term resulting from the addition of cellular metabolic activity [37,93–96] and ROD [84], g=gmet +gROD, defined later in Sections 2.2.5 and 2.2.6, respectively. 2.2.2 Geometry of the problem Tumors present a vascular network of capillaries that provide cells with oxygen and other nutrients needed to carry on their metabolic activity. The Vascular Fraction (VF) is defined as the fraction of the total tissue occupied by blood vessels. Following the work of D´ıaz-Botana [92], the VF is assumed to be homogeneous within the tumor, and oxygen distribution is described by a representative voxel (three-dimensional equivalence of a pixel) of tumor tissue. Moreover, to simplify the computational complexity of simulating oxygenation on three-dimensional vascular networks, the problem is simplified to a two-dimensional domain, which is a good surrogate of the threedimensional problem for particular vascular networks [93]. Nevertheless, it must be noted that this approximation can miss some effects derived from the chaotic nature of tumor vascularity, as the two-dimensional simplification assumes parallel capillaries of equal length. Following the previous considerations, the oxygenation problem is defined on a 1×1 mm2pixel with inner circles, representing tumor tissue and its vascular system, as shown in Figure 2.1. Capillaries are randomly located within the pixel following the next steps for the placement of each circle: Step 1: Obtain the coordinates of the center using a uniform distribution. Step 2: Assign a diameter characterized by a lognormal distribution obtained from the literature [97] and detailed in Table 2.1. 34 Chapter 2. Methods and materials Step 3: Verify that the capillary is entirely within the pixel. Otherwise, reject it and return to Step 1. Step 4: Verify that the capillary fits the minimum inter-vessel distance given in Table 2.1 with the previously positioned capillaries. Otherwise, reject it and return to Step 1. These steps are repeated until the cumulative area occupied by capillaries reaches a given VF. Figure 2.1: Example of a two-dimensional tumor domain. Geometry for the numerical simulation of the oxygenation problem. The inner circles represent blood capillaries of the tumor vascular system. This geometry corresponds to a VF vf= 0.14. Size (µm) Mean diameter 64.5 Standard deviation 58.0 Minimum diameter 0.9 Maximum diameter 161.9 Minimum inter-vessel distance 19.3 Table 2.1: Characteristics of the vascular vessels within the tumor, corresponding to experimental measurements of capillaries in the center of colorectal tumors [97]. Mean and standard deviation refer to the lognormal distribution. 35 Isabel Gonz´ alez Crespo 2.2.3 Initial and boundary conditions Imposing initial and boundary conditions, the following oxygenation problem is defined on a two-dimensional domain Ω with boundary Γ = Γcap ⊔Γout: Problem 2.1. Find p: [t0, tf]×Ω→[0, pcap]such that:                  ∂p ∂t (t, x) = DO2∆p(t, x)−g(p(t, x)),∀(t, x)∈[t0, tf]×Ω, p(t0,x) = p0(x),∀x∈Ω, p(t, x) = pcap,∀(t, x)∈[t0, tf]×Γcap, ∂p ∂n(t, x) = 0,∀(t, x)∈[t0, tf]×Γout, where p0is the partial pressure at t0. On the outer boundary of the domain, Γout, zero-flux Neumann conditions were imposed, while Dirichlet conditions with p=pcap were set on the capillary walls, Γcap. Although real tumors present heterogeneous values on the partial pressure for different capillaries, a constant pcap = 40 mmHg is the typical value used in the literature [37,84,94,95]. Similarly, oxygen diffusion properties are not homogeneous within tumor tissue, however a constant diffusion coefficient DO2= 2 ×103µm2s−1 is usually accepted [94,95]. 2.2.4 Numerical solution: FEM In this study, the oxygenation problem 2.1 was solved by the FEM using a MATLAB (The Mathworks Inc., Natick, MA) solver from the Partial Differential Equation Toolbox (version 3.9 (R2022b)). This section presents the theoretical framework for the application of the FEM. The fundaments of the use of the FEM for time-dependent problems can be found in the literature [98]. Equation (2.8) can be turned into its integral form by constructing the weak formulation. Let vbe a test function in the space V={v∈H1(Ω)/v|Γcap = 0}. Multiplying equation (2.8) by v, it follows: ∂p ∂t (t, x)v(x) = DO2∆p(t, x)v(x)−g(p(t, x))v(x). Then, integrating the above equation over Ω yields: ZΩ ∂p ∂t (t, x)v(x)dΩ = ZΩ DO2∆p(t, x)v(x)dΩ−ZΩ g(p(t, x))v(x)dΩ, 36 Chapter 2. Methods and materials where the Green’s theorem can be applied to the second term, obtaining: ZΩ ∂p ∂t (t, x)v(x)dΩ = −ZΩ DO2∇p(t, x)·∇v(x)dΩ + ZΓ DO2 ∂p ∂n(t, x)v(x)dΓ −ZΩ g(p(t, x))v(x)dΩ. Also notice that, as v= 0 on Γcap and ∂p ∂n= 0 on Γout, it follows: ZΓ DO2 ∂p ∂n(t, x)v(x)dΓ = ZΓcap DO2 ∂p ∂n(t, x)v(x)dΓcap +ZΓout DO2 ∂p ∂n(t, x)v(x)dΓout = 0. Then, for DO2constant, it results: ZΩ ∂p ∂t (t, x)v(x)dΩ = −DO2ZΩ∇p(t, x)·∇v(x)dΩ−ZΩ g(p(t, x))v(x)dΩ. From the previous considerations, the oxygenation problem 2.1 can be redefined as the following variational problem: Problem 2.2. Find a function p(t, ·)for t∈[t0, tf]on V={v∈H1(Ω)/v|Γcap = 0} such that p(t0,x) = p0(x), and ∀v∈ V, t ∈(t0, tf): ZΩ ∂p ∂t (t, x)v(x)dΩ + DO2ZΩ∇p(t, x)· ∇v(x)dΩ + ZΩ g(p(t, x))v(x)dΩ = 0. The notation can be simplified by defining p(t, ·) = p(t) and introducing the scalar product in L2(Ω), and the bilinear form a(v, w) for v, w ∈ V as: ⟨v, w⟩L2(Ω) =ZΩ v(x)w(x)dΩ, a(v, w) = ZΩ∇v(x)·∇w(x)dΩ. Moreover, as neither Ω nor vdepends on t, it follows: ZΩ ∂p ∂t (t, x)v(x)dΩ = d dt ZΩ p(t, x)v(x)dΩ. Thus, the weak formulation in Problem 2.2 is rewritten as: find p: [t0, tf]→ V such that p(t0) = p0and ∀v∈ V, t ∈(t0, tf): d dt⟨p(t), v⟩L2(Ω) +DO2a(p(t), v) + ⟨g(p(t)), v⟩L2(Ω) = 0. 37 Isabel Gonz´ alez Crespo Regarding the discretization in the spatial coordinate, to apply the FEM the infinite-dimensional space Vis approximated by a finite-dimensional subspace Vh, so it is sufficient to find ph: [t0, tf]→ Vhsuch that ph(t0) = p0,h, where p0,h ∈ Vhis an approximation of the initial solution p0, and ∀vh∈ Vh, t ∈(t0, tf): d dt⟨ph(t), vh⟩L2(Ω) +DO2a(ph(t), vh) + ⟨g(ph(t)), vh⟩L2(Ω) = 0. Let {ϕi}1≤i≤Nbe a basis for the finite-dimensional space Vh. Then, the approximate solution and initial condition can be written as: ph(t) = N X i=1 Ph i(t)ϕi, p0,h = N X i=1 P0,h i(t)ϕi, where Ph= (Ph i)1≤i≤Nand P0,h = (P0,h i)1≤i≤Nare the coordinates vectors of phand p0,h, respectively. Subsequently, values of (Ph i)1≤i≤Nand (P0,h i)1≤i≤Nmust be found, such that Ph j(t0) = P0,h jfor 1 ≤j≤N, and: N X i=1⟨ϕi, ϕj⟩L2(Ω) dPh i(t) dt +DO2 N X i=1 a(ϕi, ϕj)Ph i(t) + ⟨g N X i=1 Ph i(t)ϕi!, ϕj⟩L2(Ω) = 0. Introducing the mass matrix, Mh, and stiffness matrix, Kh, defined by: (Mh)i,j =⟨ϕi, ϕj⟩L2(Ω),(Kh)i,j =a(ϕi, ϕj), for 1 ≤i, j ≤N, the weak formulation is equivalent to the system:      Mh dPh dt (t) + DO2KhPh(t) + bh(t)=0,∀t∈(t0, tf), Ph(t0) = P0,h, being: bh(t) = ⟨g N X i=1 Ph i(t)ϕi!, ϕj⟩L2(Ω). Lastly, the problem is fully discretized in time-space by using a finite difference method, such as the forward Euler scheme (2.6), having: Mh Ph,k+1 −Ph,k ∆t+DO2KhPh(tk) + bh(tk)=0, where ∆tis the time step of the temporal discretization for the interval [t0, tf]. 38 Chapter 2. Methods and materials The FEM is characterized by three elements: a partition/triangulation of the domain Ω, τ; a set of points obtained from the partition named nodes, Σ; and the basis functions {ϕi}1≤i≤N. A partition τdivides Ω in a finite number, N, of subdomains, Ti, such as ∀Ti∈τ: •Tiis closed and its interior set ˚ Ti=∅is connected. •The boundary is Lipschitz continuous. •Ω = ∪Ti∈τTi. •If Ti=Tj, then ˚ Ti∩˚ Tj=∅. •The intersection of two elements is a common vertex or edge. The MATLAB mesh generator was used to obtain a partition of triangles. Each vertex of the triangles is a node, bi, and the pair (τ, Σ), where Σ = {bi}1≤i≤N, is the mesh of finite elements. Furthermore, the basis functions were constructed by using first-order finite Lagrange elements, verifying: •Each ϕiis continuous and bounded. •There is a basis function for each node. •ϕi(bj) = δij, being δthe Kronecker delta, i.e., 1 if i=jand 0 otherwise. •The restriction of each ϕito the element Tiis a first-order polynomial. 2.2.5 Modeling of metabolic oxygen consumption The oxygen consumption rate due to cellular metabolic activity in equation (2.8) is usually modeled following the Michaelis-Menten kinetics form [93]: gmet(p(t, x)) = gmax p(t, x) k+p(t, x),(2.9) where gmax is the maximum oxygen metabolic cosumption rate, and kis the oxygen pressure for half-maximum metabolic consumption rate. Typically, parameter values are set to gmax = 15 mmHg and k= 2.5 mmHg, based on experimental studies [94]. For the geometry depicted in Figure 2.1, the steady-state solution of the oxygenation problem 2.1 with the consumption term given by equation (2.9) is shown in Figure 2.2. 39 Isabel Gonz´ alez Crespo Figure 2.2: Steady-state solution of the oxygenation problem 2.1 with the consumption term given by equation (2.9) for the tumor geometry shown in Figure 2.1. 2.2.6 Modeling of Radiolytic Oxygen Depletion (ROD) in conv-RT and FLASH-RT As mentioned in the previous chapter, radiation promotes oxygen depletion caused by the radiolysis of water molecules. Many studies have investigated oxygen depletion patterns both with conv-RT and FLASH-RT in different solutions simulating cellular compounds [99–101]. These works suggest that the depleted oxygen depends on the initial oxygen partial pressure and shows a pronounced curvature in the depletion pattern as the oxygen pressure approaches zero. Besides, Cao et al. [78] conducted in vivo experiments and found that oxygen depletion does not depend on the radiation dose rate at UHDRs within the 50–300 Gy s−1range, but only on the radiation dose. In this context, previous works [84,100] have proposed a Michaelis-Menten kinetics to describe ROD as: gROD(p(t, x)) = D T G0p(t, x) kROD +p(t, x),(2.10) where D(Gy) is the radiation dose, T(s) the total time of irradiation, G0(mmHg Gy−1) the radiolytic consumption rate, and kROD (mmHg) the oxygen pressure for half-maximum oxygen depletion. As mentioned in the previous chapter, ROD has not been observed in in vivo experiments with conv-RT, as the tumor vascular system continuously resupplies oxygen, preventing measurable oxygen depletion at conventional dose rates. 40 Chapter 2. Methods and materials The left panel of Figure 2.3 shows oxygen depletion curves for both conventional (0.1 Gy s−1) and UHDRs (100 Gy s−1) at different oxygen statuses obtained by solving Problem 2.1 with metabolic and radiolytic consumption terms, respectively given by equations (2.9) and (2.10). As discussed, oxygen remains constant during conv-RT, while it is progressively eliminated in FLASH-RT. Moreover, the right panel of Figure 2.3 shows the transition of the oxygenation histogram belonging to Figure 2.1 from the beginning to the end of delivering 20 Gy with FLASH-RT (100 Gy s−1). 0 20 40 60 80 100 0 2 4 6 8 10 Delivered dose (%) p(mmHg) conv-RT FLASH-RT 0 10 20 30 40 0 0.02 0.04 0.06 0.08 0.1 p(mmHg) Frecuency t=t0 t=T Figure 2.3: Oxygen depletion during conventional RT (conv-RT) (0.1 Gy s−1) and FLASH-RT (100 Gy s−1) for a delivered dose of 20 Gy and the geometry depicted at Figure 2.1. The left panel shows the depletion curves for both modalities at different initial oxygen partial pressures, p. The right panel shows the change in oxygenation histograms from the beginning (t=t0) to the end (t=T) of FLASH-RT. These results were obtained by using the parameter values given in previous sections, along with G0= 0.25 mmHg Gy−1and kROD = 1 mmHg. In this context, the oxygenation problem 2.1 in FLASH-RT is as follows: Problem 2.3. Find p: [t0, tf]×Ω→[0, pcap]such that:                  ∂p ∂t (t, x) = DO2∆p(t, x)−gmet(p(t, x)) −gROD(p(t, x)),∀(t, x)∈[t0, tf]×Ω, p(t0,x) = p0(x),∀x∈Ω, p(t, x) = pcap,∀(t, x)∈[t0, tf]×Γcap, ∂p ∂n(t, x) = 0,∀(t, x)∈[t0, tf]×Γout, where gmet(p)and gROD(p)are given by equations (2.9) and (2.10), respectively. 41 Isabel Gonz´ alez Crespo Problem 2.3 is a particular case of a semilinear parabolic problem with mixed boundary conditions. The existence of a unique solution for these types of systems has been studied in the classical literature [102]. Because of the reasons mentioned before, in conv-RT, the oxygenation problem can be simplified to its steady-state form, ignoring the effect of ROD: Problem 2.4. Find p: Ω →[0, pcap]such that:          DO2∆p(x) = gmet(p(x)),∀x∈Ω, p(x) = pcap,∀x∈Γcap, ∂p ∂n(x)=0,∀x∈Γout, where gmet(p)is given by equation (2.9). 2.2.7 Implementation details Building on the work of D´ıaz-Botana [92] in conv-RT, a MATLAB code was developed to calculate the numerical solution of the oxygenation problems 2.3 and 2.4. Particularly, the original code was updated to remove deprecated functions from previous MATLAB versions (such as initmesh and refinemesh to generate the mesh, and pdenonlin to solve the problem) and adapt it to the R2022b version, and expanded to include the effect of ROD in FLASH-RT. Firstly, the geometry of the problem is created through a user-defined function createGeometry that takes as input: the VF, the side length of the squared domain, and the parameters given in Table 2.1 to generate the capillaries. This function follows the steps given in Section 2.2.2 and gives as output the decomposed geometry matrix and a vector with the radius of the capillaries. The decomposed geometry matrix contains a representation of the geometry in terms of disjointed minimal regions constructed by the decsg MATLAB algorithm. In this case, each column in the matrix belongs to a line or circle edge segment. Then, the MATLAB function createpde is used to create a model object that contains information about the PDE which defines the problem: number of equations, geometry, mesh, and boundary conditions. The generated geometry is associated with a one-equation model object through the MATLAB function geometryFromEdges, which requires as inputs the model object and the decomposed geometry matrix. Subsequently, the MATLAB function generateMesh creates the mesh to apply the FEM. In this thesis, the minimum and maximum length of mesh edges were constrained within the range 1–10 µm. Besides, the geometric order of the finite elements was specified to be linear, i.e., the nodes are the vertex of the triangles. Further analysis of the mesh suitableness to solve the oxygenation problem is presented in Appendix A. Listing 2.1 shows the code to perform the above process, where geom is 42 Chapter 2. Methods and materials a vector containing the parameter values to create de geometry, previously detailed; dl is the decomposed geometry matrix; and Ris the radius vector. 1%Create geometry 2[dl, R] = createGeometry(geom); 3%Create PDE model 4model = createpde; 5%Assign geometry to model 6geometryFromEdges(model,dl); 7%Create mesh 8generateMesh(model,'Hmin',1,'Hmax',10,'GeometricOrder','linear'); Listing 2.1: Code to create the geometry of the problem, the PDE model object, and the mesh. Once the geometry and the mesh are created, the boundary and initial conditions are defined in the PDE model object. Regarding the boundary conditions, the minimal regions in the decomposed geometry matrix must be assigned to the outer or the capillary boundary. The corresponding columns in dl are respectively stored in the vectors id out and id cap. Zero-flux Neumann conditions are imposed in line edges and Dirichlet conditions with p= 40 mmHg are set on circle segments by using the MATLAB function applyBoundaryCondition, as shown in Listing 2.2. 1%Set Neumann conditions on the outer boundary 2applyBoundaryCondition(model,'neumann','Edge',id_out,'g',0,'q',0); 3%Set Dirichlet conditions on the capillary walls 4applyBoundaryCondition(model,'dirichlet','Edge',id_cap,'u',40); Listing 2.2: Code to set the boundary conditions to solve Problems 2.3 and 2.4. To solve Problem 2.3, the initial condition is set to p0(t0) = 0through the MATLAB function setInitialConditions, as shown in Listing 2.3. 1%Set initial condition 2setInitialConditions(model, 0); Listing 2.3: Code to set the initial condition to solve Problem 2.3. Finally, the MATLAB functions specifycoefficients and solvepde are used to set the PDE coefficients and solve Problems 2.3 and 2.4, as respectively shown in Listings 2.4 and 2.5. 1%Specify coefficients: m*d^2p/d^2t + d*dp/dt - c*Lap(p) + a*p = f 2aa = @(loc,st) acoeffun(loc,st,D,G0,T,k_ROD,g_max,k); 3specifyCoefficients(model,'Face',1,'m',0,'d',1,'c',D,'a',aa,'f',0); 4sol = solvepde(model,tlist); 43 Isabel Gonz´ alez Crespo 3.1 Overview of the problem Cancer ImmunoTherapy (IT) is a therapeutic strategy that aims at boosting and exploiting the natural immune response to control and cure tumors. Particularly, Immune Checkpoint Inhibitors (ICIs), such as αCTLA4 and αPD(L)1 drugs (see Section 1.4.1), are used to treat several cancers, including melanoma, prostate, non-small cell lung cancer, and leukemia [105,106]. ICIs block different proteins, like the CTLA-4 receptor and those involved in the PD-1/PD-L1 pathway, which are well-known suppressors of the immune response against tumors. These ITs have shown promising results in preclinical experiments, and there are several monoclonal antibodies against CTLA-4 and PD-(L)1 approved to treat different types of cancer [107]. However, the efficacy of IT as a cancer treatment by itself is still limited, except for particular cases. For example, ipilimumab (αCTLA4) has shown an improvement in the survival of melanoma patients, but with low response rates, within the 10–15% range [108]. Many preclinical studies have shown that RadioImmunoTherapy (RIT) combined treatments, in particular with inhibitors of CTLA-4 and PD-(L)1, are significantly more effective than RT and IT as monotherapies [7,18,29,59,109]. Although the dominant mechanism behind the effect of RT is the direct damage on tumor cells’ DNA [110], there is evidence that radiation can trigger other indirect cell-death mechanisms, particularly at high-doses-per-fraction treatments, which may influence the synergy with IT. High radiation doses may promote an immune response against surviving tumor cells [28,29,111], increasing tumor control. These immune effects seem to be related to an increased release of tumor antigens due to radiation-mediated cell death, which enhances the activation of antitumor T-cells; and modifications of the Tumor MicroEnvironment (TME), such as the elimination of immune down-regulators like regulatory T-cells (Tregs) and Myeloid-Derived Suppressor Cells (MDSCs), facilitating T-cell infiltration in the tumor [7,50,112]. Despite the demonstrated effectiveness of RIT, how the best combination of both therapies can be achieved is still a matter of study. In this regard, validated tumor response models to RIT constitute a powerful tool to investigate the biological mechanisms behind the effect of the treatment, and may also assist in the design of optimal therapeutic strategies, potentially guiding clinicians in the selection of optimal treatments. Modeling the response of tumors to IT has been addressed in the biomathematical literature, following both phenomenological and mechanistic approaches [113– 119]. On the contrary, modeling the synergistic effect of RIT has been less studied due to the novelty of this treatment strategy, but it has become an active field of research over the last few years [34,35,120–123]. For example, Serre et al. [34,120] proposed a simple response model to RIT with αCTLA4 and αPD(L)1 ICIs that was fitted to experimental data of tumor response and rejection probability of implanted tumors. Besides, Kosinsky et al. [121] presented a more complex model, based on ODEs, which was fitted to experimental data and used to formulate hypotheses of optimal treatment strategies with RT and inhibitors of PD-L1. 50 Chapter 3. Results on radioimmunotherapy 3.2 Model of tumor response to RIT This section presents a mathematical model of tumor response to RIT with ICIs against the CTLA-4 receptor and the PD-1/PD-L1 pathway. The model has been tested by fitting it to preclinical data of tumor response to combined therapies of radiation with different fractionations along with αPDL1 or αCTLA4 [18,59], including volume dynamics and tumor control rates. The dynamics of the most relevant cell populations are described by ODEs,DDEs and IDEs. On the one hand, DDEs explicitly include biological response delays on cell dynamics, which may have importance when modeling biological systems. On the other hand, IDEs describe the effect of the treatment at the delivery times, introducing changes in cell populations through impulses. The model is presented as a continuous deterministic model to describe tumor volumes. However, when modeling TCP, it is converted into a discrete stochastic model by introducing a birth-death Markov chain. This allows the model to yield TCPs according to the clonogenic cells hypothesis, introduced in Section 1.3.4. 3.2.1 Overview of the model The proposed biomathematical model follows the cancer-immunity cycle, described in Section 1.4.1. In summary, tumor cells release antigens, both naturally and due to cell death; antigens are taken by APCs to activation sites (lymph nodes) where they trigger T-cell activation; active T-cells migrate to the tumor, infiltrate it, and eliminate tumor cells. Each RT dose fraction damages and eliminates tumor cells (contributing to the release of antigens), as well as T-cells present in the tumor. IT drugs affect either the activation of T-cells (αCTLA4) or the immune-mediated death of tumor cells by activated T-cells (αPD(L)1). A diagram of the variables of the model and the interaction between them is shown in Figure 3.1. The model includes three main species: tumor cells, T-cells, and antigens/APCs. These species are split into six compartments, which exist in two different spatial locations: in the tumor (viable cancer cells, C; tumor cells that are doomed by radiation and will eventually die, Cd; and active T-cells against tumor cells, Ta), and in the activation sites (a phenomenological compartment accounting for the effects of antigens and APCs,ˆ A; the T-cell pool, ˆ T; activated T-cells, ˆ Ta; and blocked T-cells through the CTLA-4 receptor, ˆ Tb). The compartments at the activation sites are identified with a hat ( ˆ ) to facilitate the understanding of the model. Migration between different spatial locations (antigens/APCs migrating from the tumor to activation sites and activated T-cells migrating from the activation sites and infiltrating the tumor) results in biological delays which are explicitly included in the model through DDEs (see Section 2.1.2). The effect of radiation delivery on tumor cells and irradiated T-cells, as well as the administration of IT drugs, are modeled as impulses by using IDEs (see Section 2.1.1). 51 Isabel Gonz´ alez Crespo CTa(τT) Cd RT ˆ Ta ˆ T ˆ A(τA)ˆ Tb αPD(L)1 αPD(L)1 αCTLA4 Tumor Activation site Figure 3.1: Flowchart of the model showing the interaction between different compartments: viable tumor cells (C), doomed tumor cells (Cd), active T-cells in the tumor (Ta), antigens/Antigen Presenting Cells (APCs) (ˆ A), the pool of T-cells ( ˆ T), activated T-cells ( ˆ Ta) and blocked T-cells ( ˆ Tb). The dotted line shows the separation between compartments physically located in the tumor (left) and in the activation sites (right). The notation also distinguishes between the two locations, identifying the compartments within the activation zone with a hat (ˆ). τAand τTare biological delays related to the migration of antigens and T-cells between those two physical locations. The action points of RadioTherapy (RT),anti-Programmed Death-(Ligand) 1 (αPD(L)1), and anti-CTLA-4 (αCTLA4) are also indicated. The model dynamics has the following form: dynamics C= proliferation −radiation death −immune death (PD-(L)1) dynamics ˆ A= natural release + RT-mediated release −natural elimination − T-cell activation (CTLA-4) dynamics Ta= activation/infiltration (CTLA-4)−radiation death − immune death −natural elimination The presented compartmental model attempts to describe the dynamics of the main agents involved in the tumor response to RIT. However, many simplifications were made to present a tractable problem. On the one hand, spatial coordinates were not considered, although a simple spatial dependency is indirectly introduced through biological delays simulating the migration processes between the tumor site and the T-cell activation zone. On the other hand, the process of T-cell infiltration within the tumor was not included in the model but integrated into the migration delay. Besides, immune-mediated death is overly simplified, as other cell types that participate in this process, either favoring immunity or acting as suppressors, like natural killers or Tregs, were not included. Regarding the effect of RT, the LQ model was used to account for radiation cell death, but departures from the LQ model were also studied and discussed. 52 Chapter 3. Results on radioimmunotherapy 3.2.2 Radiation cell death and tumor proliferation Dose delivery is modeled as instantaneous, which seems a good approximation, as in typical fractionated treatments it takes minutes, while the typical times of the model dynamics are days. The LQ model (1.1), described in Section 1.3.1, is typically used to model radiation-mediated cell death through the calculation of the SF. However, as previously mentioned, it is well-known that cell death can depart from the standard LQ model, especially at high-dose-per-fraction treatments. As the model was fitted to data from different fractionation schedules, including hypofractionated RT, other LQderived models were investigated and compared to provide the best fits to experimental data. In particular, the LQL model (1.7) (described in Section 1.3.2) and a simple ad hoc modification of the quadratic term in the LQ model [33]: β→β1 + c√d,(3.1) where c(Gy−1/2) is a free parameter and d(Gy) is the dose per fraction. Introducing the above modification in the parameter βof equation (1.1) results in the modified LQ (LQmod)model: SF = exp −αd −β1 + c√dd2.(3.2) Negative cvalues might be related to damage saturation effects, as sublethal damage repair in the LQL model, while positive values are linked to indirect cell death mechanisms contributing to cell killing. The LQ,LQL, and LQmod models were used to describe radiation-induced tumor cell death. While radiation-induced T-cell death was described by using the LQ model. Tumor cells lethally damaged by radiation (doomed) are not instantaneously removed from the system, but follow a given kinetics, generally a plateau (mitotic delay) followed by a progressive death as damaged cells enter mitosis and suffer mitotic catastrophe due to damage in their DNA. On the contrary, T-cells die in interphase within a few hours of irradiation (apoptosis), without intervening mitosis. The model considers a mitotic delay followed by an exponential death [37] in tumor cells and the immediate death of irradiated T-cells. The proliferation of viable tumor cells is modeled by the logistic formalism [124], as tumor growth might be limited by different factors, like physical constraints due to its location or the lack of necessary nutrients. Moreover, while doomed cells may carry some proliferating capacity (abortive divisions [125]), it should be limited and not contribute to the long-term cell population. Therefore, it was not included in the model. Being {ti, i = 1 . . . n}a collection of nradiation delivery times and {di, i = 1 . . . n} the delivered fractions of dose, the following IDE describes the proliferation dynamics 53 Isabel Gonz´ alez Crespo and radiation effect on viable tumor cells: dC dt (t) = λ1C(t)(1 −λ2Ctot(t)),∀t=ti, ∆C(ti)=(SFC(di)−1)C(ti),∀t=ti, (3.3) where λ1(days−1) is the exponential proliferation rate, λ2is a parameter related to the carrying capacity of the system, and Ctot(t) = C(t) + Cd(t) is the total number of tumor cells at time t. Moreover, SFCis the SF of tumor cells given by the LQ model (1.1), the LQL model (1.7) or the LQmod model (3.2). Each radiation fraction is considered to create new doomed cells, but not to interfere with the radiation kinetics of existing ones. Therefore, the compartment of doomed cells is split into nsubcompartments created by the ndelivered dose fractions, having: dCd,i dt (t) = −ϕω(¯ ti)Cd,i(t),∀t=ti, ∆Cd,i(ti) = (1 −SFC(di))C(ti),∀t=ti, Cd(t) = n X i=1 Cd,i. (3.4) Where, Cd,i(t) is the compartment of doomed cells created by the radiation dose fraction di,tithe delivery time of that fraction, and ¯ ti=t−ti. Notice that Cd,i(t) is defined as zero for t < ti. The parameter ϕ(days−1) is the death rate, and the function ωmodels the mitotic delay and progressive incorporation of damaged cells to cell death kinetics after a radiation fraction: ω(z) =        0,for z≤τd1, z−τd1 τd2−τd1 ,for τd1< z ≤τd2, 1,for z > τd2, (3.5) being τd1 (days) and τd2 (days) the initial and final instants of time describing the mitotic delay. Regarding the radiation-mediated death of T-cells within the tumor, it follows: ∆Ta(ti)=(SFT(di)−1)Ta(ti),∀t=ti,(3.6) where SFTis the SF of T-cells given by the LQ model (1.1). As previously mentioned, it is assumed that radiation-damaged T-cells die instantly, thus there are no kinetic terms associated with that process. 54 Chapter 3. Results on radioimmunotherapy 3.2.3 Antigen release and T-cell activation Antigens are considered to be released both naturally, due to the presence of tumor cells (with a rate proportional to the number of both viable and doomed tumor cells), and as a result of radiation-induced cell death (proportional to the rate of doomed cell elimination). Besides, the model includes a term that describes their natural removal from the system. A biological delay, τA(days), between antigen release and T-cell activation accounts for the time that APCs take to collect antigens and carry them to the activation sites (see Figures 1.5 and 3.1), having: dˆ A dt (t) = ρCtot(t−τA) + ψϕ n X i=1 ω(ti−τA)Cd,i(t−τA)−σˆ A(t),(3.7) where ρ(days−1) is the natural release rate, ψis the cell-death-mediated release rate, ti−τA=t−(ti−τA), and σ(days−1) is the natural elimination rate. The funtion ω is given by equation (3.5) and ϕis the death rate from equation (3.4). The activation of T-cells against tumor cells is modeled through four bilinear equations, which describe the generation of activated T-cells ( ˆ Ta) or blocked T-cells ( ˆ Tb) (through the CTLA-4 receptor) from a pool of blank T-cells ( ˆ T): dˆ A dt (t) = −aˆ A(t)ˆ T(t)−bˆ A(t)ˆ T(t),(3.8) dˆ T dt (t) = −aˆ A(t)ˆ T(t)−bˆ A(t)ˆ T(t) + h, (3.9) dˆ Ta dt (t) = aˆ A(t)ˆ T(t),(3.10) dˆ Tb dt (t) = bˆ A(t)ˆ T(t).(3.11) Parameters a(days−1) and b(days−1) describe the affinities for activation and inactivation, respectively [34]. Moreover, the amount of available blank T-cells decreases due to the activation/inactivation process, and it is renewed with a constant pool (due to maturation of new T-cells), h(days−1). To prevent the T-cell compartment from unlimited growth, an initial condition was set, ˆ T(0) = T0, which serves also as a carrying capacity, imposing ˆ T(t)≤T0,∀t. Active T-cells, ˆ Ta, migrate and infiltrate in the tumor with an associated biological delay τT(days), which models the time needed to move from the activation sites to the tumor zone and act on tumor cells. Once within the tumor, they become part of the compartment Ta, resulting in: dTa dt (t) = dˆ Ta dt (t−τT) = aˆ A(t−τT)ˆ T(t−τT).(3.12) 55 Isabel Gonz´ alez Crespo To test the hypothesis that vascular damage at high radiation doses [11] may reduce the effectiveness of RIT by limiting the infiltration of T-cells in the tumor, a simple modulation factor was introduced on active T-cells within the tumor. This time-dependent term accounts for vascular damage and progressive recovery. Inspired by previous works [11,95], critical vascular damage is considered for doses beyond 15 Gy, and the progressive recovery of the vascular function is modeled as: f(t) = min{0.05(t−t′),1},(3.13) where t−t′is the recovery time since the last radiation dose delivered at t=t′, and 0.05 days−1is the recovery rate [95]. This term represents the fraction of active T-cells reaching the tumor and multiplies the right side of equation (3.12). 3.2.4 Immnune-mediated cell death Interaction between active T-cells and tumor cells within the tumor results in the partial depletion of both. Following the work of de Pillis et al. [113], this interaction is modeled with a bilinear term for the compartment Ta, in addition to an exponential natural elimination: dTa dt (t) = −ιTa(t)Ctot(t)−ηTa(t),(3.14) where ι(days−1) is the depletion rate due to cell competition and η(days−1) is the natural elimination rate. On the other hand, immune-mediated death of viable tumor cells (analogous for Cd(t)) is modeled by the following rational term: dC dt (t) = −p(Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qC(t),(3.15) where p(days−1) is the maximum depletion rate, qrepresents how the elimination rate depends on the effector/target ratio, and smodules the steepness of the immunemediated-death curve. Data fitting experiments of cell lysis data showed that while a simple bilinear term is sufficient to describe immune-mediated T-cell death, the elimination of tumor cells due to the interaction with immune cells is better described by a rational term like that in equation (3.15). The rational term allows for a saturation effect which may result from antigen-specific T-cells targeting a specific tumor cell type [113]. 3.2.5 The effect of αPD(L)1 and αCTLA4 The biokinetics of αPD(L)1 and αCTLA4 ITs are included in the model through the variables p1and c4, which represent the concentration of the respective drugs in the system. The administration is modeled as instantaneous source terms at the injection 56 Chapter 3. Results on radioimmunotherapy times (respectively, {tk, k = 1 . . . l}) and {tj, j = 1 . . . m}). Besides, the drugs are progressively eliminated from the organism following an exponential clearance, having: dc4 dt (t) = −νc4(t),∀t=tj, ∆c4(tj) = ic4,∀t=tj, (3.16) and: dp1 dt (t) = −µp1(t),∀t=tk, ∆p1(tk) = ip1,∀t=tk. (3.17) Being ν(days−1) and µ(days−1) the clearance rates of αCTLA4 and αPD(L)1 drugs, and ic4and ip1the respective administered doses in arbitrary units. The pharmacokinetic modeling of IT drugs is simplified. However, equation (3.16) seems a good approximation for the kinetics of αCTLA4, as Selby et al. [126] investigated the biokinetics of different αCTLA4 drugs finding that they follow linear forms. Nevertheless, Deng et al. [127] investigated the biokinetics of αPDL1 and found a more complex non-linear behavior. Rather than considering a complex model to characterize the effect of αCTLA4 on the activation process of T-cells, a simple dependence on c4was introduced in the parameter bin equations (3.8), (3.9) and (3.11). Thus, following the work by Serre et al. [34], αCTLA4 counteracts T-cell blockade as: b→b 1 + c4(t).(3.18) On the other hand, the effect of αPD(L)1 on the interaction between tumor cells and T-cells is modeled by introducing a dependence on the parameter pin equation (3.15), having: p→p(1 + p1(t)),(3.19) that results in enhanced immune-mediated death of tumor cells. 3.2.6 The complete model Let {ti, i = 1 . . . n},{tj, j = 1 . . . m}, and {tk, k = 1 . . . l}be the delivery times of RT,αCTLA4, and αPD(L)1, respectively; and Ctot(t) = C(t) + Cd(t) with Cd(t) = Pn i=1 Cd,i. 57 Isabel Gonz´ alez Crespo Assembling the equations presented in the above sections, the complete model of tumor response to RIT is for all t=ti: dC dt (t) = λ1C(t) (1 −λ2Ctot(t)) −p(1 + p1(t)) (Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qC(t), dCd,i dt (t) = −ϕω(¯ ti)Cd,i(t)−p(1 + p1(t)) (Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qCd,i(t), dTa dt (t) = aˆ A(t−τT)ˆ T(t−τT)−ιTa(t)Ctot(t)−ηTa(t), dˆ T dt (t) = −aˆ A(t)ˆ T(t)−b 1 + c4(t)ˆ A(t)ˆ T(t) + h, dˆ Ta dt (t) = aˆ A(t)ˆ T(t), dˆ Tb dt (t) = b 1 + c4(t)ˆ A(t)ˆ T(t), dˆ A dt (t) = ρCtot(t−τA) + ψϕ n X i=1 ω(ti−τA)Cd,i(t−τA)−σˆ A(t) −aˆ A(t)ˆ T(t)−b 1 + c4(t)ˆ A(t)ˆ T(t), where the function ωis given by equation (3.5), and the biokinetics of αPD(L)1 and αCTLA4 drugs are described by: dc4 dt (t) = −νc4(t),∀t=tj, dp1 dt (t) = −µp1(t),∀t=tk. Moreover, the following impulses are introduced at treatment delivery times: ∆C(ti)=(SFC(di)−1)C(ti),∀t=ti, ∆Cd,i(ti) = (1 −SFC(di))C(ti),∀t=ti, ∆Ta(ti)=(SFT(di)−1)Ta(ti),∀t=ti, ∆c4(tj) = ic4,∀t=tj, ∆p1(tk) = ip1,∀t=tk, where SFC(d) is given by the LQ model (1.1), the LQL model (1.7) or the LQmod model (3.2); and SFT(d) is obtained from the LQ model. Besides, the initial conditions for t∈[−max{τA, τT},0] are C(t) = C0>0, ˆ T(t) = T0>0, and the remaining initial values are set to 0. Some notes on the formal analysis of the proposed model are presented in Appendix B. 58 Chapter 3. Results on radioimmunotherapy 3.3 Specific methods and materials in RIT In addition to the general methods described in Chapter 2, this section presents specific methods used in the RIT study. 3.3.1 Numerical solution The model presented in Section 3.2.6 was numerically solved by employing the forward Euler method, described in Section 2.1.3. To obtain the numerical solution, different user-defined functions were coded in MATLAB, which are available on the Dataverse repository [128]. Details about a convergence study can be found in Appendix C. 3.3.2 Experimental data To validate and test the model, experimental data were collected from published preclinical studies available in the literature. On the one hand, Dewan et al. [59] studied the response of preclinical tumors in mice populations to RT (different fractionations) and αCTLA4, either as monotherapies or combined. Tumor cells (TSA breast carcinoma cells) were implanted and let grow for 12 days, until reaching a volume of ∼32 mm3. Treatments started at that time, and the evolution of tumor volumes was monitored every 5 days until day 35 after tumor injection. The study also reported the fraction of animals where tumor control was achieved (no evidence of tumor at the time of euthanasia). Establishing the first day of the treatment as day 0, the authors investigated tumor response to the following RIT schedules (see Figure 3.2): i) control, ii) RT alone, 20 Gy single-fraction, iii) RT alone, 8 Gy×3 (days 0, 1, 2), iv) RT alone, 6 Gy×5 (days 0, 1, 2, 3, 4), v) IT alone, delivered in 3 fractions (days 2, 5, 8), vi) RIT, 20 Gy×1RT and 3 fractions IT (days 2, 5, 8), vii) RIT, 8 Gy×3RT and 3 fractions IT (days 2, 5, 8), viii) RIT, 6 Gy×5RT and 3 fractions IT (days 2, 5, 8), ix) RIT, 8 Gy×3RT and 3 fractions IT (days 0, 3, 6), x) RIT, 8 Gy×3RT and 3 fractions IT (days 4, 6, 8). 59 Isabel Gonz´ alez Crespo 0 5 10 15 20 0 200 400 600 800 1000 Time post-start of treatment (days) Volume (mm3) Control IT (0 3 6 9) RT 12 Gy×1 IT & RT 12 Gy×1 Figure 3.7: Model fitting to experimental data reported by Deng et al. [18] of tumor response to RadioTherapy (RT) (12 Gy single fraction on day 0), ImmunoTherapy (IT) with αPDL1, and combined treatments. The LQ model (1.1) accounts for the radiosensitivity of tumor cells. Notice that differences between control and IT curves are small and both of them overlap in the figure. (αC= 0.0299 Gy−1) were allowed to vary. Because this dataset only includes one dose per fraction, the LQ model (1.1) was employed to describe radiation damage on tumor cells, imposing αC/βC= 10 Gy and considering only αCas a free parameter. While there are differences in the clones and tumors that could justify using different hostrelated and tumor-related parameters, imposing such constraints on the optimization limits reaching good fits by over-fitting. 3.4.2 Study of vascular damage on RIT effectiveness Large radiation doses can seriously damage tumor vasculature, which might limit the infiltration of active T-cells in the tumor. This might explain the poorer results obtained with the 20 Gy single-fraction irradiation by Dewan et al. [59]. To test the hypothesis that vascular damage may affect the effectiveness of RIT, a dose and time-dependent T-cell infiltrating parameter was introduced in the model through equation (3.13), to account for vascular damage and recovery. This term represents the fraction of active T-cells reaching the tumor. Inspired by previous works [11,95], critical vascular damage was considered for radiation doses above 15 Gy, followed by a progressive recovery of the vascular function. Figure 3.8 shows the best fits of the model to data presented by Dewan et al. [59], considering vascular damage and the LQ model (1.1) to account for the radiationmediated damage on tumor cells. Best-fitting parameter values are reported in Table D.4, including relevant parameters related to proliferation (λ1= 0.1357 days−1), radiation damage (αC= 0.0200 Gy−1,βC= 0.0022 Gy−2) and the immune effect 66 Chapter 3. Results on radioimmunotherapy on tumor cells (p= 24.9382 days−1) and T-cells (ι= 6.2472 ×10−9days−1). For reference, the best-fitting value of the cost function was F= 25.3345. The AIC was used to compare the goodness-of-fit between the three versions of the response model used to fit the data reported by Dewan et al. [59]. The likelihood of the fit, L, in equation (2.13) was computed as the product of the individual probabilities of having the modeled volumes assuming a Gaussian distribution. For n= 50 datapoints, the LQ model with vascular damage (AIC = 436.24, k= 19) was preferred to the LQmod model (AIC = 446.02, k= 20) and the LQL model (AIC = 458.75, k= 20) according to their AIC. 0 5 10 15 20 0 200 400 600 Volume (mm3) Control IT (2 5 8) RT 20 Gy×1 IT (2 5 8) & RT 20 Gy×1 0 5 10 15 20 0 200 400 600 Time post-start of treatment (days) Volume (mm3) Control RT 8 Gy×3 IT (2 5 8) & RT 8 Gy×3 0 5 10 15 20 0 200 400 600 Control RT 6 Gy×5 IT (2 5 8) & RT 6 Gy×5 0 5 10 15 20 0 200 400 600 Time post-start of treatment (days) Control IT (0 3 6) & RT 8 Gy×3 IT (4 6 9) & RT 8 Gy×3 Figure 3.8: Model fitting to experimental data reported by Dewan et al. [59] of tumor response to RadioTherapy (RT),ImmunoTherapy (IT) with αCTLA4, and combined treatments. The LQ model characterizes tumor cell response to radiation, and vascular damage at 20 Gy per fraction was included to limit T-cell infiltration in the tumor. Radiation doses were delivered on consecutive days starting from day 0. Notice that differences between control and IT curves are small and both of them overlap in the figure. All the fits were obtained with a single set of parameters, although they are plotted separately to facilitate the visualization. The control curve is included in all the panels as a common reference. 67 Isabel Gonz´ alez Crespo Although the volumes presented in the previous figures include tumor cells and T-cells, the model should not reproduce tumor volumes by including low fractions of tumor cells and large fractions of T-cells, which would eventually lead to tumor control contradicting experimental evidence. Figures D.1 and D.2 present the contribution of tumor cells and T-cells to tumor volumes, showing that modeled tumor volumes are dominated by tumor cells. 3.4.3 Local sensitivity analysis A local sensitivity analysis was performed, as described in Section 3.3.5. Table 3.1 ranks the sensitivity of the cost function to model parameters. The model is most sensitive to parameters describing tumor cell proliferation, radiosensitivity, and immunemediated tumor cell death. Parameter Description Sensitivity index qImmune-mediated tumor cell death parameter 309.4673 λ1Proliferation rate of tumor cells (days−1) 71.6561 sSlope of the immune-mediated tumor death curve 20.7941 pImmune-mediated tumor cell death rate (days−1) 12.8123 T0Initial pool of blank T-cells in activation site 7.3442 ηRate of natural elimination of T-cells (days−1) 4.2115 hRate of production of blank T-cells (days−1) 3.5890 αCLinear parameter of LQ model for tumor cells (Gy−1) 2.1120 βCQuadratic parameter of LQ model for tumor cells (Gy−2) 1.8827 τADelay between antigen liberation and T-cell activation (days) 0.6541 Table 3.1: Analysis of the most critical model parameters. The sensitivity index given by equation (3.21) for each parameter is calculated as the difference between the cost function of fits to experimental data presented by Dewan et al. [59] (best-fitting parameters reported in Table D.4) and the cost function obtained when a 10% perturbation is applied to that particular parameter. Only the ten most critical parameters are shown. For reference, the best cost value was F= 25.3345. 3.4.4 Optimization of IT administration on RIT treatments Finding the optimal sequence of administration of RIT can improve effectiveness, as experimentally shown in results presented by Dewan et al. [59]. In that study, the authors found that the combination of 8 Gy×3RT fractions and 3 fractions of αCTLA4 lead to different tumor responses depending on the days of administration of the IT (RT was always delivered on the same days). In particular, different IT schedules lead to different tumor volumes and control rates. This may be caused by the interplay between biological mechanisms with different kinetics, such as the biological delays between the release of antigens and the activation of T-cells, as well as the migration of such T-cells to the tumor, tumor proliferation, and the progressive death of doomed cells. 68 Chapter 3. Results on radioimmunotherapy The presented model along with the best-fitting parameters reported in Table D.4 were used to investigate optimal research strategies by simulating TCP, following the procedure described in Section 3.3.6. The administration of the three 8 Gy RT fractions was fixed on days 0, 1, and 2, while the delivery schedule of the three IT doses was varied, starting from day 0 to day 7, and ending from day 2 to day 9. The evolution of tumor volumes and the TCP were evaluated from 200 simulations on each combination. 0 10 20 0 50 100 150 (a) Volume (mm3) 0 10 20 0 50 100 150 (b) Time post-start of treatment (days) 0 10 20 0 50 100 150 (c) (0 3 6) (2 5 8) (4 6 9) 0 0.2 0.4 0.6 0.8 1 (d) 3/6 13/16 1/6 112/200 122/200 24/200 Days of αCTLA4 administration TCP 0123 0 0.2 0.4 0.6 0.8 1 (e) αCTLA4 concentration on day 2 (relative units) 0123 0 0.2 0.4 0.6 0.8 1 Figure 3.9: Study of optimal schedules of RT (8 Gy×3 fractions) and IT (3 fractions of αCTLA4) obtained from the presented model and best-fitting parameters reported in Table D.4.RT fractions are delivered at days (0, 1, 2), and IT is delivered with different schedules starting from day 0 to day 7. Panels (a), (b), and (c) report the dynamics of tumor volumes and the 95% confidence intervals for the combinations investigated by Dewan et al. [59], delivering IT at days (0, 3, 6), (2, 5, 8), and (4, 6, 9) respectively. Panel (d) shows the comparative between 95% confidence intervals for Tumor Control Probability (TCP) values obtained with the model (200 simulations) and experimental controls (6 to 16 animals), for the same treatment combinations. Panel (e) presents TCP (asterisks) versus the concentration of αCTLA4 on day 2 after starting treatment, showing a positive correlation between those two variables. The solid line corresponds to the fit of a logistic function. 69 Isabel Gonz´ alez Crespo Results are reported in Figure 3.9. Panels (a)–(c) report the evolution of tumor volumes with best-fitting parameters and 95% confidence intervals (obtained from the 200 different simulations) for the combinations investigated by Dewan et al. [59]. Panel (d) shows the comparative of 95% confidence intervals for TCP values obtained from the model (200 simulations) and experimental controls (6 to 16 animals), for the same previous combinations. The combination leading to the best tumor response is that delivering IT at days 2, 5 and 8. Moreover, the relation between modeled TCP and different metrics related to the administration of IT was studied. As shown in panel (e), the concentration of αCTLA4 on day 2 is associated with TCP. 0 10 20 30 0 50 100 150 200 (a) Time post-start of treatment (days) Volume (mm3) RT 0 10 20 30 0 50 100 150 200 (b) Time post-start of treatment (days) RIT 16.25 Gy×1 10.49 Gy×2 8 Gy×3 5.58 Gy×5 3.30 Gy×10 1 2 3 5 10 0 0.2 0.4 0.6 0.8 1 (c) Number of RT fractions TCP 1 2 3 5 10 0 0.2 0.4 0.6 0.8 1 (d) Number of RT fractions Figure 3.10: Modeled responses to different RadioTherapy (RT) schedules with αCTLA4 (2, 5, 8). RT fractionations are equivalent from a classical radiobiological point of view (same BED). Tumor volume evolution is shown for RT alone (a), and RadioImmunoTherapy (RIT) with αCTLA4 (b). Tumor Control Probability (TCP) (95% confidence intervals) were obtained from the model (200 simulations) for RT alone (c), and RIT (d). The model parameters used for this study are presented in Table D.4, and include vascular damage effect at 16.25 Gy, governed by equation (3.13). 70 Chapter 3. Results on radioimmunotherapy 3.4.5 Optimization of dose fractionation on RIT treatments The effect of different RT fractionations as monotherapy or in combination with αCTLA4 (fixed to days 2, 5, and 8) was studied. The treatment with 8 Gy×3 fractions was used as a reference, and different RT fractionations (1, 2, 3, 5, or 10 fractions, one fraction per day on consecutive days) were studied. The doses of each fractionation schedule were selected to have the same BED, given by equation (1.11), using αC/βC= 9.29 Gy (corresponding to the value presented on Table D.4). For this αC/βCratio, the resulting doses per fraction were: 16.25 (×1), 10.49 (×2), 8 (×3), 5.58 (×5), and 3.30 (×10) Gy. Figure 3.10 reports the evolution of tumor volumes under different treatment regimens, as well as TCPs (replicating the experiment described in the previous section for 200 simulations and normally distributed perturbations on the reference values from Table D.4). Conventional fractionation may be sub-standard, as daily fractions can deplete active T-cells from the tumor (it is widely assumed that T-cells are radiosensitive, even though this idea may be contradicted by recent evidence [133]). In this regard, hypofractionated schedules may prove more effective, as already seen in some experimental studies, but without reaching extreme hypofractionation. In the latter case, the strategy may result disadvantageous due to two factors: on the one hand, a single fraction may fail to keep therapeutic numbers of T-cells in the tumor for a long time; on the other hand, large doses may lose some effectiveness, and can seriously damage tumor vasculature, which might limit the infiltration of active T-cells in the tumor. Figure 3.11 shows the evolution of the main populations for three different treatments. As can be seen in panels (a) and (b), in line with what was hypothesized above, the combined treatment with hypofractionated non-single-dose radiation is the most effective. Furthermore, in panels (c) and (d) it can be seen that while with single-dose treatment the number of T-cells reaching the tumor is limited, with a more fractionated treatment the radiation eliminates active T-cells in the tumor reducing its efficacy. Table 3.2 summarizes the findings and shows the strengths and weaknesses of each treatment. Activation Infiltration Death Extreme-hypofractionated RT ✓✘✓ Moderate-hypofractionated RT ✓ ✓ ✓ Conventional RT ✓ ✓ ✘ Table 3.2: Comparative between different fractionated RadioTherapy (RT) schedules regarding T-cell’s activation, infiltration, and radiation-mediated death. While extreme-hypofractionated regimens may compromise the therapeutic effect by reducing T-cell infiltration, conventional treatments might be also suboptimal by promoting the elimination of T-cells due to multiple irradiations. 71 Isabel Gonz´ alez Crespo 0 10 20 30 0 1 2 3 4 5×108 (a) Tumor cells RT 16.25 Gy×1 10.49 Gy×2 3.30 Gy×10 0 10 20 30 0 2 4 6 8×105 (c) Time post-start of treatment (days) T-cells in tumor 0 10 20 30 0 1 2 3 4 5×108 (b) RIT 0 10 20 30 0 2 4 6 8×105 (d) Time post-start of treatment (days) Figure 3.11: Evolution of tumor cells and active T-cells in the tumor zone for the different fractionations: single-dose, non-single-dose hypofractionation, and more fractionated treatment. The panels on the left correspond to RadioTherapy (RT) and those on the right to RadioImmunoTherapy (RIT) with αCTLA4 (2 5 8). Non-single-dose hypofractionated treatments appear to be more effective, as they allow T-cell infiltration into the tumor and reduce T-cell damage. 3.5 Discussion The proposed model fits experimental volume curves of tumor response to different combinations of RIT with αPDL1 and αCTLA4, as shown in Sections 3.4.1 and 3.4.2. Furthermore, it reproduces experimental TCP values, as presented in Section 3.4.4, even though these fits are favored by large experimental uncertainties obtained from trials with small population sample sizes. However, a limitation of the presented model is the risk of overfitting due to the lack of validation data compared to the number of parameters involved. To mitigate this weakness, as many parameters as possible were set to fixed values according to the existing literature, as shown in Section 3.3.4, 72 Chapter 3. Results on radioimmunotherapy and their range of variation in the optimization process was limited to consistent values, reported in Table D.1. Besides, most of the parameters (16/20) obtained from the fitting to data of RIT with αCTLA4 were kept fixed when simulating data of RIT with αPDL1, as mentioned in Section 3.4.1. These constraints strengthen the obtained results, although they should still be taken with care. On the other hand, the model fails to reproduce the effect of αPDL1 and αCTLA4 as monotherapies on the progression of tumor volumes. While this effect is small in the collected experimental data, they show a consistent influence of IT on tumor volume. However, the model shows a non-noticeable effect of αCTLA4 and αPDL1 monotherapies for the considered experiments, as reported in Sections 3.4.1 and 3.4.2. Therefore, its application to the particular case of IT as monotherapy could be debatable. In any case, it might be that the low percentage of subjects responding to IT monotherapy treatments present a particular phenotype (model parameters) that makes them more sensitive to such therapies, an effect that cannot be reproduced when fitting population-averaged data. The fits suggest very radioresistant tumors (low αvalues) but are in line with reported experimental values for tumors in mice [134]. Moreover, the best fits were obtained when the relative radiosensitivity of tumor cells decreases with increasing dose per fraction, as predicted by the LQL model. Such effect could also be explained by a limited infiltration of T-cells into the tumor due to vascular damage at large doses per fraction. It must be noticed that, for the sake of simplicity, only the effect of limited T-cell infiltration was included in the analysis, but vascular damage may cause other effects, such as proliferation arrest, starvation, and hypoxia, resulting in a complex interaction. The obtained results led to the formulation of some hypotheses regarding the effectiveness of RIT and optimal treatment combinations. On the one hand, experimental results show that delivering IT at days (2, 5, 8) leads to a better response than deliveries at days (0, 3, 6) or (4, 6, 9). Given the slow clearance rate of αCTLA4, with mean lifetime T1/2∼7 days, the elimination from day 0 to day 2, and from day 2 to day 4, is ∼20% of the drug concentration. This seems to point out that biological delays may influence the activation and/or immune effect of T-cells. The modeling study of optimal combinations (see Section 3.4.4) suggests that a better synergy could be obtained by delivering IT soon after the first fraction of RT, and then continuing IT fractions longer into the treatment, to keep therapeutic numbers of T-cells in the tumor for a longer time (according to the different kinetics of the IT drug, tumor cell death, and antigen liberation). The combined treatment starts to lose effectiveness if the administration of IT is delayed beyond 2–3 days post-start of RT. On the other hand, both conventional fractionation and extreme hypofractionation may be suboptimal in RIT treatments, as reported in Section 3.4.5. High doses per fraction used in extreme hypofractionated schedules may lead to loss of effectiveness due to saturation damage (modeled by the LQL model), and cause critical vascular damage compromising T-cell infiltration. This might also explain the poorer results obtained 73 Isabel Gonz´ alez Crespo with the 20 Gy single-fraction irradiation. On the contrary, conventional fractionations cause multiple irradiations of circulating T-cells, reducing their immune activity. Thus, moderate hypofractionation of the radiation dose may offer the best results. Direct extrapolation of the above hypotheses to clinical data must be considered with care, for they are based on the analysis of preclinical data. In human tumors, many biological processes are slower than in mice (such as cell proliferation, clearance and biokinetics of drugs, and migration) and will certainly affect the complex interplay that leads to optimal treatment combinations. The combined study of indirect tumor cell death mechanisms (due to vascular damage and radiation-triggered immune response) and direct radiation damage on tumor cells may be paramount in RT and RIT, as both processes may interfere with each other. Biomathematical models linking these effects can provide insight into the problem and help to interpret experimental results, as well as to interpret conflicting reports on the effect of indirect cell death on tumor response [11,21]. Compared to other models of IT and RIT response, this thesis presents a simpler modeling of the immune-mediated tumor cell death, including a single population of T-cells and ignoring other populations. Besides, the most relevant novelty is the introduction of biological delays (through DDEs). It is hypothesized that these biological delays play an important role in the response to RIT, and should be considered in the design of optimal combination strategies. This seems to be supported by experimental results, as argued above. Moreover, a stochastic approach was introduced to compute TCPs, and the effect of vascular damage in the response to RIT was investigated. 74 Chapter 4 Results on FLASH radiotherapy “I’m a racecar, you’re. . .a much older racecar, but under the hood, you and I are the same.” “We are not the same! Understand?” — Lightning McQueen & Doc Hudson, Cars. This chapter summarizes the outcomes achieved during the thesis on developing mathematical models for FLASH radiotherapy. This novel radiation modality has preclinically proven to decrease toxicity in healthy tissue around tumors while maintaining the anti-tumor effect of conventional radiotherapy. Building on previous works, this chapter presents a mathematical approach to replicate tumor volume evolution curves in preclinical FLASH and conventional radiotherapy studies. This allows for the examination of the mechanisms behind the effect of FLASH radiotherapy and the analysis of its presumed isoeffectiveness compared to conventional modalities. This chapter partially reproduces the content of the article: I. Gonz´alezCrespo, F. G´omez, ´ O. L´opez Pouso, and J. Pardo-Montero, “An in-silico study of conventional and FLASH radiotherapy iso-effectiveness: potential impact of radiolytic oxygen depletion on tumor growth curves and tumor control probability,” Physics in Medicine & Biology, vol. 69, p. 215016, 2024. doi:10.1088/1361-6560/ad8291. See the Copyright terms in Appendix F. Isabel Gonz´ alez Crespo Moreover, the multcompare function was used to perform multiple comparison tests and obtain the associated p-values. Two curves were considered significantly different if p-value <0.05. 4.5.6 Generation of simulated samples Following the steps given in Section 2.2.2, 100 simulated tumor geometries were generated, using random VFs sampled from a normal distribution with mean µ=0.10, and standard deviation σ= 0.04, adding a 0.04 positive cutoff to avoid extremely poor vascularizations. The steady-state problem (Problem 2.4) was solved for each of the generated tumors to obtain oxygenation profiles, having baseline mean oxygen levels in the range ∼5–30 mmHg, similar to those reported in the literature [101]. Oxygen distributions, named as the Hsample onwards, were limited to 100 simulations due to the computational cost of solving the oxygenation models. Figure E.1 illustrates the spatial oxygen distributions and the respective oxygenation histograms of two elements with different VFs,fv= 0.04 and fv= 0.11. Figure E.2 summarizes the oxygen distributions in the Hsample as a heatmap. When simulating oxygen evolution during FLASH-RT through Problem 2.3, different G0values were used for each distribution to introduce variability in oxygen depletion rates. 4.5.7 TCP estimation The Poisson-LQ model (1.14), described in Section 1.3.4, was used to calculate TCP from the SFs obtained through equations (4.4) and (4.5), respectively in conv-RT and FLASH-RT. To account for population heterogeneity on TCP calculations, a set of 1000 simulated tumors was obtained by randomly assigning different model parameters to each of them. In particular, each tumor was characterized by a radiosensitivity parameter, αox in equation (1.4), and a number of tumor cells, N0in equation (1.14), sampled from normal distributions with standard deviation equal to 20% of the mean values. Moreover, each one of the 1000 simulated tumors was associated with an oxygen distribution and radiolytic consumption rate from the Hsample, introduced in the previous section. TCP was calculated for several doses (assuming single fraction treatments) to obtain TCP-dose curves. Such curves were then characterized by calculating D50 (Gy) and D90 (Gy), the doses required to obtain TCP = 0.5 and 0.9, respectively, by interpolating the TCP-dose curve. 82 Chapter 4. Results on FLASH radiotherapy 4.6 Results This section summarizes the results on FLASH-RT regarding the use of the oxygenation model to reproduce oxygen depletion measurements from preclinical studies, the fit of the presented model of tumor evolution to preclinical volumes, and the simulation of TCP to study the assumed iso-effectiveness between conv-RT and FLASH-RT. 4.6.1 Simulation of oxygen depletion in FLASH-RT Equation (4.1) was fitted to the measurements of depleted oxygen in water and other solutions reported in Table 4.1. Using the methodology described in Section 4.5.3, best-fitting values of G0and kROD were obtained. Figure 4.1 and Table 4.4 show the fitted curves and the best-fitting parameter values with 95% confidence intervals, respectively. From this study, the radiolytic depletion rate was set to G0= 0.25 mmHg Gy−1as a reference value for the rest of the calculations, close to the mean value obtained from the fits presented in Table 4.4, and similar to values previously used by other authors [84]. The parameter kROD was set to 1 mmHg, which was also employed in previous modeling studies [84], and near to the best-fitting values reported in Table 4.4. 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 Initial p(mmHg) ∆p/D (mmHg Gy−1) i. 32.3 Gy, 201.9 Gy s−1ii. 31.9 Gy, 319 Gy s−1 iii. 32.1 Gy, 1.1×105Gy s−1iv. 30 Gy, 100 Gy s−1 v. 30 Gy, 100 Gy s−1 Figure 4.1: Oxygen depletion curves obtained by fitting equation (4.1) to the measurements reported by Jansen et al. [100] (i, ii and iii) in water, Van Slyke et al. [101] (iv) in CELL solution, and El Khatib et al. [99] (v) in BSA solution. The total amount of depleted oxygen during FLASHRT, ∆p, divided by the delivered radiation dose is presented against the initial oxygen partial pressure. 83 Isabel Gonz´ alez Crespo G0(mmHg Gy−1)kROD (mmHg) i 0.179 ±0.009 1.509 ±1.112 ii 0.202 ±0.022 2.418 ±1.547 iii 0.166 ±0.006 2.282 ±0.732 iv 0.407 ±0.019 1.303 ±0.837 v 0.365 ±0.021 0.650 ±1.120 Table 4.4: Best-fitting values of G0and kROD, and 95% confidence intervals obtained by fitting (4.1) to the oxygen depletion curves reported by El Khatib et al. [99], Jansen et al. [100], and Van Slyke et al. [101]. The fits are shown in Figure 4.1. Subsequently, the in vivo scenario was investigated. Using the simulated tumors from the Hsample described in Section 4.5.6, oxygen depletion was calculated after aFLASH-RT delivery of 30 Gy with a dose rate of 100 Gy s−1, and the results were compared with those reported by Van Slyke et al. [101]. A different radiolytic consumption rate was assigned to each simulated tumor following a normal distribution with mean µ= 0.25 mmHg Gy−1and 20% deviation from the mean value. Figure 4.2 presents the simulated and experimental data, which show a good agreement on the slope of the linear fits. 0 5 10 15 20 25 30 35 0 2 4 6 8 Initial p(mmHg) ∆p(mmHg) In vivo ∆p= (0.215 ±0.065) p Model ∆p= (0.223 ±0.069) p Figure 4.2: Amount of depleted oxygen during FLASH-RT (30 Gy, 100 Gy s−1) versus the initial mean oxygen partial pressure, ¯p, in preclinical tumors. The triangles represent in vivo data reported by Van Slyke et al. [101], and the circles represent the simulated data for 100 tumors, obtained by solving Problem 2.3. Linear fits of each dataset are presented as solid and dashed lines, respectively. 84 Chapter 4. Results on FLASH radiotherapy 4.6.2 Comparison of SF in conv-RT and FLASH-RT in vitro As mentioned in Section 4.2, the effect of ROD on the SF of irradiated cells in vitro was investigated, as a previous step to study heterogeneous oxygen distributions simulating in vivo experiments. Fixing G0= 0.25 mmHg/Gy and kROD = 1 mmHg, equation (4.1) was employed to obtain oxygen depletion curves during FLASH-RT for different radiation doses and baseline oxygen levels. On the other hand, the oxygen partial pressure remained steady during irradiation with conv-RT. Then, the SFs for conv-RT and FLASH-RT were obtained as described in Section 4.3. The study was limited to oxygenations in the range 1–40 mmHg, as it was previously reported (both from experimental [142] and modeling studies [143]) that for oxygen pressures above 30–40 mmHg the FLASH effect is not observed. Different values of the radiosensitivity parameters, αox and βox, reported in Table 4.5, were investigated. These values were set to be iso-effective for a dose of 20 Gy using the LQ model (1.1). αox/βox (Gy) αox (Gy−1)βox (Gy−2) 3 (low) 0.157 0.052 10 (medium) 0.400 0.040 20 (high) 0.600 0.030 ∞1.200 0 Table 4.5: Radiosensitivity parameters and ratios employed in the in vitro study. These values were set to be iso-effective for a dose of 20 Gy with conv-RT. 0 2 4 6 8 10 12 14 16 18 20 0 40 80 120 p(mmHg) SFF/SFC(20 Gy) αox/βox = 3 Gy αox/βox = 10 Gy αox/βox = 20 Gy αox/βox =∞Gy Figure 4.3: Ratio of Surviving Fractions (SFs) between FLASH-RT (SFF) and conv-RT (SFC) for a dose of 20 Gy versus the oxygenation status of the cells. A dose rate of 100 Gy s−1was employed to obtain SFF. 85 Isabel Gonz´ alez Crespo Figure 4.3 illustrates the difference in SF between conv-RT and FLASH-RT versus the oxygenation for a radiation dose of 20 Gy. It is observed that differences tend to zero in the limits p→ ∞ and p→0, the reason being that at large oxygen concentrations, the effect of ROD on the OERs becomes insignificant and that ROD is limited at very low oxygen concentrations. The maximum differences are observed approximately in the range 2–6 mmHg. The location of this range depends on both the modeling of ROD and the modeling of the OER, particularly in the parameter km, which was set to 3.28 mmHg, as reported in Table 4.3. 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 3 Gy Dose (Gy) SF FLASH-RT (p= 15 mmHg) conv-RT (p= 15 mmHg) 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 10 Gy Dose (Gy) 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 20 Gy Dose (Gy) 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox =∞Gy Dose (Gy) 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 3 Gy SF FLASH-RT (p= 3 mmHg) conv-RT (p= 3 mmHg) 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 10 Gy 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox = 20 Gy 0 10 20 30 10−20 10−15 10−10 10−5 100 αox/βox =∞Gy Figure 4.4: Surviving Fraction (SF) versus dose curves for cells irradiated with conv-RT and FLASH-RT at two oxygenation levels (poorly oxygenated,p= 3 mmHg, moderately-well oxygenated, p= 15 mmHg). Results are presented for αox/βox ∈ {3,10,20,∞} Gy. 86 Chapter 4. Results on FLASH radiotherapy Figure 4.4 shows the SF versus dose curves resulting from the in vitro study. The differences in SF between FLASH-RT and conv-RT match qualitatively well with previously reported experimental studies [142], showing higher SFs in FLASH-RT. Moreover, the obtained curves reveal the influence of oxygen on these differences when comparing poorly and moderately-well oxygenated statuses, respectively represented by p= 3 mmHg and p= 15 mmHg. 4.6.3 Model fitting to preclinical tumor volume curves The optimization techniques mentioned in Section 4.5.2, were used to fit the tumor response model given by equations (4.8)–(4.11) to the datasets of tumor volume dynamics after conv-RT and FLASH-RT summarized in Table 4.2. Parameters αox,λ, K,ϕand the initial volume, V0(being Cd(0) = 0), were fitted while the ratio αox/βox was set to 10 Gy. To find the best fit considering the oxygen distribution and the radiolytic consumption rate, the optimization process was performed separately for each oxygen distribution in the Hsample, using the same G0parameter values assigned for the in vivo experiment described in Section 4.6.1. Then, the best-fitting distribution and parameters were selected as those that minimize the cost function, defined in Section 4.5.2 (however, differences in the cost value between the 100 optimizations were minimal, having an increasing value of αox with decreasing oxygenation). As summarized in Table 4.2, the two experiments performed by Diffenderfer et al. [72] corresponded to the same tumor line but for different radiation doses. Thus, the same set of model parameters were used for both datasets. On the contrary, the experiments reported by Zhu et al. [140] were treated separately as belonging to different tumor lines which may present different response parameters. Additionally, three distinct initial volumes, one for each curve of control, conv-RT, and FLASH-RT, were used to emulate noticeable differences in pre-irradiation volumes observed in the experiment with Py8119 tumors. The best-fitting curves and experimental data are shown in Figure 4.5. While the best-fitting parameters for each experimental study are summarized in Table E.1. 4.6.4 Analysis of conv-RT and FLASH-RT iso-effectiveness from dose-volume curves Building on the results of the previous section, the significance of the differences between control, conv-RT and FLASH-RT tumor growth curves was investigated. The study was performed following the next steps: Step 1: A random growth curve was generated by sampling parameters from a normal distribution with the best-fitting parameters reported in Table E.1 as mean values and a 20% relative standard deviation (which qualitatively fits the 87 Isabel Gonz´ alez Crespo error bars of the experimental studies), sampling an oxygenation from the H sample, and assigning a treatment group (control, FLASH-RT,conv-RT). Step 2: Step 1 was repeated to obtain the sample size of each experimental study. Step 3: The ANCOVA methodology was used to detect significant differences between groups, as described in Section 4.5.5. Step 4: Steps 1–3 were repeated 1000 times to achieve enough statistics. Figure 4.6 illustrates one of the resulting simulations from following Steps 1 and 2 to analyze the differences between conv-RT and FLASH-RT groups. A very large percentage of the simulated experiments showed a significant difference between control and FLASH-RT or conv-RT. However, only a small fraction of the simulated experiments showed a significant difference between FLASH-RT and conv-RT, respectively 15.8% and 28.6%, as shown in Figure 4.7. The data from the experiment by Zhu et al. [140] in Py8119 tumors ((c) in Figure 4.5) was excluded from this analysis because the baseline differences in volume between groups lead to significant differences that cannot be attributed to the treatments. 0 10 20 30 0 500 1000 1500 2000 2500 (a) Volume (mm3) Control FLASH-RT conv-RT 0 10 20 30 0 500 1000 1500 2000 2500 (b) 0 5 10 15 0 500 1000 1500 2000 2500 3000 (c) Time post-RT (days) Volume (mm3) 0 10 20 30 0 100 200 300 400 500 (d) Time post-RT (days) Figure 4.5: Best fits of tumor growth curves (mean values and standard deviations) for control, FLASH-RT, and conv-RT groups, reported by Diffenderfer et al. [72] ((a) MH641905 12 Gy, (b) MH641905 18 Gy) and Zhu et al. [140] ((c) Py8119 9.5 Gy, (d) Py230 9.5 Gy). 88 Chapter 4. Results on FLASH radiotherapy 0 10 20 30 0 500 1000 1500 (a) Volume (mm3) FLASH-RT conv-RT 0 10 20 30 0 500 1000 1500 (b) 0 5 10 15 0 500 1000 1500 2000 (c) Time post-RT (days) Volume (mm3) 0 10 20 30 0 20 40 60 80 100 (d) Time post-RT (days) Figure 4.6: Example of the comparison between tumor growth curves after irradiation with FLASH-RT and conv-RT from fits to the data reported by Diffenderfer et al. [72] ((a) MH641905 12 Gy, (b) MH641905 18 Gy) and Zhu et al. [140] ((c) Py8119 9.5 Gy, (d) Py230 9.5 Gy). The solid lines represent the mean values and the shadow areas represent the standard deviation of populations with the sample size given in Table 4.2 for each experiment. (a) (b) (d) 0 20 40 60 80 100 98.8 100.0 100.0 98.1 100.0 100.0 23.0 28.6 15.8 Control vs. conv-RT Control vs. FLASH-RT conv-RT vs. FLASH-RT Figure 4.7: Percentage of simulated experiments showing a significant difference (p-value <0.05) between control, conv-RT, and FLASH-RT groups, for simulations based on the experimental data reported by Diffenderfer et al. [72] and Zhu et al. [140], summarized in Table 4.2. 89 Isabel Gonz´ alez Crespo 4.6.5 Estimation of the effect of ROD on TCP Once the effect of ROD on the SF of tumor cells in vitro was estimated, the obtained results were extrapolated to TCP versus dose curves in vivo. For that purpose, TCP was computed on simulated tumors with heterogeneous oxygen distributions, following the methodology described in Section 4.5.7. Different αox/βox values were studied and reported in Table 4.6. Fixing αox = 0.4 Gy−1and αox/βox = 10 Gy as reference values, the other radiosensitivity parameters, αox, were set to yield the same D50 with convRT. Populations of 1000 simulated tumors for each group (radiosensitivity ratio) were generated using the data on Table 4.6 and N0= 106tumor cells as mean values for parameter perturbations. αox/βox (Gy) αox (Gy−1)βox (Gy−2) 3 (low) 0.175 0.058 10 (medium) 0.400 0.040 20 (high) 0.553 0.028 Table 4.6: Radiosensitivity parameters and ratios employed in the in vivo study. These values were set to yield the same D50 (dose at which TCP is 0.5) with conv-RT. 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) TCP αox/βox = 3 Gy conv-RT FLASH-RT 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) αox/βox = 10 Gy 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) αox/βox = 20 Gy Figure 4.8: Tumor Control Probability (TCP)-dose curves for heterogeneously oxygenated tumors with different αox/βox ratios, irradiated with conv-RT and FLASH-RT. Figure 4.8 presents the resulting TCP-dose curves for each group. Using the dose D50, described in Section 4.5.7, as a metric to compare the potential differences in TCP curves between both radiation modalities, the figure shows that FLASH-RT brings higher D50 values, as reported in Table 4.7. To account for the effect of the 90 Chapter 4. Results on FLASH radiotherapy αox/βox ratio on the differences in D50 between FLASH-RT and conv-RT, the associated BED was calculated from equation (1.11). Table 4.7 shows an increasing BED with decreasing αox/βox ratio. Thus, this study suggests that TCP loss with FLASH-RT becomes bigger in tumors with low αox/βox ratios. αox/βox (Gy) DC 50 (Gy) DF 50 (Gy) ∆D50 (Gy) BED(∆D50) 3 29.17 31.52 2.35 4.19 10 29.17 31.78 2.61 3.29 20 29.16 31.98 2.82 3.22 Table 4.7: Differences in D50 between conv-RT and FLASH-RT according to the αox/βox ratio. Differences are reported in grays (∆D50), and the corresponding Biologically Effective Dose (BED) to account for the influence of the radiosensitivity ratio. 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) TCP ˜p∈(0,10] conv-RT FLASH-RT 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) ˜p∈(10,20] 0 20 40 60 0 0.2 0.4 0.6 0.8 1 Dose (Gy) ˜p∈(20,30] Figure 4.9: Tumor Control Probability (TCP)-dose curves for conv-RT (solid lines) and FLASHRT (dashed lines) in tumors with heterogeneous oxygen levels according to their median oxygen partial pressure, ˜p:poorly oxygenated (˜p≤10 mmHg), moderately oxygenated (10 <˜p≤20 mmHg), and well oxygenated tumors (˜p > 20 mmHg). The αox/βox ratio was set to 10 Gy. A similar study was performed to investigate the influence of tumor oxygenation on TCP. For that purpose, the radiosensitivity ratio was set to αox/βox = 10 Gy and the corresponding tumor population was split into three groups according to their median oxygen distribution, ˜p:poorly oxygenated (˜p≤10 mmHg), moderately oxygenated (10 <˜p≤20 mmHg), and well oxygenated tumors (˜p > 20 mmHg). TCP-dose curves were obtained separately for each group, both for FLASH-RT and conv-RT, obtaining the curves presented in Figure 4.9. Moreover, the resulting D50 and D90 are summarized in Table 4.8 for both RT modalities, along with their relative ratio. When analyzing the 91 Isabel Gonz´ alez Crespo Future work would comprehend the incorporation to the response model of other factors that can compensate TCP loss due to ROD and maintain the iso-effectiveness between conv-RT and FLASH-RT, supporting the experimental evidence. 98 Bibliography [1] S. Gianfaldoni, R. Gianfaldoni, U. Wollina, et al., “An overview on radiotherapy: From its history to its current applications in dermatology”, Open access Macedonian Journal of Medical Sciences, vol. 5, no. 4, pp. 521–5, 2017. doi: 10.3889/oamjms.2017.122. [2] E. H. Grubb´e, “Priority in the therapeutic use of X-rays”, Radiology, vol. 21, no. 2, pp. 156–62, 1933. doi:10.1148/21.2.156. [3] H. Becquerel and P. Curie, “Action physiologique des rayons du radium”, Comptes Rendus de l’Acad´emie des Sciences, vol. 132, pp. 1289–91, 1901. [4] E. O. Lawrence and M. S. Livingston, “The production of high speed light ions without the use of high voltages”, Physical Review, vol. 40, no. 1, pp. 19–35, 1932. doi:10.1103/PhysRev.40.19. [5] H. Coutard, “Principles of X ray therapy of malignant diseases”, The Lancet, vol. 224, no. 5784, pp. 1–8, 1934. doi:10.1016/S0140-6736(00)90085-0. [6] D. W. Fry, R. B. R. -S. -Harvie, L. B. Mullett, et al., “A travelling-wave linear accelerator for 4-MeV. electrons”, Nature, vol. 162, no. 4126, pp. 859–61, 1948. doi:10.1038/162859a0. [7] F. G. Herrera, J. Bourhis, and G. Coukos, “Radiotherapy combination opportunities leveraging immunity for the next oncology practice: Radiationimmunotherapy combinations”, CA: A Cancer Journal for Clinicians, vol. 67, no. 1, pp. 65–85, 2017. doi:10.3322/caac.21358. [8] C. Park, L. Papiez, S. Zhang, et al., “Universal survival curve and single fraction equivalent dose: Useful tools in understanding potency of ablative radiotherapy”, International Journal of Radiation Oncology * Biology * Physics, vol. 70, no. 3, pp. 847–52, 2008. doi:10.1016/j.ijrobp.2007.10.059. [9] S. S. Lo, A. J. Fakiris, E. L. Chang, et al., “Stereotactic body radiation therapy: A novel treatment modality”, Nature Reviews Clinical Oncology, vol. 7, no. 1, pp. 44–54, 2010. doi:10.1038/nrclinonc.2009.188. Isabel Gonz´ alez Crespo [10] B. G. Douglas and J. F. Fowler, “The effect of multiple small doses of X rays on skin reactions in the mouse and a basic interpretation”, Radiation Research, vol. 66, no. 2, pp. 401–26, 1976. doi:10.2307/3574407. [11] C. W. Song, Y. J. Lee, R. J. Griffin, et al., “Indirect tumor cell death after highdose hypofractionated irradiation: Implications for stereotactic body radiation therapy and stereotactic radiation surgery”, International Journal of Radiation Oncology * Biology * Physics, vol. 93, no. 1, pp. 166–72, 2015. doi:10.1016/j .ijrobp.2015.05.016. [12] M. Garc´ıa-Barros, F. Paris, C. Cordon-Cardo, et al., “Tumor response to radiotherapy regulated by endothelial cell apoptosis”, Science, vol. 300, no. 5622, pp. 1155–9, 2003. doi:10.1126/science.1082504. [13] H. J. Park, R. J. Griffin, S. Hui, et al., “Radiation-induced vascular damage in tumors: Implications of vascular damage in ablative hypofractionated radiotherapy (SBRT and SRS)”, Radiation Research, vol. 177, no. 3, pp. 311–27, 2012. doi:10.1667/RR2773.1. [14] C. W. Song, E. Glatstein, L. B. Marks, et al., “Biological principles of stereotactic body radiation therapy (SBRT) and stereotactic radiation surgery (SRS): Indirect cell death”, International Journal of Radiation Oncology * Biology * Physics, vol. 110, no. 1, pp. 21–34, 2021. doi:10.1016/j.ijrobp.2019.02.0 47. [15] A. Filatenkov, J. Baker, A. M. S. Mueller, et al., “Ablative tumor radiation can change the tumor immune cell microenvironment to induce durable complete remissions”, Clinical Cancer Research, vol. 21, no. 16, pp. 3727–39, 2015. doi: 10.1158/1078-0432.CCR-14-2824. [16] L. de la Maza, M. Wu, L. Wu, et al., “In situ vaccination after accelerated hypofractionated radiation and surgery in a mesothelioma mouse model”, Clinical Cancer Research, vol. 23, no. 18, pp. 5502–13, 2017. doi:10.1158/1078-0432 .CCR-17-0438. [17] Y. Wang, Z. G. Liu, H. Yuan, et al., “The reciprocity between radiotherapy and cancer immunotherapy”, Clinical Cancer Research, vol. 25, no. 6, pp. 1709–17, 2019. doi:10.1158/1078-0432.CCR-18-2581. [18] L. Deng, H. Liang, B. Burnette, et al., “Irradiation and anti–PD-L1 treatment synergistically promote antitumor immunity in mice”, The Journal of Clinical Investigation, vol. 124, no. 2, pp. 687–95, 2014. doi:10.1172/JCI67313. [19] Y. Lee, S. L. Auh, Y. Wang, et al., “Therapeutic effects of ablative radiation on local tumor require CD8+ T cells: Changing strategies for cancer treatment”, Blood, The Journal of the American Society of Hematology, vol. 114, no. 3, pp. 589–95, 2009. doi:10.1182/blood-2009-02-206870. 100 Bibliography [20] J. He, Y. Yin, T. A. Luster, et al., “Antiphosphatidylserine antibody combined with irradiation damages tumor blood vessels and induces tumor immunity in a rat model of glioblastoma”, Clinical Cancer Research, vol. 15, no. 22, pp. 6871– 80, 2009. doi:10.1158/1078-0432.CCR-09-1499. [21] E. J. Moding, K. D. Castle, B. A. Perez, et al., “Tumor cells, but not endothelial cells, mediate eradication of primary sarcomas by stereotactic body radiation therapy”, Science Translational Medicine, vol. 7, no. 278, 278ra34, 2015. doi: 10.1126/scitranslmed.aaa4214. [22] J. A. Torok, P. Oh, K. D. Castle, et al., “Deletion of ATM in tumor but not endothelial cells improves radiation response in a primary mouse model of lung adenocarcinoma”, Cancer Research, vol. 79, no. 4, pp. 773–82, 2019. doi:10.1 158/0008-5472.CAN-17-3103. [23] K. Deland, J. S. Mercer, D. M. Crabtree, et al., “Radiosensitizing the vasculature of primary brainstem gliomas fails to improve tumor response to radiation therapy”, International Journal of Radiation Oncology * Biology * Physics, vol. 112, no. 3, pp. 771–9, 2022. doi:10.1016/j.ijrobp.2021.09.047. [24] J. M. Brown, D. J. Carlson, and D. J. Brenner, “The tumor radiobiology of SRS and SBRT: Are more than the 5 Rs involved?”, International Journal of Radiation Oncology * Biology * Physics, vol. 88, no. 2, pp. 254–62, 2014. doi: 10.1016/j.ijrobp.2013.07.022. [25] R. Ruggieri, P. Stavrev, S. Naccarato, et al., “Optimal dose and fraction number in SBRT of lung tumours: A radiobiological analysis”, Physica Medica, vol. 44, pp. 188–95, 2017. doi:10.1016/j.ejmp.2016.12.012. [26] J. P. Kirkpatrick, J. J. Meyer, and L. B. Marks, “The linear-quadratic model is inappropriate to model high dose per fraction effects in radiosurgery”, Seminars in Radiation Oncology, vol. 18, no. 4, pp. 240–3, 2008. doi:10.1016/j.semra donc.2008.04.005. [27] D. J. Brenner, “The linear-quadratic model is an appropriate methodology for determining isoeffective doses at large doses per fraction”, Seminars in Radiation Oncology, vol. 18, no. 4, pp. 234–39, 2008. doi:10.1016/j.semrad onc.2008.04.004. [28] C. W. Song, H. Park, R. J. Griffin, et al., “Radiobiology of stereotactic radiosurgery and stereotactic body radiation therapy”, in Technical basis of radiation therapy: Practical clinical applications. Springer Berlin Heidelberg, 2012, pp. 51–61, isbn: 9783642115721. doi:10.1007/174_2011_264. [29] P. W. Sperduto, C. W. Song, J. P. Kirkpatrick, et al., “A hypothesis: Indirect cell death in the radiosurgery era”, International Journal of Radiation Oncology * Biology * Physics, vol. 91, no. 1, pp. 11–3, 2015. doi:10.1016/j.ijrobp.2 014.08.355. 101 Isabel Gonz´ alez Crespo [30] C. W. Song, M. S. Kim, L. C. Cho, et al., “Radiobiological basis of SBRT and SRS”, International Journal of Clinical Oncology, vol. 19, no. 4, pp. 570–8, 2014. doi:10.1007/s10147-014-0717-z. [31] M. Guerrero and M. Carlone, “Mechanistic formulation of a lineal-quadraticlinear (LQL) model: Split-dose experiments and exponentially decaying sources”, Medical Physics, vol. 37, no. 8, pp. 4173–81, 2010. doi:10.1118/1.3456927. [32] M. Guerrero and X. A. Li, “Extending the linear-quadratic model for large fraction doses pertinent to stereotactic radiotherapy”, Physics in Medicine & Biology, vol. 49, no. 20, pp. 4825–35, 2004. doi:10.1088/0031-9155/49/20/0 12. [33] A. Gago-Arias, S. Neira, M. Pombar, et al., “Evaluation of indirect damage and damage saturation effects in dose-response curves of hypofractionated radiotherapy of early-stage NSCLC and brain metastases”, Radiotherapy and Oncology, vol. 161, pp. 1–8, 2021. doi:10.1016/j.radonc.2021.05.012. [34] R. Serre, S. Benzekry, L. Padovani, et al., “Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy”, Cancer Research, vol. 76, no. 17, pp. 4931–40, 2016. doi:10.1158/0008-5472.CAN-15-3567. [35] J. Poleszczuk and H. Enderling, “The optimal radiation dose to induce robust systemic anti-tumor immunity”, International Journal of Molecular Sciences, vol. 19, no. 11, p. 3377, 2018. doi:10.3390/ijms19113377. [36] S. J. McMahon, “The linear quadratic model: Usage, interpretation and challenges”, Physics in Medicine & Biology, vol. 64, no. 1, 01TR01, 2019. doi: 10.1088/1361-6560/aaf26a. [37] A. Gago-Arias, P. Aguiar, I. Espinoza, et al., “Modelling radiation-induced cell death and tumour re-oxygenation: Local versus global and instant versus delayed cell death”, Physics in Medicine & Biology, vol. 61, no. 3, pp. 1204–16, 2016. doi:10.1088/0031-9155/61/3/1204. [38] B. J. Moeller, R. A. Richardson, and M. W. Dewhirst, “Hypoxia and radiotherapy: Opportunities for improved outcomes in cancer treatment”, Cancer and Metastasis Reviews, vol. 26, pp. 241–48, 2007. doi:10.1007/s10555-007 -9056-0. [39] A. S. E. Ljungkvist, J. Bussink, P. F. J. W. Rijken, et al., “Vascular architecture, hypoxia, and proliferation in first-generation xenografts of human headand-neck squamous cell carcinomas”, International Journal of Radiation Oncology * Biology * Physics, vol. 54, no. 1, pp. 215–28, 2002. doi:10.1016/S03 60-3016(02)02938-3. [40] P. Vaupel and A. Mayer, “Hypoxia in cancer: Significance and impact on clinical outcome”, Cancer and Metastasis Reviews, vol. 26, pp. 225–39, 2007. doi:10 .1007/s10555-007-9055-1. 102 Bibliography [41] D. R. Grimes and M. Partridge, “A mechanistic investigation of the oxygen fixation hypothesis and oxygen enhancement ratio”, Biomedical physics & engineering express, vol. 1, no. 4, p. 045 209, 2015. doi:10.1088/2057-1976/1 /4/045209. [42] B. G. Wouters and J. M. Brown, “Cells at intermediate oxygen levels can be more important than the “hypoxic fraction” in determining tumor response to fractionated radiotherapy”, Radiation Research, vol. 147, no. 5, pp. 541–50, 1997. doi:10.2307/3579620. [43] K. Sachs, P. Hahnfeld, and D. J. Brenner, “Review the link between low-LET dose-response relations and the underlying kinetics of damage production/repair/misrepair”, International Journal of Radiation Biology, vol. 72, no. 4, pp. 351–74, 1997. doi:10.1080/095530097143149. [44] J. F. Fowler, “The linear-quadratic formula and progress in fractionated radiotherapy”, The British Journal of Radiology, vol. 62, no. 740, pp. 679–94, 1989. doi:10.1259/0007-1285-62-740-679. [45] G. W. Barendsen, “Dose fractionation, dose rate and iso-effect relationships for normal tissue responses”, International Journal of Radiation Oncology * Biology * Physics, vol. 8, no. 11, pp. 1981–97, 1982. doi:10.1016/0360-3016 (82)90459-X. [46] S. Webb and A. E. Nahum, “A model for calculating tumour control probability in radiotherapy including the effects of inhomogeneous distributions of dose and clonogenic cell density”, Physics in Medicine & Biology, vol. 38, no. 6, pp. 653– 66, 1993. doi:10.1088/0031-9155/38/6/001. [47] S. Webb and A. E. Nahum, “A model for calculating tumour control probability in radiotherapy including the effects of inhomogeneous distributions of dose and clonogenic cell density”, Physics in Medicine & Biology, vol. 38, no. 6, pp. 653– 66, 1993. doi:10.1088/0031-9155/38/6/001. [48] T. L. Whiteside, S. Demaria, M. E. Rodriguez-Ruiz, et al., “Emerging opportunities and challenges in cancer immunotherapy”, Clinical Cancer Research, vol. 22, no. 8, pp. 1845–55, 2016. doi:10.1158/1078-0432.CCR-16-0049. [49] T. A. Waldmann, “Immunotherapy: Past, present and future”, Nature Medicine, vol. 9, no. 3, pp. 269–77, 2003. doi:10.1038/nm0303-269. [50] D. S. Chen and I. Mellman, “Oncology meets immunology: The cancer-immunity cycle”, Immunity, vol. 39, no. 1, pp. 1–10, 2013. doi:10.1016/j.immuni.201 3.07.012. [51] I. Mellman, G. Coukos, and G. Dranoff, “Cancer immunotherapy comes of age”, Nature, vol. 480, no. 7378, pp. 480–9, 2011. doi:10.1038/nature10673. 103 Isabel Gonz´ alez Crespo [52] E. J. Lipson and C. G. Drake, “Ipilimumab: An anti-CTLA-4 antibody for metastatic melanoma”, Clinical Cancer Research, vol. 17, no. 22, pp. 6958–62, 2011. doi:10.1158/1078-0432.CCR-11-1595. [53] K. Pang, Z. D. Shi, L. Y. Wei, et al., “Research progress of therapeutic effects and drug resistance of immunotherapy based on PD-1/PD-L1 blockade”, Drug Resistance Updates, vol. 66, p. 100 907, 2023. doi:10.1016/j.drup.2022.100 907. [54] S. C. Formenti and S. Demaria, “Systemic effects of local radiotherapy”, The Lancet Oncology, vol. 10, no. 7, pp. 718–26, 2009. doi:10.1016/S1470-2045 (09)70082-8. [55] S. Demaria, B. Ng, M. L. Devitt, et al., “Ionizing radiation inhibition of distant untreated tumors (abscopal effect) is immune mediated”, International Journal of Radiation Oncology * Biology * Physics, vol. 58, no. 3, pp. 862–70, 2004. doi: 10.1016/j.ijrobp.2003.09.012. [56] L. L. S. C. Apetoh, S. Ladoire, G. Coukos, et al., “Combining immunotherapy and anticancer agents: The right path to achieve cancer cure?”, Annals of Oncology, vol. 26, no. 9, pp. 1813–23, 2015. doi:10.1093/annonc/mdv209. [57] C. Tang, X. Wang, H. Soh, et al., “Combining radiation and immunotherapy: A new systemic therapy for solid tumors?”, Cancer Immunology Research, vol. 2, no. 9, pp. 831–8, 2014. doi:10.1158/2326-6066.CIR-14-0069. [58] M. Grapin, C. Richard, E. Limagne, et al., “Optimized fractionated radiotherapy with anti-PD-L1 and anti-TIGIT: A promising new combination”, Journal for Immunotherapy of Cancer, vol. 7, no. 160, pp. 1–12, 2019. doi:10.1186/s 40425-019-0634-9. [59] M. Z. Dewan, A. E. Galloway, N. Kawashima, et al., “Fractionated but not single-dose radiotherapy induces an immune-mediated abscopal effect when combined with anti-CTLA-4 antibody”, Clinical Cancer Research, vol. 15, no. 17, pp. 5379–88, 2009. doi:10.1158/1078-0432.CCR-09-0265. [60] R. A. Bekker, S. Kim, S. Pilon-Thomas, et al., “Mathematical modeling of radiotherapy and its impact on tumor interactions with the immune system”, Neoplasia, vol. 28, p. 100 796, 2022. doi:10.1016/j.neo.2022.100796. [61] M. Durante, E. Br¨auer-Krisch, and M. Hill, “Faster and safer? FLASH ultrahigh dose rate in radiotherapy”, The British Journal of Radiology, vol. 91, no. 1082, p. 20 170 628, 2017. doi:10.1259/bjr.20170628. [62] V. Favaudon, R. Labarbe, and C. L. Limoli, “Model studies of the role of oxygen in the FLASH effect”, Medical Physics, vol. 49, no. 3, pp. 2068–81, 2022. doi: 10.1002/mp.15129. 104 Bibliography [63] D. L. Dewey and J. W. Boag, “Modification of the oxygen effect when bacteria are given large pulses of radiation”, Nature, vol. 183, no. 4673, pp. 1450–1, 1959. doi:10.1038/1831450a0. [64] B. Lin, F. Gao, Y. Yang, et al., “FLASH radiotherapy: History and future”, Frontiers in Oncology, vol. 11, p. 1890, 2021. doi:10.3389/fonc.2021.644400. [65] E. Bogaerts, E. Macaeva, S. Isebaert, et al., “Potential molecular mechanisms behind the ultra-high dose rate “FLASH” effect”, International Journal of Molecular Sciences, vol. 23, no. 20, p. 12 109, 2022. doi:10.3390/ijms232 012109. [66] V. Favaudon, L. Caplier, V. Monceau, et al., “Ultrahigh dose-rate FLASH irradiation increases the differential response between normal and tumor tissue in mice”, Science Translational Medicine, vol. 6, no. 245, 245ra93, 2014. doi: 10.1126/scitranslmed.3008973. [67] P. Montay-Gruel, K. Petersson, M. Jaccard, et al., “Irradiation in a flash: Unique sparing of memory in mice after whole brain irradiation with dose rates above 100Gy/s”, Radiotherapy and Oncology, vol. 124, no. 3, pp. 365–9, 2017. doi:10.1016/j.radonc.2017.05.003. [68] P. Montay-Gruel, A. Bouchet, M. Jaccard, et al., “X-rays can trigger the FLASH effect: Ultra-high dose-rate synchrotron light source prevents normal brain injury after whole brain irradiation in mice”, Radiotherapy and Oncology, vol. 129, no. 3, pp. 582–8, 2018. doi:10.1016/j.radonc.2018.08.016. [69] M. C. Vozenin, P. De Fornel, K. Petersson, et al., “The advantage of FLASH radiotherapy confirmed in mini-pig and cat-cancer patients”, Clinical Cancer Research, vol. 25, no. 1, pp. 35–42, 2019. doi:10.1158/1078-0432.CCR-17-3 375. [70] P. Montay-Gruel, M. M. Acharya, K. Petersson, et al., “Long-term neurocognitive benefits of FLASH radiotherapy driven by reduced reactive oxygen species”, Proceedings of the National Academy of Sciences, vol. 116, no. 22, pp. 10 943– 51, 2019. doi:10.1073/pnas.1901777116. [71] K. Levy, S. Natarajan, J. Wang, et al., “Abdominal FLASH irradiation reduces radiation-induced gastrointestinal toxicity for the treatment of ovarian cancer in mice”, Scientific Reports, vol. 10, no. 21600, pp. 1–14, 2020. doi:10.1038 /s41598-020-78017-7. [72] E. S. Diffenderfer, I. I. Verginadis, M. M. Kim, et al., “Design, implementation, and in vivo validation of a novel proton FLASH radiation therapy system”, International Journal of Radiation Oncology * Biology * Physics, vol. 106, no. 2, pp. 440–8, 2020. doi:10.1016/j.ijrobp.2019.10.049. 105 Isabel Gonz´ alez Crespo [73] E. Liljedahl, E. Konradsson, E. Gustafsson, et al., “Long-term anti-tumor effects following both conventional radiotherapy and FLASH in fully immunocompetent animals with glioblastoma”, Scientific Reports, vol. 12, no. 12285, pp. 1–12, 2022. doi:10.1038/s41598-022-16612-6. [74] F. Gao, Y. Yang, H. Zhu, et al., “First demonstration of the FLASH effect with ultrahigh dose rate high-energy X-rays”, Radiotherapy and Oncology, vol. 166, pp. 44–50, 2022. doi:10.1016/j.radonc.2021.11.004. [75] G. Pratx and D. S. Kapp, “A computational model of radiolytic oxygen depletion during FLASH irradiation and its effect on the oxygen enhancement ratio”, Physics in Medicine & Biology, vol. 64, no. 18, p. 185 005, 2019. doi: 10.1088/1361-6560/ab3769. [76] G. Pratx and D. S. Kapp, “Ultra-high-dose-rate FLASH irradiation may spare hypoxic stem cell niches in normal tissues”, International Journal of Radiation Oncology * Biology * Physics, vol. 105, no. 1, pp. 190–2, 2019. doi:10.1016/j .ijrobp.2019.05.030. [77] K. Petersson, G. Adrian, K. Butterworth, et al., “A quantitative analysis of the role of oxygen tension in FLASH radiation therapy”, International Journal of Radiation Oncology * Biology * Physics, vol. 107, no. 3, pp. 539–47, 2020. doi:10.1016/j.ijrobp.2020.02.634. [78] X. Cao, R. Zhang, T. V. Esipova, et al., “Quantification of oxygen depletion during FLASH irradiation in vitro and in vivo”, International Journal of Radiation Oncology * Biology * Physics, vol. 111, no. 1, pp. 240–8, 2021. doi: 10.1016/j.ijrobp.2021.03.056. [79] J. Y. Jin, A. Gu, W. Wang, et al., “Ultra-high dose rate effect on circulating immune cells: A potential mechanism for FLASH effect?”, Radiotherapy and Oncology, vol. 149, pp. 55–62, 2020. doi:10.1016/j.radonc.2020.04.054. [80] D. R. Spitz, G. R. Buettner, M. S. Petronek, et al., “An integrated physicochemical approach for explaining the differential impact of FLASH versus conventional dose rate irradiation on cancer and normal tissue responses”, Radiotherapy and Oncology, vol. 139, pp. 23–7, 2019. doi:10.1016/j.radonc.2019 .03.028. [81] R. Abolfath, A. Baikalov, S. Rahvar, et al., “Differential tissue sparing of FLASH ultra high dose rates: An in-silico study”, arXiv:2210.03565, 2022. doi:10.48550/arXiv.2210.03565. [82] H. Song, Y. Kim, and W. Sung, “Modeling of the FLASH effect for ion beam radiation therapy”, Physica Medica, vol. 108, p. 102 553, 2023. doi:10.1016/j .ejmp.2023.102553. 106 Bibliography [83] R. Labarbe, L. Hotoiu, J. Barbier, et al., “A physicochemical model of reaction kinetics supports peroxyl radical recombination as the main determinant of the FLASH effect”, Radiotherapy and Oncolgy, vol. 153, pp. 303–10, 2020. doi: 10.1016/j.radonc.2020.06.001. [84] E. Taylor, R. P. Hill, and D. L´etourneau, “Modeling the impact of spatial oxygen heterogeneity on radiolytic oxygen depletion during FLASH radiotherapy”, Physics in Medicine & Biology, vol. 67, no. 11, p. 115 017, 2022. doi: 10.1088/1361-6560/ac702c. [85] D. D. Bainov and P. S. Simeonov, Impulsive differential equations: Asymptotic properties of the solutions. Singapore: World Scientific, 1995, isbn: 97898102182 32. doi:10.1142/2413. [86] B. Randjelovic, L. V. Stefanovic, and B. M. Dankovic, “Numerical solution of impulsive differential equations”, Facta Universitatis. Series Mathematics and Informatics, vol. 15, pp. 101–11, 2000. [87] N. S. B. A. Hamzah, M. Mamat, J. Kavikumar, et al., “Impulsive differential equations by using the Euler method”, Applied Mathematical Sciences, vol. 4, no. 65, pp. 3219–32, 2010. [88] I. Gy¨ori and G. Ladas, Oscillation theory of delay differential equations: With applications, 1st ed. Oxford: Clarendon Press, 1992, isbn: 9780198535829. doi: 10.1093/oso/9780198535829.001.0001. [89] L. Euler, Institutiones calculi integralis. Saint Petersburg: Academia Imperialis Scientiarum, 1768, vol. 1. [90] J. C. Butcher, Numerical methods for ordinary differential equations, 2nd ed. John Wiley &Sons, 2016, isbn: 9780470723357. doi:10.1002/9781119121534. [91] V. Wulf and N. J. Ford, “Insight into the qualitative behaviour of numerical solutions to some delay differential equations”, Proceedings of HERCMA, pp. 629–36, 1998. [92] P. D´ıaz-Botana, “Modelizaci´on de hipoxia cr´onica/aguda en tumores mediante la resoluci´on de la ecuaci´on de reacci´on-difusi´on y su efecto en tratamientos radioter´apicos”, M.S. thesis, Universidade de Santiago de Compostela, 2016. [93] I. Espinoza, P. Peschke, and C. P. Karger, “A model to simulate the oxygen distribution in hypoxic tumors for different vascular architectures”, Medical Physics, vol. 40, no. 8, p. 081 703, 2013. doi:10.1118/1.4812431. [94] A. Da¸su, I. Toma-Da¸su, and M. Karlsson, “Theoretical simulation of tumour oxygenation and results from acute and chronic hypoxia”, Physics in Medicine & Biology, vol. 48, no. 17, pp. 2829–42, 2003. doi:10.1088/0031-9155/48/1 7/307. 107 Isabel Gonz´ alez Crespo Fig. 1.7 Step 6 of the immunity cycle. T-cells detect tumor cells regulated by the Programmed Death 1 (PD-1)/Programmed Death-Ligand 1 (PD-L1) pathway and ImmunoTherapy (IT) with anti-Programmed Death-(Ligand) 1 (αPD(L)1). [Created in BioRender.com.] ..... 24 Fig. 2.1 Example of a two-dimensional tumor domain. Geometry for the numerical simulation of the oxygenation problem. The inner circles represent blood capillaries of the tumor vascular system. This geometry corresponds to a VF vf= 0.14. ..................... 35 Fig. 2.2 Steady-state solution of the oxygenation problem 2.1 with the consumption term given by equation (2.9) for the tumor geometry shown in Figure 2.1. ............................... 40 Fig. 2.3 Oxygen depletion during conventional RT (conv-RT) (0.1 Gy s−1) and FLASH-RT (100 Gy s−1) for a delivered dose of 20 Gy and the geometry depicted at Figure 2.1. The left panel shows the depletion curves for both modalities at different initial oxygen partial pressures, p. The right panel shows the change in oxygenation histograms from the beginning (t=t0) to the end (t=T) of FLASH-RT. These results were obtained by using the parameter values given in previous sections, along with G0= 0.25 mmHg Gy−1and kROD = 1 mmHg. . 41 Fig. 2.4 Flowchart of the Simulated Annealing (SA) algorithm. ........ 47 Fig. 3.1 Flowchart of the model showing the interaction between different compartments: viable tumor cells (C), doomed tumor cells (Cd), active T-cells in the tumor (Ta), antigens/Antigen Presenting Cells (APCs) (ˆ A), the pool of T-cells ( ˆ T), activated T-cells ( ˆ Ta) and blocked T-cells (ˆ Tb). The dotted line shows the separation between compartments physically located in the tumor (left) and in the activation sites (right). The notation also distinguishes between the two locations, identifying the compartments within the activation zone with a hat (ˆ). τA and τTare biological delays related to the migration of antigens and T-cells between those two physical locations. The action points of RadioTherapy (RT), anti-Programmed Death-(Ligand) 1 (αPD(L)1), and anti-CTLA-4 (αCTLA4) are also indicated. ............ 52 Fig. 3.2 Administration schedule of treatments combining different RadioTherapy (RT) fractionations and ImmunoTherapy (IT) delivery times with αCTLA4, experimentally studied by Dewan et al. [59]. ........ 60 Fig. 3.3 Administration schedule of treatments combining different RadioTherapy (RT) fractionations and ImmunoTherapy (IT) delivery times with αPDL1, experimentally studied by Deng et al. [18]. .......... 60 114 List of Figures Fig. 3.4 Sketch of the process described by the Markov birth-death model. Cis the number of clonogenic tumor cells, Ppis the proliferation probability and Pdis the death probability. .............. 63 Fig. 3.5 Model fitting to experimental data reported by Dewan et al. [59] of tumor response to RadioTherapy (RT), ImmunoTherapy (IT) with αCTLA4, and combined treatments. The LQmod model (3.2) accounts for the radiosensitivity of tumor cells. Radiation doses were delivered on consecutive days starting from day 0. Notice that differences between control and IT curves are small and both of them overlap in the figure. All the fits were obtained with a single set of parameters, although they are plotted separately to facilitate the visualization. The control curve is included in all the panels as a common reference. . . 64 Fig. 3.6 Model fitting to experimental data reported by Dewan et al. [59] of tumor response to RadioTherapy (RT), ImmunoTherapy (IT) with αCTLA4, and combined treatments. The LQL model (1.7) accounts for the radiosensitivity of tumor cells. Radiation doses were delivered on consecutive days starting from day 0. Notice that differences between control and IT curves are small and both of them overlap in the figure. All the fits were obtained with a single set of parameters, although they are plotted separately to facilitate the visualization. The control curve is included in all the panels as a common reference. . . 65 Fig. 3.7 Model fitting to experimental data reported by Deng et al. [18] of tumor response to RadioTherapy (RT) (12 Gy single fraction on day 0), ImmunoTherapy (IT) with αPDL1, and combined treatments. The LQ model (1.1) accounts for the radiosensitivity of tumor cells. Notice that differences between control and IT curves are small and both of them overlap in the figure. ........................ 66 Fig. 3.8 Model fitting to experimental data reported by Dewan et al. [59] of tumor response to RadioTherapy (RT), ImmunoTherapy (IT) with αCTLA4, and combined treatments. The LQ model characterizes tumor cell response to radiation, and vascular damage at 20 Gy per fraction was included to limit T-cell infiltration in the tumor. Radiation doses were delivered on consecutive days starting from day 0. Notice that differences between control and IT curves are small and both of them overlap in the figure. All the fits were obtained with a single set of parameters, although they are plotted separately to facilitate the visualization. The control curve is included in all the panels as a common reference. .......................... 67 115 Isabel Gonz´ alez Crespo Fig. 3.9 Study of optimal schedules of RT (8 Gy×3 fractions) and IT (3 fractions of αCTLA4) obtained from the presented model and best-fitting parameters reported in Table D.4. RT fractions are delivered at days (0, 1, 2), and IT is delivered with different schedules starting from day 0 to day 7. Panels (a), (b), and (c) report the dynamics of tumor volumes and the 95% confidence intervals for the combinations investigated by Dewan et al. [59], delivering IT at days (0, 3, 6), (2, 5, 8), and (4, 6, 9) respectively. Panel (d) shows the comparative between 95% confidence intervals for Tumor Control Probability (TCP) values obtained with the model (200 simulations) and experimental controls (6 to 16 animals), for the same treatment combinations. Panel (e) presents TCP (asterisks) versus the concentration of αCTLA4 on day 2 after starting treatment, showing a positive correlation between those two variables. The solid line corresponds to the fit of a logistic function. .................................. 69 Fig. 3.10 Modeled responses to different RadioTherapy (RT) schedules with αCTLA4 (2, 5, 8). RT fractionations are equivalent from a classical radiobiological point of view (same BED). Tumor volume evolution is shown for RT alone (a), and RadioImmunoTherapy (RIT) with αCTLA4 (b). Tumor Control Probability (TCP) (95% confidence intervals) were obtained from the model (200 simulations) for RT alone (c), and RIT (d). The model parameters used for this study are presented in Table D.4, and include vascular damage effect at 16.25 Gy, governed by equation (3.13). ....................... 70 Fig. 3.11 Evolution of tumor cells and active T-cells in the tumor zone for the different fractionations: single-dose, non-single-dose hypofractionation, and more fractionated treatment. The panels on the left correspond to RadioTherapy (RT) and those on the right to RadioImmunoTherapy (RIT) with αCTLA4 (2 5 8). Non-single-dose hypofractionated treatments appear to be more effective, as they allow T-cell infiltration into the tumor and reduce T-cell damage. ......... 72 Fig. 4.1 Oxygen depletion curves obtained by fitting equation (4.1) to the measurements reported by Jansen et al. [100] (i, ii and iii) in water, Van Slyke et al. [101] (iv) in CELL solution, and El Khatib et al. [99] (v) in BSA solution. The total amount of depleted oxygen during FLASH-RT, ∆p, divided by the delivered radiation dose is presented against the initial oxygen partial pressure. ............... 83 116 List of Figures Fig. 4.2 Amount of depleted oxygen during FLASH-RT (30 Gy, 100 Gy s−1) versus the initial mean oxygen partial pressure, ¯p, in preclinical tumors. The triangles represent in vivo data reported by Van Slyke et al. [101], and the circles represent the simulated data for 100 tumors, obtained by solving Problem 2.3. Linear fits of each dataset are presented as solid and dashed lines, respectively. ................. 84 Fig. 4.3 Ratio of Surviving Fractions (SFs) between FLASH-RT (SFF) and conv-RT (SFC) for a dose of 20 Gy versus the oxygenation status of the cells. A dose rate of 100 Gy s−1was employed to obtain SFF.. . 85 Fig. 4.4 Surviving Fraction (SF) versus dose curves for cells irradiated with conv-RT and FLASH-RT at two oxygenation levels (poorly oxygenated, p= 3 mmHg, moderately-well oxygenated,p= 15 mmHg). Results are presented for αox/βox ∈ {3,10,20,∞} Gy. ............. 86 Fig. 4.5 Best fits of tumor growth curves (mean values and standard deviations) for control, FLASH-RT, and conv-RT groups, reported by Diffenderfer et al. [72] ((a) MH641905 12 Gy, (b) MH641905 18 Gy) and Zhu et al. [140] ((c) Py8119 9.5 Gy, (d) Py230 9.5 Gy). ........ 88 Fig. 4.6 Example of the comparison between tumor growth curves after irradiation with FLASH-RT and conv-RT from fits to the data reported by Diffenderfer et al. [72] ((a) MH641905 12 Gy, (b) MH641905 18 Gy) and Zhu et al. [140] ((c) Py8119 9.5 Gy, (d) Py230 9.5 Gy). The solid lines represent the mean values and the shadow areas represent the standard deviation of populations with the sample size given in Table 4.2 for each experiment. ...................... 89 Fig. 4.7 Percentage of simulated experiments showing a significant difference (p-value <0.05) between control, conv-RT, and FLASH-RT groups, for simulations based on the experimental data reported by Diffenderfer et al. [72] and Zhu et al. [140], summarized in Table 4.2. . . . . 89 Fig. 4.8 Tumor Control Probability (TCP)-dose curves for heterogeneously oxygenated tumors with different αox/βox ratios, irradiated with convRT and FLASH-RT. ........................... 90 Fig. 4.9 Tumor Control Probability (TCP)-dose curves for conv-RT (solid lines) and FLASH-RT (dashed lines) in tumors with heterogeneous oxygen levels according to their median oxygen partial pressure, ˜p:poorly oxygenated (˜p≤10 mmHg), moderately oxygenated (10 <˜p≤20 mmHg), and well oxygenated tumors (˜p > 20 mmHg). The αox/βox ratio was set to 10 Gy. ................................ 91 117 Isabel Gonz´ alez Crespo Fig. A.1 Convergence study of the mesh to solve Problem 2.4. The left panel presents absolute error measurements for different numbers of nodes compared to the finest computed mesh (∼60 000 nodes). The right panel shows the respective execution times obtained from running the code in Listing 2.5. ............................ 132 Fig. A.2 Example of a mesh to apply the FEM to solve the oxygenation problem on the domain shown in Figure 2.1. ................... 132 Fig. B.1 Phase portrait of the system (B.4) using λ= 1, ϕ= 1, and K= 3. The critical points are shown in red. .................. 135 Fig. C.1 Tumor volume evolution with different time steps, ∆t(days), using the parameters summarized in Table D.4. In the left panel, curves for ∆t∈ {0.05,0.01,0.005,0.001}overlap, while in the right panel, curves for ∆t∈ {0.01,0.005,0.001}overlap. .................. 138 Fig. C.2 Absolute error when calculating the volume (at day 10 post-start of treatment) with different time steps (asterisks), and linear fits of the obtained values (solid line). The solution with ∆t= 0.001 days served as a reference to calculate the global error. ............... 138 Fig. C.3 Tumor volume evolution with different time steps, ∆t(days), using the parameters summarized in Table E.1). In both panels, curves for ∆t∈ {0.05,0.01,0.005,0.001}overlap. ................. 139 Fig. C.4 Absolute error when calculating the volume (at day 10 post-start of treatment) with different time steps (asterisks), and linear fits of the obtained values (solid line). The solution with ∆t= 0.001 days served as a reference to calculate the global error. ............... 140 Fig. D.1 Contribution of tumor cells and T-cells to tumor volumes in fits of the biomathematical model to experimental data of RIT with αCTLA4 reported by Dewan et al. [59]. These curves were obtained from fits reported in Figure 3.8, where the classical LQ model accounts for direct cell death, and the effect of vascular damage was taken into account. Tumor volumes are dominated by tumor cells, as expected. Numbers within parentheses indicate the delivery times of αCTLA4, if any. ................................... 145 118 List of Figures Fig. D.2 Relative contribution of T-cells (dashed lines) to tumor volumes in fits of the biomathematical model to experimental data of RIT with αCTLA4 reported by Dewan et al. [59]. These curves were obtained from fits reported in Figure 3.8, where the classical LQ model accounts for direct cell death, and the effect of vascular damage was taken into account. Tumor volumes are dominated by tumor cells, and only when the tumor is close to remission, the T-cell fraction becomes large. Numbers within parentheses indicate the delivery times of αCTLA4, if any. ................................... 146 Fig. E.1 Example of solutions to Problem 2.4 on a squared domain of 1 mm2[95] for two Vascular Fractions (VFs), fv, leading to different oxygen distributions. The top panels show the spatial distribution of the oxygenation and capillaries, while the bottom panels show the associated histograms. ................................ 148 Fig. E.2 Simulated heterogeneous oxygenations of each distribution in the H sample (rows). The numbers on each cell represent the percentage of tumor area for each oxygen partial pressure, p. Darker colors indicate higher frequencies. ............................ 149 Fig. E.3 Comparison of Tumor Control Probability (TCP)-dose curves between conv-RT and FLASH-RT for homogeneously oxygenated tumors with αox/βox = 10 Gy and oxygen partial pressure ranging from 3 mmHg to 27 mmHg. ............................... 150 119 List of Tables Tab. 2.1 Characteristics of the vascular vessels within the tumor, corresponding to experimental measurements of capillaries in the center of colorectal tumors [97]. Mean and standard deviation refer to the lognormal distribution. ........................... 35 Tab. 3.1 Analysis of the most critical model parameters. The sensitivity index given by equation (3.21) for each parameter is calculated as the difference between the cost function of fits to experimental data presented by Dewan et al. [59] (best-fitting parameters reported in Table D.4) and the cost function obtained when a 10% perturbation is applied to that particular parameter. Only the ten most critical parameters are shown. For reference, the best cost value was F= 25.3345. . . . 68 Tab. 3.2 Comparative between different fractionated RadioTherapy (RT) schedules regarding T-cell’s activation, infiltration, and radiation-mediated death. While extreme-hypofractionated regimens may compromise the therapeutic effect by reducing T-cell infiltration, conventional treatments might be also suboptimal by promoting the elimination of T-cells due to multiple irradiations. ................. 71 Tab. 4.1 Characteristic of the collected experimental data of oxygen depletion. Notice that none of the experiments were performed on cellular compounds, but over different mediums that partially mimic the intracellular milieu. Van Slyke et al. used the CELL aqueous solution composed of glycerol, glucose, glutathione, and HEPES; while El Khatib et al. employed a 5% Bovine Serum Albumin (BSA) and phosphate solution. ........................... 80 Isabel Gonz´ alez Crespo Tab. 4.2 Characteristic of the collected experimental data of tumor evolution in mice. *In the experiments by Zhu et al., the absorbed dose was measured having, for a FLASH-RT planned dose of 10 Gy, received doses of 9.75 Gy, 9.50 Gy, 9.81 Gy and 9.36 Gy (following the table order), which were approximated by 9.5 Gy. .............. 80 Tab. 4.3 List of parameters in the oxygenation model fixed to values from the literature, and references to the corresponding studies. ........ 81 Tab. 4.4 Best-fitting values of G0and kROD, and 95% confidence intervals obtained by fitting (4.1) to the oxygen depletion curves reported by El Khatib et al. [99], Jansen et al. [100], and Van Slyke et al. [101]. The fits are shown in Figure 4.1. .................... 84 Tab. 4.5 Radiosensitivity parameters and ratios employed in the in vitro study. These values were set to be iso-effective for a dose of 20 Gy with convRT. .................................... 85 Tab. 4.6 Radiosensitivity parameters and ratios employed in the in vivo study. These values were set to yield the same D50 (dose at which TCP is 0.5) with conv-RT. ............................ 90 Tab. 4.7 Differences in D50 between conv-RT and FLASH-RT according to the αox/βox ratio. Differences are reported in grays (∆D50), and the corresponding Biologically Effective Dose (BED) to account for the influence of the radiosensitivity ratio. ................. 91 Tab. 4.8 Differences in D50 and D90 between conv-RT (DC 50,DC 90) and FLASHRT (DF 50,DF 90) in tumors with heterogeneous oxygen levels according to their median oxygen partial pressures, ˜p. The αox/βox ratio was set to 10 Gy. ............................... 92 Tab. D.1 Parameter restrictions used for model fitting. Parameters not included in the table were constrained to be positive. .......... 141 Tab. D.2 Best-fitting parameters of the model to experimental data reported by Dewan et al. [59] and Deng et al. [18]. The LQmod model, given by equation (3.2), was used to account for tumor cell radiosensitivity. The symbol * indicates that the parameter value was fixed and not included in the optimization; while ** indicates that the parameter was optimized, but the same value was used for αCTLA4 and αPDL1 fittings. Initial numbers/concentrations of variables not indicated in the table (Cd,ˆ A,Ta,ˆ Ta,ˆ Tb,c4,p1) were set to zero. The parameter ris defined as r= 1 + b/a [34]. ..................... 142 122 List of Tables Tab. D.3 List of best-fitting parameters of the model to experimental data reported by Dewan et al. [59]. The LQL model (equation (1.7)) was used to account for tumor cell radiosensitivity. The symbol * indicates that the parameter value was fixed and not included in the optimization. Initial numbers/concentrations of variables not indicated in the table (Cd,ˆ A,Ta,ˆ Ta,ˆ Tb,c4,p1) were set to zero. The parameter ris defined as r= 1 + b/a [34]. ............... 143 Tab. D.4 List of best-fitting parameters of the model to experimental data reported by Dewan et al. [59]. The LQ model (equation (1.1)) was used to account for tumor cell radiosensitivity, and limited T-cell infiltration was assumed for 20 Gy irradiation due to vascular damage (equation (3.13)). The symbol * indicates that the parameter value was fixed and not included in the optimization. Initial numbers/concentrations of variables not indicated in the table (Cd,ˆ A,Ta,ˆ Ta,ˆ Tb, c4,p1) were set to zero. The parameter ris defined as r= 1 + b/a [34]. 144 Tab. E.1 Best-fitting parameters for the tumor growth curves reported by Diffenderfer et al. [72] and Zhu et al. [140]. The characteristics of the experiments are summarized in Table 4.2 ((a) MH641905 12 Gy, (b) MH641905 18 Gy, (c) Py8119 9.5 Gy, (d) Py230 9.5 Gy). ...... 147 123 Appendix A Mesh convergence study to apply the FEM Meshing to apply the FEM requires finding a good compromise between mesh quality and the resulting computational cost of solving the associated problem. A good-quality mesh is defined by uniform-shaped elements and smooth size transitions between consecutive elements. Low-quality meshes usually require fewer elements and decrease the execution time but may introduce errors and instability in the solution. Therefore, this appendix presents a mesh convergence study to solve Problem 2.4 applied to the domain shown in Figure 2.1. The FEM was run with different mesh sizes and the solutions were compared to a reference solution for the finest computed mesh (∼60 000 nodes) by calculating the absolute error in the mean partial pressure over the whole domain, p. The results are shown in the left panel of Figure A.1. Moreover, the execution time of running the code in Listing 2.5 is shown in the right panel. In this thesis, meshes with approximately 10 000 nodes were considered wellbalanced in terms of quality and computational cost. Finally, the MATLAB function meshQuality was employed to evaluate the shape quality of the elements in a representative mesh with 10 496 nodes and 19 621 elements, illustrated in Figure A.2. This function assigns a number within the range 0–1 to each element, where 1 corresponds to the optimal shape, being an equilateral triangle in the case of two-dimensional triangular meshes. The quality of an element, Q, is obtained as: Q=4√3a h2 1+h2 2+h2 3 , where ais the area of the triangle, and hi, i = 1,2,3,are the three edges. The mean quality of the mesh shown in Figure A.2 is Q= 0.98 with a typical deviation of 0.03. Only eight elements out of the total (0.04%) present Q < 0.8 and one element shows Q < 0.5. Isabel Gonz´ alez Crespo 0123 ×104 0 0.5 1 1.5 2 Number of nodes Error 0123 ×104 0 100 200 300 400 500 600 700 Number of nodes Execution time (s) Figure A.1: Convergence study of the mesh to solve Problem 2.4. The left panel presents absolute error measurements for different numbers of nodes compared to the finest computed mesh (∼60 000 nodes). The right panel shows the respective execution times obtained from running the code in Listing 2.5. Figure A.2: Example of a mesh to apply the FEM to solve the oxygenation problem on the domain shown in Figure 2.1. 132 Appendix B Formal analysis of the response models This appendix presents some aspects of the formal analysis of the response models to RIT and FLASH-RT, respectively given in Sections 3.2.6 and 4.4. B.1 RIT model Using the notation from Section 2.1, the IVP associated to the RIT response model given in Section 3.2.6 can be rewritten as:          dx(t) dt =f(x(t),x(t−τA),x(t−τT)),∀t=ti, ∆x(ti) = Ii(x(ti)),∀t=ti, x(t) = x0,∀t∈[−max{τA, τT},0], (B.1) being x1(t) = C(t), x2(t) = Cd(t), x3(t) = Ta(t), x4(t) = ˆ Ta(t), x5(t) = ˆ Tb(t), x6(t) = ˆ T(t), x7(t) = ˆ A(t), x8(t) = c4(t), and x9(t) = p1(t). As the initial condition is given by a continuous function and fis of class C∞in {x∈R9:x1>0, x2>0}× R9×R9(and particularly of class C1, thus fis locally Lipschitz continuous), then Theorems 1and 2prove that there exist a local solution to the above IVP. Isabel Gonz´ alez Crespo B.2 FLASH-RT model Using the notation from Section 2.1, the IVP associated to the FLASH-RT response model given in Section 4.4 can be rewritten as:          dx(t) dt =f(x(t)),∀t=ti, ∆x(ti) = Ii(x(ti)),∀t=ti, x(0) = x0, (B.2) being x1(t) = C(t), x2(t) = Cd(t), and: x=x1 x2,f=   λ1−x1+x2 Kx1 −λx1+x2 Kx2−ϕx2   ,∆x=(SF −1)x1 (1 −SF)x1,and x(0) = x0,1 x0,2. The function fis at least of class C1in R2, thus local Lipschitz continuity is guaranteed by the mean value theorem. Subsequently, Theorem 1proves that there exist b > 0 (for t0= 0) and a unique local solution x: (0, b)→R2to the above IVP. Notice that, as mentioned in Section 4.4, for this study the initial condition x0,2 was set to 0. Thus, in the absence of impulses, x2≡0, and the IVP is transformed into the following:      dx(t) dt =λ1−x1(t) Kx1(t),∀t > 0 x1(0) = x0,1. (B.3) The analytical global solution to the above IVP is: •If x0,1= 0, then x1≡0. •If x0,1=K, then x1≡K. •If x0,1∈(0, K), then x1=x0,1K x0,1+ (K−x0,1) exp(−λt). In the particular case that there is a single impulse at t= 0, which is the case of study in Section 4.6, the IVP (B.2) is equivalent to the following:                    dx1(t) dt =λ1−x1+x2 Kx1,∀t > 0 dx2(t) dt =−λx1+x2 Kx2−ϕx2,∀t > 0 x1(0) = SFx0,1, x2(0) = (1 −SF)x0,1, (B.4) 134 Appendix B. Formal analysis of the response models where SF ∈(0,1], thus 0 ≤x1(0) ≤x0,1and 0 ≤x2(0) < x0,1. The critical points of the above system have been studied and the respective phase portrait is shown in Figure B.1, having that: •(0,0) is a saddle. •(K, 0) is a sink. •0,−ϕK λis a source. In the context of this thesis, only positive solutions are biologically reasonable. Besides, if x1= 0 and x2= 0 for any t, it would mean that the tumor does not exist, and the presented model was not designed to describe tumor total remission. Thus, (K, 0) is the only critical point of interest. In view of the phase portrait, if the initial condition is in the first quadrant, the solution of the system tends to the sink (K, 0), as expected. Moreover, the critical point being a sink guarantees that the solutions to the system are global, i.e., they are defined in the interval (0,∞), assuming that the initial condition is near enough to (K, 0). −5−4−3−2−1012345 −5 −4 −3 −2 −1 0 1 2 3 4 5 Figure B.1: Phase portrait of the system (B.4) using λ= 1, ϕ= 1, and K= 3. The critical points are shown in red. The IVP (B.4) is a particular case of a generalized Volterra system, the dynamical behavior of which has been studied in the literature [147,148]. 135 Appendix C Forward Euler algorithm stability Given the complexity of the presented models, which involve ODEs,DDEs and IDEs, a simple forward Euler algorithm, described in Section 2.1.3, was implemented to obtain numerical solutions. This method may have stability problems for stiff differential equations, but in general, it is not possible to know a priori whether a given equation presents stiffness [90]. This appendix presents convergence studies for the numerical solution of both the RIT and FLASH-RT response models. C.1 RIT model Figure C.1 shows the evolution of tumor volumes (obtained with best-fitting parameters reported in Table D.4) in the control group (left), and RIT group with 8 Gy×3 and αCTLA4 at days 2, 5 and 8 (right), for different time steps, ranging from ∆t= 1 days to ∆t= 0.001 days. The solution for the lowest time step is considered as a reference and the best approximation to the exact solution. Figure C.2 presents the absolute error made when calculating the volume with each time step against the reference solution at day 10. The error shows a linear dependence on the discretization time step, as expected for the Euler method in well-conditioned problems. ∆t= 0.05 days was selected as a good trade-off between reducing time execution and obtaining accurate approximations of the exact solution (especially when considering the uncertainties of the experimental data). Isabel Gonz´ alez Crespo 0 5 10 15 20 0 200 400 600 800 Time post-start of treatment (days) Volume (mm3) Control ∆t= 1 ∆t= 0.5 ∆t= 0.1 ∆t= 0.05 ∆t= 0.01 ∆t= 0.005 ∆t= 0.001 0 5 10 15 20 0 10 20 30 40 50 Time post-start of treatment (days) RIT (αCTLA4 (2 5 8) & 8 Gy×3) Figure C.1: Tumor volume evolution with different time steps, ∆t(days), using the parameters summarized in Table D.4. In the left panel, curves for ∆t∈ {0.05,0.01,0.005,0.001}overlap, while in the right panel, curves for ∆t∈ {0.01,0.005,0.001}overlap. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 Time step, ∆t(days) Error Control 0 0.2 0.4 0.6 0.8 1 0 5 10 15 Time step, ∆t(days) RIT (αCTLA4 (2 5 8) & 8 Gy×3) Figure C.2: Absolute error when calculating the volume (at day 10 post-start of treatment) with different time steps (asterisks), and linear fits of the obtained values (solid line). The solution with ∆t= 0.001 days served as a reference to calculate the global error. 138 Appendix C. Forward Euler algorithm stability C.2 FLASH-RT model Figure C.3 shows the evolution of tumor volumes (obtained with best-fitting parameters reported in the first column of Table E.1) in the control group (left), and FLASHRT group (12 Gy ×1) (right) for different time steps, ranging from ∆t= 1 days to ∆t= 0.001 days. The solution for the lowest time step is considered as a reference and the best approximation to the exact solution. 0 5 10 15 20 25 0 500 1000 1500 2000 2500 Time post-start of treatment (days) Volume (mm3) Control ∆t= 1 ∆t= 0.5 ∆t= 0.1 ∆t= 0.05 ∆t= 0.01 ∆t= 0.005 ∆t= 0.001 0 5 10 15 20 25 0 400 800 1200 Time post-start of treatment (days) FLASH-RT (12 Gy×1) Figure C.3: Tumor volume evolution with different time steps, ∆t(days), using the parameters summarized in Table E.1). In both panels, curves for ∆t∈ {0.05,0.01,0.005,0.001}overlap. Figure C.4 presents the absolute error made when calculating the volume with each time step against the reference solution at day 10. The error shows a linear dependence on the discretization time step, as expected for the Euler method in wellconditioned problems. ∆t= 0.1 days was selected as a good trade-off between reducing time execution and obtaining accurate approximations of the exact solution (especially when considering the uncertainties of the experimental data). 139