It has been suggested that HIV-1 has evolved its set-point virus load to be optimized for transmission. Previous epidemiological models and studies into the heritability of set-point virus load confirm that this mode of adaptation within the human population is feasible. However, during the many cycles of replication between infection of a host and transmission to the next host, HIV-1 is under selection for escape from immune responses, and not transmission. Here we investigate with computational and mathematical models how these two levels of selection, within-host and between-host, are intertwined. We find that when the rate of immune escape is comparable to what has been observed in patients, immune selection within hosts is dominant over selection for transmission. Surprisingly, we do find high values for set-point virus load heritability, and argue that high heritability estimates can be caused by the ‘footprints’ left by differing hosts' immune systems on the virus.
Human immunodeficiency virus type 1 (HIV-1) evolves under two levels of selection. On the one hand, there is within-host selection for immune escape. On the other hand, selection on the population-level acts on infectiousness and virulence. In this paper, we explore how these two levels of selection are intertwined, keeping in mind the massive heterogeneity of the hosts with respect to their cellular immune responses. A HIV-1 infection can be separated into three phases: the acute phase, the asymptomatic phase and the symptomatic (or AIDS) phase. During the acute phase, the virus establishes high virus loads (the number of HIV-1 RNA copies per ml blood plasma) [1], until the CD4+ target cells are depleted [2], and adaptive immune responses start limiting viral reproduction. The virus load then drops to a semi-stable level called the set-point. This marks the beginning of the asymptomatic or chronic phase, during which the partially restored CD4+ T-cell count gradually drops, and at some point patients develop AIDS. The set-point virus load (spVL) differs markedly between individuals. In untreated patients, spVL ranges from 102 to 106 copies/ml. The origin of this variation is an extensively researched topic, and explanations include host and viral factors. For instance, host factors incorporate the association between the set-point and the Human Leukocyte Antigen (HLA) haplotype, which is important for cellular immunity [3]–[6]. The observation that the spVL is to some extent heritable [7]–[14], suggests that viral genetic factors sway the set-point too. The exact extent of this heritability is unknown, as estimates range from 6% to 59%. spVL is related to infectiousness and virulence. Patients with a higher spVL tend to be more infectious [15], but also develop AIDS more rapidly [16], resulting in a trade-off between infectiousness and the length of the asymptomatic phase. This life history trade-off was identified by Fraser et al. [17], and opens the door for HIV-1 adaptation with respect to transmission by means of spVL evolution. Certain spVLs (around copies/ml) allow a HIV-1 strain to cause more secondary infections than strains with lower or higher set-points. A strain that establishes on average this optimal set-point should therefore become more abundant in the population. The striking observation is that, although large variation in set-points exists, most HIV-1 infected patients show a set-point close to the transmission-optimal value [17]. Moreover, mathematical models show that this adaptation can take place within realistic time scales [18], given the heritability estimates of spVL [7], and HIV-1' s likely dates of origin [19], [20]. In such mathematical models, HIV-1' s population-level fitness (measured in terms of the basic reproduction number) is only constrained by the life history trade-off, and environment- and mutation-induced spVL-variation. It is therefore quite intuitive that in such a model evolution leads to intermediate levels of spVL [17], [18]. The inclusion of directed within-host evolution in such models introduces an extra constraint on the population-level fitness; one which dominates the evolutionary outcome, unless within-host selection is exceedingly weak. For a homogeneous host population, this has been shown recently by Lythgoe et al. [21], and they suggest that within-host evolution of traits affecting virus load must be slow. Below we argue that ‘short-sightedness’ [sensu 21,22], i. e. , the life history trade-off has no apparent effect on the evolutionary outcome, can easily be understood when the host population is homogeneous. However, in a much more realistic situation where HIV-1 needs to escape from immune responses that vary markedly between individuals, the same intuition for the effect of directed within-host evolution can no longer be applied, and needs to be revised. In this study, we explicitly incorporate such immune selection and massive host-heterogeneity with respect to immune responses in a nested epidemiological model. We investigate whether spVL evolution of HIV-1 is influenced by the virus' life history trade-off. Our model predicts that within-host immune selection has a major influence on population-wide spVL evolution. Thus, both Lythgoe' s and our model predict short-sighted spVL evolution. However, we do not agree that within-host evolution must therefore be slow. Throughout the paper, we use the term ‘between-host adaptation’ for evolutionary dynamics where HIV-1' s life history trade-off notably affects the evolution of spVL. The term ‘within-host selection’ refers to selection for immune escape and reversion of deleterious mutations. At the same time, we use our model to investigate spVL heritability. We argue that high heritability can be a result of HIV-1 rapidly escaping immune responses, and the between-individual variation of these responses. We emphasize that spVL heritability caused by such a mechanism does not provide support for between-host adaptation.
Our approach combines a caricature model for immune escape with a susceptible-infectious (SI) model for HIV-1 transmission. Both the within-host and the between-host simulations are discrete-event and individual based. The technical details are given in Methods. Here we give an intuitive exposition. Cytotoxic T-Lymphocyte (CTL) responses are arguably important for controlling HIV-1 virus load [23], [24]. Human cells notify the cellular immune system about their proteome by presenting peptides on HLA molecules. On infected cells, a subset of these peptides originate from viral proteins. If a CTL clone detects such a foreign peptide, it can kill the infected cell, and the peptide (in its proper HLA context) is called an epitope. Not all peptides can be presented by the HLA molecules of a host, and HIV-1 can escape from CTL recognition by mutating amino acids in its peptides to prevent presentation by the host' s HLA molecules [25]–[27]. Due to HLA-polymorphism, the particular subset of all peptides that can be presented by a host' s HLA molecules (the binding repertoire) differs strongly between individuals [28]. In our model we incorporate this by assuming that a wild-type virus has n peptides that can be presented in the population. A particular host can present a subset of size k of these n peptides. During infection, we assume that mutations in the n potentially recognized peptides occur according to a Markov process. Some of these mutations will result in CTL escape (escape mutations). In this case, the mutant takes over the viral population in that host. Naturally, if two hosts have have a common peptide in their binding repertoires, the mutated peptide is a CTL escape for both hosts. In line with evidence, we assume that escape mutations in HIV-1 come with a fitness cost [29], [30]. The total fitness effect of an escape mutation, resulting from immune escape and its fitness cost, must be positive before the escape mutant can replace the dominant HIV-1 strain in the host. In order to model this, we use the virus load in the asymptomatic phase as a measure for within-host fitness. An immune response causes a reduction in the log10 virus load, and a fitness cost of any mutation reduces the log10 virus load by. The total fitness effect of an escape mutation is then a increase in the log10 virus load. In the simulations, we choose and so that lies within estimated ranges [31], [32]. Qualitatively, our results do not depend on these particular choices for and, as long as (results not shown). Certain hosts have an efficient immune response to HIV-1. This can partially be explained by HLA-type. For instance, HLA-B*57, B*27, B*58 and B*18 are associated with a low spVL. HIV-1 is able to escape immune responses in hosts with these HLA-types, but the associated fitness costs tend to cripple the virus [25]. When such a crippled virus is transmitted to the next host, lacking the protective HLA-type, the virus load in this secondary host can remain low for a long time [29]. After a while, the crippled virus reverts the deleterious mutations, since the immune pressure causing these crippling mutations is not present in the secondary host [33]. We propose that this effect is not only restricted to known protective HLA-types, but holds more generally [e. g. , see 34]. We model this similar to immune escape. As a result of immune escape in previous hosts, a viral strain may carry a number of deleterious mutations. These mutations can revert to the wild-type, again according to a Markov process. In summary, our model for the log10 virus load V is [cf. 35]where Vmax is the log10 virus load of a HIV-1 strain without deleterious mutations in the absence of CTL-responses (), e. g. , the high virus load observed in a CD8+ T cell depleted individual [23], [36], [37]. The integer e represents the number of escape mutations in a host (and hence, equals the number of immune responses), and f denotes the number of deleterious mutations. In other words, f equals the number of mutated peptides outside the current host' s binding repertoire. We assume that escapes and reversions appear at a rate proportional to the number of immune responses and deleterious mutations, respectively. Hence where and are the ‘per-peptide’ rate of escape and reversion, respectively. We will refer to and as ‘mutation rates’. Keep in mind, however, that our model of escape and reversion is quite phenomenological. The rates and are a combination of many factors, such as the error rate during reverse transcriptase and the fixation rate. Moreover, the rates and should in reality depend on the virus load. We simplify this dependence by assuming that the rates differ only between disease phases. In the acute and AIDS phase the per-peptide rates are high and in the asymptomatic phase, these rates are lower. Instead of and, we therefore take distinct parameters and for the per-peptide mutation rate in the acute (), asymptomatic () and AIDS () phase. We choose with, meaning that reversion is slower than escape (see Table 1 for the exact parameterization). This is in line with the assumption that the total fitness benefit of an escape mutation is greater than the benefit of a reversion (i. e.). As mentioned earlier, a HIV-1 strain infecting a new host carries a history of mutations acquired in previous hosts [29], [38]. In the context of the new host, many of these mutations will not be beneficial. Some of them may be advantageous, because HLA molecules can share epitopes [39], [40], and individuals share HLA molecules. To keep our model simple, we assume that a random host' s binding repertoire is a random subset of size k of the set of all n possible HIV-1 epitopes. In reality, HLA haplotypes, and hence binding repertoires, are less regularly distributed. However, our simpler distribution provides us with the advantage that we only have to keep track of the number of mutated peptides. Namely, when a host transmits a virus with e escape mutations and f deleterious mutations (denoted as an -virus), then in the secondary host the virus will have phenotype with. We find the number of escape mutations e′ by choosing a new random binding repertoire of size k′. Since every peptide is part of the new binding repertoire with equal probability, the number of a priori escape mutations is drawn from the hypergeometric distribution (). An example of how a virus' phenotype can differ between hosts is given in Figure 1. By default, we choose and (see Table 1), such that about 10% of HIV-1' s peptides can serve as an epitope. The number k is chosen such that hosts have a realistic number of responses, also when many of the n peptides are mutated. We model the three phases of a HIV-1 infection based on Fraser et al. [17] and Hollingsworth et al. [41]. The acute phase has a fixed length D1, and in this phase individuals have a fixed infectiousness β1. After D1 years, the asymptomatic phase starts and infectiousness β2 (V) and the average length of the asymptomatic phase D2 (V) depend on the virus load V. The functions β2 and D2 are Hill functions with coefficients as estimated by Fraser et al. [17]. When the asymptomatic phase ends, the AIDS phase starts. This AIDS phase has, similar to the acute phase, a fixed length D3 and fixed infectiousness β3. We do not incorporate any correction for serial monogamy on infectiousness. As an illustration of the within-host model, we have simulated a large number of within-host processes for two different parameter settings (Figure 2). Stochasticity and host-heterogeneity cause large variation in the within-host evolution of the virus (the thin step-wise lines in Figure 2). As deleterious mutations are reverted and CTL responses are escaped, the virus load increases during the infection. If the mutation rate is high, almost all escapes happen during the acute phase of the infection. The cohort-average virus load (the heavy blue line in Figure 2, bottom panels) can then even decrease, since individuals with a high set-point develop AIDS more rapidly. When these fast progressors die, we exclude them from the calculation of the cohort' s mean virus load. Notice that during the acute phase, the variable V does not reflect the peak virus loads observed in patients, but is merely a measure of the virus' fitness. For the between-host model, we explicitly model a population of infected individuals (of size I), and assume a frequency-dependent contact process with susceptibles [42]. Super-infection and co-infection are ignored. We keep the total population size (N) constant, and only keep track of the susceptibles' number (S). Because of within-host evolution, an individual may transmit different viral strains during the course of an infection. When the virus load increases due to within-host adaptation, the infection rate also increases. We verified that a model with a non-constant population size does not give different results (not shown). Since in our model the virus load can increase during the asymptomatic phase, we need to specify what we mean with set-point virus load. We define the spVL (in log10 scale) as the geometric average of the log10 virus loads in the asymptomatic phase, i. e. , , where the integral is taken over the chronic phase, which lasts L years, and V (t) denotes the virus load at time t. We often write to indicate the population-wide arithmetic average spVL. Bracket notation is also used for other population-wide averages. When we choose the mutation rate low and run the agent-based model, the mean spVL converges to 4. 52log10 copies/ml; the value optimal for transmission (Figure 3A). However, this takes many centuries, depending on the maximal virus load Vmax and the initial number of mutations. By increasing the mutation rate, we make the evolutionary dynamics faster, but lose between-host adaptation (Figure 3B). In fact, the mean spVL is approximately 1. 3 log10 copies/ml higher than 4. 52. By keeping the mutation rate equally high, but lowering Vmax, the HIV-1 quasi-species can be given a population-level fitness () that is about 17% higher than what is reached in Figure 3B. Apparently, selection for spVL values that are optimal for transmission is overruled by within-host selection at high mutation rates. Both simulations in Figure 3 are approaching different steady states. Thus, to investigate between-host adaptation further, we now look at the properties of the model in population-level steady state for many different parameter combinations (Figure 4). To make the analysis computationally feasible, we stochastically approximate the next-generation matrix (NGM, see Methods). We fix all parameters except for the mutation rates (and), and the maximum virus load (Vmax). We keep the ratios and between the mutation rates constant (see Table 1 for the parameters chosen). Apart from the standard model described above, we also consider two modifications that serve as controls. In the first control, we take out the effect of population-level selection for transmission. In the second control, we make the population homogeneous. Between-host adaptation is only possible if spVL is inherited from one person to the next. If the speed of within-host adaptation is intermediate or fast, our model does not predict population-level adaptation for transmission. Here we show that the absence of between-host adaptation is not due to lack of spVL heritability (h2, see Methods). To this end, we compute heritability during an epidemic (see Figure 3, bottom panels), and in the steady state of the (standard) model for many different parameter combinations (see Figure 5C). During a simulated epidemic, we use all transmissions that take place within a time span of a year to compute heritability. This means that the sample size for computing heritability equals the (yearly) incidence. The median incidence in the simulation with a low mutation rate (, Figure 3B) equals 2335 yearly infections (2. 5th, 97. 5th percentiles: [1839,4085]). For the simulation with a faster mutation rate (, Figure 3A), more virulent viruses evolve, and the median incidence equals 3209 infections per year (2. 5th, 97. 5th percentiles: [241,3586]). Even with these large sample sizes, heritability fluctuates substantially. In Figure 3B the median of h2 is 3. 12% (2. 5th, 97. 5th percentiles: [−1. 31,8. 29]), and in Figure 3A the median of h2 is 17. 1% (2. 5th, 97. 5th percentiles: [10. 8,21. 8]). The rapid fluctuation in h2 might explain why different experimental studies to HIV-1 spVL heritability that use transmission couples [7], [11]–[14] often give quite varying results [cf. 9]. The NGM approach allowed us to produce an even larger number of transmission couples, and hence, to estimate heritability more accurately. Overall, heritability lies between 0% and 30%, and is for realistic parameter combinations. In our model, one can think of two mechanisms that cause heritability. The first mechanism applies when mutation rates are not too high. If variation in the number of mutations exists and the mutation rate is low, the spVL of transmitting and recipient hosts are correlated, although this correlation will not be perfect due to the variation in the breadth of the immune response (k). If the mutation rate increases, viruses adapt to their host more rapidly and, according to this first mechanism, the correlation vanishes. The second mechanism is related to transmission of crippled viruses. If a host controls the infection well because of a broad immune response, the virus will escape more CTL responses and, when transmitted, becomes crippled in the average new host. In the primary, controlling host, the set-point virus load is low due to good initial immune responses and the virus' fitness cost of escape, and in the average secondary host the virus load will again be low due to the high number of deleterious mutations. Vice versa, in hosts with a narrow immune response, transmitted strains will have few new escape mutations and this will lead to few deleterious mutations in the recipient. We can most clearly see the effect of the second mechanism when both mutation rate and Vmax are high (the contour in Figure 5C). In this part of the parameter space, most immune responses are escaped in the acute phase (cf. the solid graph in Figure 5B). Rapid escape causes variation in the number of deleterious mutations in the transmitted virus, because the size of the binding repertoire (k) varies among individuals. However, when, not all deleterious mutations can be reverted in the acute phase (cf. the dashed graph in Figure 5B). For high Vmax, the asymptomatic phase is short, resulting in few reversions during this phase and a ‘footprint’ of the transmitting host' s immune responses on the receiving host' s spVL [11]. Notice that the second mechanism does depend on the reasonable assumption that reversion is a slower process than escape (and, not shown), and that the size of the binding repertoire (k) differs between individuals. As the above evidence for the second mechanism—or ‘footprint effect’ as we like to call it—is only circumstantial, a quantification of this mechanism is needed. To quantify the footprint effect we analyze the simulations using a structural equation model (SEM). The model estimates heritability, and takes the fitness costs (m = e+f) and breadth of the immune response (k) into account. Heritability of spVL is the sum of two effects; one mediated by viral fitness, and the other by the breadth of the immune response of the transmitting host. Figure 7 shows a graphical representation of the model, and details of the analyses are given in the Methods section. For realistic parameter values, approximately half of the observed heritability is due to the footprint effect (Figure 5C). When we lower the rate of escape, the footprint effect, and therefore also the heritability, decreases. On the contrary, when within-host evolution is extremely fast, almost all of the heritability is due to the footprint effect, although the total heritability decreases. Our model predicts that heritability of the set-point virus load and host-heterogeneity are related. When within-host evolution is fast enough, approximately half of the observed heritability may be explained by the immunological footprint. Also, when we lower heterogeneity in our model, heritability decreases. An intuitive measure for heterogeneity in the host population is the expected similarity of hosts’ binding repertoires. This tells us how much adaptation to one host remains beneficial in the next. As a measure of the similarity of two binding repertoires K1 and K2 (of size k1 and k2, respectively) we use the Jaccard index, the overlap between binding repertoires, divided by the the number of (wild-type) epitopes that at least one of the hosts can recognize. Figure 8A shows the relation between the expected similarity between hosts () and the heritability of the set-point (h2). We modulated heterogeneity by varying n, the total number of potential epitopes, between 30 and 300, corresponding to low and high host-heterogeneity, respectively. The mutation rate equals 3y−1, such that the number of escape mutations during the acute phase lies within a realistic range. The figure shows that heritability indeed decreases when the population becomes more homogeneous, which indicates that high heritability relies on host-heterogeneity. HLA-heterogeneity differs between human populations. If our model prediction holds, then this variation could affect the heritability of the set-point measured in these populations. An unpublished study by Hodcroft et al. [44], shows that heritability in measured for HIV-1 clade C in a Sub-Saharan African population is higher than heritability for HIV-1 clade B in a European cohort [10] (30% vs. 5. 7%). Keeping our model in mind, we are able to understand this, if the European population with respect to clade B, is less heterogeneous than the Sub-Saharan population with respect to clade C. Using the peptide-MHC binding predictor NetMHCpan [45], we compared the two populations and circulating viruses (see Methods). Again, we measured heterogeneity in terms of similarity between binding repertoires. We sampled from the HLA-haplotype distributions of the European and Sub-Saharan populations, and calculated how similarity (again measured in terms of the Jaccard index) within these populations is distributed. Figure 8B shows two of these distributions. The black bars correspond to the European population, and the gray bars to the Sub-Saharan population. Although small, these populations do show a difference in heterogeneity: The Sub-Saharan population is more heterogeneous than the European population, as European binding repertoires tend to be more similar. The difference in heterogeneity is statistically significant (see Methods and Figure 8C).
In this paper, we model HIV-1 transmission and within-host adaptation by means of immune escape in a HLA-heterogeneous host population. In comparison to what data suggests, we do not find that HIV-1' s life history trade-off determines or influences spVL evolution. For realistic mutation rates, the evolutionary outcome is mostly determined by within-host selection for escape and reversion. Without HLA-heterogeneity, viruses would evolve to be within-host optimal in every host. Due to HLA-polymorphism, however, deleterious mutations accumulate, and the environment changes at each transmission. This causes virulence to evolve to intermediate levels for most hosts. Incomplete adaptation at the individual level is not exploited by the virus in order to improve it' s transmission potential. Although set-point virus loads are expected to be lower in a heterogeneous population, spVL evolution remains short-sighted. As we will point out below, our model is limited in the sense that we only incorporate immune escape and reversion as a means for within-host and between-host adaptation. Nevertheless, since population-level adaptation does occur when within-host adaptation is slow, the model' s limitations do not necessarily revoke our conclusions. In our model, we do find that spVL is heritable, even when the mutation rate is high. spVL heritability is needed for between-host adaptation. However, for realistic mutation rates, high heritability, as measured using transmission couples, is over-estimated; it mostly results from a ‘footprint’ left by the transmitter' s immune system on the receiver' s spVL. This novel explanation calls the validity of the use of high heritability as support for between-host adaptation in question. During real HIV-1 infections, immune escape sometimes requires compensatory mutations. Such escape variants need more time to revert to the wild-type in hosts lacking the escaped CTL response [46]. Such a mechanism is not incorporated in our model, but is likely to cause even higher heritability compared to what we find. Given previous results on the effects of transmitted CTL escape mutations on a receiver' s virus load [29], [33], and the sharing of HLA alleles [11], [47], we think the footprint effect provides a sound explanation for the experimentally observed high heritability of the set-point. Importantly, if this explanation were to be found true, and if spVL evolution and heritability are indeed strongly influenced by CTL escape, reversion and compensatory mutations, finding SNPs in HIV-1' s genome that control spVL might be a fool' s errand, unless this pursuit would be restricted to known CTL-epitope sites. Our claims concerning the footprint effect, and the dependency of heritability on host-heterogeneity are not just speculative. We show that the model can make testable predictions, and we give an example of how such a test can be performed, i. e. , by comparing host-heterogeneity in different human populations. In our example, we compared Sub-Saharan and European populations with respect to the viruses circulating in these populations, and showed that host-heterogeneity is higher in the African population, which is consistent with our novel explanation, and estimates of the heritability in these populations. Of course, we would not suggest that this isolated finding is evidence for the footprint effect, although we do want to stress that heritability estimates are expected to be correlated with host heterogeneity. Moreover, the heritability estimates that were used in this example were obtained using a phylogenetic analysis [10], [44], while our explanation only holds for studies that use transmission couples. In future work, we plan to investigate whether an immunological footprint can also affect heritability that has been estimated using phylogenies or pedigrees. Intuitively, the fact that within-host adaptation overrules between-host adaptation can be understood by considering that many viral generations separate the founding virus and a transmitted strain, while transmission only takes one generation. In the homogeneous model, this results in full within-host adaptation (throughout the population all epitopes are escaped), except when within-host adaptation is extremely slow. This result was also shown recently by Lythgoe et al. [21]. The intuition mentioned above works best for homogeneous populations. Adaptations to a primary host are beneficial again in a secondary host, and if within-host adaptation is fast, this leads to population-wide within-host adaptation and not between-host adaptation. This part of the intuition fails for a heterogeneous host population, where within-host adaptations in the form of immune escapes, are most likely not beneficial in the next host. Therefore, one could argue that homogeneity obstructs between-host adaptation. Here, we attempt to remove that obstruction by adding host-heterogeneity to a multi-level HIV-1 model. We find that in a heterogeneous population, HIV-1 also fails to evolve a mean spVL that maximizes the transmission potential, as shown by our sensitivity analysis and controls. Of course, when we make within-host adaptation trivial by choosing a very low mutation rate, population-level adaptation occurs. Apparently, host-heterogeneity does not solve the within- versus between-host adaptation paradox. Our models tell us that within-host adaptation overrules between-host adaptation, and yet HIV-1 appears to have adapted with respect to the life history trade-off [17], or at least is evolving its mean spVL towards the value that maximizes the transmission potential [48]. Several mechanisms that can serve as a solution for the paradoxical observation have been proposed [49]–[51]. One of these mechanisms is referred to as ‘store and retrieve’ [49]. It is hypothesized that latently infected memory CD4+ cells occasionally produce virus, and that these virions are preferentially transmitted. Preferential transmission is backed up by the observation that evolutionary rates are higher at the within-host than at the between-host level [52], and recently by a very interesting study into HIV-1' s transmission bottleneck [53]. However, transmission of CTL escape mutants within transmission couples [29], and even the spread of CTL escape mutants through populations has been observed [38], [54]–[57]. These observations indicates that ‘store and retrieve’ is not absolute, and in order for this mechanism to solve the paradox, we expect it to rely on getting the population-level evolutionary rate below a threshold; one which may not be reached. This premise could limit the robustness of the ‘store and retrieve’ model. Furthermore, when the population-level evolutionary rate is slowed down because of a mechanism like ‘store and retrieve’, the rate of between-host adaptation is also decreased, which could conflict with the short time scales at which adaptation must have been taking place for HIV-1 [18], [20]. Another possible mechanism is a heritable viral trait that influences spVL, but that is not under within-host selection. This was recently examined by Hool et al. [50]. An example of such a trait could be target cell activation rate [58], [59]. In short, if a viral trait influences target cell activation, and a mutant strain manages to increase the activation rate, then this additional activation is a ‘common good’ for the entire within-host viral population (activated cells produce more virions). Hence, the mutant does not have an advantage and will not be preferentially selected. Drift creates within-host variation in activation rate, and the transmission bottleneck leads to variation of the target cell activation rate at the population-level. This hypothesis could be challenged by other traits that affect spVL, since these may still be under within-host selection, and are likely to interfere with the within-host neutral one. We finish with a novel suggestion for solving the paradox, one which is based on our modeling formalism, and was recently also put forward by Fraser et al. [51]. One point of criticism on our model could be that we limit the evolutionary capabilities of our in-silico viruses. Strains can only evolve their number of deleterious mutations in order to approach population-level favorable spVLs. Unfortunately, in the current framework, it is not sensible to allow for mutations in other parameters, in particular Vmax, since then Vmax would only increase during within-host evolution, and hence, during the course of the epidemic. This is because we assume that no two strains simultaneously reside a single host, and that mutants with a higher fitness go to fixation rapidly. In reality, fixation of mutants within a host can take a considerable amount of time [60]. An obvious—but technically challenging—fix for this problem is to abandon the assumption that the within-host evolutionary dynamics is memoryless, and allow for multiple mutants to compete for fixation, i. e. , allow for clonal interference [61]–[65]. These mutants can then carry negative fitness effects (e. g. , Vmax decreasing mutations) along with beneficial escape mutations or reversions (genetic hitchhiking). Additionally, mutants with a small Vmax increasing effect, but that are otherwise equal to the wild-type, may have a long fixation time and can easily be out-competed by, e. g. , escape mutants. This makes within-host Vmax evolution more selectively neutral, and hence more sensible in our model. In future work, we aim to test if these speculations are valid, and whether a more detailed within-host fitness and selection model can unify within-host evolution and population-level adaptation.
In general, whenever a new event E is created during the simulation at time t, the exact moment when E will take place is unknown. Therefore, we assign to E a threshold and a load. The threshold is sampled from some probability distribution with non-negative support and mean 1. We first compute the waiting time, while conditioning on E being the first event to take place: Here, is the ‘rate’ or ‘hazard’ at which E takes place, which can depend on time. Notice that the (conditional) waiting time could be infinite (e. g. , when the number of susceptibles equals zero, the first event to take place can never be a transmission). When we perform this computation for all future events E, we find the event F that must take place first, and also the time at which it takes place, i. e. , . We then perform the following steps: First, we update the time. Then, for all future events E, we update the load as follows: Finally, we let the event F act on the current state and remove F from our event list. For instance, if F happens to be a transmission event, we should initiate a new host and decrease the number of susceptible individuals. Additionally, new transmission events should be created for the transmitter and recipient. Hereafter, we re-compute the waiting times for all events E and repeat the above steps. In most cases, the computation of wE and the updating of αE is simple. For instance, if E is a reversion event and f>0, then. For updating the load, we replace αE with. A transmission event requires more effort, because the rate of transmission varies during an individual' s infectious lifetime. The model can now be described by specifying for the events E listed in Table 2, their threshold-distribution, their ‘rates’, and the precise actions on the state (see Table 3). In order to study the steady state of the above described model, we developed a faster and more accurate method. In deterministic (e. g. ODE-based) models with multiple viral strains, one can compute the next-generation matrix (NGM), using the model' s rate equations [67]. Given a ‘generation’ (i. e. , a distribution of strains in a cohort of newly infected individuals), the NGM gives the ‘next generation’ after mutation and selection in a discrete generation-based model. The steady states of the original (continuous time) and generation-based model coincide. This steady state can be computed by finding the dominant eigenvector of the NGM. The dominant eigenvalue equals (by definition) [68]. Our model is not deterministic, but we can approximate the NGM using a Monte-Carlo method. We start with a virus that has m1 mutations. We then infect a large cohort (of size N) of individuals. These individuals may have different binding repertoires (of diverse size k), so we first sample pairs (e1, f1) with and. Then we run a within-host simulation for each of the virus-host pairs. Finally, we sample strains (e, f) that would be transmitted by the hosts at the start of an epidemic, and we count the number of transmitted stains that have mutations. The vector with approximates the m1-th column of the NGM. If the sample size N is large enough, the dominant eigenvalue and corresponding right eigenvector of the matrix approximate, respectively, and the steady state distribution of prevalent viral strains in our agent-based model. By sampling strains from the steady-state distribution, and simulating infections with these strains, we can compute statistics as in equilibrium. This method is not based on formal arguments, but below we put forward some heuristic evidence for its correctness. For the statistic heritability (h2), the above scheme is insufficient. However, we do have a cohort of potential transmitters, and hence we can create transmission couples by first sampling transmitted strains from the cohorts’ individuals, and then infecting recipients. The statistic h2 is computed as the slope of the regression between the spVL of transmitters and receivers. Classically, heritability of a trait x is defined as the proportion of variance in x that is caused by inherited genetic factors [see e. g. 18]. Hence, if we write, where is a genetic, and an environmental factor, then. The slope of the regression mentioned above is an estimator for this quantity, but only when the transmitted quantity in the recipient is independent of the the transmitter' s environmental factor. Below we will see that such an independence assumption does not hold for our model, and that the use of transmission couples results in an over-estimate of spVL heritability. To quantify the effect of the immunological footprint on heritability, we use a structural equation model (SEM), depicted as a directed, acyclic graph (DAG) in Figure 7. In our model, the actual inherited quantity is the number of mutated peptides. During an infection this quantity can change due to escapes and reversions, so we will only incorporate the number of mutations at the moment of infection (mtra for a transmitting host, and mrec for the corresponding receiver) in our statistical model. The set-point virus load of the receiver (spVLrec) depends on mrec, and the breadth of the immune response against the wild-type virus (krec). Of course, more factors determine the set-point virus load, such as the initial number of escape mutations, and stochastic effects such as mutations and progression to AIDS, but the simplified SEM only contains the variables spVL, m and k. Apart from ktra, the breadth of the transmitter' s immune response and mtra, the set-point virus load of the transmitter (spVLtra) depends also on mrec. This is because the set-point is an average over the chronic phase, and hence, the transmitted virus co-determines the set-point of the transmitter. In Figure 7, this is indicated by the arrow. During infection of the transmitter, the virus escapes a number of immune responses, and this number is dependent on ktra. This means that ktra influences the number of mutations of the transmitted virus mrec. This ‘immunological footprint’ is represented by the arrow in Figure 7. The breadth of the immune response ktra has no direct effect on mtra, since mtra corresponds to the transmitter' s founder virus. Likewise, there is no direct effect of krec on mrec. We use the the R package lavaan [69] to fit the model to (standardized) simulated data, that were produced using the NGM method and the standard model' s parameters. As an example, the result of one of such fits is given in Figure 7. In this graph, the numbers above the arrows indicate the estimated weights. The maximal virus load Vmax equals copies per ml, and the mutation rate equals 3y−1, such that the mean set-point for this population is 4. 51 log10ml−1 (cf. Figure 4A). Despite the large sample size of 25690 transmission couples, and the fact that the SEM has 4 degrees of freedom, the model describes the data quite well (the -test' s p-value equals 0. 81, and the root mean square error of approximation (RMSEA) equals 0 with a 90% CI of [0,0. 01]). In the context of our SEM, the statistic h2 equals the correlation between spVLrec and spVLtra. This correlation can be computed as the sum of the contributions of all paths that connect spVLtra with spVLrec. The contribution of each path equals the product of the coefficients along the path. The 3 paths that connect spVLtra with spVLrec are: where P3 is responsible for the immunological footprint. In the example of Figure 7, the contribution of P3 equals 0. 082, which is about half (49. 7%) of the total correlation between spVLrec and spVLtra (i. e. , of the heritability). We refer to the contribution of the path P3 as the “contribution of the immunological footprint to heritability”. We downloaded representative sequences for clades B and C from LANL' s HIV sequence database (www. hiv. lanl. gov; four sequences for each clade, as described in [70]). Then, we downloaded the HLA-A and HLA-B distributions for Europe and Sub-Saharan Africa from the NCBI database dbMHC (www. ncbi. nlm. nih. gov/projects/gv/mhc, [71]). Using the MHC binding predictor NetMHCpan (version 2. 4 [45]), we computed binding affinities of all 9-mers from the representative strains for the most common HLA alleles (covering 95% of the populations). For each HLA molecule, the binding threshold was chosen such that the top 1% of a set of 105 naturally occurring peptides would be considered a binder (as described in [72]). For our analysis, we sample pairs of HLA haplotypes from the HLA distributions of one of the populations (ignoring linkage disequilibrium), each haplotype consisting of two HLA-A alleles and two HLA-B alleles. For each two haplotypes, we then compare the similarity of the binding repertoires with respect to one of the four representative strains. As a measure of similarity, we use the Jaccard index (J): the size of the intersection, divided by the size of the union of the two binding repertoires. This gives us the distribution of similarity scores of a population with respect to a strain. Figure 8B depicts two of these distributions. The black bars correspond to the European population with respect to a clade B virus, and the gray bars to the Sub-Saharan population with respect to a clade C virus. By comparing the similarity distributions of a Sub-Saharan with a European population (Figure 8B), we can assess the difference in heterogeneity between the two populations and clades. The right panel of Figure 8C depicts the medians (one value for each representative strain). The European medians are significantly higher than the Sub-Saharan medians. For a better comparison between two distributions, we use a U-statistic, defined as, where Jeur and Jafr are distributed as the European and Sub-Saharan similarity distributions, respectively (cf. the Mann-Whitney U-test). Hence, U equals the likelihood that a random haplotype pair in the European population shows more similarity than a random pair in the Sub-Saharan population. We have four clade B strains and four clade C strains, and hence we can compute 16 probabilities U (Figure 8C, right panel). They turn out to be significantly higher than 0. 5, meaning that the European population, subject to clade B strains, is less heterogeneous than the Sub-Saharan population and clade C strains. We model within-host escape and reversion with two Markov chains: Let and denote the probability at time t that during infection phase i the host is infected by a virus with e escape mutations and f deleterious mutations, respectively, given that phase i started with an -virus at time. These probabilities satisfy the Kolmogorov forward equations [see e. g. 73]. Closed-form expressions for and are given by The probability that the host is infected with an (e, f) -virus only makes sense if we condition on the infection still being in phase i. We want to get the expected number of transmitted virus of a specific type, and in order to make the calculations possible, we take exponential distributions for the length of the phases. We tested that this assumption is not crucial by considering Erlang distributions. The rate at which phase i ends is given by. We also assume that mutation during the asymptomatic phase is slow and that the spVL is determined by the virus at the end of the acute phase (which is of type). This means that and can be kept constant. Furthermore, the fraction of susceptibles () can be kept constant, either because the population is in a steady state, or because the epidemic has just started (). Consider the probability generating function [cf. 74] for the number of transmitted virus of type during phase: Assuming that, we can write this integral in terms of the Beta function (). First we substitute the above given expressions for and and when we now assume that, we can get which equals by definition. Since one of the arguments in this Beta function is an integer, the function is rational: where we use the (rising) Pochhammer symbol. Now that we have this expression for, we can exploit the probability generating function' s useful properties. The number equals the probability that at the end of phase, the host is infected with an -virus. The expected number of transmitted -strains during phase equals. We use the following notation: If we now take into account that a transmitted -virus has a different phenotype in the receiver (with probability given by the hypergeometric distribution), we can find the NGM for the case. We verified that for this part of the parameter space (i. e. ,), the deterministic and stochastic computation give the same results (not shown). When the host population is homogeneous (), we find a threshold in the parameter space across which between-host adaptation is no longer possible. Here we will make this precise and show that this threshold is caused by a transcritical bifurcation. In a homogeneous population, we lose deleterious mutations. In the above introduced notation, we may ignore and we write for instance. Let denote the NGM, then we get the following formula in terms of and: The matrix is triangular, since the number of escape mutations, which equals the total number of mutations, can only grow during an infection. The diagonal elements of are the eigenvalues of, and the dominant eigenvalue equals (by definition) of the quasi-species. The diagonal elements can be written as If is dominant, then population-level evolution will result in strains that have escaped all CTL responses. If another eigenvalue with is dominant, then not all viruses have escaped all CTL responses and this is due to selection for transmission on the population-level. If we now fix and let approach from the right, then for high the eigenvalue is dominant. The mentioned bifurcation occurs when equals one of the (with) for the first time. We first give simple expressions for and that occur in the expression for: These expressions and the formula for enable us to (numerically) find the curves for. These curves and the resulting threshold are shown in Figure 6.
HIV-1 is a relatively young virus, being introduced in the human population somewhere between 1884 and 1924. Yet, previous studies suggest that the virus has already evolved to be efficiently transmitted among humans. Efficient transmission occurs when the set-point virus load, the semi-stable number of virus particles in the blood during the asymptomatic phase, is intermediate (approximately particles/ml). At this virus load level, individuals remain asymptomatic for a long period (7. 0 years on average), and still remain sufficiently infectious. In this study, we model the combined immunological and epidemiological dynamics of HIV-1 to explore whether population-level adaptation is feasible. We show that strong selective forces within the host are expected to dominate the much weaker population-level selection, unless the within-host dynamics of immune escape becomes exceedingly slow. Surprisingly, our analyses yield high levels of set-point virus load heritability, as observed in human populations. In the model, heritability of set-point virus load partially results from an immunological ‘footprint’ of the host-virus interaction in transmitting hosts, affecting the receiving hosts' virus load.
biotechnology
medicine and health sciences
infectious disease epidemiology
immunology
microbiology
synthetic bioengineering
immunodeficiency viruses
protein engineering
plant science
sexually transmitted diseases
population modeling
plant pathology
evolutionary adaptation
bioengineering
infectious diseases
major histocompatibility complex
theoretical biology
medical microbiology
hiv
epidemiology
aids
microbial pathogens
infectious disease modeling
immune response
macromolecular engineering
directed evolution
clinical immunology
viral pathogens
biology and life sciences
viral diseases
computational biology
evolutionary biology
evolutionary processes
Immuno-epidemiological Modeling of HIV-1 Predicts High Heritability of the Set-Point Virus Load, while Selection for CTL Escape Dominates Virulence Evolution