An Evolutionary Theory on Virus Mutation in COVID-19

1. Introduction

The novel coronavirus, SARS-CoV-2, has been continuously mutating for over three years, giving rise to different variants such as alpha, beta, gamma, delta, and omicron. The complex relationship between virus antigenicity, transmission, and virulence makes the future trajectory of COVID-19 unpredictable (Carabelli et al., 2023; Markov et al., 2023; Cao et al., 2023; Mallapaty, 2022; Callaway, 2021; Wang et al., 2023). Recent research has identified two primary forces controlling virus evolution: intrinsic transmissibility, determined by the ACE2 binding affinity of SARS-CoV-2, and immune evasion, achieved through a reduction in susceptibility to neutralizing antibodies (Ma et al., 2023). Studies have also examined the relationship between the spike protein charge and viral infectivity and transmissibility (Božič and Podgornik, 2023; Cotton and Phan, 2023).

In a prior study (Luo and Lv, 2023), we introduced a mathematical model for analyzing the dynamics of COVID-19 spread, with a specific focus on the competitive dissemination of two viruses within a given region. Building upon this work, the present article puts forth an evolutionary theory concerning virus mutation in the context of COVID-19. Our approach involves deducing an evolutionary tree based on sound assumptions regarding mutant representation, alongside the application of a well-established tree construction algorithm. Importantly, we illustrate that the theoretical tree effectively mirrors the observed evolutionary patterns of existing virus strains. Moreover, leveraging a statistical method developed in our study, we extend the utility of the evolutionary tree to predict the emergence of novel macro-lineages of viruses. These findings supply a critical theoretical framework for comprehending the evolution of SARS-CoV-2.

2. Material and Method

2.1. SARS-CoV-2 variants

We obtained a total of 25 variants from the mutation reports provided by outbreak.info (https://outbreak.info/) (Gangavarapu et al., 2023). These variants belong to the SARS-CoV-2 Pango lineages. Our analysis focused solely on mutations that occurred in the spike protein of the virus. Additionally, we only included mutations that were found in at least 75% of the sequenced SARS-CoV-2 lineages. Table 1 provides a comprehensive list of all 25 variants and their corresponding mutation sites.

2.2. ACE2 binding affinity of a single mutation

We obtained the data for single-point mutations affecting the interaction between the receptor-binding domain (RBD) (amino acid residues 331 to 531 on the spike protein) and the ACE2 receptor from reference (Starr et al., 2020). The mean and standard deviation of the affinity for the 19 mutations occurring at the i-th residue are denoted as b_i and s_i, respectively. If the affinity value m_i for the i-th residue is greater than b_i+s_i, the mutation is classified as an affinity-enhancing type. If m_i<b_i-s_i, the mutation is classified as an affinity-weakening type. If b_i-s_i≤m_i≤b_i+s_i, the mutation is considered to have little effect on the binding affinity.

2.3. The n-distance between two sequences

The evolutionary tree provides a detailed description of the relationships between different viral strains. Consistency with the tree of life is a stringent requirement for any proposed evolutionary model, making it a sensitive test for such models. Let p_a represent the probability of a letter "a" (where a can be 0, 1, 2, or 3) occurring in a sequence, and let p_ab represent the joint probability of letters "a" and "b" appearing sequentially in the sequence. In general, let σ=abc... represent an n-letter segment, and let p_σ denote the joint probabilities of the bases in σ occurring in the sequence. In the calculation of joint probabilities, all sequences are assumed to be circular. For any given n, the sum of the joint probabilities over all segments σ of length n is always equal to 1. This can be written as ${\displaystyle \sum _{\sigma }{p}_{\sigma }}=1$, where the summation is over the set of all 4ⁿ segments of length n. Given two sequences Σ and Σ' with sets of joint probabilities {p_σ} and {p'_σ} respectively, we define a distance, called an n-distance, between the two sequences based on the difference in joint probabilities as follows:

$$E_n(\Sigma ,{\Sigma }^{\prime })={\displaystyle \sum _{\sigma }\left|p_{\sigma }-{p^{\prime }}_{\sigma }\right|},n=1,2,\mathrm{...}.$$

(1)

