key: cord-0170119-8s8gezpz
authors: Chen, Jiahui; Wang, Rui; Wang, Menglun; Wei, Guo-Wei
title: Mutations strengthened SARS-CoV-2 infectivity
date: 2020-05-27
journal: nan
DOI: nan
sha: c279abafa0a8cae880e9c9fffd98f48a354b44fb
doc_id: 170119
cord_uid: 8s8gezpz

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infectivity is a major concern in coronavirus disease 2019 (COVID-19) prevention and economic reopening. However, rigorous determination of SARS-COV-2 infectivity is essentially impossible owing to its continuous evolution with over 13752 single nucleotide polymorphisms (SNP) variants in six different subtypes. We develop an advanced machine learning algorithm based on the algebraic topology to quantitatively evaluate the binding affinity changes of SARS-CoV-2 spike glycoprotein (S protein) and host angiotensin-converting enzyme 2 (ACE2) receptor following the mutations. Based on mutation-induced binding affinity changes, we reveal that five out of six SARS-CoV-2 subtypes have become either moderately or slightly more infectious, while one subtype has weakened its infectivity. We find that SARS-CoV-2 is slightly more infectious than SARS-CoV according to computed S protein-ACE2 binding affinity changes. Based on a systematic evaluation of all possible 3686 future mutations on the S protein receptor-binding domain (RBD), we show that most likely future mutations will make SARS-CoV-2 more infectious. Combining sequence alignment, probability analysis, and binding affinity calculation, we predict that a few residues on the receptor-binding motif (RBM), i.e., 452, 489, 500, 501, and 505, have very high chances to mutate into significantly more infectious COVID-19 strains.

In December 2019, an outbreak of pneumonia due to coronavirus disease 2019 (COVID- 19) was initially detected in Wuhan, China [8] , due to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It has now spread globally via travelers and breached the boundaries of 213 countries and territories, leading to more than 5.4 million infection cases and 343,000 deaths as of May 23, 2020. In the past two decades, there have been three major zoonotic disease outbreaks of betacoronaviruses: SARS-CoV in 2002, Middle East respiratory syndrome coronavirus (MERS-CoV) in 2012, and SARS-CoV-2 in 2019. Similar to SARS-CoV and MERS-CoV, SARS-CoV-2 infections were observed in hospital personnel and family clusters in the early stages of the outbreak [2, 5, 16] . Unfortunately, there is no specific antivirus drugs nor effective vaccines developed to moderate this outbreak at present. SARS-CoV-2 is an enveloped non-segmented positive-sense RNA virus and belongs to the betacoronavirus genus. Although the origin of SARS-CoV-2 remains elusive, it has undergone thousands of recorded single mutations compared to the reference genome collected on January 5, 2020 [29, 31] . Because of the lack of proofreading ability in RNA polymerases, RNA viruses are prone to random mutations. Human immune system intervention introduces viral mutations too. However, rapid global spread and transmission of COVID-19 provide the virus with substantial opportunities for the natural selection of rare-acted but favorable mutations. Therefore, although some viral mutations are benign, many mutations strengthen viral survival capability. It is of paramount importance to understand SARS-CoV-2 infectivity changes following the existing mutations and predict the future infection tendency.

It is well known that like SARS-CoV, SARS-CoV-2 enters host cells through the interaction of spike glycoprotein (S protein) and host angiotensin-converting enzyme 2 (ACE2) receptor [6, 14, 27] . In both SARS-CoV and SARS-CoV-2, the S protein receptor-binding domain (RBD) is recognized as one of two subunits, i.e., S1 and S2, to bind directly to the ACE2. Although the SARS-CoV-2 S protein harbors a furin cleavage site at the boundary between the S1/S2 subunits [31] , lessons learned from SARS-CoV are important in formulating hypotheses about SARS-CoV-2, as well as the receptor recognition when studying the host range, cross-species transmission, and pathogenesis of SARS-CoV-2. In the studies of SARS-CoV, epidemiologic and biochemical studies show that the infectivity of different SARS-CoV strains in host cells is proportional to the binding affinity between the RBD of each strain and the ACE2 expressed by the host cell [6, 15, 19, 26, 27] . Therefore, the assessment of binding affinity changes following mutations is vital for the understanding of SARS-CoV-2 infectivity evolution.

