Unrealizable target distributions

Let P_target = (P_target(¹), ? , P_target(21)) be the amino acid distribution at a target position of polypeptide, where 0 ? P_target(a) ? 1 with a = 1?20 gives the probability for amino acid a, and 0 ? P_target(21) ? 1 gives the probability of having a stop codon at the target position. These probabilities are subject to the constraint that ?_{a=1 to 21}P_target(a) = 1. Let P_n1 = (P_n1(A), P_n1(C), P_n1(T), P_n1(G)) be probability distribution over the nucleotides at the first position of the codon for the target. Similarly, P_n2 = (P_n2(A), P_n2(C), P_n2(T), P_n2(G)) is defined at the second position, and P_n3 = (P_n3(A), P_n3(C), P_n3(T), P_n3(G)) at the third position. If we assume that the three codon positions be independent from one another, the probability of having amino acid a generated from the nucleotide distributions should be

where n1, n2, and n3 stands for the nucleotides, the summation exhausts over all possible nucleotides, and ?(a | n1 n2 n3) = 1 if codon ?n1 n2 n3? encodes amino acid a and ?(a | n1 n2 n3) = 0 otherwise. The Equation (1) defines the mapping from nucleotide distribution space (a 12-dimension hypercube) to the amino acid distribution space (a 21-dimension hypercube). So, even if the mapping is injective, or one-to-one, the nature of mathematical functions makes it clear that there will be infinitely many points in the 21-dimension hypercube that cannot possibly be covered.

In practical laboratory settings, the precision for the quantities of various reagents is finite, and for the sake of discussion simply let us assume it down to 1%, as suggested in LaBean and Kauffman (⁷). In other words, these probabilities [Equation (1)] take integral values from 0 to 100 inclusive for each nucleotide and each amino acid. This then allows us to estimate the extent to which the target amino acid distributions cannot be matched by any nucleotide distributions via Equation (1). Since now P_n1(A), P_n1(C), P_n1(T) and P_n1(G) take an integer value between 0 and 100 each and sum to a total of 100, the number of possible assignments is C₁₀₃¹⁰⁰, i.e. 103 choose 100. As 3 nt form one codon, this value is then cubed yielding 5.53 ? 10¹⁵ possible codon nucleotide distributions. Similar analysis suggests that the number of possible amino acid distributions down to 1% accuracy is C₁₂₀¹⁰⁰, which is about 2.94 ? 10²². As argued above, codon distributions map into amino acid distributions, so even if the mapping is injective, or one-to-one, then at a maximum only 5.53 ? 10¹⁵ of the 2.94 ? 10²² amino acid distributions can be matched (⁷). The mapping is clearly not one-to-one because of the degeneracy of the genetic code. So, at most, 1 out of every 5.3 million amino acid distributions can have an exact match via Equation (1).

Objective functions

Because, as shown above, no exact match is achievable for many target amino acid distribution, it is a desirable compromise to find codon nucleotide distributions that match a given target amino acid distribution as closely as possible.

Objective function 1 (least squares) (¹⁰): one simple way to measure how close two amino acid distributions are is by using the difference between the calculated and the target:

where P_target (a) is the target distribution?s value for amino acid a, P_calc (a) is the calculated distribution as given in Equation (1) and wt(a) is a weighting factor for amino acid ?a? to affect how a mismatch is factored into the total measure D. Then, the task of finding the best matching distribution is choosing nucleotide distributions that minimize the difference.

This simple intuitive measure is called the first objective function in (⁵), where four other objective functions were also used, and adopted in this study as well. These four other objective functions may be selected by the user and are listed below.

Objective function 2 (cubic function) (¹⁰):

where Objective function 3 (maximum likelihood) (¹⁰): Objective function 4 (Cosine Basin) (⁸): Objective function 5 (chi square and relative entropy) (⁵): where ? is an arbitrary small constant introduced to avoid numerical instability.

Codon usage

To further enhance the yield, the codon usage information for the host organisms where the protein libraries are synthesized should be taken into account. Let k_a(n₁ n₂ n₃) be the frequency that codon ?n₁ n₂ n₃? are used for amino acid a, and then with codon usage being factored in the adjusted target distribution becomes P_target(a |n₁ n₂ n₃) = P_target(a) k_a(n₁ n₂ n₃). We can replace the target distribution in all these five objective functions with the adjusted target distribution to take into account of codon usage. For example, objective function 5 [Equation (6)] is modified as follows.

