From Predicting to Analyzing HIV-1 Resistance to Broadly Neutralizing Antibodies

Treatment with broadly neutralizing antibodies (bNAbs) has recently proven effective against HIV-1 infections in humanized mice, non-human primates, and humans. For optimal treatment, susceptibility of the patient’s viral strains to a particular bNAb has to be ensured. Since no computational approaches are so far available, susceptibility can only be tested in expensive and time-consuming neutralization experiments. Here, we present well-performing computational models (AUC up to 0.84) that can predict HIV-1 resistance to bNAbs given the envelope sequence of the virus. Having learnt important binding sites of the bNAbs from the envelope sequence, the models are also biologically meaningful and useful for epitope recognition. Additional to the prediction result, we provide a motif logo that displays the contribution of the pivotal residues of the test sequence to the prediction. As our prediction models are based on non-linear kernels, we introduce a new visualization technique to improve the model interpretability. Moreover, we confirmed previous experimental findings that there is a trend towards antibody resistance for the subtype B population of the virus. While previous experiments considered rather small and selected cohorts, we were able to show a similar trend for the global HIV-1 population comprising all major subtypes by predicting the neutralization sensitivity for around 36,000 HIV-1 sequences a scale-up which is very difficult to achieve in an experimental setting.


INTRODUCTION
To date, there is neither a vaccine nor a cure available for infection with the human immunodeficiency virus type 1 (HIV-1).With an incidence rate of around 2 million each year and 1.6 million deaths in 2012 (WHO, 2014), HIV-1 infections continue to be a major global health issue.Since humans seem not to have natural immune mechanisms to clear the infection, infected individuals need to receive lifelong antiretroviral treatment (ART).Due to the high mutation rate of the virus, drug resistances emerge frequently, often requiring a change of drug targets.However, the number of available drug target classes remains limited; this is why there is still a high demand for drugs addressing new targets.
A currently investigated treatment option is the passive transfer of a combination of broadly neutralizing antibodies (bNAbs) to HIV-1 patients.The advantage of these antibodies is that they are very broad and potent.The potency of an antibody is defined as the antibody concentration needed to inhibit HIV-1 infectivity by 50% (IC50), while the neutralization breadth of an antibody is measured by the ability of the antibody to neutralize viruses from different subtypes.Upon the advent of new single-cell based methods, an abundance of these new bNAbs has been isolated and their higher neutralization potency and breadth have been shown in several studies (Walker et al., 2009;Mouquet et al., 2012).The target of these antibodies is the HIV-1 spike, a trimeric heterodimer of two viral envelope glycoproteins, gp120 and gp41.The successful binding of an antibody to this spike blocks the two main functions of the spike, namely mediating host cell fusion and viral entry.As a consequence, the corresponding virus cannot infect any new cell.So far, there are five known epitopes on the envelope glycoprotein, which are targeted by a variety of bNAbs (given in brackets): on gp120 the CD4 binding site (e.g., VRC01, VRC-PG04, 3BNC117, NIH45-46) (Falkowska et al., 2012;Wu et al., 2010;Scheid et al., 2011), the V1/V2 region (e.g., PG9 and PG16), and the V3 GCB 2015 we built binary classifiers to distinguish between HIV-1 resistance and susceptibility to a bNAb based on the envelope sequence of the virus.Since there are differences in potency between the antibodies, it is also of interest to know how strongly a bNAb neutralizes the virus.For this reason, we also built regression models that directly predict the IC50 value from the envelope sequence of the virus, thereby enabling more fine-grained results.
For the learning process we used kernels in conjunction with large-margin based methods: support vector machines (SVMs) for the classification, and support vector regression for the regression analysis.The underlying kernel for both tasks should preferably fulfill three properties: to allow for positional uncertainty to account for the high mutation rate of the virus, to be able to identify consecutive patterns of amino acids reflecting the shape of some epitopes, and to learn multivariate signals in order to model the fact that epitopes might consist of residues that are not consecutive in the sequence.String kernels that capture positional information such as the oligo kernel (Meinicke et al., 2004) or the weighted degree kernel with shifts (WDKS) (Rätsch et al., 2005), match these requirements and therefore might lead to better performances than conventional kernels like the polynomial (Poly) or the Gaussian RBF kernel.To validate this hypothesis, we compared the performances of models based on each of these kernels.The comparison was conducted by 10 runs of a 5-fold nested cross-validation using AUC and Pearson Correlation Coefficient as performance measure for the classification and regression task, respectively.For the polynomial and Gaussian RBF kernel the amino acid sequences have to be transformed to a real-valued input.We used one-hot encoding to represent the sequence information for the polynomial kernel, i.e., each amino acid a i , i ∈ {1, . . .20} is transformed into a 20-dimensional vector, where only the i-th entry is 1, and the others are 0.For the Gaussian RBF kernel, we encoded the sequence information using physico-chemical properties based on Atchley et al. (2005) (RBF1) and Braun and Venkatarajan (2001) (RBF2).

