CMCpy: Genetic Code-Message Coevolution Models in Python.

Code-message coevolution (CMC) models represent coevolution of a genetic code and a population of protein-coding genes ("messages"). Formally, CMC models are sets of quasispecies coupled together for fitness through a shared genetic code. Although CMC models display plausible explanations for the origin of multiple genetic code traits by natural selection, useful modern implementations of CMC models are not currently available. To meet this need we present CMCpy, an object-oriented Python API and command-line executable front-end that can reproduce all published results of CMC models. CMCpy implements multiple solvers for leading eigenpairs of quasispecies models. We also present novel analytical results that extend and generalize applications of perturbation theory to quasispecies models and pioneer the application of a homotopy method for quasispecies with non-unique maximally fit genotypes. Our results therefore facilitate the computational and analytical study of a variety of evolutionary systems. CMCpy is free open-source software available from http://pypi.python.org/pypi/CMCpy/.

Introduction

Code-Message Coevolution (CMC) models were introduced to facilitate the study of co-evolutionary systems in which a large population of individuals share an evolvable genetic code coupled to many and/or long protein-coding genes (“messages”). Initial study of these models demonstrated that base mutations in protein-coding genes can significantly influence the fitness of genetic codes.1 More detailed study of CMC models has yielded formal demonstrations for how natural selection likely contributed to the origins of codon redundancy2 and non-random patterns of amino acid assignments3,4 in genetic codes. CMC models have been extended to study the effects of population structure and gene sharing on coevolving systems of codes and messages.5

In their original formulations, CMC models are deterministic evolutionary genetic models that couple together sets of large populations of asexually reproducing genotypes evolving under mutation and natural selection. One such population is called a “quasispecies”,6 a concept reviewed recently both generally7 and specifically in connection to population genetics theory.8,9 In CMC models, each of multiple quasispecies represents a large population of codons evolving to meet the same physicochemical requirements of a “site-type” in proteins under translation by a given genetic code. The coupling of multiple quasispecies through a genetic code occurs through reuse of the same codon type (or allele) in different site-types. CMC models follow deterministic trajectories that alternate between equilibration of messages to an established genetic code by mutation and selection, and locally adaptive “gradient ascent” hill-climbing of the genetic code through single codon assignments and reassignments. This “quasistatic” dynamic continues until no single codon reassignment yields higher fitness with the current message population. Published CMC models always converge to a stable local fitness optimum—a process called “code freezing.”

Because CMC models compound quasispecies models, analytical solutions to quasispecies models are relevant to their study. Previous applications of perturbation theory for approximate solutions to quasispecies models assume uniqueness of maximally fit genotypes and small perturbation parameters (small mutation rates).10–12

In this work, we present two analytical methods that relax these assumptions. The first is a perturbative method that clarifies, extends and generalizes the perturbative approach to quasispecies models. Our presentation provides a derivation to arbitrary order, which relaxes any restriction on the magnitude of the perturbation parameter. The second method is the first published application of homotopy methods to quasispecies models, which handles cases with non-unique maximally fit genotypes. Numerical tests show that, for two example problems, both algorithms converge. Notably, the error of the perturbative method decreases exponentially in the number of iterations, suggesting its potential to outperform the power method. Both approaches are flexible and can be extended to a variety of quasispecies mutation models and fitness schemes. Also, both methods have been formulated so as to facilitate future connections between quasispecies models and established results in eigenvalue perturbation theory13–15 and homotopy methods.16

In addition to these analytical results, we present a new code-base implementing CMC models. The dissemination and further study of CMC models has been sorely hindered by the lack of a modern code-base implementing them. The original C/C++ CMC code-base has been rendered obsolete by evolution of compilers and language standards from the time of its writing in the late 1990s. This has created a need for reimplementation of the original models in an easy-to-use and program, powerful, and efficient, high-level scripting language like Python.

Here we present CMCpy, a free open-source code-base that implements CMC models in an easy to use object oriented Python API. The code-base comes with a front-end command-line executable called cmc that can drive the exploration of a variety of CMC models and reproduce published results.

