Abstract

Motivation: Since the identification of human immunodeficiency virus (HIV) over twenty years ago, many mathematical models of HIV dynamics have been proposed. The purpose of this study was to evaluate intracellular and intercellular scale HIV models that best described the dynamics of viral and cell titers of a person, where parameters were determined using typically available patient data. In this case, ‘best’ was defined as the model most capable of describing experimental patient data and was determined by Bayesian-based model discrimination analysis and the ability to provide realistic results.

Results: Twenty models of HIV-1 viral dynamics were initially evaluated to determine whether parameters could be obtained from readily available clinical data from established HIV-1 patients with stable disease. Based on this analysis, three models were chosen for further examination and comparison. Parameters were estimated using experimental data from a cohort of 338 people monitored for up to 2484 days. The models were evaluated using a Bayesian technique to determine which model was most probable. The model ultimately selected as most probable was overwhelmingly favored relative to the remaining two models, and it accounted for uninfected cells, infected cells and cytotoxic T lymphocyte dynamics. The authors developed a fourth model for comparison purposes by combining the features of the original three models. Parameters were estimated for the new model and the statistical analysis was repeated for all four models. The model that was initially favored was selected again upon model discrimination analysis.

Contact:srivasta@engr.uconn.edu

INTRODUCTION

Since it was initially identified, significant progress has been made in treating human immunodeficiency virus (HIV) due to improved treatment regimens, earlier identification and increased public awareness (Fauci, 2003; Pomerantz and Horn, 2003). However, there is still disagreement and confusion concerning the disease on several fronts. For example, many researchers are trying to elucidate which cells can be infected and which contribute a significant number of viral particles to infection (Brenchley et al., 2004). Others are attempting to determine whether intermittent viral rebound arises from hidden viral reservoirs, from production of resistant viral strains or from drugs that are not completely effective (Cohen Stuart et al., 2001; Herz et al., 1996). The determination of exactly how HIV causes acquired immune deficiency syndrome (AIDS) is also still open to debate (Stebbing et al., 2004; Stevenson, 2003). Traditionally experimental data has been gathered independently of creating a kinetic model. When models are developed, they are frequently based on experimental data from only a few patients. By comparing multiple models to a large experimental dataset, it may be possible to characterize this disease more precisely. Kinetic modeling is an important tool for understanding the course of infection in an afflicted individual, as well as for optimizing treatment strategies and characterizing antivirals. The work presented here is a contribution to the global effort to provide better treatment to those infected with HIV.

A large number of deterministic and stochastic intercellular and intracellular models have been proposed to understand the dynamics of HIV (Bonhoeffer et al., 1997a,b; Nowak et al., 1997; Perelson et al., 1996; Phillips, 1996; Reddy and Yin, 1999; Ribeiro and Bonhoeffer et al., 1999; Tan and Wu, 1998; Tuckwell and Le Corfec, 1998). In addition, a large amount of experimental data has been gathered through clinical trials such as the AIDS Clinical Trials Groups of the National Institutes of Health. Although many kinetic models have been developed to simulate viral infection, little effort has been made to identify, quantitatively or qualitatively, which model can best predict the course of the disease in a specific person. One study compared three deterministic models of human immunodeficiency type 1 (HIV-1) infection, monitoring the differences in different cells and ascertaining the location of the cell populations (Snedecor, 2003). While the research yielded beneficial results on the quality of various models, no quantitative evaluation was performed among the models. The current work attempts to determine which model best describes HIV-1 dynamics based on available experimental data.

SYSTEMS AND METHODS