Visualizing the samples' interrelations in the reproducing kernel Hilbert space
Since non-linear dependencies in the data can exist, disregarding them by simply using a linear method might lead to worse performances.Transforming the data from the linear input space in to a space H , in which those dependencies are better represented, can lead to a better separability of the data.Support vector machines that only need dot products of the samples can take advantage of those non-linear dependencies using kernels, which correspond to dot products in the space H .However, the interpretation of the learnt non-linear models and the predicted results remains a challenge, which might explain why non-linear SVMs are less often used than advisable despite their good performances.So far, few methods exist that address this disadvantage of non-linear SVM classifiers using graphical representation (Caragea et al., 2001;Wang et al., 2006).These methods are usually neither generally applicable (only for certain kernels or restricting the data to be low-dimensional) nor simple (requiring additional optimization steps).
In this study, we propose a general method that displays the interrelations between the training and the test samples in the reproducing kernel Hilbert space (RKHS) without explicitly using the feature mapping function Φ. Euclidean distances between the samples in the RKHS, representing the interrelations, can be expressed with the help of the kernel function (Shawe- Taylor and Cristianini, 2004) (1) To provide a user-friendly representation, we visualize the pairwise feature space distances in a three dimensional space using multi-dimensional scaling (MDS) (Kruskal and Wish, 1978).Multidimensional scaling preserves the between-samples distances to the magnitudes of the variables' interrelationships while projecting the data into a D-dimensional space.This way, highly similar variables are spatially closer.Analyzing the stress for MDS with D in a range from 1 to 10, reveals that a three-dimensional space is sufficient to represent the data for all bNAbs (data not shown).
Visualizing the feature space distances with MDS offers new information on the interrelations between the used training data in the RKHS.Furthermore, it can be used to investigate the representation of the test sample x with respect to the training samples x i with i ∈ {1, . . ., N} in the feature space.Upon the MDS step, we include the class label and the contribution of the training points to the classifier via a color scheme into the MDS visualization (cf.Fig. 2).The two classes were colored in two different colors (blue and orange) where the color intensity was assigned by α i k(x i , x) with α i being the weight of sample i in the classifier.Thus, the color intensity of each training point increases with growing similarity to the test sample as well as with larger influence on the classification result.

PrePrints
From Predicting to Analyzing HIV-1 Resistance to Broadly Neutralizing Antibodies The labels of the close-by neighborhood of the test sample provide the user additional information on the prediction result.
We present not only a general meaningful graphical representation of the complex feature space, but also a tool to further investigate and interpret the prediction outcome.

Identifying learnt discriminant signals of the classifiers
In general, the learnt signals of a classifier can be traced back, if the kernel incorporates positional information such as the weighted degree kernel with shifts (WDKS) (Rätsch et al., 2005) or the oligo kernel (Meinicke et al., 2004).Both kernels compute the similarity between two sequences x and x of same length L by comparing the co-occurrences of their substrings within a certain distance.While the oligo kernel considers only the substrings of length l with 1 ≤ l ≤ L, the WDKS takes into account the co-occurrences of every substring up to a length l, thereby adjusting for potential overlapping signals of lower order (≤ l).For the WDKS, the significance of each oligomer can be identified using positional oligomer importance matrices (POIMs) (Sonnenburg et al., 2008).Due to better performances, we used in this study the oligo kernel to build the the prediction models (cf.Section Results and Discussion).The oligo kernel defines each sequence by the occurrences of its l-mers, which are encoded via so-called oligo functions.The oligo function µ, encoding all occurrences p ∈ x ω of a particular l-mer ω in a sequence x, is defined as with the continuous position variable t ∈ [1, L] and σ 2 controlling the positional uncertainty.As described in Meinicke et al. (2004), the corresponding learnt weight of the classifier for each oligomer ω at each position t can be retrieved by where i ∈ {1, . . ., N} denotes the i-th training sample with α i ≥ 0 and y i ∈ {−1, 1} being the learnt weight and classification label of the i-th sample.