The codon usage information can be found in various databases, such as the one at the Codon Usage Database (http://www.kazusa.or.jp/codon/). If codon usage is not specified, the default codon usage will be a uniform distribution among all the codons.

Optimization

One major improvement from previous work (⁵) is that we make use of a genetic algorithm (¹¹) to perform the optimization instead of using gradient descent method, which is known to be more susceptible to get stuck in local optima. Another major contribution of this work is detailed analysis of the constraints imposed by the genetic code and a method for circumventing these constraints by using more than one codon nucleotide distribution to match a given amino acid distribution. When multiple codon nucleotide distributions are used, the genetic algorithm, with proper encoding schemes (see below), can also optimize on the weighting factors for these distributions; these weighting factors are free parameters in the objective functions and cannot be easily optimized otherwise.

Optimization using genetic algorithm

The main evolutionary algorithm, called CodonOptima, is given by the following pseudocode.

Input: population size ps, objective function f, fitness_threshold ft, maximum time mt, maximum iterations mi, mutation rate mr, mutation density md, crossover rate cr, double crossover rate dcr

Algorithm

1. Generate a population p, of size ps, of codon nucleotide distributions at random
2. Evaluate the fitness of each member using the chosen objective function, f
3. Terminate, if the fitness of the best member is below a certain stopping fitness threshold, ft, or the maximum time, mt, or number of iterations, mi, has been exceeded.
4. Remove those with low fitness and replace with new random members. Mutate some of the population, where mutation density, md, is the probability of mutating a single member and mutation rate, mr, is the probability of mutating a single position in a member that has been selected for mutation. Crossover and double crossover also occurs at the rates specified by cr and dcr respectively. Go to step 2

Output: p, the population of codon nucleotide distributions sorted by objective function fitness and/or R-value

Encoding of distributions and implementation of mutations and crossing over

The distributions are encoded as strings representing the relative percentages of ACGT for each of the three codon positions. For the purpose of wet lab preparations, a precision of 1% for any one nucleotide is sufficient. Therefore, the percentage for each nucleotide is an integer between 0 and 100 inclusive, which takes 7 bits in binary to represent. For convenience, we use 4 bytes (1 byte each) to encode the amounts for the 4 nt at a codon position, and normalize these four integers to get the relative percentage for the corresponding nucleotides. For example, the binary string ?01001000 11101010 00001000 10000010? give the following amounts (A = 72, C = 234, G = 8 and T = 130), and relative percentages (A = 16%, C = 53%, G = 2% and T = 29%) after normalization. While the ?one-byte-per-nucleotide? representation suffices for the purposes of wet lab preparations, it can be expanded to accommodate higher precision when necessary. Since there are three positions for a codon, this representation has 12 bytes in total. Even at this 1% precision, the search space of codon distributions is already of a size of ?10¹⁶.

New distributions are generated by ?evolutionary operations? on the strings representing the existing distributions. Specifically, these evolutionary operations include point mutations and cross-overs. Point mutations are simply flips of single bits in this 12-byte string from ?0? to ?1?, or vice versa. Typically, this results in a small change in the value that a byte represents, particularly when the mutation happens at the less significant bits. These small changes in value allow for fine-tuned explorations of the codon nucleotide search space. Larger jumps through the search space are accomplished by means of cross-over, where two strings will swap their characters after a random pivot point. Double cross-over is also implemented with two pivot points to allow for ?middle? sections to be swapped between a pair of strings. Mutation rate and density values typically range from 0.1 to 0.35, with cross-over and double cross-over rates being usually 0.8 and 0.01, respectively. In Figure 1, point mutation and cross-over are illustrated on toy examples with strings of 16 bits.

Multi-test tube solutions

Not only a target distribution may be impossible to match exactly as shown in the section ?Unrealizable target distributions?, but it is often also difficult to find the exact match even when it exists. The situations are further complicated due to the constraints imposed by the genetic code and codon usage. To alleviate the difficulty, in a practical lab setting, it is acceptable to use two separate codon distributions (n1 n2 n3) and (n1? n2? n3?) in combination (i.e. two test tubes) to match a target amino acid distribution, as specified in the following equation.

We implemented a combined solution of two test tubes (one for each codon). This approach can be further generalized to using multiple codon nucleotide distributions to match any given target amino acid distribution. As shown in the ?Results? section, this relaxation has resulted in better performance as measured by relative entropy.

This computational method of multiple test tubes can be implemented using wet lab methods. Previous work on probabilistic protein design has demonstrated using automated phosphoramidite synthesizers to create target nucleotide sequence ensembles according to the calculated arbitrary nucleotide probabilities (⁶). Mixtures of custom designed oligonucleotides from commercial sources can also be utilized to match target distributions (L. Shi et al., 2009 submitted for publication). Wet lab methods may also depend on the number of positions in a polypeptide that are to be matched, the complexity of the amino acid distributions at each of these positions and how far apart the positions are. For example, if the positions are far enough apart, methods such as overlapping Polymerase Chain Reaction can be used to assemble two different regions together. In general, it is preferable for wet lab methods that the number of generated codon nucleotide distributions (i.e. test tubes) be as few as possible.

What is the minimum number of test tubes that can be used to guarantee that any given target amino acid distribution can be matched? The trivial upper bound is 20 as every amino acid that is to be matched can be created from its own codons in a separate test tube. A more stringent and realistic upper bound is 6. To explain how this bound is established, let us examine the constraints imposed by the genetic code and identify amino acids that can be grouped in one test tube and still can have any distribution over them matching exactly. It is assumed that the desired distribution will contain no stop codons. For example, by inspection, if we wanted to match a target amino acid distribution over alanine (A), aspartate (D), glutamate (E), glycine (G) and valine (V), we can see that they all lie in the last four rows of the following genetic code table (Table 1).

That is, all these five amino acids have codons with a G in the first base position. Therefore, a mixture of these five amino acids (ADEGV) can be encoded by degenerate codons GNN, where N stands for any nucleotides. The relative amounts of U, C, A and G in the second codon position control the amounts of valine, alanine, aspartate/glutamate and glycine, respectively. The third base in the codon is degenerate for valine, alanine and glycine and does not affect their amounts, but the ratio of (U + C) to (A + G) controls the relative amounts of aspartate to glutamate. As a result, with a codon distribution given by GNN (100% G in first position), we can exactly match any given distribution for ADEGV. This distribution can then be weighted accordingly to match the initial distribution containing up to 20 target amino acid frequencies.

Based on this analysis, a list of one-test tube groupings of amino acids, which can be matched perfectly and will not contain any stop codons, can be generated and is given in Table 2. For each grouping of amino acids, the corresponding degenerate codons are also listed using the nomenclature for variable bases from Cornish-Bowden (¹²): H = not G, Y = T or C, R = A or G, etc.

From these groupings, it is easy to propose the following partitioning of the 20 standard amino acids using six test tubes. One test tube for each of these four groups: IKMRT, ADEGV, HLPQR and CFSY. These four test tubes contain 18 unique amino acids, ACDEFGHIKLMPQRSTVY. Asparagine can be matched by codon AAY in one test tube (fifth) and tryptophan by UGG in its own test tube (sixth).

where in w_i, i = 1?6 are the weights to control the relative percentage of the six individual distributions from P(IKMRT) to P(W) in a target distribution. This recipe gives a total of six test tubes to exactly match any given target amino acid distribution.

Six test tubes are sufficient for exact matching, but are they necessary? The answer is yes, and the proof is outlined as follows.

1. Construct a 20 ? 20 binary valued matrix, M, whose rows and columns correspond to 20 amino acids, and each entry is either 1 if the corresponding amino acids can be exactly matched using one codon nucleotide distribution (test tube) or 0 if otherwise. An amino acid pair can be matched exactly if there is at least one codon for each amino acid in the pair such that these two codons have a Levenshtein edit distance (¹³) of exactly 1.
2. Build a graph G based on matrix M. G contains 20 nodes corresponding to the 20 amino acids. Any two nodes x and y in G are connected by an edge if and only if M(x,y) = 0. That is, edges in G correspond to pairs of amino acids that cannot be exactly matched using one codon nucleotide distribution so they must be in separate test tubes.
3. Find the maximal clique of G. All of the amino acids in this clique must be in mutually exclusive test tubes.

Two maximal cliques of size 6 are found in graph G, AFMNQW and EFMNPW. For both of these maximal cliques, the six amino acids must be placed in separate test tubes giving a lower bound of 6 on the minimum number of test tubes needed. These groupings may prove useful to wet lab researchers employing probabilistic protein design techniques using codon nucleotide distributions. While most amino acid distributions can be matched sufficiently with a few test tubes, amino acid distributions containing many or all of the amino acids in these maximal cliques may prove very difficult to match if additional test tubes are not utilized. Since 6 has been demonstrated as both a lower and upper bound on the minimum number of test tubes needed to match any given amino acid distribution, then Equation (9) represents one such optimal partitioning using a minimum number of test tubes.

[Formula ID: M10]
10