Translating HIV sequences into quantitative fitness landscapes predicts viral vulnerabilities for rational immunogen design

Abstract

A prophylactic or therapeutic vaccine offers the best hope to curb the HIV-AIDS epidemic gripping sub-Saharan Africa, but it remains elusive. A major challenge is the extreme viral sequence variability among strains. Systematic means to guide immunogen design for highly variable pathogens like HIV are not available. Using computational models, we have developed an approach to translate available viral sequence data into quantitative landscapes of viral fitness as a function of the amino acid sequences of its constituent proteins. Predictions emerging from our computationally defined landscapes for the proteins of HIV-1 clade B Gag were positively tested against new in vitro fitness measurements and were consistent with previously defined in vitro measurements and clinical observations. These landscapes chart the peaks and valleys of viral fitness as protein sequences change and inform the design of immunogens and therapies that can target regions of the virus most vulnerable to selection pressure.

Figure 1

Cartoon schematic of a viral fitness landscape. The replicative fitness of a viral strain is a function of its amino acid sequence. This information can be visualized as a topographical map where the amino acid sequence specifies a location on the landscape, and the height of the landscape prescribes viral fitness. For visualization purposes, this cartoon pertains to a virus consisting of only two residues. For multi-residue viral proteins, the fitness landscape is traced out in higher dimensions. The broken line indicates a hypothetical high-fitness mutational escape pathway from the global fitness maximum, to a nearby local maximum.

Figure 2

Comparison of our theoretical metric of fitness (E) and experimental in vitro replicative fitness data. (A) In vitro replicative fitness was measured for 19 Gag single and double mutants (cf. Computational and Experimental Procedures). A Pearson correlation coefficient of ρ = −0.52 (p = 0.02) reveals a statistically significant negative correlation between the energy difference of the engineered mutants relative to the wild-type strain computed from our model, (E − E_wt), and the logarithm of the measured relative fitness of the mutant, log (f/f_wt). (B) Replicative fitness data was compiled for 25 engineered Gag mutants containing up to five point mutations (Brockman et al., 2007 (Jurkat cell), Miura et al., 2009, Troyer et al., 2009, Schneidewind et al., 2007, Schneidewind et al., 2008). This data also exhibits a strong negative correlation ρ = −0.81 (p = 2 ×10⁻⁷). In each panel, a linear least squares fit is provided to guide the eye, and error bars delineating estimated uncertainties in the relative fitness are provided where available.

Figure 3

Clinically documented mutant strains correspond to low energy (high fitness) states within our inferred model. Circles: Energies assigned by our model to the 10 p24 single mutants, 8 p24 double mutants, and 5 p17 double mutants listed within Table S1, which identifies the particular mutants, HLA associated epitopes, and computed energies. Squares: Energies of the strains of three p17 epitopes observed to undergo sequence adaptation in longitudinal deep sequencing of an HIV infected host (Henn et al., 2012). In all cases, residues outside the epitope are treated as wild-type in the assignment of energies.

Figure 4

Comparison of point mutant energy costs with documented B*57 associated escape mutations in the TW10 (Gag_240-249) epitope in p24. The 10 residues constituting the epitope are indexed along the abscissa, and the energy cost associated with making a point mutation at each position, ΔE, along the ordinate. The greater the energy cost, the higher the fitness penalty. Mutations at position 246 were not observed within our sequence alignments, leading to the specification within our model of an infinite energy cost of a point mutation at this position, which we denote by *. The observed escape mutations at 242 and 248 (Brumme et al., 2009) occur precisely at the two positions carrying the lowest energy cost (smallest fitness penalty).

Figure 5

