Computational methods have always played an important role in protein design. This work focuses on searching the protein sequence space to find one or several protein sequences that are compatible with known structures and functions. Probabilistic computational methods provide information on the range of amino acid variations allowed within desired functional and structural constraints. The source for this experiment is "A Guide to Modern Protein Engineering Experiments" [German] K.M. Arndt, K.M. Miller, eds.
Operation method
Combinatorial protein design strategies using computational methods Move Probabilistic methods in protein design are those that provide probabilistic estimates of the probability of an amino acid (appearing) at a one-finger localization point in a given protein structure. Here, we discuss several methods for estimating these probabilities and focus on entropy-based self-consistent formulas that directly solve for these probabilities. For more product details, please visit Aladdin Scientific website.
2.1 Comparison of related sequences
Sequence variability in protein structures can be explored with sequence and structure databases. Sequences known to fold into very similar structures can be identified from protein databases or structure comparison databases [48] . If the structure of a sequence is known, sequences with sufficient sequence similarity [e.g., sequence identity greater than 40%] can be considered to share the same structure. Multiple sequence comparisons of such structurally similar proteins allow the probability of amino acid site-specificity to be estimated simply as the frequency of each amino acid at each position (occurrence) in the comparison [49] . Such a set of probabilities is often called a sequence profile. If the number of sequences is insufficient such that certain amino acids never occur at particular sites, psedocounts and other methods can be used to regularize these frequencies so that they are more representative of folding into the selected structure. Nonetheless, the probabilities obtained from such profiles would seriously bias the nature of the sequences in the database. Because there are a large number of examples of sequences with very low similarity collapsing into similar structures, we would like to get a complete understanding of sequence variability in a broader context. Profiles derived from databases are also not suitable for designing new protein structures for which there is no sequence in the database. Using a given backbone structure as a template, more general computational methods can determine amino acid probabilities from scratch.
2.2 Directed search methods for building profiles
Repeated applications of directed search methods can estimate the nature of a sequence as a whole. For this type of calculation, a target structure is selected by giving the coordinates of the main chain atoms. If a single protein structure is used, several recent direct design studies have yielded sequences that are quite similar to wild-type sequences [51~54 ]. For a given structure, multiple sequence search calculations can be run independently to obtain a set of sequences whose comparison yields site-specific probabilities.Desjarlais and co-workers, ran their sequence prediction algorithm independently for each member of a system of structures that are extremely correlated with a particular fold consistency [ 55 ]. For each structure, an optimized 'nucleation' sequence is identified, and subsequently, sequence/rotation anomalies are explored for the entire structure. This approach has been used to identify sequences that are compatible with the folding of the small β-sheet WW structural domain [4, 55]. Sequence prediction algorithms were applied to each of the 100 microstructural variants (1A rms) of a particular fold to construct computational profiles that were much more discrete than applying individual structures [ 56]. Workers at Xencor Inc. sampled an optimized sequence (in which residues near the β-lactamase active site were replaced) using Monte Carlo sampling [ 57], and found that an antibiotic that was not compatible with the small β-sheet WW structure domain [ 4, 55] was not compatible with the small β-sheet WW structure domain [ 4, 55]. The Monte Carlo sampling method [ 57 ] was used to find a sequence that increased resistance to an antibiotic by more than 1000-fold. However, these methods of constructing profiles are very computationally intensive. This is because, in order to establish amino acid site-specific frequencies, repeated directed searches are completed.
2.3 Statistical theory of sequence lineages
Statistical, entropy-based formulas have been developed to identify a set of site-specific amino acid probabilities for a given backbone structure, not just the optimal sequence [ 58, 59 ]. The theory derived from statistics is used to deal with the number and composition of sequences compatible with the main chain structure. This theory also deals with the entire space in which the composition fits, not just the small portion of the space accessible to experimental and numerical calculations and sampling. The properties of suboptimized sequences are easily tested. Large protein structures (more than 100 residues) are easy to calculate. The 'entropy' here is the number of sequences that are compatible with the target structure. This concept, derived from thermodynamics, is used to reduce the number of possible sequences: the restriction of sequences reduces the entropy and is accompanied by a decrease in energy.
The inputs in the method are the target main chain structure, and an energy function quantifying sequence-structure compatibility. For a target main chain structure, the method produces a probability for each amino acid (occurrence) at each residue site (see Note 2 ). In theory overall properties (e.g., the total energy of the sequence in that target) and local properties (e.g., the amino acids allowed at a particular site), both can be included as restrictions in the method. A collection of many amino acid probabilities is possible. The method determines the "most probable" ("most likely") such sets by maximizing the effective entropy, whereby the maximization is restricted. The method effectively provides the means for the system to reduce the volume of sequence space to be searched to experimentally accessible levels through such a restriction.
In a sequence with desired properties specified by the restriction function, let Wi denote the probability that the α amino acid appears at residue position i and such that its side chain is any one of a set of discrete conformations - rk (α) ( rotational isomorphisms; cf. [6] and [ 60] ). The total sequence-configuration entropy one Sc (here simply referred to as "conformational entropy]) can be defined as 
The summation is over every sequence site i and all possible amino acids α. For each amino acid, the summation is also over k possible rotational isomorphisms - rk(α). Wi is obtained by maximizing Sc subject to the restriction f. The maximization is done by the Lagrange multiplier method [ 61]. The variational general function V of Wi[α, rk(α)] is defined as 
In general, the restriction f is also a function of the probability Wi[α, rk (α)] . In determining the probability of a state compatible with a particular restriction, the mth restriction function fm is bounded to take the value fm0. The system of equations and the Lagrange multipliers for determining the probability take the form (see Note 3): 
This large system of coupled nonlinear equations is solved by the root- finding method. Although there are many options for such a method, we find an overall convergent method that can be widely adopted [ 62].
2.3.1 Energy Functions
Two energies, the conformational energy Ec and the environmental energy Eenv, are considered in the calculations and are used as constraints in maximizing the conformational entropy.
The conformational energy Ec is calculated using the atom-based potential energy one AMBER force field [63]. Ec includes van der Waals interactions, electrostatic interactions with a distance-dependent dielectric constant (4ery ), and a modified hydrogen bonding term. For a particular sequence ( α1, ..., αN ), where the conformational states of the amino acids are [ r1 ( α1 ), ..., rN ( αN ) ], the Ec is 
When considering the protein energy function, the one-body term εi [ α, rk(α)] includes the interactions of the main chain and side chain atoms, as well as the reference energy of the amino acid (see subsection 1.2.3.3). The two-body term is the sum of the interactions between two rotational isomers at two different sites in the structure. For a large number of sequences enjoying common energetic properties, we assume that Ec rises and falls around its mean value due to sequence changes are not significant. Then, we can write 
As another constraint, the ambient energy is introduced to account for hydrophobic effects in an equivalent way within statistical theory [59] . This potential energy takes into account the surface exposure propensity of amino acids. We can write Eenv in terms of amino acid probabilities as 
where is the local environmental energy defined in subsection 1.2.3.2. Note that this energy does not include two-body interactions and depends only on the amino acid and rotational isomerization state at each position.
2.3.2 Solvation and hydrophobicity energies
Quantifying hydrophobic interactions and other solution properties is an important parameter in protein design approaches. It is impractical to computationally test for large amounts of variation in sequence, and even calculating the solution-accessible surface area, which often correlates well with hydrophobicity, can be computationally intensive. In an effort to consider solution effects in a practical way that is consistent with statistical calculations, the environmental energy was introduced as a function of the density of β carbon atoms near each site, ρ [59] . In general, hydrophobic residues tend to localize in the buried regions of proteins, whereas hydrophilic residues tend to localize on the surface. Consequently, hydrophobic residues tend to have a higher density of β carbon atoms than hydrophilic residues. Using 500 different globular proteins of known structure, we derived a generalized "statistical" potential energy equation for calculating the effective potential energy of amino acids 
where p ( α, ρ ) is the localized β-carbon atom density of residue α observed as ρ the number of times residue α is observed in the training set; p ( α ) is the number of times residue α is observed in the training set; p ( ρ ) is the number of times, regardless of residue type, that the local density ρ is the number of times the local density ρ is observed regardless of residue type; Te is the effective temperature; density ρ is the density of β carbon atoms in the "free volume" centered on a particular orientation of the residue. The free volume is the average volume not excluded by the side chains 
where nβ is the number of β carbon atoms within a distance R (e.g., 8A) from the center of mass of the side chain. We note that the local density depends on the rotational isomerization state of the residues. This β-carbon-atom density-based local energy correlates well with other amino acid hydrophobic scales [59]. For sequence probability calculations, the Eenv restriction takes the value of a known sequence with the same structure (if there is a known value), or a representative value of a protein with the same size or chain length.
2.3.3 Reference energies
In protein design, we seek to optimize the energy of a particular sequence in the target structure relative to the unfolded state tether. To deal with the unfolded state, a reference energy is introduced into Ec for each amino acid to mimic the effect of the inactivated state [ 51, 65]. This energy is calculated as the 'free energy' of each amino acid α in the N-acetyl-α-N'-methylamine amino acid form, averaged over multiple backbone structures. This is a rough approximation of the average value for the stretched unfolded state. The reference energy consists of summing the possible rotational isomerization states and the possible main chain conformations. The values of ψ and φ for the main chain conformation are taken at 10° intervals. The reference energy for each residue can be estimated using the following formula: 
where Emf is the conformational energy of amino acid α in the N-acetyl-α-N'-methylamine form in a particular conformation determined using the molecular site energy. Here, a database of rotational isomers dependent on the main chain is always used. The reference energy is measured relative to glycine (G ) without side chains. Energy constraints on the main chain include interatomic interactions and take the following form 
2.3.4 Rotational conformation and (amino acid) identity probabilities
This theory maximizes the conformational entropy Sc to obtain the probability Wi that a particular amino acid appears at site i in side-chain conformation k. The amino acid probability Wi ( α ) is available as 
be determined.
By analogy with statistical thermodynamics, the Lagrange multiplier, which arises by limiting the conformational energy, can be considered as the equivalent inverse temperature. The corresponding heat capacity C is defined as 
Applying this theory to a specific protein, the SH 3 structural domain, serves as an example. As the effective temperature Tc decreases (i.e., Ec decreases; Figure 1.1), the conformational entropy decreases. At high energies (high Tc, low βc), there are many unfavorable interactions between residues (high energies) and a broad distribution of sequence/rotational isoforms at each site. In general, the number of possible amino acids and rotational isomers at each site decreases with energy. As shown in Figure 1.1, Cv passes a peak at Tc ( = 1/βc ) = 10 mol/kcal before reaching a trough at approximately 2 mol/kcal. At this point the internal residue types and conformational states are relatively well defined [59]; whereas the surface residues, despite being predominantly hydrophilic, have a large number of rotational isomorphisms with comparable probabilities. This is consistent with the conformational variability of surface-exposed residues of proteins in the database. Thus, for determining at which 'effective temperature' to test the probabilities, the heat capacity is helpful. In addition, direct comparison with profiles obtained by sequence comparison (e.g., protein secondary structure databases derived using homology; ref. [66]) yields very compatible results, especially in the buried region (Fig. 1.2). 
2.4 Genebanks for protein profiling
In subsection 1.1.3, we discussed how to apply structure-specific sequence profiling. If a sequence is identified by a public or targeted search sequence, it can be directly implemented by peptide synthesis or by synthesizing the gene encoding the sequence and then expressing it. Large proteins are often realized by expression. If probabilistic sequence information is to be used to build combinatorial libraries (see Note 4 ), methods are needed to transcribe protein profiles into partially randomized libraries of gene sequences. A non-uniform distribution of nucleotides is necessary to encode polypeptide sequences that are biased toward specific amino acids. The pseudo-independent nucleotide probabilities at each position of a set of partially randomized genes can be determined computationally such that the gene library encodes a protein library that best reproduces the desired amino acid profile. The calculated gene pool can then be used in standard DNA synthesis.
If the various codes for amino acid α are equally probable (no coding bias), the probability of amino acid α is the sum of the probabilities of the various codes corresponding to this amino acid: 
The objective function quantifies the difference between the desired amino acid probability distribution and the probability distribution of amino acids encoded by a given nucleotide probability distribution [ 67, 68]. In order to find the nucleotide probability that best reproduces the desired amino acid frequency and avoids termination codes, Wang et al. proposed a new objective function "69". This objective function consists of two terms, one is the X2 function which quantifies the absolute deviation between the expected and calculated amino acid probabilities, and the other is the relative entropy. Such relative entropy is widely used to quantify the "distance" between two probability distributions, and is a strong indicator that information from one distribution is not included in the other [ 50 ]: 
Figure 1.3 illustrates the nucleotide design at a specific amino acid site of a particular protein, here site 54 of the SH3 structural domain. Shown are the expected frequency of amino acids (hollow band on Figure 1.3) and the frequency of amino acids encoded by the nucleotide probabilities used for the calculation (solid band on Figure 1.3). The two match well in this example. Because of the partial parsimony of the amino acid coding, in many cases the expectation and the calculated probability distribution do not precisely match. This calculation provides excellent complete sequences (without termination codes): for 50-60 residues of a test protein that has undergone randomized selection, the output of the complete sequence would be 96% or more. The high output is particularly impressive when most or all of the gene has undergone combinatorial substitutions. 