# Hierarchy and extremes in selections from pools of randomized proteins

^{a}Laboratoire Interdisciplinaire de Physique, CNRS and Université Grenoble Alpes, 38000 Grenoble, France;^{b}Laboratoire de Biochimie, Chimie-Biologie-Innovation UMR8231, CNRS and Ecole Supérieure de Physique et Chimie Industrielles ParisTech, Paris Sciences & Lettres Research University, 75005 Paris, France

See allHide authors and affiliations

Edited by Boris I. Shraiman, University of California, Santa Barbara, CA, and approved February 5, 2016 (received for review September 8, 2015)

## Significance

Evolution by natural selection requires populations to be sufficiently diverse, but merely counting the number of different individuals provides a poor indication of the potential of a population to satisfy a new selective constraint. To achieve a more relevant characterization of this selective potential, we performed in vitro experiments of selection with populations of partially randomized proteins and analyzed the results quantitatively by high-throughput sequencing. We find that selective potentials in these populations follow simple statistical laws, which can be interpreted with extreme value theory (the mathematical theory of extreme events—here, the rare finding of a protein meeting the selective constraints). Our results provide an approach to quantitatively measure the selective potential of a population.

## Abstract

Variation and selection are the core principles of Darwinian evolution, but quantitatively relating the diversity of a population to its capacity to respond to selection is challenging. Here, we examine this problem at a molecular level in the context of populations of partially randomized proteins selected for binding to well-defined targets. We built several minimal protein libraries, screened them in vitro by phage display, and analyzed their response to selection by high-throughput sequencing. A statistical analysis of the results reveals two main findings. First, libraries with the same sequence diversity but built around different “frameworks” typically have vastly different responses; second, the distribution of responses of the best binders in a library follows a simple scaling law. We show how an elementary probabilistic model based on extreme value theory rationalizes the latter finding. Our results have implications for designing synthetic protein libraries, estimating the density of functional biomolecules in sequence space, characterizing diversity in natural populations, and experimentally investigating evolvability (i.e., the potential for future evolution).

Diversity is the fuel of evolution by natural selection, but translating this concept into quantitative measurements is not straightforward (1). A simple count of the number of different individuals in a population, for instance, fails to account for the very different responses to selection that two populations with the same number of different individuals may elicit. The problem is even acute at the molecular scale, where it also takes a very practical form: libraries of diverse proteins are routinely screened as a way to identify biomolecules of interest (binders, catalysts, etc.), and a proper “diversity” is critical for success (2, 3). However, beyond a general agreement that maximizing the number of different elements is desirable, there is no general rule for engineering and comparing diversity in these libraries.

A common design of many protein libraries is to concentrate variations at one or a few variable parts located around a fixed “framework,” which is shared by all members of the library (2, 3). The natural design of antibody repertoires, the pools of immune proteins with potential to recognize nearly every molecular target, follows this pattern. Most of the sequence variations in antibodies are, indeed, concentrated at a few loops extending from a common structural scaffold (4). This architecture has inspired the conception of artificial protein libraries built on frameworks other than the Ig fold (5).

Here, we present an approach to quantitatively characterize the selective potential of molecular libraries. To develop this approach, we designed and screened 24 synthetic protein libraries with identical sequence variations but different frameworks and analyzed their response to well-defined selective pressures by high-throughput sequencing. Between libraries, we find that selective potentials vary widely and define a hierarchy of frameworks. Within libraries, we find that selective potentials exhibit a simple scaling law, characterized by few parameters. The essence of these results is captured by an elementary probabilistic model based on extreme value theory (EVT).

Previous work has quantified the functional potential of totally or partially random biomolecules by counting the number of positive hits resulting from successive rounds of selections and amplifications of a large sample of these biomolecules (6⇓⇓⇓⇓–11). Our results lead us to propose a different approach to characterize the selective potential of a population. Compared with previous analyses, this approach does not depend strongly on the sensitivity of the experimental assay or the number of copies in which each distinct biomolecule is present in the initial population.

