Stochastic modelling of cellular populations: Effects of latency and

feedback

Daniel Sánchez Taltavull

ADVERTIMENT. La consulta d’aquesta tesi queda condicionada a l’acceptació de les següents condicions d'ús: La difusió d’aquesta tesi per mitjà del servei TDX (www.tdx.cat) i a través del Dipòsit Digital de la UB (diposit.ub.edu) ha estat autoritzada pels titulars dels drets de propietat intel·lectual únicament per a usos privats emmarcats en activitats d’investigació i docència. No s’autoritza la seva reproducció amb finalitats de lucre ni la seva difusió i posada a disposició des d’un lloc aliè al servei TDX ni al Dipòsit Digital de la UB. No s’autoritza la presentació del seu contingut en una finestra o marc aliè a TDX o al Dipòsit Digital de la UB (framing). Aquesta reserva de drets afecta tant al resum de presentació de la tesi com als seus continguts. En la utilització o cita de parts de la tesi és obligat indicar el nom de la persona autora. ADVERTENCIA. La consulta de esta tesis queda condicionada a la aceptación de las siguientes condiciones de uso: La difusión de esta tesis por medio del servicio TDR (www.tdx.cat) y a través del Repositorio Digital de la UB (diposit.ub.edu) ha sido autorizada por los titulares de los derechos de propiedad intelectual únicamente para usos privados enmarcados en actividades de investigación y docencia. No se autoriza su reproducción con finalidades de lucro ni su difusión y puesta a disposición desde un sitio ajeno al servicio TDR o al Repositorio Digital de la UB. No se autoriza la presentación de su contenido en una ventana o marco ajeno a TDR o al Repositorio Digital de la UB (framing). Esta reserva de derechos afecta tanto al resumen de presentación de la tesis como a sus contenidos. En la utilización o cita de partes de la tesis es obligado indicar el nombre de la persona autora. WARNING. On having consulted this thesis you’re accepting the following use conditions: Spreading this thesis by the TDX (www.tdx.cat) service and by the UB Digital Repository (diposit.ub.edu) has been authorized by the titular of the intellectual property rights only for private uses placed in investigation and teaching activities. Reproduction with lucrative aims is not authorized nor its spreading and availability from a site foreign to the TDX service or to the UB Digital Repository. Introducing its content in a window or frame foreign to the TDX service or to the UB Digital Repository is not authorized (framing). Those rights affect to the presentation summary of the thesis as well as to its contents. In the using or citation of parts of the thesis it’s obliged to indicate the name of the author.

Stochastic modelling of cellularpopulations: Effects of latency and

feedback

Daniel Sanchez TaltavullAdvisor: Tomas Alarcon

Tutor: Alex Haro

tesi inscrita al programa de doctorat en Matematiquesper optar al tıtol de doctor

Universitat de BarcelonaFacultat de Matematiques

Departament de Matematica Aplicada i Analisi

October 9, 2014

Acknowledgements

Aquesta ha estat un llarg camı ple de pujades i no hauria estat possible sense tota lagent que he tingut al voltant, ja sigui donant-me consells o simplement fent una birra ambmi quan ho necessitava.

Tinc que donar les gracies al meu director de tesi, Tomas Alarcon, per tota la paciencia,tot l’esforc i la gran feina feta en dirigir-me la tesi i alhora dirigir tot el grup, i en ajudar-me cada cop que li he demanat o ha vist que el necessitava. Quan vaig prendre la decisiode venir al CRM, sense saber gairebe res d’ell, tenia por de no haver comes un error, peroara se que venir aquı amb ell ha estat una de les millors decisions de la meva vida, i si hohagues de tornar a fer tornaria a triar igual.

Tambe he d’agraır al meu pare i la meva mare tot el que han fet per mi i per fer-mecostat amb totes o casi totes les meves decisions i ajudar-me a ser qui soc ara. Tambedonar-li les gracies al meu germa gran Alex per ser, a part d’un germa, un amic i ser allasempre que l’he necessitat i tambe per fer-me costat en aquesta aventura. Gracies tambea s’avia Juanita, s’avia Tonia, a en Toni, na Sili i en Victor (que arribara on vulgui almon de la poesia) per tot el que feu per mi, per ser alla i per estimar-me i ajudar-me enlo que es pot. Sense vosaltres aixo no seria el mateix.

I al Centre de Recerca Matematica pel financament i totes les facilitats que m’hanproporcionat perque pogues fer la tesi. Gracies per confiar en mi.

I com no, als membres del meu grup. A Pily por ser como un atardecer de otono,por haberme animado cada dıa que estuvo ahı y por lo bien que nos lo pasamos cadavez que quedabamos para cotillear. A l’Esther per no tirar-nos una galleda d’aigua (olejia) el dia que vam anar a cantar al seu portal. A Roberto por ser una de las pocaspersonas dispuestas a bajar las escaleras para venir a verme y hacerme companıa, y porpasar de los casinos. A Ivon porque, las pocas veces que sale del bunker en que se esconde,consigue sacarme una sonrisa. Y a las dos nuevas incorporaciones: Juan, que aun no leconozco mucho, pero que por ahora nos lo estamos pasando bien y Elisa, que encaja mejor4 parrafos mas abajo.

Agraesc molt sincerament, efusivament i de tot cor a nes meu gran amic Andreu quehagi cregut durant tants anys que π = 3, y estoy muy agradecido al chico de Movistarque me ayudo a elegir un buen movil. Aram gracias por dejarme siempre un sofa dondedormir y cuidarnos siempre. A n’Alex li agraesc que no hagi decidit comencar a tocar esviolı fins as final de sa meua tesi, ja que m’hauria provocat danys cerebrals irreversibles.Gemma, merci per ser sa millor amiga que es pot tenir, i merci per lo bona que esta sateua germana petita, lo qual me recorda que estic molt agraıt a na Maria Galeta per lobona que esta. A en Mati per no pigar-me mai. A n’Ester per totes ses ties que mus hapresentat a jo i n’Andreu. A en Sus per ses lans party que vam montar. A en Nel perhaver estat as costat meu des de fa 23 anys. a n’Afryca per fer-me colacaos es mati i perdonar-me una cullerada d’oli d’oliva quan vaig tornar borratxo de Son Bou. a na Maiteper clavar-me una forquilla i mossegar-me. A na Falsa Vegana per menjar carn. Christian,siento mucho haberte abandonado en ese pueblo... y siento que nos llevasemos tu movil,cartera, dinero y llaves. a Tere por no pedir una orden de alejamiento. A en Jordi Boschper totes ses juergues mıtiques. A na Sara Velez per tot lo que no sap sobre es porter deAlla de sempre. A n’Andrea per ser sa meua consciencia i no deixar-me fer res dolent. A

en Francis si. A en Xoric per estar amb jo durant tot es Diablo 3. a en Jake i en Kid persalvar la tierra media destruint el anillo unico.

Molta gent que m’ha acompanyat en aquest camı han estat els amics que he fet aBarcelona, aquest va per ells: Mariano, ets es millor company de pis que es pot tenir,lo unic que te dire es: clucs. Lluis, encaixaries a nes parraf de dalt pero millor posar-teaquı que es on mus hem fet amics, merci per totes ses borratxeres. Cristina que decir deti, te echare taaaanto de menos, pero se que me vendras a ver muy pronto! Y Marta teprometo que en el proximo viaje buscaremos una cancion que puedas cantar tu tambien,pero por favor, pasanos ya las fotos de Berlin. Cesar, sigo sin entender como te puedegustar Bran Stark, pero sabes que te apoyo. Eva Martinez, saps que ets la unica personaque https://www.youtube.com/watch?v=q5V01fIjN-s. Anouk m’encanta lo guay que ets!Miriam Dooooniga, si d’alguna cosa estic segur es de que mai mes tornare a trobar ningucom tu. Gerard gracies, de ve, moltes gracies per aguantar a tothom en aquella casa delocus i mantenir sa calma en tot moment. Laura Silvestre em vaig riure molt quan em vacontestar el whatsapp de Menorca. Aracelis te he perdido la pista del todo pero seguroque te va tan bien como siempre. Laura Marfill m’encanta la d’aventures que vius sempre,estic segur que en algun viatge raru que faci et trobare. Yasser porque un dıa de estos nosvamos a tomar una cerveza. Andres me encanta cuando me preparas la cena.

Otro grupo de amigo que he hecho son los E.coli & friends. David m’encanta quesempre estiguis de festa fins al final, trobo a faltar que trenquis portes i espero que not’obris el cap per California... Maria graaaaaacies per les galetes que em vas fer, en vullmes. Edgar merci pel deck de hunter, m’ho vaig passar molt be amb ell. Joan merci peltota la ginebra que ens hem begut. Roser em preguntava... ens hem vist algun cop dedia? Alba tots els records que tinc de tu son molt borrosos ;). Gabriel graaaacias por eseguacamole. Y Maggie gracias por ser la mejor! Me encanta que aunque te hayas ido tanlejos sigamos siendo igual de amigos (o mas) y perdon por lo que voy a poner ahora... YRebecca, me encanta que tengas un leon tatuado y me apetece darte un mordisco.

Mi tesis no habrıa sido la misma de no haber ido a la II Summer School del GEFENOLen Benasque (sitio donde deberıan organizarse todas las summer schools del Gefenol, nose porque se montan en otros sitios ultimamente...), donde conocı a mi grupo favorito decientıficos. Uri desde el moment que vas dir “se me ha ido de las manos” hem viscutinnumerables aventures que no canviaria per res del mon. Edgar, ¿cuando fue la ultimavez que salimos un sabado? creo que el Martes. Edoardo El Sabio, eres el protagonista delmejor momento que ha vivido nuestro grupo, siempre recordare como Juan, Edgar y yo temirabamos y vas ahı... y luego te das la vuelta y vuelves andando sin mirar atras... madremia... Y ahı con la boca abierta dos minutos y tu: “voy a por otro cubata”. Svetozarel vallekano GAAS! GAAAAAAS! A Tipo De Incognito porque eres la persona con quienmas he bebido durante la tesis... A mi higado le caes tan bien como yo a tu novia... UlaSagarra no puc posar aquı el motiu pel qual et tinc en un pedestal aixı que et guardo pelparraf de les JIPI. A Virginia por tener unos ojos tan bonitos. A Paula Villa Maravilla porser un sol, y por dejarnos quedarnos en casa de sus abuelos aka El Gran Hotel Paula Villa,y por no enfadarse cuando inundamos la casa de vecino, pasase lo que pasase: Merecio lapena. A Ruben Perez Carrasco por esa mirada que me echo en la discoteca esa... Magia!Unicornios! Que gran noche. Ana Kalt por la mejor tortilla de patatas que he probadoen mi vida. A Toni por estar siempre dispuesta a irse de fiesta. A Angel Maximilia perl’Angel. A Nicole por esas conversaciones sobre Hannah Montana. A Rafa por todas esasfiestas, he sido incapaz de encontrar el antro ese al que nos llevaste. A Alessio que aun no

me ha pagado lo de la apuesta... Alessio Paga el cubata segundo aviso. A Nacho por nomatarme todas las veces que le desperte la segunda semana en Benasque... A Eustaquiopor esas cenas en el japones de Passatge de Marimon. A Trapote per ser una taza. Ya Elisa Beltran por ensenarme a valorar la verdadera belleza de las playas. Lo cual merecuerda... Gracias Massimiliano por ensenarme como NO funciona un GPS y graciasFran por los zapatos. Y a nuestros amigos de Dawn Meats por suministrarnos siemprealimentos de la mejor calidad. Y por ultimo, aunque no estuvieron en este congreso, aquıencajan bien: Joaquin aka Marcelino Adbekunkus eres mi heroe, sere el primero en venira verte y Raquel gracias por pagarme ese vodka con limon, pero sabes que yo no bebo.

Durant la carrera de Matematiques he conegut a molta gent i molts d’ells no podenfaltar en els agraıments d’aquesta tesi. David Martı moltes gracies per tot lo que vas ferquan em van ingressar a l’hospital i per tot el suport que em vas donar, encara em ricquan penso que abans de que m’operessin et vaig dir: “tranquil no es res” i tu “pero sid’aixo et pots morir!!”. Narcis ets un crack, merci per totes les converses censurades perwhatsapp... tic tac tic tac. Arturo, gracies per totes les mates que m’has ensenyat aquestultim any, sense tu mai hauria dominat l’aproximacio semi-classica, i aviam quan et tornesa venir a fer unes birres. Jordi, encara que fas veure que ets dolent se que ets un granamic, cal ser-ho per no anar a una esquiada que ja has pagat quan un colega ho necessita.Dani, gracies per totes les classes de gnuplot i tex, i per haver estat “al pie del canon”constantment. Meri que be m’ho vaig passar amb tu i que gran va ser aquella setmana alcongres dels espais de moduli. Ari merci pels snapchats brutals que ens envies. Juan VonMalder sempre m’animes a superar-me a mi mateix. Marta Canadell moltes gracies pelvodka de Russia. Ramon, vas ser el primer amic que vaig fer a mates i m’encanten totesles aventures que hem viscut. Anna ets la millor, la meva matematica preferida. InakiPuigdollers m’ho he passat molt be en totes aquelles nits de jocs. Marcel ets el matematicmes fiestero que he tingut l’honor de coneixer, aviam quan tornes i ens bevem un altrebar. Marta Palau fas que m’agradin els chupitos de tequila. Roca merci per tot aquellvodka que ens hem fotut. Roc em fa molta gracia la frase per lligar que em vas ensenyar,algun dia la fare servir.

Tambien quiero agradecerle a Alex Haro todo lo que ha hecho por mi estos anos,sobretodo gracias por todos los papeles que has firmado y por hacer que me decidiese porla matematica aplicada. I tambe gracies als professors que mes m’han influenciat, JuanLuis, vas fer que em decidis a comencar aquest camı de les matematiques. I gracies CarlesSimo, que ets qui m’ha convertit en un matematic tant bo com soc.

I tinc que dedicar un paragraf als JIPIs (que no hippys, be n’hi ha algun que si,pero aixo es una altre historia), l’Anna Alemanai per totes aquelles birres a la TerretaValenciana i totes aquelles falses acusacions sobre que organitzo congressos per coneixerties. A la Maria Van Der Waals (aka Mireia Martı) per ser tant bona dic... per estar tantbona, no espera! ja ho havia dit be al principi: per ser tant bona. I a l’Anna May per sertant maca i tant happy i tant de tot, et vull adoptar com germana petita. I a l’Ula perlo genials que son les festes a casa seva i per la quantitat indecent d’alcohol que me begutamb ell. I al Martı per lo be que li queda la barba, la llista exhaustiva de les qualitatsque ha de tenir una noia i per les mossegades que la seva xicota, la Vero, va donar alMario o al Pol, no recordo a qui dels dos va ser. A Daniel Reta por esa serenidad que nostransmite siempre. I al Joan Camunas per aquelles xerrades sobre matematiques, ell jasap a quina matematica em refereixo ;). I al Rebled per les nits de karaoke i sauna i lesxerrades sobre com ser un gran amant. I a la Teresita de Jesus, al Chileno i a la filosofa

for fruitful discussions.Un altre congres que ha estat molt important per mi, per la gent que vaig coneixer

va ser el FisEs. Alla vaig coneixer al Pol, que espero que aquest any la convidi... I a laSara la Jabata, que apart de demostrar se una persona molt flexible... be... medallon,chuleta, fresa, frambuesa, encara espero una trucada teva. Celia me encantaron todas esasconversaciones sobre monstruos. I al David Palau, que estava per alla parlant de gossos imariners. I Mario, que es un fiestero com deu mana, i quan no em censura els agraımentsencara es mes guay. I a la Lorena, que encara em ric quan penso en la trucada del divendresa les tantes. Tambien estaban por ahı 4 chicas muy simpaticas que me ayudaron muchocuando estaba algo perdido, Patricia la de los ojos verdes, Vanessa Cobi, y las otras dosque nunca supe como se llamaban... Romualda y Eustaquia por ejemplo. Luino cada vezque dicen la palabra depravado pienso en ti. Y Jordi de Granada a ti te dire ;) ;).

I evidentment a totes les persones que he trobat al CRM. Michelle I miss you so much!!Fran Martinez, merci per tota la calma que em transmeties sempre. Marina, moltes graciesper ajudar-me a construir la urna amb l’unicorn. Carles a tu t’estic molt agraıt per allode la tassa aquella que NO era robada. Rosalba ets una persona molt especial... pero novull dir rara, que tambe, ets molt rara, pero aquı amb ningu m’ho passo tant be com ambtu. Bernat, sempre que penso en tu recordo el dia de l’incendi a l’orfanat. Isabel I jatinc ganes de saber on acabarem el proxim cop que sortim de festa. A Luis le agradezcotodas las conversaciones sobre las causas justas que apoya. Jordi gracies per fer que totfuncioni. Dani Balague et poso en aquı tambe! gracies per ajudar-me en un dels momentsque mes ho he necessitat, espero veure’t mes ara que estarem mes aprop. A Alvaro quesiempre me apoya cuando me cambio de despacho. A Tim que me hizo ganar un dineroextra con apuestas deportivas. A Andrei que monto uno de los congresos en que mejor melo he pasado. Al Joaquim, que sempre m’ha recolzat en tot lo que he volgut fer. Guillemgracies per ensenyar-me l’starlink. Mihail el poco tiempo que estuviste aquı tuve una gransuerte. Y gracias Mikel Brun por ensenarme valores, el dıa que le conocı vi un ejemplo derectitud.

Y dio la casualidad que a la primera summer school que fuı en la tesis fue el Biomaten Granada, y la ultima (tres anos mas tarde) el Biomat en Granada. Donde estaba Nico,que este si que es un fiestero como Dios manda. Y Paloma que es malvada (y lo digocomo algo bueno). Y Maite que tenıa una sonrisa contagiosa. Y Pablo, el otro Pablo, y elotro Pablo. Y el chico ese que no se como se llamaba con la camiseta de berlin que salesin cabeza en las fotos. Y tambien estaba Isabel II, que no salio... A ver si ahora que tetengo mas cerca te veo algun dıa de fiesta, gracias por ayudarme a organizar el seminario!Tambien estoy muy agradecido a Rica Schnei que le tiraba cacahuetes a Nico mientrasintentaba ligar. Y tambien a Ainize por ser tan fiestera.

A Soraya Garcia, Daniel Inglorion, Alvaro Riera Bable, Marc Parti, al vendedor degemas de Geffen (por la cantidad de negocios que acepto), al sacerdote de Elune por ser elelfo mas feo de la historia. Al agente imperial Cipher 9, a la contrabandista Pancracia. AlBounty Hunter Firex por todos esos partidos de Huttball y por su estrategia en Ironman.Al legıtimo rey de los siete reinos y protector del reino Joffrey Baratheon. Al High Priest ypalabra de Dios en Rune-Midgars, Walter Sullivan. Y al verdadero senor de Rune-Midgars,el Wizard Wallace. Special thanks to Cave Johnson and his secretary, the sweet Caroline,authors of the project GLaDOS that inspired my scientific career.

Contents

0 Resum en catala 1

0.1 Motivacio i rerefons biologic . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

0.1.1 Latencia en poblacions de cel·lules . . . . . . . . . . . . . . . . . . . 1

0.1.2 Dinamica estocastica de les cascades de diferenciacio regulades retroac-tivament . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

0.1.3 Dinamica estocastica i extincio del VIH-1 en pacients sota terapiaanti-retroviral potent . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

0.2 Marc matematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

0.3 Objectius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1 Introduction 9

1.1 Motivation and biological background . . . . . . . . . . . . . . . . . . . . . 9

1.1.1 Latency in cell population dynamics . . . . . . . . . . . . . . . . . . 9

1.1.2 Stochastic dynamics of differentiation cascades with regulatory feed-back . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.1.3 Stochastic dynamics and extinction of HIV-1 in patients under po-tent anti-retroviral therapy . . . . . . . . . . . . . . . . . . . . . . . 12

1.2 Mathematical framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.3 Aims and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Robustness of differentiation cascades with symmetric stem cell division 17

2.1 Biological Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2 Stochastic modelling of a serial differentiation cascade with symmetric SCdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.3 Extinction induced by symmetric SC division . . . . . . . . . . . . . . . . . 23

2.3.1 Delay-induced extinction in serial differentiation cascades with sym-metric SC division . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.3.2 Delayed cytokine scenario . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4 Robustness of populations with asymmetric SC division . . . . . . . . . . . 27

2.4.1 Stochastic modelling of a differentiation cascade with asymmetricstem cell division . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4.2 Competition between populations . . . . . . . . . . . . . . . . . . . . 30

2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3 Antigen-stimulation induced eradication of the HIV-1 infection in pa-tients under highly active anti-retroviral therapy 39

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.2 Stochastic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.3 Asymptotic analysis and numerical methods . . . . . . . . . . . . . . . . . . 44

3.3.1 Semi-classical approximation . . . . . . . . . . . . . . . . . . . . . . 44

3.3.2 Quasi-steady state approximation . . . . . . . . . . . . . . . . . . . . 46

3.3.3 Heteroclinic orbits and extinction time . . . . . . . . . . . . . . . . . 50

3.3.4 Lower Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.3.5 Multi-scale stochastic simulations . . . . . . . . . . . . . . . . . . . . 53

3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.4.1 Subcritical case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3.4.2 Supercritical case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.4.3 Side effects: viral blips . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4 Stochastic modelling of viral blips in HIV-1-infected patients: Effects ofinhom*ogeneous density fluctuations 61

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.2 Model formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.2.1 Compartmental Model . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.2.2 In silico blood sample analysis model . . . . . . . . . . . . . . . . . 68

4.2.3 Effects of continuous versus burst production of virions . . . . . . . 72

4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.3.1 Effect of density fluctuations on blips statistics . . . . . . . . . . . . 72

4.3.2 Effect of burst production of virions . . . . . . . . . . . . . . . . . . 73

4.3.3 Effect of post-extraction handling time . . . . . . . . . . . . . . . . . 75

4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5 Conclusions 79

5.1 Stochastic dynamics of differentiation cascades with regulatory feed-back . . 79

5.2 Stochastic dynamics and extinction of HIV-1 in patients under potent anti-retroviral therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.3 Effects of density fluctuations on the statistics of viral blips . . . . . . . . . 82

5.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Bibliography 83

A Stochastic modelling and the Gillespie algorithm 85

A.1 Basic definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

A.1.1 The Chapman-Kolmogorov Equation . . . . . . . . . . . . . . . . . . 87

A.1.2 The Master equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

A.2 classic examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

A.2.1 The Moran process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

A.2.2 The stochastic Logistic Equation . . . . . . . . . . . . . . . . . . . . 88

A.2.3 Equilibrium and absorbing States . . . . . . . . . . . . . . . . . . . . 89

A.3 Monte Carlo Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

A.4 Gillespie Stochastic Simulation Algorithm . . . . . . . . . . . . . . . . . . . 91

B Semi-Classical Approach 93B.1 Metastablility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

C Computation of Heteroclinic Connections for Hamiltonian System 101C.1 The Taylor Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101C.2 Automatic Differentiation Formulas . . . . . . . . . . . . . . . . . . . . . . . 101C.3 Heteroclinic connections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

D Multi-scale Stochastic Simulations 105D.1 Cao-Gillespie-Petzold Method . . . . . . . . . . . . . . . . . . . . . . . . . . 105

D.1.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105D.1.2 Computation of Partial equilibrium . . . . . . . . . . . . . . . . . . . 106D.1.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Chapter 0

Resum en catala

En aquesta tesi doctoral es tracten alguns problemes sobre fluctuacions en poblacionsamb una estructura complexa, en particular aquelles caracteritzades per l’existencia d’unestat latent. El nostre objectiu es entendre millor algunes questions que sorgeixen en ladinamica de les poblacions de cel·lules amb estructura jerarquica, com ara les que es trobenen els organismes multicel·lulars caracteritzats per la diferenciacio cel·lular, i la dinamicadel VIH-1 en pacients infectats que es troben amb terapia antiretroviral. Els dos sistemesestan caracteritzats per la presencia d’una poblacio en estat latent, es a dir, cel·lules mareen el primer cas, i cel·lules infectades en estat latent en el segon.

0.1 Motivacio i rerefons biologic

Els problemes biologics han estat estudiat pels cientıfics amb 4 models diferents: In vivo,In vitro, Data Driven i mechanistic models. Encara que aquests ultims son els que donenun millor enteniment sobre la dinamica del proces que s’estudia, han estat els que menysexit han tingut en biologia.

Els mechanistic models s’han utilitzat tradicionalment per matematics i fısics per estu-diar problemes biologics fent servir aproximacions continues i deterministes. No obstant,en alguns casos, aquests models poden no ser molt correctes. Si les poblacions son petitesno podem assumir que un model de camp mitja sigui una bona aproximacio de la realitat.Aixo es important en molts problemes biologics, on sovint, un es troba poblacions petites.

S’ha provat que el soroll intrınsec, es a dir, fluctuacions en el sistema degut a un tamanypetit, pot tenir una forta influencia en el comportament d’un sistema dinamic [14, 85, 42].Presentem una introduccio als processos estocastics a l’Apendix A, concretament tractasobre modelitzacio estocastica i l’algoritme de Gillespie de simulacio estocastica, que serala nostra principal eina.

0.1.1 Latencia en poblacions de cel·lules

Pels proposits d’aquesta tesi, la latencia es defineix com un estat de la cel·lula en que lamajoria de les seves funcions, prominent proliferacio, estan desactivades. Aquest mecan-isme sovint sorgeix en entorns hostils per sobreviure. Un exemple d’aquest comportaments’ha analitzat en [3], on s’ha provat que la latencia pot actuar com un mecanisme que per-met a una poblacio maligne eludir l’accio d’un medicament, fent que les cel·lules d’aquestapoblacio tinguin la capacitat d’entrar en un estat de repos, en el qual son immunes a

1

2 Resum en Catala

l’efecte de la droga. Al entrar en aquest estat, una petita poblacio latent es capac desobreviure, i, finalment, tornar a creixer un cop que el medicament ha deixat d’aplicar-se.Un cas particular en que aquest model es rellevant, es la resistencia a la quimioterapia o ra-dioterapia en els tumors amb regions amb hipoxia (falta d’oxigen). Les cel·lules hipoxiqueses sotmeten a una gran reduccio de la seva taxa de proliferacio, que les fa menys sensiblesa aquests tractaments que les cel·lules de cicle rapid, proporcionant aixı un reservori decel·lules immunes, que, quan es para la terapia, fan que aquesta poblacio torni a creixer.A mes d’escapament a la terapia, els fenomens de latencia han provat tenir un paper enel cancer per eludir els controls homeostatics del teixit normal. Roesch et al. [108] hancaracteritzat un cicle lent (latent) d’una subpoblacio dins de la poblacio del melanoma,que te una proliferacio de cel·lules molt mes rapida necessaria per al creixement tumoral.

Un altre exemple on la latencia ha provat ser crucial, es per entendre certs aspectesdel VIH i la seva resposta a la terapia. El VIH es controla amb molta eficacia ambterapia antiretroviral. No obstant, aquestes terapies no suprimeixen completament elvirus. Encara que eviten que els virus infectin cel·lules T, estudis a llarg termini de lacarrega viral en pacients tractats revela que un nivell baix del virus (per sota dels lımitsde deteccio estandard) persisteix molt de temps [110, 123]. Aquests resultats suggereixenque la infeccio persisteix en els pacients tractats en forma de compartiment insensible a ladroga. Moltes hipotesis han estat formulades sobre la naturalesa d’aquest compartiment,pero la mes consolidada es la que ho explica amb una poblacio de cel·lules infectades enestat latent. Aquest compartiment cel·lular consisteix en unes cel·lules (probablement dememoria) que proliferen lentament i guarden el virus pero no el repliquen. Com aquestescel·lules no repliquen el virus no son afectades per la terapia antiretroviral. A mes, despresde l’estimulacio amb antıgens especıfics, passen a un estat actiu i reprenen la produccio devirus. Aquest cicle mante aquest baix nivell de viremia. La presencia d’aquesta infecciolatent ha estat reconeguda com la principal barrera per a la completa eradicacio de lainfeccio, i la recerca de mitjans per combatre-la s’ha convertit en una molt activa lınia derecerca [121, 113, 65, 66].

El fenomen de la latencia tambe es part de la regulacio fisiologica normal dels teixits.En particular, les cel·lules mare, que son els components essencials dels teixits dels organ-ismes multicel·lulars, passen la major part del seu cicle de vida en un estat latent, es a dir,sense proliferar. Nomes despres de rebre els senyals apropiats, les cel·lules mare s’activeni la proliferen. Aquestes senyals son secretades en resposta a una disminucio en el nombrede cel·lules madures (totalment diferenciades) que estan al final de la corresponent cadenade diferenciacio: quan el nombre de cel·lules madures disminueix, la proliferacio cel·lulars’activa per tal de reposar les cel·lules madures. Encara que aquest mecanisme es part dela regulacio normal del teixit, pot ser subvertit per patologies com el cancer que segrestenles propietats de les cel·lules mare.

L’objectiu general d’aquesta tesi es utilitzar models estocastics, juntament amb tecniquesasimptotiques i metodes numerics per a comprendre millor els mecanismes basics que in-tervenen en la dinamica de les poblacions del cel·lules latencia retroalimentada amb lafinalitat de proposar estrategies de control que ajudin a formular millor, els enfocamentsterapeutics per situacions patologiques que resulten de la des-regulacio de les cascades dediferenciacio i del VIH-1 en estat latent.

Resum en Catala 3

0.1.2 Dinamica estocastica de les cascades de diferenciacio reguladesretroactivament

En el Capıtol 2 analitzem els factors que afecten a la robustesa de les poblacions de cel·lulesamb estructura jerarquica. Els teixits en organismes multicel·lulars com els mamıfers, man-tenen l’homeostasi per mitja de cascades de diferenciacio de cel·lules en serie [4]. Aquestescascades estan compostes per compartiments jerarquicament organitzats de diferents tipusde cel·lules. La pedra angular d’aquest sistema son les cel·lules mare (SCS) que sostenenel teixit produint mes cel·lules mare, una propietat coneguda com auto-renovacio, i lescel·lules diferenciades d’un llinatge especıfic. La cascada de diferenciacio procedeix atraves d’una serie d’etapes intermitges o compartiments de cel·lules amplificadores tran-sitories (TAC) i acaba amb el compartiment de cel·lules madures o totalment diferenciades(MCS).

Hi ha hagut molta discussio sobre la questio de com les SC mantenen tota la cadena apartir de l’auto-renovacio i la diferenciacio [91, 92, 71, 114, 72]. Una estrategia mitjancantla qual les SCS poden realitzat aquestes dues tasques es la divisio cel·lular asimetrica,d’aquesta manera, en la divisio d’una SC, una cel·lula conserva la identitat de la sevamare mentre que la seva germana es diferencia [71, 72].

Encara que la divisio asimetrica es una solucio simple i elegant, hi ha hagut moltacontroversia al voltant del que es percep com un defecte fonamental d’aquest model, esa dir, la divisio asimetrica no permet que el compartiment de les cel·lules mare creixi[91, 92]. En situacions com lesions, on la taxa de renovacio cel·lular ha d’augmentar engran mesura, amb la finalitat de regenerar el teixit afectat, el compartiment de les SC hade creixer. Aixo suggereix que la divisio asimetrica podria no ser la solucio completa.

Hi ha evidencies experimentals que les SC, en alguns teixits, com el sistema nervioscentral [50] i l’epidermis [74], es poden dividir de forma simetrica. La divisio simetrica deSC consisteix en que les dues filles son iguals, o be SC o be cel·lules diferenciades. D’aquestamanera, el nombre de cel·lules esta regulat per les frequencies de divisions simetriques queprodueixen SC o cel·lules diferenciades. Com aquest tipus de divisio permet l’expansio delcompartiment de les SC, apareixen models alternatius [91] on la majoria de les SCs potdividir-se d’ambdues formes, simetricament o asimetricament, l’equilibri entre aquests dostipus de proliferacio esta controlat per senyals ambientals per produir el nombre apropiatde SC i cel·lules diferenciades. Pero s’ha reconegut que, si be la divisio simetrica confereixla capacitat de creixement millorat i una major capacitat de regeneracio, tambe augmentala probabilitat de cancer. Aquesta idea esta recolzada pels resultats d’acord amb els qualsla maquinaria cel·lular involucrada en la divisio asimetrica ha estat conservada, perquete un paper en la supressio de tumors, els gens que promouen la divisio simetrica tambefuncionen com oncogens [91].

El nostre objectiu es analitzar la dinamica estocastica de les cascades de la diferenciaciocel·lular on trobem les dues formes de proliferacio de les SC. Estudiarem les estrategiesoptimes entre divisio simetrica (que millora la capacitat d’adaptacio) i la divisio asimetrica(que dota d’estabilitat).

0.1.3 Dinamica estocastica i extincio del VIH-1 en pacients sota terapiaanti-retroviral potent