It is very challenging to rigorously measure the relative viral infectivity of two viruses by experiments. There is a discrepancy in the literature about the relative S protein-ACE2 binding affinities of SARS-CoV and SARS-CoV-2. Wrapp et al and Shang et al reported that SARS-CoV-2 has a higher binding affinity than SARS-CoV does [22, 30] , whereas Walls et al argued that SARS-CoV-2 and SARS-CoV bind with similar affinities to ACE2. [27] The first SARS-CoV-2 genome reported on January 5, 2020 [31] has about 80% sequence identity with that of SARS-CoV. However, compared with SARS-CoV, SARS-CoV-2 S protein has 301 mutations over its 1255 residues. Their sequence identity is only 76%. Among 301 mutations on SARS-CoV-2 S protein, 50 were on the RBD, which has a total of 194 residues, suggesting that the RBD is subject to more mutations. Our recent studies using over 13,000 genome samples show that SARS-CoV-2 S protein is among the most non-conservative ones in its genome [29] . Since early January 2020, hundreds of new mutations were found on different residue positions of SARS-CoV-2 S protein. Many of them are located on the RBD [29] . The existence of so many different S protein mutations indicates that there are many different SARS-CoV-2 subtypes that might have very different infectivities. Obviously, the relatively high mutation rate at the RBD poses a real threat to the occurrence of future SARS-CoV-2 strains that might be more infectious than the current SARS-CoV-2.

Due to the continuous evolution of SARS-CoV-2, the experimental measurement of SARS-CoV-2 infectivity is extremely difficult if it is not entirely impossible. The computational estimation of mutation-induced protein-protein binding affinity changes is an important approach for understanding the impact of mutations on protein-protein interactions (PPIs). There are many standard databases available, including the AB-Bind database of mutation-induced antibody-antigen complex binding free energy changes [25] and SKEMPI for protein-protein binding affinity changes upon mutation (∆∆G) [17] . These databases have been used as a benchmark for evaluating the predictive power of various computational methods [18, 21] . To simplify the structural complexity of protein-protein complexes, we have recently introduced element-specific and site-specific persistent homology, a new branch of algebraic topology, to embed molecular mechanisms into topological invariants [28] . This approach is paired with a new deep learning algorithm called NetTree, to combined convolutional neural networks and gradient boosting trees. The resulting method, called Top-NetTree, was about 22% better than the previous best result for the AB-Bind dataset and significantly outperformed the state-of-the-art in the literature on the SKEMPI database [28] .

The objective of this work is three-fold. First, we train the TopNetTree on a set of 8338 PPI data to analyze the impacts of existing S protein RBD mutations on the binding affinity of the S protein and the ACE2. Since different SARS-CoV-2 subtypes have different mutation patterns, it is important to understand their mutation impacts accordingly. We carry our analysis based on existing six mutation clusters [29] , though more specific analysis can be easily done as well. Additionally, it is also extremely important to know whether future SARS-CoV-2 subtypes would pose an imminent danger to public health. To this end, we have conducted a systematic screening of all possible 3,686 future mutations on all 194 residues (residue IDs from 333 to 526 on S protein) on the RBD. We classify these mutations into three categories: the most likely ones which would happen by a single mutation at any one of three constitutive nucleotides; the likely mutations which would occur via two concurrent mutations at three constitutive nucleotides; and the unlikely mutations which would produce through three concurrent mutations at all of three constitutive nucleotides. Finally, we analyze how 50 mutations on the RBD of the SARS-CoV-2 S protein with respect to SARS-CoV have changed its infectivity. Figure 1 : The scatter plot of six distinct clusters in the world. The light blue, dark blue, green, red, pink, and yellow represent Cluster I, Cluster II, Cluster III, Cluster IV, Cluster V, and Cluster VI, respectively. The size of the pie charts corresponds to the number of samples. The color of the dominated cluster decides the base color of each country. N TF  US  146 829 450 2948 126 1171 1600 11963 946 6766 793 8539  CA  17  95  16  79  14  119  23  154  40  240  8  76  AU  99  555 357 3533 118 1020  115  827  69  469 122 1423  UK  852 5073 908 6031 1545 14549 124  961  3  15  544 5175  IS  145 868  89  474  89  870  70  472  15  127  17  145  ES  111 677  84  555  25  217  7  59  2  6  35  307  CN  2  8  192 910  1  13  1  7  25  58  8  69  DE  20  100  20  97  38  324  41  274  0  0  10  67  FR  64  385  14  55  12  105  98  49  0  0  46  446  IN  33  184 114 1145  9  84  5  40  0  0  69  753  RU  18  95  1  2  68  583  3  21  0  0  7  63  BE  111 647  39  151  91  815  16  114  0  0  68  599  SA  1  5  9  61  1  7  16  110  0  0  30  254  IT  18  93  5  110  17 To investigate the influences of existing S protein RBD mutations on binding affinity (BA) of S protein and ACE2, the 13752 complete SARS-CoV-2 genome samples deposited at GISAID [23] are compared with the first genome sequence of SARS-CoV-2 collected on January 5, 2020 [31] . The resulting 7525 single mutations are found in six distinct clusters as shown in Table 1 [29] and Figure 1 . There are 434 existing non-degenerated mutations on SARS-CoV-2 S protein. Among them, 55 mutations occurred on the RBD which are relevant to the binding of SARS-CoV-2 S protein and ACE2. Furthermore, 28 out of 55 mutations are on the receptor-binding motif (RBM), i.e., the region of RBD that is in direct contact with the ACE2.