Implementation

CMCpy was developed in Python 2.7. Figure 1 shows the organization of classes in CMCpy. A rich class hierarchy allows convenient specification of a wide-range of CMC models with very few lines of Python code.

Figure 1

Overview of class hierarchy and containment relationships in CMCpy.

An element not shown in Figure 1 is that the ArdellSellaEvolver class is an abstract base class for a variety of subclasses implementing different strategies to solve dominant eigenpairs. The default solver relies on the Numpy eig() function for its speed and accuracy, but for reasons of flexibility, verification, and to experiment with different methods, we include a legacy central processing unit (CPU)-based power method, a multicore CPU-based power method implementation and an experimental graphics processing unit (GPU) implementation of the power method that relies on pyCUDA.17

The pyCUDA solver is implemented in a CUDA C kernel with supportive Python code. The Python code uses NumPy for simple operations, to cast data to types acceptable by the CUDA platform, and to reshape matrices to one-dimensional data frames, for use by CUDA C. The CUDA C implementation of the power method approximately solves the eigenvector of each site-type, in parallel. Each site-type is assigned to one virtual “block” of the GPU, each of which corresponds to a physical processor when the number of blocks is below a hardware limit. This organization of the parallel workload was chosen to stay below CUDA specification limits, access GPU memory efficiently, and avoid excess complexity. Ultimately, the CUDA power method implementation is faster than the CPU power method implementation. Figure 2 shows system clock execution times of various methods on a 3.0 GHz Core 2 Duo with an Nvidia GeForce 460 GTX GPU running Ubuntu 12.04. The benchmark script in R that generated this figure is provided as supplementary data.

Figure 2

Comparison of wall-clock execution times of the cmc executable with three different eigensystem solvers using the double ring model4 with eight codons, μ = 0.1 and ϕ = 0.25.

CMCpy comes with a front-end command-line executable in Python called cmc. This executable provides users, including non-programmers, the capability to reproduce (at least qualitatively) all of the published results on CMC models1–4 as well as run individual and batch simulations of other models and parameter spaces. In Table 1 we list command-line options to the cmc executable and their corresponding model parameters.

Table 1

Options to the cmcexecutable and corresponding model parameters.

Option	Long version	Description	Default
–a	––numaas	Number of amino-acids/site-types (AA/ST)	10
–d	––numdims	Dimensionality of AA/ST space	1
–s	––seed	Seed for random initialization of AA/ST coordinates	42
–t	––numtrials	Num. trials with reinitialized AA/ST spaces	1
–c	––numcodons	Num. codons in a “double ring” model	N/A
–b	––numbases	Alphabet size for word-based codon model	4
–p	––numpositions	Length of codons for word-based codon model	N/A
–m	––mu	Base (word models) or codon mutation rate	0.1
–k	––kappa	Transition/transversion mutation bias ratio	1
–f	––phi	Missense tolerance parameter	0.25
–r	––misreading	Misreading parameters	N/A

A variety of observables and statistics are implemented including the Normalized Encoded Range2 for one-dimensional amino acid/site-type spaces. Although CMC models are deterministic, exact quantitative differences may arise in results with CMCpy depending on floating point representation differences by platform and differences in convergence thresholds using power method-based eigensolvers.

It may be useful to restate the assumptions of the “Ardell-Sella” models currently implemented and available in CMCpy. These include the following:

The numbers and/or lengths of protein-coding genes, or “messages,” translated by a common genetic code, are large with respect to every possible site-type in proteins.

For every possible site-type in proteins there corresponds a uniquely most fit amino acid.

The machineries to decode codons, and to associate any codon to any amino acid, pre-exist the evolution of the genetic code.

Fitness contributions of amino acids across sites are independent and multiplicative.

Fitness contributions of different amino acids within the same site are independent and additive.

Bases in messages mutate independently of one another.

Messages are haploid and asexually reproducing.

Genetic codes evolve much slower than messages, through discrete and independent assignments or reassignments of amino acids to codons.

Analytical Methods for Quasispecies Solutions