In the following, when there is no confusion, the arguments of E_n will be omitted. An n-distance is well-defined for sequences of different lengths that are not aligned. By repeatedly applying relations such as equation (2) below:

$$\left|p_{\sigma }-{p^{\prime }}_{\sigma }\right|=\left|{\displaystyle \sum _a(p_{\sigma a}-{p^{\prime }}_{\sigma a})}\right|\le {\displaystyle \sum _a\left|p_{\sigma a}-{p^{\prime }}_{\sigma a}\right|},$$

(2)

where σ is any n-letter segment and σa is an (n+1)-letter segment, it can be deduced that:

$$E_{n+1}\ge E_n,n=1,2,\mathrm{...}$$

(3)

2.4. The construction method of the evolutionary tree

To infer the phylogenetic tree, three main algorithms are commonly used in the literature, based on the evolutionary information contained in the sequences: distance matrix treeing, maximum parsimony analysis, and maximum likelihood method (Li, 1997). In the distance matrix method, evolutionary distances are initially computed for all pairs of taxa. Subsequently, a phylogenetic tree is constructed using an algorithm (such as the unweighted pair-group method with arithmetic mean or UPGMA, the neighbor-joining method or NJ, etc) that relies on functional relationships among the distance values. The maximum likelihood method involves studying the relationship between each sequence and the common ancestor in all possible trees. In this article, we utilized the UPGMA (Nei and Kumar, 2000) distance matrix method for constructing the phylogenetic tree. We attempted to use alternative algorithms to reconstruct the evolutionary tree of 25 virus strains but found that the UPGMA tree was superior to others (such as the NJ tree, etc.).

3. Results

3.1. Four-letter representation of the mutants

The surfaces of coronaviruses are decorated with a spike protein, about 1273 amino acids for SARS-COV2. We propose representing any SARS-CoV-2 virus strain by a 1273-long sequence of four letters: 0, 1, 2, and 3. Here, 0 represents no mutation relative to the wild type, 1 represents a mutation having little effect on the ACE2 binding affinity, 2 represents a mutation of affinity-enhancing type, and 3 represents a mutation of affinity-weakening type.

We categorize the old SARS-CoV-2 mutants (alpha, beta, gamma, delta, etc.) as the N-lineage, while the recent mutants (omicron BA.2, BA.5, etc.) are classified as the O-lineage. We consider 12 main mutants in the N-lineage and 13 mutants in the O-lineage. Figure 1 presents the representation of 25 mutants in the receptor-binding domain (RBD) of the spike protein, from site 331 to site 531. For example, the RBD of the Delta mutant is represented by (0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,2,0,0,0,0,0,0,0,0), and the RBD of the BA.1 mutant is represented by (1,0,0,1,1,3,0,0,0,0,0,0,0,0,0,0,0,0,1,2,1,0,0,1,2,1,1,2). The sequences are 27 long, corresponding to 27 mutated sites. When considering all conservative sites denoted by 0 and inserting them in the sequence, we obtain 201-long sequences to describe the RBD. The 27-long sequence is the reduced sequence of the 201-long sequence.

The above representation for the RBD domain (201 sites) can be generalized to the full spike protein (1273 sites). Outside the RBD (sites 1 to 330 and sites 532 to 1273), only two letters, 0 and 1, occur since ACE2 binding does not occur in these domains. Thus, the full-spike representation of a SARS-CoV-2 mutant is a 1273-long sequence of four letters. Among the twenty-five known mutants, there are 104 mutated sites. Therefore, the reduced sequence of the full spike protein is a 104-long sequence. Table 1 lists the 25 selected mutants and their 104 mutated sites on the spike protein.

3.2. Evolutionary tree of virus mutation deduced by UPGMA method

Using the n-distance, we can construct a distance matrix D for a set of mutants labeled by i, j, k, ..., where the matrix element D_ij is equal to the n-distance between Σ_i and Σ_j, with Σ_i being the sequence representing mutant i. Then, by applying the UPGMA algorithm, we reconstruct the evolutionary tree of 25 mutants for a given value of n. We observed that the structure of the tree changes with increasing n and becomes stable after a few steps. Partial results are shown in Figure 2. Specifically, when n≥5, the tree topology of the branch corresponding to macro-lineage O achieves stability. However, the structure of branch N slightly adjusts with increasing n. When n≥7, both branches O and N achieve stability as given in Figure 2C.