El VIH es controla de manera efectiva amb l’administracio de terapies antiretrovirals(ART). Pero a pesar del gran exit de la terapia antiretroviral de gran activacio (HAART),

4 Resum en Catala

aconseguir curar-la, en el sentit d’eliminar completament el virus, no s’ha aconseguit[105]. Analisis quantitatius de l’evolucio temporal de la carrega viral en plasma despresdel tractament HAART suggereix l’existencia de diverses fases en el decaıment induıt peltractament de la carrega viral [102, 69, 111]. Primer s’observa una fase en la que la carregaviral decau exponencialment, en aquesta fase la carrega viral es redueix d’un a dos ordresde magnitud en unes dos setmanes. Aquesta etapa reflecteix que el temps mitja de vida devirus en plasma i de les cel·lules infectades produint virus, els limfocits T CD4+, es moltcurt [57, 124, 83]. Despres d’aquesta primera fase, una segona etapa mes lenta reflecteixla contribucio a la carrega viral de cel·lules infectades amb un major temps mitja de vida,com els macrofa*gs i les cel·lules T CD4+ infectades que presenten una menor taxa dereplicacio viral [99, 58, 126]. Despres d’aquesta segona fase, la carrega viral en plasmanormalment ha caigut per sota del llindar de deteccio dels assaigs clınics estandards (∼ 50copies d’ARN / ml). No obstant, despres d’aquesta segona fase, la HAART no aconsegueixeradicar completament la infeccio. Es produeix una tercera etapa (de l’ordre d’anys [69])en que els nivells residuals (1-5 copies d’ARN / ml) persisteixen en plasma, aixi com enaltres llocs del cos com el sem*n.

La pregunta de quin es la font d’aquesta carrega viral residual ha generat molt debat is’han creat diverses hipotesis de treball. Una d’aquestes hipotesi es basa en la possibilitatde que la terapia HAART no arriba a tot el cos i hi ha uns reservoris on la terapia noafecta [102], els anomenats “drug sanctuaries”, on la infeccio persisteix [67].

Un model alternatiu, que esta molt recolzat pels estudis recents, suggereix que, tot ique la terapia HAART es totalment supressora, hi ha un reservori cel·lular que permeta la infeccio romandre en forma latent [102], i la carrega viral residual es el resultat del’activacio de les cel·lules latents [102]. Aquestes cel·lules latents s’estableixen dins de lapoblacio de les cel·lules T de memoria CD4+ infectades [25, 24] i, per tant, romanen enl’estat de repos en la presencia de la HAART durant perıodes prolongats de temps. Coma consequencia, les cel·lules infectades de forma latent son capaces d’escapar a l’efecte dela droga i al sistema immune a causa del fet que presenten nivells molt baixos de ARNde VIH-1 [111]. Pero ja que les cel·lules amb infeccio latent alliberen virus quan sonestimulades amb l’antigen apropiat, si es retira la HAART la carrega viral torna a creixeri la infeccio VIH-1 reapareix, aixo es consistent amb que la latencia serveix com a moded’escapament dels medicaments [3].

Un cop que la infeccio ha entrat en aquest estat latent, un pot observar episodistransitoris en que la carrega viral puja per sobre dels lımits estandards de deteccio (50copies d’ARNm/ml) durant un breu perıode de temps [33, 104, 54, 95]. Aquests episodises coneixem com blips. L’origen i la rellevancia medica d’aquests blips encara no esta clarai nombroses hipotesis s’han formulat [111].

Degut a la incapacitat de la HAART per acabar amb les cel·lules en estat latent,s’han proposat terapies combinades que suprimeixen el reservori latent [105]. Un camıprometedor en aquesta direccio consisteix en la combinacio de HAART amb agents queactiven especıficament les cel·lules latents, el que les fa susceptibles als atacs de la HAART.El fonament d’aquesta terapia combinada es que si s’augmenta el ritme d’activacio de lescel·lules amb infeccio latent, el numero d’aquestes ha de decaura rapidament i, degut afluctuacions estocastiques, aquestes cel·lules poden ser completament eliminades i, ambelles, tota la infeccio.

Resum en Catala 5

0.2 Marc matematic

En termes generals, els problemes adrecats en aquesta tesi es refereixen a la dinamicaestocastica de poblacions amb diferents tipus de cel·lules. El marc mes adequat per ferfront a aquests problemes es el de dinamica de poblacions. Degut a que el numero decel·lules en estat latent es molt baix, hem de modelar les nostres poblacions en termes deprocessos estocastics, en particular, dels processos de Markov [44]. El proces de Markovper la dinamica de poblacions en formulara en termes de la corresponent Equacio Mestre(Master Equation). Aquestes equacions son molt difıcils de resoldre exactament i enqualsevol cas practic, s’ha de recorrer a metodes aproximats per tractar amb elles.

En alguns casos treballem amb poblacions molt petites, en altres situacions buscaremrare events. Per aquesta rao farem servir l’aproximacio semi-classica, que es una aproxi-macio del tipus WKB [73] per calcular les probabilitats dels rare events.

Hem inclos una explicacio detallada sobre l’aproximacio semi-classica en l’Apendix B.Aquesta aproximacio transforma el nostre problema de treballar amb una equacio mestreen un problema de mecanica classica amb un Hamiltonia amb diversos graus de llibertat,tants com tipus diferents de cel·lules. Calcular les varietats invariants del sistema, que,normalment es el primer pas en l’estudi d’un sistema dinamic, no es facil i una aproximaciolineal no es suficient i requereix tecniques com el metode de la parametritzacio [17, 18] perobtenir expansions a ordres alts. Mes dificultats associades a aquests sistemes apareixendegut a la naturalesa slow-fast dels problemes. A mes, l’estudi del moviment aprop depunts hiperbolics pot ser necessari. Per integrat les equacions del moviment fem servirmetodes de Taylor [64] combinats amb tecniques de diferenciacio automatica [64], el queens permet prendre passos molt llargs i aconseguir mes precisio del que que altres metodesexplıcits, com el Runge-Kutta, ens poden donar. El metode de Taylor no es un integradorsimplectic, es a dir, no preserve l’estructura Hamiltoniana de les equacions, no obstant, aixoes irrellevant, ja que un pot demanar un error local per sota del roundof de la maquina.A l’apendix C expliquem els metodes que hem fet servir per estudiar aquests sistemesHamiltonians.

Encara que hem fet alguns progressos en l’estudi asimptotic d’aquests problemes, de-penem molt dels metodes de simulacio. La nostre principal eina per estudiar les equacionsmestre es l’Algoritme de Gillespie de Simulacio Estocastica (SSA), que es un metode deMonte Carlo basat en una reinterpretacio de l’equacio mestre. A l’Apendix A aquestametodologia esta explicada en detall. el Gillespie SSA requereix un gran temps de com-putacio, en particular quan (i) el sistema considerat es molt gran i (ii) quan la separaciod’escales temporals apareix. En aquest ultim cas, podem fer servir la separacio d’escalestemporals per formular una variant mes eficient del SSA proposat per Cao et al. [20],aquest metode esta explicat en detall a l’Apendix D.

0.3 Objectius

L’objectiu principal d’aquesta tesi doctoral es l’estudi de l’efecte de les fluctuacions enpoblacions acoplades en sistemes biologics, on cel·lules en estat latent juguen un paperimportant. Intentant trobar el significat biologic de la dinamica dels sistemes. Els puntsespecıfics que volem abordar i la organitzacio de la tesi estan explicats a continuacio.

En el Capıtol 2, estudiem el comportament de les poblacions de cel·lules amb estructurajerarquica des de el punt de vista de les propietats d’estabilitat, En particular

6 Resum en Catala

1. Divisio simetrica contra asimetrica en el compartiment de les cel·lules mare. Es-tudiem la robustesa de les poblacions amb estructura jerarquica, depenent de si lescel·lules mare es divideixen simetricament, asimetricament o de les dues maneres.Estudiem com la divisio simetrica afecta a l’estabilitat de la poblacio, ja que aixo teuna gran importancia en la progressio del cancer.

2. La competicio entre dues poblacions amb diferents tipus de divisio de les cel·lulesmare. Aixo es crucial per trobar estrategies optimes que maximitzin la robustesa(supervivencia a llarg termini, resistencia a invasions i habilitat per invadir) depoblacions amb estructura jerarquica.

3. La influencia de parametres com son la duplicacio i el ritme de mort de les cel·lulesmare, el temps de vida mitja de les cel·lules completament diferenciades, la longitudde les cadenes de diferenciacio i les fluctuacions al compartiment de les cel·lules mareen la robustesa i arquitectura optima de les cascades de diferenciacio.

En el Capıtol 3 presentem un model hom*ogeni de combinacio de HAART amb terapiesd’activacio de les cel·lules latents del VIH-1 a la sang. Estem interessats en

1. L’Efecte del ritme d’activacio de les cel·lules latents en el temps mitja de vida de lainfeccio. En particular analitzem si les terapies basades en incrementar aquest ritmeson capaces de suprimir la infeccio en un temps raonable.

2. La importancia de l’eficiencia de les terapies antiretrovirals, incloent els casos lımiten que l’eficacia es del 100%, en la quantitat de carrega viral.

3. La formulacio d’una teoria asimptotica basada en l’aproximacio semi-classica ambaproximacions quasi estacionaries per descriure la dinamica del proces. La precisiod’aquest metode asimptotic es comparat amb simulacions multi-scale proposades pelCao et al. [20].

En el Capıtol 4, estenem el model proposat pel Rong i el Perelson [110] a un modelno hom*ogeni de la dinamica del VIH-1 en el corrent sanguini, considerant que les cel·lulesi els virus no estan distribuıts de manera uniforme en la sang. Els punts especıfics quevolem estudiar son:

1. El mecanisme que fa que apareguin els episodis de viremia per sobre els lımits dedeteccio, coneguts com viral blips. En particular volem investigar si son productede fluctuacions estocastiques degudes a la inhom*ogenietat o un altre mecanisme hade ser considerat.

2. Si l’aparicio dels viral blips esta afectada pels procediments duts a terme en ellaboratori, com el temps d’espera entre les extraccions i les observacions.

3. Si la probabilitat, l’amplitud i la frequencia dels viral blips es veu afectada pelsdiferents possibles tipus de produccio viral, es a dir, continua vs burst.

En el capıtol 5 presentem i discutim els resultats obtinguts, i comparem, quan espossible, amb altres models o amb resultats experimentals, i discutim el treball que deixempel futur.

Resum en Catala 7

Els detalls relatius a questions metodologiques, aixı com una introduccio a la mod-elitzacio estocastica fent servir equacions mestres es donen en els apendixs. Per a aquellsque no estan familiaritzats amb els models basats en equacions mestres, l’autor recomanallegir primer l’apendix A que proporciona la base matematica per entendre el capıtol 2.Els Apendixs B, C i D juntament amb l’Apendix A donen la base matematica necessariaper seguir el capıtol 3 i el capıtol 4.

8 Resum en Catala

Chapter 1

Introduction

This PhD Thesis deals with some issues concerning the fluctuations in complex structuredpopulations in particular those characterised by the existence of a latent state. Our aimis to better understand some questions arising in the dynamics of hierarchically organisedpopulations of cells, such as the ones encountered in long-lived multicellular organismscharacterised by ongoing cell differentiation, and the dynamics of HIV-1 in infected pa-tients under anti-retroviral therapy. Both systems are characterised by the presence of alatent population, namely, stem cells and latently infected cells, respectively. This chapteris devoted to present the necessary biological background and motivation, summarisingprevious work done on the subject. We also give an overview of the mathematical modelsused in this PhD Thesis. Finally, we summarise our aims and objectives.

1.1 Motivation and biological background

Biological problems have been studied by scientists via four different type of models: Invivo, In vitro, Data Driven and Mechanistic Models. Although, Mechanistic Models are theones which gives a better understanding of the dynamics of the process that are studying,has also been the less successful in biology.

Mechanistic Models have been traditionally used by mathematicians and physicists tostudy biological problems using continuous and deterministic approaches. However, insome cases, this models can be not very accurate. If the population is small, we can notcan not assume that the mean-field model will be a good approximation of the reality.This is relevant in many biological problems, where often one finds small populations.

It has been shown that intrinsic noise, i.e. fluctuations to the system due to its smallsize, can have very strong influence on the behaviour of dynamical systems [14, 85, 42]. Anintroduction to the stochastic processes that we use have been included in the AppendixA about stochastic modelling and the Gillespie stochastic simulation algorithm, which willbe the main tool to deal with our problems.

1.1.1 Latency in cell population dynamics

For the purposes of this thesis, latency or quiescence is defined as a cell state such thatmost cell functions, prominently proliferation, are slowed down. This down-regulation ofcellular functions often emerges as a trade-off in exchange for longevity and survival inhostile environments. An example of this behaviour has been analysed in [3], where it

9

10 Introduction

has been shown that latency can act as a an escape mechanism that allows a malignantpopulation to elude the action of a drug, provided that cells of such population have theability to go into a quiescent state in which cells are immune to the drug. By going intoquiescence, a small latent population is able to survive, and eventually re-grow once thedrug has been cleared off. A particular instance in which this model is relevant is resistanceto chemo- or radio-therapy in tumours with hypoxic (ill-oxygenated) regions. Hypoxic cellsundergo radical reduction of their proliferation rate, making them less sensitive to thesetherapies than fast-cycling cells, thus providing a reservoir of immune cells which, uponstopping the therapy, provide a reservoir for which (partial) re-grow occurs. Besides escapefrom therapy, quiescence-like phenomena have been shown to play a role in cancer to eludenot only therapy but the homeostatic controls of normal tissue. Roesch et al. [108] havecharacterised a slow-cycling (latent) sub-population within the much faster proliferatingpopulation of melanoma cell which is necessary for continuous tumour growth.

Another example of a pathology where latency has proved crucial to understand certainaspects of HIV-1 infection and its response to therapy. HIV infection can be effectivelycontrolled by potent anti-retroviral therapies. However, such therapies do not accomplishtotal eradication of the disease. Although they effectively eliminate virus-producing in-fected T cells, early study of the long-term viral load in treated patients revealed thatlow-level (below the detection limit of standard clinical analysis), residual viremia per-sisted for very long time (months or even years) [110, 123]. These results suggested thatthe infection persisted in the treated patients in the form of a compartment that wasinsensitive to the drugs. Several hypothesis were formulated regarding the nature of thiscompartment, but the one that appears to have consolidated itself as the most likely ex-planation is that of a latently infected population of T cells. This cellular compartmentconsists of slow-proliferating cells (probably memory cells) which store the virus but donot replicate it. Since this cells do not replicate the virus they remain unaffected by theanti-retroviral drugs. Furthermore, upon stimulation with specific antigens, they becomeactive and resume virus production and proliferation, which, in turn, replenishes the latentcompartment. This feed-back cycle sustains the persistent low-level viremia. The presenceof this latent infection has been recognised as major barrier for complete eradication ofHIV-1 infection and the search for ways to tackle it has become a very active of research[121, 113, 65, 66].

Quiescent behaviour is also part of normal physiological regulation of tissues. In par-ticular, stem cells, which are the essential building blocks of tissues in multi-cellular or-ganisms, spend most of their life span in a latent, non-proliferative state. Only uponreception of the appropriate signalling cues, stem cells activate and proliferation ensues.These signalling cues are secreted in response to a decrease in the number of mature(terminally-differentiated) cells that are at the summit of the corresponding differentia-tion cascade: when mature cell numbers decline stem cell proliferation activates in orderto replenish the partially depleted mature cell compartment. Although, this mechanism ispart of normal tissue regulation, it can be subverted by pathologies such as cancer whichhighjack the properties of stem cells.

The overarching aim of this thesis is to use stochastic modelling along with asymp-totic and numerical techniques to better understand the basic mechanisms involved inthe dynamics of cell populations with feed-back regulated latency in order to proposecontrol strategies that help in formulating better, more rationally grounded therapeuticapproaches for pathological situations arising from disregulation of differentiation cascades

Introduction 11

and HIV-1 latent infection.

1.1.2 Stochastic dynamics of differentiation cascades with regulatoryfeed-back

In Chapter 2 we analyse the factors that affect the robustness of cell populations withhierarchical structure. Tissues in higher organisms such as mammals maintain homeosta-sis by means of cascades of serial cell differentiation [4]. These cascades are composedby hierarchically-organised compartments of different cell types. The cornerstone of thissystem is stem cells (SCs) which sustain the tissue by producing both more stem cells,a property known as self-renewal, and lineage-specific differentiated cells. The serial dif-ferentiation cascade proceeds through a number of intermediate stages or compartmentsof transient amplifying cells (TACs) until it terminates with the compartment of fully-differentiated mature cells (MCs).

There has been much discussion around the issue of how SCs accomplish the feat ofundergoing both self-renewal and differentiation [91, 92, 71, 114, 72]. One strategy bywhich SCs can perform these two tasks is asymmetric cell division, whereby, upon SCdivision, one daughter cell retains the cellular identity of its mother whereas its sister dif-ferentiates. Asymmetric cell division seems to proceed through segregation of the materialthat determines cell fate so that, when the SC divides, one of the daughters keep the stemcell-fate determinants. Lineage-specific determinants are thus passed onto its sister cell[71, 72].

Although asymmetric SC division is a simple and elegant solution to the stem cellpuzzle, there has been much controversy around what is perceived to be a fundamentalflaw of this model, namely, asymmetric division does not allow the stem cell compartmentto expand [91, 92]. The argument runs that in situations such as injury, where the rateof cellular turnover must be largely increased in order to regenerate the affected tissue,the size of the SC pool should undergo a marked expansion. This hints that asymmetricdivision might not be the full solution.

Experimental evidence points out to the fact that SCs in some tissues, such as thecentral neural system [50] and the epidermis [74], can divide symmetrically. SymmetricSC division consists of both daughter cells sharing the same fate: SC or differentiated.In this way, the number of cells is regulated by the frequencies of (symmetric) divisionsproducing SCs or differentiated cells. Since this division modality allows for expansion ofthe SC compartment, an alternative model arises [91] where most SCs can divide bothsymmetrically and asymmetrically with the balance between these two proliferation modescontrolled by environmental signals to produce the appropriate number of SCs and differ-entiated cells. However, it has been recognised that, whilst symmetric division confers thecapability of enhanced growth and increased regenerative capacity, it also increases thelikelihood of cancer. This view is supported by findings according to which the cellularmachinery involved in asymmetric division has been evolutionary conserved with a role intumour suppression, the gene products that promote symmetric SC division also functionas oncogenes [91].

Our aim is to analyse the stochastic dynamics of cell differentiation cascades where bothmodes of SC proliferation occur. We will study optimality strategies where by a trade-off between symmetric division (which enhances adaptability) and asymmetric division(which endows stability).

12 Introduction

1.1.3 Stochastic dynamics and extinction of HIV-1 in patients underpotent anti-retroviral therapy

HIV infection has been proved to be effectively controlled by administration of anti-retroviral therapy (ART). Despite the great success of of highly-active ART (HAART),the goal of curing, in the sense of completely eradicating, HIV-1 infection is yet to beachieved [105]. Quantitative analysis of the temporal evolution of plasma viral load uponHAART treatment suggests the existence of several phases in the therapy-induced decay ofthe viral load [102, 69, 111]. After an initial shoulder, reflecting delays associated to boththe pharmaco*kinetics and the production of virus by newly infected cells [101, 55], a firstphase of fast exponential decline of the viral load is observed where viral load is reducedby one to two orders of magnitude over a period of time of approximately two weeks. Thisfast response stage, with half-lime time of the order of days, reflects short half-life time ofplasma virus and of productively infected CD4+ T lymphocytes [57, 124, 83]. Followingthis initial phase of fast decline in plasma viral load, a second stage of slower decay startswith a half-life time between one and four weeks. This phase reflects the contribution tovirus load of infected cells with longer half-life time, such as macrophages, and infectedCD4+ T cells that exhibit a lower rate of viral replication [99, 58, 126]. After this secondstage, plasma virus load has normally fallen below the detection threshold of standardclinical assays (∼ 50 copies RNA/ml). However, following this second phase, HAARTappears to fail to completely eradicate the infection. Rather, a third stage ensues withmuch longer half-life time than the previous ones (of the order of years [69]) in whichresidual levels of viral load (1-5 copies RNA/ml detectable only by supersensitive assays)persist in plasma as well as in other bodily compartments, such as sem*n.

The question of what is the source of this residual viral load has triggered much debatewhich has materialised in several working hypothesis. One of these hypothesis invokesthe possibility that HAART is not completely suppressive thus allowing the infection tocontinue to replicate in anatomical HIV-1 reservoirs [102], in particular within the so-called”drug sanctuaries“, i.e. sites of poor drug penetration where the infection is allowed topersist [67].

An alternative model,which is abundantly supported by recent studies, suggests that,although HAART is fully suppressive, a cellular reservoir exists which allows the infectionto linger in latent form [102] with residual viral load is the result of the activation of thelatently infected cells [102]. Such a latent reservoir is established within the population ofinfected CD4+ T memory cells [25, 24] and, therefore, they remain in the resting state inthe presence of HAART for prolonged periods of time. As a consequence, latently infectedcells are able to escape the effect of the drug and immune surveillance due to the fact thatthey undergo no duplication and, consequently, exhibit very low levels of HIV-1 messengerRNA [111]. However, since latently infected cells release virus when stimulated with theproper antigen, viral rebound will eventually occur when HAART is withdrawn leading toHIV-1 infection recurrence, consistent with a wider scenario of quiescence-induced escape[3].

Once HAART has forced the infection to enter the latent stage, one can observe tran-sient episodes of viremia where the viral load raises above the standard test detectionlimit (50 copies mRNA/ml) for a brief period of time [33, 104, 54, 95]. These episodes arereferred to as blips. The origin and clinical relevance of these blips remains unclear and anumber hypothesis have been formulated [111].

Given the inability of HAART to hit latently infected cells, combined therapies sup-

Introduction 13

pressing the latent reservoir has been proposed [105]. A promising avenue in this directionconsists of combining HAART with agents that specifically activate latently effective cells,thus rendering them susceptible to attack by HAART. The rationale for this combinedtherapy is that increasing the activation rate of the latently infected cells the size of thiscompartment should decay rapidly and, due stochastic fluctuations, these cells eventuallywill be suppressed taking all the infected population with them.

1.2 Mathematical framework

Generally speaking, the problem we address in this thesis concerns the stochastic dynam-ics of populations with different cell types. The framework most suited for dealing withsuch problem is that of population dynamics. Since during the latent state numbers ofindividuals are low, we need to model our populations in terms of stochastic processes, inparticular of Markov processes [44]. The Markov process for the population dynamics willbe formulated in terms of the corresponding Master Equation. These equations are noto-riously difficult to solve exactly and in any practical case one must resort to approximatemethods to deal with it.

In some cases we are dealing with small populations, and in other situations we arelooking for rare events, for that reason we use the semi-classical approximation, which isa WKB approximation [73] to compute the probabilities of rare events.

We include a detailed explanation of the semi-classical approach in Appendix B. Thisapproximation transforms our problem of dealing with a Master Equation to a problemof classical mechanics with a Hamiltonian with several degrees of freedom, as many ascell types. To compute the invariant manifolds of this system, which usually is the firststep in the study of these dynamical system, is not easy and require non-linear approach,techniques such as parametrisation method [17, 18] to obtain high-order expansion. Moredifficulties associated to such systems appear due to their slow-fast nature. Furthermore,the study of motion near hyperbolic points may be necessary. To integrate the equa-tions of motion, Taylor integrator methods [64] combined with, automatic differentiationtechniques [64] allow us to use larger steps achieving higher accuracy than other explicitmethods, such as Runge-Kutta, can provide. The Taylor integrator method, is not a sym-plectic integrator, i.e. it does not preserve the Hamiltonian structure of the equations,however, this is irrelevant, since one can require a local error below the roundof machineerror. In Appendix C we explain the methods used to study these Hamiltonian systems.

Although we have made some progress in large-size asymptotic treatment of theseproblems, we heavily relay in simulation and numerical methods. Our main tool to studythe Master Equation will be the Gillespie Stochastic Simulation Algorithm (SSA), whichis a numerical Monte Carlo technique based on a reinterpretation of the Master Equation.In Appendix A this methodology is explained in detail. The Gillespie SSA requires largecomputational times, in particular when (i) large systems are considered and (ii) whenseparation of time scales appears. In the latter case, time-scale separation can be takenadvantage of formulate a more efficient variant of the SSA proposed by Cao et al. [20],this method is explained in detail in Appendix D.

14 Introduction

1.3 Aims and objectives

The main aim of this PhD Thesis is to study the effect of fluctuations in coupled popu-lations of biological systems where cells in latent states play a major role. In doing so,we try to put to forward some biologically meaningful aspects of the dynamics of thesesystems. The specific points we aim to cover and the organisation of the thesis are asfollows.

In Chapter 2, we study the behaviour of cell populations with hierarchical structurefrom the point of view of their stability properties. In particular,

1. Symmetric vs Asymmetric division in the stem cell compartment. We study therobustness of hierarchically-organized populations, depending on whether stem cellsdivide symmetrically, asymmetrically, or both. We assess how the symmetric divisionaffects the stability of the population, regarding the role of stem cells symmetricdivision of SC in cancer progression.

2. The competition between population with different modes of stem cell division. Thisis crucial in order to find optimal architectures to maximize the robustness (long-timesurvival, resistance to invasion and ability to invade) of hierarchical populations.

3. The influence of parameters such as the duplication and death rate of the stem cells,the average life-time of the fully differentiated cells, the length of the differentia-tion cascade and the fluctuations of population of stem cells on the robustness andoptimality of differentiation cascades.

In Chapter 3, an hom*ogeneous model of combined HAART and latently infected cellactivation therapies of the HIV-1 dynamics in the blood stream is presented. We areinterested in

1. The effect of the activation rate of the latently infected cells in the extinction timeof the infection. In particular we analyse if the therapies based on increase this rateare able to suppress the infection in a reasonable time.

2. The importance of the efficiency of the anti-retroviral therapies, including the limit-ing case of an efficiency of 100%, in the size of the persistence viral load.

3. The formulation of an asymptotic theory based on the semi-classical quasi-steadystate approximation to describe the dynamics of the process. The accuracy of thisasymptotic method is compared to the multi-scale SSA proposed by Cao et al.

In Chapter 4, we extend a model proposed by Rong and Perelson [110] to a non-hom*ogeneous model of the HIV-1 dynamics in the blood stream, we consider that cellsand virons are not uniform distributed in the blood stream. The specific points we aim toaddress are:

1. The mechanism for the emergence of transient episodes of viremia, known as viralblips. In particular, we investigate if they are the product of inhom*ogeneous densityfluctuations or another mechanism must be considered.

2. If the appearance of viral blips is affected by laboratory procedures, such as the timeelapsed between the blood sample extractions and the observations.

Introduction 15

3. Whether the likelihood, amplitude and frequency of the viral blips is affected by thedifferent of possible types of viral production, namely, continue vs burst.

In Chapter 5, we proceed to discuss our results, compare, when possible, with previousmodels or experimental results, and we give directions for future work.

Details regarding methodological issues, as well as an introduction to the stochas-tic modelling via the Master Equation, are give in the appendices. For those who arenot familiarised with models based on Master Equations, the author recommends readfirst Appendix A which provides the mathematical background to understand Chapter 2.Appendices B, C and D together with Appendix A provide the necessary mathematicalbackground to follow Chapter 3 and Chapter 4.

16 Introduction

Chapter 2

Robustness of differentiationcascades with symmetric stem celldivision

Stem cells perform the task of maintaining tissue homeostasis by both self-renewal anddifferentiation. Whilst it has been argued that stem cells divide asymmetrically, there isalso evidence that stem cells undergo symmetric division. Symmetric stem cell divisionhas been speculated to be key for expanding cell numbers in development and regenera-tion after injury. However, it might lead to uncontrolled growth and malignancies suchas cancer. In order to explore the role of symmetric stem cell division, we propose amathematical model of the effect of symmetric stem cell division on the robustness of apopulation regulated by a serial differentiation cascade and we show this may lead to ex-tinction of such population. We examine how the extinction likelihood depends on definingcharacteristics of the population such as the number of intermediate cell compartments.We show that longer differentiation cascades are more prone to extinction than systemswith less intermediate compartments. Furthermore, we have analysed the possibility ofmixed symmetric and asymmetric cell division against invasions by mutant invaders inorder to find optimal architecture. Our results show that more robust populations arethose with unfrequent symmetric behaviour.

2.1 Biological Background

There is large body of literature examining the properties of differentiation cascades, inparticular in the haematopoietic system [1, 29, 82] and in relation to several kinds ofleukaemias and haematological diseases [26, 27, 87, 77]. Marciniak et al. [82] have studieda deterministic model in which they explore the role of regulation in a system with asym-metric SC division. They have found that for their model to exhibit efficient repopulation,regulation by environmental signals must occur at the level of the fraction of self-renewal.Mackey and co-workers [26, 27, 77] have done extensive modelling and analysis of a numberof haematological diseases in which circulating cellular blood components exhibit oscilla-tory behaviour. By means of delay-differential-equation models and exhaustive parametersensitivity analysis, they have proposed which processes in the haematopoietic differen-tiation cascade are more likely to be involved in the development of disease. They have

17

18 Robustness of differentation cascades

also found that delays in the differentiation cascade are fundamentally linked to the on-set of the oscillations, a mechanism akin to the extinction mechanism we explore here.Dingli et al [32] have studied a model of the architecture and dynamics of hematopoiesisto provide estimates of the number of compartments, the size of each compartment, andthe corresponding replication rates.

There has also been a number of stochastic models of serial differentiation cascadesthat tackle evolutionary issues. The model proposed by Pepper et al. [98] deals with anumber questions arising when considering which evolutionary purpose can such a costlysystem (from the energetic point of view) possibly serve. The conclusion of their study isthat differentiation cascades provides with intrinsic protection against somatic mutationsleading to malignancies such as cancer. Evolution towards malignant behaviour and cancerhas been studied, for example, by Sun & Komarova [119] and Rodriguez-Brenes et al.[106, 107], where the mechanisms of regulation of normal tissues are studied in detail toascertain how can they be subverted in order to produce cancer. Other related work isthe symmetric stem cell division models proposed by Dingli et al. [31], were a Moranprocess is proposed to study the competition between a population of normal stem cellsand a population of cancer stem cells. They observed that a larger duplication rate in theinvader (cancer) stem cells increases the probability of invasion.

Michor and co-workers [39, 75, 87, 88, 86] have modelled the dynamics of myleoidleukemia using models consisting of hierarchical differentiation cascades. It has beenfound that cancer initiation requires only a single mutation. Their models suggest thatimatinib is a potent inhibitor of the production of differentiated leukemic cells, but doesnot deplete leukemic stem cells. Lenaerts et al [78] have shown that although imatinib,does not eradicate chronic myleoid leukemic, due to the stochastic nature of hematopoiesis,leukemic stem cells undergo extinction. This result shows the importance of stochasticeffects on the dynamics of hierarchical cell populations.

Zhao et al. [127] have proposed a linear model for colonic epithelial cells, that explicitlytakes into account the proliferation kinetics of a cell as a function cell position within thecrypt, They have observe that those cells with mitotic activities concentrated near thestem cells delay the rate of mutation accumulation in colonic stem cells.

Mutations in hierarchically structured populations have been previously studied byWerner et al. [125], they modelled a general hierarchically organized multi-compartmenttissue, allowing any number of mutations in a cell. They showed that these tissues stronglysuppress the cells that are carrying multiple mutations. Traulsen et al [120] have studiedthe effect of different kinds of mutations and showed that although neutral mutations oftenleads to clonal extinction, disease is still possible, and in this case, it may become difficultto eliminate neutral mutations with therapy.

Also related to the context of this section, a mechanism for stochastic evolutionaryescape and its consequences for drug resistance based on the presence of a quiescent,drug-resistant population has been recently described [3]. Since quiescent cancer stemcells have been recently found in intestinal crypts [16], this escape mechanism can berelevant to ascertaining the efficacy of therapies targeting cancer stem cells.

The aim of this chapter is to address some of the issues arising from the considerationof symmetric SC division in serial differentiation cascades, in particular those regardingthe evolution of uncontrolled growth and the risk of cancer. To do so, we consider astochastic mathematical model of the serial differentiation cascade and use it to addresstwo problems (i) the long-term viability (stability) of a population of cells maintained by a

Robustness of differentation cascades 19

differentiation cascade whose SCs undergo symmetric division only, and (ii) the likelihoodof invasion of a mutant population with increased symmetric SC division in competitionwith a population with asymmetric SC division (robustness). Regarding the former, weshow that serial differentiation cascades with symmetric SC differentiation exhibit anintrinsic instability mechanism which induce extinctions with high likelihood. Thus, thesecascades seem to have an in-built fail-safe mechanism against excessive symmetric SCdivision. The competition model allows us to estimate optimality conditions in, in thesense of maximum robustness, for the rate of symmetric division in terms of maximisingits chance to take over another population with less symmetric-division capability whilstminimising its chances of undergoing extinction.

