Modeling Adaptive Regulatory T-Cell Dynamics during Early HIV Infection

Regulatory T-cells (Tregs) are a subset of CD4+ T-cells that have been found to suppress the immune response. During HIV viral infection, Treg activity has been observed to have both beneficial and deleterious effects on patient recovery; however, the extent to which this is regulated is poorly understood. We hypothesize that this dichotomy in behavior is attributed to Treg dynamics changing over the course of infection through the proliferation of an ‘adaptive’ Treg population which targets HIV-specific immune responses. To investigate the role Tregs play in HIV infection, a delay differatial equation model was constructed to examine (1) the possible existence of two distinct Treg populations, normal (nTregs) and adaptive (aTregs), and (2) their respective effects in limiting viral load. Sensitivity analysis was performed to test parameter regimes that show the proportionality of viral load with adaptive regulatory populations and also gave insight into the importance of downregulation of CD4+ cells by normal Tregs on viral loads. Through the inclusion of Treg populations in the model, a diverse array of viral dynamics was found. Specifically, oscillatory and steady state behaviors were both witnessed and it was seen that the model provided a more accurate depiction of the effector cell population as compared with previous models. Through further studies of adaptive and normal Tregs, improved treatments for HIV can be constructed for patients and the viral mechanisms of infection can be further elucidated.


Introduction
Although the dynamics of CD4 + and CD8 + cells have been well characterized, in HIV infection, there is currently a lack of understanding concerning the role of regulatory T-cells, or Tregs in viral dynamics [1]. Tregs are a class of CD4 + cells which limit the activation and expansion of immune cells, including autoreactive CD4 + cells and CD8 + cells. Initial studies provided evidence that the Treg response to HIV was beneficial, limiting immune exhaustion and immune-mediated tissue damage [2][3]. Conversely, Tregs have been observed to contribute to the onset of immune dysfunction and to prevent a successful immune response [3][4][5]. Finally, there is evidence that the role of Tregs throughout infection may follow a more dynamic behavior, changing its behavior at different stages of infection [3] [6]. Divergent reports on Treg activity can be attributed in part to experimental obstacles involved in studying their dynamics in vivo. Patient data for Tregs are largely non-existent due to the absence of accurate surface markers to characterize the population. Without selective markers, computational modeling is paramount in providing insight into Tcell regulation during HIV infection, specifically examining adaptive Treg behavior.
To further exacerbate the experimental problem in HIV infection, there is evidence that multiple subsets of Tregs exist: normal Tregs, as well as a HIV-specific adaptive Tregs [3] [6][7][8]. Normal Tregs (nTreg), the body's naturally occuring regulatory T-cells, are present in the early stages of infection, but through some unknown mechanism, HIV may adapt the activity of nTregs, to the benefit of the virus; this modified class of Tregs is deemed adaptive Tregs (aTreg) [6].
Adaptive Tregs are believed to be deleterious for the patient by contributing to viral proliferation and poor immune activity. Although there are no known markers for adaptive Tregs, their existence may have a distinct effect on HIV dynamics [6]. In our model, we included both nTregs and aTregs, matching simulations to patient data and constructed an understanding of whether the two subsets could biologically exist together with differing effects. Our model focuses on the effect Tregs have on sharply declining viral load during acute infection and the examination of the importance of CD8 + activation and CD4 + population limitation as well as its subsequent steady state levels during latency.

Model Provides Realistic Effector Cell Dynamics which are Lacking in Previous Models
Through data fitting, the model was initially tested to ensure normal viral behavior could be replicated. Consistent with previous fits, two main behaviors were observed -oscillatory as well as steady state dynamics [9] (Fig. 1A, C). The model was capable of reproducing the acute viral infection and proceeding either to steady state or damped oscillations. Comparing our fits to a previously published study (Cuipe et. al), lacking normal and adaptive regulatory T-cells, reveals that similar dynamics could be obtained through the inclusion of a regulatory T-cell population. Notably however, Treg fits showed improvements in simulating effector cell populations, within the first few weeks of infection, compared to previous models (Fig. 1B, D). In comparing fits between the Ciupe model and our model, we find comparable fits; however, when we examine the ranges of effector cells, we find a significant improvement in our model. Our model is noted as staying within physiologically relevant effector population ranges for the majority of the simulation, whereas the Ciupe model never enters physiological ranges. In fact, the Ciupe model displays a sharp decline in effector populations early in simulation not seen biologically. Comparing the effector populations with data on effector populations reveals that our model provides a more biologically relevant dynamic with comparable fits to viral load [10]. We believe these improvements are due to the usage of dual equations modeling the effector cell population. Through compartmentalization of effector cells into immature and mature populations, we were able to avoid sharp drops in these cells during acute infection and construct a method of modeling effector populations within realistic physiological bounds.