Figure 2C presents a UPGMA tree of 25 virus strains, displaying two distinct large branches corresponding to the macro-lineages N and O of existing virus strains. Notably, in macro-lineage O, the sub-branch accurately reflects the bifurcation of the earliest BA.1 strain and the subsequent BA.2, BA.4, BA.5 strains, as well as the correct classification of XBB as a recombination of two BA.2 strains. In macro-lineage N, the sub-branch of the tree shows four VOCs (variants of concern) - alpha, beta, delta, and gamma, belonging to independent sub-lineages. These observations highlight the highly consistent nature of the theoretical tree with all the evolutionary peculiarities of existing virus strains. This work presents a logically simple theory of virus mutation, whereby the tree is deduced from a single basic assumption - the four-letter representation of mutants and the n-distance algorithm for tree reconstruction. It is surprising that the theory aligns so well with experimental data.

The four-letter representation involves choosing between 0 and 1 for 104 mutational sites on the spike protein, while allowing for additional choices (i.e., choices 2 and 3 in addition to 0 and 1) for 27 sites on the RBD. The first part provides 104 bits of information, while the second part contributes only log(27×19)/log2=9 bits of information. Therefore, to simplify calculations, one can approximate the four-letter representation as a two-letter representation that neglects choices 2 and 3. The validity of this approximation can be seen from the fact that the tree topology in Figure 2 remains unchanged when the two-letter approximation is used.

3.3. Generation of new strain on the tree by stochastic method

Once the evolutionary tree is reconstructed based on known lineages of virus mutants, the generation of new strains on the tree can be examined. In this section, we propose a stochastic approach to address this problem. Let A be the set of all mutated sites in the updated virus strains, with the symbol a representing the number of sites. Each strain is represented by a two-letter sequence from the set A(a). Considering stochastic mutation on the spike protein, we define the set X (including x sites) as the set of sites participating in the stochastic sampling. The union of sets A(a) and X(x) is denoted as Z, with Z=A∪X. The intersection of sets A and X is denoted as Y, which contains y sites and is represented as Y(y). Assuming that the new strain is generated by stochastic sampling in set Z, we discuss the new strain represented by a two-letter sequence from the union set Z(len), where len=a+x-y. Since the UPGMA tree provides a good description of the evolutionary peculiarities of virus mutation, we use the UPGMA method (for n=7) to reconstruct the evolutionary tree of the 25+1 mutants, where the additional strain is the predicted new strain. This approach is referred to as the A-X model.

In the A-X model, the new strain is represented by a two-letter sequence in the union set Z. Figure 3A illustrates the schematic representation of the union set Z=A∪X and its relationship with sets A, X, and Y. Figure 3B presents the coding rules for the 25 mutants and 1 new strain.

Now we aim to demonstrate the feasibility of the A-X model in predicting new strains of SARS-CoV-2.

By performing stochastic sampling M=10⁵ times, we obtained 10⁵ predictions of new strains. Interestingly, when x<17, the predicted new strain always falls within one of the two major branches, N and O, as shown in Figure 2. This indicates that the new strain belongs to either macro-lineage N or O. However, when x>17, an anomaly occurs, suggesting that the new strain may occur outside the branches N and O. In other words, the new strain does not belong to either lineage N or O for large values of M. We designate this new macro-lineage as P. Based on a typical stochastic sampling the generation of the new strain within the tree is illustrated in Figure 4. For a given pair x and y the predicted new strain is represented as Newx.y. In Figure 4A, x=16 and y=4, the new strain is located within macro-lineage N; in Figure 4B, x=16 and y=1, the new strain is located within macro-lineage O; in Figure 4C, x=35 and y=2, the new strain is located within macro-lineage O; and in Figure 4D, x=35 and y=3, the new strain is located within macro-lineage P.

3.4. The stochastic property of the intersection set