## Experimental Approach

### Library Design.

We built 24 minimal libraries with different frameworks but identical sequence diversity (*Materials and Methods*, Fig. 1, *SI Appendix*, Fig. S1, and Dataset S1). Twenty frameworks consist of single-domain antibodies taken from natural heavy-chain genes of diverse origins (*SI Appendix*, Fig. S2); they originate from maturated antibodies, which are mutated relative to their germ-line form, except for the S1 framework that comes from a germ-line (naïve) antibody. Three additional frameworks are more closely related and correspond to the germ-line and two maturated forms of the same human antibody, with the maturated frameworks sharing 65% and 85% sequence identity with the germ line. Finally, one framework consists exclusively of glycines to serve as a control. Diversity is limited to four consecutive amino acids at the complementarity determining region 3 (CDR3), the part of antibody sequences most critical for specificity (12). Structurally, the CDR3 forms one of three loops that define the binding pocket of a

### Selection.

We screened our libraries by phage display for binding to one of two targets: a neutral synthetic polymer, polyvinylpyrrolidone (PVP), and a short DNA loop of 9 nt (*Materials and Methods*). Two previous studies established the capacity of antibody phage display to select binders for these targets (15, 16). Phage display is a standard high-throughput screening technique (17). It is based on the fusion of each antibody sequence to the sequence of the pIII surface protein of the filamentous bacteriophage M13, a natural virus of the bacterium *Escherichia coli* with the shape of a 1-μm-long and 10-nm-wide cylinder (17). The engineered phage encapsulates the DNA sequence of an antibody and displays the corresponding polypeptide at its surface. Populations of up to *SI Appendix*). Our analysis and its interpretation, however, do not rely on the precise nature of the selective pressure.

### High-Throughput Sequencing.

We sequenced samples of *i* in the population at each round *Materials and Methods*).

Provided that a sequence *i* is present in many copies *i* as*a* so that *t* of selection; we explain below how our conclusions depend on this choice. We compare the frequencies between rounds

Previous studies have applied next generation sequencing to the outcome of phage display screens as a way to identify a large number of binders (20, 21) or infer sequence–function relationships (19) but have not investigated the statistical properties of the distribution of the relative selectivities of these binders.

### Reproducibility and Specificity.

Several observations based on the frequencies and amino acid patterns of the sequences in populations under selection validate our experimental approach. (*i*) Screening the same library against the same target in separate experiments yields reproducible frequencies *SI Appendix*, Fig. S3). (*ii*) Screening the same library against different targets yields target-specific amino acid patterns (*SI Appendix*, Fig. S4). (*iii*) Screening two libraries against the same target yields library-specific amino acid patterns (*SI Appendix*, Fig. S4). Taken together, these results show that enrichment of some of the sequences is reproducible and that it arises from selection for specific binding to the targets.

We note that one feature of our experiments is critical for reproducibility: the initial populations have a large degeneracy (the number of copies of each sequence) and not just a large diversity (the number of distinct sequences). For a sequence *i* with probability

## Results

### Hierarchy Between Libraries.

To compare the selective potentials of libraries built around different frameworks, we performed experiments in which the initial population of sequences consists of a mixture of libraries with distinct frameworks—a metalibrary. The results of these experiments reveal a striking hierarchy. Diverse members of the same library (i.e., sequences sharing a common framework) typically dominate. When repeating the experiment with an initial mixture of libraries that excludes the dominating library, another library dominates (Fig. 2). Libraries not selected when mixed with other libraries, nevertheless, do contain sequences with detectable selectivities as shown by screening them in isolation (*SI Appendix*, Fig. S4). These results are not explained by uneven representations of the libraries in the initial population (because the distribution of frequencies at round 2 is remarkably different from the distribution at round 1) or framework-specific differences during amplification (*SI Appendix*, Fig. S5).