Sensitivity of Normal Tregs on Viral Load
Partial rank correlation coefficient (PRCC) analysis was performed in order to examine the importance of the individual parameters on viral load. PRCC isolates a single parameter for analysis independent of others (further explained in the methods) and examines its behavior in the system. This analysis is important in systems where parameters are highly sensitive or not well known in the literature, as individual contributions must be taken into account when making predictions. Our model considers welldefined parameters for viral infection as well as newly introduced parameters such as the rate at which normal Tregs downregulate CD4 + cells (a 1 ): a 1 accounts for the rate of change of CD4+ cells in response to normal regulator T-cells. Previous models have not considered this CD4+ cell response directly but recent research has suggested regulatory T-cells do effect CD4+ cells throughout the infection process [3] [5]. We noticed a clear bifurcation in viral load behavior when isolating (a 1 ) ( Fig. 2A). Two clusters were observed (cluster 1 and 2) with divergent viral behaviors. Parameter sets within cluster 1 were simulated and complete viral load clearance was seen within the first two weeks of infection for all points found within the cluster, thus signifying a largely nonbiological dynamic (Fig. 2B). Parameter sets within cluster 2 revealed high steady state viral loads representative of the patient fits (Fig. 2C).

Sensitivity of Treg Downregulation of CD4 + on Viral Load
Bifurcated clustering was further examined through numerically sampling, physiologically bound parameter sets defined in Table 1 (a 1~0 :0{1:0(day) 21 ). Since a 1 has values that are currently unknown in the literature, we first varied it over three orders of magnitude. All of our results, however, suggest that a 1 is bounded between (0:0,1:0). Viral loads at steady state (t~200 days) were recorded for each parameter set over the a 1 range. Simulations reveal a 1 parameters to be (w0:6) and (0:0{0:8) for clusters 1 and 2, respectively (Fig. 3A). Steady state behavior is relatively constant among the various curves; specifically, viral load seems unchanged over small values of a 1 , however, a clear bifurcation to zero occurs when a 1 is sufficiently large (w0:6 day 21 ). These results suggest a physiological range in normal Treg behavior.

Sensitivity of Transition Rate of normal Tregs on Viral Load
Ranging c (the transition rate of normal Tregs into adaptive Tregs) revealed surprising dynamics in viral load, particularly during early infection (Fig. 4). It was seen that viral load was highly dependent on c in the acute stage of infection (t~80 days), but was to a lesser degree dependent during steady state (t~200 days). For some parameter sets, it was even found that c~0 would result in viral clearance; thus for these simulations it was necessary to have an adaptive population for simulation to achieve steady state behavior. In addition, it was revealed that a larger value forc, i.e., a larger rate of nTregs becoming aTregs, results in, on average, a 20% increased viral load at steady state.

Impact of Treg Parameters on Viral Load and Cellular Populations
To elucidate the dynamics of the system, different populations were examined over increasing values of the Treg-related parameters c and a i : The best fit for Patient 8 was utilized as the starting parameter set, and the respective parameters were varied over several orders of magnitude (Fig. 5). T-cell concentration, viral load, as well as effector cell concentration (measured as the sum of immature and mature effectors) were graphed against a i and c (Fig. 5). Through the plots, it was seen that steady state concentrations of the three populations were relatively unchanged. The exception was in the case of adaptive transition rate (c) where low values of c resulted in a moderate decrease in steady state viral load. However, it was generally seen that the main differences in dynamics occurred in the acute stage of infection, i.e. for tv100 days.
An important finding was that a 1 was found to be inversely proportional to peak viral load, suggesting that reducing the target population can limit the proliferation of virus during the early stages of infection. a 2 did not seem to have much of an effect on Tcell and virus dynamics; however for large values of a 2 , there was an early effector population drop which then rebounded to a higher level. When a 3 was increased, it was found that viral load would increase, as expected due to the fact that a 3 controls the rate of downregulation of mature effectors by adaptive Tregs. However, it was also found that large a 3 was inversely proportional to T-cells, which gave a surprising result (Fig. 5D). It can be seen that viral load can remain elevated during T-cell depletion (Fig. 5D). It was additionally found that c~0 resulted in low viral load while cw0 on the same order as the data fitting would result in an elevated viral load, suggesting adaptive Tregs can have a deleterious effect on the patient (Fig. 5A).

Discussion
Although the behavior of regulatory T-cells in HIV infection is largely unknown, recent literature contrasts heavily on whether such T-cell populations are helpful or deleterious to the patient [3]. It has been theorized that there are two subsets of Tregs, a basal population of normal Tregs, and a population of adaptive Tregs derived from interactions between normal populations and viral particles [3] [6]. While the explicit mechanism for how such a population can emerge is still unknown, it has been hypothesized to occur at the transcriptional level resulting in differential regulatory T-cell distribution in various regions of the body [6]. Here, a physiologically relevant model of normal and adaptive Tregs was utilized to model patient viral loads and examine the effect of normal Treg parameter boundaries and an adaptive Treg subgroup.
Normal regulatory T-cells, specifically the efficiency of their clearance of CD4 + cells, were characterized throughout the model by the parameter a 1 : Using the analysis as discussed, a 1 was observed to range from 0.0-0.7 day 21 due to a bifurcation in viral steady state at higher efficiencies. This suggests a defined range on normal Treg action and acts as a first step in studying the dynamics of this regulatory population within the framework of HIV infection. These bounds additionally suggest that nTreg efficiency in down-regulating CD4 + cells can drastically affect viral load behavior. Specifically, if efficiency is relatively high, the number of target cells for HIV is reduced enough to prevent proliferation of virus due to a clearance of the target cell population. The limits to which nTreg dynamics can be manipulated in vivo will elucidate the dynamics of regulatory Tcells within HIV dynamics and give light to potential therapeutic techniques against the infection.
Adaptive Tregs are suggested in our model to counteract the immune system's defenses against HIV through depletion of the effector cell populations. We reproduce normal viral load behavior in the presence of aTregs (Fig. 1) and through the sensitivity dynamics of the aTreg parameter c observe that, in some instances, the presence of an adaptive population is essential to produce realistic viral behavior. Additionally, the adaptive transition parameter, c, appears to influence viral loads; specifically, larger c values result in higher steady state viral loads and increased T-cell counts (Fig. 5). Distinct biological markers will be necessary for further study of adaptive T-cell populations; however, the possibility of two classes of regulatory T-cells is plausible within the context of mature effector cell depletion. Additional mechanisms of aTreg action, such as depletion of other immune cell populations or other events would be of great interest to further examine Treg-HIV dynamics. . Models correspond to best fits to the viral load data from two separate representative patients. (inset) Patient data is shown (circular data points), and viral load is shown (red). Parameter values for simulations can be found in Table 1. (B,D) Comparative plots of effector populations using Treg model as compared to a previously published model [9] as well as physiological ranges for effector cells [10] as seen by the red horizontal lines in B and D. In all cases, we notice that our model does not dip to unrealistic ranges initially as does the Ciupe model. However, the viral load fits between both models are nearly the same suggesting a degree of robustness between both models. doi:10.1371/journal.pone.0033924.g001  The use of two effector cell populations, immature and mature, proved successful in maintaining realistic cell population concentrations and while producing comparable viral fits to previous models. Sensitivity analyses reveal that varying parameters related to regulatory T-cells largely influence the acute stages of infection but are relatively robust at long times (tw100 days). Also, we observe regulatory T-cells to be pivotal in early viral peaks and early depletion of CD4 + T-cells.
While production and death rates for the regulatory T-cells could be extrapolated from similar parameter values approximated in previous research, little was known about a i and c, which reflect the downregulatory capabilities of the nTregs and aTregs as well as the rate at which normal Tregs migrate into the adaptive population. Parameter fitting gives bounds on these values and desirable visualization to differentials in viral load dynamics. Through further studies of Tregs and their different classes, new therapeutic techniques can be developed to potentially limit normal Treg action and reflect optimal viral load reduction providing additional insights into viral and immune system's capabilities during infection.

Ethics Statement
There was no need for consent in the work presented here as no work with patients was performed.
To examine the dynamics of HIV viral load during primary infection, we propose a system of delay differential equations consisting of seven variables: normal Tregs (R), adaptive Tregs (R a ), target CD4 + cells (T), productively infected CD4 + T cells (T * ), free virus (V), immature CD8 + cells (E i ) and mature CD8 + cells (E m ).

Normal and Adaptive Regulatory T-Cells
Here, we introduce normal Tregs as well as the hypothesized adaptive Tregs. It has been suggested that HIV interacts with normal Tregs and induces them into a secondary adaptive population which can have deleterious impacts on the patient and assist the virus. To model this dynamic, normal Tregs (R) were given a production rate of p R and a death rate of d R : The death rate of adaptive Tregs (R a ) was defined by d R : The massaction term cRV was utilized, in which c is the proportion of interactions between R and V which result in adaptive Tregs.  The following are the equations for the R and R a populations, dR a dt~c RV {d R R a : ð2Þ For viral load, we use the same form as seen in the Perelson, Nelson 1999 model, For the CD4 + population, an almost identical form as the Perelson, Nelson 1999 model is utilized, with a proliferation term s, a half-life of 1/d, and a bikinetic expression kVT to express the notion that interaction of virus with T-cells at some probability k results in the infection of T-cells. In the literature, it is well characterized that Tregs the proliferation rate of T-cells [1] [3][4]. Thus, the term a 1 R was introduced, which serves as the basal downregulation of CD4 + cells by the normal population of Tregs. The uninfected T-cell population equation is as follows, Infected T-cells were produced from the bikinetic interaction (kVT) in which uninfected T-cells interact with HIV to become infected. The half-life for these infected cells is 1=d with a clearance term d x E m T Ã ; thus infected cells are cleared also through a bikinetic interaction with mature CD8 z effectors at a rate d x ,

Effector Cell Equations
Various models of HIV include or exclude cytotoxic CD8 + Tcells; depending on the specific dynamics of interest [9] [11]. In the laboratory, it has been observed that Tregs impact the activity of these cells, thus we include them as a key element in our model. Some previous models [9], use a single equation to represent effector cells, however, this population would often trend to nonphysiological conditions. To adjust for this issue, we separated the effector population into two subsets: an immature class of effectors and a mature class. Since effector cells need to be presented with HIV-antigen in order to specifically combat infected cells, we write our equations with immature effectors interacting with HIV in order to become mature effectors.
For immature effectors, we assume similar dynamics to the uninfected T-cell equation. There is a constant rate of production p Ei , a death rate d Ei and a loss of cells due to conversion into mature effector cells, dependent on virus interaction at some rate k E : Additionally, there is a downregulation of the immature population by normal Tregs, expressed through the a 2 R term, Mature effectors are produced through viral interaction with the immature effector population (kVE i ) and die at a rate d Em : A time delay was used in the production term of the mature effectors to represent the biological lag that occurs between immature effector recognition of virus and the production of mature effectors. Adaptive Tregs are hypothesized to specifically downregulate the mature effector response through the expression a 3 R a , thus suppressing normal immune responses, Parameter Fitting Viral load patient data were used from data collected in [12]. Fitting was accomplished through running 2,000 iterations using a Simplex optimization algorithm through JSim Software [13]. All parameters were varied with known parameters from [9] used as start values in the optimization; parameter descriptions, ranges and best fits can be found in Table 1. Parameters which were relatively unknown, for instance production rate of nTregs, were extrapolated from similar parameters and were varied over an order of magnitude to ensure that their starting points were biologically feasible. A bootstrapping algorithm was used in order to find confidence intervals for the parameters related to Treg dynamics (Table 1).

PRCC Analysis and other Sensitivity Analysis
In the analysis of our model, we were interested in determining which parameters have a strong effect on key outputs such as viral load. This is especially important for parameters which are not well-characterized in the literature, which included terms related to normal and adaptive regulatory T-cell function; if these have a strong influence on model outputs, we must be careful to qualify any conclusions drawn from model analysis with the disclaimer that such conclusions may be very sensitive to rough parameter estimates. The results of a sensitivity analysis can be used as a guide in determining which parameters are most important, whether to re-express the model in terms of better-understood parameters, and potentially which parameters can be measured or investigated more thoroughly.
The Partial-Rank Correlation Coefficient (PRCC) measures the association between two variables, with the effects of other variables removed. Given a model parameter of interest and a model output at a particular time point, the PRCC can be used as a sensitivity metric of the output to the parameter, independent of the other parameters, as described in Marino et al. [14].
The first step of computing the PRCC of parameter x j and output y is to rank-transform the data in variablesx x j andŷ y: The PRCC is then given by the correlation between the residuals (y{ŷ y) and (x j {x x j ), where, wherex x j andŷ y are given by [14], y y~b 0 z X k p~1,p=j b p x p : ð9Þ A high correlation between these residuals suggests that output y is sensitive to parameter x i : We can apply this technique to our HIV model by setting the model output to, for example, viral load at time 80 days after primary infection. PRCC analysis can give us a snapshot of the influential parameters on viral load at this time point.
In addition to PRCC analysis, a second sensitivity analysis was performed for the parameters a i and c: One of the parameters was isolated, and the other parameters were chosen through a Latin Hypercube Sampling (LHS). The isolated parameter was ranged over several orders of magnitude and for each value of the parameter, viral load was simulated and recorded at t~80 days as well as t~200 days. Plots of viral load vs. parameter value were constructed to depict the dynamics. Analysis was performed for 100 distinct parameter sets, with representative plots in (Fig. 3, 4). In both cases, we find using the PRCC analysis that a i and c are robust parameters and viral dynamics are not highly sensitive to small perturbations in either parameter, suggesting that our new model is more robust then previous models.