We examine the free energy changes following the existing site-specific mutations. Our studies are based on the X-ray crystal structure of SARS-CoV-2 S protein and ACE2 (PDB 6M0J) [11] (see Fig. 19 ), whose S protein gene sequence is consistent with that of the reference SARS-CoV-2 [31] . The BA change following mutation (∆∆G) is defined as the subtraction of the BA of the mutant from the BA of wild type, ∆∆G = ∆G W − ∆G M where ∆G W is BA of the wild type and ∆G M is BA of mutant type. Therefore, a positive BA change means that the mutation increases affinities, making the mutant more stable and more infective.

We present the overall BA changes ∆∆G of SARS-CoV-2 S protein RBD in Figure 2 . Most mutations have a small number of changes in their binding affinities, while some of them have large changes. There are 56% mutations on RBD having positive BA changes (i.e., 31 over 55) including V367F and V483A, which have the most frequencies, 13 and 23 respectively. This statistic implies that the evolution of SARS-CoV-2 is mostly driven by selection and COVID-19 evolves toward more infectious. It is noted that many mutations on the RBM, such as N439K, L452R, and T478I, have significant free energy changes. The mutations on the RBM take 51% (28 over 55) of all mutations on the RBD, which potentially increases the complexity of antiviral drug and vaccine development. This global analysis indicates that mutations on RBD strengthen the binding of S protein and ACE2, leading to more infectious SARS-CoV-2.

The SARS-CoV-2 genotypes are clustered into six clusters or subtypes based on their single nucleotide polymorphism (SNP) variants [29] . Accordingly, a more detailed analysis of mutation impacts on the BA changes can be carried out on each cluster, which reveals the diversity of COVID-19 infection rates and provides evidence for transmission pathways and spread dynamics across the world.

It is worth noting that residue 414 has three mutations, Q414P, Q414E, and Q414R, due to mutations at two adjacent nucleotides 22802 and 22803: 22803A>C, Q414P; 22802C>G, Q414E; and 22803A>G, Q414R. At the protein level, some or all of these mutations show up in different clusters. Similarly, residues 354 and 521 have two existing mutations. Figure 3 depicts the binding affinity changes ∆∆G of Cluster I. Four out of sixth mutations have positive binding affinity changes which indicate increasing infectivity on Cluster 1. Particularly, mutation N439K, which has a higher frequency and a larger free energy amplitude than the rest, attained a very positive binding affinity change to increase the affinity between S protein and ACE2. Therefore, Cluster I has a moderate increase in its infectivity. Cluster I is associated with COVID-19 in most countries except for Japan and South Korea [29] . Figure 4 illustrates the binding affinity changes following the mutations of Cluster II. As shown in the figure, there are many mutations on the RBD. However, most mutations are associated with small free energy changes. When only considering the absolute value of BA change greater than 0.5 kcal/mol, five mutations have positive BA changes whereas two mutations have negative BA changes. The mutation D364Y has the highest frequency and the highest free energy change, indicating the increase in infectivity. Therefore, Cluster II has a minor increase in infectivity. Note that Cluster II COVID-19 is found in every country that has submitted SARS-CoV-2 genome samples. 

From Figure 5 , a significant decreasing trend of affinity is observed such that the largest change is also the most frequent mutation with a negative BA change while the rest of the mutation-induced BA changes are negligible. Interestingly, the mutation T478I changes from amino acid with polar uncharged side chains, Threonine, to amino acids with hydrophobic side chains, Isoleucine, which significantly decreases affinity between the S protein and the ACE2 receptor. Another observation is that mutations that happened in the same residues have close changes such as P384L and P384S, S477R and S477N, or Q414P and Q414R.

Cluster III is only one which has a decreasing trend of binding affinities as SARS-CoV-2 evolving from this cluster. This cluster involves genome samples from all countries except for South Korea. Notably, most samples are submitted by the UK. Figure 6 shows the binding affinity changes of mutations in Cluster IV. Among all mutations in Cluster 4, mutation L452R has the largest free energy change and directly connects the ACE2 receptor. Although the most frequent mutation Q414E has a negative change, the BA change ∆∆G is -0.055 kcal/mol which is negligible compared with others. The overall trend of this cluster is considered as increasing the COVID-19 infectibility. Note that most genome samples in Cluster IV are submitted by the US. 