Differences in frameworks are, thus, generally more significant than differences between variable parts, although these parts are clearly under selection for binding (different CDR3s have different selectivities) (Fig. 3 *B* and *D*). This result may not be surprising for very dissimilar frameworks, but our frameworks are all expected to share the same structural fold, and some frameworks have few sequence differences. In particular, the dominating framework when selecting the mixture of all 24 libraries against the DNA target (Fig. 2) is a germ-line human *SI Appendix*, Fig. S6), the dominating framework is the only other germ-line framework of the mixture (the S1 framework). As noted previously, differences between frameworks also appear in the patterns of amino acids that are selected at the level of CDR3s (*SI Appendix*, Fig. S4 *C–E*).

### Scaling Within Libraries.

To compare the selectivities of sequences sharing a common framework and therefore, differing by, at most, four amino acids (Fig. 1), we rank these sequences in decreasing order of their selectivity *r*, then for the sequences with top ranks,*A* shows an example where the exponent is *B*), and for several experiments, a power law cannot be justified (Fig. 3*D*).

Both the power law and its various deviations can, however, be rationalized under an elementary mathematical model. This model rests on two assumptions. First, it assumes that the selectivity of each sequence in a library is drawn independently at random from a common probability density

The model is, thus, probabilistic, although—barring experimental noise—the experiments have no inherent stochastic element. To the extent that selectivity reflects binding at thermodynamic equilibrium, the selectivity *i* is, indeed, determined by its binding free energy *T* represents the temperature, and

EVT, indeed, indicates that random variables *s* independently drawn from the tail of a common probability density have themselves a probability density of the form (24)*a* introduced in Eq. **1**), and *a*), which defines the universality class to which the distribution of selectivities belongs: the probability densities *SI Appendix*, Fig. S18).

As suggested by the notations, when **2**. Mathematically, when considering a large number *N* of samples, the rank *s* as predicted by Eq. **4**, for **2**. In other words, the power law seen in Fig. 3*A* corresponds to the expected relationship between the rank and the values of random variables drawn from the tail of a probability density when this density belongs to a class associated with

To precisely assess the ability of our model to describe all of the different cases, we followed the point over threshold approach, a standard method in applications of EVT to empirical data (24). This approach consists of fitting the data *A*, an illustration is provided in Fig. 4*A*, with error bars indicating 95% confidence intervals (*SI Appendix* discusses the analyses of other experiments). In this example, we observe that

Given *s*) against the cumulative **5**, a straight line *B*, *Inset* shows to be nearly the case in this example. The Q-Q plot makes a similar comparison but by representing *s* against *s* above which a fraction *x* of the data is located. This representation has two advantages over the P-P plot: it relies only on the estimation of κ, and it displays more clearly the contribution of the most extreme values. A straight line is expected if the fit is perfect but this time, with a slope τ and a *y*-intercept *B* indicates, again, a very good fit in the illustrated case.

Performing the same analysis on results of selections of various libraries against various targets, we find that the model is able to describe all of the experiments (*SI Appendix*, Figs. S8–S12). Different values of κ are obtained with differences that are statistically significant (*SI Appendix*, Table S1). In particular, the three cases,

Although many models can lead to a power law (25), our probabilistic model has the merit of explaining the various deviations from this behavior that the data exhibit. First, when *SI Appendix*, Fig. S7), and sequences of smaller selectivities, where **4** can provide an excellent fit well beyond the point where the power law applies (*B* and *D* and *SI Appendix*, Fig. S19). Second, EVT predicts behaviors differing from a power law if the probability density *D* and *SI Appendix*, Fig. S10).

## Discussion

We presented a quantitative analysis of in vitro selections of multiple libraries of partially randomized proteins with variations limited to four consecutive amino acids. The distribution of selectivities of the top-ranked sequences is described by few parameters, with an interpretation provided by an elementary probabilistic model based on EVT.