Understanding the classification result
Additional knowledge on the classifiers such as the provided interrelationships of the training and test data in the kernel feature space or the discriminant signals learnt by the classifiers, leads to interpretable prediction models.The classification decision for each test sample remains however elusive.In this paper, we offer for each classification of a test sequence, a motif logo -a representation of the test sequence -that shows those residues in the test sequence that contributed the most to the classification result.
Using the kernel feature representation of the oligo kernel, we retrieve the contribution for each residue of the test sequence x * to the classification result.Since the residue at a certain location t of the test sequence is fixed, there exists only one oligomer ω containing this residue as starting point whose contribution is calculated as with µ * ω being the oligo function of l-mer ω of the test sequence.For l-mers > 1 the computed contribution is assigned to all amino acids of the oligomer.Since showing the contribution of each residue of the test sequence might rather be confusing than improving the interpretability of the classification result, we limit the motif logos to only represent the most discriminant residues of the test sequence.This is possible, since we could show that models using only the strongest p% signals with p ∈ {1, 3, 5, 7, 10, 15, 20, 25} do not have a significantly worse performance compared to models using the full envelope sequence as information (data not shown).To generate the motif logos we used Weblogo 3.0 (Crooks et al., 2004).

Kernel preselection
To validate the hypothesis that string kernels such as the oligo kernel or the weighted degree kernel with shifts might be more suitable for this application than commonly used kernels such as the Gaussian RBF or the polynomial kernel, we compared the corresponding prediction performances with each other for both, the classification and the regression task.Here, we show the performances of the classifiers for the bNAbs PG9, PG16, 10-996, 10-1074, PGT121, VRC01, and VRC-PG04 as listed in Table 1.
In general, the classifiers using the string kernels yield not a better performance compared to the Gaussian RBF or the polynomial kernel.Only for the VRC-PG04 bNAb the classifier using the oligo kernel performed better than all the other kernels.A reason for this might be the characteristic binding site of VRC-PG04 on the envelope protein.The binding site of VRC-PG04 is rather a large consecutive sequence than single residues as for the others bNAbs.While all kernels are good at identifying single residues, only the oligo and the weighted degree kernel with shifts are able to capture substrings of length l (l-mers).In contrast to the oligo kernel, the weighted degree kernel with shifts counts all co-occurrences of substrings of length ≤ l, thereby, adding a lot of noise the higher the parameter l is set, if only the long k-mers are informative.Parameter settings for the models Upon kernel comparison, the oligo kernel was selected for all bNAbs to predict the neutralization susceptibility to each bNAb for unseen viral strains.Table 2 presents the final parameters settings for the classifiers for the bNAbs PG9, PG16, 10-669, 10-1074, PGT121, VRC01, and VRC-PG04 fitted by a 5-fold cross-validation (cf.Table 2) and the corresponding performances (cf. Figure 1).
For the PGT121 and VRC-PG04 classifier an l-mer of length 6 led to the best performance whereas the l-mer length for the other antibodies was comparatively small (2-mers for VRC01 and single positions for the remaining antibodies).The length differences of the l-mers for different epitope classes supports the knowledge gained from experimental findings.For the N-glycan dependent antibodies, a single glycan site is the most important residue for a successful binding.The N332linked (V3-loop-directed) antibodies PGT121, 10-1074, and 10-996 need in the first instance an asparagine at position 332 for successful binding (Julien et al., 2013).The N160-linked antibodies PG9 and PG16 bind in a hammerhead-like way to the virus, building contacts with two glycans (160 and 156 or 171) (Louder et al., 2011).For the CD4 binding site (CD4bs), which forms a cavity, it is only known that it is sterically not easy to bind to for antibodies (West et al., 2014) .Longer l-mers led to the best prediction results for the CD4bs classifiers which is likely due to the fact that the CD4bs-directed bNAbs target a larger epitope compared to the other bNAbs.