The binding affinity changes in the last cluster is shown in Figure 8 . Obviously, most mutations on Cluster VI have enhanced the binding affinity of the S protein and ACE2 receptor except T478I. The most significant positive free energy change is caused by mutation K417N. Overall, Cluster VI has strengthened infectivity. This cluster involves genome samples submitted from all countries except for Japan and South Korea. The US submitted most samples. Cluster VI is a new cluster of SARS-CoV-2 according to its low frequency of each mutation.

In summary, five of six clusters, Clusters I, II, IV, V, and VI, have moderate or minor positive binding affinity changes which indicate that the evolution of the SARS-CoV-2 trend to increase infectability by increasing its binding affinity with ACE2 receptor. One cluster, Cluster III, reduced its binding affinity, indicating a decrease in its infectivity.

In this section and the next section, we analyze the impacts of all of 3686 possible mutations on 194 S protein RBD of 194 residues. On each amino acid, we classify all 19 possible mutations into most likely future mutations, likely future mutations, and unlikely future mutations. Here, most likely, likely, and unlikely future mutations are defined by the protein mutations induced by only one, simultaneous two, and simultaneous three of genetic changes on three underlying nucleotides on a codon. Based on the codon analysis of all 194 amino acid residues on the RBD, we have 1149 most likely, 1912 likely, and 625 unlikely mutations.

We compute the ∆∆Gs following most likely future mutations on the RBD. Figure 9 depicts 20 most likely future mutations that can have the highest adversarial impacts on COVID-19 infectivity. First, it is noted that mutation Y495N on the RBM has the highest free energy change and if it occurs, it will make the virus significantly more infectious. Additionally, mutation Y489H on the RBM would incur another large infectivity strengthening. It is worthy to note that residue 489 is a potentially hot spot, where 5 possible mutations, Y489H, Y489D, Y489F, Y489C, and Y489N, will lead to the strengthened S protein-ACE2 binding. The other potentially hot spot is residue 423 with Y423C, Y423F, and Y423S being infectivitystrengthening mutations. Residue 452 on the RBM has been proven to hot spot as it has an existing mutation L452R (see Fig. 6 ) and another infectivity strengthening mutation, L452P. In general, the highest free energy changes are due to mutations on the RBM. However, mutations away from the RBM can have a considerable impact on the infectivity as well.

The above analysis considers only BA strengthening mutations. To have a global view of how future Binding affinity changes (kcal/mol) Figure 10 : An illustration of the average and variance of ∆∆G (kcal/mol) for most likely mutation types on the RBM. y-axes: wild type residues; x-axises: mutant type residues. Colors on the axes indicate residue types. Variance (kcal/mol) Figure 11 : An illustration of the average and variance of ∆∆G (kcal/mol) for most likely mutation types away from the RBM. y-axes: wild type residues; x-axes: mutant type residues. Colors on the axes indicate residue types.

mutations would change the COVID-19 infectivity, we analyze the general trend of the free energy changes of most likely mutations according to 400 possible mutation types. The ∆∆G values following mutations on each amino acid are predicted and averaged by their mutation types. Figure 10 shows the average and variance of ∆∆G (kcal/mol) of each mutation type for most likely mutations on the RBD. Here, y-axes stand for wild type residues, and x-axes are mutant type residues. The colors on the axes are the residue types, such as charged, polar uncharged, hydrophobic, and special cases. The colors in the heat maps indicate the binding affinity changes strengths and directions. It is worthy to note that there are more positive binding affinity changes (green cubes) than negative changes (pink cubes), showing a trend of more infectious COVID19 strains due to most likely future mutations. For example, if a wild type mutation takes place from wide type K, T, N, Q, L, F, or Y to any other residue type except for W, it will end up with a more infectious COVID19 strain. However, mutations from R to T, or from A to many other residue types might lead to a less infectious COVID-19. The large values on the variance map indicate where the above average values might not be reliable. It is seen that the variances are general small. Figure 11 shows a similar trend for the most likely mutations away from the RDM.