The organisation of the section is as follows. Section 2.2 is devoted to presenting ourmodel of a differentiation cascade with symmetric SC division in detail and discussing itsunderlying hypothesis. In Section 2.3, we present simulation results for this model anddiscuss the emergence of a mechanism of extinction induced by symmetric SC division. InSection 2.4, we analyse the robustness of the different division strategies by formulating amodel for the competition between populations using different such strategies and extractconclusions regarding their optimality. Finally, in Section 2.5, we summarise our results.

2.2 Stochastic modelling of a serial differentiation cascadewith symmetric SC division

In this section we proceed to present and discuss our stochastic model of a hierarchically-organised cell population. Our aim is to formulate a regulated stochastic model which, byincorporating symmetric SC division, can lead to changes in the size of SC compartment.Regulation is achieved by means of a cytokine which is produced by the MC compartmentwhich controls the rate of apoptosis of the stem cells.

The stochastic approach is justified by the small number of stem cells, typically presentin their host organism. For example, according to [91], in the haematopoietic system onlya 0.003 % of the cells are SCs. Given these small numbers in the SC compartment, weexpect random effects to be important if not dominant.

2.2.1 Model description

We consider a compartmental model of the serial differentiation cascade with symmetricSC division. As shown in Fig 2.1 each compartment corresponds to a stage in the cascade:from the SC compartment to the MC compartment with a number of intermediate TACcompartments. We further introduce in our model the secretion of a cytokine by the MCswhich regulates the size of the SC compartment. We study the dynamics of the populationof cells in each compartment. In particular, we examine the behaviour of a differentiationcascade of length n+1, with stem cells, n−1 partially-differentiated, transient amplifyingcells, and mature cells: x0 represents the number of SCs, xi with i = 1, ..., n − 1, thenumber of TACs at differentiation stage i, and xn, the number of MCs. Furthermore,xn+1 represents the amount of a cytokine which regulates the size of the SC compartmentthrough a negative-feedback [81, 26, 27, 118]. Such cytokine is secreted at a rate wichis proportional to the number of fully-mature cells. Therefore, this negative feedbackaccounts for the mature-cell regulation of the size of the SC compartment which has beenincorporated, for example, in models of haematopoeitisis [82].

20 Robustness of differentation cascades

Figure 2.1: Diagrammatic representation of the hierarchic differentiation cascade. x0

represents the population of the SC compartment, xi the population of TAC cells at istage, xn the population of the MC compartment and xn+1 represents the amount ofcytokine. Plot (a) corresponds to the well-stirred scenario, whereas plot (b) correspondsto the delayed scenario (see main text for details).

Robustness of differentation cascades 21

Transition rate rj = (∆xj0,∆xj1, ...,∆x

jn+2) Description

W1 = px0 (1, 0, ..., 0) SC self-renewal

W2 = d0x0 (−1, 2, 0, ..., 0) SC differentiation

W3 = λ0x0xn+1 (−1, 0, ..., 0,−1, 0) SC apoptosis

W2+2i = λixi (0, ...,−1, ..., 0) TAC apoptosis

W3+2i = dixi (0, ...,−1, 2, 0, ..., 0) TAC differentiation

W2+2n = λnxn (0, 0, ...,−1, 0, 0) MC apoptosis

W3+2n = sxn (0, 0, ..., 0, 1, 0) Cytokine secretion

W4+2n = λn+1xn+1 (0, 0, ..., 0,−1, 0) Cytokine clearance (well-stirred scenario)

W5+2n = s2xn+1 (0, 0, ..., 0,−1,+1) Cytokine transport

W6+2n = λn+2xn+2 (0, 0, ..., 0,−1) Cytokine clearance (delayed scenario)

W7+2n = λ0x0xn+2 (−1, 0, ..., 0,−1) SC apoptosis

Table 2.1: Transition rates corresponding to the stochastic model of the serial differentia-tion cascade with symmetric stem cell division. A description of the corresponding elemen-tary population-dynamical processes is given in Section 2.2.1. We choose d0 = 1−p day−1,to have that stem cells divide, in average, once per day, λi = 0.01 day−1 and di = 1− λiday−1 for i = 1, ..., n−1. s is fixed in the following way: from equation (2.5), at the steady

state we have 0 = sxn − λn+1xn+1 − λ0x0xn+1, this implies s = λn+1xn+1+λ0x0xn+1

xn.

Our stochastic model is formulated in terms of the corresponding master equation[122, 44]:

∂P (X, t)

∂t=∑j

(Wj(X − rj , t)P (X − rj , t))−Wj(X, t)P (X, t)) (2.1)

where X = (x0, . . . , xn, xn+1) and P (X, t) is the probability of the population vector tohave value X at time t. For full specification of our stochastic models, we need to prescribethe transition rates, Wj(X, t), corresponding to the probability per unit time of the processj. These rates are modelled in terms of the law of mass action [47]. rj is the change in Xwhen elementary process j occurs, i.e. X(t+ ∆t) = X(t) + rj with probability Wj∆t.

We will consider two different scenarios schematically represented in Fig. 2.1. Thefirst one, corresponding to Fig. 2.1.(a) and which will be referred to as the well-stirredscenario, the cytokine which controls the size of stem cell compartment is assumed to beinstantaneously transported to the it. By contrast, the second scenario, corresponding toFig. 2.1.(b) and to be referred to as the delayed scenario, the cytokine requires a finiteamount of time to reach the stem cell compartment. We have modelled this by introducinga compartment where, upon secretion, the cytokine is transported to at a certain rate.The elementary processes involved in our model are:

• SCs divide symmetrically, so that they can undergo:

1. Self-renewal. Both daughter cells share the SC fate of the mother with transi-tion rate W1 (see Table 2.1)

2. Differentiation. Both daughter cells differentiate to first-stage TAC with tran-sition rate W2 (see Table 2.1)

22 Robustness of differentation cascades

3. Apoptosis with cytokine-regulated transition rate W3 if we are considering thewell-stirred scenario or W7+2n, otherwise (see Table 2.1)

• Stage-i, i = 1, . . . , n− 1, TACs can undergo:

1. Differentiation. Both daughter cells differentiate to either the next TAC stage(if i < n − 1) or to the MC compartment (if i = n − 1) with transition rateW2+2i as given in Table 2.1

2. Apoptosis with transition rate W3+2i as given in Table 2.1

• MCs can undergo:

1. Release cytokine with transition rate W3+2n as defined in Table 2.1

2. Apoptosis with transition rate W2+2n as defined in Table 2.1

• Cytokines can undergo:

1. Clearance with transition rate W4+2n if we are considering the well-stirred sce-nario, or W6+2n, otherwise (see Table 2.1)

2. Transport to the intermediate compartment (xn+2 in Fig. 2.1.(b)) with transi-tion rate W5+2n if we are considering the delayed scenario (see Table 2.1).

We will analyse this system by means of numerical Monte-Carlo simulations usingthe Gillespie stochastic simulation algorithm [47, 48]. The mean-field limit of the modelalso provides useful information. For the model described by Eq. (B.1) and Table 2.1,corresponding to the well-stirred scenario , the mean-field model is given by the followingsystem of ODEs:

x0 = px0 − d0x0 − λ0x0xn+1, (2.2)

xi = 2di−1xi−1 − λixi − dixi, i = 1, ..., n− 1, (2.3)

xn = 2dn−1xn−1 − λnxn, (2.4)

xn+1 = sxn − λn+1xn+1 − λ0x0xn+1, (2.5)

which has two fixed points. The trivial (unstable) equilibrium (0, 0, ..., 0), which corre-sponds to the extinction of the system, a positive equilibrium fixed point, PS

PS =

(λ1x1 + d1x1

2d0, ...,

λj+1xj+1 + dj+1xj+1

2dj, ...,

λnxn2dn−1

, xn,p− d0

λ0

)(2.6)

We will fix the number of xn cells at the metastable state. We denote as xi, the numberof cells in the i compartment at the metastable state.

Our model assumes that the signal is produced by the mature cells, we do this as-sumption instead of consider that the signal is produced by, for example, the transientamplifying cells in position k. In this case the signal producing rate s, should be of theform:

s =λn+1xn+1 + λ0x0xn+1

xk. (2.7)

Robustness of differentation cascades 23

Figure 2.2: Gillespie simulation of a system with some parameter values that kill in ashort time the population. The black line corresponds to the x0 population and red andto xn and xn+1 respectively. p = 0.6 day−1, λn = 0.03 day−1, λn+1 = 1 day−1, λ0 = 1

4000day−1, the Length of the differentiation cascade is 4.

The signal producing rate necessary to the cells in each compartment to keep the amountof signal necessary to control the population. We see that the only efficient possibility isthat the signal is generated by the mature cells. However, this possibility of TAC cellsproducing signal is not out of our study. If the signal is produced by the xk cells, the newsystem is equivalent to a system described above with length k taking as a new λ′k (thedeath rate of the new mature cells) the value of λk + dk, since the xk+1, ...xn cells has noeffect on the system.

In the scenario (b), that is, the cytokine is not transported immediately to the stem cells

compartment, we consider another compartment, xn+2, and the reactions xn+1s2−→ xn+2

and xn+2λn+2−→ ∅. changing x0 + xn+1 → ∅ by x0 + xn+2 → ∅. In the Section 2.3 we will

see the effect of this change in function of the transport rate s2.

2.3 Extinction induced by symmetric SC division

We now proceed to analyse the stochastic model described in Section 2.2.1 and describe anextinction mechanism intrinsic to the structure of a serial differentiation cascade with sym-metric cell division. We interpret this mechanism as an anti-cancer mechanism inherentto the hierarchical structure of serial differentiation cascades.

24 Robustness of differentation cascades

2.3.1 Delay-induced extinction in serial differentiation cascades withsymmetric SC division

In order to proceed with our analysis, we perform numerical simulations of our model bymeans of the Gillespie stochastic simulation algorithm [47, 48]. This algorithm is basedon an exact reformulation of the stochastic process described by the master equation Eq.(B.1) which allows us to generate exact sample paths or realisations of the underlyingprocess. An example of such a realisation which illustrates the extinction mechanism inour model with symmetric SC division is illustrated in Fig. 2.2.

Fig. 2.2 shows simulation results where the population has been initially set to beequal to the (stable) positive equilibrium of the mean-field limit Eqs (2.2)-(2.5). This par-ticular example is illustrative of the mechanism through which extinction occurs. Initiallya fluctuation in the system induces an increase in the number of SCs. This increase inSC population propagates through the differentiation cascade and it eventually reachesthe MC compartment which starts expanding. Such expansion in MC numbers inducesan increase in cytokine concentration, which, in turn, leads to decay in the size of SCcompartment. As a consequence, oscillatory behaviour ensues. However, since the sys-tem, specially the SC compartment, is subject to fluctuations and these are amplifiedby the length of the differentiation cascade (an increase of one cell in the size of the SCcompartment induces an increase O(2n) in the MC compartment), the amplitude of theoscillations grow until the SC compartment size hits the absorbing barrier at which pointthe population undergoes extinction.

Since delays appear to be instrumental for this extinction pathway, we are interestedin the effects of altering the balance between the xn+1/xn cells life-time and the replenishrate via x0 on the average extinction time. To carry on this analysis, we compute theextinction time as a function of the SC proliferation rate, p, and the death rate of theMCs, λn. Our results are shown in Figure 2.3.

Fig. 2.3 shows that there is a region in the λn-p-plane where the average extinctiontime, TE , is such that TE > 104 days. For all other values of these two parameters, thewaiting time for extinction to occur becomes much shorter. It is relatively straightforwardto see that the effect of these two parameters within the region of the parameter spacewhere TE > 104 is to cancel the effect of the delay between the SC and MC compartmentsdue to the length of the differentiation cascade. Large values of the MCs death rate, λn,stabilises the system via a dissipative-like effect, namely, it impairs the ability of the MCcompartment to expand in response to fluctuations in the SC compartment size. It alsoweakens the effect of the feed-back between the MC and SC compartments, as it effectivelyreduces the concentration of cytokine.

The effect of varying the value of the self-renewal rate, p, is slightly different. Whenpd0

& 1, the SC compartment cannot grow, on the contrary, under the effect of SC apoptosisand random fluctuations, it will dwindle and the population will become extinct. As pgrows, its effect is strongly coupled to λn. If λn is small, i.e. the MC life expectancy isincreased, the effects of fluctuations and delay are reinforced and the population undergoesrapid extinction. If, on the contrary, λn is large, we observe two regimes: A regime ofintermediate values of p where the effects of fluctuations and delay are cancelled by therapid death of the MCs (see above), and a regime of high values of p where MC death is notfast enough to compensate for the increased effects of fluctuations in the SC compartmentwhich take over and produce short survival times.

Another model parameter that has an obvious effect on the delay is the length of the

Robustness of differentation cascades 25

Figure 2.3: This plot represents the extinction time of the system depending on theparameters λn and p. The length of the differentiation cascade is n = 5, the x-axis is thedeath rate of the mature cells, the y-axis is the duplication probability of the stem cells,the z-axis is the days needed to get an extinction, the white color means that more than10000 days (∼ 30 years) are needed to get an extinction. λn+1 = 1 day−1, The number ofmature cells at the metastable state is xn = 5000, λ0 = 1

2500 day−1.

Figure 2.4: Sample paths of the dynamics of our serial differentiation cascade model ob-tained by means of the stochastic simulation algorithm. Plot (a) corresponds two differentcascades with length n = 4 and length n = 6 with the same number of MCs. Plot (b),idem fixing the number of SCs, instead. We observe that, in both cases, shorter cascades(n = 4) are more stable that longer ones (n = 6). λn+1 = 1 day−1, λ0 = 1

6000 day−1,λn = 0.25 day−1, p = 0.6.

26 Robustness of differentation cascades

Figure 2.5: Simulation results corresponding to the dependence on the length of the dif-ferentiation cascade, n, of the probability of extinction within a time window (0, TE),PE(TE), with TE = 1000 days. p = 0.53 day−1, λn = 0.25, λn+1 = 1, xn = 26000 ,λ0 = 1

26000 day−1.

differentiation cascade, n. Longer chains promote the effect of the delay and thereforefavours extinctions. This is illustrated in Fig. 2.4, where we show sample paths of ourmodel for different values of n. We see that indeed longer differentiation cascades de-stabilise the system and favour extinction. This is further explored in Fig. 2.5, wherewe show the extinction probability within the time window of time of duration (0, TE),PE(TE), as a function of n. Fig. 2.5 shows a sharp increase in this quantity as n increases.

2.3.2 Delayed cytokine scenario

In this section we consider the delayed cytokine scenario where:

• The secreted cytokine, xn+1, is transported to xn+2, that is, xn+1s2−→ xn+2, with

transition rate transition rate W5+2n = s2xn+1

• The cytokine in xn+2 is degraded, that is, xn+2λn+2−→ ∅, with transition rate transition

rate W6+2n = λn+2xn+1

• The size of the SC compartment is controlled by xn+2, that is x0 + xn+2 → ∅, withtransition rate transition rate W7+2n = λ0x0xn+2

We start by analysing the mean-field behaviour of the system:

x0 = px0 − d0x0 − λ0x0xn+1, (2.8)

xi = 2di−1xi−1 − λixi − dixi, i = 1, ..., n− 1, (2.9)

xn = 2dn−1xn−1 − λnxn, (2.10)

xn+1 = sxn − λn+1xn+1 − s2xn+1, (2.11)

xn+2 = s2xn+1 − λn+2xn+1 − λ0x0xn+2, (2.12)

Robustness of differentation cascades 27

Figure 2.6: Extinction time in function of the transport rate s2. p = 0.6 day−1, λn = 0.03,λn+1 = 1, λn+2 = 1, xn = 10000 , λ0 = 1

10000 day−1, n = 5.

which implies that the secretion rate, s, is now determined by

s =(s2 + λn+1)(λn+2xn+2 + λ0x0xn+2)

xns2, (2.13)

and the positive equilibrium is of the form:

PS =

(λ1x1 + d1x1

2d0, ...,

λj+1xj+1 + dj+1xj+1

2dj, ...,

λnxn2dn−1

, xn,λn+2xn+2 + λ0x0xn+2

s2,p− d0

λ0

).

(2.14)Note that if λn+2 = λn+1, we recover the previous model when s2 →∞.

We are interested in how the dependence of the behaviour of population on the cytokinedelay, which is controlled by s2. Fig 2.6 shows the extinction time as a function of s2. InFig 2.7, we have plotted s as a function of s2, as given by Eq 2.13. Taken together, thesetwo results show that a decrease in s2 contributes to destabilise the population, since thiscontributes to increase the delay between cytokine secretion and SC response. Moreover,more cytokine needs to be produced to maintain the SC compartment size.

2.4 Robustness of populations with asymmetric SC division

So far, we have focused on the study of differentiation cascades with symmetric stem celldivision. We have uncovered that symmetric stem cell division may induce extinction ofthe population, specially in long differentiation cascades where the delays in the regulationof stem cells by mature cells are more important. This result points out that utilisationof symmetric stem cell alone may not be an optimal strategy.

In normal situations, stem cell division is commonly believed to be mostly asymmetricwith symmetric division being a rare event [71]. However, it has also been put forwardthat symmetric SC division could also play an important role, specially in situations, suchas wound repair, where the SC compartment must be allowed to expand in order to com-pensate for an increased demand in cellular turnover [74]. In view of our results regardingthe induction of instabilities by symmetric SC division, one could ask the question of howfrequent symmetric SC division can be before anomalies in the population start to appear.

28 Robustness of differentation cascades

Figure 2.7: Value of s according to Eq. 2.13 in function of s2. p = 0.6 day−1, λn = 0.03,λn+1 = 1, λn+2 = 1, xn = 10000 , λ0 = 1

10000 day−1.

The aim of this section is to address these two issues, namely, whether symmetric SCdivision is an optimal strategy and the frequency of symmetric SC division, as well as as-certaining which are the generic features of an optim, i.e. maximally robust, differentiationcascade.

2.4.1 Stochastic modelling of a differentiation cascade with asymmetricstem cell division

Our first step is to present our stochastic model for a differentiation cascade with asym-metric stem cell division. We formulate our model in terms of the transition rates that gointo the corresponding master equation Eq. (B.1), which are shown in Table 2.2. Modelformulation is very similar to the one for the differentiation cascade with symmetric SCdivision: we consider the same n + 2 different cell and chemical species, i.e. SC, n − 1TACs, MCs, and the MC-secreted cytokine.

The dynamics of this system is essentially the same as in the symmetric SC divisioncase (see Table 2.1) except for the SC cell dynamics. Following Knoblich [71], we assumethat SCs divide asymmetrically with occasional symmetric cell division is a rare event. Wethus consider that SCs divide asymmetrically with probability 1− e. We further considerthan the rate of asymmetric division is regulated (inhibited) by a cytokine produced bythe mature cells [82]. These two factors are summarised in the transition rate W1, Table2.2. With probability e, SCs behave as a symmetrically-dividing cell, whose behaviour ischaracterised by transition rates W1, W2, and W3, Table 2.2, i.e. a SC which behave as asymmetrically-dividing cell can (i) self-renew, (ii) differentiate, or (iii) undergo apoptosis,with probability rates W1, W2, and W3, respectively. In our study, the parameters a andb will be chosen to have the same division rate as in the symmetric division model at themetastable state, that is, (1−e)a

b+yn+1+ ep+ ed0 = p+ d0.

The mean-field behaviour of the model is described by the following system of law-of-mass action ODEs:

Robustness of differentation cascades 29

Reacting rate rj = (∆yj0,∆yj1, ...,∆y

jn+1) Description

W1 = (1−e)ab+yn+1

y0 (0, 1, ..., 0) asymmetric division–differentiation

W2 = epy0 (1, 0, ..., 0) symmetric self-renewal

W3 = ed0y0 (−1, 2, 0, ..., 0) symmetric differentiation

W4 = eλ0y0yn+1 (−1, 0, ..., 0,−1) SC apoptosis

W3+2i = λiyi (0, ...,−1, ..., 0) TAC apoptosis

W4+2i = diyi (0, ...,−1, 2, 0, ..., 0) TAC differentiation

W3+2n = λnyn (0, 0, ...,−1, 0) MC apoptosis

W4+2n = syn (0, 0, ..., 0, 1) Cytokine secretion

W5+2n = λn+1yn+1 (0, 0, ..., 0,−1) Cytokine clearance

Table 2.2: Reaction rates corresponding to the stochastic dynamics of a differentiationcascade with both symmetric and asymmetric division. Symmetric division is consideredto be a rare event which occurs with probability e [72]. We choose d0 = 1 − p day−1,λi = 0.01 day−1 and di = 1− λi day−1 for i = 1, ..., n− 1. s is fixed in the following way:from equation (2.19), at the steady state we have 0 = syn − λn+1yn+1 − eλ0y0yn+1, this

implies s = λn+1yn+1+eλ0y0yn+1

yn.

Figure 2.8: Relative size of the stationary average size SC compartment at the metastablestate in function of e, y0(e)/y0(e = 0), as determined by Eq. 2.20. We have normalisedwith y0(e = 0). Parameter values: a = 1000, b = 700, λ0 = 1

10000 day−1, yn = 10000,p = 0.55 day−1, λn = 0.05 day−1, λn+1 = 1 day−1, the Length of the differentiationcascade is 4.

30 Robustness of differentation cascades

y0 = epy0 − ed0y0 − eλ0y0yn+1 (2.15)

y1 =(1− e)ab+ yn+1

y0 + 2ed0y0 − λ1y1 − d1y1 (2.16)

yi = 2di−1yi−1 − λiyi − diyi i = 2, ..., n− 1 (2.17)

yn = 2dn−1yn−1 − ynλn (2.18)

yn+1 = syn − λn+1yn+1 − eλ0y0yn+1 (2.19)

Where yi is the number of cells in the compartment i. The variables and parametersare described in the Table 2.2. Eqs. (2.15)-(2.19) have a positive equilibrium given by:

PA =

λ1y1 + d1y1

(1−e)ab+yn+1

+ 2ed0

, ...,λj+1yj+1 + dj+1yj+1

2dj, ...,

λnyn2dn−1

, yn,p− d0

λ0

(2.20)

Fig. 2.8 shows the average stationary size of the stem cell compartment as predicted byEqs. (2.15)-(2.19). We observe that, whilst an increase in the size of the SC compartmentis obtained as e is let to grow, the magnitude of such increase, for the chosen set ofparameter values, of the order of a 10 %.

Regarding the behaviour of the stochastic model, the main difference with respect tothe symmetric cell division model (Table 2.1) involves the lack of delay-induced extinctionbehaviour exhibited by differentiation cascade model with pure symmetric SC division.This behaviour is illustrated in Fig. 2.9 where two typical sample paths are shown fordifferent values of the probability of symmetric division, e. If the probability of symmetricdivision, e, is small, the variation in size of the SC compartment is bounded, which rendersthe system resilient to the effect of the delay-induced oscillations that, in the pure symmet-ric SC division system, eventually lead to extinction of SC compartment followed by theextinction of the whole population. Fig. 2.10, we have plotted the extinction probabilityof a differentiation with asymmetric SC division as a function of e and show that, indeed,differentiation cascades with small probability of symmetric division are impervious to thedelay-induced extinction mechanism.

2.4.2 Competition between populations

We now move on to study of the robustness of the asymmetric population, that is, whetherit can withstand attempts of invasion by other populations using different strategies. Inwhat follows, such strategies will refer to variations of two parameters: The probability ofsymmetric division, e and the rate of symmetric self-renewal p. Our aim is to ascertain cri-teria of evolutionary optimality based on our results regarding the outcome of competitionbetween populations characterised by different values of these parameters.

Asymmetric vs symmetric population

We start by studying the competition between a resident population generated by differen-tiation cascade with asymmetric SC division, hereafter referred to as A, and a symmetric-SC-division invader, referred to as S, whose population dynamics is given by the modeldescribed in Section 2.2.

Robustness of differentation cascades 31

Figure 2.9: Sample paths showing the time evolution of the number of SCs for differentvalues of e. Panel (a) corresponds to e = 0.05 whereas panel (b) shows a realisation withe = 1. Parameter values: a = 1000, b = a − xn+1, λ0 = 1

5000 day−1, p = 0.55 day−1,λn = 0.03 day−1, λn+1 = 1 day−1 n = 4, xn = 5000 . We can see that the population withsymmetric stem cells division becomes extinct much faster than the asymmetric.

Figure 2.10: Simulation results corresponding to the dependence on the non asymmetricdivision probability of stem cells e, of the probability of extinction within a time window(0, τ), PE(τ), with τ = 1000 days. p = 0.55 day−1, λn = 0.03, λn+1 = 1, xn = 6000,λ0 = 1

6000 day−1 a = 1000, b = a− xn+1.

32 Robustness of differentation cascades

Mean-field behaviour. In order to proceed further, we start by giving some resultsregarding the behaviour of the mean-field limit using qualitative theory of ODEs. Forsimplicity, we will assume that both populations of mature cells are producing the samecytokine.

Let y0 (x0) be the size of SC compartment of the asymmetric (symmetric) divisionpopulation. The mean-field behaviour of the model of competition between the populationsA and S is described by the following system of law-of-mass action ODEs:

y0 = epy0 − ed0y0 − eλ0y0xn+1 (2.21)

y1 =(1− e)ab+ yn+1

y0 + 2ed0y0 − λ1y1 − d1y1 (2.22)

yi = 2di−1yi−1 − λiyi − diyi i = 2, ..., n− 1 (2.23)

yn = 2dn−1yn−1 − λnyn (2.24)

x0 = e′p′x0 − e′d′0x0 − e′λ′0x0xn+1 (2.25)

x1 =(1− e′)a′

b′ + xn+1x0 + 2ed′0x0 − λ′1x1 − d′1x1 (2.26)

xi = 2d′i−1xi−1 − λ′ixi − d′ixi i = 2, ..., n− 1 (2.27)

xn = 2d′n−1xn−1 − λ′nxn (2.28)

xn+1 = syn + s′xn − λn+1xn+1 − eλ0xn+1y0 − e′λ′0xn+1x0 (2.29)

Let the quantities r and r′ be defined as:

r =p− d0

λ0

r′ =p′ − d′0λ′0

(2.30)

i.e. r and r′ are the ratios between the net SC growth rate and the SC differentiationrate for each of the competing populations, which coincide with steady-state value of themean-field cytokine concentration for each of the populations.

From Eq. (2.21) we see that y0 > 0, i.e. y0 grows in time, if xn+1 ∈ (0, r). On thecontrary y0 decreases as time progresses if xn+1 ∈ (r,∞). Similarly, Eq. (2.25) impliesthat x0 is increasing with time when xn+1 ∈ (0, r′) and decreasing when xn+1 ∈ (r′,∞).Moreover, according to Eqs. (2.21) and (2.25), when t → ∞, the stationary level ofcytokine is given by xn+1 → max(r, r′). We thus consider three possible scenarios for themean-field limit:

• r > r′. The S population becomes extinct: as t→∞ xn+1 → r which implies that,as t→∞, y0 = 0 and x0 ≤ 0 with the equality in the latter holding for x0 = 0 only.

• r = r′. Both populations coexist: as t → ∞, xn+1 → r = r′ which implies that, ast→∞, both y0 = 0 and x0 = 0 hold.

• r < r′. By the same reasoning as in the first case, the A population becomes extinct.

Regarding the behaviour of the stochastic system, there are some important differenceswith respect to its mean-field counterpart. In Figure 2.11, we show that the probability of

Robustness of differentation cascades 33

invasion as a function of the self-renewal rate, p′, and death rate λ′0, of the invader stemcells. Instead of the 3 scenarios that characterise the mean-field limit, we have:

• r ≥ r′. In this case, extinction of S is the most likely scenario: the population S isforced to evolve under higher on-average cytokine levels which induces upregulationof apoptosis of the S-stem cells. As r′ → r, this competitive advantage of A overS weakens. However, at that point, the delay-induced extinction mechanism for Stakes over. Since we are considering e to be small, delays have little effect on thepopulation A, which leads to eventual extinction of S. This behaviour is illustratedin Fig. 2.11, which shows that the probability of S to invade A vanishes in thisregime.

Figure 2.11: Probability invasion of the competition of two populations, one with asym-metric division with e = 0.1, and the other with pure symmetric division, e′ = 1, and weplot the probability of invasion in function of the duplication rate of the invader populationand for 2 different values of invader stem cells death. a = 1000, b = a − yn+1, λ0 = 1

6000day−1, yn = 6000, p = 0.55 day−1, λn = 0.03 day−1, λn+1 = 1 day−1.When λ′0 = 1

6000 ,r′ = r for p′ = 0.55, λ′n = 0.03 day −1, λ′n+1 = 1 day −1, a′ = 1000, b′ = a − xn+1 whenλ′0 = 1

3000 , r′ = r for p′ = 0.6. x0(0) = 1, xi(0) = 0 ∀i 6= 0.

• r < r′. Under these conditions, the population A is the one which evolves underexcess of cytokine conditions, as the average steady-state concentration of cytokineis higher than the one the A population would have in the absence of the populationS. As a consequence, the probability of S invading A is now positive, as shown inFig. 2.11. We also observe that the probability of invasion as a function of r′ reachesa maximum value and then starts decreasing. This behaviour is again due to delay-induced extinctions: as the relative net growth rate of the S-stem cells increases theprobability of delay-induced extinction of the S population also increases (as shownin Fig. 2.3).

34 Robustness of differentation cascades

In order to further explore the behaviour of the stochastic competition model in ther < r′ regime, we have done simulations to compute the probability of invasion and theprobability of extinction of both populations as a function of e, i.e. the probability ofasymmetric behaviour of A-stem cells, and the rate of (symmetric) self-renewal of S-stemcells. The results are shown in Fig. 2.12. We observer (Fig. 2.12(a)) that there is athreshold value of e below which the resident A-type population cannot be invaded bya S-type invader, regardless of the value of p′. If the A-type resident has an e abovethis threshold, then two outcomes are possible: invasion, for moderately high values ofp′, or extinction of the symmetric invader (see also Fig. 2.11), for even higher values ofp′. Furthermore, as e grows, Fig. 2.12(b) shows that the more likely outcome is bothpopulations be extinct. This is due to the fact that for larger values of e the symmetriccomponent becomes dominant in population A as well as S, and consequently the A-population starts to becoming affected by the effects of the corresponding delays. Theprobability of long term coexistence is negligible, according to our simulations.

Optimal differentiation cascades

Finally, we tackle the issue of examining the optimality properties of differentiation cas-cades. In this context, we define optimality of a differentiation cascade in terms of itsresilience to invasion, namely, optimal differentiation cascades are those which are maxi-mally resilient to invasion.

In order to address this issue, we again consider the competition between two popula-tions generated by differentiation cascades characterised by different values of the proba-bility of symmetric behaviour, e, and the self-renewal rate, p. We study the ability of theresident population to withstand invasion as a function of the values of these parametersfor both populations. We have chosen these two particular quantities as the focus of ouranalysis, as the are directly related to the ability for adaptive behaviour of differentiationcascades, which has been argued to be at the root of symmetric cell division [91].

We start by analysing the optimal strategy regarding the probability of symmetricbehaviour, e. To this end, we have ran simulations of the competition between a residentwith probability of symmetric behaviour, e, and an invader with probability of symmetricbehaviour, e′. All the other parameters are taken to be equal for both populations. Resultsare shown in Fig. 2.13.

In the light of our results regarding delay-induced extinction for higher values of e, i.e.when the symmetric part of the dynamics dominates, we expect that there exists a trade-off between resilience to invasion and long-term survival and elevated frequency of SCsymmetric behaviour (which favours adaptive responses). This expectation is confirmedby the simulations results shown in Figs. 2.13(a), which show the probability of invasion.We can observe that resident populations with e 1 are very stable and resilient toinvasion. As e increases, so that a faster and more efficient adaptive response be possible,an invader with e′ ≤ e is likely to wipe out the resident population (see Fig. 2.13(a)).

We have also investigated the effect of the rate of self-renewal on the resilience againstinvasion of an A-type population. Our results are reported in Fig. 2.14 where we have plot-ted the probability of invasion (Fig. 2.14(a)), the probability of extinction (Fig. 2.14(b)),and the probability of coexistence after a time TF has elapsed (2.14(c)). These simulationshave been done for e = 0.1 and e′ = 0.05, which correspond to the region of parametervalues where, according to Fig. 2.13, the probability of invasion is positive. The resultsshown in Fig. 2.14 indicate that if p > p′, i.e. if the rate self-renewal of the resident is

Robustness of differentation cascades 35

Figure 2.12: Simulation results for the competition between a resident with mixed asym-metric division, i.e. 0 < e < 1, and an invader with pure symmetric division, i.e. e′ = 1.Plot (a) and plot (b) show, respectively, the probability of invasion and the probabilityof extinction of both populations as a function of the duplication rate of the symmetricdivision and the probability of symmetric-like division of the asymmetric population, e.a = 1000, b = a − yn+1, λ0 = 1