Parameters were estimated based on data provided by the VA hospital in West Haven, Connecticut, for a cohort of 338 people monitored for up to 2484 days (Hamilton et al., 1992; O'Brien et al., 1996). This data set is suitable for the current study because of the large number of patients included and the length of time for which they were monitored. To qualify for the VA study, patients had to have symptoms or signs of HIV infection as well as a CD4+ titer between 200 and 500 cells/mm3. In addition, patients with AIDS-defining illnesses were screened out of the study, ensuring that no one had progressed to AIDS. Patients were generally stable in their disease, as those classified as unstable were omitted from observation. Additional characteristics of the patients involved are summarized by Hamilton et al. (1992). The patients were treated with AZT either immediately upon entering the study or when their CD4+ counts dropped <200 cells/mm3. The experimental data includes virus, CD4+ cell and white blood cell (WBC) titers, the percentage of WBCs that are lymphocytes and the number of days from baseline that the measurements were taken. Baseline was defined as the time point when the first reading was taken for a given study participant.

Twenty intercellular models, ranging from simple to complex, were initially examined for the current study (Bonhoeffer et al., 1997a,b; De Boer and Perelson, 1998; Essunger and Perelson, 1994; Funk et al., 2001; Nelson et al., 2001; Perelson, 2002; Perelson et al., 1996; Stafford et al., 2000; Wein et al., 1998; Wodarz and Nowak, 2002) with several papers proposing more than one model. Table 1 lists the features of HIV infection that were accounted for in each model. The initial criterion for model selection was the ability to estimate parameters from typical patient data. Models that compensate for drug treatment and the corresponding efficacy of reverse transcriptase and protease inhibitors provide useful detail; however, for a number of reasons for this study, only data collected from untreated patients were used. By studying untreated patients, it was possible to focus on the details of HIV dynamics without external interference from medication. Another reason for not evaluating drug-dependent models was that they often monitored variables such as non-infectious virus particles. While it is true that these particles are created during treatment, there is no simple way to determine the number of defective particles from available experimental data. It would also be difficult to compare the treatment regimen used in the VA study to present courses of therapy because monotherapy has been replaced by highly active antiretroviral therapy (HAART) and drug cocktails (Flint et al., 2004; Goldsby et al., 2003). Further, AZT resistance develops quickly thus limiting its use (Flint et al., 2004; Goldsby et al., 2003) and various HAART treatments rarely use the same combinations of reverse transcriptase inhibitors (RTI) and protease inhibitors (PI). As a result, drug-dependent models were excluded. Other models were not included for analysis because they monitored variables that could not be quantified without significant effort or genome sequencing. These variables included active or latent infection, infection with different viral strains and virus particles susceptible or resistant to therapy.

Given that infected and uninfected cells are subpopulations of the total CD4+ count, it was necessary to approximate what proportion of a patient's CD4+ count could be attributed to infected cells. Many reports in the literature attempt to determine what fraction of the CD4+ population is infected with HIV-1, with results ranging from one infected cell per 150 000 CD4+ cells (6.7× 10−6%) to 69 infected cells per 100 CD4+ cells (69%) (Bagasra et al., 1992,1993; Bieniasz et al., 1993; Brinchmann et al., 1991; Chun et al., 1997; Harper et al., 1986; Ho et al., 1989; Hsia and Spector, 1991; Macatonia et al., 1990; Poznansky et al., 1991; Re et al., 1994; Schnittman et al., 1889,1990; Simmonds et al., 1990). In the present analysis an intermediate value of 2% infected CD4+ cells was assumed, which may be high for newly infected people and low for someone nearing fully-developed AIDS.

A further assumption was necessary to determine what percentage of WBCs were cytotoxic T lymphocytes (CTL), as shown in Figure 1. CTLs are the mature form of cytotoxic T (Tc) cells that mediate killing of infected cells. For each experimental measurement of WBCs, there was a corresponding percentage of lymphocytes. Between 48 and 67% of lymphocytes are T cells (Goldsby et al., 2003). For this study 58% of lymphocytes were assumed to be T cells. Further, ∼30% of T cells have CD8+ on their surface (Goldsby et al., 2003), indicative of a Tc cell; however, to the knowledge of the authors, there is no known relation between the number of Tc cells and the number of active CTL effector cells at any given point in time. To simulate a minimal and a complete immune response, CTL cells were assumed to be either 1 or 100% of the Tc population in two different scenarios.

Three literature models of cellular and viral interactions were chosen for detailed analysis. Table 2 summarizes the parameter values used for each model. Parameters having the same function in different models are referred to by the same name. For example, k1 refers to the rate constant characterizing uninfected cell production in all models. The value of the parameter, however, will differ from model to model based on what is and what is not accounted for. By presenting the parameters in this way, the authors felt a better reference for comparison was established.

The first model, referred to as the basic model, monitors uninfected cells (X), infected cells (Y) and virus (V) 

(1)
\[\hbox{ d }X/\hbox{ d }t={k}_{1}-{k}_{2}X-{k}_{3}XV\]
 
(2)
\[\hbox{ d }Y/\hbox{ d }t={k}_{3}XV-{k}_{4}Y\]
 
(3)
\[\hbox{ d }V/\hbox{ d }t={k}_{5}Y-{k}_{6}V,\]
where k1 is the supply rate of uninfected cells by the thymus, k2X is the death rate of uninfected cells, k3XV is the rate of infection, k4Y is the death rate of infected cells, k5Y is the rate of virus production by infected cells and k6V is the clearance rate of virus (Bonhoeffer et al., 1997b).

The second model is a simpler version of the basic model and will be referred to as the simplified model. The number of uninfected cells is assumed to be constant, 

(4)
\[\hbox{ d }Y/\hbox{ d }t={k}_{3}XV-{k}_{4}Y\]
 
(5)
\[\hbox{ d }V/\hbox{ d }t={k}_{4}{k}_{7}Y-{k}_{6}V,\]
where k7 is the number of new virus particles produced per infected cell (Perelson et al., 1996).

In the third model, or the immune model, a constant number of virions is assumed. The model also includes the response of the immune system (Z) 

(6)
\[\hbox{ d }X/\hbox{ d }t={k}_{1}-{k}_{2}X-{k}_{8}XY\]
 
(7)
\[\hbox{ d }Y/\hbox{ d }t={k}_{8}XY-{k}_{4}Y-{k}_{9}YZ\]
 
(8)
\[\hbox{ d }Z/\hbox{ d }t={k}_{10}Y-{k}_{11}Z,\]
where Z is a CTL, k8XY is the rate of infection, k9YZ is the rate that CTLs kill infected cells, k10Y is the rate of CTL stimulation by infected cells and k11Z is the death rate of CTLs (Bonhoeffer et al., 1997a).

A fourth model, referred to as the composite model, was developed during this study to incorporate aspects of the first three models presented. This was done with the ambition of creating a more comprehensive model that simultaneously accounted directly for the immune system response and for the viral dynamics. 

(9)
\[\hbox{ d }X/\hbox{ d }t={k}_{1}-{k}_{2}X-{k}_{3}XV\]
 
(10)
\[\hbox{ d }Y/\hbox{ d }t={k}_{3}XV-{k}_{4}Y-{k}_{9}YZ\]
 
(11)
\[\hbox{ d }V/\hbox{ d }t={k}_{5}Y-{k}_{6}V\]
 
(12)
\[\hbox{ d }Z/\hbox{ d }t={k}_{10}Y-{k}_{11}Z.\]

One should note that CTLs do not remove free virus particles; they only remove infected cells as indicated in Equation (10) by the rate k9YZ. The actions of CTLs are generally only against cells displaying viral antigen on their surface, not virus particles themselves (Goldsby et al., 2003).

To provide a common frame of reference, all parameters for each model were estimated from the VA data. Parameter estimation was carried out using non-linear least squares regression analysis. Minimization of the objective function was accomplished by using the Global Optimization package v.4.2 from Loehle Enterprises. Once values of rate constants were estimated, the posterior probability share for each model was calculated by using the Bayesian technique described in the Algorithm section.

ALGORITHM

Model discrimination was carried out using the method of Stewart et al. (1998). For convenience, a brief description of the method is provided here. Measurements of experimental data can be expressed mathematically as 

(13)
\[{y}_{iu}={f}_{ji}\left({\underset{\bar}{\xi }}_{u},{\underset{\bar}{\theta }}_{j}\right)+{\epsilon }_{iu},\]
where the observation y is a function f of the independent variable vector ξ and the parameter vector θ plus some error ε. There are i = 1, …, m responses for each of the u = 1, …, n events of j = 1, …, J models. To account for variation in the precision of observations with experimental conditions, Equation (13) may be weighted by the square root of the inverse variance of yiu (Stewart et al., 1998). 
(14)
\[{Y}_{iu}={F}_{ji}\left({\underset{\bar}{\xi }}_{u},{\underset{\bar}{\theta }}_{j}\right)+{E}_{iu}.\]
The variables Y, F and E are the weighted values of y, f and ε, respectively.

The experimental dataset provided for discrimination analysis was highly irregular, as shown in Table 3, meaning that many responses were missing and there was no identifiable pattern of data collection. A ‘response’ was the reading of one patient for a specified event. An ‘event’, for the purposes of this study, was defined as any day on which a measurement was made relative to the baseline. Although methods exist for comparing models with highly irregular data, also described by Stewart et al. (1998) the methods require a known covariance matrix. The covariance matrix was not calculable from the data obtained for this study, so the only methods applicable are those for regular data. To force the data into a regular, rectangular form, as shown in Table 4, responses that had a missing variable, such as CD4+ or virus titer, were not included in the analysis. Events that had only one response were ignored, because calculating the variance requires a minimum of two responses.

For full data with unknown covariance, the posterior probability for a given model is calculated with a Bayesian- based technique. The probability of the j-th model, given a set of observations, is 

(15)
\[p\left({M}_{j}\right|\underset{\bar}{\underset{\bar}{Y}})\propto p({M}_{j}){2}^{-{p}_{j}/2}{\left|{\underset{\bar}{\underset{\bar}{v}}}_{j}\right|}^{-{v}_{e}/2}\]
where Mj is model j;
\(\underset{\bar}{\underset{\bar}{Y}}\)
is the matrix of weighted observations; pj is the number of parameters in model j;
\({\underset{\bar}{\underset{\bar}{v}}}_{j}\)
is the matrix of products of the deviation of Yiu from the predicted value for model j, evaluated at the maximum likelihood θj; and ve is the number of degrees of freedom. p(Mj) is the prior probability of model j, and is assumed to be equal for all models unless there is reason to give a higher probability to a particular model. The ik-th element of
\({\underset{\bar}{\underset{\bar}{v}}}_{j}\)
is calculated according to (Stewart et al., 1992, 1998; Stewart and Sørenson, 1981). 
(16)
\[{v}_{ik}\left({\underset{\bar}{\theta }}_{j}\right)={\displaystyle \sum _{u=1}^{n}}\left[{Y}_{iu}-{F}_{ji}\left({\underset{\bar}{\xi }}_{u},{\underset{\bar}{\theta }}_{j}\right)\right]\left[{Y}_{ku}-{F}_{jk}\left({\underset{\bar}{\xi }}_{u},{\underset{\bar}{\theta }}_{j}\right)\right],\]

Once the posterior probability is calculated for each model, normalization to a sum of one gives the posterior probability share, π, of model Mj: (Stewart et al., 1996, 1998). 

(17)
\[\pi \left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)=\frac{p\left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)}{{\sum }_{k}p\left({M}_{k}|\underset{\bar}{\underset{\bar}{Y}}\right)}.\]
The model with the largest value of π is considered the most probable. Although the methods detailed above will allow for the selection of the best model among a set of potential models, the model with the highest probability share is not necessarily a good model. If all models characterize the system poorly, calculating π will simply choose the least poor model (Stewart et al., 1998). A good model must be able to use clinical data to accurately describe viral dynamics.