Longitudinal deep sequencing within a single infected host identified three p17 CTL Gag epitopes as undergoing sequence adaptation during the first four years of HIV infection (Henn et al., 2012). In panels A–C, the title provides the location of each CTL epitope within Gag (e.g., Gag_28-36), the HLA association (e.g., A24), and the name of the epitope (e.g., KW9). On the left side we present bar charts illustrating our inferred h parameter (c.f. Computational and Experimental Procedures) at each position in the epitope. On the right, we list the strains observed by deep sequencing at the six time points (Henn et al., 2012), the fraction of the deep sequencing reads corresponding to each strain (Henn et al., 2012), and the energy of each strain assigned by our model. Red letters indicate point mutations relative to our MSA consensus; all residues outside the epitope are treated as wild-type in the computation of energies from our model. In panel D, we present the temporal adaptation courses tracking the energy assigned by our model to the mutant strains of each epitope sequenced at each time point.(A) The infecting strain contains two mutations in the KW9 epitope (Gag_28-36) at positions 28 and 34, corresponding to the positions with the second (h₂₈ = 1.50) and third (h₃₄ = 2.12) lowest h values and a compensatory J coupling (J_28,34 = −0.15, c.f. Computational and Experimental Procedures), lowering the energy of the double mutant relative to the two independent point mutations. By day 476, 61.5% of the viral population has reverted to the lower energy single mutant (h₂₈ = 1.50). By day 1543, 79.5% of viral strains have reverted to the lowest energy (E=0) wild-type state. The remaining 20.5% occupy the third least costly singly mutated state (h₃₄ = 2.12). The lowest (h₃₀ = 0.55), second lowest (h₂₈ = 1.50) and third lowest (h₃₄ = 2.12) singly mutated states are similar in energy.(B) The LY9 epitope (Gag_78-86) enters with a single mutation in the infecting strain at the single lowest energy position (h₈₂ = 0.17). No sequence adaptation is observed until day 1543, at which time 31.3% of the population has reverted to the lowest energy (E=0) wild-type. 57.1% of the population occupy the second lowest energy singly mutated state (h₇₉ = 0.97), while 6.3% occupy a low energy doubly mutated state (h₇₉ = 0.97, h₈₂ = 0.17, J_79,82 = 0.25), and 5.4% remain in the infecting state. (C) The GY9 epitope (Gag_71-79) entered with a single mutation at the second least costly position in the epitope (h₇₆ = 1.30), nearly equi-energetic with the least costly point mutation (h₇₉ = 0.97). By day 165, the wild-type strain transiently emerged in 30.8% of the population, but has vanished by day 476. At day 1543, 37.8% of the population remains in the infecting state. The remaining 62.2% possess mutations at the second least costly (h₇₉ = 0.97) and least costly (h₇₆ = 1.30) positions, and a small compensatory coupling (J_76,79 = −0.04). This double mutant is of very low energy, lying within the lowest 1% of all 8,646 possible p17 double mutants. Of the three epitopes considered, only for GY9 is a significant fraction of the day 1543 population not observed in the lowest energy wild-type state, although it transiently emerges in day 165. Of the three epitopes, GY9 was the only one against which a CTL response was reported (Henn et al., 2012), consistent with a situation in which the wild-type state is effectively less fit than a competing mutant strain capable of evading immune pressure.(D) Temporal adaptation courses follow high fitness (low energy) routes for each of the three epitopes: KW9 (red circles), LY9 (green squares), and GY9 (blue triangles).

Figure 6

Rank ordered bar chart of 121 p24 class I HLA epitopes according to the computed energy penalty, Δ, imposed upon the viral ensemble. The particular epitopes considered are listed in Table S2. Epitopes associated with a reported immunodominant response within a particular HLA class I allele (Streeck et al., 2009) are designated by colored bars and labeled with the epitope location and HLA association. Red bars indicate that the corresponding HLA allele has been linked with enhanced HIV control, whereas blue bars denote those that have not (Miura et al., 2009; Pereyra et al., 2010; Trachtenberg & Erlich, 2001). The immunodominance data pertaining to the B*57 KF11 epitope Gag_162-172were assumed to also apply to the epitope formed from its nine residue subset, Gag_164-172. Six of the seven immunodominant epitopes leading to the greatest fitness cost are associated with protective HLA alleles; the dashed line at Δ<E> = 1.54 represents the cutoff above which all immunodominant responses are associated with protective alleles.

Figure 7

Evaluation of immunogen candidates. (A–C) Scatter plots of (A) p6, (B) p17 and (C) p24 vaccine candidates in the three-dimensional design space spanned by: (i) the weighted average fitness impact in the target population, $\overline{\mathrm{\Delta }<\text{E}>}$, (ii) fraction of the target population that respond to at least one epitope in the immunogen (fractional coverage), and (iii) the number of components in the vaccine. The target population comprised the 21 most prevalent haplotypes in North Americans of European ancestry, accounting for 44.6% of this population. Respectively, 1, 1023 and 1,048,575 immunogen candidates were tested for p6, p17 and p24 (black crosses), of which 1, 25 and 44 candidates were located on the Pareto frontier (red circles). The composition of these Pareto efficient immunogens is listed in Tables S3–5. (D–F) Scatter plots analogous to those in panels A–C for (D) p6, (E) p17, and (F) p24, in which criterion (ii) was modified to the fraction of the target population that respond to at least two epitopes in the immunogen; 1, 31 and 62 candidates are located on the Pareto frontier. (G) Scatter plot of the 2,340 combination Gag vaccine candidates in the three-dimensional design space spanned by: (i) the weighted average fitness impact in the target population, $\overline{\mathrm{\Delta }<\text{E}>}$, (ii) fraction of the target population that respond to at least one epitope in the immunogen (fractional coverage), and (iii) the number of components in the vaccine. The compositions of the 95 Pareto efficient candidates are listed in Table S6.