3000 day−1, xn = 6000, p = 0.55 day−1, λn = 0.03 day−1,λn+1 = 1 day−1. We run simulations until a final time= 6000 days. e′ = 1, λ′0 = 1

3000day−1, xn = 6000, p′ = 0.55 day−1, λ′n = 0.03 day−1, λ′n+1 = 1 day−1. x0(0) = 1,xi(0) = 0 ∀i 6= 0.

36 Robustness of differentation cascades

Figure 2.13: Simulation results for the competition between two A-type populations. Plot(a) and (b) show the probability of invasion and extinction of both populations, respec-tively. We have plotted these three quantities as a function of the probability of symmetricbehaviour of the resident, e, and the invader e′. a = 1000, b = a− yn+1, λ0 = 1

3000 day−1,yn = 6000, p = 0.55 day−1, λn = 0.05 day−1, λn+1 = 1 day−1, a′ = 1000, b′ = a′ − xn+1,λ′0 = 1

3000 day−1, xn = 6000, p′ = 0.55 day−1, λ′n = 0.05 day−1, λ′n+1 = 1 day−1, We runsimulations until a final time= 6000 days. x0(0) = 1, xi(0) = 0 ∀i 6= 0.

Robustness of differentation cascades 37

Figure 2.14: Simulation results for the competition between two A-type populations. Plot(a), (b), and (c) show the probability of invasion, extinction of both populations, and prob-ability of coexistence, respectively. We have plotted these three quantities as a functionof the probability of the stem cells duplication resident, p, and the invader p′. a = 1000,b = a − yn+1, λ0 = 1

3000 day−1, yn = 6000, e = 0.1 day−1, λn = 0.05 day−1, λn+1 = 1day−1, a′ = 1000, b′ = a′ − xn+1, λ′0 = 1

3000 day−1, xn = 6000, e′ = 0.05 day−1, λ′n = 0.05day−1, λ′n+1 = 1 day−1, We run simulations until a final time TF = 6000 days. x0(0) = 1,xi(0) = 0 ∀i 6= 0.

38 Robustness of differentation cascades

larger than that of the invader, the resident population is resilient to invasion.Taken together, the results shown in Figs. 2.13 and 2.14 show that the optimal trade-

off between resilience to invasion and capability for adaptive behaviour is realised bypopulations with small probability of symmetric behaviour, e, and large value of the self-renewal rate.

2.5 Conclusions

In this chapter, we have presented an analysed stochastic models of differentation cas-cades to study the effect of symmetric and asymmetric stem cell division. Our model isa compartmental model in which each compartment corresponds to a different differenti-ation stage, we have considered stem cells, transient amplifying cells, fully differentiatedcells and a cytokine, produce by the FD cells which regulates the size of the stem cellcompartment. In particular, we have deduced from our model some protection mecha-nism of the hierarchically structured cell populations. We have proposed a new extinctionmechanism in which, the delay between the changes due noise fluctuations in the stemcell compartment, and the cytokine compartments, unstablises the whole population andeventually leads the system to extinction. Moreover, we have investigated a more realisticsituation in which both, symmetric and asymmetric stem cell division are present, and wehave looked for optimal strategies to be robust against invasion.

Chapter 3

Antigen-stimulation inducederadication of the HIV-1 infectionin patients under highly activeanti-retroviral therapy

HIV-1 patients under potent anti-retroviral develop a small reservoir of latent cells whichpersists and that anti-retroviral therapies are unable to eliminate, which hinders the totaleradication of the infection. Recent research activity has been focused on formulatingstrategies aimed at removing the latent infection [121, 113, 65, 71]. One such proposalconsists of elevating the rate of activation of the latently infected cells. We present astochastic model of dynamics of the HIV-1 infection and study the effect of the rate oflatently infected cell activation on the expected extinction time of the infection. We analyseour model by means of an asymptotic approximation using the semi-classical quasi steadystate approximation (QSS). We also present a lower bound of the extinction reducingour stochastic problem to an easier 1-dimensional system which, in some cases, is a goodapproximation of the original one. We test the accuracy of our asymptotic results bymeans of a hybrid multi-scale stochastic simulation algorithm.

3.1 Introduction

HIV-1 infected patients are effectively treated with highly active anti-retroviral therapy(HAART). Whilst HAART is successful in keeping the disease at bay with average levelsof viral load well below the detection threshold of standard clinical assays, it fails tocompletely eradicate the infection [63]: the infection persist in form of a latent a reservoirwith a half-life time of years. This implies that life-long administration of HAART is, atthe moment, necessary for HIV-1-infected patients, which is prone to drug resistance andcumulative side effects as well as imposing a considerable financial burden on developingcountries, those more afflicted by HIV, and public health systems [121].

The latent reservoir consists of a population of latently infected T cells (probably,memory cells) which host, but do not produce the virus. Since virion production is notoccurring within latently infected cells, they are immune to effects of HAART [110]. Duringtheir life-cycle, latently infected cells can undergo activation thus becoming active infected

39

40 Antigen-stimulation

cells with the ability to produce virus. The latently infected cell population is replenishedby both active infected cell proliferation and also by slow, density-dependent homeostaticproliferation of the memory CD4+ T memory cells (and, therefore, of the latently infectedcell compartment), which, according to [23], drives persistence and determines the size ofthe latent reservoir. This cycle is able to maintain positive, albeit small levels of viral loadfor very long time.

The presence of this latent infection has been recognised as major barrier for completeeradication of HIV-1 infection and the search for combination therapies, i.e. HAART plusspecific agents that tackle the latent reservoir, has become a very active field of research[121, 113, 65, 66]. One of such potential approaches consist using agents that activate la-tently infected cells: by activating the cells in the latent reservoir, they would be renderedsensitive to HAART. This would, at least theoretically, clear the latent reservoir and,eventually, the infection. Early studies regarding reactivation of the latent reservoir usedinterleukins (IL-2, IL-3, and IL-7) in combination with specific antibodies (CD3). Theseattempts failed because it increased the absolute count of T cells [121, 113]. More recentstudies have focused on the use of small molecules that reactivate latent virus produc-tion without inducing global T cell activation, in particular several histone deacetylases(HDAC) and other chromatin modifiers [121]. Here, we focus on modelling the efficiencyof the latter type of therapy in combination of HAART. In particular, our aim is to analysehow the average extinction time of the infected cell population (both active and latent)changes as the activation rate of latently infected cells is increased.

We perform a WKB asymptotic analysis of the partial differential equation for the prob-ability generating function associated to the Master Equation which describes stochasticpopulation dynamics. This method relays on the solution of a variational problem for theoptimisation of an action functional, which effectively reduces the problem to the analy-sis of the associated Hamilton equations [73, 35]. In order to make further progress, wetake advantage of the multiplicity of time scales in the system, to perform a quasi-steadystate approximation (QSSA) on the Hamilton equations, thus reducing the dimensional-ity of the problem [2]. We further used advanced numerical techniques to analyse themulti-dimensional Hamiltonian system: computation of invariant manifolds [53] and Tay-lor integration method [64] together with automatic differentiation techniques [64], whichare used to obtain high accuracy and larger time steps. Our asymptotic results are testedby means of multi-scale stochastic simulations [20].

The organisation of the chapter is as follows. Section 3.2 is devoted to presenting ourmodel of the HIV-1 dynamics under HAART in detail and discussing its underlying hy-pothesis. In Section 3.3 we present an asymptotic analysis of our stochastic process, whichtransforms it into a classical mechanics problem via the semi-classical approximation. Weperform a time-scale analysis of the resulting Hamilton-Jacobi equations to use a Quasisteady state approximation to reduce the dimension of the problem. Moreover, we presenta simplification of the model which allows us to reduce the system to an easy handle 1-dimensional system which gives us a lower bound of the extinction time. In Section 3.4we present the results obtained via the different methods presented and a comparison ofthem. Finally, in Section 3.5 we summarise our results.

Antigen-stimulation 41

Transition rate rj = (∆T,∆L,∆T ∗,∆V ) Description

W1 = λΩ (1, 0, 0, 0) Recruitment of T cells

W2 = dTT (−1, 0, 0, 0) Death of T

W3 = η(1− ε)kV TΩ−1 (−1, 1, 0,−1) Infection T + V =⇒ L cell

W4 = (1− η)(1− ε)kV TΩ−1 (−1, 0, 1,−1) Infection T + V =⇒ T ∗ cell

W5 = rL (0, 1, 0, 0) Proliferation L =⇒ 2L

W6 = rLmax

L∗(L−1)2 Ω−1 (0,−2, 0, 0) L+ L =⇒ ∅

W7 = d0L (0,−1, 0, 0) Death of LiW8 = aLL (0,−1, 1, 0) Activation L→ T ∗

W9 = δT ∗ (0, 0,−1, 0) Death of T ∗

W10 = cV (0, 0, 0,−1) Clearance of virus

W11 = pvT∗ (0, 0, 0, 1) Production of virus by T ∗

W12 = εkV TΩ−1 (0, 0, 0,−1) Infection failed

Table 3.1: Transition rate corresponding to the stochastic model of HIV dynamics. Adescription of the corresponding elementary population-dynamical processes is given inSection 3.2.

3.2 Stochastic Model

Our model is a stochastic generalisation of the model proposed by Rong and Perelson[110]. It takes into account 4 different variables: T represents the number of CD4+ T-cells that are susceptible to HIV-1 infection, L, which stands for the number of latentlyinfected cells that cannot produce virus. L cells can be activated by stimulation with theirrecall antigens. Finally, T ∗ and V represent productively infected cells, i.e. those that canproduce virus particles, and the total viral load, respectively.

Our stochastic model is formulated in terms of the corresponding master equation[122, 44]:

∂P (X, t)

∂t=∑j

(Wj(X − rj , t)P (X − rj , t))−Wj(X, t)P (X, t)) (3.1)

where X = (x1, x2, x3, x4) = (T, L, T ∗, V ) and P (X, t) is the probability of the populationvector to have X particles at time t. For full specification of our stochastic model, weneed to prescribe the transition rates, Wj(X, t), corresponding to the probability per unittime of the elementary processes j to occur. These rates are modelled in terms of thelaw of mass action [47]. rj is the change in X when elementary process j occurs, i.e.X(t+ ∆t) = X(t) + rj with probability Wj∆t: P (X(t+ ∆t) = X(t) + rj |X(t)) = Wj∆t.

The elementary processes involved in our model are:

• Recruitment of new CD4+T cells with transition rate W1 (see Table 3.1).

• CD4+T cells can undergo:

42 Antigen-stimulation

Parameter Description Value

T0 CD4+T cells 599999 cells ml−1

at the metastable state

L0 Latently infected cells 0.934778 cells ml−1

at the metastable state

T ∗0 Productively infected cells 0.115067 cells ml−1

at the metastable state

V0 Viral load 10 copies ml−1

at the metastable state

T1 CD4+T cells 600000 cells ml−1

at the extinction

L1 Latently infected cells 0 cells ml−1

at the extinction

T ∗1 Productively infected cells 0 cells ml−1

at the extinction

V1 Viral load 0 copies ml−1

at the extinction

λ Recruitment rate of T cells 10000 ml−1 day−1

dT Death rate of T cells 0.0166 day−1

k infection rate 2.4 · 10−8 ml day−1

ε Drug efficacy 0.85

η Fraction resulting in latency 0.001

dL Death rate of latently infected cells 0.001 day−1

asL Standard rate of transition from latently to productively 0.1 day−1

aL rate of transition from latently to productively varied

σ death rate of productively infected cells 1 day−1

c clearance rate of free virus in blood stream 23 day−1

pv Viral production rate 2000 day−1

r proliferation rate of activated cells 0.2 day−1

Lmax Carrying capacity density of latent cells 1.888 cells ml−1

Ω System size 5000 ml

Table 3.2: Parameter values of the different rates

Antigen-stimulation 43

1. Apoptosis with transition rate W2 (see Table 3.1).

• Latently infected cells Li can undergo:

1. Homeostatically balanced proliferation. Following [110], we account for ho-moeostatic control of proliferation by means of a combination of branching(Li → 2Li) with binary annihilation (Li+Li → ∅). It has been shown (see e.g.[35]) that this combination is an stochastic counterpart of the standard logisticgrowth. The corresponding transition rates are W3 for branching and W4 forbinary annihilation (see Table 3.1).

2. Death. We assume a simple linear decay with transition rate W5 as shown inTable 3.1.

3. Activation. By means of this process a latently infected cell becomes an activeinfected cell Li → T ∗i . The corresponding transition rate is W6, see Table 3.1.

• Active infected cells, T ∗i , cell are subjected to:

1. Apoptosis with transition rate W7 (see Table 3.1).

2. Virion production. Contrary to latently infected cells, active infected cellssynthesise and release new virions. In general, viral production can occur ina continuous fashion over the life span of an infected cell or in a burst whichkills the cell. For the HIV infection both modes have been proposed [97].Transition rate W9 (Table 3.1) corresponds to continuous production. Later inthe chapter, we consider an additional scenario, in which both continuous andburst production occur.

• Finally, virus, Vi, can:

1. Infect a healthy T cell producing a latently infected cell. In patients underHAART, the infection process is hindered by the presence of an anti-retroviraldrug. The efficiency of such drug is measured by a parameter, ε, which takesvalues between 0 and 1, the latter (former) corresponding to a maximally(in)efficient drug. (1 − ε) is interpreted as the proportion of virions capableof infection under HAART treatment. We also assume that, upon infection,the cell can become latently infected with probability η or active with proba-bility (1 − η). Therefore the corresponding transition rate W8 is proportionalto η(1− ε) as shown in Table 3.1.

2. Infect a healthy T cell producing an active infected cell. In this case, thecorresponding transition rate W2 is proportional to (1 − η)(1 − ε) (see Table3.1).

3. Undergo clearance. Virions are removed from the blood, and we model thisprocess as a simple linear decay with transition rate W8 as per Table 3.1.

4. Fail to infect and being eliminated by the drug with transition rate W10 (seeTable 3.1).

The parameter values used in our analysis are based on current estimates available inthe literature. They are summarised in Table 3.2.

44 Antigen-stimulation

As a first approach to the dynamics, one can study the mean-field behaviour of thissystem, which is given by the set of ordinary differential equations (3.2), (3.3), (3.4) and(3.5).

d

dtT (t) = λ− dTT − (1− ε)kV T (3.2)

d

dtL(t) = η(1− ε)kV T + rL(1− L

Lmax)− d0L− aLL (3.3)

d

dtT ∗(t) = (1− η)(1− ε)kV T − δT ∗ + aLL (3.4)

d

dtV (t) = NδT ∗ + pvT

∗ − cV − kV T (3.5)

The mean-field system has two fixed points. The point with coordinates (T1, 0, 0, 0),i.e. the disease-free equilibrium associated to extinction of the infection, is a repeller fixedpoint. The other fixed point, of the form (T0, L0, T

∗0 , V0) with L0, T

∗0 , V0 6= 0, is a global

attractor. When noise is considered, the latter equilibrium becomes a metastable state,whereas the former becomes an absorbing state [35]. From Eqs. (3.2)-(3.5), we obtainthat this metastable state exists as long as aL < a∗L, where:

a∗L =r − d0

1−(

η(1−ε)kT1δ(c+kT1)pv−1−(1−η)(1−ε)kT1)

) (3.6)

3.3 Asymptotic analysis and numerical methods

In this section we study the stochastic dynamics described by Eq. (3.1). We follow themethodology presented in [2] and we use the semi-classical approximation, but performinga time-scale analysis on the resulting Hamilton equations to simplify to reduce the dimen-sion of the resulting problem. A multi-scale stochastic simulation method, which takesadvantage of the separation of time scales to produce fast simulations, is used to test ourresults.

3.3.1 Semi-classical approximation

In Appendix B we have included all details of the semi-classical approximation, in thissection we apply that methodology to our specific problem.

Consider the generating function associated to the probability distribution Eq. (3.1).

G(p1, p2, p3, p4, t) =∞∑n=0

px11 px22 p

x33 p

x44 P (x1, x2, x3, x4, t). (3.7)

The evolution of the generating function is determined by a partial differential equation

Antigen-stimulation 45

(PDE), which is derived from the Master Equation (3.1):

∂G

∂t= −λ(p1 − 1)− dT (1− p1)

∂G

∂p1− (1− η)(1− ε)k(p3 − p1p4)

∂G

∂p1

∂G

∂p4

−η(1− ε)k(p2 − p1p4)∂G

∂p1

∂G

∂p4− r(p2

2 − p2)∂G

∂p2− r

2Lmax(1− p2

2)∂2G

∂p22

−dL(1− p2)∂G

∂p2− aL(p3 − p2)

∂G

∂p2− pv(p3p4 − p3)

∂G

∂p3− δ(1− p3)

∂G

∂p3

−c(1− p4)∂G

∂p4− εk(p1 − p1p4)

∂G

∂p1

∂G

∂p4. (3.8)

Which can be interpreted as a Schrodinger equation,

∂

∂tG = −HG, (3.9)

where the Hamiltonian operator, H, in the p (momentum) representation is

H(p1, p2, p3, p4, q1, q2, q3, q4) = λ(p1 − 1) + dT (1− p1)q1 + (1− η)(1− ε)k(p3 − p1p4)q1q4

+η(1− ε)k(p2 − p1p4)q1q4 + r(p22 − p2)q2 +

r

2Lmax(1− p2

2)q22

+dL(1− p2)q2 + aL(p3 − p2)q2 + pv(p3p4 − p3)q3 + δ(1− p3)q3

+c(1− p4)q4 + εk(p1 − p1p4)q1q4 (3.10)

where qi ≡ − ∂∂pi

. The operators pi and qi satisfy the commutation relation [pi, qj ] = δi,j .From Eq. (3.9), we formulate an Ansatz where G = exp (−S). When the system

size, Ω, is big, one may employ the WKB approximation and consider a S with S of theform S = −ΩS expanding up to order O(Ω−1) (neglecting O(Ω−2) terms), one obtains theclassical Hamilton-Jacobi equations

∂S

∂t= H

(p,∂S

∂p

), (3.11)

Instead of directly tackling with the Hamilton-Jacobi equation, which, in general, wedo not know how to solve, we exploit the analogy of the Shrodinger equation and use aFeynman path-integral representation we obtain a solution of (3.8), that is

G(p, t) =

∫ t

0(exp (−S(p, q))Dq(s)Dp(s)) ds+ px0 , (3.12)

with

S(p, q, t) =

∫ t

(−H(p, q) +

n∑i=1

qi(s)pi(s)

)ds. (3.13)

The semi-classical approximation consists of approximating the path integral (B.13)by

G(p, t) ≈ exp (−S(p, t)), (3.14)

where pi are the solution of the Hamilton equations of the (classical) Hamiltonian

H(p, q) = λ(p1 − 1) + dT (1− p1)q1 + (1− η)(1− ε)k(p3 − p1p4)q1q4

+η(1− ε)k(p2 − p1p4)q1q4 + r(p22 − p2)q2 +

r

2Lmax(1− p2

2)q22

+dL(1− p2)q2 + aL(p3 − p2)q2 + pv(p3p4 − p3)q3 + δ(1− p3)q3

+c(1− p4)q4 + εk(p1 − p1p4)q1q4. (3.15)

46 Antigen-stimulation

This implies, in particular (B.6) corresponds to the classical Lagrangian action related to(3.15). The equations of motion are given by dqi/dt = ∂H/∂pi, dpi/dt = −∂H/∂qi [6],that is

dp1

dt= dT (p1 − 1) + (1− η)(1− ε)k(p1p4 − p3)q4 + η(1− ε)k(p1p4 − p2)q4 + εk(p1p4 − p1)q4,

dp2

dt= r(p2 − p2

2) +r

Lmax(p2

2 − 1)q2 + dL(p2 − 1) + aL(p2 − p3),

dp3

dt= pv(p3 − p3p4) + δ(p3 − 1),

dp4

dt= (1− η)(1− ε)k(p1p4 − p3)q1 + η(1− ε)k(p1p4 − p2)q1 + c(p4 − 1) + εk(p1p4 − p1)q1,

dq1

dt= λ− dT q1 − (1− ε)kp4q1q4 + εk(1− p4)q1q4,

dq2

dt= η(1− ε)kq1q4 + r(2p2 − 1)q2 −

r

Lmaxp2q

22 − dLq2 − aLq2,

dq3

dt= (1− η)(1− ε)kq1q4 + aLq2 + pv(p4 − 1)q3 − δq3,

dq4

dt= −kp1q1q4 + pvp3q3 − cq4. (3.16)

These equations are formally solved with boundary conditions

qi(0) = xi(0),

pi(t) = pi. (3.17)

3.3.2 Quasi-steady state approximation

In this section, we proceed to formulate a quasi-steady approximation, that allows us tosimplify the system Eq. (3.16). We proceed by re-scaling our variables in such a way tomake explicit the separation of time scales and simplify the model according to the quasi-steady state approximation. For clarity, we perform this analysis in two steps. We startby defining the following set of re-scaled (dimensionless) quantities:

In this section, we proceed to re-scale variables so that we can analyse the intrinsictime scales of each of the variables involved and simplify the model according to the quasi-steady state approximation. For clarity, we perform this analysis in two steps. We startby defining the following set of re-scaled (dimensionless) quantities:

q1 =q1

T0,

q2 =q2

L0,

q3 =q3

T ∗0,

q4 =q4

V0,

s = kV0t. (3.18)

Antigen-stimulation 47

Parameter Order of magnitude

κ1 = λkT0V0

O(105)

κ2 = dTkV0

O(105)

κ3 = rL0kV0T0

O(100)

κ4 = dLL0kV0T0

O(10−2)

κ5 = aLL0kV0T0

O(100)

κ6 =pvT ∗0kV0T0

O(103)

κ7 = δT0kV0T0

O(100)

κ8 = ckT0

O(103)

Table 3.3: Re-scaled parameters following the re-scaling of variables shown in Eq. 3.18

In this re-scaled variables, the Hamiltonian can be written as

H(p, x) = kT0V0Hκ(p, x), (3.19)

where Hκ(p, q) is given by

Hκ(p, q) = κ1(p1 − 1) + κ2(1− p1)q1 + (1− η)(1− ε)(p3 − p1p4)q1q4

+η(1− ε)(p2 − p1p4)q1q4 + κ3(p22 − p2)q2 + κ3

L0

2Lmax(1− p2

2)q22

+κ4(1− p2)q2 + κ5(p3 − p2)q2 + κ6(p3p4 − p3)q3 + κ7(1− p3)q3

+κ8(1− p4)q4 + ε(p1 − p1p4)q1q4, (3.20)

and where, for simplicity of the notation, we have dropped the hats of the re-scaledvariables qi. The parameters κi are given in Table 3.3.

The Hamiltonian equations for the re-scaled variables read as follows

kV0T0dq1

ds= kV0T0

∂Hκ

∂p1, kV0

dp1

ds= −kV0T0

T0

∂Hκ

∂q1,

kV0L0dq2

ds= kV0T0

∂Hκ

∂p2, V0

dp2

ds= −kV0T0

L0

∂Hκ

∂q2,

kV0T∗0

dq3

ds= kV0T0

∂Hκ

∂p3, kV0

dp3

ds= −kV0T0

T ∗0

∂Hκ

∂q3,

kV 20

dq4

ds= kV0T0

∂Hκ

∂p4, kV0

dp4

ds= −kV0T0

V0

∂Hκ

∂q4,

(3.21)

and, therefore,dq1

ds=

∂Hκ

∂p1,

dp1

ds= −∂Hκ

∂q1,

L0

T0

dq2

ds=

∂Hκ

∂p2,

L0

T0

dp2

ds= −∂Hκ

∂q2,

T ∗0T0

dq3

ds=

∂Hκ

∂p3,

T ∗0T0

dp3

ds= −∂Hκ

∂q3,

V0

T0

dq4

ds=

∂Hκ

∂p4,

V0

T0

dp4

ds= −∂Hκ

∂q4.

(3.22)

48 Antigen-stimulation

In order to proceed further with our time-scale analysis, we note that the dynamicalequations for q1 and q4 are dominated by the terms in κ1 and κ2, which, according toTable 3.3, are O(105), and κ6 and κ8, which are O(103), respectively. On the contrary,the dominant terms in the equations for q2 and q3 are O(1). In view of this, we re-scalethe dimensionless time s as

T = κ1s. (3.23)

In terms of the re-scaled time T , our Hamilton equations are given by

dq1

dT=

1

κ1

∂Hκ

∂p1,

dp1

dT=

1

κ1

∂Hκ

∂q1,

ε1dq2

dT=

∂Hκ

∂p2, ε1

dp2

dT=

∂Hκ

∂q2,

ε2dq3

dT=

∂Hκ

∂p3, ε2

dp3

dT=

∂Hκ

∂q3,

ε3dq4

dT=

1

κ6

∂Hκ

∂p4, ε3

dp4

dT=

1

κ6

∂Hκ

∂q4,

(3.24)

where ε1 = κ1L0T0

= O(10−1), ε2 = κ1T ∗0T0

= O(10−2), and ε3 = κ1κ6

L0T0

= O(10−3), and all theright hand sides of the Hamilton equations are now O(1). In view of this separation oftime scales, we can assume that q1 is a slow variable which will remain roughly constantin relation to q2, and q3 (and similarly q4) is a fast variable which can be considered tobe at equilibrium with the rest of the system. Therefore, the quasi-steady state Hamilton(QSSH) equations read

dq1

dT=

∂Hκ

∂p1,

dp1

dT= −∂Hκ

∂q1,

ε1dq2

dT=

∂Hκ

∂p2, ε1

dp2

dT= −∂Hκ

∂q2,

ε2dq3

dT=

∂Hκ

∂p3, ε2

dp3

dT= −∂Hκ

∂q3,

ε3dq4

dT=

1

κ6

∂Hκ

∂p4, ε3

dp4

dT= − 1

κ6

∂Hκ

∂q4.

(3.25)

In following sections we consider two different approximations. In both approximationswe assume that the healthy cells are constant, that is p1 = 1 and q1 = 1, this hypothesisholds since if we project on the plane (p1, q1), the stochastic extinction and the metastablestate are very close, the distance is O(10−5), and that the equations for p1 and q1 areessentially linear in p1 and q1 plus a small perturbation in the other variables. The firstapproximation consists in considering ε2 = 0 and ε3 = 0, we refer this as Strong QSS. Theother option consists in considering only ε3 = 0, this approximation is refereed as WeakQSS. In this particular problem, the Strong QSS and the Weak QSS lead us to a verysimilar solutions, see Figure 3.4. However, in other cases their results could be different,for this reason we give tools to study the Weak QSS and to compare them.

Antigen-stimulation 49

Strong QSS

The so-called strong QSSA consists of applying QSS conditions to the equations forq3, p3, q4 and p4:

dq1dT = 0, dp1

dT = 0dq2dT = ∂H

∂p2, dp2

dT = − ∂H∂q2

0 = ∂H∂q3, 0 = ∂H

∂q4,

0 = ∂H∂p3

, 0 = ∂H∂p4

.

(3.26)

That is,

dp2

dT=

1

ε1

(−κ3(p2

2 − p2)− κ3L0

Lmax(1− p2

2)q2 − κ4(1− p2)− κ5(p3 − p2)

),

dq2

dT=

1

ε1

(η(1− ε)q1q4 + κ3q2(2p2

2 − 1)− κ3L0

Lmaxp2q

22 − κ4q2 − κ5q2

),

(3.27)

where

p1 = 1,

q1 = 1,

p3 =κ6 − κ7 − Cκ6 −Bκ6p2 ±

√(κ6 − κ7 − Cκ6 −Bκ6p2)2 − 4κ7κ6A

2κ6A,

q3 =−κ5

(1− η)(1− ε)q1Dp3 + κ6(Ap3 +Bp2 + C)− κ7q3q2,

p4 =−(1− η)(1− ε)q1p3 − η(1− ε)p2q1 − κ8 − εp1q1

−(1− η)(1− ε)q1p1 − η(1− ε)p1q1 − κ8 − εp1q1

q4 =κ6

p1q1 + κ8p3q3,

with

A =−(1− η)(1− ε)q1

−(1− η)(1− ε)p1q1 − η(1− ε)p1q1 − κ8 − εp1q1

B =−η(1− ε)p2q1

−(1− η)(1− ε)p1q1 − η(1− ε)p1q1 − κ8 − εp1q1

C =−κ8 − εp1q1

−(1− η)(1− ε)p1q1 − η(1− ε)p1q1 − κ8 − εp1q1

D =κ6

p1q1 + κ8.

We can discard the positive branch of p3 for continuity. A phase portrait of the StrongQSSA is shown in Figure 3.1.

50 Antigen-stimulation

Figure 3.1: Phase diagram of the equations of motions with the Strong QSS assumptionin the supercritical case (left) and the subcritical case (right).

Weak QSS

The Weak QSS equations in which quasi-steady state conditions are applied only to theequations for p4, q4 reads

dp2

dT=

1

ε1

(−κ3(p2

2 − p2)− κ3L0

Lmax(1− p2

2)q2 − κ4(1− p2)− κ5(p3 − p2)

),

dq2

dT=

1

ε1

(η(1− ε)q1q4 + κ3q2(2p2

2 − 1)− κ3L0

Lmaxp2q

22 − κ4q2 − κ5q2

),

dp3

dT=

1

ε2(−κ6(p3p4 − p3)− κ7(1− p3)) ,

dq3

dT=

1

ε2((1− η)(1− ε)q1q4 + κ5q2 + κ6(p4 − 1)q3 − κ7q3) , (3.28)

where

p1 = 1,

q1 = 1,

p4 =−(1− η)(1− ε)q1p3 − η(1− ε)p2q1 − κ8 − εp1q1

−(1− η)(1− ε)q1p1 − η(1− ε)p1q1 − κ8 − εp1q1

q4 =κ6

p1q1 + κ8p3q3.

3.3.3 Heteroclinic orbits and extinction time

As we have seen in Section 3.2, this system has two different regimes: aL ≥ a∗L, in whichonly the trivial fixed point exists, and aL < a∗L, where in addition to the trivial fixedpoint, a positive fixed point (associated to a persistent infection exists. If this metastable

Antigen-stimulation 51

state exists, the mean field behaviour, which is exactly the dynamics along pi = 1, hastwo different fixed points: an attractor (which corresponds to the metastable state) and arepeller (which corresponds to the extinction).

As time goes to infinity, the trajectories approach the 0 level energy. The action isS = 0 on pi = 1 and on x2 = x3 = x4 = 0. However, it is positive on the heteroclinicconnection between the metastable state and the stochastic extinction, which is of theform x2 = x3 = x4 = 0, 0 < pi < 1.

Since the trivial fixed point corresponds to the only absorbing state of the stochasticdynamics, and the population can not explode, extinction takes place with probability 1.If the metastable state is of the form xi 0, that is the system size Ω 0, the extinctionprobability P0(t) = G(0, t) can be approximated by:

G(0, t) ≈ 1− exp (− tτ

), (3.29)

We look for τ of the formτ = AΩBexp(ΩC), (3.30)

where C is the integral of the re-scaled action, S, along the heteroclinic connection betweenthe metastable state and the stochastic extinction. The semi-classic approximation isvalid when Ω 0. We will perform a numerical fit of A and B using the exact solution(computed with Gillespie simulations) of the system for large values of Ω. A full derivationof this approximation is included in Appendix B.

Numerical fit of A and B

It is clear that as Ω → ∞ the solution of the semi-classical approximation tends to theexact solution of Eq. (3.1), and we can use this fact in our advantage to determine A andB. We also have to notice that we can determine analytically C, as it is the integral ofthe action along the heteroclinic connection. A direct fit of a function of the form

f(Ω) = AΩB exp (CΩ)

can be very unstable. We rather analyse

g(x) =1

Ωlogf(Ω) = Ax+ Bx log (x) + C,

with x =1

Ω, which is more favourable, since first we are interested in C, therefore we want

to have C alone with small terms suppressed by Ω. Our procedure is as follows,

1 Given a aL we select an initial Ω and we compute τ using stochastic simulations (themulti-scale stochastic simulation employed to perform this simulations is explainedlater).

2 Save εΩ =1

Ωand δΩ = εlog

(1

τ

). Increase Ω and go to 1 until we have some values

of (εΩ, δΩ).

3 using an implementation of the non-linear least-squares (NLLS) Levenberg–Marquardtalgorithm [96] fit g(x) = Ax + Bxlog(x) + CS via A,B,CS where x = log(δΩ) andthe right-hand term is εδΩ with some of the higher values of Ω we have.

52 Antigen-stimulation

4 If |CS − C| > tol, where C is the one obtained analytically, increase Ω and go to 1.

5 If |CS − C| < tol we choose CS = C and we fit again via A and B.

6 As we can expect B of the form B = ±n2 with n ∈ N, we choose B as the closest

number of this form to B.

7 finally we fit again just via A.

8 As we are interested in Ω = 5000ml, the mean extinction time, obtained with thesemi-classical approximation is, τs = −

(g(

15000

)5000

).

Note that we expect B to take this form, since this system exhibits the so-called Stokesphenomenon [6]. Therefore it is natural to expect that there is a singularity in the complexspace which is a pole. Moreover, for the binary-Annihilation-decay process, it has beenproven that B = 1

2 [8, 9, 10]. Numerically B approaches 12 (independently of aL). To

obtain analytic expressions for A and B requires a careful analysis of the related Stokesphenomenon [10]. This is postponed for future work.

3.3.4 Lower Bound

We now present another approximation which provides a lower bound to the extinctiontime. The variable x2, associated to latently infected cells, is the most important variablesince the drug is unable to clear it. It is also the slowest variable and it is affected bythe other variables only via the reaction W3, which is the slowest reaction in the system.Numerical simulations suggest the ratio W3/W5 is negligible (see Fig. 3.2). If we considerW3 = 0, only the reactions W5, W6, W7 and W8 affect x2. In this case, we can reduce our4-dimensional stochastic system to a much simpler one dimensional system. The extinctiontime, τLB, of the x2-cells in this reduced system is a lower bound to the extinction timeof the x2-cells in the 4-dimensional system, and provided η is small or ε is large enoughit is a good approximation. Therefore, we can use it to get a first approximation of theaverage extinction time.

Note that, to study the extinction probability of the latently infected cells, this lowerbound is equivalent to consider either a perfect efficacy of the drug (ε = 1) or η = 0, sincein these cases W3 = 0.

Under this assumption, our system is described by the following master equation

∂P (L, t)

∂t= (W5(L− rj , t)P (L− rj , t))−W5(L, t)P (L, t)) +