RESULTS

After 17 of the original 20 models were excluded from consideration, non-linear least squares regression analysis was used to determine parameter values for the remaining models. Parameter values, units and definitions are displayed in Table 2. Figures 25 show the predictions of each model versus average experimental data with respect to the four variables: uninfected cells, infected cells, virus particles and CTL response. Results presented in all figures and tables were based on the assumption that CTLs comprised 1% of Tc cells because it is unlikely that every Tc cell would be an effector at any point in time.

The three models that accounted for uninfected cells each reached a steady state shortly after baseline. A noticeable difference among the models was the steady state value, where the basic model predicted 180 cells/mm3 after 150 days, and the immune model predicted the largest value of 267 cells/mm3 after 100 days (Fig. 2). To be included in the VA Cooperative Study Group, patients had to have CD4+ concentrations ranging from 200–500 cells/mm3, which indicates that they were far along in the course of infection. The basic model was slightly different from the immune and composite models in that it predicted an initial decline and resurgence of uninfected cells, which is accurate for initial infection, but not for individuals that have had the disease for some time. The cause of the transient behavior may be the result of an elevated initial condition, which was calculated by averaging all baseline readings. On day zero 167 individuals had CD4+ readings, but few subsequent days included more than 10 readings on average. This fact may indicate that as the number of days from baseline increased, those with higher virus titers discontinued the study. Fewer recordings per day can also be attributed to not having all patients return on the same day from baseline for follow-up; many days had only one response and were thus not included in model discrimination.

All four models accounted for the change in the number of infected cells but their trajectories were very different. In Figure 3, the basic model generated an infected cell trajectory that had an initial abrupt increase followed shortly thereafter by a sharp decrease and maintenance of a steady state. This result correlated with the decrease and increase of uninfected cells, shown in Figure 2. As HIV proliferated, the pool of healthy cells decreased and the number of infected cells increased. Both the immune and composite models predicted a very low value of infected cells beginning at baseline and continuing throughout the course of infection. Because the number of infected cells went to zero after a short time, the uninfected cells reached steady state in both models. In the immune model k8XY dropped out of Equation (6) and k1 balanced k2X. Once infected cells had been eliminated in the composite model, the number of virus particles declined because of k6V of Equation (11) which then caused a decrease in the rate of change of X as represented by Equation (9). When virus was not present, k1 and k2X of Equation (9) balanced each other to achieve steady state. The simplified model trajectory decreased monotonically and appeared to best fit the infected cell data. From Equation (4) it is apparent that the decrease in infected cells is due to the fact that more infected cells die than are produced after each round of infection.

The basic, simplified and composite models had very different trajectories for viral titer. The simplified model most accurately fits the viral experimental data upon visual inspection. The slow decrease mimicked the slow decline in the number of infected cells. The basic model predicted a dramatic increase in concentration followed by decline and achievement of steady state. The number of virions increased dramatically because the number of infected cells also rose, yielding more successful infections. It is easy to see from Figure 4 that the steady state of the basic model was significantly higher than the other model predictions (note the secondary axis on the right), and not very representative of the experimental data. A large production rate of new virions, k5, and low clearance rate, k6, most likely caused this trend. The composite model decreased from the initial condition until the eventual prediction of viral clearance. Steady-state analysis, shown in Table 5, indicated that although the simplified model also eventually predicted viral clearance, this was not reached until after approximately 30 000 days.