In this section we develop two different analytical methods to solve for the equilibrium growth rate and genotype distributions for a wide range of quasispecies models. The Matlab/Octave code used to implement these solutions are provided as supplementary data. A version of the homotopy method is also implemented in CMCpy for ring models. Full implementions of both methods will be incorporated into CMCpy at a later date.

Perturbative method: quasispecies with unique fittest genotype

We start with a ring mutation model discussed in prior work.1 Let μ denote the N × N mutation matrix where blank spaces should be interpreted as zeros. Throughout this derivation and the following one, we assume 0 ≤ μ < 1 and N > 1. Let w denote the N × N fitness matrix

We assume that w1 is the unique maximum of the finite set {w1, . . . , w}. We also assume wi > 0 for each i, consistent with prior work.1

Our goal in this derivation is to determine the leading eigenpair, consisting of the largest eigenvalue together with its corresponding eigenvector, of the matrix

We write Q̃(μ) to emphasize the fact that since μ is a function of μ, so is Q̃. The diagonal matrix w is constant with respect to μ. Next we define

Since the matrices Q and Q̃ are related by a similarity transformation, their eigenpairs are closely related. One can check that

We focus our attention on Q because it is symmetric, i.e., Q(μ) = Q(μ) for all μ.

For Q(μ), the leading eigenpair (λ(μ), v(μ)) must satisfy

When μ = 0, the matrix μ reduces to the N × N identity matrix. Therefore, Q(0) = w, and the leading eigenpair of Q(0) is given by

Here e denotes the j-th basis vector in N-dimensional space, i.e., the vector with all zeros except for 1 in the j-th slot.

Guiding principle

Since Q(μ) is symmetric for all real μ, standard theoretical results in eigenvalue perturbation theory14 guarantee that both the eigenvalue λ(μ) and eigenvector v(μ) are analytic functions of μ. This means that there exists M > 0 such that for μ ∈ (−M, M), the following power series expansions converge:

In this notation, (6) can be written as λ0 = w1 and v0 = e1. Note also that we have arranged the coefficients so that

Our strategy now will be to derive from (5) a recursive set of equations for the coefficients λ, μ. Once we have these coefficients for j = 0, 1, ..., J, we have an approximation to the leading eigenpair.

Let I denote the n × n identity matrix. Then μ = I + μs, where

This implies that

By (3), we have

Perturbative solution: part I

Armed with the above facts, we differentiate (5) with respect to μ on both sides:

We then set μ = 0 and use (8), (3), (9), and (6) to obtain

Multiply this equation on the left by the row vector :

Note that , which implies . Using this equality in (11), we have

This shows that if we already know (λ0, v0), we can determine λ1. To determine v1, we return to (10) except now we treat λ1 as known. Rearranging the equation, we have

We substitute our definition of w on the left-hand side to obtain

Suppose that v1 = (v1(1), v1(2), ..., v1()). We set v1(1) = 0. For the remaining components, we solve the above matrix-vector system to obtain

We have completed the loop, showing how to proceed from the zeroth-order eigenpair (λ0, v0) to the first-order eigenpair (λ1, v1).

In the next section, we show how to iterate this procedure to generate the j-th order eigenpair (λ, v) from the previously obtained eigenpairs.

Perturbative solution: part II

We return to (5) and take j derivatives with respect to μ on both sides. Using the general Leibniz rule, we have

After differentiating, we set μ = 0 and use (8) to obtain

Applying (3) and (9) yields

We peel off the m = 0 term from the right-hand side:

Multiplying on the left by and using , we see that the first term on the left-hand side cancels the first term on the right-hand side. We are left with

Now on the right-hand side, we peel off the m = j term—note that this the only term in which λ appears. Hence

We use (6) and , and we solve for λ:

We have therefore shown that if we already know (λ, v) for m = 0, 1, …, j – 1, we can solve for λ. Now, using this λ, we can solve for v in the same way as before. We go back to (12), isolate all terms involving v, and apply (6) to derive

The right-hand side is clearly valid only for j ≥ 1. The matrix on the left-hand side is the same one that appeared earlier:

We again set . For the remaining components, we have for k = 2, 3, …, N.

Algorithmic improvements

Equations (13) and (14) complete the step of deriving both (λ, v) using only the previously derived eigenpairs (λ, v) for m = 0, 1, …, j – 1, giving us a recursive solution procedure. Once we have determined (λ, v) for j = 0, 1, …, J, we can use these coefficients in (7) and thereby obtain approximations to the the leading eigenpair (λ(μ), v(μ)). Hence we view (13) and (14) as an algorithm for computing the leading eigenpair.

Turning to the numerical implementation, we now describe two improvements to the algorithm given by (13) and (14).

First, we note that in (14), we always set . This implies that for all j ≥ 1. This means that all terms under the summation symbol in (13) vanish, yielding the simplified update formula:

Second, we note that (14) contains a binomial coefficient that becomes prohibitively large to compute for large j. A natural question is whether these large coefficients are compensated by the inverse factors of j! in (7). To quantify this, we define

We then substitute λ=j!λ̂ and v = j!v̂ in (14) and derive

Dividing through by j! and using (16), we derive, again for j ≥ 1,

Applying the same substitutions in (15), we derive

Using (16) in (7), we obtain

Examining the equations in the λ̂ and v̂ variables, we see that all large binomial coefficients and factorials have disappeared. For this reason, in our Octave implementation of the perturbative method, we use the recursive system given by (18) and (17), together with the summation formula (19). Finally, applying (4), our answer for the leading eigenpair of the matrix Q̃(μ) defined by (2) is (λ(μ), w−1/2v(μ)).

Example

Let us give an example of the perturbative method in practice. We set N = 5, μ = 0.01, and w equal to a diagonal matrix whose entries along the diagonal are (ϕ, ϕ, ϕ2, ϕ2, ϕ) where d = 0.2 and ϕ = 0.32768.

Let J denote the total number of iterations we run the perturbative method. Starting from λ̂0 = w1 and v̂0 = e1 as in (6), we iterate using (18) and (17) from j = 1 up to j = J. We then evaluate the solution using (19) truncated at j = J, giving us an approximate solution that we denote (λ(μ), w−1/2v(μ)).

For an N-dimensional vector x = (x1, x2, ..., x), let ||x||∞ = max1≤ |x|, the infinity-norm of x, with respect to which the error of the approximate solution after J iterations is

In Figure 3, we plot (in circles) the log10 of the error as a function of the number of iterations J, and (in solid black) the least-squares line of best fit to the data. From J = 1 to J = 14, the log10 errors show a strongly linear trend, confirmed by the R2 = 0.9958 value for the regression line. The slope of the line is approximately −0.9944, implying

Figure 3

For a particular eigenvalue problem, we plot (in circles) the log10 of the error committed by the the perturbative method after J iterations, where J goes from 1 to 14.

Notes: We also fit and plot (in solid black) a least-squares regression line to the log10 errors; for this line, R2 = 0.9958. The plot shows that the perturbative method converges exponentially to an eigenpair with a final error of 2.2590 × 10−16.

Machine epsilon in Octave is approximately 2.2204 × 10−16, and the error after J = 14 iterations is 2.2590 × 10−16. Therefore, for this particular example, Figure 3 shows that the perturbative method converges exponentially to a solution with error on the order of machine epsilon.

Extension to base/codon/word mutation models

Consider now a matrix μB defined as follows: with constant parameters 0 < μ < 1 and k ≥ 1, and the indexing of both rows and columns corresponding to bases in the ordered set B = (A, G, C, T). This matrix was employed in prior work3 and corresponds to the Kimura two-parameter base mutation model.18 When k = 1, μ represents the Jukes-Cantor mutation model.19

Matrix μB shares with matrix μ the properties of linearity in parameter μ, and reduction to the identity matrix when μ = 0. The methods of this section therefore apply to μB directly.

More biologically realistic CMC models2–4 employ codon mutation models. Let C = B, the p-th Cartesian product of the set B. A codon c ∈ C is a string of bases b1b2...bp of pre-specified length p with bi ∈ B.