+ (W6(L− rj , t)P (L− rj , t))−W6(L, t)P (L, t)) +

+ (W7(L− rj , t)P (L− rj , t))−W7(L, t)P (L, t)) +

+ (W8(L− rj , t)P (L− rj , t))−W8(L, t)P (L, t)) ,

(3.31)

whose mean field behaviour is described by the equation

x2 = rx2(1− x2

Lmax)− d0x2 − aLx2. (3.32)

We can proceed as in Section 3.3.3,to obtain an asymptotic approximation of theextinction time for the reduced system Eq. (3.31). For simplicity x and p will denote x2

and p2 respectively.

Antigen-stimulation 53

Figure 3.2: Average of the ratio W3/W5 computed with 1000 Gillespie simulations, fordifferent values of aL.

HLB(p, x) = r(p2 − p)x+r

2Lmax(1− p2)x2 + (aL + dL)(1− p)x (3.33)

The associated Hamilton equations are:

dx

dt=

∂HLB

∂p= r(2p− 1)x− r

Lmaxpx2 − (aL + dL)x

dp

dt= −∂HLB

∂p= −r(p2 − p)− r

Lmax(1− p2)x− (aL + dL)(1− p) (3.34)

3.3.5 Multi-scale stochastic simulations

In order to test the accuracy of our asymptotic analysis, we need to carry out numericalsimulations of the stochastic model. The direct simulation method proposed by Gillespiebecomes very inefficient for systems with separation of time scales. To address this issue,we use the methodology proposed by Cao et al. [20] to perform faster simulations. Adetailed explanation of this method is presented in Appendix D.

Before proceeding with the multi-scale stochastic simulations, we notice that, in thisparticular case, the number of healthy cells, x1, is big enough so that their dynamics areessentially determined by the mean-field dynamics (i.e. we can assume there is no noise).

The main idea of the multi-scale simulation method is to produce stochastic simulationsfor the slow variables and assume that the fast variables are in “stochastic equilibrium”,that is, in an state in which birth and death rates are balanced.

How to decide which are the fast and the slow reactions, in general, is an open question,and usually the best way to determine them is to run few simulations of the Gillespiealgorithm and decide. However, in this case we can use the results of Section 3.3.2. In thissection we consider that x2 is the slow variable, x3, x4 are fast variables and, as before, x1

follows the ordinary differential equation (2.2). The algorithm is summarised as follows.

54 Antigen-stimulation

1. Set initial conditions at t = t0, x1 = x1(t0), x2 = x2(t0), x3 = x3(t0), x4 = x4(t0).

2. Choose the next step size, ∆t, and reaction channel, rj , according to the Multi-scalesimulation algorithm for x2, x3, x4.

3. Update x1, which is evolving according to Equation (3.2), with x2, x3, x4 constantsfrom t to t+ ∆t.

4. Update t = t+ ∆t and x2, x3, x4 according to rj .

5. If t < tf go to 2.

Note that x2, x3 and x4 are considered constant in step 3 because the next reactionchannel does not occur until t+ ∆t.

3.4 Results

In this section we show how the probability of extinction and the average extinction timeof the infection changes as the activation rate of the latently infected cells varies. Weanalyse the system behaviour in two different regimes (subcritical and supercritical), viathe different methods explained in the previous section. Finally, we study a possible sideeffect of the therapy.

3.4.1 Subcritical case

In this section we show the results for the subcritical case, that is, when a metastablestate exists, that is aL < a∗L, using the different methods explained in Section 3.3.

Stochastic Simulation Methods

Performing simulations of this system by means of the regular stochastic simulation algo-rithm [47] requires a huge amount of computational time. In order to reduce the compu-tation time, we propose two different numerical methods. One method consists of using amulti-scale stochastic simulation method. An alternative method is to perform Gillespiesimulations in a reduced 1-dimensional system, in which we have neglected the effect of thevirus infection in the latently infected cells. This second method provides a lower bound ofthe average extinction time (see Section 3.3.4). In Fig 3.3 we compare the average extinc-tion time obtained from these two methods. Good agreement between both methods isobserved, which provides evidence of the suitability of the multi-scale simulation methodto perform simulations in this problem. We use this method, as it allows us to reach awider range of values of aL. The extinction times obtained from the lower bound appearsto underestimate slightly the results from the other methods. This difference increases asaL is close to a∗L, since the metastable state for the reduced systems disappears earlier.

Weak vs Strong QSSA

The probability of extinction at a given time is given by G(0, t), in previous section wehave derived an expression for the semi-classical approximation of G, using it we canapproximate G(0, t) by

G(0, t) = 1− exp (−t/τ) (3.35)

Antigen-stimulation 55

100

1000

10000

100000

0.193 0.194 0.195 0.196 0.197 0.198

Ω C

aL

SSA

multi-scale SSA

LB

Figure 3.3: Average over 1000 simulations of the extinction time of the system describedby Eq. (3.1) computed with the Gillespie SSA, the multi-scale SSA and Gillespie SSA forthe Lower Bound system (Eq. (3.31)) for different values of aL.

with

τ = A√

ΩeΩC . (3.36)

where ΩC is the integral of the re-scaled Lagrangian action, S, along the heteroclinicconnection between the stochastic extinction and the metastable state. Computing thisheteroclinic connection is, in general, very hard. To reduce the dimensionality of theproblem, we have performed a QSSA on the Hamilton equations. In Figure 3.4 we compareΩC for different values of aL computed with the Weak and the Strong QSSA, as aLincreases the difference between both tends to 0, then for large values of aL we can use theStrong QSSA. In Figure 3.5 we show the heteroclinic connection, projected on the (p2, q2)plane, computed with the Strong QSSA and the Weak QSSA.

Lower Bound

Using the reduced system described by Eq (3.31), the computation of the heteroclinic istrivial. In this case the metastable state has the form

(p, x) =

(1, Lmax

r − aL − dLr

), (3.37)

and the stochastic extinction is

(p, x) =

(aL + dL

r, 0

). (3.38)

56 Antigen-stimulation

200

400

600

800

1000

1200

1400

0.1 0.11 0.12 0.13 0.14 0.15 0.16

Ω C

aL

Strong QSSA Weak QSSA

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.1 0.11 0.12 0.13 0.14 0.15 0.16

Diffe

ren

ce

aL

Difference Weak QSSA and Strong QSSA

Figure 3.4: ΩC computed with the Strong QSSA and the Weak QSSA (left), a differencebetween them (right), for different values of aL.

500

1000

1500

2000

2500

3000

3500

4000

4500

5000

0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1

q2

p2

Strong QSSA Weak QSSA

Figure 3.5: Projection on the (p2, q2)−plane of the heteroclinic connection between themetastable state and the stochastic extinction, computed with the Strong QSSA and theWeak QSSA. aL = asL.

Antigen-stimulation 57

1000

10000

100000

1e+06

1e+07

0.19 0.191 0.192 0.193 0.194 0.195 0.196

τ (

days)

aL

Multi-scale SSAsemiclassical Strong QSSA

semiclassical LB

Figure 3.6: Mean extinction time for different values of aL obtained via the semi-classicalapproximation with the Strong QSSA, the semi-classical approximation for the reduced(Lower Bound) system. Average is done over 1000 realisations of the multi-scale stochasticsimulation algorithm

From Eq. (3.33), using the fact that the energy is preserved (H = 0), we obtain ananalytic expression for the heteroclinic connection

x(p) = 2Lmaxr(p2 − p) + (aL + dL)(1− p)

r(p2 − 1), (3.39)

therefore

ΩCLB =

∫ aL+dLr

1x(p)dp =

∫ aL+dLr

12Lmax

r(p2 − p+ (aL + dL)(1− p)r(p2 − 1)

dp

= 2Lmaxrp− log(p+ 1)(r + aL + dL)

r

∣∣∣∣aL+dL

r

1

. (3.40)

In Figure 3.6 we show a comparison between the mean extinction time, and the proba-bility of extinction at a given time computed via stochastic simulations, the semi-classicalapproximation with the Strong QSSA and the semi-classical approximation for the re-duced (lower bound) system. We observe good agreement between the semi-classical andthe exact solution, except in a neighbourhood of a∗L where the number of latently infectedcells is given in a small number and, in this case, the system size is not large enough,therefore the semi-classical approximation is not accurate.

58 Antigen-stimulation

3.4.2 Supercritical case

We are interested in the probability of extinction in a certain time, t, which is given byG(0, t). Of particular interest is the case that the initial condition is the metastable statewith aL = asL.

For aL ≥ a∗L the phase diagram is sketched in Figure 3.1 left. The main feature is thatthere is no metastable state. We approximate G(p, t) ≈ exp(−S(p, t)). A full derivationof this approximation is given in Appendix B.

Given a time Tf , the action is given by

S(p, Tf ) = HTf −∫ Tf

0(q1, q2, q3, q4)(p1, p2, p3, p4)dt (3.41)

−(n1(0), n2(0), n3(0), n4(0))(log (p1(0)), log (p2(0)), log (p3(0)), log (p4(0))),

where ni(0) = pi(0)qi(0) is the initial condition. To obtain the probability of extinction,P0(t) = G(0, t), we proceed as follows:

• Choose initial conditions (p(0), q(0)) satisfying pi(0)qi(0) = ni(0),

• Integrate the equations (3.26) up to p2(T ) = 0, from this we will obtain a time Tand a trajectory γ.

• Compute Sγ the action S along that trajectory γ.

• Then, approximate P0(Tf ) = G(0, Tf ) ≈ exp(−Sγ(0, Tf )).

A different way to approximate G(0, t) should be to study only the extinction timeof the latently infected cells, and assume that the dynamics of the other variables aredescribed by the mean-field behaviour. This assumption is equivalent to consider p1 =p3 = p4 = 1 and look for G(1, 0, 1, 1, t).

The same procedure is applicable for the 1-dimensional system described by Eq. (3.31).Figure 3.7 shows the probability of extinction as a function of time, for different values

of aL, computed according to these methods. For comparison purposes, we also representthe results obtained by multi-scale stochastic simulations. In this case, the semi-classicalapproximation provides values of the extinction probability slightly below the numericalvalues.

3.4.3 Side effects: viral blips

Any increment of aL reduces the viral load of the metastable state, however this does notgive us information on the short time. In this subsection we study the short time effectof the antigen stimulation on the viral load and we will show that this stimulation couldproduce viral blips [110].

We observe that for a short time after the antigen stimulation, the viral load increases(Fig. 3.8). This is because the more x2 cells are activated, the more x3 cells appear in thesystem, therefore the amount of virus, x4, increases as well.

Starting at the metastable state for aL = asL = 0.1 day−1, we perform different simu-lations for different values of aL > asL. We observe, see results in Fig 3.8 (left), that forany aL there is a peak of the viral load peak and then decays. The mean field behaviouralso shows the same evolution pattern, as can be seen in Fig. 3.8 (right), and it gives usa good approximation of the solution of (3.1). In all our tests cases the peak in viral loadstays within a few times the detection limit of standard assays (50 RNA copies/ml).

Antigen-stimulation 59

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 500 1000 1500 2000 2500 3000

Pro

bab

ility

time (days)

SSASemi-classial Strong QSSASemi-classical + meanfield

Figure 3.7: Probability of extinction as a function of time comparing three different ap-proximations: assuming p1 = p3 = p4 = 1(blue), Multi-scale stochastic simulations (red),semi-classical approximation for the QSSA (green)

(a) (b)

50

100

150

200

250

0 2 4 6 8 10 12 14 16 18 20

Viral lo

ad

days

aL=0.1 aL=1 aL=100 aL= ∞

50

100

150

200

250

300

0 2 4 6 8 10

Viral lo

ad

time (days)

SSA mean-field

Figure 3.8: Gillespie simulations of the viral dynamics for different values of aL. As initialcondition we take the value of the variables at the metastable state before the antigenstimulation, that is aL = asL (left). A comparison between a Gillespie simulation and themean-field behaviour of the viral dynamics. As initial condition we take the value of thevariables at the metastable state before the antigen stimulation, that is aL = asL, and wechoose aL = 100. (right)

60 Antigen-stimulation

3.5 Conclusions

In this chapter we have presented and analysed a stochastic model of HIV-1 dynamics inthe blood stream under anti-retroviral therapy to study of the survival probability in acertain time of the viral load. We have investigated how the extinction probability dependon the activation rate of the latently infected cells using stochastic simulations and analyticapproximations. Furthermore we have presented some reductions of the system which hasgiven a lot of tractability to the problem without losing accuracy in the solution. Froma medical point of view, according to this model, this therapy can reduce the averageextinction of the infection to a few years.

Chapter 4

Stochastic modelling of viral blipsin HIV-1-infected patients: Effectsof inhom*ogeneous densityfluctuations

An stochastic model of HIV-1 infection dynamics under HAART is proposed in order toanalyse the origin and dynamics of the so-called viral blips, i.e. episodes of transientviremia that occur in the phase of where the disease remains in a latent state duringwhich the viral load raises above the detection limit of standard clinical assays. Based onprior work in this subject, we consider an infection model in which latently infected cellcompartment sustains a residual (latent) infection over long periods of time and, unlikeprevious models, we include the effects of inhom*ogeneities in cell and virus concentration inthe blood stream. We further consider the effect of burst virion production. By comparingwith the experimental results obtained during an study in which intensive sampling wascarried out on HIV-1-infected patients undergoing HAART over a long period of time,we conclude that our model supports the hypothesis that viral blips are consistent withrandom fluctuations around the average viral load. We further conclude that agreementbetween our simulation results and the blip statistics obtained in the aforementioned studyimproves when burst virion production is considered. We also study the effect of samplemanipulation artifacts on the results produced by our model, in particular, that of thepost-extraction handling time, i.e. the time elapsed between sample extraction and actualtest. Our results support the notion that the statistics of viral blips can be criticallyaffected by such artifacts.

4.1 Introduction

An early hypothesis regarding the origins of viral blips, whose emergence was, at thetime, suspected of heralding imminent virological failure, was that they are due to theappearance of new, drug-resistant viral variants [80]. However, a number of studies havecompiled evidence against viral blips being correlated with virological failure [54, 89, 95],thus weakening the case of viral blips being early warnings related to viral evolution leadingto drug resistant strains of the virus. Further to the hypothesis of blips being originated by

61

62 Viral Blips

viral evolution, other mechanisms have been proposed which include antigen-driven CD4+T cell activation due to vaccination or secondary infections [37, 52, 41, 40, 43, 61, 62]. It hasalso been shown that activation of latently infected cells may play a role in the emergenceof viral blips. For example, Jones & Perelson have proposed a model in which increasedactivation of latently infected cells can lead to a burst in viral load [63]. The possibilitythat asymmetric division of activated latently infected cells may help to explain the decaykinetics of the latent compartment and intermittent viral blips has been in explored in[109]. Recently, Rong & Perelson [111] have formulated a model in which stochasticactivation of latently infected cells can maintain viral blips without completely depletingthe latent reservoir, thus maintaining long-term, low-level viremia. They also developeda model which incorporated density-dependent homeostatic proliferation of the memoryCD4+ T memory cells (and, therefore, of the latently infected cell compartment), which,according to [23], drives persistence and determines the size of the latent reservoir.

Alongside all these models which postulate that viral blips have a physiological origin,there is an alternative school of thought which claims that most viral blips are randomoccurrences of probabilistic origin related to the small number of virus mRNA in the latentphase as well as being partly produced by laboratory artefacts during the processing of thesamples [76, 94]. According to this view, viral blips are mostly uncorrelated with virologicalfailure, virus evolution, vaccination, or non-adherence to treatment regimen [95, 76, 94]and, therefore, only a small fraction of such occurrences are of clinical significance. Conway& Coombs [28] have proposed a model to analyse the stochastic viral dynamics in treatedpatients. This model treats viral blips as random events occurring every time the viral loadreaches the standard detection limit (i.e. 50 mRNA/ml). Although this model reproducesmany of the features of HIV-1 viral dynamics and provides a detailed description of itsstochastic dynamics, there are several properties of the statistics of viral blips in which[28] appears to depart from experimental observation. One such departure refers to thefrequency of blips. According to the model formulated in [28], blip frequency decreasesexponentially in time which disagrees with the results by DiMascio et al. [30] where aconstant blip rate is reported. Conway & Coombs [28] claim that this difference is dueto the fact that, whilst theirmodel accounts only for blips of small amplitude causedby fluctuations, [30] inlcude in their statistics blips of both small and large amplitude.Therefore, Conway & Coombs implicitly assume that large blips must have a biologicalorigin. Similarly, the data collected by Nettles et al. [95], where intensive sampling wascarried out on HIV-1-infected patients with samples collected every 2-3 days over a periodof 3 to 4 months, reported that blips were observed in 9 out of 10 patients with an averageof two blips per patient. The blip statistics produced by the model formulated in [28]appears to predict that blips be too rare to be comparable to the experimental data, evenin conditions of elevated virus rate production and latent cell activation. It can be arguedthat the rationale for this discrepancy is the same as before: Nettles et al. [95] measurethe full range of blips (both small and large amplitudes) whilst Conway & Coombs [28]account only for small-amplitude blips. Therefore, according to this model [28], smallamplitude blips are consistent with an stochastic model whereas large amplitude blipsmust be produced by causes other than random fluctuations.

One common feature shared by all of the models discussed in the Introduction, bothstochastic and deterministic, is that they all assume that the system is well-mixed, i.e.all the cellular and molecular species are evenly distributed over the volume of bloodso that only the number of each species determines the state of the system. For this

Viral Blips 63

modelling assumption to hold, the numbers of each species must be large enough so thattheir densities remain approximately uniform. The system we are dealing with in thischapter has at least two species, virus and latently infected cells, which are present invery small numbers. As a consequence, their densities can show large fluctuations whichwould lead to large inhom*ogeneities in the local numbers of such species. The effect ofthis inhom*ogeneity can be rather sizeable, specially since measurements of viral load areperformed by extracting small samples of blood which are then analysed.

The aim of this chapter is to ascertain whether density fluctuations affect the stochasticdynamics of the viral load in HAART-treated patients beyond the predictions of reference[28] and investigate if an stochastic model which includes density fluctuations is capable ofa more faithfull reproduction of the experimental results reported in [30, 95]. To this endwe formulate an stochastic compartmental model [13] where volume of blood is divided insmall compartments (whose volume is assumed of the order of the volume of one bloodsample taken for analysis). Within each of this compartment, the system is assumedto be well-mixed but we consider that differences may arise between compartments thusaccounting for density fluctuations [46, 34, 35, 36, 79, 117]. Furthermore, in order toaddress the issue of whether laboratory artefacts alter the statistics of viral blips we havedesigned an in silico blood extraction procedure that allows us to study, at least partially,this aspect of the problem.

This chapter is organised as follows. Section 4.2 is devoted to a detailed explanation ofthe model assumptions and formulation. In Section 4.3 we conduct extensive simulationsof our model and present our results. Finally, in Section 4.4 we proceed to discuss ourresults.

4.2 Model formulation

We now proceed to present our compartmental model of stochastic viral dynamics. Ourmodel is based on the modelling approach whereby spatially inhom*ogeneous systems aredealt with by compartmentalising the domain where the system lives into small compart-ments [13]. Within each of these compartments all of the species participating in thedynamics of the system are assumed to be well-mixed, so that our resolution to measurespatial variations is of the order of the compartment size. The local, well-mixed, within-compartment stochastic dynamics is supplemented by stochastic transitions that allow fortransport of species between compartments. This approach has been widely applied tomany situations in which spatial heterogeneity is essential, in particular, reaction-diffusionsystems [46, 34, 35, 36, 79, 117]. This setting will allow us to ascertain the effect of densityfluctuations on stochastic viral dynamics.

Further to our compartmental model, we also formulate an in silico blood analysismodel, which intends to follow as closely as possible the intensive sampling proceduredesigned by Nettles et al. [95] to analyse the dynamics and statistics of viral blips inHIV-1-infected patients under HAART.

4.2.1 Compartmental Model

In order to account for density fluctuations, we consider a compartmental model, schemat-ically shown in Fig. 4.1. We will assume that, in an attempt to roughly imitate thestructure of the blood stream, compartments are arranged on a one dimensional, closed

64 Viral Blips

Figure 4.1: This figure depicts an schematic diagram of our compartmental model. Ourmodel represents the circulatory system in a very simplistic way as a one dimensional closedcircuit of adjacent compartments. Within each compartment we model the (local) cellulardynamics of an HIV-1-infected T-cell population with latency [110, 28]: We consider thedynamics of the HIV-1 virus load (V ), active infected T-cells (T ∗), and latently infectedcells (L). Virus capsides, V , are released by active infected cells and infect uninfectedT-cells. We consider that the number of uninfected T-cells is so much larger than thatof infected cells that we consider them as a reservoir. Latently infected cells can becomeactive upon estimulation with proper antigens. In addition to this within-compartmentcellular dynamics, we consider random movement of both latently and active infected cellsand virus capsides between compartments.

Viral Blips 65

Transition rate rj+(i−1)RL= (∆Li,∆T

∗i ,∆Vi) Description

W1+(i−1)RL= η(1− ε)kViTiΩ−1 (1, 0,−1) Latent infection Ti + Vi → Li cell

W2+(i−1)RL= (1− η)(1− ε)kViTiΩ−1 (0, 1,−1) Active infection Ti + Vi → T ∗i cell

W3+(i−1)RL= rLi (1, 0, 0) Latent cell proliferation Li → 2Li

W4+(i−1)RL= r

Lmax

Li∗(Li−1)2 Ω−1 (−2, 0, 0) Latent cell binary annihilation Li + Li → ∅

W5+(i−1)RL= d0Li (−1, 0, 0) Latent cell clearance Li

W6+(i−1)RL= aLLi (−1, 1, 0) Activation Li → T ∗i

W7+(i−1)RL= δT ∗i (0,−1, 0) Active cell death

W8+(i−1)RL= cVi (0, 0,−1) Virion clearance

W9+(i−1)RL= pvT

∗i (0, 0, 1) Continuous virion production

W10+(i−1)RL= εkViTiΩ

−1 (0, 0,−1) Failed infection

Table 4.1: Transition rates corresponding to the stochastic model of viral blip generation.See text (Section 4.2) for details. RL = 10. Ti, the number of uninfected cells is taken tobe constant [28]. Ω is the compartment size.

(i.e. with periodic boundary conditions) lattice (see Fig. 4.1). In each compartment, weconsider three types of interacting species, namely, active infected cells, latently infectedcells, and virus. The number of each of these species in compartment i is referred to as T ∗i ,Li, and Vi, respectively, where i = 1, . . . , Nc with Nc is the number of compartments. Wefurther introduce the vector X(t) ≡ (T ∗1 (t), L1(t), V1(t), . . . , T ∗Nc

(t), LNc(t), VNc(t)) whosecomponents correspond to the number of cells of of each species at site i = 1, . . . , Nc andtime t. We will find convenient to define xi(t) ≡ (T ∗i (t), Li(t), Vi(t)), so that the vectorX(t) can be written as X(t) = (x1(t), . . . , xNc).

The stochastic dynamics of our system is described by the corresponding Master Equa-tion [122, 79, 44]:

∂

∂tP (X, t) =

Nc∑i=1

RL∑j=1

(Wj+(i−1)RL

(X − rj+(i−1)RL)P (X − rj+(i−1)RL

)−Wj+(i−1)RL(X)P (X)

)Nc∑i=1

RT∑j=1

(ωj+(i−1)RT

(X − ρj+(i−1)RT)P (X − ρj+(i−1)RT

)− ωj+(i−1)RT(X)P (X)

)where we have splitted the Master Equation in a purely local part, corresponding to thepopulation dynamics within each compartment, and a transport part, which accounts fortransport of species between adjacent compartments. RL and RT are the number of localand transport events, respectively, that can affect the state, xi, of each compartment. Thetransition rates are defined in Tables 4.1 and 4.2.

In order to facilitate later formulation of our in silico blood extraction protocol and itscomparison with [95], we will consider compartments of volume Vc = 8.5 ml each, which isthe volume of blood sampled extracted for analysis in the study reported in [95]. We willconsider that an average individual has 5 litres of blood so we need to consider Nc = 588compartments.

The within-compartment stochastic viral dynamics model is based on previous worksby Rong & Perelson [110] and Conway & Coombs [28] and it is schematically shown inFig. 4.1. Following the latter, we assume that, due to the large numbers in which CD4+

66 Viral Blips

Transition rate ρj+(i−1)RTDescription

ω1+(i−1)RT= µL+Li (0, ..., 0,−1, 0, 0, 1, 0, ...) Li + Li+1 → Li − 1 + Li+1 + 1

ω2+(i−1)RT= µL−Li (0, ..., 0, 1, 0, 0,−1, 0, ...)) Li + Li−1 → Li − 1 + Li−1 + 1

ω3+(i−1)RT= µT ∗+T

∗i (0, ..., 0,−1, 0, 0, 1, 0, ...)) T ∗i + T ∗i+1 → T ∗i − 1 + T ∗i+1 + 1

ω4+(i−1)RT= µT ∗−T

∗i (0, ..., 0, 1, 0, 0,−1, 0, ...) T ∗i + T ∗i−1 → T ∗i − 1 + T ∗i−1 + 1

ω5+(i−1)RT= µV+Vi (0, ..., 0,−1, 0, 0, 1, 0, ...) Vi + Vi+1 → Vi − 1 + Vi+1 + 1

ω6+(i−1)RT= µV−Vi (0, ..., 0, 1, 0, 0,−1, 0, ...) Vi + Vi−1 → Vi − 1 + Vi−1 + 1

Table 4.2: Transition rate corresponding to the stochastic model of viral blip generation.A detailed is given in the main text (Section 4.2). RT = 6.

T uninfected cells are present, they are unaffected by fluctuations and their number istaken as constant. We assume that two types of infected cells exist, latenly infected cells,Li, which carry the virus but do not synthesize new virions, and active infected cells, T ∗i ,which produce and release new virions. As a consequence, active infected cells are targetedby HAART whereas latently infected cells are immune to its effects. Both types of infectedcells are produced by infection of uninfected T cells. Furthermore, latently infected cellscan become active, for example, by estimulation with appropriate antigens. Both latentlyand active infected cells are assumed to die at a certain rate and blood-borne virions areassumed to be cleared off at a constant rate. We also account for experimental resultsshowing that the latently infected cells are maintained by homeostatic proliferation [23],by means of a phenomenological model which includes branching and binary annihilationof latently infected cells [110, 35].

The elementary events involved in our spatial model of stochastic viral dynamics areas follows.

• Latenly infected cells Li can undergo:

1. Homeostatically balanced proliferation. Following [110], we account for ho-moeostatic control of proliferation by means of a combination of branching(Li → 2Li) with binary annihilation (Li+Li → ∅). It has been shown (see e.g.[35]) that this combination is an stochastic counterpart of the standard logisticgrowth. The corresponding transition rates are W3+(i−1)RL

for branching andW4+(i−1)RL

for binary annihilation (see Table 4.1).

2. Death. We assume a simple linear decay with transition rate W5+(i−1)RLas

shown in Table 4.1.

3. Activation. By means of this process a latently infected cell becomes an activeinfected cell Li → T ∗i . The corresponding transition rate is W6+(i−1)RL

, Table4.1.

4. In addition to the processes described above, which have to do with the (within-compartment) population dynamics, we assume that latently infected cells canmove between neighbouring compartments. We will assume a coupling betweencompartments where transitions occur at constant per cell rate, νL. In periodic,one dimensional setting, cells in compartment i can move to compartment i+1at rate ω1+(i−1)RT

or to compartment i − 1 at rate ω2+(i−1)RTas defined in

Table 4.2.

• Active infected cells, T ∗i , cell are subjected to:

Viral Blips 67

1. Apoptosis with transition rate W7+(i−1)RL(see Table 4.1).

2. Virion production. Contrary to latently infected cells, active infected cellssynthesise and release new virions. In general, viral production can occur in acontinuous fashion over the life span of an infected cell or in a burst which killsthe cell. For the HIV infection both modes have been proposed [97]. Transitionrate W9+(i−1)RL

(Table 4.1) corresponds to continuous production. Later inthe chapter, we consider an additional scenario, in which both continuous andburst production occur.

3. Transport between neighbouring compartments. Just as latently infected cellsdo, active infected cells can move between compartments at rates ω3+(i−1)RT

and ω4+(i−1)RTas per Table 4.2.

• Finally, virus, Vi, can:

1. Infect a healthy T cell producing a latently infected cell. In patients underHAART, the infection process is hindered by the presence of an antiretrovi-ral drug. The efficiency of such drug is measured by a parameter, ε, whichtakes values between 0 and 1, the latter (former) corresponding to a maximally(in)efficient drug. (1− ε) is interpreted as the proportion of virions capable ofinfection under HAART treatment. We also assume that, upon infection, thecell can become latently infected with probability η or active with probability(1−η). Therefore the corresponding transition rate W1+(i−1)RL

is proportionalto η(1− ε) as shown in Table 4.1.

2. Infect a healthy T cell producing an active infected cell. In this case, thecorresponding transition rate W2+(i−1)RL

is proportional to (1− η)(1− ε) (seeTable 4.1).

3. Undergo clearance. Virions are removed from the blood, which we model as asimple linear decay with transition rate W8+(i−1)RL

as per Table 4.1.

4. Fail to infect and being eliminated by the drug with transition rate W10+(i−1)RL

(see Table 4.1).

5. Transport between neighbouring compartments. Virions can move betweencompartments at rates ω5+(i−1)RT

and ω6+(i−1)RTas per Table 4.2).

Numerical methodology We will perform numerical simulations of the Master Equa-tion 4.1 using the methodology proposed by Bernstein [13], which is a straight forwardgeneralisation of the original stochastic simulation algorithm proposed by Gillespie [47].

Parameter values The (default) parameter values used in our simulations are givenin Table 4.3. The parameter values corresponding to the model of cellular dynamics ofan HIV-1-infected T-cell population with latency are based on estimates available in theliterature on the subject.

Random motility between compartments has been modelled as simple diffusion. Thetransition probability between compartments is therefore given by µJ+ = DJ

h2, where J =

V, L, T ∗, DJ is the diffusion coefficient of species J , and h = 3√Vc, where Vc is the

compartment volume. The magnitude of the diffusion coefficient for several types of viralcapsides has been estimated between 1.6 and 30 µm2s−1 [93]. We consider a typical

68 Viral Blips

Parameter Rate Description Reference

λ 10000 ml−1 day−1 Recruitment rate of T cells [19]

dT 0.0166 day−1 Death rate of T cells [90]

k 2.4 · 10−8 ml day−1 infection rate [100]

ε 0.85 Drug efficacy [110]∗

η 0.001 Fraction resulting in latency [63]

d0 0.001 day−1 Death rate of latently infected cells [19]

aL 0.1 day−1 Rate of transition from latently to productively [110]

σ 1 day−1 death rate of productively infected cells [83]

c 23 day−1 clearance rate of free virus in blood stream [104]

pv 2000 day−1 Viral production rate [59]

r 0.2 day−1 proliferation rate of activated cells [110]

Lmax See caption Carrying capacity density of latent cells [110]

Vc 8.5 ml Compartment volume [95]

Table 4.3: Parameter values used in our numerical simulations. The value of Lmax is letto vary depending on the average viral load we impose on our simulations. For the valuesof the average viral load considered in Section 4.3, i.e. 10, 12.5, 20, 30, and 40 copies/ml,the corresponding values of Lmax = 1.89, 2.36, 3.78, 5.66, and 7.55 cells/ml, respectively.

