BACKGROUND
antibody-protein interactions play a major role in infectious diseases, autoimmune diseases, oncology, vaccination and therapeutic interventions.
antibodies present in human blood interact with antigens  with different affinities and in a sequence- and structure-specific manner.
when studying protein-antibody interactions, two types of epitopes are to be distinguished:  conformational and  linear epitopes.
in this study we focus on linear epitopes; see a recent review  <cit>  for a discussion of conformational epitopes.
all potential linear epitopes of a protein can be represented by short peptides derived from the primary amino acid sequence.
the binding site of an epitope covered by an antibody typically includes a minimal stretch of  <dig> to  <dig> amino acids.
if peptides of  <dig> amino acids in length are incubated with one specific antibody, that antibody will bind to its epitope independently of the physical position of the binding motif within the peptide.
motifs running from position  <dig> to position  <dig> up to motifs running from position  <dig> to position  <dig> would be possible.
this uncertainty results in difficulties for determining consensus binding sites as well as meaningful position weight matrices .
individual amino acids within epitope binding sites may have different impact on antibody recognition not only due to the nature of amino acids involved in binding  but also because of the specific position of the amino acid within the whole peptide sequence .
here, we present a method, pythia-design, for designing novel peptides with a desired binding affinity .
this method is built upon a successful, novel discriminative classifier called pythia  that can accurately label a given peptide as either a high- or low-affinity binder.
to test the quality of the designs that pythia-design produces, we experimentally constructed our designed peptides  and tested their binding affinity.
we show that pythia-design more accurately designs such peptides than barbarini et al.
<cit> .
we further show that pythia-design produces a more diverse set of designed peptides, which is important for generating a varied set for experimental construction.
finally, we show that the two methods of pythia-design and barbarini et al.
<cit>  can be combined, exploiting the relative strengths of both, to achieve even higher accuracy in epitope design.
while there is less prior work on epitope design , much previous work has focused on the task of predicting binding affinity of a given peptide to various target molecules  <cit> , e.g.
antibodies  <cit> , to mhc class i and class ii complexes alone or in concert with t cell receptor binding .
machine learning classifiers such as artificial neural networks  <cit> , hidden markov models  <cit> , and support vector machines  <cit>  and other approaches have been explored in tackling the problem of predicting human leukocyte antigen  binding peptides  <cit> .
much work has also focused on the prediction of t-cell and b-cell binding peptides .
zhao et al.
<cit>  explore various classifiers to predict peptide t-cell binding.
using a 10-dimensional feature vector to represent each amino acid, they discover that svms provide the best classification performance in their task.
huang and dai  <cit>  also explore the classification of peptide binding to t-cells using a support vector machine classifier.
they present a novel peptide feature based on combining a 20-dimensional indicator vector with amino acid similarity information encoded by the blosum <dig>  <cit>  matrix.
nanni and lumini  <cit>  introduced the mpps system that relies on an ensemble of support vector machines, trained on various physicochemical properties, to classify peptide binding to hiv-protease and t-cells.
they use sequential floating forward selection  <cit>  to select a subset of features, and combine the individual classifier predictions using the max rule  <cit> .
more recently, nanni and lumini  <cit>  have explored the use of a novel peptide-encoding scheme that relies on the use of nonlinear dimensionality reduction to extract the information encoded across a large number of physicochemical properties.
they demonstrate that this novel feature representation, when used in conjunction with a support vector machine classifier, exhibits state-of-the-art performance in predicting peptide t-cell binding.
wang et al.
<cit>  also showed that combining multiple classifiers using a consensus approach improved the classification of mhc class ii peptide binding predictions.
others  <cit>  have used motif mining for mhc i and mhc ii peptides.
recently, a flexible t cell receptor docking algorithm achieved near-native predictions for  <dig> % of the tcr/pmhc cases  <cit> .
zhang et al.
<cit>  use 3d features and a random forest classifier to predict b-cell epitopes.
lin et al.
<cit>  provide a method for b-cell epitope prediction that exploits phylogenetic information.
yao et al.
<cit>  introduce an svm approach called svmtrip for b-cell linear epitope prediction.
yao et al.
<cit>  compare various methods to find conformational epitopes of b-cells.
el-manzalawy et al.
<cit>  used support vector machines in combination with a subsequence kernel reaching a auc of  <dig>  and an accuracy of  <dig>  % to predict peptide/epitope-antibody-binding  <cit> .
our novel discriminating classifier upon which our design method is based uses an ensemble of support vector machines  to classify design candidates.
this classifier is broadly similar to that of nanni and lumini  <cit> .
however, we use probabilistic svms with platt’s extension  <cit> , along with a different set of features.
in addition, no other previous work deals with such a wide variety of paratopes — the regions of antibodies which recognize antigens — as is found in intravenous immunoglobulin fractions, as we do here.
as our experiments below show, a diverse set of ivig-binding peptides can be computationally designed using pythia-design with high accuracy.
in addition, our random-walk strategy for ensuring diversity is general and could be applied to any accurate discriminative classifier.
these computational techniques, and the collection of designed and validated ivig binders and non-binders, will be useful both for gaining a more thorough understanding of ivig binding properties and the diversity of possible epitopes.
methods
overview
the pythia-design method has two main parts.
the first part is a machine learning classifier that is trained to predict whether a given peptide is a high-affinity or low-affinity binder.
we designed this classifier, which we refer to as pythia, to use an ensemble of probabilistic support vector machines  trained on various sequence, chemical, and structural features.
the features used in the classifiers are described in .
for features where we did not compute the kernel directly, we used a radial basis function kernel.
the second part of pythia-design is a method for using random walks to generate candidate peptides with novel sequences to feed into the classifier.
this is described in section “generating novel peptides that bind to ivig”.
discriminative classifier for predicting binding and non-binding epitopes
we here describe a method for predicting whether a given peptide  is going to be a low-affinity or high-affinity binder.
motivated by the success of previous work in various protein-related prediction tasks  <cit> , we use support vector machines  as our classifiers within the ensemble.
we trained these svms via the libsvm software  <cit> .
for all features, the optimal svm parameters were discovered via a grid search  and 5-fold cross-validation.
the parameters were selected based entirely on a training set without the availability of the subsequent testing set, which was held hidden from method developers until after predictions were made.
for a given svm model, the cross-validation accuracy for the optimal set of parameters is used as a weight to combine the corresponding model’s predictions with the others from the ensemble.
during standard svm classification, instances are assigned a hard label, as belonging to the negative  or positive class , denoted by c− and c+ respectively.
such a hard labeling poses no problem when only a single classifier is used to label test data.
however, when an ensemble of classifiers is used, it is useful to have extra information about the degree to which the label assigned by each individual classifier should be trusted.
for this reason, we chose to use platt’s extension  <cit> , which provides probabilistic outputs for a support vector machine’s classifications.
instead of receiving a 0/ <dig> label, each instance is given an a posteriori estimate of the probability with which it belongs to the positive class.
thus, we expect that instances that clearly belong to the negative class will be given a value close to  <dig>  while instances that belong to the positive class will be given values close to  <dig>
having a probabilistic interpretation of the classification for data instances makes it possible to combine the output of different classifiers.
we used a variant of the sum rule, where the predictions of the individual classifiers are summed and normalized to yield the prediction of the ensemble.
specifically, the prediction of the ensemble for a particular instance xi was computed using
p+ensxi=1a∑j=1majp+jxi,
where xi is a feature vector representing the ith peptide, constructed using some subset of the features described in section “features used in the classifiers” , p+jxi is the a posteriori probability output by classifier j that the peptide with features xi is a high-affinity binder, and aj is classifier j’s cross-validation accuracy.
a is a normalization factor equal to ∑j=0maj.
we can then take p+ens to be the probability with which the ensemble predicts xi to belong to the positive class, or we can use it to obtain a discrete class prediction with the decision rule:
xi∈c+ifp+ensxi≥τ,c−otherwise.
in our experiments, we set τ= <dig> , but other values may be reasonable.
in fact, one may even learn the value of τ which yields the best performance by using a held-out subset of the training data, though we do not explore that here.
each svm model will yield a prediction for each peptide in the testing set.
we combined the predictions for all of the classifiers in the ensemble using a variation on the approach presented by nanni and lumini  <cit> , which is itself an extension of the sum-rule.
we normalized the predictions for each classifier to have a standard deviation of  <dig>  next, we combined the predictions from each of the j classifiers according to eq.
<dig>  by sorting the peptides in the testing set according to this value, we can produce a rank ordered list of the peptides in order of the likelihood that they belong to the positive  class.
features used in the classifiers
numerically encoded sequence features
there are two distinct types of sequence features that we encode numerically.
first, we used a simple variation on the peptide encoding scheme presented by huang and dai  <cit> .
we encoded each amino acid in the peptide by replacing its single letter code with its corresponding row in the blosum <dig> matrix.
the blosum <dig> matrix contains empirically derived log-odds scores that encode the frequency of different amino acid substitutions and is commonly used to measure the similarity between different amino acids.
let the peptide of length d be given as p=, where ai is the amino acid in the ith position of the peptide.
further, let row map the amino acid a to its corresponding row in the blosum <dig> matrix.
we encoded the peptide as enc=.
for the length d peptide p, enc will be a 20d dimensional feature vector.
in addition to blosum <dig>  we use the same type of encoding with matrices nlf and sa introduced by nanni and lumini  <cit> .
these matrices are derived by performing dimensionality reduction on a large, rectangular  matrix where each row corresponds to an amino acid and each column to some physicochemical property.
the goal of the dimensionality reduction is to decorrelate the physicochemical properties, reducing the column space of the matrix significantly.
the nlf matrix is a 20× <dig> matrix obtained using a nonlinear fisher transform, while the sa matrix is a 20× <dig> matrix obtained using a combination of clustering and principal component analysis.
the second type of sequence feature that we encode numerically involves various physicochemical properties of the constituent amino acids of each peptide.
we analyze the amino acid properties present in the amino acid index   <cit> .
each aaindex property provides a mapping from each of the  <dig> amino acids to a numerical scale measuring some physicochemical attribute .
the aaindex listed  <dig> different amino acid properties.
we use an approach based on a sliding window and histograms to turn each aaindex property into a numerical feature vector for a peptide.
consider a single aaindex property aaij, and let aaij represent the numerical value to which the amino acid a is mapped under aaindex property j. to form a representation for the entire peptide p under the property aaij, we used a window of length w sliding across the peptide to produce a -dimensional vector where entry i in this vector is the average value of the aaindex property over the window starting at position i. by varying w, we can change the coarseness of this representation.
through a process of experimenting with different values of w for this classification task , we computed these features for w∈ .
string kernel features
string kernels are used to evaluate the sequence similarity between peptides.
there are many different varieties of string kernels, ranging from the somewhat simple k-spectrum kernel, which essentially counts the occurrence of all length k substrings in each peptide, to the more complex substring-mismatch kernel  <cit> , which considers all shared subsequences between two peptides, allowing for gaps and mismatches.we use the k-spectrum string kernel  <cit>  for k= <dig> , <dig> , the sssk kernel  <cit>  with parameter d= <dig>  and the bounded range substring kernel  <cit>  with parameter r= <dig>  the output of each of these methods is a matrix, known as the kernel matrix, in which the entry at row i and column j is the result of the kernel evaluation between peptides i and j. to train a svm model for each of these string kernels, we simply compute the kernel matrix, and then make use of the ability of libsvm to train a model using a precomputed kernel.
structure features
we use rigid docking to estimate the binding affinity between ivig and a candidate peptide.
we computed a hypothesized 3-dimensional structure for each peptide using the biochemical algorithms library   <cit> .
we built a starting model for each peptide p by positioning the side chains for each amino acid by choosing the most frequently occurring rotamer position from a rotamer library.
we then optimized this initial structure by performing an energy minimization using the amber  <cit>  force field.
this relaxes the structure until a  energy minimum is achieved.
we also obtained an experimentally measured 3d structure for igg <dig>  <cit> , the most prevalent class of igg antibody present in intravenous immunoglobulin.
we measure the conformational complementarity of each of our hypothesized peptide structures with the immunoglobulin structure.
to compute this complementarity, we performed a protein-protein docking simulation for each of the constructed peptides against igg <dig> using the zdock software  <cit> .
each zdock run produces a list of the  <dig> top-ranked  docking predictions for each peptide.
the zdock score provides a measure of the complementarity of the peptide and immunoglobulin conformation in the docking region and is used as a proxy for the overall quality of the docking.
for each peptide, we formed a histogram from the  <dig> zdock scores, and use this histogram as a feature vector with which to train the svm model.
intuitively, we expect peptides whose zdock score distributions are skewed toward high scores to have better shape complementarity and, therefore, to be more likely binders than peptides whose zdock score distributions are skewed toward low scores.
generating novel peptides that bind to ivig
pythia-design builds upon the classifier described above by first generating many peptide sequences and then assigning them a reactivity category high , medium , or low  according to the predictions of our classifier.
we generated the de novo sequences using a sampling approach that corresponds to a seeded random walk in sequence space.
to obtain a sequence for reactivity class c, we choose a random seed sequence s∈c  from the training set, and randomly mutate its constituent amino acids until it adheres to several required sequence diversity rules, in order to ensure that the novel peptide sequences that are generated are sufficiently different from those in the training set.
specifically, the designed sequences predicted to be in the high- or low-reactivity category could not share any 4-mer, or exhibit a sequence identity of greater than  <dig> amino acids in any subsequence of length  <dig>  with any sequence in the same reactivity category in the training set.
further, the peptides generated by the pythia-design method were required to adhere to the same set of constraints with respect to the peptides in the original experiment’s testing set.
we seeded our peptide generator with  <dig> sequences from the training set —  <dig> sequences with the highest experimentally measured reactivity and  <dig> sequences with the lowest experimentally measured reactivity.
running the random walks then produced  <dig> candidate peptides that we classified using pythia and sorted according to their probability of belonging to the positive class.
there were  <dig> peptides with an probability greater than or equal to  <dig> , and  <dig> with an probability less than  <dig> .
the  <dig> sequences with the highest probabilities were predicted belong to the high-affinity class h, while the  <dig> sequences with the lowest probabilities were predicted to belong to the low-affinity class l. the remaining  <dig> peptides, with probabilities closest to  <dig> , were predicted to belong to the medium-affinity class m.
training and testing sets of peptides for classifier
to train and test our discriminative classifier, we used the data set of luštrek et al.
<cit>  , in which  <dig>  peptides were incubated with commercially available intravenous immunoglobulin  fractions, which was originally presented as part of the dream <dig> challenge  <dig> .
ivig is a mixture of naturally occurring human antibodies isolated from up to  <dig>  healthy individuals.
from this dataset, high-confidence negative and positive pools of recognized peptides were determined based on epitope-antibody-reactivities  of more than  <dig>  different peptides subjected to a huge number of structurally different antibodies present in ivig.
the training and test datasets for the discriminative classifier were assembled from these peptide pools .
from the collection of all the peptides incubated with human ivig, a pool of  <dig>  epitope-containing peptide sequences reactive with human immunoglobulins  was experimentally identified.
this was called the positive set.
from the same original collection of peptides  <dig>  peptides were identified that showed no antibody binding activity in any of the triplicate assays .
this peptide set was called the negative set.
the training set was formed by random sampling of  <dig>  peptides from the positive set and  <dig>  peptides from the negative set.
the training set thus created contained  <dig>  peptides and their respective binding reactivities.
the test set was created by joining together the remaining  <dig>  peptides from the positive set and the remaining  <dig>  peptides from the negative set, for a total of  <dig>  peptides.
selection of designed peptides for experimental validation
pythia-design was used to generate  <dig> peptides predicted to be reactive  as well as  <dig> predicted not to be reactive.
the method of barbarini et al.
<cit>  was also used to generate  <dig> peptides of each class.
to select a subset to experimentally construct and validate, the designed peptides were re-categorized by selecting  <dig> high binders and  <dig> non-binders as follows.
a stratified instead of a randomized sampling procedure was chosen in particular to investigate whether the designed peptides are robust to existing classification methods  <cit>  as well as position weight matrix  analysis.
the initial peptide set were subgrouped by using pwm and svm analyses as described in  <cit> ; see categories in legend of additional file 4: table s <dig> .
categories  were determined from a linear scale representing pwm and lustrek classification scores of the original data sets  from minmax sampling of the classifiers was used.
the stratified quota sampling was restricted to  <dig> peptide sequences in case of binding peptides .
in total,  <dig> peptides predicted to be bound and  <dig> predicted not being bound by ivig were taken from each of pythia-design and barbarini et al.
<cit>  and subjected to experimental testing.
RESULTS
overview
the test set of peptides was withheld from the algorithm designers until after the algorithm was finalized.
only the training set was used in the design and initial evaluation of the pythia classifier.
once finalized, pythia was then evaluated on the held-out test set of peptides to validate the classifier component of pythia-design.
finally, pythia-design was used to generate a collection of new likely binder and non-binder peptides, a subset of which were then experimentally validated.
validation of the pythia affinity classifier
although our goal is to produce a peptide design method, we first validate that our peptide classifier is accurate.
this classifier is based on an ensemble of learners that aggregates the prediction of many individual classifiers, each of which was trained on a set of features.
analyzing the auc for the individual classifiers in the pythia ensemble  on the test peptides, we observe that many of the classifiers show similar performance, with the exception of the structural classifier which displays significantly lower classification performance.
the ensemble, however, yields superior performance compared to any of its constituent classifiers with an auroc of  <dig> .
at any given false positive rate, the ensemble classifier will obtain a higher true positive rate than any of the other classifiers.
the precision-recall curve presents a related view of classifier performance to the roc curve.
it measures how the precision changes as the recall is increased.
for very small recall values , the sparse spatial sample and k-spectrum string kernels yield the best  precision.
however, for the vast majority of recall values, the ensemble classifier yields the highest precision.
just as was the case with the roc curves, the ensemble again achieved the maximum area under the pr curve.
while the auprs were generally lower than the aurocs, we did observe that the benefit of the ensemble was larger with respect to the pr curves than the roc curves.
local compositiona
<dig>
<dig>
athe best single classifier under both the auroc and aupr metrics
boldface indicates the best solution
because the test and training sets of peptides were chosen randomly, it is possible that overlapping or shared sequences between the test and training set partially leads to this high performance.
nevertheless, our ultimate aim is to design novel peptides, and this analysis suggests that pythia has reasonable performance identifying low- and high-affinity binding peptides.
accuracy of predicted affinity of computationally designed peptides
pythia-design was used to generate a number of probable high-affinity binding peptides , probable low-affinity binding peptides , and medium-affinity binding peptides .
in addition, the method of barbarini et al.
<cit>  was used to generate the same number of peptides.
a subset of designed peptides of both the pythia-design and barbarini et al.
<cit>  methods were experimentally constructed .
the binding affinities of these designed peptides were then experimentally measured.
additional file 4: table s <dig> gives the designed peptides and their measured affinities.
fig.
<dig> quality of designed peptides from two approaches.
the distribution of measured affinities for the designed peptides predicted to belong to the low  and high  binding affinity classes for the  pythia-design method, and  method of barbarini et al.
<cit> .
the horizontal line at  <dig>  indicates the binding affinity cutoff above which a peptide is considered to have a high binding affinity.
both methods produce a statistically significant separation of high- and low- binders , but pythia-design is much better at generating high-affinity binders
to quantify the degree to which the methods are able to generate binders and non-binders effectively, we use the method of ojala and garriga  <cit>  to compute a p-value indicating the probability of the observed high- and low-binding separation.
to do this, we compute the standard f1-score of the predictions, then randomly permute the labels  <dig> times to get a distribution of f1-scores.
for both pythia-design and pavia, the true f1-score is always better than the randomized score, meaning that for each method p< <dig> .
this indicates that the approaches are both truly designing peptides better than random.
using specific binding defined by a factor  <dig> above control measures  plus a minimum signal of above  <dig>  the designed peptides were categorized as truly or falsely predicted as outlined in additional file 4: table s <dig>  using these categorizations, we computed precision, recall, and accuracy for the two methods.
since we do not have a hard threshold to define high or low affinity measurements, and such threshold cannot be detected from the distribution of measurements itself, we used a two-valued cutoff centered at  <dig>  the middle of the excluded range of intensity measurements.
the cutoff was chosen to be 5500±δ, , to define the high  and low  affinity thresholds, thus leaving out of consideration those predicted peptides with measurements that fall in the “grey” zone of .
when δ= <dig>  we take into account all the measurements, when δ= <dig> we have a grey zone of peptides with signal intensities between  <dig> and  <dig>  figure  <dig> shows the performance  as well as the fraction of excluded peptides for both methods as well as the aggregate prediction as a function of δ. from fig.
<dig> top-left panel, we can see that the precision of the method of barbarini et al.
<cit>  at δ= <dig> is around  <dig> , much lower than the precision of close to  <dig>  of pythia-design.
roc curves showing the same conclusion are presented in fig.
<dig>
fig.
<dig>
a precision, b recall and c accuracy of peptide design.
performance of pythia-design, the method of  <cit>  , and the aggregate formed by taking the positives from pythia-design and the negatives from  <cit>  is shown.
in , the lines for pythia and the aggregate overlap.
the x-axis δ is a measure of the activities from the high- and low-binding affinities , and  shows the fraction of peptides excluded for a given δ
fig.
<dig> performance of the pythia-design and barbarini et al.
<cit>  method  for designing peptides with desired reactivities.
roc curves were determined from predicted peptides incubated with ivig  diluted to 1: <dig>  1: <dig> and 1: <dig>  and epitope-antibody reactivities  determined as described by lustrek et al.
2013
diversity of designed peptides
it is important that a design method is able to generate many diverse peptides rather than a small number of very similar peptides.
this allows for a greater sampling of the space of binding peptides and allows for more effective screening based on downstream criteria.
the set of peptides produced using the pythia-design approach is much more diverse than the set generated by the barbarini et al.
<cit>  approach .
to quantify the diversity, we create a graph g from the set of predicted high and low binding affinity peptides, where each peptide is a vertex in g and two different peptides are connected by an edge if they have fewer than a specified number c of differences under the hamming distance .
for each vertex v, we compute a maximal independent set that contains v. an independent set is a subset of vertices such that no two are connected by an edge.
this yields n maximal independent sets, where n is the order of the graph.
the size of the independent set containing v is a measure of how dissimilar v is to the other designed peptides.
we compute the average size of these maximal independent sets, and observe how this value changes as we vary the cutoff parameter c defining the edges of the graph.
at a given cutoff level c, if the average size of the maximal independent sets is larger, then there are more independent peptides and these peptides have, by construction, c or more differences.
pythia-design peptides exhibit substantially more diversity  than the designs of barbarini et al.
<cit> .
until a distance cutoff of  <dig>  almost all of the pythia-design peptides  belong to a single independent set that spans the entire graph, meaning that nearly all the peptides are dissimilar in at least  <dig> of their positions.
the designed peptides of barbarini et al.
<cit> , however, share a great deal of sequence similarity, representing a very dense sampling of the sequence space near only a few particular points.
fig.
<dig> diversity of the designed, predicted high- and low-binding affinity peptides.
the sequence diversity among the pythia-design peptides is significantly higher than the approach of  <cit> .
the y-axis gives a measure of diversity of a set of designed peptides  under a particular hamming-distance threshold defining similar peptides .
almost all the pythia-designed peptides differ in at least  <dig> of their  <dig> possible positions
aggregation of two strategies for constructing binding and non-binding peptides
we also created an aggregate approach for de novo peptide design to assess how pythia-design complements previous methods.
pythia-design and the method of barbarini et al.
<cit>  take different approaches to the design problem.
pythia-design generates dissimilar random peptides and classifies them, while the barbarini et al.
<cit>  method extracts some motifs from clusters of positive and negative peptides and then generates new peptide sequences.
the method of barbarini et al.
<cit>  produces a high number of peptides that, even though they were predicted positives, are actually non-reactive.
on the other hand pythia has very few non-reactive peptides in its positive set.
a confusion matrix  comparing the binding by ivig dilutions of peptides from pythia-design and the method of barbarini et al.
<cit> , showing dilutions of 1: <dig>  1: <dig>  1: <dig>  underline that the pythia-design method does preferentially select peptides that show high affinity binding to ivig.
six hundred peptides representing  <dig> binders and  <dig> non-binders of each of pythia-design and the method of barbarini et al.
were incubated with ivig.
the confusion matrix below indicates that the peptides selected by pythia-design bind antibodies with higher affinity than the peptides designed by barbarini et al.
to combine the strengths of each method, we used the strategy of barbarini et al.
<cit>  for generating low-affinity peptides and pythia-design for generating high-affinity predictions.
the accuracy of the aggregate method, with an auroc of  <dig> – <dig> , is better than either of the two methods .
one caveat, however, is that this combined method produces a set of non-binders with lower diversity than that produced by pythia-design in isolation.
fig.
<dig> performance of a method that combines pythia-design and barbarini et al.
<cit>
inclusion of citrulline, cysteine, and tryptophan in designed peptides
pythia-design was allowed to include the citrulline amino acid  in its designed peptides.
many of the designed peptides included this non-standard amino acid .
citrulline was included with high prevalence in both the tested high and low binders .
since both high and low binders included them, we do not see a significant effect on binding in general with the inclusion of these amino acids.
it is likely they were included with such high prevalence because such peptides are able to satisfied the imposed peptide diversity constraints.
similarly, peptides with cysteine  are over-represented in the pythia-design peptides , likely for partially similar reasons.
this interpretation is further supported by the fact that c and z were under-represented in the training and test sets used to train the pythia discriminative classifier on which pythia-design is based .
it is also consistent with the implementation of the classifier, which omitted features for which no training data  were available.
in particular, all of the propensity to include or exclude z in pythia-design comes from the string kernel and structural features extracted from the training set.
figure  <dig> shows that the true positive designed peptides often include tryptophan  and cysteine , which are relatively uniformly over-represented along the entire length of the designed peptides.
fig.
<dig> position-specific peptide propensity within true positive pythia-design peptides  divided by the pwm of the negative set of peptides.
pwm segments in red indicate amino acids that are predicted to interfere with antibody binding.
green highlights amino acids that favor binding of antibodies present in ivig.
an over-representation of cysteine  and tryptophan  in all positions is seen
the prevalence of citrulline and cysteine is likely due to the fact that citrulline and cystiine were less represented in the peptides used in the training data sets, allowing these designed peptides to more easily satisfy the imposed diversity requirements
discussion
one caveat of epitope mapping experiments is the nature of how peptides are presented e.g.
in solution or fixed to a support and/or whether the peptides are fixed amino- or carboxyterminally.
the peptides that we analyzed were coupled via the aminoterminal end to glass slides  <cit> .
the nature of how peptides are presented in solution or fixed to a platform might influence the binding affinities obtained.
this experimental restriction may affect their binding properties, and if peptides were bound via the carboxyterminal end, their binding affinities may change.
hecker et al.
<cit>  have demonstrated that the epitopes found via the experimental protocol used here represent epitopes that are functional in other assay systems as well, so it is likely that the predictive methods and their results will be robust to the experimental assay.
in addition, naturally the experiments reported here are specific to the ivig sets selected.
performance on additional sets is a promising direction for future work.
another caveat of the experimental validation performed here is the non-random selection procedure used for choosing which peptides were experimentally validated.
because stratified minmax sampling was used, the designed peptides were chosen so that they represented both high- and low- predicted binders of two other computational methods.
we find post-selection that the peptides designed by pythia-design fell relatively uniformly across these categories, so the sampling represents a mostly unbiased sampling of the designed peptides .
the binders designed by barbarini et al.
<cit>  displayed more bias toward appearing in only a few sampling categories , and so it is possible that this has led to a bias under-representing that method’s overall performance.
another caveat is our interpretation of the use of citrulline in our high-binding designed peptides.
although it appears that citrulline is primarily included only to increase peptide diversity, it is possible that the predicted peptides are only reactive to the subset of antibodies within the ivig serum that bind to citrullinated peptides.
in this case, the predicted results are still of interest, since they represent a high-binding set of peptides, but to this more restrictive class of antibodies.
despite its high performance, there is room for improvement of the pythia method for affinity binding prediction that is at the core of pythia-design.
in particular, zdock  <cit> , which is used to compute the structural features, considers only rigid docking of the peptides to immunoglobulin.
the structural features may change significantly if we use a non-rigid docking procedure, where conformational changes in the paratope, epitope, or both are allowed.
it is actually quite surprising that, using only a single igg <dig> structure model and using only perfectly rigid docking, the structure-based classifier obtained such respectable performance.
this indicates that improving the computation of the structural features is a promising way to increase prediction accuracy.
the integration of additional structural information might guide and improve computational processes studying epitope binding.
not all immunogobulins necessarily use similar or comparable binding modes  <cit> .
as such, knowing the structural heterogeneity of immunoglobulins found in one human individual might lead to the description of different types of epitope-antibody recognition modes  <cit> .
another potential source of improvement is the inclusion of phylogenetic information as was done in lin et al.
<cit> .
an important direction for future work is the determination in greater detail, from the svm weights, which sequence features are particularly indicative of binding or non-binding.
for a machine learning method to be reasonable, the testing and training peptides should be drawn from the same distribution, as we have done here.
the peptides used to train and test the pythia classifier were randomly chosen  from a large set of bound and non-bound peptides, and randomly divided into testing and training.
given the very high level of performance on the random set, it is unlikely that performance is driven by any small number of sequence patterns.
further, and most importantly, the generalizability of the pythia classifier was tested in a particularly strong way: it was used to design novel, sequence-dissimilar peptides and the accuracy on that task is very good, and provides strong evidence that the developed models are not specific to any single sequence pattern.
CONCLUSIONS
we have provided a method, pythia-design, for the design of peptides with specific reactivity properties , showing experimentally that the designs accurately exhibit the desired affinities.
in addition to producing more accurate designs than previous approaches, pythia-design is able to sample the space of possible designed peptides more completely, creating many dissimilar designs rather than variations on a few similar peptides.
again, combination of two dissimilar methods shows improved performance.
understanding antibody binding patterns is crucial for understanding the immune response for many human diseases.
we show that diverse sets of peptides that exhibit the desired binding properties can be computationally designed.
this work moves us closer to understanding the interplay and interactions between human antibodies and the targets to which they bind.
additional files
additional file  <dig>
table s <dig>  binding training data.
training set of measured binding affinities.
this collection of peptides was used to train the pythia classifier.
additional file  <dig>
table s <dig>  classification test binders.
set of binders in the test set including computational analysis.
the peptides in column a represent the test data set, sorted by column b. the highest measured values  are given in column b, in column c  the mean of secondary antibody control, and in column d  the mean of all ivig measures.
in column e, the sum of all computational methods are summed up whose predictions were correct as outlined in the remaining columns, where the el-manzalawy, luštrek, pwm, pythia, barbarini et al.
<cit>  , and pwm <dig> methods are given.
additional file  <dig>
table s <dig>  classification test non-binders.
set of non-binders in the test set including computational analysis.
the peptides in column a represent the test data set, sorted by column b. the highest measured values  are given in column b, in column c  the mean of secondary antibody control, and in column d  the mean of all ivig measures.
in column e, the sum of all computational methods are summed up whose predictions were correct as outlined in the remaining columns, as in additional file 2: table s <dig>
additional file  <dig>
table s <dig>  designed peptides and their experimental validation.
data set of the designed peptides including computational analysis.
the peptides in column a were further categorized as outlined in material and methods, see column e and f. the columns b, c, d show the computational assessment based on methods used in  <cit> .
the experimental data are given in column j, k, and l using dilutions 1: <dig>  1: <dig> and 1: <dig> of ivig purchased from omrix.
column m gives the background values to the binding of the secondary antibody to the peptides.
columns g, h, and i outlines whether binding occurred  or not .
the first excel data sheet shows the potential binders, the second excel sheet the non-binders.
additional file  <dig>
figure s <dig>  distribution of amino acids within training and testing sets.
figure s <dig>  selection bias of validated peptides.
competing interests
hjt is ceo of indymed, which studies epitope-antibody recognition in ivig sera and from patients that suffer from autoimmune disease.
although improved general understanding of epitope-antibody recognition may benefit indymed, the study and data reported here are not directly related to the business interests of indymed.
authors’ contributions
rp, rjp, rsr, gs, ck conceived and designed the experiments.
rp, rjp, pl performed the experiments.
rp, rn, rjp, fs, bz, nb, at, rb, hjt, gs, and ck analyzed the data.
rn contributed analysis tools.
rp, rn, ml, hjt, gs, and ck contributed to the writing of the manuscript.
all authors read and approved the final manuscript.

 antibody-protein interactions play a major role in infectious diseases, autoimmune diseases, oncology, vaccination and therapeutic interventions. antibodies present in human blood interact with antigens with different affinities and in a sequence- and structure-specific manner. when studying protein-antibody interactions, two types of epitopes are to be distinguished: conformational and linear epitopes. in this study we focus on linear epitopes; see a recent review for a discussion of conformational epitopes. all potential linear epitopes of a protein can be represented by short peptides derived from the primary amino acid sequence. the binding site of an epitope covered by an antibody typically includes a minimal stretch of to amino acids. if peptides of amino acids in length are incubated with one specific antibody, that antibody will bind to its epitope independently of the physical position of the binding motif within the peptide. motifs running from position to position up to motifs running from position to position would be possible. this uncertainty results in difficulties for determining consensus binding sites as well as meaningful position weight matrices. individual amino acids within epitope binding sites may have different impact on antibody recognition not only due to the nature of amino acids involved in binding but also because of the specific position of the amino acid within the whole peptide sequence. here, we present a method, pythia-design, for designing novel peptides with a desired binding affinity. this method is built upon a successful, novel discriminative classifier called pythia that can accurately label a given peptide as either a high- or low-affinity binder. to test the quality of the designs that pythia-design produces, we experimentally constructed our designed peptides and tested their binding affinity. we show that pythia-design more accurately designs such peptides than barbarini et al.. we further show that pythia-design produces a more diverse set of designed peptides, which is important for generating a varied set for experimental construction. finally, we show that the two methods of pythia-design and barbarini et al. can be combined, exploiting the relative strengths of both, to achieve even higher accuracy in epitope design. while there is less prior work on epitope design, much previous work has focused on the task of predicting binding affinity of a given peptide to various target molecules, e.g. antibodies, to mhc class i and class ii complexes alone or in concert with t cell receptor binding. machine learning classifiers such as artificial neural networks, hidden markov models, and support vector machines and other approaches have been explored in tackling the problem of predicting human leukocyte antigen binding peptides. much work has also focused on the prediction of t-cell and b-cell binding peptides. zhao et al. explore various classifiers to predict peptide t-cell binding. using a 10-dimensional feature vector to represent each amino acid, they discover that svms provide the best classification performance in their task. huang and dai also explore the classification of peptide binding to t-cells using a support vector machine classifier. they present a novel peptide feature based on combining a 20-dimensional indicator vector with amino acid similarity information encoded by the blosum matrix. nanni and lumini introduced the mpps system that relies on an ensemble of support vector machines, trained on various physicochemical properties, to classify peptide binding to hiv-protease and t-cells. they use sequential floating forward selection to select a subset of features, and combine the individual classifier predictions using the max rule. more recently, nanni and lumini have explored the use of a novel peptide-encoding scheme that relies on the use of nonlinear dimensionality reduction to extract the information encoded across a large number of physicochemical properties. they demonstrate that this novel feature representation, when used in conjunction with a support vector machine classifier, exhibits state-of-the-art performance in predicting peptide t-cell binding. wang et al. also showed that combining multiple classifiers using a consensus approach improved the classification of mhc class ii peptide binding predictions. others have used motif mining for mhc i and mhc ii peptides. recently, a flexible t cell receptor docking algorithm achieved near-native predictions for % of the tcr/pmhc cases. zhang et al. use 3d features and a random forest classifier to predict b-cell epitopes. lin et al. provide a method for b-cell epitope prediction that exploits phylogenetic information. yao et al. introduce an svm approach called svmtrip for b-cell linear epitope prediction. yao et al. compare various methods to find conformational epitopes of b-cells. el-manzalawy et al. used support vector machines in combination with a subsequence kernel reaching a auc of and an accuracy of % to predict peptide/epitope-antibody-binding. our novel discriminating classifier upon which our design method is based uses an ensemble of support vector machines to classify design candidates. this classifier is broadly similar to that of nanni and lumini. however, we use probabilistic svms with platt’s extension, along with a different set of features. in addition, no other previous work deals with such a wide variety of paratopes — the regions of antibodies which recognize antigens — as is found in intravenous immunoglobulin fractions, as we do here. as our experiments below show, a diverse set of ivig-binding peptides can be computationally designed using pythia-design with high accuracy. in addition, our random-walk strategy for ensuring diversity is general and could be applied to any accurate discriminative classifier. these computational techniques, and the collection of designed and validated ivig binders and non-binders, will be useful both for gaining a more thorough understanding of ivig binding properties and the diversity of possible epitopes. methods overview the pythia-design method has two main parts. the first part is a machine learning classifier that is trained to predict whether a given peptide is a high-affinity or low-affinity binder. we designed this classifier, which we refer to as pythia, to use an ensemble of probabilistic support vector machines trained on various sequence, chemical, and structural features. the features used in the classifiers are described in. for features where we did not compute the kernel directly, we used a radial basis function kernel. the second part of pythia-design is a method for using random walks to generate candidate peptides with novel sequences to feed into the classifier. this is described in section “generating novel peptides that bind to ivig”. discriminative classifier for predicting binding and non-binding epitopes we here describe a method for predicting whether a given peptide is going to be a low-affinity or high-affinity binder. motivated by the success of previous work in various protein-related prediction tasks, we use support vector machines as our classifiers within the ensemble. we trained these svms via the libsvm software. for all features, the optimal svm parameters were discovered via a grid search and 5-fold cross-validation. the parameters were selected based entirely on a training set without the availability of the subsequent testing set, which was held hidden from method developers until after predictions were made. for a given svm model, the cross-validation accuracy for the optimal set of parameters is used as a weight to combine the corresponding model’s predictions with the others from the ensemble. during standard svm classification, instances are assigned a hard label, as belonging to the negative or positive class, denoted by c− and c+ respectively. such a hard labeling poses no problem when only a single classifier is used to label test data. however, when an ensemble of classifiers is used, it is useful to have extra information about the degree to which the label assigned by each individual classifier should be trusted. for this reason, we chose to use platt’s extension, which provides probabilistic outputs for a support vector machine’s classifications. instead of receiving a 0/ label, each instance is given an a posteriori estimate of the probability with which it belongs to the positive class. thus, we expect that instances that clearly belong to the negative class will be given a value close to while instances that belong to the positive class will be given values close to having a probabilistic interpretation of the classification for data instances makes it possible to combine the output of different classifiers. we used a variant of the sum rule, where the predictions of the individual classifiers are summed and normalized to yield the prediction of the ensemble. specifically, the prediction of the ensemble for a particular instance xi was computed using p+ensxi=1a∑j=1majp+jxi, where xi is a feature vector representing the ith peptide, constructed using some subset of the features described in section “features used in the classifiers”, p+jxi is the a posteriori probability output by classifier j that the peptide with features xi is a high-affinity binder, and aj is classifier j’s cross-validation accuracy. a is a normalization factor equal to ∑j=0maj. we can then take p+ens to be the probability with which the ensemble predicts xi to belong to the positive class, or we can use it to obtain a discrete class prediction with the decision rule: xi∈c+ifp+ensxi≥τ,c−otherwise. in our experiments, we set τ=, but other values may be reasonable. in fact, one may even learn the value of τ which yields the best performance by using a held-out subset of the training data, though we do not explore that here. each svm model will yield a prediction for each peptide in the testing set. we combined the predictions for all of the classifiers in the ensemble using a variation on the approach presented by nanni and lumini, which is itself an extension of the sum-rule. we normalized the predictions for each classifier to have a standard deviation of next, we combined the predictions from each of the j classifiers according to eq. by sorting the peptides in the testing set according to this value, we can produce a rank ordered list of the peptides in order of the likelihood that they belong to the positive class. features used in the classifiers numerically encoded sequence features there are two distinct types of sequence features that we encode numerically. first, we used a simple variation on the peptide encoding scheme presented by huang and dai. we encoded each amino acid in the peptide by replacing its single letter code with its corresponding row in the blosum matrix. the blosum matrix contains empirically derived log-odds scores that encode the frequency of different amino acid substitutions and is commonly used to measure the similarity between different amino acids. let the peptide of length d be given as p=, where ai is the amino acid in the ith position of the peptide. further, let row map the amino acid a to its corresponding row in the blosum matrix. we encoded the peptide as enc=. for the length d peptide p, enc will be a 20d dimensional feature vector. in addition to blosum we use the same type of encoding with matrices nlf and sa introduced by nanni and lumini. these matrices are derived by performing dimensionality reduction on a large, rectangular matrix where each row corresponds to an amino acid and each column to some physicochemical property. the goal of the dimensionality reduction is to decorrelate the physicochemical properties, reducing the column space of the matrix significantly. the nlf matrix is a 20× matrix obtained using a nonlinear fisher transform, while the sa matrix is a 20× matrix obtained using a combination of clustering and principal component analysis. the second type of sequence feature that we encode numerically involves various physicochemical properties of the constituent amino acids of each peptide. we analyze the amino acid properties present in the amino acid index. each aaindex property provides a mapping from each of the amino acids to a numerical scale measuring some physicochemical attribute. the aaindex listed different amino acid properties. we use an approach based on a sliding window and histograms to turn each aaindex property into a numerical feature vector for a peptide. consider a single aaindex property aaij, and let aaij represent the numerical value to which the amino acid a is mapped under aaindex property j. to form a representation for the entire peptide p under the property aaij, we used a window of length w sliding across the peptide to produce a -dimensional vector where entry i in this vector is the average value of the aaindex property over the window starting at position i. by varying w, we can change the coarseness of this representation. through a process of experimenting with different values of w for this classification task, we computed these features for w∈. string kernel features string kernels are used to evaluate the sequence similarity between peptides. there are many different varieties of string kernels, ranging from the somewhat simple k-spectrum kernel, which essentially counts the occurrence of all length k substrings in each peptide, to the more complex substring-mismatch kernel, which considers all shared subsequences between two peptides, allowing for gaps and mismatches.we use the k-spectrum string kernel for k=,, the sssk kernel with parameter d= and the bounded range substring kernel with parameter r= the output of each of these methods is a matrix, known as the kernel matrix, in which the entry at row i and column j is the result of the kernel evaluation between peptides i and j. to train a svm model for each of these string kernels, we simply compute the kernel matrix, and then make use of the ability of libsvm to train a model using a precomputed kernel. structure features we use rigid docking to estimate the binding affinity between ivig and a candidate peptide. we computed a hypothesized 3-dimensional structure for each peptide using the biochemical algorithms library. we built a starting model for each peptide p by positioning the side chains for each amino acid by choosing the most frequently occurring rotamer position from a rotamer library. we then optimized this initial structure by performing an energy minimization using the amber force field. this relaxes the structure until a energy minimum is achieved. we also obtained an experimentally measured 3d structure for igg, the most prevalent class of igg antibody present in intravenous immunoglobulin. we measure the conformational complementarity of each of our hypothesized peptide structures with the immunoglobulin structure. to compute this complementarity, we performed a protein-protein docking simulation for each of the constructed peptides against igg using the zdock software. each zdock run produces a list of the top-ranked docking predictions for each peptide. the zdock score provides a measure of the complementarity of the peptide and immunoglobulin conformation in the docking region and is used as a proxy for the overall quality of the docking. for each peptide, we formed a histogram from the zdock scores, and use this histogram as a feature vector with which to train the svm model. intuitively, we expect peptides whose zdock score distributions are skewed toward high scores to have better shape complementarity and, therefore, to be more likely binders than peptides whose zdock score distributions are skewed toward low scores. generating novel peptides that bind to ivig pythia-design builds upon the classifier described above by first generating many peptide sequences and then assigning them a reactivity category high, medium, or low according to the predictions of our classifier. we generated the de novo sequences using a sampling approach that corresponds to a seeded random walk in sequence space. to obtain a sequence for reactivity class c, we choose a random seed sequence s∈c from the training set, and randomly mutate its constituent amino acids until it adheres to several required sequence diversity rules, in order to ensure that the novel peptide sequences that are generated are sufficiently different from those in the training set. specifically, the designed sequences predicted to be in the high- or low-reactivity category could not share any 4-mer, or exhibit a sequence identity of greater than amino acids in any subsequence of length with any sequence in the same reactivity category in the training set. further, the peptides generated by the pythia-design method were required to adhere to the same set of constraints with respect to the peptides in the original experiment’s testing set. we seeded our peptide generator with sequences from the training set — sequences with the highest experimentally measured reactivity and sequences with the lowest experimentally measured reactivity. running the random walks then produced candidate peptides that we classified using pythia and sorted according to their probability of belonging to the positive class. there were peptides with an probability greater than or equal to, and with an probability less than. the sequences with the highest probabilities were predicted belong to the high-affinity class h, while the sequences with the lowest probabilities were predicted to belong to the low-affinity class l. the remaining peptides, with probabilities closest to, were predicted to belong to the medium-affinity class m. training and testing sets of peptides for classifier to train and test our discriminative classifier, we used the data set of luštrek et al., in which peptides were incubated with commercially available intravenous immunoglobulin fractions, which was originally presented as part of the dream challenge. ivig is a mixture of naturally occurring human antibodies isolated from up to healthy individuals. from this dataset, high-confidence negative and positive pools of recognized peptides were determined based on epitope-antibody-reactivities of more than different peptides subjected to a huge number of structurally different antibodies present in ivig. the training and test datasets for the discriminative classifier were assembled from these peptide pools. from the collection of all the peptides incubated with human ivig, a pool of epitope-containing peptide sequences reactive with human immunoglobulins was experimentally identified. this was called the positive set. from the same original collection of peptides peptides were identified that showed no antibody binding activity in any of the triplicate assays. this peptide set was called the negative set. the training set was formed by random sampling of peptides from the positive set and peptides from the negative set. the training set thus created contained peptides and their respective binding reactivities. the test set was created by joining together the remaining peptides from the positive set and the remaining peptides from the negative set, for a total of peptides. selection of designed peptides for experimental validation pythia-design was used to generate peptides predicted to be reactive as well as predicted not to be reactive. the method of barbarini et al. was also used to generate peptides of each class. to select a subset to experimentally construct and validate, the designed peptides were re-categorized by selecting high binders and non-binders as follows. a stratified instead of a randomized sampling procedure was chosen in particular to investigate whether the designed peptides are robust to existing classification methods as well as position weight matrix analysis. the initial peptide set were subgrouped by using pwm and svm analyses as described in; see categories in legend of additional file 4: table s. categories were determined from a linear scale representing pwm and lustrek classification scores of the original data sets from minmax sampling of the classifiers was used. the stratified quota sampling was restricted to peptide sequences in case of binding peptides. in total, peptides predicted to be bound and predicted not being bound by ivig were taken from each of pythia-design and barbarini et al. and subjected to experimental testing.  overview the test set of peptides was withheld from the algorithm designers until after the algorithm was finalized. only the training set was used in the design and initial evaluation of the pythia classifier. once finalized, pythia was then evaluated on the held-out test set of peptides to validate the classifier component of pythia-design. finally, pythia-design was used to generate a collection of new likely binder and non-binder peptides, a subset of which were then experimentally validated. validation of the pythia affinity classifier although our goal is to produce a peptide design method, we first validate that our peptide classifier is accurate. this classifier is based on an ensemble of learners that aggregates the prediction of many individual classifiers, each of which was trained on a set of features. analyzing the auc for the individual classifiers in the pythia ensemble on the test peptides, we observe that many of the classifiers show similar performance, with the exception of the structural classifier which displays significantly lower classification performance. the ensemble, however, yields superior performance compared to any of its constituent classifiers with an auroc of. at any given false positive rate, the ensemble classifier will obtain a higher true positive rate than any of the other classifiers. the precision-recall curve presents a related view of classifier performance to the roc curve. it measures how the precision changes as the recall is increased. for very small recall values, the sparse spatial sample and k-spectrum string kernels yield the best precision. however, for the vast majority of recall values, the ensemble classifier yields the highest precision. just as was the case with the roc curves, the ensemble again achieved the maximum area under the pr curve. while the auprs were generally lower than the aurocs, we did observe that the benefit of the ensemble was larger with respect to the pr curves than the roc curves. local compositiona athe best single classifier under both the auroc and aupr metrics boldface indicates the best solution because the test and training sets of peptides were chosen randomly, it is possible that overlapping or shared sequences between the test and training set partially leads to this high performance. nevertheless, our ultimate aim is to design novel peptides, and this analysis suggests that pythia has reasonable performance identifying low- and high-affinity binding peptides. accuracy of predicted affinity of computationally designed peptides pythia-design was used to generate a number of probable high-affinity binding peptides, probable low-affinity binding peptides, and medium-affinity binding peptides. in addition, the method of barbarini et al. was used to generate the same number of peptides. a subset of designed peptides of both the pythia-design and barbarini et al. methods were experimentally constructed. the binding affinities of these designed peptides were then experimentally measured. additional file 4: table s gives the designed peptides and their measured affinities. fig. quality of designed peptides from two approaches. the distribution of measured affinities for the designed peptides predicted to belong to the low and high binding affinity classes for the pythia-design method, and method of barbarini et al.. the horizontal line at indicates the binding affinity cutoff above which a peptide is considered to have a high binding affinity. both methods produce a statistically significant separation of high- and low- binders, but pythia-design is much better at generating high-affinity binders to quantify the degree to which the methods are able to generate binders and non-binders effectively, we use the method of ojala and garriga to compute a p-value indicating the probability of the observed high- and low-binding separation. to do this, we compute the standard f1-score of the predictions, then randomly permute the labels times to get a distribution of f1-scores. for both pythia-design and pavia, the true f1-score is always better than the randomized score, meaning that for each method p<. this indicates that the approaches are both truly designing peptides better than random. using specific binding defined by a factor above control measures plus a minimum signal of above the designed peptides were categorized as truly or falsely predicted as outlined in additional file 4: table s using these categorizations, we computed precision, recall, and accuracy for the two methods. since we do not have a hard threshold to define high or low affinity measurements, and such threshold cannot be detected from the distribution of measurements itself, we used a two-valued cutoff centered at the middle of the excluded range of intensity measurements. the cutoff was chosen to be 5500±δ,, to define the high and low affinity thresholds, thus leaving out of consideration those predicted peptides with measurements that fall in the “grey” zone of. when δ= we take into account all the measurements, when δ= we have a grey zone of peptides with signal intensities between and figure shows the performance as well as the fraction of excluded peptides for both methods as well as the aggregate prediction as a function of δ. from fig. top-left panel, we can see that the precision of the method of barbarini et al. at δ= is around, much lower than the precision of close to of pythia-design. roc curves showing the same conclusion are presented in fig. fig. a precision, b recall and c accuracy of peptide design. performance of pythia-design, the method of, and the aggregate formed by taking the positives from pythia-design and the negatives from is shown. in, the lines for pythia and the aggregate overlap. the x-axis δ is a measure of the activities from the high- and low-binding affinities, and shows the fraction of peptides excluded for a given δ fig. performance of the pythia-design and barbarini et al. method for designing peptides with desired reactivities. roc curves were determined from predicted peptides incubated with ivig diluted to 1: 1: and 1: and epitope-antibody reactivities determined as described by lustrek et al. 2013 diversity of designed peptides it is important that a design method is able to generate many diverse peptides rather than a small number of very similar peptides. this allows for a greater sampling of the space of binding peptides and allows for more effective screening based on downstream criteria. the set of peptides produced using the pythia-design approach is much more diverse than the set generated by the barbarini et al. approach. to quantify the diversity, we create a graph g from the set of predicted high and low binding affinity peptides, where each peptide is a vertex in g and two different peptides are connected by an edge if they have fewer than a specified number c of differences under the hamming distance. for each vertex v, we compute a maximal independent set that contains v. an independent set is a subset of vertices such that no two are connected by an edge. this yields n maximal independent sets, where n is the order of the graph. the size of the independent set containing v is a measure of how dissimilar v is to the other designed peptides. we compute the average size of these maximal independent sets, and observe how this value changes as we vary the cutoff parameter c defining the edges of the graph. at a given cutoff level c, if the average size of the maximal independent sets is larger, then there are more independent peptides and these peptides have, by construction, c or more differences. pythia-design peptides exhibit substantially more diversity than the designs of barbarini et al.. until a distance cutoff of almost all of the pythia-design peptides belong to a single independent set that spans the entire graph, meaning that nearly all the peptides are dissimilar in at least of their positions. the designed peptides of barbarini et al., however, share a great deal of sequence similarity, representing a very dense sampling of the sequence space near only a few particular points. fig. diversity of the designed, predicted high- and low-binding affinity peptides. the sequence diversity among the pythia-design peptides is significantly higher than the approach of. the y-axis gives a measure of diversity of a set of designed peptides under a particular hamming-distance threshold defining similar peptides. almost all the pythia-designed peptides differ in at least of their possible positions aggregation of two strategies for constructing binding and non-binding peptides we also created an aggregate approach for de novo peptide design to assess how pythia-design complements previous methods. pythia-design and the method of barbarini et al. take different approaches to the design problem. pythia-design generates dissimilar random peptides and classifies them, while the barbarini et al. method extracts some motifs from clusters of positive and negative peptides and then generates new peptide sequences. the method of barbarini et al. produces a high number of peptides that, even though they were predicted positives, are actually non-reactive. on the other hand pythia has very few non-reactive peptides in its positive set. a confusion matrix comparing the binding by ivig dilutions of peptides from pythia-design and the method of barbarini et al., showing dilutions of 1: 1: 1: underline that the pythia-design method does preferentially select peptides that show high affinity binding to ivig. six hundred peptides representing binders and non-binders of each of pythia-design and the method of barbarini et al. were incubated with ivig. the confusion matrix below indicates that the peptides selected by pythia-design bind antibodies with higher affinity than the peptides designed by barbarini et al. to combine the strengths of each method, we used the strategy of barbarini et al. for generating low-affinity peptides and pythia-design for generating high-affinity predictions. the accuracy of the aggregate method, with an auroc of –, is better than either of the two methods. one caveat, however, is that this combined method produces a set of non-binders with lower diversity than that produced by pythia-design in isolation. fig. performance of a method that combines pythia-design and barbarini et al. inclusion of citrulline, cysteine, and tryptophan in designed peptides pythia-design was allowed to include the citrulline amino acid in its designed peptides. many of the designed peptides included this non-standard amino acid. citrulline was included with high prevalence in both the tested high and low binders. since both high and low binders included them, we do not see a significant effect on binding in general with the inclusion of these amino acids. it is likely they were included with such high prevalence because such peptides are able to satisfied the imposed peptide diversity constraints. similarly, peptides with cysteine are over-represented in the pythia-design peptides, likely for partially similar reasons. this interpretation is further supported by the fact that c and z were under-represented in the training and test sets used to train the pythia discriminative classifier on which pythia-design is based. it is also consistent with the implementation of the classifier, which omitted features for which no training data were available. in particular, all of the propensity to include or exclude z in pythia-design comes from the string kernel and structural features extracted from the training set. figure shows that the true positive designed peptides often include tryptophan and cysteine, which are relatively uniformly over-represented along the entire length of the designed peptides. fig. position-specific peptide propensity within true positive pythia-design peptides divided by the pwm of the negative set of peptides. pwm segments in red indicate amino acids that are predicted to interfere with antibody binding. green highlights amino acids that favor binding of antibodies present in ivig. an over-representation of cysteine and tryptophan in all positions is seen the prevalence of citrulline and cysteine is likely due to the fact that citrulline and cystiine were less represented in the peptides used in the training data sets, allowing these designed peptides to more easily satisfy the imposed diversity requirements discussion one caveat of epitope mapping experiments is the nature of how peptides are presented e.g. in solution or fixed to a support and/or whether the peptides are fixed amino- or carboxyterminally. the peptides that we analyzed were coupled via the aminoterminal end to glass slides. the nature of how peptides are presented in solution or fixed to a platform might influence the binding affinities obtained. this experimental restriction may affect their binding properties, and if peptides were bound via the carboxyterminal end, their binding affinities may change. hecker et al. have demonstrated that the epitopes found via the experimental protocol used here represent epitopes that are functional in other assay systems as well, so it is likely that the predictive methods and their results will be robust to the experimental assay. in addition, naturally the experiments reported here are specific to the ivig sets selected. performance on additional sets is a promising direction for future work. another caveat of the experimental validation performed here is the non-random selection procedure used for choosing which peptides were experimentally validated. because stratified minmax sampling was used, the designed peptides were chosen so that they represented both high- and low- predicted binders of two other computational methods. we find post-selection that the peptides designed by pythia-design fell relatively uniformly across these categories, so the sampling represents a mostly unbiased sampling of the designed peptides. the binders designed by barbarini et al. displayed more bias toward appearing in only a few sampling categories, and so it is possible that this has led to a bias under-representing that method’s overall performance. another caveat is our interpretation of the use of citrulline in our high-binding designed peptides. although it appears that citrulline is primarily included only to increase peptide diversity, it is possible that the predicted peptides are only reactive to the subset of antibodies within the ivig serum that bind to citrullinated peptides. in this case, the predicted results are still of interest, since they represent a high-binding set of peptides, but to this more restrictive class of antibodies. despite its high performance, there is room for improvement of the pythia method for affinity binding prediction that is at the core of pythia-design. in particular, zdock, which is used to compute the structural features, considers only rigid docking of the peptides to immunoglobulin. the structural features may change significantly if we use a non-rigid docking procedure, where conformational changes in the paratope, epitope, or both are allowed. it is actually quite surprising that, using only a single igg structure model and using only perfectly rigid docking, the structure-based classifier obtained such respectable performance. this indicates that improving the computation of the structural features is a promising way to increase prediction accuracy. the integration of additional structural information might guide and improve computational processes studying epitope binding. not all immunogobulins necessarily use similar or comparable binding modes. as such, knowing the structural heterogeneity of immunoglobulins found in one human individual might lead to the description of different types of epitope-antibody recognition modes. another potential source of improvement is the inclusion of phylogenetic information as was done in lin et al.. an important direction for future work is the determination in greater detail, from the svm weights, which sequence features are particularly indicative of binding or non-binding. for a machine learning method to be reasonable, the testing and training peptides should be drawn from the same distribution, as we have done here. the peptides used to train and test the pythia classifier were randomly chosen from a large set of bound and non-bound peptides, and randomly divided into testing and training. given the very high level of performance on the random set, it is unlikely that performance is driven by any small number of sequence patterns. further, and most importantly, the generalizability of the pythia classifier was tested in a particularly strong way: it was used to design novel, sequence-dissimilar peptides and the accuracy on that task is very good, and provides strong evidence that the developed models are not specific to any single sequence pattern.  we have provided a method, pythia-design, for the design of peptides with specific reactivity properties, showing experimentally that the designs accurately exhibit the desired affinities. in addition to producing more accurate designs than previous approaches, pythia-design is able to sample the space of possible designed peptides more completely, creating many dissimilar designs rather than variations on a few similar peptides. again, combination of two dissimilar methods shows improved performance. understanding antibody binding patterns is crucial for understanding the immune response for many human diseases. we show that diverse sets of peptides that exhibit the desired binding properties can be computationally designed. this work moves us closer to understanding the interplay and interactions between human antibodies and the targets to which they bind. additional files additional file table s binding training data. training set of measured binding affinities. this collection of peptides was used to train the pythia classifier. additional file table s classification test binders. set of binders in the test set including computational analysis. the peptides in column a represent the test data set, sorted by column b. the highest measured values are given in column b, in column c the mean of secondary antibody control, and in column d the mean of all ivig measures. in column e, the sum of all computational methods are summed up whose predictions were correct as outlined in the remaining columns, where the el-manzalawy, luštrek, pwm, pythia, barbarini et al., and pwm methods are given. additional file table s classification test non-binders. set of non-binders in the test set including computational analysis. the peptides in column a represent the test data set, sorted by column b. the highest measured values are given in column b, in column c the mean of secondary antibody control, and in column d the mean of all ivig measures. in column e, the sum of all computational methods are summed up whose predictions were correct as outlined in the remaining columns, as in additional file 2: table s additional file table s designed peptides and their experimental validation. data set of the designed peptides including computational analysis. the peptides in column a were further categorized as outlined in material and methods, see column e and f. the columns b, c, d show the computational assessment based on methods used in. the experimental data are given in column j, k, and l using dilutions 1: 1: and 1: of ivig purchased from omrix. column m gives the background values to the binding of the secondary antibody to the peptides. columns g, h, and i outlines whether binding occurred or not. the first excel data sheet shows the potential binders, the second excel sheet the non-binders. additional file figure s distribution of amino acids within training and testing sets. figure s selection bias of validated peptides. competing interests hjt is ceo of indymed, which studies epitope-antibody recognition in ivig sera and from patients that suffer from autoimmune disease. although improved general understanding of epitope-antibody recognition may benefit indymed, the study and data reported here are not directly related to the business interests of indymed. authors’ contributions rp, rjp, rsr, gs, ck conceived and designed the experiments. rp, rjp, pl performed the experiments. rp, rn, rjp, fs, bz, nb, at, rb, hjt, gs, and ck analyzed the data. rn contributed analysis tools. rp, rn, ml, hjt, gs, and ck contributed to the writing of the manuscript. all authors read and approved the final manuscript. 
