{'fraction_non_alphanumeric': 0.03925073282333595, 'fraction_numerical': 0.015192678916136413, 'mean_word_length': 4.322488584474886, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 redundancy 2 and non-random patterns of amino acid assignments 3,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 generally 7 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][11][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 theory [13][14][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 codebase 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 easyto-use and program, powerful, and efficient, highlevel scripting language like Python.
Here we present CMCpy, a free open-source codebase 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.
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.
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 models 1-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.
A variety of observables and statistics are implemented including the Normalized Encoded Range 2 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: 1. 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. 2. For every possible site-type in proteins there corresponds a uniquely most fit amino acid. 3. The machineries to decode codons, and to associate any codon to any amino acid, pre-exist the evolution of the genetic code.

7.
Messages are haploid and asexually reproducing. 8. 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 w 1 is the unique maximum of the finite set {w 1 , . . . , w N }. We also assume w i . 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(µ) T = 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 j 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 theory 14 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 = w 1 and v 0 = e 1 . 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 λ j , µ j . 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 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 .
Note that e we This shows that if we already know (λ 0 , v 0 ), we can determine λ 1 . To determine v 1 , 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 We have completed the loop, showing how to proceed from the zeroth-order eigenpair (λ 0 , v 0 ) to the first-order eigenpair (λ 1 , v 1 ).
In the next section, we show how to iterate this procedure to generate the j-th order eigenpair (λ j , v j ) 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: Now on the right-hand side, we peel off the m = j term-note that this the only term in which λ j appears. Hence We use (6) and e e T 1 1 1 = , and we solve for λ j : We have therefore shown that if we already know (λ m , v m ) for m = 0, 1, …, j -1, we can solve for λ j . Now, using this λ j , we can solve for v j in the same way as before. We go back to (12), isolate all terms involving v j , 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: for k = 2, 3, …, N.
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 v j This implies that e T j 1 0 v = 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 !
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 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 (φ, φ d , φ 2d , φ 2d , φ d ) where d = 0.2 and φ = 0.32768. Let J denote the total number of iterations we run the perturbative method. Starting from 0 λ = w 1 and 0 v = e 1 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 (λ J (µ), w -1/2 v J (µ)).
For an N-dimensional vector x = (x 1 , x 2 , ..., x N ), let ||x|| ∞ = max 1#i#N |x i |, the infinity-norm of x, with respect to which the error of the approximate solution after J iterations is In 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 work 3 and corresponds to the Kimura two-parameter base mutation model. 18 When k = 1, µ B 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 models 2-4 employ codon mutation models. Let C = B p , the p-th Cartesian product of the set B. A codon c ∈ C is a string of bases 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 c 1 ∈ C to another codon c 2 ∈ 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 µ B , then λ C = λ p is the leading eigenvalue of µ C . Similarly, if v is the eigenvector corresponding to the leading eigenvalue of µ B , then v v C i p = ⊗ =1 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 F (∈) = µ + ∈(w 1/2 µw 1/2 − µ). (20) 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.
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 P = w 1/2 µw 1/2 − µ, (23) 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(∈) T F(∈) = λ(∈) v(∈) T . Thus, after multiplying (25) through on the left by v(∈) T , 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 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: Fix an integer number of steps n steps and then set ∆∈ = 1/n steps . Also set ∈ = 0, initially. The algorithm will terminate in n steps 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

Example
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) T 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 n steps ; specifically, we take n steps = 10 j , where j = 1, 2, 3, 4, 5, 6. For each value of n steps , we compute an approximation to the leading eigenpair, which we denote by (λ j , v j ). We then evaluate the residual error of this approximation using In Figure 4, we plot (in circles) the log 10 of the error as a function of the log 10 of the number of steps. We also plot (in solid black) the least-squares line of best fit to the data; for this line, R 2 = 0.9899 and the slope is approximately -1.0115, implying

∆ ∈.
Note that with n steps = 10 3 , 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 n steps as the default value in the CMCpy implementation of the homotopy method.
Further note that when n steps = 10 6 , 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  Figure 4. For a particular eigenvalue problem, we plot (in circles) the log 10 of the error committed by the the homotopy method using a number of steps given by n steps = 10 j , 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, R 2 = 0.9899. The plot shows that the homotopy method's error is linear in ∆∈ = 1/n steps . sophisticated methods to solve the system of nonlinear ODEs (28), an issue we leave for future work.

Extension to base/codon/word mutation models
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][21][22][23][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.

competing Interests
Author(s) disclose no potential conflicts of interest.

Disclosures and ethics
As a requirement of publication author(s) have provided to the publisher signed confirmation of compliance with legal and ethical obligations including but not limited to the following: authorship and contributorship, conflicts of interest, privacy and confidentiality and (where applicable) protection of human and animal research subjects. The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material. Any disclosures are made in this section. The external blind peer reviewers report no conflicts of interest.

supplementary Materials Supplementary 1: derivation of complete eigendecomposition of circulant µ matrices
The matrix µ is circulant, an important property because the eigenvalues and eigenvectors of any circulant matrix can be written in closed form. Let i = −1 and define for j = 0, 1, …, N -1 the roots of unity Then the N eigenvalues of µ are, for j = 0, 1, …, N -1, given by Let us now show that the maximum eigenvalue is always equal to 1 and that this eigenvalue is unique.
Consider the function , then x = 0 will still be the unique global maximum of f(x). Therefore, the unique maximum eigenvalue s j occurs at j = 0 and is always equal to s 0 = 1.

Eigenvectors
Let v * denote v -T , the conjugate transpose of the vector v. The eigenvector of µ that corresponds to the eigenvalue s j is written here as a column vector. The 1/ N factor is included to normalize the eigenvalue so that y y j j * . = 1 Note that the eigenvector corresponding to the eigenvalue s 0 = 1 is particularly simple: Let us now show that (30), as defined, is an orthonormal system of eigenvectors. We take the Hermitian inner product of an eigenvector y j with an eigenvector y k , assuming j ≠ k: It is clear that each row of µ B sums to 1. This, by itself, implies that v 1 = (1, 1, 1, 1) T is an eigenvector of µ B with corresponding eigenvalue s 1 = 1.
In what follows, we use the fact that The four column vectors that appear on the righthand side are the four columns of µ B .
It is clear that if we subtract the fourth column of µ B from the third column of µ B , then the first two entries cancel, while the last two entries are exact opposites of one another: Hence we have shown that v 2 = (0, 0, 1, -1) T is an eigenvector of µ B with corresponding eigenvalue In a very similar way, we can see that Note that s 2 = s 3 . Since v 2 and v 3 are clearly orthogonal, we see that {v 2 , v 3 } span a two-dimensional eigenspace, and that the geometric multiplicity of s 2 equals its algebraic multiplicity. This is consistent with the fact that µ B is, by definition, a real symmetric matrix and therefore, by the spectral theorem, must possess four real eigenvalues together with an orthonormal basis of eigenvectors.
As we have found three orthogonal eigenvectors {v 1 , v 2 , v 3 }, we appeal to orthogonality to find v 4 . Let v 4 = (1, x, y, z) T . Orthogonality with v 2 and v 3 respectively imply 1 -x = 0 and y -z = 0, so v 4 = (1, 1, y, y) T . Now, orthogonality with v 1 implies that 2 + 2y = 0, so y = -1, and we have v 4 = (1, 1, −1, -1) T . Multiplying µ B by v 4 , we see that Normalizing the eigenvectors, we obtain the set where (y j )*y j = 1 for j = 1, 2, 3, 4, and each eigenvector y j corresponds to the eigenvalues s j given above. Indeed, examining the expressions of the eigenvalues s j , and using the fact that 0 , µ , 1 and k $ 1, it is clear that s 1 = 1 is the unique maximal eigenvalue of µ B .
Supplementary 3: Benchmark script in R to generate Figure 2.
Supplementary 4: Matlab/Octave implementations of the perturbative and homotopy methods presented in this manuscript.