The codon mutation models studied by Ardell and Sella assume independence of mutation of bases within codons and that all bases mutate according to the same model of evolution μB. With these assumptions, mutation from any codon c1 ∈ C to another codon c2 ∈ C is represented by a matrix μC that is the p-th Kronecker power of a matrix μB as follows: where (· ⊗ ·) is the Kronecker product, and codons are indexed in both rows and columns in lexicographic order.

If λ is the leading eigenvalue μ, then λ = λ is the leading eigenvalue of μC. Similarly, if v is the eigenvector corresponding to the leading eigenvalue of μB, then is the eigenvector corresponding to the leading eigenvalue of μC.

Therefore, the methods of this section allow calculation of the leading eigenpairs of matrices μC as Kronecker powers of leading eigenpairs of μB.

Homotopy method: quasispecies with multiple most fit genotypes

We now give a second method for finding the leading eigenpair of the matrix Q defined in (3). This method, which we call the homotopy method, is motivated by the desire to handle a fitness matrix w that does not have a unique maximal element along its diagonal. The homotopy method produces accurate approximations of the leading eigenpair for such problems.

Problem formulation

Our goal is still to find the leading eigenpair of Q̃ defined in (2). As before, we will instead focus our attention on the symmetric matrix Q defined by (3). We define the matrix-valued function

Note that F(0) = μ and F(1) = Q. The function F smoothly deforms μ into Q—such a function is often called a homotopy in the mathematical literature.16 Note that for all ∈ ∈ [0,1], the leading eigenpair (λ(∈), v(∈)) of F must satisfy

The basic idea behind the homotopy method is to use F to form a bridge between μ, a matrix whose leading eigenpair we already know, and Q, a matrix whose leading eigenpair we seek.

When ∈ = 0, (21) reduces to μv(0) = λ(0)v(0). By the results provided in supplementary materials, we know that the leading eigenvalue of F(0) = μ is 1 with corresponding eigenvector 1⃗ = (1, 1, ..., 1), the column vector of N ones. This implies that

When ∈ = 1, (21) reduces to Qv(1) = λ(1)v(1). Thus the question is how we can use our knowledge of (λ(0), v(0)) and the function F(∈) to derive (λ(1), v(1)), the leading eigenpair of Q.

Unlike the perturbative solution, at no point will we assume that the entries of w have a unique maximum. To make the derivation easier to read, we define so that F(∈) = μ + ∈P and

Homotopy solution

We differentiate both sides of (21) once with respect to ∈ and obtain

Since F(∈) is symmetric for all ∈, we see that the transposition of (21) can be written v(∈)(∈) = λ(∈) v(∈). Thus, after multiplying (25) through on the left by v(∈), the second term on the left-hand side cancels the second term on the right-hand side, leaving

Now let us substitute (26) back into (25). Solving for v′(∈), we have

We now recognize (26) and (27) as a system of ordinary differential equations (ODEs) with ∈ playing the role of a time-like independent variable:

We also recognize (22) as the initial conditions for this system of ODEs. Let us now describe an elementary algorithm for solving this system:

Set λ = 1 and v = 1⃗. Fix an integer number of steps nsteps and then set Δ∈ = 1/nsteps. Also set ∈ = 0, initially.

While ∈ < 1:

Compute λ′ using (26), i.e., λ′ = (v)/(v).

Set λ ⇆ λ + (Δ∈)λ′.

Using this updated λ, compute v′ using (27), i.e., v′ = −[μ + ∈P − λI]−1Pv.

Set v ← v + (Δ∈)v′.

Set ∈ ← ∈ + Δ∈.

The algorithm will terminate in nsteps steps, yielding an approximation to the leading eigenpair of Q that is stored in λ and v.

To obtain the leading eigenvector of Q̃, we compute ṽ = w−1/2v. The leading eigenpair of Q̃ is then (λ, ṽ).

We now give test results for the homotopy method applied to a particular problem. We set N = 8, μ = 0.1, and w equal to a diagonal matrix whose entries along the diagonal are (a, b, b, b, b, b, b, b) with a = 0.63631836 and b = 0.73306514. We also set the parameter, φ = 1.