As discussed above, likely and unlikely future mutations require two and three concurrent nucleotide mutations on each codon to happen, respectively. Figure 12 presents the top 20 likely future mutations that will strengthen the COVID-19 infectivity. The most energetic adversarial mutation is Y505G. Note that residue 505 has is a most likely mutation Y505H shown in Fig. 9 . Therefore, residue 505 is a potentially hot spot on the RBM. The next few energetic adversarial mutations are away from the RBM. Among them, N423P and N422G are hot-spot mutations. Figure 9 shows that residue 423 has three most likely energetic mutations while in Fig. 12 , it has the other three likely energetic mutations. Similarly, residue 489 on the RBM has 5 most likely energetic mutations (see Fig. 9 ) and 4 likely energetic mutations as shown in Fig. 12 . It is on our top surveillance list for the next generation of infectious COVID-19 strains. Another potentially hot spot is residue 495. Figure 13 shows the average and variance of ∆∆G (kcal/mol) of all likely mutations on the RBM where values of most likely mutations and unlikely mutations are excluded. About the same amount of mutations have positive and negative binding affinity changes. In Figure 14 , similar results are shown for the RBD excluding the RBM. Interestingly, mutations on the RBM have larger magnitude changes rather than out of this region for second potential mutations. It again shows that the RBM is the most important region to study. Binding affinity changes (kcal/mol) Figure 13 : An illustration of the average and variance of ∆∆G (kcal/mol) for likely mutation types on the RBM. y-axes: wild type residues; x-axises: mutant type residues. Colors on the axes indicate residue types. Variance (kcal/mol) Figure 14 : An illustration of the average and variance of ∆∆G (kcal/mol) for likely mutation types away from the RBM. y-axes: wild type residues; x-axises: mutant type residues. Colors on the axes indicate residue types. 

2.0 R H K D E S T N Q A V I L M F Y W C G P R H K D E S T N Q A V I L M F Y W C G P R H K D E S T N Q A V I L M F Y W C G

To further understand the evolution trend and potential infectivity changes of the COVID-19, we carry out S protein sequence alignment analysis to examine residue conservativeness. Figure 17 presents the alignment analysis of SARS-CoV-2 S protein sequence and those of the other four closely related species, namely, SARS-CoV, bat coronavirus RaTG13, bat coronavirus BM48-31, and bat coronavirus CoVZC45. We note that among the residues we discussed in the last section, 414, 422, 423, 492, and 495 are very conservative. They have not undergone any mutations among five related species. In contrast, RBM residues 452, 489, Variance (kcal/mol) 415 S1/S2 Figure 17 : Sequence alignments of SARS-CoV-2 S protein with those of closely related species, including SARS-CoV [12] , bat coronavirus RaTG13 [32] , bat coronavirus BM48-31 [3] , and bat coronaviruse CoVZC45 [7] . Detailed numbering is given according to SARS-CoV-2. Residue 364 Ala (A) of bat coronavirus BM48-31 is omitted. 500, 501, and 505 have a history of mutations and are non-conservative. Therefore, the predicted infectivitystrengthening mutations on these residues are more likely to happen.

As mentioned earlier, there are inconsistent assessments about relative infectivity between SARS-CoV and SARS-CoV-2 in the literature [22, 27, 30] . Our validated computational method can be employed to resolve this discrepancy. Based on the sequence alignment shown in Figure 17 , we conduct binding affinity change calculations for all relevant mutations on the SARS-CoV S protein RBD. Figure 18 illustrates the S protein-ACE2 binding affinity changes following the mutations from SARS-CoV to SARS-CoV-2. The SARS-CoV S protein and ACE2 complex 3D0G [13] is used as the wide type in our predictions. It is interesting to note that overall, there are more infectivity-strengthening mutations than infectivity-weakening mutations on the RBD. This is particularly true for mutations on the RBM. This result indicates that the SARS-CoV-2 sample collected on January 5, 2020, [31] is slightly more infectious than SARS-CoV found in 2003 [12] .

For a comparison between various SARS-CoV-2 subtypes and SARS-CoV of 2003, our results indicate that SARS-CoV-2 Clusters, I, II, IV, V, and VI become more infectious than SARS-CoV whereas SARS-CoV-2 Cluster III may have a similar infective rate with that of SARS-CoV.

Compared with SARS-CoV, SARS-CoV-2 has four extra residues, PRRA from 681 to 684, as shown in Figure 17 . It is believed that these extra residues might change SARS-CoV-2's behavior in ACE2 assisted entry of host cells [27] . However, this speculation has no qualitative nor quantitative validation at present.