value DV = 5 µm2s−1 [56, 12]. Similarly, we take a generic estimate for the diffusioncoefficient of a cell DL = DT ∗ = 0.05 µm2s−1 [15]. Furthermore, in order to account forthe directionality of blood flow, we assume that µJ− = 3µJ+.

4.2.2 In silico blood sample analysis model

Further to the compartmental model of Section 4.2.1 which allows us to account for densityfluctuations in stochastic viral dynamics, we have formulated an In silico blood sampleanalysis model. This additional model has the aim of trying to reproduce the experi-mental procedure described in [95] as closely as possible. Nettles et al. [95] designed anexperimental protocol in which ten HIV-1 infected patients under HAART were inten-sively sampled (every 2 or 3 days) for a period of 3 to 4 months. Each time blood wasextracted, two samples per patient were taken and sent to two different laboratories foranalysis. They concluded that blips had short duration (less than 3 days on average) andlow amplitude (79 copies/ml on average) with an average of 1.8 blips per patient.

Our procedure is as follows. After running the compartmental dynamics until timetc, which is chosen to be long enough so that the average properties of the system reachan steady state, we choose two compartments at random i and j among the Nc com-partments that compose our system. We then record the corresponding state xi(tc) =(T ∗j (tc), Li(tc), Vi(tc)) and xj(tc) = (T ∗j (tc), Lj(tc), Vj(tc)). Recall that we are assumingthat the number of uninfected T cells is assumed to be constant and uniform [28]. In or-der to account for possible delays between the time of extraction and the actual analysis,we assume that the extracted samples continue to evolve subject to the local (within-compartment) dynamics determined by the rates in Table 4.1 with the transition ratescorresponding to between-compartment transitions ωl(X) = 0. This post-extraction dy-namics is ran for a duration tw.

Viral Blips 69

Figure 4.2: Probability of observing at least one blip in our in silico blood sample analysismodel as a function of the average virus load for different values of the burst size. Nettleset al. [95] in their intensive sampling study showed that 9 out 10 patients underwent atleast one viral blip. We observe that, as the average viral load increases the probability ofobserving at least one blip tends to one. We also observe that as the burst size increasesthe observation of at least one blip becomes more likely. Note that the rate of continuousvirion production has been chosen for each value of the burst size so that the averagenumber of virions produced per active infected cell during its life-time is kept constant.We have fixed an average viral load of 12.5 copies/ml [95]. Parameter values as per Table4.3. p′v is taken such that p′v + δN = pv, where N is the burst size and pv is given in Table4.3.

70 Viral Blips

Figure 4.3: Comparison between blips statistics obtained in the study by Nettles et al.[95] and those obtained from our simulations for different burst sizes. This plot showsthe probability of observing a blip of a given duration (measured in days). The rate ofcontinuous virion production has been chosen for each value of the burst size so that theaverage number of virions produced per active infected cell during its life-time is keptconstant. We have fixed an average viral load of 12.5 copies/ml [95]. Parameter values asper Table 4.3. p′v is taken such that p′v + δN = pv, where N is the burst size and pv isgiven in Table 4.3.

Viral Blips 71

Figure 4.4: This plot shows the average number of blips and median amplitude as afunction of the burst size obtained with our model. According to Nettles et al. [95] ,theaverage number of blips per patient is 1.7 and the median of the amplitude is 79 copies/ml.The rate of continuous virion production has been chosen for each value of the burst sizeso that the average number of virions produced per active infected cell during its life-timeis kept constant. We have fixed an average viral load of 12.5 copies/ml [95]. Parametervalues as per Table 4.3. p′v is taken such that p′v + δN = pv, where N is the burst size andpv is given in Table 4.3.

72 Viral Blips

4.2.3 Effects of continuous versus burst production of virions

The issue of whether virions are secreted by active infected cells in a continuous way duringtheir lifespan or, on the contrary, they are released in a single burst which, actually, kills thecell remains controversial. Standard models based on mean-field, deterministic systemsof ordinary differential cannot distinguish between both production modes, since theypredict exactly the same dynamics as long as the total number of virions produced overtheir lifespan is the same [97]. A recent analysis by Pearson and co-workers has concludedthat stochastic models of early infection lead to significantly different results dependingon whether virion production is considered as continuous or bursty [97].

Here, we consider this issue regarding its possible influence on the statistics of viralblips. To this end, we consider two models, namely, one in which only continuous virionproduction is taken into account (given in Table 4.1) and another one in which bothmodes of virion production, continuous and bursty, are considered. The latter requires analteration of the process described by Table 4.1. The reaction corresponding to active celldeath now reads:

W7+(i−1)RL= δT ∗i and r7+(i−1)RL

= (0,−1, N) (4.1)

which implies that upon active cell death, one active cell death is lost and N virions arerealeased. N is the so-called burst size. A further modification needs to be done, this timeaffecting the rate of continuous virion production:

W9+(i−1)RL= p′vT

∗i and r9+(i−1)RL

= (0, 0, 1) (4.2)

where p′v, i.e. the rate of continuous virion production, is now a function of the burstsize, N . Given a value of N , we fix p′v so that the average number of virion producedduring the lifespan of am active infected cell is the same as for the purely continuousvirion production model (i.e. N = 0 as per Table 4.1).

The remaining reactions are left unmodified with respect to those shown in Table 4.1.

4.3 Results

In this section we will use the models presented in section 4.2 to explore the phenomenologyassociated to our model and obtain statistics of viral blips which can then be comparedto the experimental data reported by Nettles et al. [95].

4.3.1 Effect of density fluctuations on blips statistics

We have carried out extensive numerical simulations of our model using the methodologyexplained in Section 4.2. We start by focusing in our inhom*ogeneous infection modelwith continuous virion production as defined by the transition rates given in Tables 4.1,accounting for the local infection dynamics, and 4.2, accounting for diffusion of cells/virionsbetween compartments. Results are included in Figs. 4.2, 4.3 and 4.4.

Nettles et al. [95] reported that 9 out 10 patients who underwent their intensive sam-pling protocol were observed to exhibit at least one blip. In other words, the probabilityof detecting at least one blip during the course of their experiments can be estimated tobe 90%. By means of numerical simulations of our in silico blood sample model, we havecomputed the probability of detecting at least one blip as a function of the average viral

Viral Blips 73

load (see Fig. 4.2, burst size N = 0). Our results show that, as the average viral loadgrows, so does the probability of detecting at least one blip. According to our resultsfor viral loads bigger that 20 copies/ml, this probability approaches 1, in good agreementwith figure reported in [95].

Nettles et al. [95] reported further statistical information, in particular, regardingblip duration (reproduced in Fig. 4.3) by recording the frequency with which blips of acertain duration are observed. They observed the blips with duration of 3 days or shorteraccounted for over 80% of the observations. Blips of duration between 3 and 5 days, 5and 7 days and of 7 days or longer were observed in less than 10% of the cases each.Our simulation results are in good agreement with these observations (see Fig. 4.3). Ourmodel, however, appears to overestimate the frequencies of the shorter blips (duration of3 days or shorter and between 3 and 5 days) at the expense of the longer ones.

Further information that can be extracted from the data reported in [95] concerns theaverage number of blips and their average amplitude. According to their experiments, theaverage number of blips is 1.8 and their average amplitude 79 copies/ml. Regarding thesetwo metrics, our simulation results reported in Fig. 4.4 for burst size equal to zero showthat our inhom*ogeneous infection model with continuous virion production underestimatesthese quantities by a significant amount.

4.3.2 Effect of burst production of virions

In order to improve the quantitative agreement between the data reported in [95] andour model predictions, we have explored the effect of including both continuous and burstvirion production (see Section 4.2.3). In particular, we have analysed the behaviour ofour system as a function of the burst size, keeping constant the average number of virionsproduced by active infected cells during their lifespan.

We have studied the effect of burst production on the statistics concerning viral blipsshown in Figs. 4.2, 4.3 and 4.4. Regarding the probability of detecting at least one blip,although its behaviour for positive burst size is qualitatively similar to the burst size N = 0case, our model predicts an increase in the probability of observing at least one blip forlower values of the average viral load (see Fig. 4.2). This increase is proportional to theburst size: The bigger the burst size, the larger the increase in the probability of observingat least one blip.

The discrepancies between experimental data and simulation results concerning thefrequency of blip duration tend to be corrected when a positive burst size is considered(see Fig. 4.3). As we have discussed previously, the inhom*ogeneous infection model withcontinuous virion production overestimates the frequency of shorter blips at the expenseof longer ones. This trend is compensated when we consider burst production of virions:As the burst size is let to increase, shorter blips become less frequent at the expense oflonger ones, whose frequency is observed to increase as shown in Fig. 4.3.

Further discrepancies between the data reported in [95] and simulation results of ourinhom*ogeneous infection model with continuous virion production regarding the averagenumber of blips and their average amplitude. Our model with no burst virion productionunderestimates these quantities. Fig. 4.4 shows that, as burst size increases, the valuepredicted by our model for these two quantities approaches those reported in [95].

To summarise, although our inhom*ogeneous infection model with continuous virionproduction appears to qualitatively capture the blip statistics obtained by Nettles et al.

74 Viral Blips

Figure 4.5: Probability of observing n blips within a sample, P (n), as a function of thepost-extraction handling time, tw, and the virus clearance rate, c. We have also considereddifferent value of the average viral load: The first, second, third, and fourth correspondto 10, 20, 30, and 40 copies/ml, respectively. We have considered the continuous virionproduction only.

Viral Blips 75

[95], better quantitative agreement is obtained when burst production of virions is consid-ered.

4.3.3 Effect of post-extraction handling time

A controversial aspect regarding the origin and clinical significance of viral blips is theclaim that their statistics are affected by artifacts related to the sample manipulationprocess [76, 94]. In order to partially address this issue within the context of our model,we have devised a set of simulations in which we consider the effect of the post-extractionhandling time in our in silico blood test model, tw (see Section 4.2.2). This time isassumed to correspond to the time elapsed between extraction of the blood sample andactual analysis at the laboratory. We have considered that tw is of the order several hours.During this time lapse, the infection dynamics stays running but since the sample is nowextracted and therefore isolated there is no between-compartments cell/virion motility(i.e. the dynamics is described by a purely local infection model determined by the ratesgiven in Table 4.1). Extraction of a compartment (sample) and its consequent isolationfrom the blood stream has the likely effect of reducing the rate at which virus particlesare cleared, since it has been isolated from usual removal mechanisms. We have thereforeconsidered the effect of reducing the virus clearance rate with respect to current estimatesin the literature (see Table 4.3). Note that all the simulation results presented in Sections4.3.1 and 4.3.2 correspond to tw = 0.