The immune system response was a feature of only the immune and composite models. All results presented in Figure 5 are for the assumption that CTLs were 1% of Tc cells. However, it was interesting to note that, except for CTLs themselves, results were essentially the same when CTLs were assumed to be 100% of the Tc population (data not shown). This behavior suggested the system was robust with respect to a wide range of CTL values. The composite model appeared to reflect the data well, with a slight decrease in titer over time. This effect could be explained by recent evidence that CD4 on the surface of CD8+ cells causes infection of those cells by HIV (Kitchen et al., 2004), thus decreasing the population through the course of infection. The immune model also predicted values close to experimental titers, but the decrease in concentration of CTL cells was more dramatic and larger at the beginning of the experiment. The effect was an increasing discrepancy between the two models over the course of time. Consistent decline was seen because once infected cells were removed from the system, the only term present in Equation (8) was death of CTLs.

A summary of the steady-state values for all four models is presented in Table 5. A number of the steady states are reached beyond the life expectancy of a patient, and others are achieved beyond the life expectancy of a healthy human being. Although the simplified model predicts a trivial state of viral and infected cell clearance, this occurs after a very long time. On the other hand, though the immune and composite models take a long time to reach their steady-state values, it is easy to see that after a short time the predicted values are already approaching zero in Figures 3 and 4, which is known to be unrealistic.

The immune model was overwhelmingly favored as most likely, based on model discrimination analysis. Table 6 presents the probability share of the four models considered in the present analysis. While the results of Table 1 are for the entire range of data, the immune model still had a probability of ∼1.00 when all models were fit to the first 25, 50, 100, 250 and 500 days (data not shown). It was peculiar that with the parameters used in this study, the immune model predicted clearance of infected cells, which is not a biological reality. The immune model may have had the largest π because it had a variance several orders of magnitude lower than that of all other models, and therefore a larger determinant in Equation (15). One can see from this equation that having a large variance will penalize the model. The end result was a probability share many orders of magnitude larger than that of the other three models considered.

If the immune model was ignored and the probability share was recalculated over the entire course of the infection the simplified model was the most favored. If the immune model was disregarded and models were again fitted and evaluated for the first 25, 50, 100, 250 or 500 days of experimental data, the simplified model still retained the highest probability share (data not shown). This was expected because the composite model predicted clearance of virus and infected cells, both of which are incorrect in vivo. In addition, the basic model poorly fits the data in the first 100 days of all variables considered. As a result, the simplified model best fits the data over the range considered, yielding a probability share of 1.00. Although the simplified model eventually predicted the elimination of the virus, this did not occur within the time frame of the data presented.

DISCUSSION

Three steps were employed to determine the best model of HIV-1 intercellular kinetics: estimating parameters from clinical data, using model discrimination analysis to mathematically determine the best model and utilizing steady-state analysis to conclude if cell titer and viral load predictions were logical and sensible. A prediction was considered logical and sensible if it agreed with clinical data. Twenty deterministic models of HIV-1 intercellular dynamics were narrowed to three based on the variables they accounted for and the available experimental data. Upon using model discrimination analysis, the immune model had the largest probability share of the three models evaluated, implying it was the most favorable. However, steady-state analysis made it clear that this result was unrealistic because, as shown in Figure 3, infected cells were falsely predicted to be cleared from the body. Interestingly, although the immune model was very similar to the composite model, and analogous parameters were assigned values within one order of magnitude to each other, the probability shares were very different. Both models incorporated immune-mediated cell killing in addition to infected and uninfected cells, and both erroneously predicted clearance of infected cells. A significant difference was that the immune model explicitly assumed a constant number of virions. However, if at steady state Equation (11) is set equal to zero, as in the immune model, Y = k6V/k5. Substituting this into the k8XY of Equations (6) and (7), one can see that both models are equivalent at steady state. The composite model's additional prediction of viral clearance was a possible reason why it had a lower probability share than the immune model. While it is true that being exposed to HIV will not definitely cause infection, once an individual has contracted the disease, it is not eliminated from the body.

Although it was initially surprising that the immune model was the most favored, this result can be explained mathematically. It is important to realize the immune model only accounted for infected/uninfected CD4+ T cells and CTLs. It did not explicitly account for viral dynamics and was not fitted to viral data; rather, it assumed a constant level of virus. Because of this assumption, a steady state was reached where virus and uninfected cells coexisted without the presence of infected cells. The constant viral load assumption may have been an oversimplification because although a steady state is reached within a single person, a uniform steady-state concentration is not achieved within a population (Bonhoeffer et al., 2003). More suitable data may have avoided such assumptions, and may have given different probability shares for the models considered.

Caution should be exercised if applying the immune model to the entire course of infection because the viral titer is not always at steady state. In the initial stage of the disease, the amount of virus has a precipitous increase and then declines to a steady-state value (Stafford et al., 2000). At the later stages immediately preceding and after the on-set of AIDS, virus titer increases, causing a decline in the number of healthy CD4+ T cells. Assuming that the immune model is most likely, it would be applicable to the experimental data because many of the patients in this study had HIV for a long enough time to have reached a steady viral titer. However, the number of infected cells would not be accurately accounted for, as the model predicts eventual clearance.

Even though the immune model was determined to be the most favorable under the conditions studied based on the Bayesian analysis, the model may in fact not be the best in general. The immune model ignored changes in virus concentration, and poorly predicted the number of infected cells (Fig. 3). Perhaps with the appropriate data this model would forecast titers more precisely. However, given the available data, the best model to use for clinical studies may be the simplified model because it only had infected cells and virus in measured quantities, and did not predict clearance of either quantity in the time frame of interest. In addition, this model was assigned the highest probability share when the immune model was disregarded.

It is apparent that as a model becomes more complex it is more difficult to accurately account for everything occurring in the system, which is evidenced by the fact that the most complex model had a low probability share. The composite model was an extension of the immune model because it explicitly accounted for viral dynamics. Yet when compared, the composite model had a significantly lower probability share then the immune model. Furthermore, the parameters estimated for the composite model not only predicted clearance of infected cells, but also predicted clearance of virus particles. Similarly, the basic model was an extension of the simplified model, with the basic model explicitly accounting for uninfected cells. Yet the simplified model was considered the better candidate based on the discrimination analysis, although the difference in the probability shares was not as dramatic.