analysis of the experimental validation of pythia-design peptides indicates that binding of ivig is favored by epitopes that contain trypthophan and cysteine. understanding the interactions between antibodies and the linear epitopes that they recognize is an important task in the study of immunological diseases. we present a novel computational method for the design of linear epitopes of specified binding affinity to intravenous immunoglobulin.
 understanding the interactions between antibodies and the linear epitopes that they recognize is an important task in the study of immunological diseases. we present a novel computational method for the design of linear epitopes of specified binding affinity to intravenous immunoglobulin.  we show that the method, called pythia-design can accurately design peptides with both high-binding affinity and low binding affinity to ivig. to show this, we experimentally constructed and tested the computationally constructed designs. we further show experimentally that these designed peptides are more accurate that those produced by a recent method for the same task. pythia-design is based on combining random walks with an ensemble of probabilistic support vector machines classifiers, and we show that it produces a diverse set of designed peptides, an important property to develop robust sets of candidates for construction. we show that by combining pythia-design and the method of: 2011), we are able to produce an even more accurate collection of designed peptides. analysis of the experimental validation of pythia-design peptides indicates that binding of ivig is favored by epitopes that contain trypthophan and cysteine.  our method, pythia-design, is able to generate a diverse set of binding and non-binding peptides, and its designs have been experimentally shown to be accurate. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords protein bindingmachine learningantibodiesprotein designnational science foundation1053918kingsford carl national science foundation0849899kingsford carl national science foundation0812111kingsford carl national science foundation1256087kingsford carl national institutes of healthai085376kingsford carl national institutes of healthhg006913kingsford carl state of mecklenburg-vorpommern mv-excellenceur08051thiesen hans-jürgen state of mecklenburg-vorpommern mv-excellenceur09012thiesen hans-jürgen bmbf0313692athiesen hans-jürgen bmbf0313692bthiesen hans-jürgen bmbf0315450dthiesen hans-jürgen bmbf0315450gthiesen hans-jürgen innovative medicine initiativeimi/115006barbarini nicola issue-copyright-statement© the author 2016 