In order to assess the effects of both the handling time, tw, and the virus clearancerate, c, we have done numerical simulations to compute the Probability of observing nblips within a sample, P (n), for different values of the average viral load. Our resultsare reported in Fig. 4.5. For lower values of the average viral load (see Fig. 4.5, 10 and20 copies/ml), we observe that there is a strong dependence on the average number ofblips with tw: Due to accumulation of virus capsides as time progresses (recall that weare assuming that active infected cells continue to release virions after extraction untilsuch time as the actual count is carried out), the average number of blips increases, asindicated by the fact that the peak of P (n) shifts to larger values of n as tw increases.Likewise, the position of this peak shifts to larger values of n as the virus clearance rate, c,decreases. For larger values of the average viral load ((see Fig. 4.5, 30 and 40 copies/ml),our simulation results show that the peak of P (n) shifts to smaller values of n as twincreases. This result is actually an artifact of how blips are defined and counted: A blipis the event whereby the viral load, initially below 50 copies/ml, transiently grows abovethe detection threshold. When the average viral load is 30 or 40 copies/ml, the durationof these events becomes very long, which has the by-product of reducing the blips count.To illustrate this phenomenon, we show in Fig. 4.6 the number of observations which areabove 50 RNA copies/ml. When the viral load is 10 copies/ml, the difference between thenumber of observations above the detection limit and the number of viral blips is small.However, in the case of 40 copies/ml we see strong discrepancies between both pictures.Although the viral load is most of time above 50 copies/ml, we are considering most ofthat observations as the same viral blip and, in this case, the number of viral blips is givenby the number of observations below the detection limit.

76 Viral Blips

Figure 4.6: Probability of having n observations above the detection limit within a sample,P (n), as a function of the post-extraction handling time, tw, and the virus clearance rate,c. We have also considered different value of the average viral load: The first, and secondcorrespond to 10 and 40 copies/ml, respectively. We have considered the continuous virionproduction only.

Viral Blips 77

4.4 Conclusions

In this chapter we have studied if the inhom*ogeneous density fluctuations, together withexperimental artifacts could be the origin of the viral blips, or if an extra assumptions hasto be included in the model of chapter 3. To study this we have presented a stochasticmodel of HIV-1 dynamics in the blood stream, very similar to the model presented inchapter 3 but supplemented with a spatial dependence. To do this we have considereda compartmental model, in each compartment the same dynamics of the model of theprevious chapter are considered, supplemented with a random diffusion between compart-ments. To determine how the experimental protocols can affect the appearance of viralblips we have presented an in silico blood test model, to reproduce in a realistic way theexperimental protocols and to study if the experimental artifacts can affect the results.Moreover we have shown that the dynamics of the infection in the blood samples, duringthe time between the extraction and the observation, can play a major role in the appear-ance of the viral blips. And we have studied how the different virion production (continuevs burst) affect the appearance, the duration and the magnitude of the viral blips.

78 Viral Blips

Chapter 5

Conclusions

The this PhD Thesis we have studied the effect of fluctuations in complex structuredpopulations characterised by the existence of a latent state. In particular we have put lighton some questions regarding the dynamics of hierarchically organissed cell populations,and we have investigated the dynamics of HIV-1 in infected patients under anti-retroviraltherapy. In this chapter we expose the conclusions of our study.

5.1 Stochastic dynamics of differentiation cascades with reg-ulatory feed-back

The aim of Chapter 2 is to study the robustness of the cell populations with hierarchicalstructure. To study this we have presented and analysed stochastic models of hierarchicalpopulations governed by differentiation cascades to study the effect of symmetric stem-cellself-renewal on their long-time stability properties as well as their resilience to invasion.The behaviour of these models has given us some insight on intrinsic protection mecha-nisms, in particular against invasions by malignant, potentially tumorigenic, populations.

Our model of a hierarchically-organised tissue consists of a stochastic compartmentalmodel, where each compartment corresponds to a differentiation stage, with negative-feedback regulating the size of the stem cell compartment. This feedback is assumed to bemediated by a cytokine whose secretion rate is, in turn, modulated by the number of fully-mature cells. The general framework we have used in our model formulation, i.e. in termsof cellular compartments with feedback between the SC and the MC compartments, hasbeen used in a wide variety of studies, ranging from general models of robustness againstmutations leading to cancer [98, 106, 107, 119] to more specific models of particular tissuessuch as the haematopoeitic system [1, 26, 27, 29, 82], epidermis [112] or neurogenesis [7].Our model formulation is therefore rather generic and the resulting model should beapplicable to a wide variety of tissues where hierarchical structure is present.

However, the issue of the cytokine-mediated feedback has more supporting evidence inthe haematopoeitic system, where several of these cytokines have been identified which,upon variation in the size of the differentiated cell compartment, regulate the size of the SCcompartment by means of negative feedback loop. For example, in the production of whiteblood cells, the granulocyte colony-stimulating factor (G-CSF) has been demonstrated tomediate a negative feedback [26]. Another example is thrombopoietin, which is involvedin negative-feedback regulation in the production of platelets [26]. Erythropoiesis, i.e. the

79

80 Conclusions

production of red blood cells, is yet another example in which a negative feed-back ismediated by a cytokine, namely, erithropoietin [81, 118].

In order to investigate the effects of symmetric self-renewal on the robustness of hier-archical populations, we have first formulated a model in which stem cells undergo onlysymmetric self-renewal. Regarding this model, we have identified a new extinction mech-anism in which the coupling between fluctuations that increase the pool stem cells andthe delay in stem cell regulation leads to the extinction of the stem cell population and,therefore, to the eventual extinction of the whole: fluctuations in the size of the stemcell compartment propagate along and are amplified by the differentiation cascade in anamount which typically scales as 2n, where n is the length of the cascade. Since the reg-ulation of stem cell proliferation by the fully differentiated cells has a delay proportionalto n, the negative feed-back between both populations may lead to extinction of the stemcells and therefore of the whole population.

The extinction time, TE , in such symmetric population has been shown to be sensitiveto variations of the death rate of the mature cells, λn, and the rate of self-renewal, p. Wehave shown that increasing λn leads longer extinction times, whereas populations withlarger values of the self-renewal rate p exhibit shorter extinction times. Biologically, theseresults imply that symmetric stem-cell self-renewal is more likely to be observed in thosewith shorter delays between the stem cell and mature cells and those where the life-time ofthe mature cells is shorter and therefore larger cellular turn over is necessary to maintaintissue homeostasis.

We have also explored the possibility that the cytokine is secreted in response tovariations in the size of some of TACs compartments. We observe that the behaviour issimilar to that of shorter differentiation cascades: If the cytokine is produced at a rateproportional to the number of cells in compartment i, the behaviour of the stem cells is,for all intents and purposes, independent of compartments i + 1, ...m which therefore donot affect the system. We have also seen that if the secretion rate of the cytokine wereproportional to the number of some of the TACs, the corresponding secretion rate wouldhave to be much larger to maintain the same population of stem cells, as the number ofTACs is smaller than the number of MCs (results not shown). These two results suggestthat is more likely that the cytokine is produced in response to variations in the size ofthe mature cells compartment.

In order to proceed further, we have turned our attention to a model which correspondsto a more biologically-realistic situation: a stochastic model of a differentiation cascadeswhere stem cells can proliferate both symmetrically or asymmetrically. In our model, thechoice between these two division modes is a random event. Here, we have focused onthe robustness of the population generated by a mixed (i.e. with both symmetric andasymmetric stem cell division) differentiation cascade in terms of its resilience againstinvasion (for example, invasion by a potential tumorigenic differentiation cascade). Inparticular, we have focused on optimal strategies for a differentiation cascade to be robustagainst invasions as defined by the values of two parameters: the probability of symmetricbehaviour, e, and the rate of symmetric self-renewal, p. Our result show that thereexist a trade-off between robustness and the capability of the differentiation cascade foradaptive behaviour in response to some sort of stressful situation (e.g. healing responseto injury): Whereas the smallest the value of e, the more robust the population, we havealso found that upregulation of symmetric self-renewal, i.e. increased value of p, increasesthe likelihood of the population to withstand invasion attempts. However, we found that

Conclusions 81

populations with big values of p are not robust, in the sense that lead the population toextinction due fluctuations and the mechanisms studied in this section.

The issue of the optimal (i.e. maximally anti-tumour) architecture of has been recentlyaddressed Rodrıguez-Brenes et al. [107]. Comparisons between our results and theirs isdifficult as the models are based on slightly different premises. For example, the modelconsidered in [107] does not take regulation into account. In turn, our models do notaccount for self-renewal of the transient-amplifying cells. In spite of these differences,some of the results of Rodrıguez-Brenes et al. [107] are compatible with ours. For example,Rodrıguez-Brenes et al. [107] conclude that anti-tumour mechanisms should favour theevolution of shorter (few intermediate transient amplifying stages) differentiation cascadesover longer ones. The delay-induced extinction mechanism described in this section alsofavours the evolutionary prevalence of shorter cascades. The effect of self-renewal in theintermediate differentiation stages has not been explored in this section and it is left asthe subject for future investigation.

Related models has been proposed by Traulsen et al. [31, 78, 120, 125], where they usea Moran process to study the competition between two populations as a function of theduplication probability to find that the probability of invasion is an increasing function ofthe duplication probability of the invader population. We observe similar results, only forsmall values of the duplication rate p. For larger values of p the probability of invasiondecreases, as the invader population gets extinct due the extinction mechanisms thatwe have characterised in this section. This difference is due to the fact that we areconsidering a model with regulated-SC compartment size. Similarly, the models presentedin [70, 87, 39, 75, 88, 86, 127] also lack feed-back, and, therefore do no exhibit any ofphenomenology described in this work.

We also observe significant differences between deterministic and stochastic models.Whereas the mean-field behaviour model presents a critical point, characterised in termsof the net growth rate of the stem cell compartment, beyond which the resident populationstops being robust, the stochastic model exhibits non-monotonic behaviour in which prob-ability: The invasion probability is virtually zero below the deterministic critical point,however beyond the critical point it reaches a maximum and starts decreasing (due to theextinction mechanism studied in the symmetric division model).

5.2 Stochastic dynamics and extinction of HIV-1 in patientsunder potent anti-retroviral therapy

We have presented and analysed a stochastic model of HIV-1 dynamics in the blood streamunder anti-retroviral therapy. The analysis of the behaviour of this model has providedus with valuable insight on the ability of a therapy based on stimulating the activation oflatently infected cells through an study of the extinction of the infection.

In order to investigate the effect of the activation rate of the latently infected cells, aL,on the average extinction time of the latent reservoir we have used several methods. Wehave studied the mean-field behaviour of our system to determine the range of parametervalues in which a positive equilibrium (i.e. a steady-state associated to a persistent in-fection) exists. We have further performed Gillespie simulations of our stochastic modelto study the extinction average time of the system. This simulation results allow us todetermine if increasing the activation rate induces a reduction of the average extinctiontime. We have showed that, in the range of values of aL that we have analysed, the

82 Conclusions

extinction time is dramatically reduced from decades to 2-3 years. This result supportsthe viability of increasing the rate of activation of latently infected cells by a possibletherapeutic approach for total eradication of HIV-1 infection.

This problem presents a very particular multi-scale nature, each variable has a differenttime scale. This property has been used to reduce the dimensionality of the problem. Wehave performed a time-scale analysis of the Hamilton Equations to produce two differentQSS approximations. We further assume that the generalised coordinate correspondingto the healthy cells, q1, and its momentum, p1 are constant.

Separation of time-scales have also been used to justify the uses of a multi-scale stochas-tic algorithm [20], which is excellent for our problem. We have used the semi-classical ap-proximation [10, 73, 35] together with QSS approximations in the Hamilton equations [2]to approximate the generating function G(p, t) and we have observed excellent agreementbetween the results of the semi-classical approximation and the ones obtained via Gillespiesimulations, except in a small region near a∗L, the value of aL where the metastable stateceases to exist. In that region Gillespie simulations must be performed to compute theextinction probability. A better understanding of the Stokes phenomenon and analyticexpressions for all the terms in the semi-classical approximation are left for future work.

5.3 Effects of density fluctuations on the statistics of viralblips

The aim of Chapter 4 is to contribute to the discussion of whether viral blips in HIV-1 in-fected patients treated with HAART are symptomatic of ensuing drug failure or if, on thecontrary, they are random occurrences uncorrelated to drug-efficiency status and thereforelacking clinical significance. To this end, we have presented an inhom*ogeneous stochas-tic HIV-1 infection dynamics model, based on previous mean-field and stochastic models[110, 28], which accounts for local density fluctuations and burst virion production. Previ-ous models [110, 28] assume that the system is well-stirred and, therefore, hom*ogeneous.Given the small number of individuals involved in our model, in particular, regarding thenumber of virus copies during the latent infection phase, we have hypothesised that spatialinhom*ogeneities may play a role.

We have proposed an stochastic infection model which extends the model proposedin [28] by accounting for density inhom*ogeneities and bursty virion production. We haveavoided the consideration of factors such as random activation of latently infected cellsor upregulated virion production rate which have been introduced in previous models[110, 28]. We have further designed an in silico blood sample model, in order to reproduceas faithfully as possible the experimental protocol of Nettles et al. [95]. By doing so, wehave been able to show that our model closely reproduces the blip statistics obtained in[95], in particular when burst virion production is considered (see Figs. 4.2, 4.3 and 4.4),thus supporting the hypothesis that viral blips in HAART-treated patients are randomevents unrelated to drug failure [76, 94].

Previous modelling approaches have made different claims regarding the origin of blips.Rong & Perelson [110] have postulated that random activation of latently infected cellsby exposure to the corresponding antigens account for the presence of blips in HAART-treated patients. Conway & Coombs [28] have formulated a model in which a heightenedrate of virion production is necessary for their model to reproduce realistic values of blipfrequency, which could be interpreted as blips being the consequence of an evolutionary

Conclusions 83

process where new variants of the virus have emerged with an increased virion productionrate. Whilst our model does not disprove any of these mechanisms as well as othersproposed in the literature [37, 52, 41, 40, 43, 61, 62, 80], our model points out that viralblips may be purely random occurrences uncorrelated with factors of clinical significance.

Furthermore, our approach allows us to (partially) address the issue of whether lab-oratory and sample manipulation artifacts affect the observation of blips [76, 94]. Wehave investigated the effect that the post-extraction handling time, i.e. the time elapsedbetween sample extraction and the actual analysis, has on the statistics of the numberof blips. According to our model, this factor contaminates the statistics of the numberand duration of blips (see Fig. 4.5), which supports the position regarding the effectslaboratory artifacts on viral blip observation [76].

5.4 Future work

The models and results presented in this thesis open a number of interesting avenues forfuture research.

In Chapter 2 we have shown how the tissues are protected against mutations thatlead to symmetric division. However, according to our model, the optimal strategy isasymmetric division, since the number of stem cells remains constant and, therefore, theextinction probability is equal to zero, but this does not reflect the reality. The next stepis to formulate a new model including more reactions such as mutations in the stem cellcompartment (e.g. SC

µ−→ SCm, with SCm a mutant stem cell with another properties).This improvement would lead us to a more realistic situation in which the asymmetricdivision is not the best strategy.

In Chapter 3, using the semi-classical approximation we have obtained the mean extinc-tion time τ = AΩBeCΩ. Although in simple problems such as the branching-annihilation[35] or branching-annihilation-decay [68] A and B can be determined. In more generalsituations the determination of A and B requires a careful analysis of the associated Stokesphenomenon [10]. Such analysis, together with the development of asymptotic methodsto address this problem in more general settings is postponed for future work.

In Chapter 3 we have analysed the combination therapy of increasing the activation rateof the latently infected cells while continuously administering HAART, just by changingthe value of aL. A more realistic way to study this, is to formulate an accurate model ofhow the levels of histone deacetylases (HDAC) and other chromatin modifiers change theactivation of the latently infected cells [121].

In this thesis we have used the semi-classic approximation in the Chapter 3 to computethe average extinction time of a population. But this methodology can be applied tocompute the average extinction time in Chapter 2 or to count the number of blips inChapter 4.

In Chapter 4, we have modelled the dynamics of the HIV-1 infection in the bloodsamples, assuming that there is no replenishment of new T cells and the clearance rate isreduced. However, post-extraction dynamics require a more careful study. It is not wellknown which changes in the dynamics are produced, if the samples are refrigerated, thiscan slow the viral replication. We left for future work a complete study of the dynamicsin the blood samples.

84 Conclusions

Appendix A

Stochastic modelling and theGillespie algorithm

This appendix is devoted to give a brief introduction to stochastic processes, stochasticmodelling and a computational tool, called the Gillespie Stochastic Simulation Algorithm.

A.1 Basic definitions

Definition A.1. Let Ω be a set, let A be a collection of subsets of Ω. Then A is called aσ-algebra and the ordered pair (Ω, A) is called a measurable space if A satisfies:

1. Ω ∈ A

2. If B ∈ A, then, Bc := b : b ∈ Ω, b /∈ B ∈ A.

3. If An, n ∈ N ⊂ A, the union ∪∞n=1An ∈ A.

Definition A.2. Let (Ω, A) be a measurable space. Let P be a real valued set functiondefined on the σ-algebra A. The function P : A→ [0, 1] is called a probability measure ifit satisfies:

1. P (Ω) = 1.

2. (σ-additivity) If An, n ≥ 1 ⊂ A and Ai ∩Aj = ∅ for i 6= j, then

P (∪∞n=1An) =∑∞

n=1 P (An).

The ordered triple (Ω, A, P ) define a probability space.

Definition A.3. Let (Ω, A, P ) a probability space. Let A1, A2 ∈ A. The conditionalprobability of A1 given A2 is

P (A1|A2) = P (A1∩A2)P (A2) .

Theorem A.1.

P (A|B) = P (A∩B)P (B) = P (B|A)P (A)

P (B)

85

86 Stochastic modelling

Definition A.4. Let (Ω, A, P ) be a probability space. Two A1, A2 ∈ A are said to beindependent if

P (A1 ∩A2) = P (A1)P (A2).

The σ-algebra generated by the open sets of R is denoted as B and it is called Borel’sσ-algebra.

Definition A.5. A Random variable is an application X : Ω→ R such that

∀B ∈ B, X−1(B) ∈ A.

Definition A.6. Let X be a random variable taking values in R. The probability thatX ∈ (x, x + dx) is given by p(x)dx where p(x) is the Probability Distribution Function(PDF) of X.

Some properties of the PDF are:

•∫∞−∞ p(x)dx = 1

• P (x1 < X < x2) =∫ x2x1p(x)dx

Definition A.7. Let (Ω, A, P ) be a probability space. A stochastic process is a collectionof random variables Xt(s) : t ∈ T, s ∈ Ω where T is some index.

Such systems are described in in terms of an infinite set of joint probability densities:P (xn, tn, ...x1, t1) or, equivalently, by a set of joint conditional probability densities:

P (xn, tn, ...xk+1, tk+1|xk, tk, ...x1, t1). Obviously, such a system is not possible to deal within practice, therefore additional conditions must be imposed in order to make the systemtractable.

The Markov property is an additional condition whereby we assume that the systemloses memory of all its past history except for the most recent event.

Definition A.8. An stochastic process is Markovian, or has the Markov Property, if it issatisfied the following property:

P (xn, tn|xn−1, tn−1, ...x1, t1) = P (xn, tn|xn−1, tn−1).

Remark:This statement is equivalent to claiming that the waiting time between twosuccessive events is exponentially distributed [47].

This property allows us to write any joint PDF in the infinite hierarchy in terms ofjust one: the one-step PDF P (xn, tn|xn−1tn−1):

P (xntn, xn−1, tn−1, ...x2, t2|x1, t1) =∏ni=2 P (xi, ti|xi−1, ti−1)

Stochastic modelling 87

A.1.1 The Chapman-Kolmogorov Equation

The Chapman-Kolmogorov equation (CKE) is a direct consequence of the Markov prop-erty and provides a first step towards writing an equation for the time evolution of theprobability density. Consider the identity

P (x3, t3|x1, t1) =

∫P (x3, t3, x2, t2|x1, t1)dx2 (A.1)

Now

P (x3, t3, x2, t2|x1, t1) =P (x3, t3, x2, t2, x1, t1)

P (x1, t1)

=P (x3, t3, x2, t2, x1, t1)

P (x2, t2|x1, t1)

P (x2, t2, x1, t1)

P (x1, t1)

= P (x3, t3|x2, t2, x1, t1)P (x2, t2|x1, t1).

Applying the Markov property we get

P (x3, t3, x2, t2|x1, t1) = P (x3, t3|x2, t2)P (x2, t2|x1, t1).

This lead us to the Chapman-Kolmogorov Equation (CKE)

P (x3t3|x1t1) =

∫P (x3, t3|x2, t2)P (x2, t2|x1, t1)dx2 (A.2)

A.1.2 The Master equation

The Master equation is a reformulation of the CKE that is easier to handle and moredirectly related to physical models. Consider P (x3, t3|x2t2) and let dt = t3 − t2. Then

P (x3, t3|x2, t2) = (1− a0(x2)dt)(x3 − x2) +W (x3|x2)dt. (A.3)

W (x3|x2) is the probability per unit time (transition rate) of the transition x2 → x3, and1− a0dt is the probability of no transition occurring. Therefore

a0(x2) =

∫W (x3|x2)dx3. (A.4)

According to the Chapman Kolmogorov Equation

P (x3, t2 + dt|x1t1)− P (x2, t2|x1, t1dt

=

∫W (x3|x2)P (x2, t2|x1, t1)dx2−a0(x3)P (x3t2|x1, t1)

(A.5)Using the definition of a0(x3) and taking the limit dt→ 0, we obtain

dP (x3, t|x1, t1)

dt=

∫(W (x3|x2)P (x2, t|x1, t1)dx2 −W (x2|x3)P (x3t|x1, t1))dx2. (A.6)

If only a discrete set of transitions Wi(X), is possible, and, upon occurrence of reactioni, X → X + ri, then, the corresponding Master equation reads:

dP (X, t)

dt=∑i

(Wi(X − ri)P (X − ri)−Wi(X)P (X, t)). (A.7)

88 Stochastic modelling

A.2 classic examples

A.2.1 The Moran process

The Moran process, named after the australian statistician Pat Moran, is a widely-usedvariant of the Wright-Fisher model and is commonly used in population genetics.

• N individuals of two types. N is keep fixed.

• n: number of normal individuals. m: number of mutant individuals. N=m+n.

• At each time step:

– n→ n+ 1 and m→ m− 1 with probability rate W+(n) = nN

(1− n

N

).

– n→ n− 1 and m→ m+ 1 with probability rate W−(n) = nN

(1− n

N

).

Note that W (n + 1) = W (n − 1), i.e. E[∆n] = E[n(t + ∆t) − n(t)] = 0. This impliesthat with m = E[n]:

dm

dt= 0 (A.8)

The system has two absorbing states: W+(n = 0) = W−(n = 0) = 0 and W+(n =N) = W−(n = N) = 0. This means that

limt→∞

P (n(t) = 0 ∪ n(t) = N) = 1 (A.9)

However, this behaviour is not at all captured by the deterministic equation, whichpredicts that the population will stay constant, as can be seen in Fig. A.1

A.2.2 The stochastic Logistic Equation

As a second example we choose the logistic equation A.10,

dm

dt= m

(1− m

K

), (A.10)

has two steady states: m = 0 unstable and m = K stable, i.e. regardless of thevalue of K and for any initial condition such that m(t = 0) > 0, m(t) will asymptoticallyapproach K.

Consider now a continuous-time Markov process nt whose dynamics are given by thefollowing transition rate:

Stochastic modelling 89

0 10 20 30 40 50 60 70 80 90 1000

5

10

15

20

25

30

35

40

45

50

Time

n

Figure A.1: Moran realizations

• n→ n+ 1 with probability rate n.

• n→ n− 2 with probability rate n(n−1)2

1K .

This stochastic process has a unique absorbing state: n = 0. This, together with thepopulation can not grow to infinity (a linear term increasing population against quadraticterm decreasing the population) leads to

limt→∞

P (n = 0, t) = 1. (A.11)

We see in Fig A.2 that the stochastic dynamics show strong discrepancies with Equa-tion A.10 when randomness is dominant. We observe that for small K fluctuations domi-nate the behaviour of the system. Extinctions are common for small K , in contradictionto the behaviour predicted by the logistic equation Eq A.10, and become rarer as K isallowed to increase.

A.2.3 Equilibrium and absorbing States

There are several definitions of equilibrium states in stochastic systems, we will explainour with an easy example. Consider again the stochastic logistic growth, i.e. a process ntsuch that:

• n→ n+ 1 with probability rate W+(n) = n

• n→ n− 2 with probability rate W−(n) = n(n−1)2

1K

90 Stochastic modelling

0 10 20 30 40 50 60 70 80 90 1000

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time

n/K

Figure A.2: Red line K = 10, green K = 50, blue K = 100, black K = 1000

Consider a state of system, ns, is, roughly speaking, a state of the process such thatW+(ns) = W−(ns). W+(ns) is the number of births within a population of ns individuals,likewise, W−(ns) = the number of deaths within a population ns individuals. So an steadystate of our population dynamics is reached when nt = ns, since death rate is balanced bybirth rate and therefore the population stays roughly constant. It is easy to derive thatns = K, which coincides with the deterministic stable fixed point. The idea is to notethat W+(n)−W−(n) > 0 if n < ns and W+(n)−W−(n) < 0 if n > ns.

Another important concept is the absorbing state. An absorbing state, n0, is char-acterised by Wi(n0) = 0 i.e. once the system has reached the absorbing state, it cannotleave anymore.

Consider again, the stochastic logistic growth rate, we have:

• Steady states are in general not absorbing states.

• W+(ns) 6= 0 and W+(ns) 6= 0.

• If n = 0 then W+(0) = W−(0) = 0 therefore n = 0 is an absorbing state.

ns belongs to the set of accessible states of n = 0, that means there is at least oneconsecutive set of transition that connects ns and n0. For Example:

K → K − 1→ K − 2→ · · · → 1→ 0.

However, if K 1 the probability of such a chain of events is vanishingly small.ns is an steady state in the sense that births and deaths are balanced. Moreover,

W+(n) −W−(n) > 0 if n < ns and W+(n) −W−(n) < 0 if n > ns. This is essentiallyequivalent to what happens in the deterministic logistic growth model. However, ns isnot an absorbing state of the stochastic dynamics. The only absorbing state is n =0. Stochastic extinctions are relatively rare provided K is big, if this is the case, then,the deterministic system provides a reasonable approximation to the behaviour of the

Stochastic modelling 91

model. If, on the contrary, K is small stochastic extinctions are relatively common andthe deterministic description is not an accurate one. In general, we should expect non-trivial random effects for small populations or dynamics with absorbing states.

A.3 Monte Carlo Methods

There is no consensus on how Monte Carlo method should be defined, we use the followingdefinition:

Definition A.9. a Monte Carlo Method is an algorithm which uses for any reason apseudorandom number.

Example A.1. Approximate the value of π

Consider a random point P in the unit square, which is the probability that P is inthe unit circle?

Clearly, the probability has to be the quotient of the areas, that isπ

4. Using this we

will use the following algorithm to compute the value of π:

1. Set A = 0.

2. Generate two random numbers x, y ∈ [−1, 1] if x2 + y2 < 1, then A = A+ 1.

3. Repeat (2) N times. When N →∞, then AN →

π4 , this lead us to π ≈ 4AN .

A.4 Gillespie Stochastic Simulation Algorithm

Consider a Chemical reacting system, that is, molecules of N chemical species S1, ..., SN ,inside some volume Ω, at some temperature T . Interacting through M elemental reactionchannels W1, ...,WM , where Wj is assumed to describe a single instantaneous physicalevent which changes the population of at least one species.

Each reaction channel Wj can then be characterized by two entities [47]:

• Its propensity function aj(x): If the system is currently in state x,

aj(x) · dt := probability that one Wj event will occur in the next dt.

• Its state change vector vj ≡ (v1j , ..., vNj):

vi,j = the change in Xi caused by one Wj event.

Wj induces the state change x→ x+ vj .

Example A.2. Combustion of methane

Consider 4 different chemical species: Methane (CH4), oxygen (O2), carbon dioxide(CO2) and water (H2O). And the reaction

CH4 + 2O2 → CO2 + 2H2O

92 Stochastic modelling

The state change vector of this reaction is v1 = (−1,−2,+1,+2), and its propensity

function a1(CH4, O2, CO2, H2O) = r1CH4O2(O2−1)

2 .The Master equation describing this process has the form:

P (n1, n2, n3, n4, t) = r1(n1 + 1)(n2 + 2)(n2 − 1)

2P (n1 − 1, n2 + 2, n3 − 1, n4 − 2)

−r1n1n2(n2 − 1)

2P (n1, n2, n3, n4)

Note that, both, vj do not determine the reaction, as con be seen in the followingexample

Example A.3. S1 → S2: vj = (−1, 1, 0, ..., 0); aj(x)dt = (cjdt)x1 =⇒ aj(x) = cjx1

S1 + S2 → 2S2: same vj ; aj(x)dt = (cjdt)x1x2 =⇒ aj(x) = cjx1x2

The main idea of the Gillespie Stochastic Simulation Algorithm is to produce samplepaths or realizations of X(t). To do this, the method generates properly distributedrandom numbers for the time τ to the next reaction, and the index j of that reaction.

Let p(τ, j|x, t) · dτ ≡ be the probability that, given X(t) = x, the next reaction willoccur in [t+ τ, t+ τ + dτ), and will be Wj .Let

a0(x) ≡M∑k=1

ak(x).

Then

p(τ, j|x, t)dτ =(

1− a0(x)τ

n

)naj(x)dτ → e−a0(x)τaj(x)dτ

Therefore,

p(τ, j|x, t)dτ = e−a0(x)τaj(x)dτ = a0(x)e−a0(x)τ · aj(x)

a0(x)dτ .

Note that τ is the exponential r.v. with mean 1a0(x) , j is the integer r.v. with probability

massaj(x)a0(x) and j and τ are statistically independent.

The following scheme is a Method of implementing the SSA:

1. In state x at time t, evaluate a1(x), ..., aM (x) and a0(x) =

M∑k=1

ak(x).

2. generate r1 and r2 random numbers in [0, 1], and compute τ and j according to

τ =1

a0(x)ln

(1

1− r1

),

j = the smallest integer satisfying

j∑k=1

ak(x) > r2a0(x).

3. t = t+ τ and x = x+ vj

4. Record (x, t). Go to 1 or end the simulation.

Appendix B

Semi-Classical Approach

This appendix is devoted to explain the semi-classical approximation used in Chapter 3to approximate the solution of a Master equation.

The basic idea of the semi-classical approximation, is to transform a stochastic processdescribed by a Master Equation.

Consider a stochastic process described by the following Master equation

∂P (n, t)

∂t=∑j

(Wj(n− rj , t)P (n− rj , t))−Wj(n, t)P (X, t)) . (B.1)

Let us define the generating function as Equation

G(p, t) =∞∑n=0

pnP (n, t) (B.2)

Clearly, if G is known, the probability distribution function Pn, can be recovered [60]

P (n, t) =1

n!

∂nG

∂pn

∣∣∣∣p=0

. (B.3)

G satisfies linear partial differential equation

∂G

∂p(p, t) = H

(p,

∂

∂p

)G (B.4)

satisfying G(1, t) = 1. This equality means the sum of the probabilities is equal to 1. Theorder of this partial linear differential equation is n if and only if, at most, reactions with nspecies reacting in the left hand side of the reaction appears, that isA1+...+An → anything[35]. If the order is bigger than one, in general we can not solve this PDE. However, wecan consider this as a Schrodinger Equation,

∂

∂tG = −HG. (B.5)

Equation (B.5) tell us that there exists an S such that G = exp (−S). When the systemsize, Ω, is big, one may employ the WKB approximation [73] and consider the re-scaledaction S with S = ΩS, expanding up to order O(Ω−1) (neglecting O(Ω−2) terms), oneobtains the classical Hamilton-Jacobi equations [6]

∂S

∂t= H

(p,∂S

∂p

). (B.6)

93

94 Semi-Classical Approach

Example B.1. Consider a stochastic process the binary annihilation reaction A+Ak−→ ∅.

This process is described by the Master equation

∂

∂tP (n, t) =

(n+ 2)(n+ 1)

2kP (n+ 2, t)− (n)(n− 1)

2kP (n, t), (B.7)

consider the generating function

G(p, t) =∑n≥0

pnP (n, t), (B.8)

then, Eq (B.7) can be written as

∂G

∂t=

∑n≥0

pn(n+ 2)(n+ 1)

2kP (n+ 2, t)−

∑n≥0

pnn(n− 1)

2kP (n, t)

=∑n≥2

pn−2n(n− 1)

2kP (n, t)− p2

∑n≥0

pn−2n(n− 1)

2kP (n, t)

=∑n≥0

pn−2n(n− 1)

2kP (n, t)− p2

∑n≥0

pn−2n(n− 1)

2kP (n, t)

=k

2(1− p2)

∂2G

∂p2. (B.9)

We consider this equation as a Schrodinger equation,

∂

∂tG = −HG, (B.10)

with

H(p, q) =k

2(p2 − 1)q2, (B.11)

where q = − ∂

∂p.

Using the Ansatz G(p, t) = exp (−ΩS(p, t)), Eq. (B.10) reads

−Ω∂S

∂texp (−ΩS(p, t)) = −k

2(p2 − 1)

∂2

∂p2exp (−ΩS(p, t))

= −k2

(p2 − 1)

(−Ω

∂

∂p

∂S

∂pexp (−ΩS(p, t)) + Ω2

(∂S

∂p

)2

exp (−ΩS(p, t))

),

expanding up to order Ω−1, we obtain the Hamilton-Jacobi equation

∂S

∂t= H

(p,∂S

∂p

)=k

2(p2 − 1)

(∂S

∂p

)2

, (B.12)

with k = Ω−1k.

Semi-Classical Approach 95

Instead of dealing directly with the Hamilton-Jacobi equations, which, in general, wedo not know how to solve, we exploit the analogy with the Schrodinger equation to use aFeynman path-integral representation further to obtain the solution of Eq. (B.5)

G(p, t) =

∫ t

0(exp (−S(p, q))Dq(s)Dp(s)) ds+ px0 , (B.13)

with

S(p, q, t) =

∫ t

(−H(p, q) +

n∑i=1

qi(s)pi(s)

)ds. (B.14)

Too see this note that, in the Fourier representation [22],

δ(p− p0) =

∫ ∞−∞

dπe−i(p−p0)πdπ. (B.15)

The probability of pass from a p′′ to a p′ can be seen as

(p′| exp

(−Ω∆tH(x′,Ω−1 ∂

∂p′, t)

)|p′′) = (B.16)

=

∫ ∞−∞

dk exp(−(Ω∆t)H(p′, iεk, t) + ik(p′ − p′′) +O(∆t2)

)=

∫ ∞−∞

dπΩ exp

(−(Ω∆t)H(p′, iπ, t)− iπ p

′ − p′′

∆t

)+O(∆t2).

for a short time ∆t. Defining exp← as an ordered exponential, the solution of Eq. (B.5)can be written as

G(pn, tn|p0, t0) = exp←

[−Ω

∫ t

t0

dsH

(p,Ω−1 ∂

∂p, s

)]δ(p− p0)

= lim∆t→0

n∏j=1

exp

(−Ω∆tH(p,Ω−1 ∂

∂pj, t)

)δ(pj − pj−1)dpj−1.(B.17)

Writing Eq. (B.17) in the representation of Eq. (B.16) we obtain

G(p, t|p0, t0) = limn→∞

∫· · ·∫dp1dπ1 · · · pndπnΩn

= exp

Ωn∑j=1

∆t

[H(pj , iπj , tj)− iπ

pj − pj−1

∆t

]. (B.18)

In the limit ∆t→ 0 we recover Eq. (B.13). Note that, in fact, the integral of the Equation(B.13) represents all possible trajectories, we just have composed the solution at eachinfinitesimal time step.

The semiclassic approximation consists on taking just the most probable path, that is,to approximate the path integral (B.13) by

G(p, t) ≈ exp (−S(p, t)), (B.19)

where pi are the solution of the Hamilton equations of the (classic) Hamiltonian H, sincethese equations are the ones that maximizes exp (−S(p, q)).

96 Semi-Classical Approach

The approximation of this integral B.13 is performed via a Laplace method or a steepestdescent method, which is valid if Ω is large enough. The Laplace method consists in thefollowing, if we have an integral of the form∫ b

aexp (Ωf(x))dx,

with Ω big enough and f(x) has a unique maximum at x0 ∈ [a, b], then we can approximate∫ b

aexp (Ωf(x))dx ≈ eΩf(x0)

∫ b

ae−Ω|f ′′(x)|(x−x0)2 . (B.20)

For the values of x far from x0 this decays exponentially fast.And in Equation (B.13) the maximum is reached when S is the solution of the Hamilton-

Jacobi equations. The path that maximizes the action integral, IS :

IS =

∫S(p, q, s)ds, (B.21)

is determined by the Euler Equation

d

ds

∂S

∂q− ∂S

∂p= 0, (B.22)

or by the canonical equations of motion

p(s) = −∂H∂q

,

q(s) =∂H

∂p,

with boundary conditions

qi(0) = xi(0),

pi(t) = pi.

Remark: We have used the Generating function formulation for convenience, all thisderivation can be done using directly the master equation, see for instance [73]

B.1 Metastablility

Equation (B.5) can, in general, be solved, only in few cases (e.g. if H is of order 1. In thatcase, one can use the characteristic method to solve that partial differential equation). Ingeneral we cannot solve the Hamilton-Jacobi PDE. However, taking advantage of the factthat H is linear to expand G in their eigenvalues λk and eigenmodes φk, we obtain

G(p, t) = Gst(p)−∑k≥1

φ(p)e−λkt (B.23)

where Gst(p) is the stationary solution (corresponding to λk = 0), λk are the eigenvaluesof −H, φk(p) the eigenmodes.

Semi-Classical Approach 97

Assume the system has a metastable state, which is an attractor of the mean-fielddynamics. In the semi-classical approach we can think it as a positive fixed point onpi = 1, or as an hyperbolic point of the Hamiltonian with stable manifold pi = 1. Wecan expect the dynamics evolve according the mean-field behaviour until they reach themetastable state. However, we have to note that, although the metastable state is anstable fixed point of the mean-field behaviour, it is not thermodynamically stable. If itwere, we would have reactions, that require energy, at infinite time, in an isolated system,and this contradicts the second-law of thermodynamics [38]. Thus, after a long timenear the metastable the system will escape to the stochastic extinction or to a differentmetastable state in the case of multi-stable systems. If the system size is large enough, wecan expect that the dynamics will follow a trajectory close to the heteroclinic connectionas this trajectory is the optimal one. The probability of observing a different trajectoryis exponentially small [35].

In this case, as we are interested in large values of t, we can to approximate theextinction probability by truncating Eq. B.23 at k = 1, since the exponential terms withnegative eigenvalues (−λk < −λ1 < 0) decay very fast and we obtain

G(0, t) ≈ 1− exp (−λ1t). (B.24)

with

τ =1

λ1(B.25)

as the mean extinction time.In Chapter 3 this mean extinction time is of the form τ = AΩB exp(CΩ), where C is

the integral of the re-scale action, S, along the heteroclinic connection. In the particularcase of the Lower bound 3.31, this can be solved analytically [10, 68], and this example isvery illustrative.

Example B.2. Consider the stochastic process described by 3.31, that is:A

r−→ A+A.

AaL+dL−−−−→ ∅.

A+Ar

Lmax−−−−→ ∅.The Master equation reads

dP (n, t)

dt= r(n− 1)P (n− 1, t)− rnP (n, t) + (aL + dL)(n+ 1)P (n+ 1, t)− (aL + dL)nP (n, t)

+r

LmaxΩ

(n+ 2)(n+ 1)

2P (n+ 2, t)− r

LmaxΩ

(n)(n− 1)

2P (n, t), (B.26)

and can be rewritten as

∂G

∂t= r(p− p2)

∂G

∂p+ (aL + dL)(p− 1)

∂G

∂p+

r

2LmaxΩ(p2 − 1)

∂2G

∂p2. (B.27)

As we have explained, we consider this as Schrodinger equation,

∂G

∂t= −HG, (B.28)

whereH(p, q) = r(p2 − p)q + (aL + dL)(1− p)q +

r

2LmaxΩ(1− p2)q2. (B.29)

98 Semi-Classical Approach

Making an eigenmode expansion of G one gets

G(p, t) = Gst(p)−∑k≥1

φ(p)e−λkt (B.30)

and, for large values of t, G(p, t) can be approximated by

G(p, t) ≈ 1− φ(p)e−λ1t. (B.31)

Note that if G(p, t) = 1− φ(p)e−λ1t, then

∂G

∂p= −φ′(p)e−λ1t,

∂2G

∂p2= −φ′′(p)e−λ1t,

∂G

∂t= φ(p)e−λ1tλ1.

This allows us to reduce our PDE to into an Sturm-Liouville problem of the form

− Hφ+ λ1φ = 0 (B.32)

To solve this problem we proceed as in [10] and we look for a solution in the momentumspace. However, it was previously solved via a WKB approximation by Kessler and Shnerb[68]. Eq. (B.32) is

r

2LmaxΩ(1− p2)φ′′ + (p− 1)(rp− (aL + dL))φ′ + λ1φ = 0, (B.33)

We proceed as follows, we divide the problem in three regions: the first region isA := p | p ≈ 0, B := p | p ≈ 0 and C := [0, 1]−A−B. We look for a solution nearp ≈ 0, where from Eq. (B.31) we have φ(p) ≈ 1, then a solution p ≈ 1, where, from Eq.(B.31) we have φ(p) ≈ 0, and finally we match both solutions in the three regions.

We use a WKB Ansatz φ(p) = a(p)exp(−Ωs(p)), then

φ(p) = a(p)e−ΩS(p)

φ′(p) = a′(p)e−ΩS(p) + ae−ΩS(p)(−ΩS′)

φ′′(p) = a′′(p)e−ΩS(p) + 2a′(p)e−ΩS(p) + a(p)e−ΩS(p)(−ΩS′)2.

so, our Sturm-Liouvulle problem reads,

r

2LmaxΩ

(a′′(p)e−ΩS(p) + 2a′(p)e−ΩS(p) + a(p)e−ΩS(p)(−ΩS′)2

)+(p− 1)(rp− (aL + dL)

(a′(p)e−ΩS(p) + ae−ΩS(p)(−ΩS′)

)= 0. (B.34)

The terms in Ω, gives us a differential equation for S(p)

S′(p) = 2Lmaxrp− (aL + dL)

r(p− 1)(B.35)

note that, the right-hand term is exactly the heteroclinic connection between the metastablestate and the stochastic extinction.

Semi-Classical Approach 99

The terms without Ω gives as a differential equation for a(p)

a′(p) = Lmax(p− 1)[rp− (aL + dL)]

r(1− p2)S(p))(B.36)

To study the non-quasi-stationary region we add the term λ1φ to the equation pertur-batively, φ(p) = 1 + δφ(p) with δφ(p) 1, and we get

r

2LmaxΩ(1− p2)φ′′ + (p− 1)(rp− (aL + dL))φ′ = −λ1, (B.37)

which gives us a differential equation for φ′(p)

φ′′(p) =2Lmaxr(1− p2)

Ω− λ1 − (p− 1)(rp− (aL + dL)φ′, (B.38)

in this case, the solution has been derived in [10] and is of the form

φ′(p) = −2Ωλ1e2Ω

(p−

(1+

aL+dLr

)log(1+p)

) ∫ p

−1

exp(

2Ω(s−

(1 + aL+dL

r

)log(1 + s)

))1− s2

ds

(B.39)Eq. (B.39) can be evaluated by the saddle point approximation with a saddle point ataL + dL

r

φ′(p) ≈ −2λ1

√LmaxπΩ

(r

aL+dL

)3/2

√r

aL+dL+ 1

(r

aL+dL− 1)e−2ΩLmaxI(p), (B.40)

with

I(p) =aL + dL

r−(

1 +aL + dL

r

)log

(1 +

aL + dLr

)−(p−

(1 +

aL + dLr

)log(1 + p)

).

(B.41)

The quasi-stationary solution is valid as long as pr

aL + dL− 1 1√

Ω, and the pertur-

bative solution holds when 1− p 1√Ω

[9, 8], so both solutions are valid when

1√Ω p

r

aL + dL− 1 1.

Finally, we find the eigenvalue λ1 matching both solutions.

λ1 =

√Lmax

4π

(r

aL + dL+ 1

)( raL+dL

− 1)2

( raL+dL

)5/2

√Ωe−ΩS0 , (B.42)

with

S0 = Lmaxr aL+dL

r − log(aL+dL

r + 1)

(r + aL + dL)

r− Lmax

r − log(2)(r + aL + dL)

r.

(B.43)

100 Semi-Classical Approach

Appendix C

Computation of HeteroclinicConnections for HamiltonianSystem

This chapter is devoted to give a brief introduction to the methods used in Chapter 3 forthe computation of the invariant manifolds

C.1 The Taylor Method

Consider the initial value problem

x(t) = f(x, t),

x(t0) = x(0) (C.1)

One step of the taylor method of order m, to compute the solution x(t) of the differ-ential equation (C.1) from time t to t+ h, of order m with step size h is

xN = xN−1 + x′(tN−1)h+ · · ·+ x(p)(tN−1)

m!hm,

starting with x0 = x(0).

This high order method allow us to obtain high accuracy and large step-sizes h alongthe integration (in the implementation m is choosen large and h is choosen to keep thelocal error of the method beyond the precision required, for example m = 24 is a commonvalue if double precision is required). Although this method is very simple, computingm derivatives of a vectorial field can be very hard. This problem is solve by using theso-called Automatic Differentiation Formulas.

C.2 Automatic Differentiation Formulas

Automatic Differentiation is a recursive algorithm to evaluate the derivatives of a closedexpression on a given point [64, 51].

Definition C.1. Let f be a smooth function. The normalized j−th derivative of f at t is

101

102 Computation of heteroclinic connections

f [j](t) =1

j!f (j)(t).

Some of the differentiation formulas are as follows:

1. If f(t) = g(t)± h(t), then f [n](t) = g[n](t)± h[n](t).

2. If f(t) = g(t)h(t), then f [n](t) =∑n

i=0 g[n−i](t)h[i](t).

3. If f(t) =g(t)

h(t), then f [n](t) =

1

h[0](t)

(g[n](t)−

∑ni=1 h

[i](t)f [n−i](t)).

There are automatic differentiation formulas for exponentials, logarithms, trigonomet-ric functions, and so on.

Example C.1. In Chapter 3, the Weak QSS approximation produces a reduced systemof 4 ordinary differential equations

dp2

dt=

1

ε1

(−κ3(p2

2 − p2)− κ3L0

Lmax(1− p2

2)q2 − κ4(1− p2)− κ5(p3 − p2)

),

dq2

dt=

1

ε1

(η(1− ε)q1q4 + κ3q2(2p2

2 − 1)− κ3L0

Lmaxp2q

22 − κ4q2 − κ5q2

),

dp3

dt=

1

ε2(−κ6(p3p4 − p3)− κ7(1− p3)) ,

dq3

dt=

1

ε2((1− η)(1− ε)q1q4 + κ5q2 + κ6(p4 − 1)q3 − κ7q3) , (C.2)

where

p1 = 1,

q1 = 1,

p4 =−(1− η)(1− ε)q1p3 − η(1− ε)p2q1 − κ8 − εp1q1

−(1− η)(1− ε)q1p1 − η(1− ε)p1q1 − κ8 − εp1q1

q4 =κ6

p1q1 + κ8p3q3,

for simplicity p4 = Ap3 +Bp2 + C and q4 = Dp3q3.In this problem, computation of the derivatives lead to very complicated expressions.

For that reason, we define the following auxiliar variablese1 = Dp3q3,e2 = Ap3 +Bp2 + C,e3 = p2p2,e4 = e3 − p2,e5 = 1− e3,e6 = e5q2,e7 = 1− p2,e8 = p3 − p2,

e9 =1

ε1

(−κ3e4 − κ3

L0

Lmaxe5 − κ4e7 − κ5e8

),

e10 = (e3 − 1),

Computation of heteroclinic connections 103

e11 = e10q3,

e12 =1

ε1((1− η)(1− ε)q1e1 + κ5q2 + κ6e11 − κ7q3),

e13 = p3e2,

e14 = e13 − p3,

e15 = 1− p3,

e16 =1

ε2(−κ6e14 − κ7e15),

e17 = e2 − 1,

e18 = e17q3,

e19 =1

ε2((1− η)(1− ε)q1e1 + κ5e2 + κ6e18 − κ7q3).

Using the Automatic Differentiation formulas, for e[n]i , we obtain

p[n+1]2 =

1

n+ 1e

[n]9 ,

q[n+1]2 =

1

n+ 1e

[n]12 ,

p[n+1]3 =

1

n+ 1e

[n]16 ,

q[n+1]3 =

1

n+ 1e

[n]19 .

In this particular case we only need the formulas for the sum and the product of twovariables.

C.3 Heteroclinic connections

In Chapter 3, in Fig. 3.4 we compared two different approximations of the exponentialgrowth of the mean relaxation time (3.30) for the semi-classical approximation. The esti-mates of ΩC are obtained after performing either a one-step (weak) or a two-step (strong)quasi-stationary approximation of the classical equations of motion and computing theintegral of the classical (Lagrangian) action S along the only non-trivial path connectingthe metastable state and the stochastic extinction of the system. Looking at the phasespace of the mean field system, one realizes that this non-trivial path is given by theheteroclinic trajectory γ(t) from the hyperbolic fixed point p1 (with pi = 1 and xi > 0)and the fixed point x0 located on xi = 0 with 0 < pi < 1 (here i = 2 for the strong QSSapproximation and i = 1, 2 for the weak QSS case).

It is a simple exercise to obtain an expression of the heteroclinic connection in theStrong QSS approximation, since, under this approximation, the system is reduced to onedegree of freedom, and, therefore, the heteroclinic trajectory is simply determined by con-servation of energy. However, the computation of γ(t) requires numerical techniques forthe Weak QSSA case since one has to deal with a 2 degrees of freedom (2-dof) Hamilto-nian system. Moreover, the Weak QSSA reduction exhibits an slow-fast structure, henceone has to deal with two quite different time-scalings when performing the numerics. Weshould note that all the computations have been performed using multiprecision arith-metics (around 100 digits where enough for most of the aL values considered) and thecodes have been implemented in PARI/GP [11].

Let us consider, from now on, the 2-dof Hamiltonian case whose related equations ofmotion, given by (3.28), are expressed in the coordinates x2, x3, p2, p3. The hyperbolic-hyperbolic fixed point p1 has a 2-dimensional stable and a 2-dimensional unstable invariant

104 Computation of heteroclinic connections

manifolds, to be denoted W s(p1) and W u(p1) respectively. Similarly, W s/u(x0) willdenote the 2-dimensional stable/unstable manifold associated to the fixed hyperbolic-hyperbolic point x0. The heteroclinic connection γ(t) corresponds to the intersectionW u(p1) ∩W s(x0).

To compute W u(p1) (and also W s(x0)) we use the so-called parametrization method[17, 18]. Basically the idea is the following. We represent the local invariant manifoldaround p1 as a vector series G(s1, s2), being G : R2 → R4. That is, a point X =(x2, x3, p2, p3) ∈ R4 will be considered to be on W u(p1) if X = G(s1, s2) =

∑i,j≥0 ai,js

i1sj2.

The coefficients ai,j ∈ R can be order by order computed by imposing the so-called invari-ance condition. This condition requires that the dynamics within the invariant manifold,expressed in the s1, s2 coordinates, must be conjugated to the linear dynamics aroundthe hyperbolic-hyperbolic point. Denote by λ1, λ2 > 0 (resp. < 0) the real eigenval-ues associated to p1 (resp. to x0, below we denote by Gp1 and Gx0 the correspondingparametrizations). Then, the linear dynamics is s1 = λ1s1, s2 = λ2s2. If X = f(X) refersto the equations 3.28, then the invariance condition requires

∂G

∂sν(s1, s2)λνsν = f(G(s1, s2)), ν = 1, . . . , 2.

By imposing this equality order by order, being k = i + j ≥ 1 the total order, one getsa sequence of linear systems for the coefficients ai,j , k > 1 (for k = 0 one gets the fixedpoint condition, and for k = 1 the eigenvectors system). See [115] for further details onthis procedure.

In practice, we truncate the series representation to suitable order N and we requirethat the invariance condition holds up to a tolerance tol. Then, the invariance conditionwill hold for X in a domain of radius rp1 around the fixed point p1 (similarly, for x0

we obtain rx0). Typical values used in the computations are N = 150 and tol = 10−40.The local domain size depends on the parameter aL in the system. As an example, foraL = 0.15 one obtains rp1 = 265 and rx0 = 1470.

Once the local representations Gp1 and Gx0 , of W u(p1) and W s(x0) respectively, areobtained we extend the W u(p1) up to Σ = p2 = Gx03 (rx0 , 0), where Gx03 denotes the 3rdcomponent of Gx0 (the one corresponding to the p2-variable). This is done by integratingthe equations (3.28) using a Taylor method, which turns out to be an appropiate time-stepper for the high precision computations required. Assume that Gp1(z1, z2) ∈ R4 is apoint such that the transported one Gp1(z1, z2) ∈ Σ is close to a zero of F (s1, s2, z1, z2) =Gx0(s1, s2) − Gp1(z1, z2), F : R4 → R4. Then we refine this initial condition, using aNewton method and we obtain a point of the heteroclinic orbit as a zero of F . Note thatthis requires to integrate the first variational equations. See [115] for further details.

Similar computations for a 2-dof system as well as for a 4D map were carried out in[45], where further technical details on the parametrization method and the computation ofhom*o/heteroclinic trajectories can be found. For a recent overview of the parametrizationmethod for computation of invariant manifolds in different contexts see [53].

Appendix D

Multi-scale Stochastic Simulations

The Gillespie Stochastic Simulation Algorithm (SSA) is exact. Its procedurally simple,even when the CME is intractable. There exist several implementations of the SSA, suchas, the direct method [47], the first reaction method [47], the next reaction method [46],thefirst-family method [49], the modified direct method [21], the sorting direct method [84],the modified next reaction method [5], the composition-rejection SSA method [116] and thepartial propensity direct method [103] among others. However, the stochastic simulationalgorithm remains too slow for most practical problems: simulating every reaction eventone at a time usually takes too much time if any reactants are present in very largenumbers. For example, in the problem presented in Chapter 4 and Chapter 3, the healthycells are given in a extremely big number, O(109) cells, if we consider them in our stochasticmodel, the Gillespie algorithm takes very small time steps. In this appendix we describea different method to avoid this problem.

D.1 Cao-Gillespie-Petzold Method

Cao et al. [20] presented a multi-scale stochastic simulation algorithm, based on a stochas-tic partial equilibrium assumption, which is much more efficient than the Gillespie stochas-tic simulation algorithm, particularly in the presence of very small populations or fastreactions.

D.1.1 Definitions

Assume we have N different species S1, ..., SN , Xi(t) denotes the number of Si particles attime t. Assume there are M different reactions W1, ...,WM . Each reaction channel, Wi,changes state vector of the system by ri. The process is described by the following masterequation:

dP (X, t)

dt=∑i

(Wi(X − ri)P (X − ri)−Wi(X)P (X, t)). (D.1)

Let W fj = (W 1

j , ...,Wmf

j ) be the set of fast reactions and W sj = (W 1

j , ...,Wmsj ), the set

of slow reactions.

Similarly, we define Xs = (xs1, ..., xsns

) as the set of the slow species and Xf =

(xf1 , ..., xfnf ), the set of the fast species. A species is considered to be fast if it is altered

by, at least, one of the fast reactions of the system. Otherwise it is a slow species.

105

106 Multi-scale Stochastic Simulations

We define two different propensity functions: The slow propensity functions

asi (x) = asi (xf , xs), i = 1, ...,ms

and the fast propensity functions,

afi (x) = afi (xf , xs) i = 1, ...mf .

We also define the fast and slow vector change as,

vfi =(

(vff1 , ..., vffnf, 0, ..., 0)

),

vsi =(vfs1 , ..., vfsnf

, vss1 , ..., vssns

)),

where the superindex ff denotes the fast reaction affecting a fast species, fs the slowreaction affecting a fast species, and ss the slow reactions affecting a slow species.

In Appendix A, we have defined an equilibrium state in stochastic systems, as a statein which, roughly speaking, the births and deaths are balanced. Consider now just thereduced stochastic system described by the fast species with the fast variables. If thereduced stochastic system has an equilibrium state, Xf , we say that the full stochasticprocess has a partial equilibrium, which is (Xs, Xf ).

We can assume the existence of partial stochastic equilibria if two conditions are sat-isfied: the existence of an equilibrium state for the fast systems, and the relaxation timeof the fast system being small in comparison with the time scale of the slow reactions.

Moreover we define the following new propensity density functions, in which the fastvariables are in the equilibrium.

asj(X) = E(asi (Xf )). (D.2)

D.1.2 Computation of Partial equilibrium

How to decide which are the fast and the slow reactions, in general, is an open question,and usually the best way to determine them is to run few simulations of the SSA simu-lation algorithm and decide. Another option is to try to obtain some information of thepropensity functions aj and the different time scales of the system.

Compute asj(Xs) is the most difficult part. There are 5 different cases [20]:

1. If asj(X) is independent of Xf , then asj(Xs) = asj(X).

2. If asj(X) = csjXfi , then asj(X

s) = csj(Xfi ).

3. If asj(X) = csjXfi x

sk, then asj(X

s) = csjXskE(Xf

i ).

4. If asj(X) = csjXfi (Xf

i − 1), then asj(Xs) =

csj2E(Xf

i (Xfi − 1)) ≈

csj2

(E(Xfi ))2.

5. If asj(X) = csjXfi X

fk , then asj(X) = csjE(Xf

i Xfk ) ≈ csjE(Xf

i )E(Xfk ).

Multi-scale Stochastic Simulations 107

The computation of the partial equilibrium requires write and solve equations for thebalance of the propensity functions. Consider that the fast system has the form,

A1 +A2c1−→ A1,

A2c2−→ A2 +A2.

In this case a1(X) = c1X1X2 and a2(X) = c2X2, and the conservation law is

c1X1X2 = c2X2. (D.3)

Remark: In general we do not obtain a linear system, however, if the partial equilib-rium exists, we can solve the equations of the partial equilibrium via the Newton method[96].

However, there problematic cases. Consider the system described by

A1 +A1c1−→ A2

A2c2−→ A1 +A1.

In this case The equilibrium law given by the deterministic equations is

c1

2X2

1 = c2X2 (D.4)

however, the MSSA requiresc1

2X1(X1 − 1) = c2X2. (D.5)

For a large X1 this is not problematic, however, when X1 is small this is no longervalid.

Example D.1. To illustrate this methodology, we use as example the stochastic processdescribed in Chapter 3. That is, the stochastic process described by the Eq. 3.1. As inChapter 3 we assume that the variable T (the number of healthy cells) follows an ordinarydifferential equation.

From our time-scale analysis we already know that the variable V , corresponding tothe number of virions, is the fast variable. But we can obtain the same running the SSAand counting the time that each variable occurs. From a few realizations, we observe thateither the reaction corresponding to the virion clearance or the reaction correspondingto the virus production are chosen around 99.9% of times. Therefore we choose as fastreactions:

T ∗pv−→ T ∗ + V

Vc−→ ∅

The fast variables are the variables altered by the fast reactions, in this case the onlyvariable which is altered by these two reactions is V , which account the number of thevirus. The equation for the partial equilibrium come from the balance of this two ratios:

pvT∗ = cE[V ]. (D.6)

Thus, partial equilibrium is V =pvT

∗

c.

108

D.1.3 Algorithm

The main idea of the algorithm is to produce an stochastic simulation for the slow speciesXs, assuming that the fast species, Xf , evolves faster and have reached the partial stochas-tic equilibrium.

The Multi-scale Stochastic Simulation Algorithm steps are

1. Compute the equilibrium state for the fast system and update Xf according to thecurrent state.

2. Compute asj , as0 =

∑asj .

3. generate r1 and r2 random numbers in [0, 1], and compute τ and j according to

• τ =1

as0(x)ln

(1

1− r1

),

• j = the smallest integer satisfying

j∑k=1

ask(x) > r2as0(x).

4. t = t+ τ and Xs = Xs + vsj

5. go to 1 or stop.

Bibliography

[1] M. Adimy and F. Crauste. Mathematical model of hematopoiesis dynamics withgrowth factor-dependent apoptosis and proliferation regulations. Mathematical andComputer Modelling, 49(11–12):2128 – 2137, 2009.

[2] T. Alarcon. Stochastic quasi-steady state approximations for asymptotic solutionsof the chemical master equation. J. Chem. Phys., 140:184109, 2014.

[3] T. Alarcon and H. J. Jensen. Quiescence: a mechanism for escaping the effectsof drug on cell populations. Journal of The Royal Society Interface, 8(54):99–106,2011.

[4] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter. MolecularBiology of the Cell. Garland Publishing, New York, 5th edition, 2008.

[5] D. F. Anderson. A modified next reaction method for simulating chemical systemswith time dependent propensities and delays. The Journal of chemical physics,127(21):214107, 2007.

[6] V. I. Arnol’d. Mathematical methods of classical mechanics, volume 60. Springer,1989.

[7] J. Ashbourn, J. Miller, V. Reumers, V. Baekelandt, and L. Geris. A mathematicalmodel of adult subventricular neurogenesis. Journal of The Royal Society Interface,9(75):2414–2423, 2012.

[8] M. Assaf and B. Meerson. Spectral theory of metastability and extinction in birth-death systems. Physical review letters, 97(20):200602, 2006.

[9] M. Assaf and B. Meerson. Spectral theory of metastability and extinction in abranching-annihilation reaction. Physical Review E, 75(3):031122, 2007.

[10] M. Assaf and B. Meerson. Extinction of metastable stochastic populations. Phys.Rev. E, 81:021116, 2010.

[11] C. Batut, K. Belabas, D. Bernardi, H. Cohen, and M. Olivier. Users’ guide toPARI/GP. http://pari.math.u-bordeaux.fr/.

[12] H. C. Berg. Random walks in biology. Princeton University Press, Princeton, NJ,USA, 1993.

[13] D. Bernstein. Simulating mesoscopic reaction-diffusion systems using the Gillespiealgorithm. Phys. Rev. E, 71:041103, 2005.

109

110

[14] R. P. Boland, T. Galla, and A. J. McKane. Limit cycles, complex floquet multipliers,and intrinsic noise. Physical Review E, 79(5):051131, 2009.

[15] D. Bray. Cell movements: From molecules to motility. Garland Publishing, NewYork, NY, USA, 2001.

[16] S. Buczacki, R. Davies, and D. Winton. Stem cells, quiescence and rectal carci-noma: an unexplored relationship and potential therapeutic target. British journalof cancer, 105(9):1253–1259, 2011.

[17] X. Cabre, E. Fontich, and R. de la Llave. The parameterization method for invariantmanifolds. I. Manifolds associated to non-resonant subspaces. Indiana Univ. Math.J., 52(2):283–328, 2003.

[18] X. Cabre, E. Fontich, and R. de la Llave. The parameterization method for invari-ant manifolds. II. Regularity with respect to parameters. Indiana Univ. Math. J.,52(2):329–360, 2003.

[19] D. S. Callaway and A. S. Perelson. HIV-1 infection and low steady state viral loads.Bulletin of Mathematical Biology, 64(1):29–64, 2002.

[20] Y. Cao, D. Gillespie, and L. Petzold. Multiscale stochastic simulation algorithm withstochastic partial equilibrium assumption for chemically reacting systems. Journalof Computational Physics, 206(2):395–411, 2005.

[21] Y. Cao, H. Li, and L. Petzold. Efficient formulation of the stochastic simulation algo-rithm for chemically reacting systems. The journal of chemical physics, 121(9):4059–4067, 2004.

[22] J. Cerda. Linear functional analysis, volume 116. American Mathematical Soc.,2010.

[23] N. Chom*ont, M. El-Far, P. Ancuta, L. Trautmann, F. A. Procopio, B. Yassine-Diab,G. Boucher, M.-R. Boulassel, J. M. Brenchley, T. W. Schacker, B. J. Hill, D. C.D. G. Ghattas, J.-P. Routy, E. K. Haddad, and R.-P. Sekaly. HIV reservoir size andpersistence are driven by T cell survival and homeostatic proliferation. Nature Med.,15:893–901, 2009.

[24] T.-W. Chun, L. Carruth, D. Finzi, X. Shen, J. A. DiGiuseppe, H. Taylor, M. Her-mankova, K. Chadwick, J. Margolick, T. C. Quinn, Y.-H. Kuo, R. Brookmeyer,M. A. Zeiger, P. Barditch-Crovo, and R. F. Siliciano. Quantification of latent tissuereservoirs and total body viral load in HIV-1 infection. Nature, 387:183–188, 1997.

[25] T. W. Chun, D. Finzi, J. Margolick, K. Chadwick, D. Schwartz, and R. F. Siliciano.In vivo fate of HIV-1-infected T cells: Quantitative analysis of the transition tostable latency. Nat. Med., 1:1284–1290, 1995.

[26] C. Colijn and M. C. Mackey. A mathematical model of hematopoiesis–I. Periodicchronic myelogenous leukemia. Journal of Theoretical Biology, 237(2):117–132, 2005.

[27] C. Colijn and M. C. Mackey. A mathematical model of hematopoiesis–II. Cyclicalneutropenia. Journal of theoretical biology, 237(2):133–146, 2005.

111

[28] J. M. Conway and D. Coombs. A Stochastic Model of Latently Infected Cell Re-activation and Viral Blip Generation in Treated HIV Patients. PLoS Comput Biol,7:e1002033, 04 2011.

[29] F. Crauste. Delay model of hematopoietic stem cell dynamics: asymptotic stabilityand stability switch. Mathematical Modelling of Natural Phenomena, 4(02):28–47,2009.

[30] M. DiMascio, M. Markowitz, M. Louie, C. Hogan, A. Hurley, C. Chung, D. D.Ho, and A. S. Perelson. Viral Blip Dynamics during Highly Active AntiretroviralTherapy. J Virol, 77(22):12165–12172, 2003.

[31] D. Dingli, A. Traulsen, and F. Michor. (a) symmetric stem cell replication andcancer. PLoS computational biology, 3(3):e53, 2007.

[32] D. Dingli, A. Traulsen, and J. M. Pacheco. Compartmental architecture and dy-namics of hematopoiesis. PLoS One, 2(4):e345, 2007.

[33] G. Dornadula, H. Zhang, B. VanUitert, and et al. Residual HIV-1 RNA in bloodplasma of patients taking suppressive highly active antiretroviral therapy. JAMA,282(17):1627–1632, 1999.

[34] J. Elf and M. Ehrenberg. Spontaneous separation of bi-stable biochemical systemsinto spatial domains of opposite phases. Syst. Biol., 1:230–235, 2004.

[35] V. Elgart and A. Kamenev. Rare event statistics in reaction-diffusion systems. Phys.Rev. E, 70:041106, 2004.

[36] R. Erban, J. Chapman, and P. Maini. A practical guide to stochastic simulations ofreaction-diffusion processes. arXiv preprint arXiv:0704.1908, 2007.

[37] N. M. Ferguson, F. deWolf, A. C. Ghani, C. Fraser, C. A. Donnelly, P. Reissi,J. M. A. Langei, S. A. Danneri, G. P. Garnett, J. Goudsmit, , and R. M. Anderson.Antigen-driven CD41 T cell and HIV-1 dynamics: Residual viral replication underhighly active antiretroviral therapy. Proc. Natl. Acad. Sci., 96:15167–15172, 1999.

[38] E. Fermi. Thermodynamics, 1936. Thermodynamics, 1936, 1956.

[39] J. Foo, M. W. Drummond, B. Clarkson, T. Holyoake, and F. Michor. Eradica-tion of chronic myeloid leukemia stem cells: a novel mathematical model predictsno therapeutic benefit of adding G-CSF to imatinib. PLoS computational biology,5(9):e1000503, 2009.

[40] C. Fraser, N. M. Ferguson, and R. M. Anderson. Quantification of intrinsic residualviral replication in treated HIV infected patients. Proc. Natl. Acad. Sci., 98:15167–15172, 2001.

[41] C. Frazer, N. M. Ferguson, F. de Wolf, and R. M. Anderson. The role of antigenicstimulation and cytotoxic T cel activity in regulating the long-term immunopatho-genesis of HIV: Mechanisms and clinical implications. Proc. R. Soc. London B.,268:2085–2095, 2001.

112

[42] M. I. Freidlin and A. D. Wentzell. Random perturbations of dynamical systems, vol-ume 260 of grundlehren der mathematischen wissenschaften [fundamental principlesof mathematical sciences], 1998.

[43] L. M. Frenkel, Y. Wang, G. H. Learn, J. L. McKernan, G. M. Ellis, K. M. Mohan,S. E. Holte, S. M. D. Vange, D. M. Pawluk, A. J. Melvin, P. F. Lewis, L. M.Heath, I. A. Beck, M. Mahalanabis, W. E. Naugler, N. H. Tobin, and J. I. Mullinhs.Multiple viral genetic analyses detect low-level human immunodeficency virus typeI replication during effective highly active antiretroviral therapy. J. Virol., 77:5721–5730, 2003.

[44] C. W. Gardiner. Stochatic methods. Springer-Verlag, Berlin, Germany, 2009.

[45] V. Gelfreich, C. Simo, and A. Vieiro. Dynamics of 4D symplectic maps near a doubleresonance. Physica D, 243(1):92–110, 2013.

[46] M. A. Gibson and J. Bruck. Efficient exact stochastic simulation of chemical systemswith many species and many channels. J. Phys. Chem. A, 104:1876–1889, 2000.

[47] D. T. Gillespie. A general method for numerically simulating the stochastic timeevolution of coupled chemical reactions. Journal of computational physics, 22(4):403–434, 1976.

[48] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. Thejournal of physical chemistry, 81(25):2340–2361, 1977.

[49] D. T. Gillespie. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem.,58:35–55, 2007.

[50] M. Gotz and W. B. Huttner. The cell biology of neurogenesis. Nature reviewsMolecular cell biology, 6(10):777–788, 2005.

[51] A. Griewank. Evaluating Derivatives: Principles and Techniques of AlgorithmicDifferentiation. Number 19 in Frontiers in Appl. Math. SIAM, Philadelphia, PA,2000. Technical report, ISBN 08-987-1451-6.

[52] H. F. Gunthard, J. Wong, C. A. Spina, C. Ignacio, S. Kwok, C. Christopherson,J. Hwang, R. Haubrich, D. Havlir, and D. D. Richman. Effect of influenza vaccinationon viral replication and immune response in persons infected with human immun-odeficency virus receiving potent antiretroviral therapy. J. Infect. Dis., 181:522–531,2000.

[53] A. Haro, M. Canadell, J. Figueras, A. Luque, and J. Mondelo. The parameterizationmethod for invariant manifolds: from rigorous results to effective comptutations.http://www.maia.ub.es/ alex/review/review.pdf/.

[54] D. V. Havlir, R. Bassett, D. Levitan, P. Gilbert, P. Tebas, A. C. Collier, M. S.Hirsch, C. Ignacio, J. Condra, H. F. Gunthard, D. D. Richman, and J. K. Wong.Prevalence and Predictive Value of Intermittent Viremia With Combination HIVTherapy. JAMA, 286:171–179, 2001.

113

[55] A. V. M. Herz, S. Bonhoefer, R. M. Anderson, R. M. May, and M. A. Nowak. Viraldynamics in vivo: limitations on estimates of intracellular delay and virus decay.Proc. Natl. Acad. Sci., 93:7247–7251, 1996.

[56] B. Hille. Ionic channels of excitable membranes. Sinauer Associates, Sunderland,MA, USA, 1992.

[57] D. D. Ho, A. U. Neumann, A. S. Perelson, W. Chen, J. M. Leonard, andM. Markowitz. Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature, 373:123–126, 1995.

[58] D. D. Ho, T. R. Rota, and M. S. Hirsch. Infection of monocyte/macrophages byhuman T lymphotropic virus type III. J. CLin. Invest., 77:1712–1715, 1986.

[59] R. D. Hockett, J. M. Kilby, C. A. Derdeyn, M. S. Saag, M. Sillers, K. Squires,S. Chiz, M. A. Nowak, G. M. Shaw, and R. P. Bucy. Constant mean viral copynumber per infected cell in tissues regardless of high, low, or undetectable plasmaHIV RNA. The Journal of experimental medicine, 189(10):1545–1554, 1999.

[60] M. S. i Sole. Edicions universitat de barcelona, 1999.

[61] L. E. Jones and A. S. Perelson. Modelling the effects of vaccination on chronicallyinfected HIV-positive patients. J. Acquir. Immune Defic. Syndr., 31:369–377, 2002.

[62] L. E. Jones and A. S. Perelson. Opportunistic infections as a cause of transientviremia in chronically infected HIV patients under treatment with HAART. Math.Biosci., 67:1227–1251, 2005.

[63] L. E. Jones and A. S. Perelson. Transient viremia, plasma viral load and reservoirreplenishment in HIV-infected patients on antiretroviral therapy. J. Acquir. ImmuneDefic. Syndr., 45:483–493, 2007.

[64] A. Jorba and M. Zou. A software package for the numerical integration of ODEs bymeans of high-order Taylor methods. 14(1):99–117, 2005.

[65] C. Katlama, S. G. Deeks, B. Autran, J. Martinez-Picado, J. van Luzen, C. Rouzioux,M. Miller, S. Vella, J. E. Schmitz, J. Ahlers, D. D. Richman, and R. P. Sekaly.Barriers to a cure for HIV: new ways to target and eradicate HIV-1 reservoirs.Lancet, 381:2109–2117, 2013.

[66] S. J. Kent, J. C. Reece, J. Petravic, A. Martyushev, M. Kramski, R. D. Rose,D. A. Cooper, A. D. Kelleher, S. Emery, P. U. Cameron, S. R. Lewin, and M. P.Davenport. The search for an HIV cure: tackling latent infection. Lancet Infect.Dis., 13:614–621, 2013.

[67] T. B. Kepler and A. S. Perelson. Drug concentration heterogeneityof drug resistancefacilitates the evolution of drug resistance. Proc. Natl. Acad. Sci., 95:11514–11519,1998.

[68] D. A. Kessler and N. M. Schnerb. Extinction rates for fluctuation-induced metasta-bilities: A real-space WKB approach. J. Stat. Phys., 127:861–886, 2007.

114

[69] H. Kim and A. S. Perelson. Viral and latent reservoir persistence in HIV-1-infectedpatients on therapy. PLoS Comp. Biol., 2:e135, 2006.

[70] A. Klinakis, C. Lobry, O. Abdel-Wahab, P. Oh, H. Haeno, S. Buonamici, I. vanDe Walle, S. Cathelin, T. Trimarchi, E. Araldi, et al. A novel tumour-suppressorfunction for the notch pathway in myeloid leukaemia. Nature, 473(7346):230–233,2011.

[71] J. A. Knoblich. Mechanisms of asymmetric stem cell division. Cell, 132(4):583–597,2008.

[72] J. A. Knoblich. Asymmetric cell division: recent developments and their implicationsfor tumour biology. Nature reviews Molecular cell biology, 11(12):849–860, 2010.

[73] R. Kubo, K. Matsou, and K. Kitahara. Fluctuations and relaxation of macrovari-ables. J. Stat. Phys., 9:51–96, 1973.

[74] T. Lechler and E. Fuchs. Asymmetric cell divisions promote stratification and dif-ferentiation of mammalian skin. Nature, 437(7056):275–280, 2005.

[75] K. Leder, J. Foo, B. Skaggs, M. Gorre, C. L. Sawyers, and F. Michor. Fitnessconferred by bcr-abl kinase domain mutations determines the risk of pre-existingresistance in chronic myeloid leukemia. PloS one, 6(11):e27682, 2011.

[76] P. K. Lee, T. L. Kieffer, R. F. Siciliano, and R. E. Nettles. HIV-1 viral load blipsare of limited clinical significance. J. Antimicrob. Chemother., 57:803–805, 2006.

[77] J. Lei and M. C. Mackey. Multistability in an age-structured model of hematopoiesis:Cyclical neutropenia. Journal of theoretical biology, 270(1):143–153, 2011.

[78] T. Lenaerts, J. M. Pacheco, A. Traulsen, and D. Dingli. Tyrosine kinase inhibitortherapy can cure chronic myeloid leukemia without hitting leukemic stem cells.haematologica, 95(6):900–907, 2010.

[79] C. Lugo and A. J. McKane. Quasicycles in a spatial predator-prey model. Phys.Rev. E, 78:051911, 2008.

[80] J. Macias, J. C. Palomares, J. A. Mira, M. J. Torres, J. A. Garcia-Garcia, J. M.Rodriguez, S. Vergera, and J. A. Pineda. Transient rebounds of HIV plasma viremiaare associated with the emergence of drug resistance mutations in patients of highlyactive antiretroviral therapy-. J. Infection, 5:195–200, 2005.

[81] J. M. Mahaffy, J. Belair, and M. C. Mackey. Hematopoietic model with movingboundary condition and state dependent delay: applications in erythropoiesis. Jour-nal of Theoretical Biology, 190(2):135–146, 1998.

[82] A. Marciniak-Czochra, T. Stiehl, A. D. Ho, W. Jager, and W. Wagner. Modelingof asymmetric cell division in hematopoietic stem cells-regulation of self-renewal isessential for efficient repopulation. Stem cells and development, 18(3):377–386, 2009.

[83] M. Markowitz, M. Louie, A. Hurley, E. Sun, M. D. Mascio, A. S. Perelson, and D. D.Ho. A Novel Antiviral Intervention Results in More Accurate Assessment of HumanImmunodeficiency Virus Type 1 Replication Dynamics and T-Cell Decay In Vivo.J. Virol., 77:5037–5038, 2003.

115

[84] J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson, and N. F. Samatova.The sorting direct method for stochastic simulation of biochemical systems withvarying reaction execution behavior. Computational biology and chemistry, 30(1):39–49, 2006.

[85] A. J. McKane and T. J. Newman. Predator-prey cycles from resonant amplificationof demographic stochasticity. Physical review letters, 94(21):218102, 2005.

[86] F. Michor. Quantitative approaches to analyzing imatinib-treated chronic myeloidleukemia. Trends in pharmacological sciences, 28(5):197–199, 2007.

[87] F. Michor, T. P. Hughes, Y. Iwasa, S. Branford, N. P. Shah, C. L. Sawyers, and M. A.Nowak. Dynamics of chronic myeloid leukaemia. Nature, 435(7046):1267–1270, 2005.

[88] F. Michor, Y. Iwasa, and M. A. Nowak. The age incidence of chronic myeloidleukemia can be explained by a one-mutation model. Proceedings of the NationalAcademy of Sciences, 103(40):14931–14934, 2006.

[89] J. A. Mira, J. Macias, C. Nogales, J. Fernandez-Rivera, J. A. Garcia-Garcia,A. Ramos, and J. A. Pineda. Transient rebounds of low-level viraemia among HIV-infected patients under HAART are not associated with virological or immunologicalfailure. Antivir. Ther., 7:251–256, 2002.

[90] H. Mohri, S. Bonhoeffer, S. Monard, A. S. Perelson, and D. D. Ho. Rapid turnoverof T lymphocytes in SIV-infected rhesus macaques. Science, 279(5354):1223–1227,1998.

[91] S. J. Morrison and J. Kimble. Asymmetric and symmetric stem-cell divisions indevelopment and cancer. Nature, 441(7097):1068–1074, 2006.

[92] S. J. Morrison and A. C. Spradling. Stem cells and niches: mechanisms that promotestem cell maintenance throughout life. Cell, 132(4):598–611, 2008.

[93] A. G. Murray and G. A. Jackson. Viral dynamics: a model of the effects of size shape,motion and abundance of single-celled planktonic organisms and other particles.Mar. Ecol. Prog. Ser., 89:103–116, 1992.

[94] R. E. Nettles and T. L. Kieffer. Update on HIV-1 viral load blips. Current Opinionin HIV and AIDS, 1:157–161, 2006.

[95] R. E. Nettles, T. L. Kieffer, P. K. andDaphne Monie, Y. Han, T. Parsons,J. Cofrancesco, J. E. Gallant, T. C. Quinn, B. Jackson, C. Flexner, K. Carson,S. Ray, D. Persaud, and R. F. Siliciano. Intermittent hiv-1 viremia (blips) and drugresistance in patients receiving haart. JAMA, 293:817–829, 2005.

[96] J. Nocedal and S. J. Wright. Numerical optimization 2nd. 2006.

[97] J. E. Pearson, P. Krapivsky, and A. S. Perelson. Stochastic theory of early viral infec-tion: Continuous versus burst production of virions. PLoS Comp. Biol., 7:e1001058,2011.

[98] J. W. Pepper, K. Sprouffske, and C. C. Maley. Animal cell differentiation patternssuppress somatic evolution. PLoS computational biology, 3(12):e250, 2007.

116

[99] A. S. Perelson, P. Essunger, Y. Cao, M. Vesanen, A. Hurley, K. Saksela,M. Markowitz, and D. D. Ho. Decay characteristics of HIV-1-infected compartmentsduring combination therapy. Nature, 387:188–191, 1997.

[100] A. S. Perelson, D. E. Kirschner, and R. De Boer. Dynamics of HIV infection ofCD4+ T cells. Mathematical biosciences, 114(1):81–125, 1993.

[101] A. S. Perelson, A. U. Neumann, M. Markowitz, J. M. Leonard, and D. D. Ho. HIV-1dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generationtime. Science, 271:1582–1586, 1996.

[102] T. Pierson, J. McArthur, and R. F. Siciliano. Reservoirs for HIV-1: Mechanismsfor viral persistence in the presence of antiviral immune response and antiretroviraltherapy. Annu. Rev. Immunol., 18:665–708, 2000.

[103] R. Ramaswamy, N. Gonzalez-Segredo, and I. F. Sbalzarini. A new class of highlyefficient exact stochastic simulation algorithms for chemical reaction networks. TheJournal of chemical physics, 130(24):244104, 2009.

[104] B. Ramratnam, S. Bonhoeffer, J. Binley, A. Hurley, L. Zhang, and et al. Rapidproduction and clearance of HIV-1 and hepatitis C virus assessed by large volumeplasma apheresis. Lancet, 354:1782–1785, 1999.

[105] D. D. Richman, D. M. Margolis, M. Delaney, W. C. Greene, D. Hazuda, and R. J.Pomerantz. The challenge of finding a cure for HIV infection. Science, 323:1304–1307, 2009.

[106] I. A. Rodriguez-Brenes, D. Wodarz, and N. L. Komarova. Asymmetric cell division:Recent developments and their implications for tumour biology. Front. Oncol., 3:82,2013.

[107] I. A. Rodriguez-Brenes, D. Wodarz, and N. L. Komarova. Minimizing the risk ofcancer: tissue architecture and cellular replication limits. Journal of The RoyalSociety Interface, 10(86):20130410, 2013.

[108] A. Roesch, M. f*ckunaga-Kalabis, E. C. Schmidt, S. E. Zabierowski, P. A. Brafford,A. Vultur, D. Basu, P. Gimotty, T. Vogt, and M. Herlyn. A temporarily distinct sub-population of slow-cycling melanoma cells is required for continuous tumour growth.Cell, 141:583–594, 2010.

[109] L. Rong and A. S. Perelson. Asymmetric cell division of activated latently infectedcells may explain the decay kinetics of the HIV-1 latent reservoir and intermittentviral blips. Math. Biosci., 217:77–87, 2009.

[110] L. Rong and A. S. Perelson. Modeling Latently Infected Cell Activation: Viral andLatent Reservoir Persistence, and Viral Blips in HIV-infected Patients on PotentTherapy. PLoS Comput. Biol., 5:e1000533, 10 2009.

[111] L. Rong and A. S. Perelson. Modelling HIV persistence, the latent reservoir, andviral blips. J. Theor. Biol., 260:308–331, 2009.

[112] N. J. Savill. Mathematical models of hierarchically structured cell populations underequilibrium with application to the epidermis. Cell proliferation, 36(1):1–26, 2003.

117

[113] L. Shan, K. Deng, N. S. Shroff, C. M. Durand, S. A. Rabi, H.-C. Yang, H. Zhang,J. B. Margolick, J. N. Blankson, and R. F. Siciliano. Stimulation of HIV-1-specificcytolityc T lymphocytes facilitates elimination of latent viral reservoir after virusreactivation. Immunity, 36:491–501, 2012.

[114] H. Shenghui, D. Nakada, and S. J. Morrison. Mechanisms of stem cell self-renewal.Annual Review of Cell and Developmental, 25:377–406, 2009.

[115] C. Simo. On the analytical and numerical approximation of invariant manifolds.In Les Methodes Modernes de la Mecanique Celeste. Modern Methods in CelestialMechanics, volume 1, pages 285–329, 1990.

[116] A. Slepoy, A. P. Thompson, and S. J. Plimpton. A constant-time kinetic montecarlo algorithm for simulation of large biochemical reaction networks. The journalof chemical physics, 128(20):205101, 2008.

[117] F. Spill, P. Guerrero, T. Alarcon, P. K. Maini, and H. M. Byrne. Mesoscopic andcontinuum modelling of angiognesis. J. Math. Biol., page To appear., 2014.

[118] T. Stiehl and A. Marciniak-Czochra. Characterization of stem cells using mathe-matical models of multistage cell lineages. Mathematical and Computer Modelling,53(7):1505–1517, 2011.

[119] Z. Sun and N. L. Komarova. Stochastic modeling of stem-cell dynamics with control.Mathematical biosciences, 240(2):231–240, 2012.

[120] A. Traulsen, T. Lenaerts, J. M. Pacheco, and D. Dingli. On the dynamics of neutralmutations in a mathematical model for a hom*ogeneous stem cell population. Journalof The Royal Society Interface, 10(79):20120810, 2013.

[121] D. Trono, C. V. Lint, C. Rouziux, E. Verdin, F. Barre-Sinoussi, T.-W. Chun, andN. Chom*ont. HIV persistence and the prospect of long-term drug-free remissionsfor HIV-infected individuals. Science, 329:174–180, 2010.

[122] N. G. Van Kampen. Stochastic processes in Physics and Chemistry. Elsevier, TheNetherlands, 2007.

[123] S. Wang and L. Rong. Stochastic population switch may explain the latent reservoirstability and intermittent viral blips in HIV patients on suppresive. J. Theor. Biol.,360:137–148, 2014.

[124] X. Wei, S. K. Ghosh, M. E. Taylor, V. A. Johnson, E. A. Emini, P. Deutsch, J. D.Lifson, S. Bonhoeffer, M. A. Nowak, B. H. Hahn, M. S. Saag, and G. M. Shaw. Viraldynamics in human immunodeficiency virus type 1 infection. Nature, 373:117–122,1995.

[125] B. Werner, D. Dingli, and A. Traulsen. A deterministic model for the occurrenceand dynamics of multiple mutations in hierarchically organized tissues. Journal ofThe Royal Society Interface, 10(85):20130349, 2013.

118

[126] Z.-Q. Zhang, T. Schuler, M. Zupancic, S. Wietgrefe, K. A. Staskus, K. A. Reimann,T. A. Reinhart, M. Rogan, W. Cavert, C. J. Miller, R. S. Veazey, D. Notermans,S. Little, S. A. Danner, D. D. Richman, D. Havlir, J. Wong, H. L. Jordan, T. W.Schacker, P. Racz, K. Tenner-Racz, N. L. Letvin, S. Wolinsky, and A. T. Haase.Sexual Transmission and Propagation of SIV and HIV in Resting and ActivatedCD4+ T Cells . Science, 286:1712–1715, 1999.

[127] R. Zhao and F. Michor. Patterns of proliferative activity in the colonic crypt de-termine crypt stability and rates of somatic evolution. PLoS computational biology,9(6):e1003082, 2013.