Identified discriminant signals
Using the oligo kernel properties as described in the Methods section, we examined the 15% strongest learnt signals for each classifier.Among this set of signals, several residues were learnt by the classifiers that are also supported by literature (Lacerda et al., 2013;West et al., 2012).In Table 3 we present the confirmed learnt signals of the classifiers for the bNAbs PG9, PG16, 10-669, 10-1074, PGT121, VRC01, and VRC-PG04 as an example.
Most of the found discriminant signals for the N-glycan dependent antibodies, that is, for the V1/V2-loop-and V3-loop-directed antibodies, contain the amino acids asparagine (N), serine (S) and threonine (T).These amino acids are also part of the pattern N-X-[S or T], which defines potential N-glycosylation sites (PNGS) (Marshall, 1974).The classifiers for the CD4bs antibodies identified known required residues for CD4-binding as reported in West et al. (2012).The fact that all classifiers learnt some known discriminant position, further support the reliability of the prediction models in addition to the provided prediction performances.Additionally to the already known epitope sites, we also found other discriminant residues that might be interesting for follow-up structural studies.

Application of the visualization methods
To demonstrate the introduced visualization methods, we retrieved several HIV-1 envelope sequences from the Los Alamos HIV sequence database (Foley et al., 2013) serving as test input for the classifiers.
We present here the test case for the sequence with the GenBank ID HM469973 and the PG9 classifier, which classified the sequence as susceptible.

Visualizing the samples' interrelations in the reproducing kernel Hilbert space (RKHS)
In general, the training samples form dense agglomerations (clusters) with respect to existent interrelations in the RKHS; exemplary shown for the PG9 classifier and the test sequence HM469973 in Figure 2 (fitted to 2-dimensional space for better representation).By adding subtype information to the plot, we observed that the clusters comprise mainly the sequences of the same subtype but not exclusively (data not shown).This is an expected finding, since the arrangement of the points is only based on the kernel similarities, which consider the whole sequence and not only the discriminant positions.Due to the coloring scheme, the most visible close-by point to the test sample (+) is also the most similar and influencing training sequence in the feature space offering the user more information about the representation of the training samples in the classifier as well as a better understanding of the classification result.

Motif logo for the test sequence
To provide the influence of each residue of the test sequence on its classification outcome, we proceeded as described in the Methods Section.In Figure 3 such a motif logo is presented for the test sequence HM469973 and the PG9 classifier.It can be observed that the test sequence has an asparagine (N) at position 160, which is known to be decisive for a successful binding of the PG9 bNAb.In all three logos, this pivotal residue has the most influence on the classification outcome compared to the other contributions.Considering the contribution of the 1, 3 and 5% learnt discriminant signals to the classification outcome of the test sequence, it can be seen that more discriminant signals occurred in the test sequence that are linked with neutralization susceptibility than with neutralization resistance.