Within a library with members that share a common framework, this distribution is characterized by a shape parameter κ, which may be positive, negative, or zero. This parameter is independent of the factor *a* in Eq. **1** and has several interpretations. For instance, it controls the relative spacing between selectivities: ranking the sequences from best to worst, the expected difference of selectivity between sequences at rank *r* and *SI Appendix*). The shape parameter also provides a statistical answer to the following question: if sampling *N* sequences yields a top-ranked sequence of selectivity *SI Appendix*, Fig. S13); as a consequence, multiplying by a factor of 1,000 the number of sequences when

Other than the shape parameter κ, the other parameters are the scaling parameter τ, the threshold of selectivities parameter *SI Appendix*). Within our experimental setup, where the selectivities are determined only up to a multiplicative factor (Eq. **1**), the values of *SI Appendix*, Table S1).

Based on these results, we propose these parameters as general descriptors of the selective potential of a population of random variants facing a given selective constraint. In particular, these descriptors could be applied to revisit the fundamental problem of estimating the density of functional proteins or RNA in sequence space. Previous studies have estimated this density by counting the number of different sequences enriched in in vitro selections (6, 7). The results of such experiments depend on experimental noise, which sets a lower limit

Power laws are seemingly ubiquitous in distributions of protein features (26, 27). Most closely related to our work, the distribution of abundances of distinct antibody sequences in zebrafish has been shown to follow a power law with exponent *n* times the same selection leads to

Although many models can be consistent with a power law, our model based on EVT covers without additional assumption the deviations from a power law observed in the data; in particular, it can fit the data over a wider range of selectivities and account for nonpower law behaviors. Our work is, however, not the first application of EVT to the description of biological variation: Gillespie (31, 32) first introduced it in models of evolutionary dynamics as a way to constrain the distribution of beneficial effects obtained when mutating a wild-type individual. Gillespie (31, 32) assumed

Several experimental studies have recently investigated the value of κ applicable to the distribution of beneficial effects in viral or bacterial populations (35, 36). The sample sizes available in these studies are, however, insufficient to conclusively validate or invalidate the EVT hypothesis. In these experiments, the number of mutants found in the tail has, indeed, been so far very low (of the order of a dozen); estimating the sign of the shape parameter κ can be attempted (37), but assessing the validity of the fit using Q-Q plots as in Fig. 4 is not possible with such limited data. Our rich dataset provides a thorough test of the applicability of EVT to the analysis of biological diversity.

Comparable datasets are now being increasingly produced. In particular, several groups have characterized the phenotype of every single-point mutant of a protein (38). Our model may be viewed as a mathematical formalization of the concept of a random library, from which single-point mutants may deviate. We note, however, that selectivities from nonrandom subsets of one of our libraries do follow the same model as the full library (*SI Appendix*, Fig. S14). In any case, significant deviations will have to be quantified against our null model.

Beyond protein libraries, the model is relevant to the screening of synthetic chemical libraries, including the combinatorial libraries of small molecules developed in the pharmaceutical industry for drug discovery (39, 40). In this context, one previous study was performed with enough data points to possibly discriminate between different universality classes but considered only the exponential case

Finally, our work raises a question for future studies: if the selective potential of a partially randomized library is captured by few parameters and if these parameters can vary from library to library, what controls them? More simply, what features of the framework define a universality class? For instance, how does extending the variable parts to other sites change κ? The patterns of amino acids forming the sequences, which we have analyzed here only to confirm the reproducibility of the experimental results and their specificity with respect to the targets and libraries, may provide valuable insights (29).

The question may also be asked at another level: can we or natural evolution control these parameters to optimize the selective potential of a population? This question relates to the debated “evolution of evolvability” (42, 43), cast here into a concrete conceptual and experimental setting. Antibodies potentially define an excellent model system to experimentally study this question, because they are subject to selection and maturation toward a diversity of targets as part of their natural function. The approach and concepts introduced in this work provide the means to address the problem with quantitative experiments.

## Materials and Methods

### Phage Display.