We repeatedly run the algorithm given above for different values of nsteps; specifically, we take nsteps = 10, where j = 1, 2, 3, 4, 5, 6. For each value of nsteps, we compute an approximation to the leading eigenpair, which we denote by (λ, v). We then evaluate the residual error of this approximation using

In Figure 4, we plot (in circles) the log10 of the error as a function of the log10 of the number of steps. We also plot (in solid black) the least-squares line of best fit to the data; for this line, R2 = 0.9899 and the slope is approximately −1.0115, implying

Figure 4

For a particular eigenvalue problem, we plot (in circles) the log10 of the error committed by the the homotopy method using a number of steps given by nsteps = 10, where j goes from 1 to 6.

Notes: We also fit and plot (in solid black) a least-squares regression line to the log 10 errors; for this line, R2 = 0.9899. The plot shows that the homotopy method’s error is linear in Δ∈ = 1/nsteps.

Note that with nsteps = 103, the error is approximately 2 × 10−5. This level of error is acceptable if we seek to use the homotopy method to reproduce, for example, Figure 4 in a previously published paper.4 Hence we use this value of nsteps as the default value in the CMCpy implementation of the homotopy method.

Further note that when nsteps = 106, even the elementary algorithm described above to solve the system of nonlinear ODEs (28) is capable of producing a residual error of approximately 2 × 10−8.

For this example, the homotopy method displays convergence that is linear in Δ∈—this relatively slow rate can be improved dramatically by using more sophisticated methods to solve the system of nonlinear ODEs (28), an issue we leave for future work.

In order to apply the homotopy method to models where the mutation matrix is given by μB as defined earlier, there is one requirement: we must be sure that μB has a unique maximal eigenvalue of 1. For the μB matrix, it turns out that we can explicitly derive all eigenvectors and eigenvalues, for general values of both μ and k. With the constraints 0 < μ < 1 and k ≥ 1, the derivation that we give in supplementary materials below proves that the μB matrices have a unique largest eigenvalue of 1. Matrices μB therefore fulfill the minimum requirements for applicability of the methods of this section, after which the leading eigenpair for corresponding matrices μC may be calculated using Kronecker powers as before.

Future outlook

We succeeded in programming a GPU-based power method implementation that is faster than its CPU analogue; however, the performance gain that we obtained with it was not as great as we had hoped. Furthermore, even though our implementation is correct, we could not completely eliminate divergence in evolutionary trajectories in power method implementations arising from differences in machine number representations and precision across platforms. We believe that this arises from deviations in the way double precision floating point numbers are represented, computed on, and rounded in CUDA compute capability 1.3. Perhaps utilization of the cuBLAS library in the future would make performance closer to Numpy and better conform to IEEE standards.

Furthermore, while the power method converges linearly, our new perturbation method provides exponential convergence with the same accuracy. However, this method cannot handle the case of non-unique most fit genotypes that occurs in CMC models at their initialized state. On the other hand, our new homotopy method does handle this case, yet as currently implemented it also converges linearly, rather slower than the power method (results not shown). Incorporation of more sophisticated methods to solve systems of nonlinear ordinary differential equations should dramatically improve performance of our new homotopy method application. We leave further development of both methods and their implementation in CMCpy for future work; perhaps other CPU or GPU implementations of them will compete with Numpy. More generally, our analytical results greatly expand the domain of quasispecies models that can be accurately solved using analytical approaches, particularly multi-site models with biologically realistic mutation parameters.

A variety of open problems remain concerning CMC models and in the field of the evolution of the genetic code.20–24 CMCpy can easily be extended to implement the model studied by Vetsigian et al. (2006)5 with variations, or alternative observables, such as the “evenness” of amino acids.23 Current models of the genetic code have not yet integrated a theory for the origin of translation per se.5,21 We believe that extensions to CMC models will better address such fundamental questions and hope that CMCpy and our analytical solutions to quasispecies models will play a role in that work.