Amino acid sequences and mutant data of the S protein used in the analysis were obtained from NCBI GenBank and GISAID [23] . After the first complete genome sequence of SARS-CoV-2 released on NCBI GenBank (Access number: NC 045512.2) [31] , there has been a large number of genome sequences. Other sequences from GenBank are as follows: bat coronavirus RaTG13 (MN996532.1) [32] , bat coronavirus BM48-31 (NC 014470.1) [3] and bat coronavirus CoVZC45 (MG772933.1) [7] . The mutant information of 13752 whole-genome sequences of S protein with high coverage of SARS-CoV-2 strains from the infected individuals around the world was obtained from the GISAID database [23] (https://www.gisaid.org/). Data without the exact submission date in GISAID were not considered.

Beta-CoV S protein structures were obtained from the RCSB Protein Data Bank: SARS-CoV RBD with ACE2 (PDB 3D0G) [13] and SARS-CoV-2 RBD with ACE2 (PDB 6M0J) [11] . The structures were presented by using PyMOL [20] . Sequences alignments were performed on SARS-CoV-2 S protein sequences by using MAFFT v7.388 [10] and on SARS-CoV-2 genome by using the Clustal Omega multiple sequence alignment with default parameters [24] .

The topology-based network tree (TopNetTree) is constructed by an innovative integration between the topological representation and NetTree for predicting protein-protein interaction binding affinity changes following mutation ∆∆G. In this work, TopNetTree is applied to predict the binding affinity changes of mutations that happened on RBD with ACE2 of SARS-CoV-2 after January 5, 2020. As shown in Fig. 19 , topology-based feature generation is the first step followed by a convolutional neural network (CNN)-assisted model. The topological representation uses element-and site-specific persistent homology to simplify the structural complexity of protein-protein complexes and encode vital biological information into topological invariants. NetTree is a recently developed deep learning algorithm that integrates the advantages of convolutional neural networks and gradient-boosting trees (GBT). In this section, we briefly describe the topology representation for machine learning training and prediction. Details can be found in the literature [28] .

The topology-based feature generation is built upon from persistence homology starting with simplicial complex and filtration. As a type of algebraic topology, persistence homology studies simplicial complex on discrete datasets under various settings. Among the many constructions, two that are widely used for point clouds are the Vietoris-Rips (VR) complex and alpha complex [4] which are applied in our approach. Built upon a simplicial complex, the topological invariants of a point-cloud dataset can be identified, such as the set of atoms in protein-protein interactions. Meanwhile, topological invariants (separated components, rings, and cavities) can be enumerated by counting the numbers referred to as Betti-0, Betti-1, and Betti-2, respectively. Thus taxing and uninformative features or calculations are fully abandoned, whereas geometric and topological characteristics persevere as data representation. Moreover, using persistent homology, the original 3D point-cloud data are characterized by topological barcodes that record the "birth" and "death" of each topological invariants and simplifying complicated structural representations of a PPI-complex. Although topology data presentation much simplifies the problem in many directions, better construction is required for it to extract patterns of different biological or chemical aspects. Before talking about more detailed feature generations, we first preset the constructions for a PPI complex into various subsets.

A m : atoms of the mutation sites.

A mn (r): atoms in the neighbourhood of the mutation site within a cut-off distance r. 

A ele (E): atoms in the system that has atoms of element type E. The distance matrix is specially designed such that it excludes the interactions between the atoms form the same set. For interactions between atoms a i and a j in set A and/or set B, the modified distance is defined as

where D e (a i , a j ) is the Euclidian distance between a i and a j .

In algebraic topology, molecular atoms can be treated as points presented by v 0 , v 1 , v 2 , ..., v k as k +1 affinely independent points. Simplicial complex, the essential building blocks, is a finite collection of sets of points K = {σ i }, and σ i are called linear combinations of these points in R n (n ≥ k). For instance, a 0-, 1-, 2-, or 3-simplex in geometry representation is a vertex, an edge, a triangle, or a tetrahedron, respectively. A simplicial complex K is valied if a face τ of a k-simplex σ i of K is also in K, such that τ ⊆ σ i and σ i ∈ K imply τ ∈ K and the non-empty intersection of any two simplieces is a face for both. Given a simplicial complex K, a k-chain is a finite formal sum of k-simplices; that is, i α i σ k i . The set of all k-chains of the simplicial complex K equipped with an algebraic field (typically, Z 2 ) forms an abelian group C k (K, Z 2 ), which means the coefficients a i are chosen from Z 2 . A boundary operator ∂ k :

Consequently, a important property of boundary operator, ∂ k−1 ∂ k = ∅, follows from that boundaries are boundaryless. Moreover the kth cycle group Z k = ker∂ k = {c ∈ C k | ∂ k c = ∅} is deffined to be the kernel of ∂ k , whose elements are called k-cycles; and the kth boundary group is the image of ∂ k+1 denoted as

The algebraic construction to connect a sequence of complexes by boundary maps is called a chain complex

−→ 0 and the kth homology group is the quotient group defined by

The key property of boundary operators implies B k ⊆ Z k ⊆ C k . The Betti numbers are defined by the ranks of kth homology group H k which counts k-dimensional holes, especially, β 0 = rank(H 0 ) reflects the number of connected components, β 1 = rank(H 1 ) reflects the number of loops, and β 2 = rank(H 2 ) reveals the number of voids or cavities. Together, the set of Betti numbers {β 0 , β 1 , β 2 , · · · } indicates the intrinsic topological property of a system.

Persistent homology is devised to track the multiscale topological information over different scales along a filtration [4] . A filtration of a topology space K is a nested sequence of subspaces {K t } t=0,...,m of K such that ∅ = K 0 ⊆ K 1 ⊆ K 2 ⊆ · · · ⊆ K m = K. Moreover, on this complex sequence, we obtain a sequence of chain complexes by homomorphisms: C * (K 0 ) → C * (K 1 ) → · · · → C * (K m ) and a homology sequence: H * (K 0 ) → H * (K 1 ) → · · · → H * (K m ), correspondingly. The p-persistent kth homology group of K t is defined as H t,p

. Intuitively, this homology group records the homology classes of K t that are persistent at least until K t+p . Under the filtration process, the persistent homology barcodes can be generated. Then the feature vectors can be constructed from these sets of intervals for machine learning models.

In a variety of vectorization methods, one discretizes the filtration parameter interval into bins and model the behavior of the barcodes in each bin [1] . To make use of advanced machine learning algorithms, we subdivide a filtration interval into bins of length. Then the numbers of persistence intervals are counted for each bin, such that birth events and death events can be represented. This approach gives us three feature vectors for each topological barcode for the machine learning method. Note for different discretizations, the characterization of birth and death might not be stable so that only Betti-0 (H 0 ) barcodes obtained from the VR filtration are applied in this approach. Intuitively, features generated by binned barcode vectorization can reflect the strength of atom bonds, van der Waals interactions, and can be easily incorporated into a CNN, which captures and discriminates local patterns. Another method of vectorization is to get the statistics of bar lengths, birth values, and death values, such as sum, maximum, minimum, mean, and standard derivation. This method is applied to vectorize Betti-1 (H 1 ) and Betti-2 (H 2 ) barcodes obtained from alpha complex filtration based on the facts that higher-dimensional barcodes are sparser than H 0 barcodes.

Prediction of binding affinity changes following mutation for PPIs is very challenging due to the complex dataset and 3D structures. A hybrid machine learning algorithm that integrates a CNN and GBT is designed to overcome difficulties. Briefly speaking, partial topologically simplified descriptions, specifically vectorized H 0 barcode feature, are converted into concise features by the CNN module. Then a GBT module is trained on the whole feature set for a robust predictor with effective control of overfitting.

TopGBT model. The gradient boosting tree (GBT) method produces a prediction model as an ensemble method which is a class of machine learning algorithms. It builds a powerful module for regression and classification problems from weak learners. By the assumption that the individual learners are likely to make different mistakes, the method using a summation of the weak learners to eliminate the overall error. Furthermore, a decision tree is added to the ensemble depending on the current prediction error on the training dataset. Thus this method (a topology-based GBT or TopGBT) is relatively robust against hyperparameter tuning and overfitting, especially for a moderate number of features. The GBT is shown for its robustness against overfitting, good performance for moderately small data sizes, and model interpretability. The current work uses the package provided by scikit-learn (v 0.23.0).

TopCNN model. CNN is a class of deep neural networks and is considered as the most successful architectures. CNN is a regularized case of a multilayer connected neural network, such that each neuron is connected locally to neurons in the next convolution layers and the weights are shared across different locations. To prepare the integration of CNN and GBT, CNN is treated as an intermediate model that converts vectorized H 0 features into a higher-level abstract feature for the downstream model.

TopNetTree model. A supervised CNN model with the PPI ∆∆G as labels is trained for extracting high-level features from H 0 barcodes. Once the model is set up, the flatten layer neural outputs of CNN are feed into a GBT model to rank their importance. Based on the importance, and ordered subset of CNN-trained features is combined with features constructed from high-dimensional topological barcodes, H 1 and H 2 into the final GBT model as shown in Fig. 19 . As for the parameters of the GBT model, 10 times 10-fold experiments are done for searching the optimal parameter setting.

The proposed TopNetTree method is trained on the SKEMPI 2.0 dataset [9] , which has 4,169 variants in 319 different complexes. A set "S8338" with 8,338 variants was derived from SKEMPI 2.0 dataset by setting the reverse mutation energy changes to the negative values of its original energy changes. To address the reliability of the TopNetTree method, we did the tenfold cross-validation on the SKEMPI 2.0 dataset with the averaged training accuracy, Pearson correlation coefficients R p , Kendall's τ , and the root mean square error (RMSE), being 0.98, 0.89, and 0.37 kcal/mol. As shown in Table 2 , these metrics are based on the average of ten ten-fold cross-validations which indicate TopNetTree is well trained. The performance test of tenfold cross-validation on dataset gives as R p = 0.84, τ = 0.60, and RMSE = 1.06 kcal/mol, which is of the same level of accuracy as the best in the literature [28] . 

The infectivity of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a vital factor for preventive measurements against coronavirus disease 2019 (COVID-19) and reopening the global economy [22, 27, 30] . However, it is very challenging to rigorously determine the viral infectivity experimentally and quantitatively compare the relative infectivity between SARS-CoV and SARS-CoV-2. These challenges are deteriorated by the continuous evolution of SARS-CoV-2 due to its existing 13752 SNP variants in six distinct clusters [29] . In the present work, we develop an advanced TopNetTree method based on algebraic topology and deep learning to predict the spike glycoprotein (S protein) and the host angiotensin-converting enzyme 2 (ACE2) binding affinity changes induced by mutations. Based on binding affinity changes, we reveal that mutations have made five out of six clusters of SARS-CoV-2 more infectious than the original virus found in Wuhan [31] . Additionally, based on sequence alignment and mutation-induced binding affinity changes, we show that SARS-CoV-2 [31] is slightly more infectious than SARS-CoV found in 2003 [12] . Finally, we systematically compute the binding affinity changes of all possible 3686 future mutations to unveil that the most likely mutations will further strengthen SARS-CoV infectivity. We predict that residues 452, 489, 500, 501, and 505 on the receptor-binding motif (RBM) have high chances to mutate into significantly more infectious COVID-19 strains.

Supporting Material is available for tables of (1) Six clusters of SARS-CoV-2 mutations on the RBD and predicted BA changes; (2), (3), and (4) Predicted BA changes for most likely, likely, and unlikely future mutations, respectively; and (5) Predicted BA changes for all mutations from SARS-CoV to SARS-CoV-2.

thank The IBM TJ Watson Research Center, The COVID-19 High Performance Computing Consortium, NVIDIA, and MSU HPPC for computational assistance. RW thanks Dr. Changchuan Yin for assistance.

Representability of algebraic topology for biomolecules in machine learning based scoring and virtual screening

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster

Genomic characterization of SARS-related coronavirus in european bats and classification of coronaviruses based on partial RNA-dependent RNA polymerase gene sequences

Topological persistence and simplification

Return of the coronavirus: 2019-nCoV

SARS-CoV-2 cell entry depends on ACE2 and tmprss2 and is blocked by a clinically proven protease inhibitor

Genomic characterization and infectivity of a novel SARS-like coronavirus in Chinese bats

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

SKEMPI 2.0: an updated benchmark of changes in protein-protein binding energy, kinetics and thermodynamics upon mutation

Mafft multiple sequence alignment software version 7: improvements in performance and usability

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

A major outbreak of severe acute respiratory syndrome in Hong Kong

Structural analysis of major species barriers between humans and palm civets for severe acute respiratory syndrome coronavirus infections

Structure of SARS coronavirus spike receptor-binding domain complexed with receptor

Bats are natural reservoirs of SARS-like coronaviruses

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

SKEMPI: a structural kinetic and energetic database of mutant protein interactions and its use in empirical models

mCSM-Lig: quantifying the effects of mutations on protein-small molecule affinity in genetic disease and emergence of drug resistance

Identification of two critical amino acid residues of the severe acute respiratory syndrome coronavirus spike protein for its variation in zoonotic tropism transition via a double substitution strategy

The PyMOL molecular graphics system, version 1.8

The FoldX web server: an online force field

Structural basis of receptor recognition by SARS-CoV-2

Gisaid: Global initiative on sharing all influenza data-from vision to reality

Clustal omega, accurate alignment of very large numbers of sequences

AB-Bind: antibody binding mutational database for computational affinity predictions

Cross-host evolution of severe acute respiratory syndrome coronavirus in palm civet and human

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

A topology-based network tree for the prediction of protein-protein binding affinity changes following mutation

Decoding SARS-CoV-2 transmission, evolution and ramification on COVID-19 diagnosis, vaccine, and medicine

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

A new coronavirus associated with human respiratory disease in China

A pneumonia outbreak associated with a new coronavirus of probable bat origin

This work was supported in part by NIH grant GM126189, NSF Grants DMS-1721024, DMS-1761320, and IIS1900473, Michigan Economic Development Corporation, Bristol-Myers Squibb, and Pfizer. The authors