Mathematical Models for Immunodeficiency Virus, Post-treatment and Memory Activation by Alejandra Donaji Herrera Reyes A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2012 c© Alejandra Donaji Herrera Reyes 2012 Abstract Nowdays, HIV infection can be controlled by anti-retroviral drug therapy (ART). However, a persistent viral reservoir in treated patients prevents the eradication of HIV infection. H-iART is an innovator treatment that con- sists of regular ART and the drugs Maraviroc and Darunavir, and H-iART was enforced with Auranofin. The drug Maraviroc (MRV) was proved to be a good CCR5 inhibitor, which is a HIV correceptor. The drug Auranofin has been shown to accelerate the activation rate of latent cells and also alters the kinetics of viral rebound when drug treatment is interrupted. Recent studies on monkeys infected with SIV have shown a complete suppression of the viral load during H-iART with Auranofin treatment and a persis- tent suppression of it in the absence of ART. Motivated by the results of the experiments I present deterministic and stochastic models of HIV after treatment interruption. For H-iART treatment, the ODE models were used as a start point to create three different continuous time multi-type branch- ing process. From equations for the probability generating function we use analytic solutions, numerical approximations, or numerical simulations to extract the probability of observable viral blips. We compare our results with the data of two rhesus macaques. We find that more than one latent cell needs to activate in order to observe the data blips, and that the net re- productive number of virus must be very close to one. Since this is unlikely, these results suggest that the viral dynamics must be more complex than our model allows for. For the ART+Auranofin treatment, I will present an ODE model of HIV population dynamics including drug treatment and the immune response to model the viral rebound at treatment interruption. ii Preface Chapter 2 is based on data by I.L. Shytaj, A. Savarino and others. The data was published on [19] and [18] and was given to us by personal communi- cation. I performed all the mathematical modeling described in the cited chapter. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Post treatment model . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Data summary of P177 and P157 monkeys . . . . . . . . . . 4 2.3 ODE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Introduction to the stochastic approach . . . . . . . . . . . . 15 2.5 T ∗ stochastic model . . . . . . . . . . . . . . . . . . . . . . . 15 2.5.1 Kolmogorov’s forward differential equations . . . . . . 17 2.5.2 Kolmogorov’s backward differential equations . . . . 18 2.5.3 Probability generating function (pgf) definition . . . 19 2.5.4 Applying the probability generating function to the backward equation (2.9) . . . . . . . . . . . . . . . . 20 2.5.5 The probability generating function applied to for- ward equation (2.6) . . . . . . . . . . . . . . . . . . . 21 2.5.6 Transition probabilities . . . . . . . . . . . . . . . . . 23 2.5.7 Extinction probability . . . . . . . . . . . . . . . . . . 25 2.5.8 Limit analysis of the transition probabilities . . . . . 26 2.5.9 First passage time to extinction . . . . . . . . . . . . 28 iv Table of Contents 2.5.10 Numerical simulations . . . . . . . . . . . . . . . . . . 31 2.5.11 Comparing results with data of P157 and P177 mon- keys . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.5.12 Discussion of the model . . . . . . . . . . . . . . . . . 38 2.6 Two-compartment stochastic model . . . . . . . . . . . . . . 40 2.6.1 Kolmogorov’s forward differential equations . . . . . . 41 2.6.2 Kolmogorov’s backward differential equations . . . . 41 2.6.3 New properties of probability generation function in more than one dimension . . . . . . . . . . . . . . . . 42 2.6.4 The probability generating function for the Kolmogo- rov forward equation . . . . . . . . . . . . . . . . . . 43 2.6.5 The probability generating function for the Kolmogo- rov backward equation . . . . . . . . . . . . . . . . . 45 2.6.6 Calculating the individual probabilities . . . . . . . . 46 2.6.7 Numerical simulations . . . . . . . . . . . . . . . . . . 47 2.6.8 Comparing the viral load from simulation with data . 51 2.6.9 Discussion of the model . . . . . . . . . . . . . . . . . 52 2.7 Latent cells stochastic model . . . . . . . . . . . . . . . . . . 53 2.7.1 Numerical simulations . . . . . . . . . . . . . . . . . . 56 2.7.2 Comparing with data . . . . . . . . . . . . . . . . . . 61 2.7.3 Numerical simulation with fast activation rate . . . . 63 2.7.4 Discussion of the model . . . . . . . . . . . . . . . . . 69 2.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3 Immune response and memory activation . . . . . . . . . . 72 3.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . 72 3.2 Analytical work . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.3 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . 76 3.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . 78 3.5 Other parameter dynamics . . . . . . . . . . . . . . . . . . . 87 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 v List of Tables 2.1 Summary of viral blip data . . . . . . . . . . . . . . . . . . . 6 2.2 Parameters to the classic ODE model . . . . . . . . . . . . . 11 2.3 Parameters to the latent ODE model . . . . . . . . . . . . . . 14 2.4 Summary of viral blip data with T ∗ converted values . . . . . 34 2.5 Probabilities of blips with parameters as in Table 2.2 . . . . . 35 2.6 Probabilities of blips with new parameters . . . . . . . . . . . 39 2.7 Probabilities of small blips with parameters as in Table 2.2 . 52 3.1 Parameter description of model [3.1] . . . . . . . . . . . . . . 75 3.2 Parameters of the model as we have found them in literature 79 vi List of Figures 2.1 P177 viral load per ml, log10 scale . . . . . . . . . . . . . . . 4 2.2 P157 viral load per ml, log10 scale . . . . . . . . . . . . . . . 4 2.3 P177 CD4 counts per µl . . . . . . . . . . . . . . . . . . . . . 5 2.4 P157 CD4 counts per µl . . . . . . . . . . . . . . . . . . . . . 5 2.5 Zoomed-in plots of the blips observed on P177 . . . . . . . . 7 2.6 Zoomed-in plots of the blips observed on P157 . . . . . . . . 8 2.7 Duration of each viral blip on P177 . . . . . . . . . . . . . . . 9 2.8 Duration of each viral blip on P157 . . . . . . . . . . . . . . . 9 2.9 High of each viral blip on P177 . . . . . . . . . . . . . . . . . 9 2.10 High of each viral blip on P157 . . . . . . . . . . . . . . . . . 9 2.11 T cells solution of classic model . . . . . . . . . . . . . . . . . 12 2.12 Infected cells solution of classic model . . . . . . . . . . . . . 12 2.13 Viral load per monkey solution of classic model . . . . . . . . 12 2.14 Viral load per ml solution of classic model . . . . . . . . . . . 12 2.15 Viral load solution of latent model with IC L = 1 . . . . . . . 14 2.16 Viral load solution of latent model with IC T ∗ = 1 . . . . . . 14 2.17 T ∗V simplified model . . . . . . . . . . . . . . . . . . . . . . 15 2.18 T ∗ simplified model . . . . . . . . . . . . . . . . . . . . . . . . 15 2.19 Probability distribution of model T ∗ at time 7 days . . . . . . 24 2.20 Probability to have 5 infected cells as function of time . . . . 24 2.21 Extinction probability function starting with T ∗(0) = 1, for T ∗ model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.22 Probability that T ∗(7) = 0 for different initial condition for T ∗ model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.23 Probability density function of first time to extinction . . . . 29 2.24 Different paths from the T ∗ model Gillespie simulation . . . . 32 2.25 Maximum number of infected cells at different times . . . . . 32 2.26 Mean of the paths and the analytic solution of the expectation 32 2.27 Variance of the simulations . . . . . . . . . . . . . . . . . . . 32 2.28 First time to extinction from simulated data . . . . . . . . . . 33 vii List of Figures 2.29 Different paths from the T ∗ model numerical simulation, new parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.30 The longest 5 paths of the T ∗ model with new parameters . . 36 2.31 Mean of the paths with new parameters . . . . . . . . . . . . 37 2.32 Variance of the simulations with new parameters . . . . . . . 37 2.33 Maximum with new parameters but excluding longest paths . 37 2.34 Maximum number of infected cells with new parameters . . . 37 2.35 First time of extinction from numerical simulations . . . . . . 38 2.36 Different paths of T ∗ from the T ∗V model numerical simulation 48 2.37 Different paths of V from the T ∗V model numerical simulation 48 2.38 Maximum of T ∗ from the T ∗V model numerical simulation . 48 2.39 Maximum of V from the T ∗V model numerical simulation . . 48 2.40 Mean values of T ∗ from the T ∗V model numerical simulation 49 2.41 Mean values of V from the T ∗V model numerical simulation . 49 2.42 Variance of T ∗ from the T ∗V model numerical simulation . . 50 2.43 Variance of V from the T ∗V model numerical simulation . . . 50 2.44 First time to extinction of the system from the T ∗V model numerical simulation . . . . . . . . . . . . . . . . . . . . . . . 50 2.45 Histogram of viral loads in simulation at t = 20 days . . . . . 50 2.46 Probability of having more than 20000 virus per monkey as a function of time . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.47 Latent cells stochastic model . . . . . . . . . . . . . . . . . . 53 2.48 Different paths of T ∗ from the LT ∗V model Gillespie simula- tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.49 Different paths of V from the LT ∗V model Gillespie simulations 57 2.50 Different paths of L from the LT ∗V model Gillespie simulations 57 2.51 Maximum of T ∗ from the LT ∗V model Gillespie simulations . 58 2.52 Maximum of V from the LT ∗V model Gillespie simulations . 58 2.53 Maximum of L from the LT ∗V model Gillespie simulations . 58 2.54 Mean of T ∗ from the LT ∗V model Gillespie simulations . . . 59 2.55 Mean of V from the LT ∗V model Gillespie simulations . . . . 59 2.56 Mean of L from the LT ∗V model Gillespie simulations . . . . 59 2.57 Variance of T ∗ from the LT ∗V model Gillespie simulations . 60 2.58 Variance of V from the LT ∗V model Gillespie simulations . . 60 2.59 Variance of L from the LT ∗V model Gillespie simulations . . 60 2.60 First passage time to extinction of LT ∗V model . . . . . . . . 61 2.61 Cumulative distribution of the FPT to extinction of LT ∗V model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.62 LT ∗V model paths showing special dynamics, T ∗ plotted . . 62 2.63 LT ∗V model paths showing special dynamics, V plotted . . . 62 viii List of Figures 2.64 LT ∗V model paths showing special dynamics, L plotted . . . 62 2.65 Probability of having more than 20000 virus per monkey as a function of time . . . . . . . . . . . . . . . . . . . . . . . . . 63 2.66 Different paths of T ∗ from the LT ∗V model simulations with a = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.67 Different paths of V from the LT ∗V model Gillespie simula- tions with a = 10 . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.68 Different paths of L from the LT ∗V model simulations with a = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.69 Maximum of T ∗ from the LT ∗V model simulations with a = 10 65 2.70 Maximum of V from the LT ∗V model simulations with a = 10 65 2.71 Maximum of L from the LT ∗V model simulations with a = 10 65 2.72 Mean of T ∗ from the LT ∗V model simulations with a = 10 . 66 2.73 Mean of V from the LT ∗V model simulations with a = 10 . . 66 2.74 Mean of L from the LT ∗V model simulations with a = 10 . . 67 2.75 Variance of T ∗ from the LT ∗V model simulations with a = 10 67 2.76 Variance of V from the LT ∗V model simulations with a = 10 67 2.77 Variance of L from the LT ∗V model simulations with a = 10 68 2.78 First passage time to extinction of LT ∗V model . . . . . . . . 69 2.79 Cumulative distribution of the FPT to extinction of LT ∗V model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 2.80 Probability of having more than 20000 virus per monkey as a function of time . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.1 Normal dynamics of the viral load after treatment interruption 73 3.2 Observed treatment interruption dynamics of the viral load when Auranofin is also included in the treatment regimen . . 73 3.3 Graphic description of immune system model . . . . . . . . . 74 3.4 Dynamics of the model with parameters ρ = 0.049, ǫ = 0.9592, p = 20000 . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.5 Close look at viral load dynamics . . . . . . . . . . . . . . . . 81 3.6 Dynamics of the improved model with parameter q = 0.1 . . 83 3.7 Close look at viral load dynamics with improve model q=0.1 84 3.8 Dynamics of the second improve model with aM∗ . . . . . . . 85 3.9 Close look at viral load using second improved model . . . . . 86 3.10 Viral load with oscillation for special parameters . . . . . . . 87 ix Acknowledgments I would like to thank my advisor Daniel Coombs for all the help while we complete this work. Thanks to Jessica Conway and Bernhard Konrad for all the helpful discussions. Thanks to Luca Shytaj and Andrea Savarino for the data and the con- versations. Thanks to Eric Cytrynbaum for all the helpful comments and final revi- sion. x To my family. xi Chapter 1 Introduction Humman Immunodeficiency Virus (HIV) infection is considered pandemic by the World Health Organization. From its discovery in 1981 to 2009, AIDS killed nearly 30 million people [23]. However, since the advent of anti-retroviral therapy (ART) in 1996, HIV infection is no longer the death sentence it once was. ART drugs are capable of suppressing HIV to undetectable levels. How- ever, there remains a persistent viral reservoir in infected patients that pre- vents the complete eradication of HIV infection [4], [17]. The HIV reservoir is not sensitive to current therapies but there are some possible medical approaches that are under investigation. Macaques infected with simian immunodeficiency virus (SIV) represent the preferred animal model to study HIV. HIV and SIV have similar patho- genesis and genetic organization [14]. Recent studies on rhesus macaques analyze the combination of ART drugs with the drug Auranofin and the CCR5 inhibitor Maraviroc (MRV) [19] [13]. CCR5 is an essential coreceptor that HIV uses to infect immune cells. Auranofin has been shown to accelerate the activation rate of latent cells and alter the kinetics of viral rebound when drug treatment is inter- rupted [13], [1]. The obtained viral rebound was previously observed only to early days of infection [14]. Adding Auranofin to treatment also seems to decrease the resulting viral set point. On the other hand, MRV improves ART by inducing control of viremia and decreasing the viral reservoir. Such treatment improvement is called H-iART. The combination of Auranofin and H-iART substantially reduces the numbers of long lived memory T- cells restricting the viral reservoir, suppress the viral load levels completely, and leads to a persistent suppression of the viremia in absence of ART [19]. There remain observable transient episodes of viremia above the detection limit after the interruption of treatment, called viral blips. In order to understand the mechanisms of this innovative treatment, I developed mathematical models to study the post-treatment dynamics of these related immunodeficiency virus (HIV and SIV) infections. As a first step, I propose new mechanisms to analyze viral dynamics 1 Chapter 1. Introduction assuming a drug therapy that can cure the infection (assumed to be H- iART with Auranofin). I use both deterministic and stochastic approaches to investigate possible viral resurgence from the viral reservoir. The ODE models express the average dynamics but since we assume small numbers of viremia and infected cells we need more than the average dynamics. Three stochastic models are then developed. They will be assumed to be continu- ous time Markov chains. The first two models are used to study whether the activation of one infected cell can lead to the observable viral blips. From these models we found that more than one latent cell must be activated in order to have simulated blips comparable with the data. The third model includes multiple latent cell activation, with the finding that with multiple latent cells the probability of having a viral blip over the detection limit is greater than with one latent cell. It is also shown that a single latent cell would produce more than one (possible undetectable) viral blip on a path. All this is developed in Chapter 2. As a second part of the project, I develop a reasonable model that repro- duces the observed viral rebound in the presence and absence of Auranofin at treatment interruption. I believe that the faster activation of memory cells is not only clearing the latent reservoir but also returning the immune system to preinfection state. This depletion of memory immunity along with lower levels of CD4 T cells could explain the rebound observed after treatment interruption. This is described in Chapter 3. However, using a novel model, I was not able to reproduce the rebound of the viral blip after regular ART interruption. Without that characteristic in the model, I was unable to test the possible mechanism of Auranofin presence. 2 Chapter 2 Post treatment model 2.1 Introduction In this and the next chapter, we will study chronic immunodeficiency virus infection with successful ART treatment. In this chapter we will examine the data from two rhesus macaques (Macacca mulatta). Those monkeys were infected with SIVmac251. They were treated with a new treatment that included the CCR5 blocker maraviroc (MRV) as well as the memory-cell depleting drug Auranofin, in combination with antiretrovirals [19]. CCR5 is an essential coreceptor that HIV and SIVmac251 use to infect immune cells. MRV is an antagonist to this coreceptor and also inhibits SIVmac251 replication in the nanomolar concentration range. Auranofin is a central memory T cell depleting drug and we will talk about it in more detail in next chapter. The monkeys P177 and P157 presented observable viral blips after treat- ment interruption. However, the viral load doesn’t return to the viral set point after the blip. By viral set point we refer to the value of viral load when the virus has stabilized in the infected individual. Rather the mea- sured viral load returns to zero following each episode of detectable viremia. This observation suggest the possibility of a functional cure of SIV in these animals. Here, we will seek to analyze the characteristics and probabilities of the observed blips. Assuming that the activation of the remaining latent cells is driving the viral blips we begin with an ODE model. This model will include the latent pool and the classic interactions. Numerical simulation of this model will give us an general overview of the possible dynamics. However, since we are assuming a really small viral load, the stochastic fluctuation play a major role in shaping the dynamics. Therefore, I will construct a stochastic model based on the ODE model. I will start with an oversimplified stochastic model assuming the viral load is in quasi steady state. This model can be solved analytically. Then I will proceed to examine more complex models by numerical simulation. A crucial assumption underlying the stochastic models is that the system 3 2.2. Data summary of P177 and P157 monkeys is homogeneously mixed at all times in the body. Therefore in all this chapter I will modify the parameters to cover the total blood of a macaque, which according to [25] is around .5L in average for a rhesus macaque. Usually the parameters are in concentration units because the ODE model assume average numbers, while the stochastic models counts every single cell or virus in the system. 2.2 Data summary of P177 and P157 monkeys The specific drug treatments used are described below. ART Tenofovir and emtricitabine, two nucleosidic/nucleotidic reverse tran- scriptase inhibitors; and raltegravir an integrase inhibitor. Intensified ART (iART) ART with the addition of the protease inhibitor darunavir boosted with ritonovir. Highly intensified ART (H-iART) iArt reinforced with the CCR5 blocker maraviroc (MRV). Both monkeys were infected 523 days before treatment. They had an initial therapeutic cycle of 46 days of ART, followed by 80 days of iART and ending with 164 days of H-iART. [Table S1 from [19]] The P177 monkey was treated with H-iART and Auranofin for 95 days. Following therapy suspension, P177 was subjected to a short H-iART cycle for 41 days at viral rebound. −150 −100 −50 0 50 100 150 200 101 102 103 104 Days VL (lo g 1 0) VL p177 Figure 2.1: P177 viral load per ml, log10 scale −200 −100 0 100 200 300 101 102 103 104 105 Days VL (lo g 1 0) VL p157 Figure 2.2: P157 viral load per ml, log10 scale 4 2.2. Data summary of P177 and P157 monkeys Monkey P157 was treated with Auranofin and H-iART during other 100 days [18]. The plots of the data mark zero as the time of treatment interruption. The viral load is plotted on Figures 2.1 and 2.2. Observing the viral load data, the maximum high is around 104 copies/mL. The P157 monkey showed more variability on its blip’s maximum varying from 102 to more than 104 copies/mL. This monkey has the highest blip and also the lowest blip ob- served in both monkeys. On the other hand, the monkey P177 has a more ho- mogeneous maximum viral load that varies between 103 and 104 copies/mL. The CD4 counts of both monkeys are shown in Figures 2.3 and 2.4. There is not any trend observable from the data. While P177 CD4 counts oscillates with a clearly increased average after suspension of treatment, P157 monkey looks mode like a linear decreasing function. −150 −100 −50 0 50 100 150 200 400 500 600 700 800 900 1000 1100 Days CD 4 co un ts CD4 p177 Figure 2.3: P177 CD4 counts per µl −200 −100 0 100 200 300 600 800 1000 1200 1400 1600 1800 Days CD 4 co un ts CD4 p157 Figure 2.4: P157 CD4 counts per µl Plots of every of blip are on Figures 2.5 and 2.6, and a summary of the data is on Table 2.1. The maximum viral load is 13500 RNA copies/ml and goes for 40 days. Plots of the duration of the blips are shown in Figures 2.7 and 2.8. The mean time of blip duration is 26 days. Except for one blip from each monkey, the duration of the blips is around 15 to 30 days over both monkeys. On monkey P177 data, the longest one is the third blip. This one could actually be thought as two viral blips (Figure 2.5(c)), and in that case all blip’s duration are on the range of 20 to 30 days for this monkey. For macaque P157, the fist blip is the longest one. The maximum viral load obtain on each blip is plotted on Figures 2.9 and 2.10. While the highest viral load is obtained on the first blip of P157, 5 2.2. Data summary of P177 and P157 monkeys monkey P177 has the highest on the last blip reported. The average high of blip is 3914 RNA copies/ml. In both animals, there was no obvious trend in the size or duration of blips over time. Summary of viral blip data VL Days to Days from max Total data reach max to extinction duration P157 13500 12 28 40 120 6 7 13 4540 11 13 24 4000 14 7 21 160 8 8 16 1260 14 7 21 P177 2600 8 13 21 2340 6 18 24 3180 23 29 52 7440 14 14 28 Table 2.1: Summary of viral blip data. The duration of the blips are from undetected levels to reach the maximum and then to go undetected levels again 6 2.2. Data summary of P177 and P157 monkeys −10 −5 0 5 10 15 0 500 1000 1500 2000 2500 3000 Days Vi ra l l oa d (a) 45 50 55 60 65 70 0 500 1000 1500 2000 2500 Days Vi ra l l oa d (b) 90 100 110 120 130 140 150 0 500 1000 1500 2000 2500 3000 3500 Days Vi ra l l oa d (c) 160 165 170 175 180 185 190 0 1000 2000 3000 4000 5000 6000 7000 8000 Days Vi ra l l oa d (d) Figure 2.5: Zoomed-in plots of the blips observed on P177. Note the linear scale. The level of detection level at 40 cells/mL is also plotted. 7 2.2. Data summary of P177 and P157 monkeys 910 920 930 940 950 0 2000 4000 6000 8000 10000 12000 Days Vi ra l l oa d (a) 978 980 982 984 986 988 990 992 0 20 40 60 80 100 120 140 Days Vi ra l l oa d (b) 1000 1005 1010 1015 1020 1025 1030 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Days Vi ra l l oa d (c) 1075 1080 1085 1090 1095 1100 0 500 1000 1500 2000 2500 3000 3500 4000 Days Vi ra l l oa d (d) 1125 1130 1135 1140 1145 0 20 40 60 80 100 120 140 160 180 Days Vi ra l l oa d (e) 1145 1150 1155 1160 1165 1170 0 200 400 600 800 1000 1200 Days Vi ra l l oa d (f) Figure 2.6: Zoomed-in plots of the blips observed on P157. Note the linear scale. The detectable level at 40 is also plotted. 8 2.2. Data summary of P177 and P157 monkeys 1 2 3 4 20 25 30 35 40 45 50 55 Blip number D ur at io n (da ys ) Figure 2.7: Comparing the dura- tion of each viral blip on P177. 1 2 3 4 5 6 10 15 20 25 30 35 40 Blip number D ur at io n (da ys ) Figure 2.8: Comparing the dura- tion of each viral blip on P157. 1 2 3 4 2000 3000 4000 5000 6000 7000 8000 Blip number H ig h (V L) Figure 2.9: Comparing the high of each viral blip on P177. 1 2 3 4 5 6 102 103 104 105 Blip number H ig h (V L) Figure 2.10: Comparing the high of each viral blip on P157, log10 scale. 9 2.3. ODE model 2.3 ODE model The treatment described above has two extra components beyond iART: an anti-memory drug and a CCR5 blocker. This suggests that the memory pool is an important compartment in the evolution of the viral load of these macaques. In this section, I will examine simple models of HIV including latent cells. The ODE model will focus on each viral blip, not the complete post- treatment dynamics. I assume that the blips occur because a latent cell activates and that there are no other latent cells activated on the body. The ODE system is given by (2.1) below. T , T ∗ and V represent the T-cells, infected cells and virus compartments respectively. This is known as the classic HIV model. The parameters are similar to Table 3.1 but all of them are constant in this case. dT dt = λ− dT − kTV, dT ∗ dt = kTV − δT ∗ (2.1) dV dt = pT ∗ − cV. The classic model (2.1) can be generalized by adding the latent infected cell compartment L. The new system is given by (2.2) and it’s assuming that: (1) a fraction η of the infected cells goes to latency, (2) latent cells can then reproduce at rate ρ, die at rate δL, and reactivate at rate a. dT dt = λ− dT − kTV, dT ∗ dt = k(1− η)TV − δT ∗ + aL (2.2) dL dt = ηkV T − δLL− aL+ ρL, dV dt = pT ∗ − cV. The net reproductive number R0 represents the number of new infections caused by one infected cell. For the classic HIV model (2.1) R0 is given by equation (2.4) and in the case of the latent model (2.2) by equation (2.5). It has been proved that there is viral extinction if and only if R0 < 1. 10 2.3. ODE model Then I modify the parameter λ such that the extinction of the viral blip approximate the duration of the data blips as explained below. For the simulation, I used the parameters from [24] and [3] as summarized on Table 2.2 but with a change on the units. I assume the total blood volume of a macaque is .5L as an approximate value of the ranges given on [25]. Then I rescale all the parameters to mean rates on the number of cells and not of concentration of cells. In order to estimate λ, I set the mean rate of the infection r be on the timescale of 7 days. This rate is represented by the smallest eigenvalue of the linearized of [2.1]. The new λ is one order of magnitude less than the estimate in [24]. At the end, I calculate R0 for the new set of parameters. Analytic expression of λ is given by equation (2.3). λ = ( (2r + δ + c)2 − (δ − c)2 ) d 4kp (2.3) R0 = pkλ cδd (2.4) Parameter values λ 6.2707 × 104 l−1day−1 k (4.8 × 10−8)/500 l day−1 [24] d 0.011 day−1 [24] δ 1.34 day−1 [3] p 50000 day−1 [24] c 23 day−1 [24] r −17 days −1. R0 0.8878 Table 2.2: Parameters of the classic ODE model with units per .5L. Pa- rameter λ estimated as marked on the text. To solve the system I use the MATLAB function ode45 starting at the activation of one infected cell, with initial conditions T = λ d and T ∗ = 1 cells per liter, and V = 0 copies/l. The viral load is set to be zero since as shown on [19] there were less than 2 RNA copies/ml at treatment interruption. The solution shows a almost constant T cell population and a fast decreasing on the infected cells (Figures 2.11 and 2.12). On the viral load, there is an observable blip. The duration of the viral blip is approximately correct but the maximum viral count perml is too low 11 2.3. ODE model 0 10 20 30 40 50 60 106.75593 106.75593 106.75593 106.75593 Days T ce lls (lo g 1 0) Figure 2.11: T cells solution of classic model with parameters cor- responding to the complete blood per monkey. 0 10 20 30 40 50 60 10−4 10−3 10−2 10−1 100 Days T* ce lls (lo g 1 0) Monkey parameters, R0 = 0.88784 Figure 2.12: Infected cells solu- tion of classic model with parame- ters corresponding to the complete blood per monkey. 0 10 20 30 40 50 60 100 101 102 103 Days Vi ra l l oa d (lo g 1 0) Monkey parameters, R0 = 0.88784 Figure 2.13: Viral load solution of classic model with parameters cor- responding to the complete blood per monkey on logarithmic scale. 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Days Vi ra l l oa d m l Figure 2.14: Viral load solution of classic model with parameters on ml units, linear scale. 12 2.3. ODE model in comparison with the data (Figures 2.13 and 2.14). I carried out similar simulation for different values of r but there were not significant difference on the blip maximum. Now I use this approximate value of R0 to estimate λ for the latent model. I calculated R0 for the latent model using the next generation matrix approach [6]. This method can be used for models when more than one class of infectives is involved [11]. In the next generation method, R0 is defined as the spectral radius of the next generation operator. The next generation operator is obtained as following: determine F1 the matrix of rates of appearance of new infections, and V1 the matrix of the rate of incoming minus the rate of outgoing from each infectious state; then take partial derivatives and evaluate on the diseases free equilibrium (new matrix F and V ); and then the next generator operator is FV −1. In this case, matrix F and V are given by: F =  0 0 ηk λd0 0 (1− η)k λ d 0 0 0   V =  dI + a− ρ 0 0−a δ 0 0 −p c   Therefore, the analytic expression of R0 is: R0 = pkλ cδd ( 1 + η(ρ− δL) δL + a− ρ ) (2.5) The parameters used for the simulations are given on Table 2.3. The parameters for the latent population are assumed to be the same as in human and they are taken from [17]. The simulations were performed using ode45. I tested two sets of initial conditions. In both cases I used T = λ d and V = 0. One option assumes that you have one latent cell not activated yet, then L = 1 and T ∗ = 0. The second option assumes that the latent cell activated at the moment of beginning of the blip, then L = 0 and T ∗ = 1. The first set of IC shows a lower and longer viral blip than the second option as shown on Figures 2.15 and 2.16. The second option shows a solution closer to the one of the classic model. This difference highlights a problem with the use of ODE models for single-cell behavior. Essentially ODE models are appropriate only for large populations. 13 2.3. ODE model Parameter values λ 6.2647 × 104 l−1 day−1 k (4.8 × 10−8)/500 l day−1 [24] d 0.011 day−1 [24] a 0.1 day−1 [17] δ 1.34 day−1 [3] ρ 0.049 day−1 [5] η 0.001 [17] δL 0.001 day −1 [17] p 50000 day−1 [24] c 23 day−1 [24] R0 0.8878 Table 2.3: Parameters of the latent ODE model. λ estimated as described on the text. The latent cells rates are assumed to be the same as described on Rong and Perelson paper with human parameters [17]. 0 10 20 30 40 50 60 70 0 100 200 300 400 500 600 700 800 900 Days Vi ra l l oa d Monkey parameters, R0 = 0.8878 Figure 2.15: Viral load solution of latent model with initial condition T = λ d , V = 0, L = 1 and T ∗ = 0 0 10 20 30 40 50 60 70 0 500 1000 1500 2000 2500 Days Vi ra l l oa d Monkey parameters, R0 = 0.8878 Figure 2.16: Viral load solution of classic model with initial con- dition T = λ d , V = 0, L = 0 and T ∗ = 1 14 2.4. Introduction to the stochastic approach 2.4 Introduction to the stochastic approach Under the assumption of low number of cells and viruses, a stochastic ap- proach is suitable to analyze the system. I will proceed to obtain stochastic approximations of the previous latent and classic models. To start with, we assume that the T pool is constant and approximately equal to λ d . This is justified from the solution of the ODE system but this is a simplification because the experimental data shows variation in CD4 count (Figures 2.3 and 2.4). Assuming again that the activation of one latent cell is driving the observed viral blip and that there is only one latent cell in the system, the classic model simplifies as illustrated in Figure 2.17. This model now includes only infected cells and virus. I will analyze this model in detail in Section 2.6. T ∗ V cδ p kT Figure 2.17: T ∗V simplified model The two compartment model remains difficult to work with analytically. Therefore, I oversimplify the T ∗V model assuming that the important dy- namics occurs on the infected cells population. Thus we assume a quasi steady state for the viral population. This system is illustrated in Figure 2.18. This model will be called the T ∗ model. It is possible to find an analytic solution to this model: see Section 2.5. T ∗ δ pkλ cd Figure 2.18: T ∗ simplified model 2.5 T ∗ stochastic model Assuming that the most important dynamics occur at the infected cell pop- ulation, the resulting system can be model as a birth death process (Figure 15 2.5. T ∗ stochastic model 2.18). The death rate of the process is given as the death rate of the in- fected cells, since we assume that this is the only way to decrease the T ∗ cell population. The birth rate is a bit more complicated. A new infection occurs when a susceptible cell encounters a virus with rate k. Then the rate of infection is given by k times the mean number of susceptible cells λ c times the mean number of virus. From the quasi steady state of the T ∗V system, the mean number of virus is assumed to be the ratio between the new produced virus and the clearance, p c . Then the birth date is given by pkλ cd . We need to make a few technical (and natural) assumptions regarding the system. It is assumed to be a stationary Markov chain, i.e. that the future is independent of the past, given the present. Also, it is homogeneous in time, i.e. P (T ∗(t) = n|T ∗(r) = ñ) = P (T ∗(t − r) = n|T ∗(0) = ñ). This last assumption implies that we only need to calculate the probability distribution assuming that the process starts at time 0. The possible events in the system are listed as follows: An uninfected T cell becomes infected ∅ pkλ cd =ξ −−−−→ T ∗, An infected T cell dies T ∗ δ −→ ∅. The probability of having n infected cells at time t1 given that we start with ñ at time t2 is called the transition probability. Simplifying the nota- tion, let Pñ,n(t) = P (T ∗(t) = n|T ∗(0) = ñ). Now assume that h is small enough such that the probability of two or more events happening in the time interval (0, h) is in o(h). Then there are only three possibilities on (0, h): (1) one infected cell dies, (2) one infected cell is born, or (3) nothing happens. The rate that one cell dies is δ, then the rate that one of a total of n + 1 cells dies is δ(n + 1). Therefore the probability of going from n+ 1 infected cells at time zero to n at time h is given by (n + 1)δh + o(h), which is exactly the definition of the transition probabilities. The calculation of the other possible changes are analogous. Then we have that the transition probabilities on that interval are given by: Pn+1,n(h) =(n+ 1)δh + o(h) Pn−1,n(h) =(n− 1)ξh+ o(h) Pn,n(h) =1− (nδ + nξ)h+ o(h). 16 2.5. T ∗ stochastic model 2.5.1 Kolmogorov’s forward differential equations Let T ∗(0) = ñ be the given initial number of infected cells. In order to find the distribution of T ∗(t+ h) at a later time t+ h, we split the time interval (0, t + h) into the long time step (0, t) and the short time step (t, t + h). Then we use Kolmogorov’s forward equation: Pñ,n(t+ h) = ∞∑ r=0 Pñ,r(t)Pr,n(h). to obtain a relationship on the transition probabilities. Thus, it follows that: Pñ,n(t+ h) = Pñ,n+1(t)Pn+1,n(h) + Pñ,n−1(t)Pn−1,n(h)+ Pñ,n(t)Pn,n(h) + o(h) = (n+ 1)δhPñ,n+1(t) + (n− 1)ξhPñ,n−1(t) + (1− (nδ + nξ)h)Pñ,n(t)) + o(h) Then dividing by h and taking the limit when h goes to zero we get an equation for the derivative of the transition probabilities called the Kol- mogorov’s Forward differential equation (KFDE). dPñ,n(t) dt = δ((n + 1)Pñ,n+1(t)− nPñ,n(t)) + ξ((n − 1)Pñ,n−1(t)− nPñ,n(t)) I.C. Pñ,n(0) = δñ,n Since the initial condition does not change, we can simplify the notation to be Pñ,k(t) = pk(t) where k is the ending point. Then the forward equation remains as in (2.6). dpn(t) dt = δ((n + 1)pn+1(t)− npn(t)) + ξ((n − 1)pn−1(t)− npn(t)) (2.6) I.C. pn(0) = δñ,n 17 2.5. T ∗ stochastic model Now notice that the mean of the process satisfies: dE(T ∗(t)|T ∗(0) = ñ) dt = ∞∑ n=0 n ( dPñ,n dt ) = ∞∑ n=0 nδ((n + 1)Pñ,n+1(t)− nPñ,n(t)) (2.7) + ∞∑ n=0 nξ((n− 1)Pñ,n−1(t)− nPñ,n(t)) = ∞∑ n=0 ξn(n+ 1− n)pñ,n + ∞∑ n=0 δn(n− 1− n)pñ,n = (ξ − δ)E(T ∗(t)|T ∗(0) = ñ) Then the mean is the solution of the ordinary differential equation rep- resenting an exponential growth or decay of the infected population [7]. dT ∗ dt = (ξ − δ)T ∗ (2.8) 2.5.2 Kolmogorov’s backward differential equations As in the previous section we will find a differential equation for the tran- sition probabilities. This time we will split the time interval (0, t + h) into the short time step (0, h) and the long time step (h, t+ h). Then: Pñ,n(t+ h) = ∞∑ r=0 Pñ,r(h)Pr,n(t). = Pñ,ñ−1(h)Pñ−1,n(t) + Pñ,ñ+1(h)Pñ+1,n(t) + Pñ,ñ(h)Pñ,n(t) + o(h) = ñδhPñ−1,n(t) + ñξhPñ+1,n(t) + (1− (ñδ + ñξ)h)Pñ,n(t)) + o(h) Again dividing by h and taking the limit as h goes to zero we get a differential equation for the transition probabilities but this time depending on the initial condition. It is called the Kolmogorov’s backward differential equation (KBDE). Again we simplify the notation since now the ending time is not changing. Then equation (2.9) is the backward differential equation where Pk,n(t) = pk(t) for k the starting point. 18 2.5. T ∗ stochastic model dpñ(t) dt = δñ(pñ−1(t)− pñ(t)) + ξñ(pñ+1(t)− pñ(t)) (2.9) I.C. pñ(0) = δñ,n 2.5.3 Probability generating function (pgf) definition We would gain a lot of insight into our model if we could solve Kolmogorov’s forward or backward differential equation. Unfortunately, this is not so easy as both formulations describe a system of infinitely many linear ODEs. Fur- thermore, we do not have a natural first equation we can solve to proceed inductively from there, as the equations are highly interconnected. A popu- lar way to still make progress is to proceed using the probability generating function defined below. The multi-variable function G := G(x; t) := Gñ(x; t), defined as Gñ(x; t) = E[x T ∗ ] = ∞∑ n=0 Pñ,n(t)x n, is called probability generating function (pgf). Note that for x ∈ D, where D = {z ∈ C : |z| ≤ 1}, and each time t ∈ [0,∞) the right hand side converges by the triangle inequality, since probabilities are non-negative and ∑ ∞ n=0 Pñ,n(t) = 1 . In particular, we have that G(1; t) = 1. The probability generating function has many useful properties: 1. We can recover the probability of having n infected T cells at time t by taking the derivatives, P (T ∗(t) = n|T ∗(0) = ñ) = 1 n! dnG(x; t) dxn ∣∣∣∣ x=0 . In particular, the probability of the infected cells to be extinct at time t is given by G(0; t). 2. We can calculate the mean value of the number of infected cells at time t by evaluating E[T ∗] := dG(x; t) dx ∣∣∣∣ x=1 respectively. The calculation of higher order moments of the distribution can be calculated also using higher order derivatives of the pgf. 19 2.5. T ∗ stochastic model All these properties are easily generalized to higher dimension random variables. 2.5.4 Applying the probability generating function to the backward equation (2.9) In order to obtain an equation for the pgf we will multiply equation (2.9) by xn and sum over all the possibles values of n, ∞∑ n=0 dpñ(t) dt xn = ∞∑ n=0 δñ(pñ−1(t)− pñ(t))x n + ∞∑ n=0 ξñ(pñ+1(t)− pñ(t)))x n Then using the definition of the pgf and expanding over the sums we get: dGñ(t) dt = δñ(Gñ−1(t)−Gñ(t)) + ξñ(Gñ+1(t)−Gñ(t)) (2.10) with I.C. Gñ(0) = x ñ Furthermore we can assume that each infected cell acts as an independent individual. In that case, the pgf has a branching property Gñ(t) = (G1(t)) ñ [21] and then the system of equations given by (2.10) reduces to (2.11). dG1(t) dt = δ(G0(t)−G1(t)) + ξ(G2(t)−G1(t)) dG1(t) dt = δ(1 −G1(t)) + ξ((G1(t)) 2 −G1(t)) dG1(t) dt = (ξG1(t)− δ)(G1(t)− 1) (2.11) with I.C. G1(x; 0) = x Now in order to solve the (2.11) we make a change of notation y = G1 and proceed by separation of variables. 20 2.5. T ∗ stochastic model y′ = (ξy − δ)(y − 1) ⇒ 1 = y′ (ξy − δ)(y − 1) ⇒ ∫ 1dt = ∫ 1 (ξy − δ)(y − 1) dy = ∫ 1 (ξ − δ)(y − 1) dy + ∫ −ξ (ξ − δ)(ξy − δ) dy ⇒ t+A = 1 ξ − δ ln(y − 1)− 1 ξ − δ ln(ξy − δ) = 1 ξ − δ ln ( y − 1 ξy − δ ) Then y = 1−Aδe(ξ−δ)t 1−Aξe(ξ−δ)t with I.C. y(0) = x ⇒ A = x− 1 ξx− δ Therefore the pgf for the stochastic process assuming independent bran- ching with initial condition one cell is given by (2.12). G1(x; t) = ξx− δ − (x− 1)δe(ξ−δ)t ξx− δ − (x− 1)ξe(ξ−δ)t (2.12) The general solution for any initial condition can be obtained by applying the branching property to equation (2.12). 2.5.5 The probability generating function applied to forward equation (2.6) Now we do a similar approach using the forward differential equation given by (2.6). Multiplying by xn and summing over all the possibles values of n 21 2.5. T ∗ stochastic model we get: ∞∑ n=0 dpn(t) dt xn = ∞∑ n=0 δ(n + 1)pn+1(t)x n − ∞∑ n=0 δnpn(t)x n + ∞∑ n=0 ξ(n− 1)pn−1(t)x n − ∞∑ n=0 ξnpñ(t)))x n ⇒ ∂Gñ ∂t = δ ∂Gñ ∂x − δx ∂Gñ ∂x + ξx2 ∂Gñ ∂x − ξx ∂Gñ ∂x In order to find the pgf, it will be needed to solve the partial differential equation given by (2.13). ∂Gñ ∂t = (ξx− δ)(x − 1) ∂Gñ ∂x (2.13) with I.C. Gñ(x; 0) = x ñ The solution of this pde is obtained using the method of characteristics [22]. The characteristic equation is give by (2.14). dx dt = −(ξx− δ)(x − 1) (2.14) with I.C. x(0) = x0 t > 0 The characteristic equation is exactly the negative of the equation that we solve for the backward case. The solution is give on (2.15). x(t) = ξx0 − δ − (x0 − 1)δe −(ξ−δ)t ξx0 − δ − (x0 − 1)ξe−(ξ−δ)t (2.15) Now, since we assume that the solution over the characteristics is con- stants on time we just need to look for x0. Thus the solution of the pgf coincides with the solution of the backward equation. Gñ(x(t); t) = G(x0; 0) = x ñ 0 ∴ Gñ(x; t) = ( δ(x− 1)− (ξx− δ)e−(ξ−δ)t ξ(x− 1)− (ξx− δ)e−(ξ−δ)t )ñ (2.16) 22 2.5. T ∗ stochastic model 2.5.6 Transition probabilities From the pgf we can calculate the transition probabilities defined before as P (T ∗(t) = n|T ∗(0) = ñ). They are given by the derivatives of the pgf over a factorial number as introduced previously. I compute the exact solution for the case when ñ = 1, ie, when there is only one infected cell to start with. Since cells are independent, in the model, this case gives us a complete picture of the dynamics. P (T ∗(t) = n|T ∗(0) = 1) = 1 n! dnG1(x; t) dxn ∣∣∣∣ x=0 = (−bc+ ad)cn−1(−1)n−1 dn+1 (2.17) For n ≥ 1, where a = δ − ξe−(ξ−δ)t b = −δ + δe−(ξ−δ)t c = ξ − ξe−(ξ−δ)t d = −ξ + δe−(ξ−δ)t Notice that this formula does not apply for the case when n = 0. For that case see subsection 2.5.7. Since this is a probability it has to be positive. The factor (−1)n makes us wish to go further and be sure it is positive. To prove that we proceed by two cases. Case 1: ξ < δ Since ξ < δ then e−(ξ−δ) > 0 ∀t > 0. Thus, b = −δ + δe−(ξ−δ)t > 0 c = ξ − ξe−(ξ−δ)t < 0 d = −ξ + δe−(ξ−δ)t > 0 And, ad = −ξδ + δ2e−(ξ−δ)t + ξ2e−(ξ−δ)t − ξδ(e−(ξ−δ)t)2 bc = −ξδ(e−(ξ−δ)t)2 − ξδ + 2ξδe−(ξ−δ)t ⇒− cb+ ad > 0 ∴ P (T ∗(t) = n|T ∗(0) = 1) > 0 23 2.5. T ∗ stochastic model Case 2: ξ > δ This case is analogous to the previous one. Using the parameters listed in Table 2.2, I calculated examples of the transition probabilities with an initial condition of one infected cell. Figure 2.19 represents the probability distribution function of the process at day 7. It seems to follow a power law decay. Figure 2.20 illustrates the probability of having 5 infected cells at different times. It is observable an exponential decay as time goes on. The observations on decay from these plots are formally proven in subsection 2.5.8. Figure 2.19: Probability distribu- tion of model T ∗ at time 7 days Figure 2.20: Probability to have 5 infected cells as function of time In general, the derivative of the pgf is given below (2.18). This derivative can be used to calculate the expectation and variance. G (n) 1 (x; t) = (−bc+ ad)cn−1(−1)n−1n! (cx+ d)n+1 (2.18) (For n ≥ 1 and a, b, c and d as before.) 24 2.5. T ∗ stochastic model 2.5.7 Extinction probability From the pgf we easily can calculate the probability that the process becomes extinct at time t. P (T ∗(t) = 0|T ∗(0) = ñ) = (G1(0; t)) ñ = ( −δ + δe−(ξ−δ)t −ξ + δe−(ξ−δ)t )ñ (2.19) The extinction probability as a function of time is an increasing function since zero is an absorbing state (once you get there you never leave). This is illustrated when there is only one infected cell at time zero in Figure 2.21 using parameters from Table 2.2. Notice that the probability of extinction by day 10 is approximately 0.969. Figure 2.21: Extinction probabil- ity function starting with T ∗(0) = 1. Notice that the probability of extinction is greater than .9 after 10 days. Figure 2.22: Probability that at day 7 there are no infected cells as a function of the initial condi- tions. Using the same parameters I plotted the probability of extinction at time 7 for different initial conditions (Figure 2.22). This probability is decreasing. With more infected cells at the beginning, the probability of extinction at one week is smaller. Starting with 15 infected cells will give us less than .5 probability of extinction after a week. Starting with more than one infected cell could be interpreted as the simultaneous activation of latent cells. 25 2.5. T ∗ stochastic model The extinction probability is given by taking the limit of (2.19) when t goes to infinity. That limit depends on the parameters. We calculate the two possible cases. lim t→∞ P (T ∗(t) = 0|T ∗(0) = ñ) =   1 iff ξ ≤ δ( δ ξ )ñ iff ξ > δ. (2.20) 2.5.8 Limit analysis of the transition probabilities We now formally look at the behavior of the transition probabilities for the case where the extinction is guaranteed. This is when ξ < δ. How does the probability of having T ∗(t) = n, at a fixed time, decay as you increase n? We have P (T ∗(t) = n|T ∗(0) = 1) = (−bc+ ad)cn−1(−1)n−1 dn+1 Then for n ≥ 1, where b > 0 c < 0 d > 0 − cb+ ad > 0 and ξ < δ ⇒ − c d < 1 we find P (T ∗(t) = n|T ∗(0) = 1) = (−bc+ ad) d2 ( − c d )n−1 Therefore the probability decays as a power function of n− 1. At a fixed number of T ∗, how does the probability P (T ∗(t) = n|T ∗(0) = 1) decay as a function of time? Here I need another factorization. Remember from previous calculations that −bc+ ad = (ξ − δ)2e−(ξ−δ)t, thus: 26 2.5. T ∗ stochastic model P (T ∗(t) = n|T ∗(0) = 1) = (ξ − δ)2e−(ξ−δ)t(ξ − ξe−(ξ−δ)t)n−1(−1)n−1 (−ξ + δe−(ξ−δ)t)n+1 = (ξ − δ)2(−ξ)n−1 ( 1− e−(ξ−δ)t −ξ + δe−(ξ−δ)t )n−1 × ( e−(ξ−δ)t −ξ + δe−(ξ−δ)t ) 1 −ξ + δe−(ξ−δ)t Using L’Hopital’s rule it is easy to prove that when t goes to infinity: ( 1− e−(ξ−δ)t −ξ + δe−(ξ−δ)t )n−1 → ( − 1 δ )n−1 e−(ξ−δ)t −ξ + δe−(ξ−δ)t → 1 δ Therefore the probability decays over time as 1 −ξ+δe−(ξ−δ)t . Notice that this does not depend on the number of T ∗. Finally, we calculate the probability of exceeding a particular number of T ∗ at any point during the blip. P (T ∗(t) > n|T ∗(0) = 1) = 1− P (T ∗(t) ≤ n|T ∗(0) = 1) = 1− n∑ i=0 P (T ∗(t) = i|T ∗(0) = 1) = 1− b d − n∑ i=1 (−bc+ ad)ci−1(−1)i−1 di+1 = 1− b d − (−bc+ ad) d2 ( 1− ( − c d )n 1− ( − c d ) ) = 1− b d − (−bc+ ad)(dn − (−c)n) (d+ c)dn+1 27 2.5. T ∗ stochastic model = 1 + dn+1(−b− a) + (−bc+ ad)(−c)n (d+ c)dn+1 = 1− b+ a d+ c + ( −c d+ c ) P (T ∗(t) = n|T ∗(0) = 1) = ξ(e−(ξ−δ)t − 1) (−ξ + δ)e−(ξ−δ)t P (T ∗(t) = n|T ∗(0) = 1) 2.5.9 First passage time to extinction Now let’s find a relationship between the transition probabilities and the time at which the number of cells reaches a certain value for the first time. We will call this the first passage time (FPT). The generalization of the next method to any initial condition is simple but I will do the calculation starting with one infected cell first. Let Fn be the random variable representing the FPT to have n infected cells. Notice that the probability of having n infected cells at time t is the same as the probability of first have n infected cells at some time t1 < t times the probability of going back to n at time t. This relation leads directly to a solution if we assume that time is discrete but in this case the FPT is a continuous random variable so we will work with the probability density function. The equation (2.21) represents this relationship where fFn is the density function of the FPT at n: P (T ∗(t) = n|T ∗(0) = 1) = ∫ t 0 fFn(t1)P (T ∗(t) = n|T ∗(t1) = n)dt1 (2.21) = ∫ t 0 fFn(t1)P (T ∗(t− t1) = n|T ∗(0) = n)dt1 The case of the first time of extinction is when n = 0. And since zero is a stationary state we can calculate explicitly the cumulative distribution function and the density function of the F (0). P (T ∗(t) = 0|T ∗(0) = 1) = ∫ t 0 fF0(t1)P (T ∗(t) = 0|T ∗(t1) = 0)dt1 = ∫ t 0 fF0(t1)dt1 = P (F0 ≤ t) 28 2.5. T ∗ stochastic model From the previous section, we know that P (F0 ≤ t) = −δ + δe−(ξ−δ)t −ξ + δe−(ξ−δ)t (2.22) so fF0(t) = dP (F0 ≤ t) dt (2.23) = δ(ξ − δ)2e−(ξ−δ)t (−ξ + δe−(ξ−δ)t)2 . The cumulative function gives the probability of having zero infected cells at time t and the density function gives the probability of having one infected cell at time t times δ, i.e, the probability of having one infected cell times the probability that it dies. The cumulative function is illustrated on Figure 2.21 and the density function is represented on Figure 2.23. Figure 2.23: Probability density function of first time to extinction Since we are interested in the case when the process goes extinct with probability one, we take ξ − δ < 0. We can integrate the density function 29 2.5. T ∗ stochastic model fF0 from zero to infinity to check we are in the correct track:∫ ∞ 0 fF0(t)dt = ∫ ∞ 0 δ(ξ − δ)2e−(ξ−δ)t (−ξ + δe−(ξ−δ)t)2 dt = δ(ξ − δ) ∫ −∞ 0 e−u (−ξ + δe−u)2 du = −δ(ξ − δ) ∫ ∞ 1 dv (−ξ + δv)2 = −(ξ − δ) ∫ ∞ −ξ+δ dw w2 = 1 We now calculate the mean passage time from the definition. E[F0] = ∫ ∞ 0 tfF0(t)dt = ∫ ∞ 0 t δ(ξ − δ)2e−(ξ−δ)t (−ξ + δe−(ξ−δ)t)2 dt = −δ ∫ 0 −∞ ue−u (−ξ + δe−u)2 du = u 1 −ξ + δe−u ∣∣∣∣ 0 −∞ − ∫ 0 −∞ 1 −ξ + δe−u du = 1 ξ ln ( δ δ − ξ ) (2.24) Using the definition of R0 = pkλ cδd as in the classic ODE model, the mean first time to extinction is related to it in the following way, E[F0] = 1 ξ ln ( δ δ − ξ ) = − 1 ξ ln(1−R0) The generalization of the initial condition is given as follows. Define Fñ,0 as the first passage time to extinction given that we start with ñ infected cells and fFñ,0 the corresponding density function. Analogous to the previous case, we have: 30 2.5. T ∗ stochastic model P (T ∗(t) = 0|T ∗(0) = ñ) = ∫ t 0 fF0(t1)P (T ∗(t) = 0|T ∗(t1) = 0)dt1 = ∫ t 0 fF0(t1)dt1 = P (F0 ≤ t) We also know from previous section that P (F0 ≤ t) = ( −δ + δe−(ξ−δ)t −ξ + δe−(ξ−δ)t )ñ (2.25) so fF0(t) = dP (F0 ≤ t) dt (2.26) = ñ ( −δ + δe−(ξ−δ)t −ξ + δe−(ξ−δ)t )ñ−1( δ(ξ − δ)2e−(ξ−δ)t (−ξ + δe−(ξ−δ)t)2 ) . 2.5.10 Numerical simulations There is an exact algorithm to simulate a path of the process known as Gillespie’s algorithm. Gillespie’s Algorithm became popular in 1977 when Gillespie published what became a landmark paper for stochastic simula- tions. [8]. Gillespie himself was motivated by chemical reactions of spatially uniformly mixed species in a fixed volume. Therefore a possible reaction is sometimes referred to as reaction channel. The algorithm works as follows. At time t calculate the probability of each possible reaction. The probabilities depend on the number of individu- als ai in each compartment i = 1, 2, · · · , and on the total number of individ- uals a0 := ∑ i ai. To decide which reaction occurs next, generate a random number r1 and find the parameter µ such that ∑µ−1 i=1 ai < r1a0 ≤ ∑µ i=1 ai. Interpret the parameter µ as the number of the next reaction and update the compartments accordingly. Next, generate a number r2 which is expo- nentially distributed with updated parameter a0, and interpret this as the time when the next reaction occurs. Repeat this procedure at t = t+ r2. We now use this algorithm to simulate the T ∗ model, using the same parameters as in the ODE model (Table 2.2) where δ = 1.34 day−1, ξ = 1.18 day−1 and R0 = 0.8878. The simulation starts with one infected cell and 31 2.5. T ∗ stochastic model 0 2 4 6 8 10 0 5 10 15 20 25 30 35 40 Days T* p er .5 L Figure 2.24: Different paths from the T ∗ model Gillespie simulation. δ = 1.34 day−1, ξ = 1.18 day−1 and R0 = 0.8878. 0 2 4 6 8 10 0 10 20 30 40 50 60 Days T* p er .5 L Figure 2.25: Maximum number of infected cells at different times continues until the population is extinct. This process was repeated for 10000 trajectories. All the paths were extinct after 80 days and only 3% of the trajectories went for more than 10 days. Some of the resulting paths are shown on Figure 2.24. The majority of the paths go extinct really quickly in line with the FPT calculation (above). In order to compare with the maximum blip from the data, the maximum of the 10000 simulations at each 0.01 time is plotted on Figure 2.25. 0 2 4 6 8 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Days T* p er .5 L Simulation Analytic Figure 2.26: Mean of the paths and the analytic solution of the ex- pectation 0 2 4 6 8 10 0 0.5 1 1.5 2 2.5 3 3.5 4 Days T* p er .5 L Figure 2.27: Variance of the sim- ulations 32 2.5. T ∗ stochastic model We also compute the mean and the variance of the paths. The mean follows a slow exponential decay as expected from analysis of the analogous continuous problem (Figure 2.26). The variance increases as the time in- creases (Figure 2.27). This is expected since a lot of the paths go extinct fast but a few of them go higher and longer. 0 2 4 6 8 10 0 20 40 60 80 100 120 140 Days # of e xt s im ul at io ns Mean of simulations = 1.7833 Mean FPT = 1.8389 Figure 2.28: First time to extinction from simulated data We can also calculate the first time to extinction plotted on Figure 2.28. The mean of the simulation is 1.7833 that is comparable with the analytic solution for those parameters, 1.839. At day 10, 3% of the paths were not extinct which is also in agreement with the analytic calculation of 3.1%. 2.5.11 Comparing results with data of P157 and P177 monkeys So far, the numerical simulations have been in agreement with the analytic solutions. Now I will compare the model with the data using both the simulation and the calculations. The first thing to do is change the data to infected cells. The data is viral load RNA per ml. Assuming that a T infected cell produce virus at rate p and the virus dies at rate c, each T cells produce in average p c virus, i.e. V ≈ p c T ∗. Then we can estimate the maximum viral load per blip with an approximate number of infected cells. This is summarized on Table 2.4 where the T ∗ conversion is already on a per monkey basis, i.e. per .5L of blood. From the simulation the maximum number of infected cells is around 50 infected cells per monkey. Comparing that with Table 2.4, it’s clear that only two of the experimentally measured viral blips are consistent with the simulated viral blips; the rest are much larger than the simulations predict. 33 2.5. T ∗ stochastic model Summary of viral blip data Max VL T ∗ Days to Days from max Total data converted reach max to extinction duration P157 13500 1552.5 12 28 40 120 13.8 6 7 13 4540 522.1 11 13 24 4000 460 14 7 21 160 18.4 8 8 16 1260 144.9 14 7 21 P177 2600 299 8 13 21 2340 269.1 6 18 24 3180 365.7 23 29 52 7440 855.6 14 14 28 Table 2.4: Conversion of the maximum of the data to be used with the T ∗ model. The Viral load is change to T infected cells using the quasi steady state approximation. The duration of the blips are from undetected levels to reach the maximum and then to go undetected levels again Therefore the two possible observed viral blips have a probability greater than .0001 of happen with this choise of parameters. In order to compare the probabilities of the duration of the blips, the blip was divided in to time intervals. The first interval is taken from the last undetectable level to the time t1 where the maximum high of the blip appear. Table 2.4 shows this as days to reach max. The second interval is from the maximum high of the blip until the first undetectable level measurement. Notice that then we have to calculate probabilities of the form P (T ∗(t1) = N |T ∗0 = 1) and P (T ∗(t2) = 0|T ∗ 0 = N), for which we already have analytic expressions. The results of those probabilities are on Table 2.5. It is clear from Table 2.5 that only two of the observed blips are even remotely probable in a set of 10000 blips. These is in agreement with the observation we already made about the maximum of the Gillespie simula- tions. The probabilities of higher blips are small starting with one infected cells points. Then, there are other to two possibilities to get the blipsfrom this model: the number of latent cells activated to get the observable blip is 34 2.5. T ∗ stochastic model Table of probabilities of blips N t1 t2 P (T ∗(t1) = N |T ∗ 0 = 1) P (T ∗(t2) = 0|T ∗ 0 = N) P157 1553 12 28 3.35E-098 0.07 14 6 7 0.0010 0.44 522 11 13 4.97E-036 6.91E-005 460 14 7 3.09E-030 2.11E-012 18 8 8 0.0004 0.43 145 14 7 7.57E-012 0.0002 P177 299 8 13 2.30E-024 0.0041 269 6 18 4.55E-025 0.12 366 23 29 1.54E-023 0.59 856 14 14 2.36E-053 1.79E-006 Table 2.5: Probabilities of having N infected cells after t1 days post infection starting with one infected cell at time zero and probabilities of going extinct after t2 days starting with N infected cells. bigger than one; or the R0 is even closer to one. I will analyze the activation of multiple latent cells first and then modify the parameters to have an R0 closer to one. If we consider the independent activation of the latent cells, then it can be seen as multiple independent processes each starting with one infected cell, which is exactly the case that has been described so far in this section. From the Gillespie simulation the maximum number of infected cells obtained was around 50. Then if we set the maximum of a blip starting with one infected cell to be approximately 50 infected cells we can find a lower bound on the number of latent cells initially needed to get the observed blips. This will be equivalent to dividing N from Table 2.5 by 50. In that case, we will need the activation of at least 3 to 31 latent cells at the same time to observe the data blips (3 for the shortest and 31 for the largest data blip). The average number of latent cells need is 10. Now this is a lower bound assuming independent latent cell activating at the same time, but this event seems to be rare. Having more than one latent cell could be better modeled assuming the latent cell divides into one or more latent cells and activate. In that case, the latent cells are not independent and then a new way of modelling is needed. In section 2.7 I 35 2.5. T ∗ stochastic model analyze the stochastic version of the latent ODE model, where the division of the latent cells is also included. The second idea to reconcile the model with the data is to change the parameters so that R0 is closer to 1. The idea used to modify the parameters was to approximate the mean first passage time to extinction (equation (2.24)), MFPT, to be approximate the average time of data blip duration, which is 24 days. While I was trying to calculate this I realize that the precision on the parameters has to be more than 10−10 to get a MFPT = 24. I decide to approximate to a big enough MFPT within a 10−4 accuracy on the param- eters. The parameters that changed were p and c, the production and the clearance of virus. p =50072.6059 c =20.45 The production of virus increases and the clearance decrease so that in average the birth rate of infected cells increases. For these parameters, MFPT = 16.94 and R0 = 0.99999999986. I run 10000 simulations of the model with the new parameters. Again only 3% of the paths go for more than 10 days but now some of them go for more than 500 days. The longest one goes for 16798 days. This is showing how now the tale of the distribution is heavier for longer times. Examples of the paths for the first 20 days are showed on Figure 2.29 and the longer 5 paths are plotted for 5000 days on Figure 2.30. 0 5 10 15 20 0 5 10 15 20 25 30 35 40 Days T* p er .5 L Figure 2.29: Different paths from the T ∗ model numerical simula- tion, new parameters 0 1000 2000 3000 4000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Days T* p er .5 L Figure 2.30: The longest 5 paths of the T ∗ model with new param- eters 36 2.5. T ∗ stochastic model For these parameters the mean and variance are bigger than for previ- ous ones since now the exponential decay is really slow and the tail of the distribution is heavier. This is observable on Figures 2.31 and 2.32. 0 5 10 15 20 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 Days T* p er .5 L Figure 2.31: Mean of the paths with new parameters 0 5 10 15 20 0 5 10 15 20 25 30 35 40 45 Days T* p er .5 L Figure 2.32: Variance of the sim- ulations with new parameters The maximum values go really high for some cases, but every path that goes extinct in no more than 1000 days has no more than 1500 infected cells. For longer times, the paths reach values over 5000 and the maximum of all simulations was around 8300 infected cells. Figures 2.33 and 2.34 show the maximum amount of infected cells not including the longest paths and including them. Notice that the longest path are only 3 paths on 10000 simulations but it makes a huge difference on the maximum values. 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 1400 Days T* p er .5 L Figure 2.33: Maximum with new parameters but excluding longest paths 0 200 400 600 800 1000 0 1000 2000 3000 4000 5000 6000 Days T* p er .5 L Figure 2.34: Maximum number of infected cells with new parameters When I calculate the mean first passage to extinction over the simulated 37 2.5. T ∗ stochastic model data, the result was different from the analytic result. This is because this mean value is high sensitive to the accuracy of the parameters and the computer error is affecting the approximation. Also there will be needed more simulation with these parameters given the heavy of the distribution tail. Figure 2.35 shows the first time to extinction of the data for the first 20 days from numerical simulations. The mean value obtained numerically was 6.8247. 0 5 10 15 20 0 20 40 60 80 100 120 140 Days # of e xt s im ul at io ns Mean of simulations = 6.8247 Mean FPT = 16.94 Figure 2.35: First time of extinction from numerical simulations Now there a really high blips that can be compared with the data. But these blips are long blips and the blips in the data do not last more than 52 days. The maximum value for the first 30 days is around 300 infected cells. Then it could be possible to observe 5 blips similar to the data with a probability greater than 0.0001. This statement can not be completely proved since the duration of the blips is now longer than the observed data. On the other hand, using the analytic solutions I can compute the prob- abilities as I did for the previous parameters. With these parameters the probabilities decrease a lot with respect to the previous ones but still at really low levels, (Tables 2.5 and 2.6 respectively). 2.5.12 Discussion of the model The T ∗ model presents an insight into possible dynamics of the infection. It was a good model to start with because there was possible to find analytic expressions for the transition probabilities and FPT. This model helped to develop an intuition on how much the infected cells drives the infection. This model also proves that the simple dynamics of infected cells is not enough to explain the observed post treatment blips. It was able to 38 2.5. T ∗ stochastic model Table of probabilities of blips with new parameters N t1 t2 P (T ∗(t1) = N |T ∗ 0 = 1) P (T ∗(t2) = 0|T ∗ 0 = N) P157 1553 12 28 2.85E-039 1.81E-016 14 6 7 0.0034 0.2965 522 11 13 2.55E-016 5.65E-012 460 14 7 1.69E-012 1.12E-018 18 8 8 0.0019 0.24 145 14 7 3.32E-006 2.11E-006 P177 299 8 13 4.33E-013 3.77E-007 269 6 18 9.34E-015 6.07E-005 366 23 29 3.17E-008 0.0003 856 14 14 1.95E-020 7.23E-018 Table 2.6: Probabilities of having N infected cells after t1 days post infection starting with one infected cell at time zero and probabilities of going extinct after t2 days starting with N infected cell using the new parameters. reproduce only the small data blips from monkey P157 maybe because there is more than one latent cell producing those blips. It was also observed that the model is sensitive to the parameters when R0 is really close to one. This is observed in the differences between the Gillespie simulations and the analytic results. In the case when R0 is not close to one, the results were consistent. From this model we can conclude that smaller viral blips can be observed in the timescale of days and that bigger blips are also possible but they will be detectables for longer times, even years for certain parameters. This is different than observations, where there are higher blips in short timescale. The natural extension of this model is to include the viral dynamics of the disease and that is described and simulated in the next section. 39 2.6. Two-compartment stochastic model 2.6 Two-compartment stochastic model In the previous section we assumed that the viral dynamics were controlled by the infected cells, or equivalently, that the viral load was in quasi-steady- state with the infected cells. In this section, that assumption is removed, leading to a two dimensional branching process. The diagram of the process is given in Figure 2.18. The possible reactions are: An infected cell dies T ∗ δ −→ ∅ A new virus particle is produced T ∗ p −→ T ∗ + V An uninfected T cell becomes infected V kλ d−−→ T ∗ A virus particle is cleared V c −→ ∅ Analogous to the previous section, this process will also be assumed to be an stationary continuous time Markov chain with discrete states. Define the transition probabilities to be Pñ,ṽ;n,v(t) = P (T ∗(t) = n, V (t) = v |T ∗(0) = ñ, V (0) = ṽ), In this and the next section we will need the concept of the marginal distributions. The marginal distributions are defined as: Pñ,ṽ;n = P (T ∗(t) = n|T ∗(0) = ñ, V (0) = ṽ) = ∞∑ v=0 Pñ,ṽ;n,v(t), Pñ,ṽ;v = P (V (t) = v|T ∗(0) = ñ, V (0) = ṽ) = ∞∑ n=0 Pñ,ṽ;n,v(t). Also assuming that the probability of two or more events more events in the time interval (0, h) is in o(h), the following postulates hold: Pn+1,v;n,v(h) = δ(n + 1)h+ o(h), Pn,v−1;n,v(h) = pnh+ o(h), Pn−1,v+1;n,v(h) = ( kλ(v + 1) d ) h+ o(h), Pn,v+1;n,v(h) = c(v + 1)h+ o(h). The next few subsections will be following similar methods as in previous section. The Chapman-Kolmogorov equality will be used to find backward and forward equations of the transition probabilities. Then we will derive equations for the probability generating function. 40 2.6. Two-compartment stochastic model 2.6.1 Kolmogorov’s forward differential equations Let T ∗(0) = ñ and V (0) = ṽ be the initial condition of the process. Then, using the same arguments as in Subsection 2.5.1, we derive that Pñ,ṽ;n,v(t+ h) = Pñ,ṽ;n+1,v(t)Pn+1,v;n,v(h) + Pñ,ṽ;n,v−1(t)Pn,v−1;n,v(h) + Pñ,ṽ;n−1,v+1(t)Pn−1,v+1;n,v(h) + Pñ,ṽ;n,v+1(t)Pn,v+1;n,v(h) + Pñ,ṽ;n,v(t)Pn,v;n,v(h) + o(h) = δ(n + 1)hPñ,ṽ;n+1,v(t) + pnhPñ,ṽ;n,v−1(t) + kλ(v + 1) d hPñ,ṽ;n−1,v+1(t) + c(v + 1)hPñ,ṽ;n,v+1(t) + (1− bh)Pñ,ṽ;n,v(t) + o(h), which leads to p′n,v(t) = δ ((n+ 1)pn+1,v(t)− npn,v(t)) + p (npn,v−1(t)− npn,v(t)) + kλ d ((v + 1)pn−1,v+1(t)− vpn,v(t)) + c ((v + 1)pn,v+1(t)− vpn,v(t)) , pn,v(0) = δn,ñδv,ṽ . Now the marginal expectations of the system are solutions of the fol- lowing ODE system. This is easily proved with arguments similar to those given in previous section. dE(T ∗(t)) dt = kλ d E(V (t))− δE(T ∗(t)) dE(V (t)) dt = pE(T ∗(t))− ( c+ kλ d ) E(V (t)) 2.6.2 Kolmogorov’s backward differential equations The backward Kolmogorov equations can be found just as quickly. Pñ,ṽ;n,v(t+ h) = Pñ,ṽ;ñ−1,ṽ(h)Pñ−1,ṽ;n,v(t) + Pñ,ṽ;ñ,ṽ+1(h)Pñ,ṽ+1;n,v(t) + Pñ,ṽ;ñ+1,ṽ−1(h)Pñ+1,ṽ−1;n,v(t) + Pñ,ṽ;ñ,ṽ−1(h)Pñ,ṽ−1;n,v(t) + Pñ,ṽ;ñ,ṽ(h)Pñ,ṽ;n,v(t) + o(h) = δñhPñ−1,ṽ;n,v(t) + pñhPñ,ṽ+1;n,v(t) + kλṽh d Pñ+1,ṽ−1;n,v(t) + cṽhPñ,ṽ−1;n,v(t) + (1− bh)Pñ,ṽ;n,v(t) + o(h). 41 2.6. Two-compartment stochastic model It follows that the backward KDE is given by P ′ñ,ṽ;n,v(t) = δñ ( Pñ−1,ṽ;n,v(t)− Pñ,ṽ;n,v(t) ) + pñ ( Pñ,ṽ+1;n,v(t)− Pñ,ṽ;n,v(t) ) + kλṽ d ( Pñ+1,ṽ−1;n,v(t)− Pñ,ṽ;n,v(t) ) + cṽ ( Pñ,ṽ−1;n,v(t)− Pñ,ṽ;n,v(t) ) , Pñ,ṽ;n,v(0) = δñ,nδṽ,v. 2.6.3 New properties of probability generation function in more than one dimension The next step will be handle the infinite number of ODE’s from the backward and forward equation using the probability generating function. The pgf G := G(x, y; t) := Gñ,ṽ(x, y; t) is defined (equivalently to the one dimensional model) as Gñ,ṽ(x, y; t) = E[x T ∗yV ] = ∞∑ v=0 ∞∑ n=0 Pñ,ṽ;n,v(t)x nyv. It has the same properties plus some new ones related to the marginal distributions. The new properties are: 1. We can recover the marginal probability of the number of infected cells, by computing G(x, 1; t) = ∞∑ n=0 ∞∑ v=0 Pñ,ṽ;n,v(t)x n. So defining T ∗n(t) as the probability of having n infected T cells at time t, i.e. T ∗n(t) := ∑ ∞ v=0 Pñ,ṽ;n,v(t), we obtain that T ∗n(t) = 1 n! ∂nG(x, 1; t) ∂xn ∣∣∣∣ x=0 . In particular, the probability that the infected cells are extinct at time t is given by ∑ ∞ v=0 Pñ,ṽ;0,v(t) = T ∗ 0 (t) = GT ∗(0, t) = G(0, 1; t). Similarly for the probability of the number of infected virus particles we get G(1, y; t) = ∞∑ v=0 ∞∑ n=0 Pñ,ṽ;n,v(t)y v. 42 2.6. Two-compartment stochastic model Again, letting Vv(t) be the probability of having v virus particles at time t, i.e. Vv(t) := ∑ ∞ n=0 Pñ,ṽ;n,v(t), we obtain that Vv(t) = 1 v! ∂vG(1, y; t) ∂yv ∣∣∣∣ y=0 . In particular, the probability of the virus to be extinct at time t is given by ∑ ∞ n=0 Pñ,ṽ;n,0(t) = V0(t) = GV (0, t) = G(1, 0; t). 2. We can calculate the mean value of the number of infected cell and virus particles, by evaluating E[T ∗] := ∞∑ n=0 nT ∗n(t) = ∂G(x, 1; t) ∂x ∣∣∣∣ x=1 and E[V ] := ∞∑ v=0 vVv(t) = ∂G(1, y; t) ∂y ∣∣∣∣ y=1 . 2.6.4 The probability generating function for the Kolmogorov forward equation In this case the definition of the pgf becomes slicker as we can abbreviate pn,v(t) = Pñ,ṽ;n,v(t). In order to find the pgf for the KFDE, we need to express all the terms in the KFDE with respect to the pgf G. To do so, let us multiply both sides by xnyv and sum over all n and v, where x, y ∈ D = {z ∈ C : |z| ≤ 1}, and n, v ∈ N0: ∞∑ n=0 ∞∑ v=0 p′n,v(t)x nyv = ∞∑ n=0 ∞∑ v=0 xnyv [ δ ((n+ 1)pn+1,v(t)− npn,v(t)) + p (npn,v−1(t)− npn,v(t)) + kλ d ((v + 1)pn−1,v+1(t)− vpn,v(t)) + c ((v + 1)pn,v+1(t)− vpn,v(t)) ] . We now rewrite this term-by-term, making use of the fact that the probabil- ity vanishes whenever one of the species has a negative number of individuals. 1. ∑ ∞ n=0 ∑ ∞ v=0 p ′ n,v(t)x nyv = ∂G/∂t. 2. δ ∑ ∞ n=0 ∑ ∞ v=0(n+ 1)pn+1,v(t)x nyv = δ ∑ ∞ n=0 ∑ ∞ v=0 npn,v(t)x n−1yv = δ(∂G/∂x). 43 2.6. Two-compartment stochastic model 3. δ ∑ ∞ n=0 ∑ ∞ v=0 npn,v(t)x nyv = δx(∂G/∂x). 4. p ∑ ∞ n=0 ∑ ∞ v=0 npn,v−1(t)x nyv = py ∑ ∞ n=0 ∑ ∞ v=0 npn,v(t)x nyv = pxy(∂G/∂x). 5. p ∑ ∞ n=0 ∑ ∞ v=0 npn,v(t)x nyv = px(∂G/∂x). 6. kλ d ∑ ∞ n=0 ∑ ∞ v=0(v + 1)pn−1,v+1(t)x nyv = kλ d ∑ ∞ n=0 ∑ ∞ v=0 vpn,vx n+1yv−1 = kλx d (∂G/∂y). 7. kλ d ∑ ∞ n=0 ∑ ∞ v=0 vpn,v(t)x nyv = kλy d (∂G/∂y). 8. c ∑ ∞ n=0 ∑ ∞ v=0(v + 1)pn,v+1(t)x nyv = c ∑ ∞ n=0 ∑ ∞ v=0 vpn,vx nyv−1 = c(∂G/∂y). 9. c ∑ ∞ n=0 ∑ ∞ v=0 vpn,v(t)x nyv = cy(∂G/∂y). Using these facts, we can eventually arrive at the following partial dif- ferential equation for the pgf: ∂G ∂t = [δ(1 − x) + px(y − 1)] ∂G ∂x + [ kλ d (x− y) + c(1− y) ] ∂G ∂y , G(x, y; 0) = xñyṽ. Let’s try to solve the pde by the method of characteristics [9], and sim- plify the notation as: x→ x1 y → y2 t→ x3 G→ u Then, the pde in the new notation becomes: F (x1, x2, x3, u, ux1 , ux2 , ux3) = − [δ(1 − x1) + px1(x2 − 1)] ux1 − [ kλ d (x1 − x2) + c(1− x2) ] ux2 + ux3 = 0, u(x1, x2, 0) = x ñ 1x ṽ 2. Then assuming that each of the variables on F are function of (s1, s2, t) and using the constraint t = 0, the characteristic equations become: 44 2.6. Two-compartment stochastic model dx1 dt = Fux1 = − [δ(1− x1) + px1(x2 − 1)] dx2 dt = Fux2 = − [ kλ d (x1 − x2) + c(1− x2) ] dx3 dt = Fux3 = 1 with x3(s1, s2, 0) = 0 ⇒ x3 = t du dt = 0 with u(s1, s2, 0) = s ñ 1s ṽ 2 ⇒ u(s1, s2, t) = s ñ 1s ṽ 2 This is a nonlinear ODE system of two equation to solve and then look for the values of s1 and s2 w.r.t. x1 and x2 in order to solve for u. That system has to be solved numerically so we stop here. Let’s then proceed to find the equation from the backward approach. 2.6.5 The probability generating function for the Kolmogorov backward equation Similarly to the section above we multiply both sides of the KBDE equation with xn, yv and sum over all n and v, where x, y ∈ D and n, v ∈ N0: ∞∑ n=0 ∞∑ v=0 dPñ,ṽ;n,v dt (t)xnyv = ∞∑ n=0 ∞∑ v=0 xnyv [ δñ ( Pñ−1,ṽ;n,v(t)− Pñ,ṽ;n,v(t) ) + pñ ( Pñ,ṽ+1;n,v(t)− Pñ,ṽ;n,v(t) ) + kλṽ d ( Pñ+1,ṽ−1;n,v(t)− Pñ,ṽ;n,v(t) ) + cṽ ( Pñ,ṽ−1;n,v(t)− Pñ,ṽ;n,v(t) )] . Using the definition of Gñ,ṽ and interpreting this a function that only de- pends on the time t, the above can be rewritten as dGñ,ṽ dt (t) = δñ ( Gñ−1,ṽ(t)−Gñ,ṽ(t) ) + pñ ( Gñ,ṽ+1(t)−Gñ,ṽ(t) ) + kλṽ d ( Gñ+1,ṽ−1(t)−Gñ,ṽ(t) ) + cṽ ( Gñ,ṽ−1(t)−Gñ,ṽ(t) ) Gñ,ṽ(0) = x ñyṽ, 45 2.6. Two-compartment stochastic model which is an infinite system of ODEs. Note that the ODE formulation is only possible because all involved probability functions are with respect to the same final data. This is not the case for the KFDE. A crucial step is to exploit the fact that in a branching process individuals act individually [21] so that Gñ,ṽ = (G1,0) ñ(G0,1) ṽ , for all ñ, ṽ ∈ N0. It therefore suffices to solve the ODE above for the func- tions G1,0 and G0,1. The corresponding equations are dG1,0 dt (t) = δ(1 −G1,0(t)) + p(G1,0(t)G0,1(t)−G1,0(t)), dG0,1 dt (t) = kλ d (G1,0(t)−G0,1(t)) + c(1−G0,1(t)), G1,0(0) = x, G0,1(0) = y. This system is actually the negative of the one obtained from the forward approach, as expected. Then I will need to find numerical solution of the system using, for example, the function ode45 from MATLAB. 2.6.6 Calculating the individual probabilities We will now use the properties of the pgf to explicitly calculate the proba- bility of having v virus particles at some time t. The same analysis can be done for the number of infected T cells at a certain time and the results are analogous. Recall from the definition of the pgf that Vv(t) := ∞∑ n=0 Pñ,ṽ;n,v(t) and that further Vv(t) = 1 v! ∂vG(1, y; t) ∂yv ∣∣∣∣ y=0 . As accurate numerical differentiation is difficult to perform numerically, in particular when we want to calculate hundreds or thousands of derivatives. Therefore, we turn towards complex analysis to find a different way to cal- culate the desired derivatives [5]. By Cauchy’s integral formula it holds that dn f dxn ∣∣∣∣ x=a = n! 2πi ∮ C f(z) (z − a)n+1 dz, 46 2.6. Two-compartment stochastic model whenever (i) C is a simple positively oriented close curve that encloses an open region which contains a, (ii) f is holomorphic on that region and (iii) f is continuous on its closure. In our case we have that f = G(1, y; t), as a function of y, is holomorphic on the interior of D and continuous on D, where D = {z ∈ C : |z| ≤ 1}. As we are interested in a = 0 we thus choose C to be the circle with radius 1 around the origin, and it follows that Vv(t) = 1 v! ( v! 2πi ∮ C G(1, z; t) zv+1 dz ) z=eiθ = 1 2πi ∫ 2pi 0 G(1, eiθ ; t) (eiθ)v+1 ieiθ dθ = 1 2π ∫ 2pi 0 G(1, eiθ ; t)e−ivθ dθ = 1 2π (∫ pi 0 G(1, eiθ ; t)e−ivθ dθ + ∫ 2pi pi G(1, eiθ ; t)e−ivθ dθ ) = 1 2π (∫ pi 0 G(1, eiθ ; t)e−ivθdθ + ∫ 0 pi G(1, ei(2pi−α) ; t)e−iv(2pi−α) (−dα) ) = 1 2π (∫ pi 0 G(1, eiθ ; t)e−ivθdθ + [∫ pi 0 G(1, eiα; t)e−ivα dα ] ∗ ) , where we substituted α = 2π − θ as well as G(1, ei(2pi−α) ; t) = ∞∑ v=0 Vv(t)e iv(2pi−α) = ∞∑ v=0 Vv(t)e −ivα = [ G(1, eiα; t) ] ∗ Here the superscript ∗ denotes the complex conjugate. Substituting this identity we finally arrive at Vv(t) = 1 π Re (∫ pi 0 G(1, eiθ ; t)e−ivθ dθ ) . We used the built-in ODE solver ode45 and the integrator trapz from MATLAB to compute the probability Vv(t) of having v virus particles at time t. 2.6.7 Numerical simulations In order to compare with the data, this section will presents results from two numerical approaches. The first one consists of the simulation of the pro- cess using Gillespie’s algorithm as before. The second one is the numerical solution of the pgf and the integrals described in the previous section. 47 2.6. Two-compartment stochastic model The Gillespie’s algorithm was again coded in C. There was one infected cell and zero virus as an initial condition. The simulation was repeated 10000 times. Only 75 of the 10000 trajectories were not extinct after 20 days. Some of the trajectories of the number of infected cells and virus are plotted on Figures 2.36 and 2.37 respectively. Figures 2.38 and 2.39 show the maximum values of each populations. 0 5 10 15 20 0 5 10 15 20 25 Days T* p er .5 L Figure 2.36: Different paths of T ∗ from the T ∗V model numeri- cal simulation 0 5 10 15 20 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 Days VL p er .5 L Figure 2.37: Different paths of V from the T ∗V model numerical simulation 0 5 10 15 20 0 10 20 30 40 50 60 Days T* p er .5 L Figure 2.38: Maximum of T ∗ from the T ∗V model numerical simula- tion 0 5 10 15 20 0 2 4 6 8 10 12 14 x 10 4 Days VL p er .5 L Figure 2.39: Maximum of V from the T ∗V model numerical simula- tion The maximum values of the infected cells are similar to the ones observed on previous model (Figures 2.25 and 2.38. It is also observable that the total 48 2.6. Two-compartment stochastic model viral load reached maximum values of 12 × 104 per monkey in less than 20 days (Figure 2.39). This represents a viral load of 240 per ml. Comparing the maximum observed with the values of the data (Table 2.4) there are two observable viral blips with a probability greater than 0.0001. Those viral blips are the same that were observables with previous model. It is also valid to compare with the data even if the simulation goes only for 20 days because only one of the data blips has it maximum after more than 20 days ( third viral blip of P177 at day 23). The rest of the data have the maximum before 15 days. Similar to the previous section, plots of the mean and the variance were made for this process. The mean and variance of the infected cells are similar to the one observed previously. The ones from this model are plotted on Figures 2.40 and 2.42. It is remarkable that the variance of the virus is also comparable to the one of the infected cells (Figure 2.43). The mean of the virus presents a rapid increase followed by a more slow decay (Figure 2.41). This is in agreement with the solution of the corresponding ODE model. 0 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Days T* p er .5 L Figure 2.40: Mean values of T ∗ from the T ∗V model numerical simulation 0 5 10 15 20 0 500 1000 1500 2000 2500 Days VL p er .5 L Figure 2.41: Mean values of V from the T ∗V model numerical simulation The first time to extinction was also calculated. In this process, the first time to extinction is defined as the first time that T ∗ = 0 and V = 0. In the simulations, the mean at 20 days was 2.0518 days but there were 75 not extinct paths at that time. Therefore a lower bound for the mean time of extinction is 2.0518 + 20 ∗ 7510000 = 2.2018 days. From that mean value and Figure 2.44 it is noticed that the first time to extinction function is slower compared with the single-compartment T ∗ model. At time 20, there were 75 non extinct paths. The histogram of those is on Figure 2.45. There is an observable bimodal distribution. This is because 49 2.6. Two-compartment stochastic model 0 5 10 15 20 0 0.5 1 1.5 2 2.5 3 3.5 4 Days T* p er .5 L Figure 2.42: Variance of T ∗ from the T ∗V model numerical simula- tion 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 7 Days VL p er .5 L Figure 2.43: Variance of V from the T ∗V model numerical simula- tion 0 5 10 15 20 0 20 40 60 80 100 120 Days # of e xt s im ul at io ns Mean of simulations = 2.0518 # Not extinct = 75 Figure 2.44: First time to extinc- tion of the system from the T ∗V model numerical simulation 0 2 4 6 8 x 104 0 5 10 15 20 25 30 Number of virus # of n on e xt s im ul at io ns Figure 2.45: Histogram of the non extinct viral load in simulations at t = 20 days 50 2.6. Two-compartment stochastic model zero is the attractor of the system and all the paths will go extinct, but some of them can run for longer times. In general, those paths that were running for longer times would present higher values of virus and then it will take them longer to go back to zero. The longer paths will then have a quasi stationary mean different than zero. Then there is a mode at zero and a mode at the quasi stationary mean. As time goes to infinity, the distribution will converge to a unimodal distribution centered at zero. 2.6.8 Comparing the viral load from simulation with data In order to compare the results with the data, we calculate the probability of observing a viral blip from the 10000 Gillespie simulations. The detectable level on monkey basis is given by the conversion from 40 RNA copies/ml to 20000 copies/l. The normalized frequency of having more than 20000 copies/L on the virus trajectories over the 10000 simulations is an estimate the probability of the having a blip (Figure 2.46). 0 5 10 15 20 0 0.005 0.01 0.015 Days P( V> 20 00 0) Figure 2.46: Probability of having more than 20000 virus per monkey as a function of time After 20 days of simulation there remains a probability of more than 0.002 of observing a blip. The maximum at time 20 is around 7 ∗ 104 which is equivalent to around 140 copies per ml. That number of virus is enough to fit only one of the data blips but with longer duration than the observed in the data, were it goes only for 13 days (Table 2.4). Using the marginal distribution, it is possible to estimate numerically the probability of having N virus at time t1. Remember that the marginal distribution can be found as derivatives from the pgf and the derivatives are calculated using the Cauchy integral formula. Since it is needed to calculate 51 2.6. Two-compartment stochastic model the total viral load per monkey, the maximum viral blip from the data is again transformed from copies per ml to copies per monkey. The number of derivatives needed to find the smallest data blips are of order 104. The number of derivatives needed to calculate the other data blips are more than 6 × 105, then there will take at least 10 times more time than the smaller ones. It is computationally intensive to calculate many derivatives even with the use of the Cauchy integral formula. Therefore we decided to calculate only the probability of the two smallest data blips and the results are in Table 2.7. Table of probabilities of blips N t1 P (V (t1) = N |T ∗ 0 = 1, V (0) = 0) P157 120*500 6 7.6781e-08 160*500 8 1.8596e-08 Table 2.7: Probabilities of having N infected cells after t1 days post infection starting with one infected cell at time zero using the T ∗V model. This probabilities are smaller comparing with the observed with the T ∗ model (Table 2.5). This is because the faster dynamic of the viral load compared with the infected cells. The viral load fluctuates a lot in the timescale of one day. Then the probability of observe a certain number of virus at a fixed time is smaller since it has more possible results. On the other hand, the T ∗ model assumes a constant number of virus per infected cell. 2.6.9 Discussion of the model The T ∗V model has more realistic assumptions than the ones from the T ∗ model. It shows how the viral dynamics could lead to a longer extinction times. However it presents similar results to the ones obtained with the past model and it does not have an analytic solution. Those similarities a long with the results from simulations enforce the hypothesis of a close relation between the infected cells and the viral load. From the path’s figures it is clear that qualitatively the virus and the infected cells have the same behavior. As observed in the previous subsection, there is only one possible ob- servable data blip from the 10000 Gillespie simulation. Also the marginal probability of the virus is giving really small probabilities of observing the 52 2.7. Latent cells stochastic model small viral blips. In future work we will calculate the probability of having the viral load in a certain range instead of at a specific number. This is not a bad idea, since the experimental data has some margin of error. This model also supports the idea that there must be more than one latent cell active to generate the data blips. That will be modeled in the next section. 2.7 Latent cells stochastic model From the previous section we realize that one latent cell is probably not enough to produce the data blips. We will extend those models by including the latent cells directly, ie, including the latent cell in an inactive stage. We will then do analogous analysis. Figure 2.47 illustrated the model. T ∗ Vp δL ηkλd L a c (1− η)kλ d ρ δ Figure 2.47: Latent cells stochastic model The possible reactions for the new process are: An infected cell dies T ∗ δ −→ ∅ A new virus particle is produced T ∗ p −→ T ∗ + V An uninfected T cell becomes infected V (1−η) kλ d−−−−−→ T ∗ A virus particle is cleared V c −→ ∅ A uninfected T cell becomes infected V η kλ d−−→ L and goes to latent state. A latent cell activates L a −→ T ∗ A latent cell reproduces L ρ −→ 2L A latent cell dies L δL−→ ∅ 53 2.7. Latent cells stochastic model In this section, as in previous ones, this process will also be assumed to be a stationary continuous time Markov chain with discrete states. Define the transition probabilities to be P ñ,ṽ,l̃;n,v,l(t) = P (T ∗(t) = n, V (t) = v, L(t) = l|T ∗(0) = ñ, V (0) = ṽ, L(0) = l̃) Similarly to before, let’s assume that the probability of more than one event in the time interval (0, h) is in o(h). Then the following postulates hold: Pn+1,v,l;n,v,l(h) = δ(n + 1)h+ o(h), Pn,v−1,l;n,v,l(h) = pnh+ o(h), Pn−1,v+1,l;n,v,l(h) = (1− η) ( kλ(v + 1) d ) h+ o(h), Pn,v+1,l:n,v,l(h) = c(v + 1)h + o(h) Pn,v+1,l−1;n,v,l(h) = η ( kλ(v + 1) d ) h+ o(h), Pn−1,v,l+1;n,v,l(h) = a(l + 1)h + o(h) Pn,v,l−1;n,v,l(h) = ρ(l − 1)h + o(h) Pn,v,l+1;n,v,l(h) = δL(l + 1)h+ o(h) Now following similar methods to the previous sections, we find backward and forward equations of the transition probabilities and derive equations for the probability generating function. Let T ∗(0) = ñ, V (0) = ṽ and L(0) = l̃ be the initial condition for the pro- cess. TheKolmogorov forward equation for the transition probabilities, where pn,v,l(t) = Pñ,ṽ,l̃;n,v,l(t), is as follows: p′n,v,l(t) = δ ((n+ 1)pn+1,v,l(t)− npn,v,l(t)) + p (npn,v−1,l(t)− npn,v,l(t)) + (1− η) kλ d ((v + 1)pn−1,v+1,l(t)− vpn,v,l(t)) (2.27) + c ((v + 1)pn,v+1,l(t)− vpn,v,l(t)) + a((l + 1)pn−1,v,l+1 − lpn,v,l) + η kλ d ((v + 1)pn,v+1,l−1(t)− vpn,v,l(t)) + ρ((l − 1)pn,v,l−1 − lpn,v,l) + δL((l + 1)pn,v,l+1 − lpn,v,l) pn,v,l(0) = δn,ñδv,ṽδl,l̃. 54 2.7. Latent cells stochastic model As in previous sections, the marginal expectations of the process are solutions of the equivalent ODE system. This can be proved using the forward equation. Defining the pgf as G(x, y, z) = E(xT ∗ yV zL) and using the forward equa- tion and similar techniques as before we get the following equation for the pgf: ∂G ∂t = [δ(1 − x) + px(y − 1)] ∂G ∂x + [ (1− η) kλ d (x− y) + η kλ d (z − y) + c(1− y) ] ∂G ∂y + [a(x− z) + ρz(1− z) + δL(1− z)] ∂G ∂z , G(x, y; 0) = xñyṽ. On the other hand, theKolmogorov backward equation for the tran- sition probabilities is given by: P ′ ñ,ṽ,l̃;n,v,l (t) = δñ ( P ñ−1,ṽ,l̃;n,v,l(t)− Pñ,ṽ,l̃;n,v,l(t) ) + pñ ( P ñ,ṽ+1,l̃;n,v,l(t)− Pñ,ṽ,l̃;n,v,l(t) ) + ( 1− η )kλṽ d ( P ñ+1,ṽ−1,l̃;n,v,l(t)− Pñ,ṽ,l̃;n,v,l(t) ) + cṽ ( Pñ,ṽ−1,l̃;n,v,l(t)− Pñ,ṽ,l̃;n,v,l(t) ) , + η kλṽ d ( Pñ,ṽ−1,l̃+1;n,v,l(t)− Pñ,ṽ,l̃;n,v,l(t) ) + al̃ ( Pñ+1,ṽ,l̃−1;n,v(t)− Pñ,ṽ,l̃;n,v(t) ) + ρl̃ ( P ñ,ṽ,l̃+1;n,v(t)− Pñ,ṽ,l̃;n,v(t) ) + δL l̃ ( P ñ,ṽ,l̃−1;n,v(t)− Pñ,ṽ,l̃;n,v(t) ) Pñ,ṽ,l̃;n,v(0) = δñ,nδṽ,vδl̃,l The corresponding system of equations for the pgf is given by the fol- lowing equations. Here we have used again the branching property, G ñ,ṽ,l̃ = (G1,0,0) ñ(G0,1,0) ṽ(G0,0,1) l̃. 55 2.7. Latent cells stochastic model dG1,0,0 dt (t) = δ(1 −G1,0,0(t)) + p(G1,0,0(t)G0,1,0(t)−G1,0,0(t)), dG0,1,0 dt (t) = (1− η) kλ d (G1,0,0(t)−G0,1,0(t)) + η kλ d (G0,0,1(t)−G0,1,0(t)) + c(1−G0,1,0(t)), dG0,0,1 dt (t) = a(G1,0,0 −G0,0,1) + δL(1−G0,0,1) + ρ(G 2 0,0,1 −G0,0,1) G1,0,0(0) = x, G0,1,0(0) = y, G0,0,1(0) = z. This system cannot be solved analytically but it is possible to find nu- merical solutions. We can calculate the marginal probabilities by taking derivatives of the pgf. Now, analogous to subsection 2.6.6, we simplify the computational effort by integrating instead of taking derivatives. Then the probability of having v viruses at time t is: P (V (t) = v) = 1 π Re (∫ pi 0 G(1, 1, eiθ ; t)e−ivθ dθ ) . 2.7.1 Numerical simulations For the numerical simulations we will first simulate the model with just one latent cell and similar parameters as previous models. The second step will be to modify those parameters and the initial conditions to evaluate possible mechanisms to observe the data blips. I performed 10000 Gillespie simulations of the process with parameters as in Table 2.3. I assumed initial conditions: L = 1, T ∗ = V = 0. Examples of the trajectories of each of the populations are on Figures 2.48. 2.49 and 2.50. There is an observable close relation between the infected cells and the viruses, as observed with the previous model. We also plotted the maximum values, mean and variance of the latent cells, infected cells, and virus particles over the 10000 simulations. The maximum values are illustrated on Figures 2.51, 2.52, and 2.53. The maximum value of possible latent cells after 20 days is 10. The infected cells were no more than 50 after 20 days. We obtain similar infected cell’s maximum in previous models. The maximum number of viruses is around 56 2.7. Latent cells stochastic model 0 5 10 15 20 0 5 10 15 20 25 Days T* p er .5 L Figure 2.48: Different paths of T ∗ from the LT ∗V model Gillespie simulations 0 5 10 15 20 0 1 2 3 4 5 6 x 10 4 Days VL p er .5 L Figure 2.49: Different paths of V from the LT ∗V model Gillespie simulations. The detection limit is plotted as a dotted line. 0 5 10 15 20 0 1 2 3 4 5 6 7 Days L pe r . 5L Figure 2.50: Different paths of L from the LT ∗V model Gillespie simulations 105 which is smaller than the one from the T ∗ model. Then again we can only observe the two smaller data blips with a probability greater than 0.0001 (Table 2.4). From the figures it is observable that both the infected cells and the virus are increasing. This is because the reservoir of latent cells is making the dynamics slower than the previous models. The mean plots correspond to Figures 2.54, 2.55, and 2.56. Comparing the mean of the infected cells here with the results from previous models 2.40, they are completely different. All this is because the latent cells. They are creating an extra source pool for virus and infected cells. That is why they first increase and then they will decrease, but this decreasing is not 57 2.7. Latent cells stochastic model 0 5 10 15 20 0 5 10 15 20 25 30 35 40 45 50 Days T* p er .5 L Figure 2.51: Maximum of T ∗ from the LT ∗V model Gillespie simula- tions 0 5 10 15 20 0 1 2 3 4 5 6 7 8 9 10 x 10 4 Days VL p er .5 L Figure 2.52: Maximum of V from the LT ∗V model Gillespie simula- tions 0 5 10 15 20 1 2 3 4 5 6 7 8 9 10 Days VL p er .5 L Figure 2.53: Maximum of L from the LT ∗V model Gillespie simulations 58 2.7. Latent cells stochastic model observable over the first 20 days. The latent cells do present a decreasing mean behavior. 0 5 10 15 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Days T* p er .5 L Figure 2.54: Mean of T ∗ from the LT ∗V model Gillespie simulations 0 5 10 15 20 0 100 200 300 400 500 600 700 800 900 Days VL p er .5 L Figure 2.55: Mean of V from the LT ∗V model Gillespie simulations 0 5 10 15 20 0.4 0.5 0.6 0.7 0.8 0.9 1 Days L pe r . 5L Figure 2.56: Mean of L from the LT ∗V model Gillespie simulations The variance plots of the simulations are on Figures 2.57, 2.58 and 2.59. They show relatively similar results to the observed before. The first time to extinction in this case correspond to the time where all the latent cells, infected cells and virus are extinct. The mean time to extinction at time 20 was 5.9624 days, but there remain 2718 paths that were no extinct. A lower bound for the time of extinction will be 5.9624 + 20× 271810000 , thus the process will take at least 11.3984 days on average to go extinct. This model shows the longest mean to extinction comparing with the previous ones, where the mean time of extinction was around 2 days. Figure 2.60 shows the fist passage time to extinction from the simulation. We observe a lot of fluctuation on the distribution. This is because the 59 2.7. Latent cells stochastic model 0 5 10 15 20 0 0.5 1 1.5 2 2.5 3 3.5 4 Days T* p er .5 L Figure 2.57: Variance of T ∗ from the LT ∗V model Gillespie simula- tions 0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 x 10 6 Days VL p er .5 L Figure 2.58: Variance of V from the LT ∗V model Gillespie simula- tions 0 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Days L pe r . 5L Figure 2.59: Variance of L from the LT ∗V model Gillespie simulations 60 2.7. Latent cells stochastic model extinction is slower than before. The latent cells die slowly, thus once you get one, even if it dies without activating, it will take longer for the system to go extinct. Figure 2.61 shows the cumulative distribution of the first time to ex- tinction. From here we can observe that the extinction is increasing but for some small interval of time near zero it is not increasing. That is the reason we observe the big variance and fluctuations in Figure 2.60 0 5 10 15 20 0 2 4 6 8 10 12 14 Days # of e xt s im ul at io ns Mean of simulations = 5.9624 # Not extinct = 2718 Figure 2.60: First passage time to extinction of LT ∗V model 0 5 10 15 20 0 1000 2000 3000 4000 5000 6000 7000 8000 Days # of e xt s im ul at io ns Mean of simulations = 5.9624 # Not extinct = 2718 Figure 2.61: Cumulative distribu- tion of the FPT to extinction of LT ∗V model 2.7.2 Comparing with data This model includes activation, reproduction and death of latently infected cells, changing deeply the results and dynamics of the system. In this model we can observe that there are multiple viral blips without initiating the system again because at least one latent cell remain inactive while the rest of the infected active cells and virus goes extinct. This is what we observe on the data. There is also the possibility that during 20 days the first latent cell never activates or dies, and then no virus is generated. These different dynamics are illustrated in Figures 2.62, 2.63, and 2.64. As in the previous model, we estimate the probability of having an ob- servable blip for each time from the frequency of simulated trajectories with more than 20000 virus. The normalized frequency of having more than 20000 viruses on the virus trajectories over the 10000 is on Figure 2.65. After 20 days there is a probability of approximately 0.008 of having a viral blip. This is bigger than the one observed on T ∗V model (Figure 2.46). On the other hand the probability of having a viral blip after 5 days 61 2.7. Latent cells stochastic model 0 5 10 15 20 0 1 2 3 4 5 6 7 Days T* p er .5 L Figure 2.62: LT ∗V model paths showing special dynamics, T ∗ plotted. The red and blue paths presents more than one viral blip. The green path shows the on acti- vation of L before 20 days. 0 5 10 15 20 0 2000 4000 6000 8000 10000 12000 14000 Days VL p er .5 L Figure 2.63: LT ∗V model paths showing special dynamics, V plot- ted. The red and blue paths presents more than one viral blip. The green path shows the on acti- vation of L before 20 days. 0 5 10 15 20 0 1 2 3 4 5 6 Days L pe r . 5L Figure 2.64: LT ∗V model paths showing special dynamics, L plotted. The red and blue paths presents more than one viral blip. The green path shows the on activation of L before 20 days. 62 2.7. Latent cells stochastic model is around 0.003 for the LT ∗V model while in the T ∗ model is around 0.013. This is again a consequence of the long lived latent cells. 0 5 10 15 20 0 0.002 0.004 0.006 0.008 0.01 0.012 Days P( V> 20 00 0) Figure 2.65: Probability of having more than 20000 virus per monkey as a function of time From the observation of the Gillespie simulations, we can conclude that one latent cell is not enough to produce the observable blips from the date. Thus we will modify the model a bit to test other hypothesis. 2.7.3 Numerical simulation with fast activation rate We were using a slow activation which drives the system to slower dynamics. In this section we will increase that activation rate. This assumption could be caused because secondary infections, a non independent activation of latent cell or other factors. Non independent activation of latent cells could happen when they activate as a group or when one latent cell activates it starts a chain reaction that activate a few other latent cells. From the T ∗ model we estimated that the average number of infected cells needed to observe the present blips is 10. Then we will use this estimate as our initial number of latent infected cells. We wish that all of them activate in around one day, then a = 10 day−1 will be our new activation rate. In order to simulate just one viral blip we will neglect the reproduction and death parameters. This is in order to compare with previous models, where we were analyzing one viral blip per trajectory. Thus, we let ρ = 0 and δL = 0. The rest of the parameters remain as in Table 2.3. We simulated 10000 paths using Gillespie’s algorithm. A few of the trajectories of each of the populations are on Figures 2.66. 2.67 and 2.68. 63 2.7. Latent cells stochastic model As showed with the slow activation parameters, there is a close relation between the number of productive infected cells and the viral load. 0 5 10 15 20 0 5 10 15 20 25 30 Days T* p er .5 L Figure 2.66: Different paths of T ∗ from the LT ∗V model simulations with fast activation rate. 0 5 10 15 20 0 1 2 3 4 5 6 x 10 4 Days VL p er .5 L Figure 2.67: Different paths of V from the LT ∗V model simulations with fast activation rate. The de- tection limit is plotted as a dotted line. 0 5 10 15 20 0 1 2 3 4 5 6 7 8 9 Days L pe r . 5L Figure 2.68: Different paths of L from the LT ∗V model simulations with fast activation rate. Notice that many of the simulations go extinct before 10 days. Also notice from Figure 2.67 that exceeding the detectable level is more common than in the previous models. In order to understand the model better we calculate the mean, variance and maximum values for the fast activation parameters. The maximum values as a function of time are in Figures 2.69, 2.70, and 64 2.7. Latent cells stochastic model 0 5 10 15 20 0 10 20 30 40 50 60 70 Days T* p er .5 L Figure 2.69: Maximum of T ∗ from the LT ∗V model simulations with fast activation rate. 0 5 10 15 20 0 2 4 6 8 10 12 14 x 10 4 Days VL p er .5 L Figure 2.70: Maximum of V from the LT ∗V model simulations with fast activation rate. 0 5 10 15 20 0 1 2 3 4 5 6 7 8 9 Days VL p er .5 L Figure 2.71: Maximum of L from the LT ∗V model simulations with fast activation rate. 65 2.7. Latent cells stochastic model 2.71. It can be observed that fast latent cell activation immediately reduces the latent cell to low levels as expected in this model. It was also expected that the maximum values of the infected cells and viruses should be higher than in previous models. In this case we have a maximum of infected cells around 65 per monkey. Comparing Figure 2.70 with Table 2.4 we find that there is again only the two small data blips with probability greater than 0.0001. The third small blip will consist in 1260 = 63 × 104 viral, which is more than 4 times bigger than the maximum observed from simulations. The mean values of the process are on Figures 2.72, 2.73 and 2.74. In this case the mean behaviors present a rapidly decreasing population of the latent cells. There are almost no latently infected cells after 1 day, as expected from the definition of the activation rate. Notice that even from Figure 2.67 we observe viral load above the detection limit but the mean of the viral load remains under detectable levels (20000 copies per monkey). 0 5 10 15 20 0 1 2 3 4 5 6 7 8 9 Days T* p er .5 L Figure 2.72: Mean of T ∗ from the LT ∗V model simulations with fast activation rate. 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 4 Days VL p er .5 L Figure 2.73: Mean of V from the LT ∗V model simulations with fast activation rate. The variance of the process is illustrated in Figures 2.75, 2.76, and 2.77. The variance of the infected cells and virus has similar dynamics as in pre- vious models. The variance of the latent infected cells is really small after one day, which again reflects the fact that the probability of activating all the latent cell in one day is almost one. This also implies that the mean behavior of the latent cells describes closely the dynamics of the latent cell population. Thus, after one day, we are again in a situation with almost no latent cell on the body of the monkey. The first passage time to extinction is calculated in Figure 2.78 and its cumulative distribution is in Figure 2.79. As in the simulation of previous 66 2.7. Latent cells stochastic model 0 5 10 15 20 0 1 2 3 4 5 6 7 8 9 Days L pe r . 5L Figure 2.74: Mean of L from the LT ∗V model simulations with fast activa- tion rate. 0 5 10 15 20 0 5 10 15 20 25 30 35 40 Days T* p er .5 L Figure 2.75: Variance of T ∗ from the LT ∗V model simulations with fast activation rate. 0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 x 10 7 Days VL p er .5 L Figure 2.76: Variance of V from the LT ∗V model simulations with fast activation rate. 67 2.7. Latent cells stochastic model 0 5 10 15 20 0 0.5 1 1.5 2 2.5 Days L pe r . 5L Figure 2.77: Variance of L from the LT ∗V model simulations with fast activation rate. 68 2.7. Latent cells stochastic model set of parameters, the latent cell population creates a more fluctuating first passage distribution. The cumulative distribution shows that at time 20, the probability of being extinct is greater than 0.9 while the slow activation rate gave only a .75 probability of extinction after 24 days (Figure 2.61). 0 5 10 15 20 0 5 10 15 20 25 Days # of e xt s im ul at io ns Mean of simulations = 7.0672 # Not extinct = 696 Figure 2.78: First passage time to extinction of LT ∗V model 0 5 10 15 20 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Days # of e xt s im ul at io ns Mean of simulations = 7.0672 # Not extinct = 696 Figure 2.79: Cumulative distribu- tion of the FPT to extinction of LT ∗V model The mean time to extinction from the simulation at t = 20 days was 7.0672 days but 696 paths were not extinct by this time. A lower bound for the mean time of extinction will be 7.0672 + 696 ∗ 20/10000, which corre- sponds to 8.4592 days. This is faster than the expected from the data where the mean duration time of a viral blip is 20 days. Then we estimate the probability to observe a blip from the frequency of observing more than 20000 virus particles. The normalized distribution is in Figure 2.80. For this process, the probability of observing a viral blip at day one is greater than 0.4 compared with that obtained from the T ∗V model (around 0.003 (Figure 2.46)). At the same time there is no difference in the probability of observing a viral blip at day 20 between the fast activation latent model and the T ∗V model. The fast activation latent model has an increased the probability of ob- serving a blip. It also accelerates the time to observe a viral blip. These two were the objective of the fast activation modification. 2.7.4 Discussion of the model From the simulation with parameters from Table 2.3, we observed that the reproduction of the latent cells along with the slow activation rate is gener- ating multi-blips in a same trajectory. This allows us to simulate multiple 69 2.7. Latent cells stochastic model 0 5 10 15 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Days P( V> 20 00 0) Figure 2.80: Probability of having more than 20000 virus per monkey as a function of time sequential blips as observed in the monkey data. The disadvantage of this model is that the low activation rate prevents us from observing viral blips at the beginning of the process. The low activation rate model also was not generating enough viruses to compare with the data. Thus the assumption of more than one latent cell was tested. Assuming a fast activation rate, no reproduction of latent cells and mul- tiple latent cells present at time zero, we could no longer find multi-blips in a single trajectory. On the other hand, the new parameter assumptions lead to a greater probability of observing a viral blip. This was still not enough to reproduce the observed maximun viral load. The fast activation model also showed a mean time to extinction of around 8 days, while the data showed a mean duration of 24 days. A mix of the two previous assumptions will need to be tested. From the first set of parameters, the mix between the slow activation and the reproduction of latent cells produces multiple blips in a path. From the second assumptions, starting with one infected cell that has been activated along with the fast activation rate produce immediate viral blips. Then a mix of a medium activation rate, a normal reproduction of latent cell rate, an infected activated cell with several latent cells as initial conditions, should lead to a better description of the data. 70 2.8. Conclusions 2.8 Conclusions In this chapter we modeled in detail one possible mechanism driving a vi- ral blips. From all models, we conclude that more than one latent cell is probably needed to observe viral blips as high as the ones from the data. From the T ∗ model we estimated that the mean number of latent cells to be activated should be around 10 per monkey but when analyzing the fast activation latent model we realized than more than 10 latent cells are needed. At the same time, modifications of the mean time to extinction were ob- served. The T ∗ model shows that in order to have a mean time to extinction of 24 days we need to modify the parameters such that R0 is almost one. But such a model is very sensitive to its parameters. The first latent model present a mean time to extinction closer, than the other models, to the data values, and without strong parametric sensitivity. Using the characteristics of the Auranofin drug, a fast activation of the latent cell was tested. In this case, the probability of having a viral blip over the detection limit is greater than with the other models. The latent model seems to be a good model to describe the data but improvements need to be done. One general improvement to all of the models would be to modify another parameter to make R0 less than one. Currently, I modify the parameter λ. This is an artificial assumption and there are other more natural assumptions using the characteristics of the H-iART treatment. The T ∗ model works as an start point to imagine the dynamics of more complicated models. The T ∗V model was in closed agreement with the T ∗ model results. With the latent model we proved that more than one latent cell is needed to observe the data blips. It was also proved from the latent model that one latent cell can drive multiple blips as observed in the data. 71 Chapter 3 Immune response and memory activation 3.1 Model description As mentioned before, a persistent viral reservoir in patients receiving potent ART prevents the eradication of HIV infection. This is known because on treatment interruption there is a smooth recovery of the viral load to pre- treatment levels, as in Figure 3.1. However, the drug Auranofin has been shown to accelerate the activa- tion rate of latent cells and alter the kinetics of viral rebound when drug treatment is interrupted [13]. This drug is an oral drug adopted for treat- ment of rheumatoid arthritis, a disease that is particularly characterized by increased central memory T cells. Auranofin action is partially understood. It has been observed that: • It impairs the proliferative capacity of T lymphocytes in vitro. • It decreases the production of pro-inflammatory cytokines in macrophages and T-cells. • It induces apoptosis in the Jurkat T-cell line. • It affects T cell proliferation and survival in vivo. • There are insights that it does not affect T cell regeneration. • It contributes to apoptosis and differentiation of leukemic cells in vitro in combination with anticancer therapies. This effects are making Auranofin a valuable candidate to reduce the reservoir of latently infected cells. According to [13], adding Auranofin be- fore stopping the treatment modified the rebound of the viral load as shown in Figure 3.2. 72 3.1. Model description Figure 3.1: Normal dynamics of the viral load after treatment in- terruption Figure 3.2: Observed treatment interruption dynamics of the vi- ral load when Auranofin is also in- cluded in the treatment regimen We want to develop reasonable models that reproduce the observed viral rebounds in the presence and absence of Auranofin. We hypothesize that the faster activation of the memory cells is decreasing the latent reservoir and is almost erasing the memory of the immune system. This depletion of memory immunity along with lower levels of CD4 T cells could explain the rebound observed after treatment interruption with Auranofin included. In order to achieve that, in a simplistic model, we introduce memory cell compartments (infected and non-infected) and a black-box immune activa- tion compartment. In the current version of the model, immune activation increases T cell reproduction, the rate of killing of infected cells and the rate of activation of memory cells. The immune activation is assumed to be driven by the presence of productively infected cells. The current model is designed to reproduce two important observations. The first one is that without Auranofin the viral load goes back smoothly and monotonically to the pre-treatment steady state as in Figure 3.1. The second one is that adding Auranofin to the treatment results in the number of memory cells being reduced substantially. As a consequence, the immune response after treatment interruption responds as slowly as primary infection with an early infection’s characteristic viral peak illustrated on Figure 3.2. It is not clear if this model will reproduce this effect with reasonable parameter values. 73 3.1. Model description c ρ dm λI dm km(1− ǫ) p δ(I) fa(T ∗) M∗ dI λ(I) M I V T ∗ a(I) a(I) λm k(1− ǫ) ρ d λmT Figure 3.3: Graphic of the model. Solid lines represent constant parameters. Dashed lines represent function of I. Dotted lines represent functions of T ∗ The corresponding system of differential equations reads: dT dt =λ(I)− k(1− ǫ)V T − dT + a(I)M − λmT dT ∗ dt =k(1− ǫ)V T − δ(I)T ∗ + a(I)M∗ − λmT ∗ dM dt =− a(I)M − dmM + λmT + ρM − km(1− ǫ)VM dM∗ dt =km(1− ǫ)VM − dmM ∗ − a(I)M∗ + λmT ∗ + ρM∗ dV dt =pT ∗ − cV dI dt =λI + fa(T ∗)− dII (3.1) The compartments T cells (T), infected T cells (T∗), memory T cells (M), infected memory cells (M∗), virus (V) and the action of the immune response (I). Here the parameters are described on Table 3.1. 74 3.1. Model description λ(I) Regeneration of T-cells is assumed to be function of I. k Infection rate of T cells. ǫ Drug efficacy. d Rate of death of T cells. a(I) Activation rate of memory cells to T cells. Assumed to be the same on infected and uninfected cells. λm Rate at which T cells convert to memory cells. Assumed to be the same for infected and uninfected cells. δ(I) Death of infected T cells is assumed to be function of I. dm Death rate of the memory cells. Assumed to be the same for infected and uninfected cells. km Infection rate of memory cells. ρ Reproduction of memory cells. p Viral production rate per infected cell c Viral clearance rate. λI Rate of regeneration of immunity not depending on the infection. Representing the immune system before infection. fa(T ∗) Generation of immunity driven by the infected T cells. dI Death rate of the immunity. Table 3.1: Parameter description of model [3.1] All the functions are assumed to be Hill functions: λ(I) =λ+ (λmax − λ)(I − λI d1 )2 (I − λI d1 )2 + 1 δ(I) =δ + (δmax − δ)(I − λI d1 )2 (I − λI d1 )2 + 1 a(V ) =a+ (amax − a)(I − λI d1 )2 (I − λI d1 )2 + 1 fa(T ∗) = fmaxT ∗ 2 T ∗2 + β2f (3.2) 75 3.2. Analytical work 3.2 Analytical work I tried to calculate the analytical expression for the steady state values. I was not able to find all of them. The disease free equilibrium is given as in (3.3). This equilibrium will be used to estimate some parameters further on. T̂ = λ(Î)(a(Î) + dm − ρ) (λm + d)(dm − ρ) + da(Î) = λ(a+ dm − ρ) (λm + d)(dm − ρ) + da T̂ ∗ =0 M̂ = λmλ(Î) (λm + d)(dm − ρ) + da(Î) = λmλ (λm + d)(dm − ρ) + da M̂∗ =0 V̂ =0 Î = λI dI (3.3) The R0 value for this model is given in (3.4). It will be used as well on the estimation of some parameter values. It was calculated using the next generation matrix with similar calculations as in previous chapter. R0 = p(1− ǫ) c kmM̂a(Î) + kT̂ (dm + a(Î)− ρ) δ(Î)(dm + a(Î)− ρ) + λm(dm − ρ) = p(1− ǫ)λ(Î) c(δ(Î)(dm + a(Î)− ρ) + λm(dm − ρ)) × kmλma(Î) + k(dm + a(Î)− ρ) 2 d(dm + a(Î)− ρ) + λm(dm − ρ) = p(1− ǫ)λ(kmλma+ k(dm + a− ρ) 2) c(δ(dm + a− ρ) + λm(dm − ρ))(d(dm + a− ρ) + λm(dm − ρ)) (3.4) 3.3 Parameter estimation The parameter values used to run the simulation are describe in Table 3.2. The values of the parameters that were found in experimental papers are λ, λmax, k, d, a, λm, δ, δmax, dm, p, and c. The following parameters were not available in the literature but were estimated using the R0 value, the 76 3.3. Parameter estimation disease free equilibrium and the quasi steady state: ρ, dI , λI , βf , and fmax. The rest of the parameters were tested on numerical simulations or were set as zero for simplicity (ǫ, amax, and km). We calculate ρ using the disease free equilibrium. We assume initial values of the total T cell population to be 106 cells/ml, and one third of that to be memory cells [2]. From the equations of T̂ and M̂ (3.3), we get: ρ = M̂ad− λmλ (λm + d)M̂ + dm ρ = T̂ ad− λ(a+ dm − ρ) (λm + d)T̂ + dm making those equal we get: a+ dm − ρ T = λm M ρ =a+ dm − λm T̂ M̂ =0.1 + 0.001 − 0.026 ∗ 2 =0.049 ∴ ρ = 0.049 day−1 (3.5) Let’s proceed to estimate the death rate of the immune response dI . We will let the average lifespan of the immune response be 3 weeks [20]. Then: dI = 1 lifespan = 0.0476 ≈ 0.05 ∴ dI = 0.05 days −1 (3.6) Now, it is acceptable to assume that the immune system level is small the beginning of infection. It could not be zero because of the inate immune system and other cells. We set it to be less than one, ie, Î >> 1 and using (3.3) we have dI >> λI . From that, we fix λI to be 0.005 λI = 0.005 days −1 (3.7) We also want the function fa to have its inflection point at the number of infected T cells at the set point of infection. In order to estimate this value, we use the steady state before treatment. 77 3.4. Numerical simulations From [12], we find that the viral load value at initiation of treatment is 103.5 RNA copies/ml. Then the count of infected cells is given by (3.8). T ∗0 = cV0 p = (1/20000) ∗ (23 ∗ .5) ∗ 103.5 = 1.82 (3.8) ⇒ βf = 1.82 cells/ml Now we fix the immune action to be 5 at this steady state. Using the equation of the steady state of the immune action we look for a fmax esti- mator. 0 =λI + fa(T ∗ 0 )− dI ∗ I ⇒ fa(T ∗ 0 ) = (dI ∗ I − λI) ⇒ fmax =(dI ∗ I − λI) ( β2f T ∗20 + 1 ) =2 (dI ∗ I − λI) ∴ fmax = 2 ∗ 0.245 = 0.49 (3.9) All values are summarized in Table [3.2]. The table also indicates where the parameters were found if they were not calculated in this section. The parameters found in previous work were taken from Rong and Perelson 2009 [17], Conway 2011 [5], Ribeiro et al 2010 [16], Ribeiro et al 2002 [15], and Hazemberg et al 2000 [10]. 3.4 Numerical simulations First of all, we will perform the simulations without Auranofin in the treat- ment. We want to observe a smooth recovery after treatment interruption as well as an early infection peak of the viral load. All the numerical simulations were conducted using MATLAB ode45 function. They assume start of infection at day one. The treatment starts approximately one year later. The interruption is after almost 3.8 years of treatment. The initial conditions used were: T = (2/3) × 106, T ∗ = 0, M = (1/3) ∗ 106, M∗ = 0, V = 10 and I = .1. Figure [3.4] shows the numerical solutions of the system of equations given by (3.1). From Figure [3.4(a)] we can observe that this model presents 3 stages after treatment initiation. There is one fast decay for a few days followed 78 3.4. Numerical simulations Parameter values λ 104 ml−1day−1 [17] λmax 10 5 ml−1day−1 [10] k 2.4 ∗ 10−8 ml×day−1 [17] ǫ Tested on numerical simulation: 0.9592 d 0.1 day−1 [17] a 0.1 day−1 [17] amax 0.14 day −1 (Fixed) λm 0.026 day −1 [15] δ 0.6 day−1 [16] δmax 1 day −1 [16] dm 0.001 day −1 [17] km 0 (for simplicity) ρ 0.049 day−1 p Using R0 > 1: 20000day −1 [5] c 23 day−1 [17] λI 0.005 βf 1.82 cells/ml fmax 0.49 dI 0.05 days −1 Table 3.2: Parameters of the model as we have found them in literature 79 3.4. Numerical simulations 0 500 1000 1500 2000 100 102 104 106 108 1010 Days Vi ra l l oa d (lo g 1 0) ← Starts → Interruption (a) 0 500 1000 1500 2000 0 1 2 3 4 5 6 7 8 9 10 Days Im m un e ac tio n ← Starts → Interruption (b) 0 500 1000 1500 2000 103 104 105 106 107 Days T ce lls (lo g 1 0) ← Starts → Interruption (c) 0 500 1000 1500 2000 10−6 10−4 10−2 100 102 104 106 Days T* ce lls (lo g 1 0) ← Starts → Interruption (d) 0 500 1000 1500 2000 103 104 105 106 107 Days M em or y ce lls (lo g 1 0) ← Starts → Interruption (e) 0 500 1000 1500 2000 10−15 10−10 10−5 100 105 Days M * ce lls (lo g 1 0) ← Starts → Interruption (f) Figure 3.4: Dynamics of the model with parameters ρ = 0.049, ǫ = 0.9592, p = 20000. (a) Viral load V (t). (b) Immune action I(t). (c) Susceptible T-cells T (t). (d) Infected T-cells T ∗(t). (e) Susceptible memory cells M(t). (f) Infected memory cells M∗(t). 80 3.4. Numerical simulations by a medium decay for a few weeks. Then the viral load starts decreasing slowly for the rest of the time. This agrees with experimental observations. On the other hand, the viral load shows a high rebound after treatment interruption (Figure [3.5(a)]). Actually, we can observe from Figure[3.5(b)] that the peak after treatment interruption is even bigger than the one ob- served at early infection. This is not observed in reality. Therefore, we need to reconsider the model construction. 1650 1700 1750 1800 107 108 109 Days Vi ra l l oa d (lo g 1 0) (a) 0 500 1000 1500 2000 107 108 109 Days Vi ra l l oa d (lo g 1 0) ← Starts (b) Figure 3.5: (a) Close look at the viral load after treatment interruption. (b) Comparison of the early infection and the post interruption peak. Parameters of the model ρ = 0.049, ǫ = 0.9592, p = 20000. We consider that one of the reasons the model is not working is because the immune action has no memory. In order to incorporate memory to the immune action, we change the equation for I to be: dI dt = λI + (1− q) ∗ fa(T ∗) + q ∗ fa(M ∗)− dII (3.10) where the parameter q will represent the relative proportion of the immunity driven by the memory cells. We used the rest of the parameters as before and let q be 0.1. Figure [3.6] shows the numerical solution of the new model. There is not much difference to the previous model. Hence we realize that the immune action is not driven by the entire memory cell pool, but only by the activated cells. Thus we change the equation 3.10 by dI dt = λI + (1− q) ∗ fa(T ∗) + q ∗ fa(a ∗M ∗)− dII (3.11) 81 3.4. Numerical simulations As observed in figure 3.9(b), the peak after treatment interruption is smaller than the one at early infection. On the other hand, if we observe the complete simulation for the viral load (Figure3.8(a)), we notice a fourth decrease stage on treatment interruption. This last stage is faster than the one before and this is not observed on experiments. Moreover, this last stage seems to drive the viral load to zero. 82 3.4. Numerical simulations 0 500 1000 1500 2000 100 102 104 106 108 1010 Days Vi ra l l oa d (lo g 1 0) ← Starts → Interruption (a) 0 500 1000 1500 2000 0 1 2 3 4 5 6 7 8 9 10 Days Im m un e ac tio n ← Starts → Interruption (b) 0 500 1000 1500 2000 103 104 105 106 107 Days T ce lls (lo g 1 0) ← Starts → Interruption (c) 0 500 1000 1500 2000 10−6 10−4 10−2 100 102 104 106 Days T* ce lls (lo g 1 0) ← Starts → Interruption (d) 0 500 1000 1500 2000 103 104 105 106 107 Days M em or y ce lls (lo g 1 0) ← Starts → Interruption (e) 0 500 1000 1500 2000 10−15 10−10 10−5 100 105 Days M * ce lls (lo g 1 0) ← Starts → Interruption (f) Figure 3.6: Dynamics of the improved model with parameter q = 0.1. (a) Viral load V (t). (b) Immune action I(t). (c) Susceptible T-cells T (t). (d) Infected T-cells T ∗(t). (e) Susceptible memory cells M(t). (f) Infected memory cells M∗(t). 83 3.4. Numerical simulations 1650 1700 1750 1800 107 108 109 Days Vi ra l l oa d (lo g 1 0) (a) 0 500 1000 1500 2000 107 108 109 Days Vi ra l l oa d (lo g 1 0) ← Starts (b) Figure 3.7: (a) Close look at the viral load after treatment interruption. (b) Comparison of the early infection and the post interruption peak. Parameters of the model q = 0.1 84 3.4. Numerical simulations 0 500 1000 1500 2000 10−2 100 102 104 106 108 1010 Days Vi ra l l oa d (lo g 1 0) ← Starts → Interruption (a) 0 500 1000 1500 2000 0 1 2 3 4 5 6 7 8 9 10 Days Im m un e ac tio n ← Starts → Interruption (b) 0 500 1000 1500 2000 103 104 105 106 107 Days T ce lls (lo g 1 0) ← Starts → Interruption (c) 0 500 1000 1500 2000 10−6 10−4 10−2 100 102 104 106 Days T* ce lls (lo g 1 0) ← Starts → Interruption (d) 0 500 1000 1500 2000 103 104 105 106 107 Days M em or y ce lls (lo g 1 0) ← Starts → Interruption (e) 0 500 1000 1500 2000 10−15 10−10 10−5 100 105 Days M * ce lls (lo g 1 0) ← Starts → Interruption (f) Figure 3.8: Dynamics of the model with parameters ρ = 0.049, ǫ = 0.9592, p = 20000, q = 0.1 and evaluating the function of influence of the memory on the immune action when a memory cell activates. These parameters present a four stage decay of the viral load (a). The last stage is faster than the previous two. The four stage decay with fast ending is not observed on experiments. (b) Immune action I(t). (c) Susceptible T-cells T (t). (d) Infected T-cells T ∗(t). (e) Susceptible memory cells M(t). (f) Infected memory cells M∗(t). 85 3.4. Numerical simulations 1650 1700 1750 1800 107 108 109 Days Vi ra l l oa d (lo g 1 0) (a) 0 500 1000 1500 2000 107 108 109 Days Vi ra l l oa d (lo g 1 0) ← Starts (b) Figure 3.9: (a) Close look at the viral load after treatment interruption. (b) Comparison of the early infection and the post interruption peak. Parameters of the model ρ = 0.049, ǫ = 0.9592, p = 20000, q = 0.1 and evaluating the function of influence of the memory on the immune action when a memory cell activates. From here, we can observe that the post interruption peak is smaller than the early infection one. 86 3.5. Other parameter dynamics 3.5 Other parameter dynamics While I was simulating model3.1 I found some interesting results. The viral load can present oscillations for some values of drug efficacy, maximum of the activation rate and rate of reproduction of memory cells. If the memory cells reproduce faster and become active to a greater level while the drug efficacy is less efficient, we obtain a limit cycle. This is illustrated on Figure 3.10. 0 500 1000 1500 2000 10−2 100 102 104 106 108 1010 Days Vi ra l l oa d (lo g 1 0) ← Starts → Interruption Figure 3.10: Dynamics of the viral load of model (3.1) with parameters ρ = 2 ∗ 0.049, ǫ = 0.98, amax = 0.3. The existence of those parameters suggest a Hopf bifurcation. We don’t go further in this result. It is left for possible future analysis of the model. 3.6 Summary This model was not able to reproduce the viral load rebound when treat- ment without Auranofin is interrupted (Figure 3.1). This was the first step on the modeling, even before trying to test different possible mechanism of Auranofin. From here it is not possible to prove something about the differences in dynamics with and without Auranofin therapy. From the numerical solutions of the model it is observable that the acti- vation of the memory cell is affecting the fast return to the viral set point. Maybe looking for new values of the activation rate will improve the results. Another weak part of this modeling is the generalized way to model the immune action. The immune system has a complex reaction in the 87 3.6. Summary presence of external agents, such as virus. It is influenced by more than the T cell cycle. It has innate and adaptive responses that were not included in this model and that are related on the new generation of T cells and on the depletion of infected cells and virus. It will be needed to improve this compartment to test the hypothesis that: first the immune action is different at normal ART interruption than before infection; and second that when Auranofin is added, the solution at time interruption has similarities to the one at the start of infection. 88 Bibliography [1] J.B. Angel, J.P. Routy, C. Tremblay, D. Ayers, R. Woods, J. Singer, N. Bernard, C. Kovacs, F. Smaill, S. Gurunathan, et al. A randomized controlled trial of HIV therapeutic vaccination using ALVAC with or without Remune. AIDS, 25(6):731, 2011. [2] T.P. Arstila, A. Casrouge, V. Baron, J. Even, J. Kanellopoulos, and P. Kourilsky. A direct estimate of the human αβ T cell receptor diver- sity. Science, 286(5441):958–961, 1999. [3] H.Y. Chen, M. Di Mascio, A.S. Perelson, D.D. Ho, and L. Zhang. Determination of virus burst size in vivo using a single-cycle SIV in rhesus macaques. Proceedings of the National Academy of Sciences, 104(48):19079, 2007. [4] T.W. Chun, L. Carruth, D. Finzi, X. Shen, J.A. DiGiuseppe, H. Taylor, M. Hermankova, K. Chadwick, J. Margolick, T.C. Quinn, et al. Quan- tification of latent tissue reservoirs and total body viral load in HIV-1 infection. Nature, 387(6629):183–188, 1997. [5] J.M. Conway and D. Coombs. A stochastic model of latently infected cell reactivation and viral blip generation in treated HIV patients. PLoS Computational Biology, 7(4):e1002033, 2011. [6] O. Diekmann, JAP Heesterbeek, and JAJ Metz. On the definition and the computation of the basic reproduction ratio R0 in models for infec- tious diseases in heterogeneous populations. Journal of mathematical biology, 28(4):365–382, 1990. [7] CW Gardiner. Handbook of stochastic methods(for physics, chemistry and the natural sciences). Springer series in synergetics. [8] D.T. Gillespie. Exact stochastic simulation of coupled chemical reac- tions. The journal of physical chemistry, 81(25):2340–2361, 1977. [9] R. Haberman. Applied partial differential equations with fourier series and boundary value problems. Pearson Education, 2004. 89 Bibliography [10] M.D. Hazenberg, J.W.T.C. Stuart, S.A. Otto, J.C.C. Borleffs, C.A.B. Boucher, R.J. de Boer, F. Miedema, and D. Hamann. T-cell division in human immunodeficiency virus (HIV)-1 infection is mainly due to immune activation: a longitudinal analysis in patients before and during highly active antiretroviral therapy (HAART). Blood, 95(1):249–255, 2000. [11] JM Heffernan, RJ Smith, and LM Wahl. Perspectives on the basic reproductive ratio. Journal of the Royal Society Interface, 2(4):281– 293, 2005. [12] M.M. Kitahata, S.J. Gange, A.G. Abraham, B. Merriman, M.S. Saag, A.C. Justice, R.S. Hogg, S.G. Deeks, J.J. Eron, J.T. Brooks, et al. Effect of early versus deferred antiretroviral therapy for HIV on survival. New England Journal of Medicine, 360(18):1815–1826, 2009. [13] M.G. Lewis, S. DaFonseca, N. Chomont, A.T. Palamara, M. Tardugno, A. Mai, M. Collins, W.L. Wagner, J. Yalley-Ogunro, J. Greenhouse, et al. Gold drug auranofin restricts the viral reservoir in the monkey AIDS model and induces containment of viral load following ART sus- pension. AIDS, 25(11):1347, 2011. [14] M.A. Nowak, A.L. Lloyd, G.M. Vasquez, T.A. Wiltrout, L.M. Wahl, N. Bischofberger, J. Williams, A. Kinter, A.S. Fauci, V.M. Hirsch, et al. Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection. Journal of virology, 71(10):7518– 7525, 1997. [15] R.M. Ribeiro, H. Mohri, D.D. Ho, and A.S. Perelson. In vivo dynamics of T cell activation, proliferation, and death in HIV-1 infection: Why are CD4+ but not CD8+ T cells depleted? Proceedings of the National Academy of Sciences, 99(24):15572, 2002. [16] R.M. Ribeiro, L. Qin, L.L. Chavez, D. Li, S.G. Self, and A.S. Perel- son. Estimation of the initial viral growth rate and basic reproductive number during acute HIV-1 infection. Journal of Virology, 84(12):6096– 6102, 2010. [17] L. Rong and A.S. Perelson. Modeling latently infected cell activation: viral and latent reservoir persistence, and viral blips in HIV-infected pa- tients on potent therapy. PLoS computational biology, 5(10):e1000533, 2009. 90 Bibliography [18] A. Savarino. Eradication trials in SIVmac251-infected macaques. In Fifth International Workshop on HIV Persistence during Therapy, vol- ume 7, pages 3–4. Global Antiviral Journal, 2011. [19] I.L. Shytaj, S. Norelli, B. Chirullo, A. Della Corte, M. Collins, J. Yalley- Ogunro, J. Greenhouse, N. Iraci, E.P. Acosta, M.L. Barreca, et al. A highly intensified ART regimen induces long-term viral suppression and restriction of the viral reservoir in a simian AIDS model. PLoS Pathogens, 8(6):e1002774, 2012. [20] L.M. Sompayrac. How the Immune System Works. Wiley-Blackwell, 2008. [21] H.M. Taylor and S. Karlin. An introduction to stochastic modeling. Academic Pr, 1984. [22] A. Tveito and R. Winther. Introduction to partial differential equations: a computational approach, volume 29. Springer Verlag, 1998. [23] UNADIS. Global report fact sheet, 2010. [24] N.K. Vaidya, R.M. Ribeiro, C.J. Miller, and A.S. Perelson. Viral dy- namics during primary simian immunodeficiency virus infection: Effect of time-dependent virus infectivity. Journal of virology, 84(9):4302– 4310, 2010. [25] S. Wolfensohn and M. Lloyd. Handbook of laboratory animal manage- ment and welfare. Wiley Online Library, 2003. 91