PVP plates were prepared as described in ref. 15. The DNA target was prepared by self-assembly of a hairpin DNA, labeled with biotin at its 5′ end (5′-biotin: AAAAGACCCCATAGCGGTCTGCGT), and purchased from Eurogentec. *E. coli* TG1-competent cells were purchased from Lucigen Lt. Phage production, phage display screens based on the pIT2 phagemid vector, and helper phage KO7 production were performed following the standard protocol from Source BioScience (lifesciences.sourcebioscience.com/media/143421/tomlinsonij.pdf) and our own previous work (15, 44), with some modifications as specified in *SI Appendix*.

### Sequencing Data.

Library phagemids were purified from *E. coli* stocks after each selection round using Midiprep Kits from Macherey-Nagel. v3 Illumina MiSeq sequencing was performed by Eurofins Genomics. The MiSeq paired-end technology was used. Frameworks were recovered on the forward read, and only the reads having all of the expected restriction sites and less than four errors on 126 base pairs were kept. The CDR3s were accessible on the reverse read, and only the reads having all of the expected restriction sites and an average value of quality read *Q* > 30 on 12 base pairs defining the CDR3 were kept (*SI Appendix*, Table S2 has an estimation of sequencing errors). Datasets (Datasets S2–S19) are provided, with the identity of the framework in the column 1, the CDR3 sequence in column 2, and the count in column 3.

### Computational Analysis.

We infer the selectivity *i* by Eq. **1** with *i* present in the sample. Given sampling errors, estimated as *SI Appendix*, Table S2), the estimation of *SI Appendix*, Table S3). With *i*.

### Extreme Value Statistics.

We followed the standard approach for modeling threshold excesses (24). The parameters κ and τ were estimated by maximum likelihood, and the 95% confidence intervals shown in Fig. 4*A* were obtained under the hypothesis of normality by calculating the inverse of Fisher’s information. To ensure that the data allow us to discriminate between *P* value was calculated by a likelihood ratio test, whose distribution was estimated by numerical simulations. Maximum likelihood estimations are calculated on at least 50 data points. Codes in the format of an ipython notebook are provided in *SI Appendix* to facilitate similar analyses with other datasets.

## Acknowledgments

We thank S. Girard, B. Houchmandzadeh, T. Mora, R. Ranganathan, and A. Walczak for discussions, and K. Reynolds for help with sequencing. This work was supported by an AXA Research Fund Postdoctoral Grant (to D.B.) and Agence Nationale de la Recherche Grant ANR-10-PDOC-004-01 (to O.R.).

## Footnotes

↵

^{1}D.B., A.K.S., and N.S. contributed equally to this work.- ↵
^{2}To whom correspondence should be addressed. Email: olivier.rivoire{at}ujf-grenoble.fr.

Author contributions: S.B., C.N., and O.R. designed research; A.K.S. set up phage display; A.K.S. and O.R. designed the libraries; D.B. constructed the libraries; S.B. performed the selections; N.S. and C.N. supervised the experiments; S.B. and O.R. analyzed data; and S.B., C.N., and O.R. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1517813113/-/DCSupplemental.

## References

- ↵.
- Magurran AE

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Fellouse FA,
- Wiesmann C,
- Sidhu SS

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Ravn U, et al.

- ↵.
- Mehta ML

- ↵.
- Gümbel EJ

- ↵.
- Coles S

- ↵
- ↵
- ↵
- ↵.
- Weinstein JA,
- Jiang N,
- White RA 3rd,
- Fisher DS,
- Quake SR

- ↵.
- Mora T,
- Walczak AM,
- Bialek W,
- Callan CG Jr

- ↵.
- Desponds J,
- Mora T,
- Walczak AM

- ↵
- ↵.
- Gillespie JH

- ↵.
- Lancet D,
- Sadovsky E,
- Seidemann E

- ↵
- ↵.
- Beisel CJ,
- Rokyta DR,
- Wichman HA,
- Joyce P

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Young SS,
- Sheffield CF,
- Farmen M

- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics

- Biological Sciences
- Biophysics and Computational Biology