Supplementary research can apply the same methods of model discrimination to systems that include terms, such as drug efficacy, that account for changes in the viral network as a result of drug treatment. Such research may provide a definitive decision of which model is most acceptable.

Although the immune model had the highest probability share, the simplified model appears to be the most accurate for the present. In general, the models considered were most applicable for the intermediate stages of HIV-1 infection, as that was the range of data used in this study. During the analysis it was necessary to make several assumptions, including how many Tc cells were CTLs. Better assays for measuring the number of CTLs in HIV patients would allow a more comprehensive analysis of the models. However, the feasibility of such assays remains to be determined. Appropriate CTL data may give a more accurate representation of titers over time and may lead to the definitive decision that the immune model is the best.

APPENDIX

  • ε = Error

  • E = Weighted error

  • fji = Prediction of response i by model j

  • Fji = Weighted prediction of response i by model j

  • i, k = Response number

  • j = Model number

  • M = Model

  • p = Probability

  • pj = Number of parameters in model j

  • π = Posterior probability share

  • θj = Vector of parameters for model j

  • u = Event number

  • ve = Degrees of freedom

  • \({\underset{\bar}{\underset{\bar}{v}}}_{j}\)
    = Maximum likelihood value of v(θj)

  • V = Virus particle

  • ξu = Vector of independent variables for event u

  • X = Uninfected CD4+ T cell

  • yiu = Observation of response i in event u

  • Yiu = Weighted observation of response i for event u

  • Y = Infected CD4+ T cell

  • Z = Cytotoxic T lymphocyte (CTL) effector

Table 1

Characterization of evaluated models

Variable or feature monitored Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997b) Bonhoeffer et al. (1997b) De Boer and Perelson (1998) De Boer and Perelson (1998) Essunger and Perelson (1994) Essunger and Perelson (1994) Funk et al. (2001) Nelson et al. (2001) Nelson et al. (2001) Nelson et al. (2001) Perelson et al. (1996) Perelson et al. (1996) Perelson (2002) Stafford et al. (2000) Wein et al. (1998) Wodarz and Nowak (2002) Wodarz and Nowak (2002) Composite 
Uninfected CD4+ cells (no further distinction)           
Uninfected quiescent T cells                     
Uninfected active T cells                  
Uninfected memory T cells                    
Uninfected virgin T cells                    
Uninfected long-lived cells (such as macrophages)                     
Infected cells (no further distinction)        
Infected active T cells                   
Latently infected cells                     
Defectively infected cells                     
Persistently infected cells                     
Transiently infected virgin T cells                     
Transiently infected memory T cells                     
Infected memory T cells                    
CD4+ cells infected by different virus strains                    
Long-lived cells infected by different virus strains                     
Variable death rate of infected cells                     
Virus (no distinction of sensitivity or strain)            
Infectious or wild-type virus               
Noninfectious, drug sensitive, or mutant virus               
Multiple, unspecified virus strains                     
CTL response                    
CTL precursors                     
CTL effectors                   
Immune response against various strains                     
Reverse transcriptase inhibitor efficacy                 
Protease inhibitor efficacy                 
Time delay                    
Antigen                     
Probability or rate of mutation                    
Variable or feature monitored Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997a) Bonhoeffer et al. (1997b) Bonhoeffer et al. (1997b) De Boer and Perelson (1998) De Boer and Perelson (1998) Essunger and Perelson (1994) Essunger and Perelson (1994) Funk et al. (2001) Nelson et al. (2001) Nelson et al. (2001) Nelson et al. (2001) Perelson et al. (1996) Perelson et al. (1996) Perelson (2002) Stafford et al. (2000) Wein et al. (1998) Wodarz and Nowak (2002) Wodarz and Nowak (2002) Composite 
Uninfected CD4+ cells (no further distinction)           
Uninfected quiescent T cells                     
Uninfected active T cells                  
Uninfected memory T cells                    
Uninfected virgin T cells                    
Uninfected long-lived cells (such as macrophages)                     
Infected cells (no further distinction)        
Infected active T cells                   
Latently infected cells                     
Defectively infected cells                     
Persistently infected cells                     
Transiently infected virgin T cells                     
Transiently infected memory T cells                     
Infected memory T cells                    
CD4+ cells infected by different virus strains                    
Long-lived cells infected by different virus strains                     
Variable death rate of infected cells                     
Virus (no distinction of sensitivity or strain)            
Infectious or wild-type virus               
Noninfectious, drug sensitive, or mutant virus               
Multiple, unspecified virus strains                     
CTL response                    
CTL precursors                     
CTL effectors                   
Immune response against various strains                     
Reverse transcriptase inhibitor efficacy                 
Protease inhibitor efficacy                 
Time delay                    
Antigen                     
Probability or rate of mutation                    

Each row is either a variable monitored by one or more models or a feature that sets models apart from each other. Each column represents one model. Note that some papers presented more than one model. Boxes with an ‘X’ indicate that the model included the particular variable. With the exception of the ‘composite’ model, the name at the top of each column is the reference publication for the model, as listed in the References section. Details regarding the composite model may be found in System and Methods section of this paper.

Fig. 1

Delineation of the number of white blood cells that are Tc, lymphocytes. Lowest tier values come from assumptions made by the authors. All other values come from Goldsby et al. (2003).

Fig. 1

Delineation of the number of white blood cells that are Tc, lymphocytes. Lowest tier values come from assumptions made by the authors. All other values come from Goldsby et al. (2003).

Table 2

Parameter values used in each model, as determined through non-linear least squares analysis