Since M=100,000 corresponds to 16.6 bits of information, which is only 1.3% of the 1273 bits for the entire spike protein, we employ a supplementary method to expand the range of stochastic sampling by randomizing the parameter y. This means that both the intersection set Y(y) and the union set Z(a+x-y) are altered. In Figure 4, we discovered that the change in the parameter y for a given x significantly affects the tree topology. Therefore, randomizing the parameter y is expected to be an efficient approach to enhance the stochastic information. To gain further insights into y-randomization, we performed six randomizations and treated y as a fluctuating quantity in each group of given x, calculating the statistical distribution of y. The resulting distributions are presented in Table 2.

Upon examining Table 2, we observe that the mean value of y remains constant within each group, specifically y=(a/1273)x, where a represents the number of sites in set A(a=104). This constancy arises from the fact that the probability of a stochastic-produced site occurring in set A is a/1273, and the total occurrence is given by (a/1273)x for x sites. Additionally, the skew value remains approximately constant for a given x and decreases as x increases. Consequently, the parameter y tends towards a normal distribution as x becomes larger. It is worth noting that extending the randomization of y from six times to a greater number of times does not alter the aforementioned conclusion. Thus, randomizing y for a given x to expand the range of stochastic sampling is a viable approach. This technique effectively increases stochastic information and demonstrates the reliability of the stochastic method in the A-X model.

3.5. The prediction of macro-lineage of SARS-CoV-2

By employing the A-X model, we can make predictions regarding the emergence of a new macro-lineage P, along with the already known macro-lineages N and O, when the number of sites (x) in X is sufficiently large. In Figure 4, we have illustrated how this new macro-lineage can arise through stochastic sampling under the condition of large x. To further investigate this phenomenon, let us consider the following example: a=104 for set A, x=77 and y=6 for set X. The resulting UPGMA tree is presented in Figure 5. Our analysis reveals that the newly identified strain, New77.6, is not classified under macro-lineages N or O, but instead falls into a distinct macro-lineage, P. Furthermore, by conducting 10⁵ rounds of stochastic sampling for x and 6 rounds of randomization for y, we consistently observe the inevitable generation of the novel macro-lineage P when x is greater than or equal to 77. Additionally, the tree structure exhibits similarities to that depicted in Figure 5 across various samplings and randomizations for such large values of x.

To delve into the formation of the new macro-lineage P, let us examine the generation of 10⁵ new variants for each specific x value. The probability of the N, O, or P-lineage at a given x is determined by dividing its occurrence count by the total number of variants, which is 10⁵. The findings of this calculation are succinctly presented in Figure 6, where the probability of the j-lineage (j=N, O, or P) is plotted against x. From the diagram, we can observe that the N-lineage prevails when x is below approximately 22, the O-lineage dominates between roughly 22 and 56, and the P-lineage takes over when x exceeds about 56. Hence, when the number of mutated sites is less than 22 (as x represents the count of stochastic-mutational sites per variant), there is a high likelihood that the variant belongs to macro-lineage N. Similarly, if the number of mutated sites falls within the range of approximately 22 to 56, the variant is likely part of macro-lineage O. However, if the number of mutated sites surpasses 56, the variant probably belongs to macro-lineage P. Currently, the maximum number of mutated sites for a single variant is 44 (as EG.5.1), indicating that the newly generated strain still falls within macro-lineage O. The emergence of the P-lineage will probably occur only after a variant with a greater number of mutated sites appears.

In our A-X model, the selection of set X is stochastic, resulting in the stochastic occurrence of new variants on the tree. Nonetheless, we have demonstrated that the generation of multiple macro-lineages and the subsequent tree topology bifurcation into several major branches are both certain and unavoidable. Furthermore, we have proved that the assignment of a macro-lineage to a new variant is statistically dependent on x, and there are specific demarcation values of x that differentiate between different macro-lineages. By performing 10⁵ iterations of stochastic production of the x-long sequence and randomizing the parameter y six times, we discovered that the demarcation values of x are as follows:

Newx.y∈N when x≤13; Newx.y∈N or O when 13<x≤ 16; Newx.y∈N or O or P when 16<x≤ 39; Newx.y∈O or P when 39<x≤ 76; and Newx.y∈P when x>76.