PrePrints
German Conference on Bioinformatics 2015 (GCB'15) GCB 2015 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −1.

HIV-1 resistance trend analysis
Apart from predicting neutralization sensitivity of unseen viral strains, we used our models to investigate whether neutralization sensitivity of HIV-1 to bNAbs has changed over time.For subtype B variants, a continuous trend towards resistance has been already confirmed in certain cohorts of the French and Dutch HIV-1 population (Bunnik et al., 2008;Bouvin-Pley et al., 2014, 2013).Since evolving resistance towards antibody neutralization in the HIV-1 species would have major implications on the antibody selection for current vaccine development, it is important to know whether such a drift towards resistance also exists in the global HIV-1 population for all subtypes.In contrast to an experimental setting, where expensive neutralization assays need to be performed for a large number of viral strains, we can use our learnt prediction models to examine this question.
To model the global HIV-1 population over time, we used all available envelope sequences from the Los Alamos database (∼36,000) comprising viral isolates from all major subtypes over a time interval from 1981 to 2013.We divided the given time interval into 5 time periods to account for changes in HIV-1 treatment strategies : 1981: -1986 before ART; before ART;1987-1991 ART monotherapy; ART monotherapy;1992-1995 ART combination therapy (cART); 1996-1999 cART with protease inhibitors; 2000-2005 cART with Lopinavir/Ritonavir; 2006-2013 cART with Maraviroc/Raltegravir.The neutralization sensitivity of the samples to the 11 considered bNAbs was determined using our support vector regression models to predict directly the IC50 value.
Taking all available samples into consideration, we observed a general trend to resistance over time to all bNAbs except PG9 and PG16, for which the virus seems to become more susceptible (cf.Predicted neutralization sensitivity (log scale) for the ∼36,000 HIV-1 sequences from the Los Alamos database to the PG9 and NIH-4546 bNAbs using the regression models.While there is a trend towards antibody resistance for the CD4 binding site targeting bNAb NIH-4546, the neutralization susceptibility seems to increase for PG9.all bNAbs, but PG9 and PG16 again; for PGT121 and PGT128 the trend was not significant anymore.By predicting the coreceptor usage for all sequences with geno2pheno[coreceptor] (Lengauer et al., 2007), we detected an increasing ratio of R5/X4-capable viruses over the time periods.Together with the known R5-bias of PG9 and PG16 (Pfeifer et al., 2014) (better against R5, worse against X4), this might lead to the observed trend towards susceptibility.In general, we could confirm that HIV-1 variants of subtype B show a trend towards antibody resistance (Bunnik et al., 2008;Bouvin-Pley et al., 2014, 2013).By using our prediction models, we extended the analysis to the world wide HIV-1 population considering all major subtypes and thus, validating the hypothesis on a large scale, which is very difficult to achieve in an experimental setting.

CONCLUSION
In this study, we showed that neutralization sensitivity of new HIV-1 variants to broadly neutralizing antibodies (bNAbs) is predictable using existing neutralization assays.The performance of the prediction models for the 11 considered bNAbs motivate their use in the selection of a bNAb combination therapy as a recommendation tool.The credibility of the models were enhanced by the finding that the prediction models learnt important binding sites for the bNAbs only based on the envelope sequence.Hence, additional information such as structural binding site information is unlikely to boost the performance.We increased the interpretability of the models, by offering the user more information on the prediction outcome in form of a motif logo where the logo displays the contribution of the pivotal residues of the test sequence to the prediction.In addition, we introduced a new general method that enables to visualize the feature space interrelations of the SVM models, providing thereby more information on the SVM classifiers.
Apart from their potential use as recommendation tool, the models can be used to analyze the change in the neutralization sensitivity of HIV-1 over time.We could confirm previous results suggesting a trend towards antibody resistance in the subtype B population (Bunnik et al., 2008;Bouvin-Pley et al., 2014, 2013).Moreover, we scaled up the analysis to the global HIV-1 population, showing that there is a general drift towards antibody resistance in the world-wide HIV-1 population.These findings are relevant for the selection of suitable vaccine candidates; a combination of bNAbs is however still very potent in neutralizing HIV-1 (Bouvin-Pley et al., 2014).

Figure 1 .
Figure 1.The AUC performances for the best parameter setting for each bNAb classifier using the oligo kernel.

Figure 2 .Figure 3 .
Figure 2. Visualization of the interrelationships between the training and test samples in the PG9 classifier with the test input sample HM469973 projected into two dimensions using MDS.Training samples that could be positively neutralized by PG9 are colored blue, while neutralization resistant training samples are colored orange.The test input sample is displayed as a black cross.The color intensities of the training points increase with growing similarity to the test sample as well as with larger importance of the training sample to the classification outcome.The arrangement of the feature space distances is based only on the distances of the support vectors.
Figure 4. Predicted neutralization sensitivity (log scale) for the ∼36,000 HIV-1 sequences from the Los Alamos database to the PG9 and NIH-4546 bNAbs using the regression models.While there is a trend towards antibody resistance for the CD4 binding site targeting bNAb NIH-4546, the neutralization susceptibility seems to increase for PG9.

Table 1 .
Tested parameter ranges for each kernel together with the average AUC performance for each bNAb classifier in the nested cross-validation: For VRC-PG04 (in bold), the oligo kernel based classifier performed better than the others.

Table 2 .
The selected parameter settings for the oligo kernel classifier for each bNAb.

Table 3 .
Discriminant signals for each bNAb classifier that are supported by literature.Positive signals are colored in blue, negative in orange.