Parameter Parameter description Units Basic model Simplified model Immune model Composite model 
k1 Supply of uninfected cells by the body cells/(L*day) 1 × 107 — 2 × 107 1 × 107 
k2 Uninfected cell death day−1 0.05 — 0.075 0.05 
k3 Infection of uninfected cells L/(cells*day) 5 × 10−10 5 × 10−10 — 5 × 10−10 
k4 Infected cell death day−1 0.4 0.025 0.01 0.5 
k5 Production of new virus from infected cells day−1 40 — — 
k6 Clearance of virus day−1 2.95 — 0.008 
k7 Number of viruses produced by an infected cell virus/cell — 16.75 — — 
k8 Infection of uninfected cells, independent of virus L/(cells*day) — — 5 × 10−10 — 
k9 Clearance of infected cells by CTLs L/(cells*day) — — 0.01 0.001 
k10 Stimulation of CTLs by infected cells day−1 — — 10 
k11 CTL death day−1 — — 0.001 0.0001 
Parameter Parameter description Units Basic model Simplified model Immune model Composite model 
k1 Supply of uninfected cells by the body cells/(L*day) 1 × 107 — 2 × 107 1 × 107 
k2 Uninfected cell death day−1 0.05 — 0.075 0.05 
k3 Infection of uninfected cells L/(cells*day) 5 × 10−10 5 × 10−10 — 5 × 10−10 
k4 Infected cell death day−1 0.4 0.025 0.01 0.5 
k5 Production of new virus from infected cells day−1 40 — — 
k6 Clearance of virus day−1 2.95 — 0.008 
k7 Number of viruses produced by an infected cell virus/cell — 16.75 — — 
k8 Infection of uninfected cells, independent of virus L/(cells*day) — — 5 × 10−10 — 
k9 Clearance of infected cells by CTLs L/(cells*day) — — 0.01 0.001 
k10 Stimulation of CTLs by infected cells day−1 — — 10 
k11 CTL death day−1 — — 0.001 0.0001 

Many parameters were used in more than one model, but they did not retain the same value. For example, more than one model may have had a term for ‘Uninfected cell death’ or k2. However, the exact numerical value may have differed from model to model depending on what else was accounted for.

Table 3

An irregular structure of responses Yu1, Yu2 and Yu3