The demarcation values of x align with Figure 6 and provide a more detailed classification of the occurrence of the three macro-lineages.

4. Discussions

4.1. On mutational robustness

A genotype is represented by a four-letter (two-letter) sequence. The evolutionary phenotype refers to the evolution of a sequence that exhibits two or more macro-lineages. A group of mutants is assigned to a specific branch on the tree, such as the alpha-beta-gamma-delta branch (macro-lineage N) or the omicron branch (macro-lineage O). The clear distinction between the two branches of the tree signifies the existence of mutational robustness and the occurrence of the evolutionary phenotype. The occurrence of the evolutionary phenotype depends on the effective organization of different genotypes. What is the principle behind the organization of genotypes? We observed that the number of mutated sites in the sequence is a crucial factor. Mutants in the N branch exhibit 34 mutated sites, whereas those in the O branch have more than 80 sites. The evolution from an earlier branch to a later branch occurs through an increase in the number of mutated sites. Thus, the evolutionary theoretical study provides evidence for the presence of mutational robustness in virus evolution. Furthermore, our calculations indicate that when the number of mutated sites is sufficiently large, there is a nonzero probability of generating a new macro-lineage (macro-branch P on the tree). The concept of mutational robustness, which maps genotypes to a specific phenotype, was proposed in (Mohanty et al., 2023), where the organization of genotypes was also discussed. This study provides an additional example of mutational robustness in virus evolution.

4.2. Electric charge dependence of virus mutation

Apart from immune evasion, it is generally believed that the primary force driving virus mutation is intrinsic transmissibility, which is mainly determined by the ACE2 binding affinity. However, the ACE2 binding affinity only affects the receptor-binding domain (RBD). Are there any other physical forces that influence the entire spike protein domain? The literature has highlighted the significant role of electrostatic interaction in biological systems (Zhou and Pang, 2018). The population of conformational states can be determined by examining the free energy change during conformational transitions of the spike protein. Electrostatic interaction plays a crucial role in determining this free energy change (Luo and Lv, 2023). As commonly known, the twenty amino acids (AAs) can be classified into three types: D, E, and Y are acidic residues (negatively charged), H, K, and R are basic residues (positively charged), and the remaining amino acids are classified as neutral residues. The viewpoint presented above regarding the relationship between conformational change and electrostatic interaction suggests that variations in amino acid charge may potentially trigger mutations in the spike structure. Therefore, changes in amino acid charge represent structural potential in virus mutations, similar to the ACE2 binding affinity but affecting the entire spike protein. To evaluate the influence of electric charge on virus mutation, we calculated and presented the results in Table 3 (only two mutants for omicron are listed since the cases for other mutants in macro-lineage O are similar).

Analysis of Table 3 reveals that, for the majority of mutants, there is a tendency towards a positive charge in the alteration of amino acid charge during mutation. Previous works (Cotton and Phan, 2023; Pawłowski, 2021) indicate that positively charged amino acid types are generally observed during mutations. However, our calculations indicate that this trend is more pronounced in macro-lineages O. In macro-lineages N, the trend is less apparent, and there are counterexamples. The fundamental principle of virus mutation is the presence of electric charge dependence, which exhibits distinct characteristics in different lineages and macro-lineages.

4.3. Conclusions drawn from the current study

1. We propose representing the SARS-CoV-2 virus strain with a 4-letter sequence, based on amino acid mutations within the spike protein.

2. The application of an n-distance algorithm in UPGMA effectively yields a variant phylogenetic tree that closely aligns with experimental data on virus evolution.

3. We introduce the A-X model for generating new strains on the phylogenetic tree. By combining set A (comprising existing mutated sites) and set X (including x randomly generated sites), detailed studies on the emergence of new strains can be conducted.

4. Expanding the scope of stochastic sampling can be achieved through two means: increasing the size of x and randomizing the parameter y within the intersection set Y=A∩X.

5. When x reaches a sufficient scale, a new macro-lineage of SARS-CoV-2 is predicted, alongside existing lineages. We have derived demarcation values of x that differentiate between various macro-lineages.

6. Therefore, we have proposed a logically simple theory of virus mutation to enhance our understanding of the evolution of SARS-CoV-2.