Event (uYu1 Yu2 Yu3 
 
  
 
  
 
Event (uYu1 Yu2 Yu3 
 
  
 
  
 

Yiu is the weighted observation y for response i and event u. An ‘X’ represents collected data. Many data are missing, thus no pattern is present.

Table 4

Data are present for all events, revealing a rectangular pattern

Event (uYu1 Yu2 Yu3 
Event (uYu1 Yu2 Yu3 

Fig. 2

Model predictions compared to experimental data for uninfected cells. Open circles are average titers for a specific day from baseline. All models achieved steady state in the first 150 days. The basic model oscillated initially before reaching steady state while the immune and composite models had a small initial decrease. The inset diagram provides details of the transient dynamics of uninfected cells.

Fig. 2

Model predictions compared to experimental data for uninfected cells. Open circles are average titers for a specific day from baseline. All models achieved steady state in the first 150 days. The basic model oscillated initially before reaching steady state while the immune and composite models had a small initial decrease. The inset diagram provides details of the transient dynamics of uninfected cells.

Fig. 3

Model predictions compared to experimental data for infected cells. Open circles are average titers for a specific day from baseline. The immune and composite models predicted immediate clearance of infected cells, while the simplified model predicted a gradual decline. A secondary axis on the right has an appropriate scale for the basic model, which showed a strong initial peak followed by a sharp decrease in cell concentration before a steady state was reached. Note that composite and immune models followed each other with a dramatic decline and clearance of cells.

Fig. 3

Model predictions compared to experimental data for infected cells. Open circles are average titers for a specific day from baseline. The immune and composite models predicted immediate clearance of infected cells, while the simplified model predicted a gradual decline. A secondary axis on the right has an appropriate scale for the basic model, which showed a strong initial peak followed by a sharp decrease in cell concentration before a steady state was reached. Note that composite and immune models followed each other with a dramatic decline and clearance of cells.

Fig. 4

Model predictions compared to experimental data for virus. Open circles are average titers for a specific day from baseline. The composite model predicted clearance of virions after approximately 500 days, while the simplified model predicted a gradual decline and clearance at a much later time. A secondary axis on the right has an appropriate scale for the basic model. This model predicted a precipitous increase followed by a sharp decrease in virus concentration before a steady state was reached.

Fig. 4

Model predictions compared to experimental data for virus. Open circles are average titers for a specific day from baseline. The composite model predicted clearance of virions after approximately 500 days, while the simplified model predicted a gradual decline and clearance at a much later time. A secondary axis on the right has an appropriate scale for the basic model. This model predicted a precipitous increase followed by a sharp decrease in virus concentration before a steady state was reached.

Fig. 5

Model predictions compared to experimental data for Tc lymphocytes. Open circles are average titers for a specific day from baseline. For the analysis, it was assumed that one percent of Tc cells were active CTL cells. The composite model slowly decreased over time, and the immune model had a much larger decline.

Fig. 5

Model predictions compared to experimental data for Tc lymphocytes. Open circles are average titers for a specific day from baseline. For the analysis, it was assumed that one percent of Tc cells were active CTL cells. The composite model slowly decreased over time, and the immune model had a much larger decline.

Table 5

Steady-state values of models

Model X Y V Z 
Basic 180 2.50 1.11 × 104 — 
Simplified — 0* 0* — 
Immune 267 2.47 × 10−9* — 1.23 × 10−5* 
Composite 200 6.20 × 10−7* 0.387* 6.20 × 10−2* 
Model X Y V Z 
Basic 180 2.50 1.11 × 104 — 
Simplified — 0* 0* — 
Immune 267 2.47 × 10−9* — 1.23 × 10−5* 
Composite 200 6.20 × 10−7* 0.387* 6.20 × 10−2* 

X, uninfected cells; Y, infected cells; V, free virus; Z, cytotoxic T lymphocytes. X, Y, Z are cells/mm3, V is particles/ml. Asterisk indicates that the steady-state value is reached beyond the life expectancy of a patient.

Table 6

Posterior probability share (π) for the four models under investigation

Model Number of equations/Number of parameters Variables 
\(\pi \left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)\)
with immune model 
\(\pi \left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)\)
without immune model 
Basic 3/6 X, Y, V 4.14 × 10−6696 6.14 × 10−2360 
Simplified 2/4 Y, V 6.75 × 10−4337 1.00 
Immune 3/7 X, Y, Z 1.00 — 
Composite 4/9 X, Y, V, Z 7.25 × 10−6067 1.07 × 10−1730 
Model Number of equations/Number of parameters Variables 
\(\pi \left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)\)
with immune model 
\(\pi \left({M}_{j}|\underset{\bar}{\underset{\bar}{Y}}\right)\)
without immune model 
Basic 3/6 X, Y, V 4.14 × 10−6696 6.14 × 10−2360 
Simplified 2/4 Y, V 6.75 × 10−4337 1.00 
Immune 3/7 X, Y, Z 1.00 — 
Composite 4/9 X, Y, V, Z 7.25 × 10−6067 1.07 × 10−1730 

‘Model’ is the name of the model, as defined in System and Methods section. The number of equations and number of parameters in each model are also listed. Variables listed are those monitored by the model. X, uninfected cells; Y, infected cells; V, virus particles; Z, CTL effector cells. The final two columns are the posterior probability share with and without the immune model. When all four models were compared with each other, the immune system model was overwhelmingly favored. When the immune model was not considered, the simplified model had the highest probability share.

We thank Dr Pamela Hartigan and Tenya Economou of the West Haven VA for providing the experimental data used in this research, as well as Dr Warren Stewart for helpful discussions regarding model discrimination. We also thank Alex Shutak for comments regarding the manuscript.

REFERENCES

Bagasra, O., Hauptman, S.P., Lischner, H.W., Sachs, M., Pomerantz, R.J.
1992
Detection of human immunodeficiency virus type 1 provirus in mononuclear cells by in situ polymerase chain reaction.
N. Engl. J. Med.
 
326
1385
– 1391
Bagasra, O., Seshamma, T., Oakes, J.W., Pomerantz, R.J.
1993
High percentages of CD4-positive lymphocytes harbor the HIV-1 provirus in the blood of certain infected individuals.
AIDS
 
7
1419
– 1425
Bieniasz, P.D., Ariyoshi, K., Bourelly, M.A., Bloor, S., Foxall, R.B., Harwood, E.C., Weber, J.N.
1993
Variable relationship between proviral DNA load and infectious virus titre in the peripheral blood mononuclear cells of HIV-1-infected individuals.
AIDS
 
7
803
– 806
Bonhoeffer, S., Coffin, J.M., Nowak, M.A.
1997
Human immunodeficiency virus drug therapy and virus load.
J. Virol.
 
71
3275
– 3278
Bonhoeffer, S., May, R.M., Shaw, G.M., Nowak, M.A.
1997
Virus dynamics and drug therapy.
Proc. Natl Acad. Sci. USA
 
94
6971
– 6976
Bonhoeffer, S., Funk, G.A., Gunthard, H.F., Fischer, M., Muller, V.
2003
Glancing behind virus load variation in HIV-1 infection.
Trends Microbiol.
 
11
499
– 504
Brenchley, J.M., Hill, B.J., Ambrozak, D.R., Price, D.A., Guenaga, F.J., Casazza, J.P., Kuruppu, J., Yazdani, J., Migueles, S.A., Connors, M., et al.
2004
T-cell subsets that harbor human immunodeficiency virus (HIV) in vivo: implications for HIV pathogenesis.
J. Virol.
 
78
1160
– 1168
Brinchmann, J.E., Albert, J., Vartdal, F.
1991
Few infected CD4+ T cells but a high proportion of replication-competent provirus copies in asymptomatic human immunodeficiency virus type 1 infection.
J. Virol.
 
65
2019
– 2023
Chun, T.W., Carruth, L., Finzi, D., Shen, X., DiGiuseppe, J.A., Taylor, H., Hermankova, M., Chadwick, K., Margolick, J., Quinn, T.C., et al.
1997
Quantification of latent tissue reservoirs and total body viral load in HIV-1 infection.
Nature
 
387
183
– 188
Cohen Stuart, J.W., Wensing, A.M., Kovacs, C., Righart, M., de Jong, D., Kaye, S., Schuurman, R., Visser, C.J., Boucher, C.A.
2001
Transient relapses (‘blips’) of plasma HIV RNA levels during HAART are associated with drug resistance.
J. Acquir. Immune Defic. Syndr.
 
28
105
– 113
De Boer, R.J. and Perelson, A.S.
1998
Target cell limited and immune control models of HIV infection: a comparison.
J. Theor. Biol.
 
190
201
– 214
Essunger, P. and Perelson, A.S.
1994
Modeling HIV infection of CD4+ T-cell subpopulations.
J. Theor. Biol.
 
170
367
– 391
Fauci, A.S.
2003
HIV and AIDS: 20 years of science.
Nat. Med.
 
9
839
– 843
Flint, S.J., Enquist, L.W., Racaniello, V.R., Skalka, A.M.
Principles of Virology: Molecular Biology, Pathogenesis, and Control of Animal Viruses
 
2004
2nd edn , Washington, DC ASM Press
Funk, G.A., Fischer, M., Joos, B., Opravil, M., Gunthard, H.F., Ledergerber, B., Bonhoeffer, S.
2001
Quantification of in vivo replicative capacity of HIV-1 in different compartments of infected cells.
J. Acquir. Immune Defic. Syndr.
 
26
397
– 404
Goldsby, R.A., Kindt, T.J., Osborne, B.A., Kuby, J.
Immunology
 
2003
5th edn , New York W. H. Freeman and Company
Hamilton, J.D., Hartigan, P.M., Simberkoff, M.S., Day, P.L., Diamond, G.R., Dickinson, G.M., Drusano, G.L., Egorin, M.J., George, W.L., Gordin, F.M., et al.
1992
A controlled trial of early versus late treatment with zidovudine in symptomatic human immunodeficiency virus infection. Results of the Veterans Affairs Cooperative Study.
N. Engl. J. Med.
 
326
437
– 443
Harper, M.E., Marselle, L.M., Gallo, R.C., Wong-Staal, F.
1986
Detection of lymphocytes expressing human T-lymphotropic virus type III in lymph nodes and peripheral blood from infected individuals by in situ hybridization.
Proc. Natl Acad. Sci. USA
 
83
772
– 776
Herz, A.V., Bonhoeffer, S., Anderson, R.M., May, R.M., Nowak, M.A.
1996
Viral dynamics in vivo: limitations on estimates of intracellular delay and virus decay.
Proc. Natl Acad. Sci. USA
 
93
7247
– 7251
Ho, D.D., Moudgil, T., Alam, M.
1989
Quantitation of human immunodeficiency virus type 1 in the blood of infected persons.
N. Engl. J. Med.
 
321
1621
– 1625
Hsia, K. and Spector, S.A.
1991
Human immunodeficiency virus DNA is present in a high percentage of CD4+ lymphocytes of seropositive individuals.
J. Infect. Dis.
 
164
470
– 475
Kitchen, S.G., Jones, N.R., LaForge, S., Whitmire, J.K., Vu, B.A., Galic, Z., Brooks, D.G., Brown, S.J., Kitchen, C.M., Zack, J.A.
2004
CD4 on CD8+ T cells directly enhances effector function and is a target for HIV infection.
Proc. Natl Acad. Sci. USA
 
101
8727
– 8732
Macatonia, S.E., Lau, R., Patterson, S., Pinching, A.J., Knight, S.C.
1990
Dendritic cell infection, depletion and dysfunction in HIV-infected individuals.
Immunology
 
71
38
– 45
Nelson, P.W., Mittler, J.E., Perelson, A.S.
2001
Effect of drug efficacy and the eclipse phase of the viral life cycle on estimates of HIV viral dynamic parameters.
J. Acquir. Immune Defic. Syndr.
 
26
405
– 412
Nowak, M.A., Bonhoeffer, S., Shaw, G.M., May, R.M.
1997
Anti-viral drug treatment: dynamics of resistance in free virus and infected cell populations.
J. Theor. Biol.
 
184
203
– 217
O'Brien, W.A., Hartigan, P.M., Martin, D., Esinhart, J., Hill, A., Benoit, S., Rubin, M., Simberkoff, M.S., Hamilton, J.D.
1996
Changes in plasma HIV-1 RNA and CD4+ lymphocyte counts and the risk of progression to AIDS. Veterans Affairs Cooperative Study Group on AIDS.
N. Engl. J. Med.
 
334
426
– 431
Perelson, A.S.
2002
Modelling viral and immune system dynamics.
Nat. Rev. Immunol.
 
2
28
– 36
Perelson, A.S., Neumann, A.U., Markowitz, M., Leonard, J.M., Ho, D.D.
1996
HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time.
Science
 
271
1582
– 1586
Phillips, A.N.
1996
Reduction of HIV concentration during acute infection: independence from a specific immune response.
Science
 
271
497
– 499
Pomerantz, R.J. and Horn, D.L.
2003
Twenty years of therapy for HIV-1 infection.
Nat. Med.
 
9
867
– 873
Poznansky, M.C., Walker, B., Haseltine, W.A., Sodroski, J., Langhoff, E.
1991
A rapid method for quantitating the frequency of peripheral blood cells containing HIV-1 DNA.
J. Acquir. Immune Defic. Syndr.
 
4
368
– 373
Re, M.C., Furlini, G., Gibellini, D., Vignoli, M., Ramazzotti, E., Lolli, E., Ranieri, S., La Placa, M.
1994
Quantification of human immunodeficiency virus type 1-infected mononuclear cells in peripheral blood of seropositive subjects by newly developed flow cytometry analysis of the product of an in situ PCR assay.
J. Clin. Microbiol.
 
32
2152
– 2157
Reddy, B. and Yin, J.
1999
Quantitative intracellular kinetics of HIV type 1.
AIDS Res. Hum. Retroviruses
 
15
273
– 283
Ribeiro, R.M. and Bonhoeffer, S.
1999
A stochastic model for primary HIV infection: optimal timing of therapy.
AIDS
 
13
351
– 357
Schnittman, S.M., Greenhouse, J.J., Psallidopoulos, M.C., Baseler, M., Salzman, N.P., Fauci, A.S., Lane, H.C.
1990
Increasing viral burden in CD4+ T cells from patients with human immunodeficiency virus (HIV) infection reflects rapidly progressive immunosuppression and clinical disease.
Ann. Intern. Med.
 
113
438
– 443
Schnittman, S.M., Psallidopoulos, M.C., Lane, H.C., Thompson, L., Baseler, M., Massari, F., Fox, C.H., Salzman, N.P., Fauci, A.S.
1989
The reservoir for HIV-1 in human peripheral blood is a T cell that maintains expression of CD4.
Science
 
245
305
– 308
Simmonds, P., Balfe, P., Peutherer, J.F., Ludlam, C.A., Bishop, J.O., Brown, A.J.
1990
Human immunodeficiency virus-infected individuals contain provirus in small numbers of peripheral mononuclear cells and at low copy numbers.
J. Virol.
 
64
864
– 872
Snedecor, S.J.
2003
Comparison of three kinetic models of HIV-1 infection: implications for optimization of treatment.
J. Theor. Biol.
 
221
519
– 541
Stafford, M.A., Corey, L., Cao, Y., Daar, E.S., Ho, D.D., Perelson, A.S.
2000
Modeling plasma virus concentration during primary HIV infection.
J. Theor. Biol.
 
203
285
– 301
Stebbing, J., Gazzard, B., Douek, D.C.
2004
Where does HIV live?
N. Engl. J. Med.
 
350
1872
– 1880
Stevenson, M.
2003
HIV-1 pathogenesis.
Nat. Med.
 
9
853
– 860
Stewart, W.E. and Sørenson, J.P.
1981
Bayesian estimation of common parameters from multiresponse data with missing observations.
Technometrics
 
23
131
– 141
Stewart, W.E., Caracatsios, M., Sørenson, J.P.
1992
Parameter estimation from multiresponse data.
AIChE J.
 
38
641
– 650
Stewart, W.E., Henson, T.L., Box, G.E.P.
1996
Model discrimination and criticism with single-response data.
AIChE J.
 
42
3055
– 3062
Stewart, W.E., Shon, Y., Box, G.E.P.
1998
Discrimination and goodness of fit of multiresponse mechanistic models.
AIChE J.
 
44
1404
– 1412
Tan, W.Y. and Wu, H.
1998
Stochastic modeling of the dynamics of CD4+ T-cell infection by HIV and some Monte Carlo studies.
Math. Biosci.
 
147
173
– 205
Tuckwell, H.C. and Le Corfec, E.
1998
A stochastic model for early HIV-1 population dynamics.
J. Theor. Biol.
 
195
451
– 463
Wein, L.M., D'Amato, R.M., Perelson, A.S.
1998
Mathematical analysis of antiretroviral therapy aimed at HIV-1 eradication or maintenance of low viral loads.
J. Theor. Biol.
 
192
81
– 98
Wodarz, D. and Nowak, M.A.
2002
Mathematical models of HIV pathogenesis and treatment.
Bioessays
 
24
1178
– 1187

Comments

0 Comments