Prediction of infectivity of SARS-CoV2: Mathematical model with analysis of docking simulation for spike proteins and angiotensin-converting enzyme 2.

Objectives: : Variants of a coronavirus (SARS-CoV-2) have been spreading in a global pandemic. Improved understanding of the infectivity of future new variants is important so that effective countermeasures against them can be quickly undertaken. In our research reported here, we aimed to predict the infectivity of SARS-CoV-2 by using a mathematical model with molecular simulation analysis, and we used phylogenetic analysis to determine the evolutionary distance of the spike protein gene (S gene) of SARS-CoV-2.
Methods: : We subjected the six variants and the wild type of spike protein and human angiotensin-converting enzyme 2 (ACE2) to molecular docking simulation analyses to understand the binding affinity of spike protein and ACE2. We then utilized regression analysis of the correlation coefficient of the mathematical model and the infectivity of SARS-CoV-2 to predict infectivity.
Results: : The evolutionary distance of the S gene correlated with the infectivity of SARS-CoV-2 variants. The calculated biding affinity for the mathematical model obtained with results of molecular docking simulation also correlated with the infectivity of SARS-CoV-2 variants. These results suggest that the data from the docking simulation for the receptor binding domain of variant spike proteins and human ACE2 were valuable for prediction of SARS-CoV-2 infectivity.
Conclusion: : We developed a mathematical model for prediction of SARS-CoV-2 variant infectivity by using binding affinity obtained via molecular docking and the evolutionary distance of the S gene.

S gene

ACE2

RBD

MD

PDB

Introduction

Variants of the novel coronavirus (SARS-CoV-2 [severe acute respiratory syndrome coronavirus 2]) that are responsible for the worldwide pandemic known as COVID-19 have led to difficulties in enacting countermeasures against infection, because many variants have occurred in succession (Bar-On et al., 2020; Cao et al., 2021; Tao et al., 2021). Infection control methods for SARS-CoV-2 have mainly included wearing face masks, avoiding close contact with other people (such as via lockdowns in urban areas), and providing multiple injections of vaccines (Adil et al., 2021).

The infectivity of SARS-CoV-2 variants has been rapidly changing, so understanding the infectivity of new variants is important for effective responses to the pandemic (Singh et al., 2021). These changes have resulted in higher infectivities of new SARS-CoV-2 variants compared with past variants (Allen et al., 2022; Araf et al., 2022; Chen and Wei, 2022; Davies et al., 2021; Faria et al., 2021; Radvak et al., 2021). These alterations in infectivity indicate that the infectivity of new coronavirus variants may be estimated by utilizing the evolutionary distance of the spike protein gene (S gene) between the wild type and the variants. We previously described a mathematical model in which we used docking simulation results to predict UDP-glucuronosyltransferase 1A1 conjugation capacity (Takaoka et al., 2019). By using the molecular docking simulation analyses, we found a plant leaf extract that inhibited the binding of SARS‑CoV‑2 spike protein to angiotensin-converting enzyme 2 (ACE2) (Hagiyama et al., 2022). Similarly, in silico docking data for a variant of SARS-CoV-2 spike protein and ACE2 may be used to estimate the infectivity of the new variant. In this research here, we choose these two approaches with mathematical models to achieve rapid and better understanding of the infectivity of new SARS-CoV-2 variants.

Materials and methods

Determination of the evolutionary distance between wild-type and variant s genes and infectivities of the variants

We chose six variants of SARS-CoV-2—alpha (Claro et al., 2021), beta (Tegally et al., 2021), gamma (Loconsole et al., 2021), delta (Cherian et al., 2021), omicron BA.1 (Dejnirattisai et al., 2021), and omicron BA.2 (Christensen et al., 2022)]—for this research. Table 1 shows the amino acid substitution of spike protein variants. We used the multiple sequence alignment method for the S gene of the SARS-CoV-2 variants, whose nucleotide sequences were obtained from NCBI (Accession no: MN908947) (National Library of Medicine US, 2022), to perform an analysis via the ClustalW program (Thompson et al., 1994). FastTree (Price et al., 2009) with default parameters then provided the phylogenetic tree and evolutionary distances between each mutant and wild type (Kimura and Ota, 1972). The infectivity of each SARS-CoV-2 variant was obtained from previous research, and we summarized these data and developed an infectivity index for our research (Allen et al., 2022; Araf et al., 2022; Chen and Wei, 2022; Davies et al., 2021; Faria et al., 2021; Radvak et al., 2021).

Table 1

Amino acid substitutions of spike proteins of SARS-CoV-2 variants.

SARS-CoV-2 strain	Mutations in spike protein
Alpha (B.1.1.7)	H69-V70del, Y144del, N501Y, A570D, D614G, P681H, T716I, S982A, D1118H
Beta (B.1.351)	D80A, D215G, L241-A243del, K417N, E484K, N501Y, D614G, A701V
Gamma (P.1)	L18F, T20N, P26S, D138Y, R190S, K417T, E484K, N501Y, D614G, H655Y, T1027I, V1176F
Delta (B.1.617.2)	T19R, E156-F157del, R158G, L452R, T478K, D614G, P681R, D950N
Omicron BA.1 (B.1.1.529/BA.1)	A67V, H69-V70del, T95I, G142D, V143-Y145del, N211del, L212I, ins214EPE, G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681H, N764K, D796Y, N856K, Q954H, N969K, L981F
Omicron BA.2 (B.1.1.529/BA.2)	T19I, L24-P26del, A27S, G142D, V213G, G339D, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, S477N, T478K, E484A, Q493R, Q498R, N501Y, Y505H, D614G, H655Y, N679K, P681H, N764K, D796Y, Q954H, N969K

Bold, amino acids in RBD. The information of amino acid substitutions are obtained from the following sources: Alpha, Beta and Gamma, https://covdb.stanford.edu/variants/; Delta, https://covariants.org/variants/21A.Delta; Omicron BA.1, https://www.cdc.gov/coronavirus/2019-ncov/science/science-briefs/scientific-brief-omicron-variant.html; Omicron BA.2, Perumal et al. (J Med Virol, DOI: 10.1002/jmv.27601, 2022).

Analysis of the three-dimensional structures of SARS-CoV-2 spike protein, analysis of docking of the receptor-binding domain (RBD) and ACE2 protein, and development of a mathematical model

The amino acid sequence of wild-type SARS-CoV-2 spike protein was obtained from UniProt (UniProt ID: P0DTC2) (UniProt Consortium, 2018). The three-dimensional structure of SARS-CoV-2 spike protein (Protein Data Bank [PDB] ID: 6ZGG) (Berman et al., 2000), which is an open state trimer, was used as the template for homology modeling of wild-type and variant SARS-CoV-2 spike proteins. With the amino acid sequences and the template structure, we used Molecular Operating Environment software (Chemical Computing Group, Montreal, Quebec, Canada) to perform the homology modeling. The model structures were then subjected to structural optimization with molecular dynamics (MD) simulation, as described in our previous report (Takaoka et al., 2021). After 10,000-ps MD simulation of each spike protein, docking analysis with ACE2 was performed as follows: The ACE2 structure was obtained from the PDB (PDB ID: 6M0J), an X-ray crystal RBD and ACE2 complex with a high quality among the reported ones (other PDB ID: 6VW1, 6LZG). To reduce the search space for docking, a partial structure of the RBD in the “up” conformation was cut from our simulated structure of the SARS-CoV-2 spike protein trimer. Docking sites were defined on the basis of the crystal structure of SARS-CoV-2 spike protein RBD bound with ACE2 (PDB ID: 6M0J). Two thousand docking runs were performed for each variant by using ZDOCK software (University of Massachusetts Medical School, Worcester, MA) with the following parameters as the docking site: amino acid no 401, 430, 433, 437, 439, 440, 459, 470, 471, 473, 477, 480, 482, 484–486, 489 for RBD and amino acid no 24, 27, 28, 30, 31, 34, 35, 37, 38, 41, 42, 79, 82, 83, 330, 353–355, 357, 393 for ACE2, respectively. To determine the stable complex group among the docking runs, the ZDOCK scores of the resulting complexes were clustered by using the group average clustering algorithm. The upper tail rule Mojena (1977) was applied to determine the number of clusters. We used the most stable ZDOCK score in the docking runs and the average ZDOCK score of the most stable complex group except the one with the most stable score as the binding affinity to derive a mathematical model according to our previous research (Takaoka et al., 2019).

To verify our mathematical model, we used leave-one-out cross-validation for infectivity of SARS-CoV-2 as follows: We removed one of the SARS-CoV-2 strains, determined the constants of the mathematical model by using the rest of the five strains, and then performed calculations to predict the infectivity for the excluded strain. We then confirmed the correlation between the reported infectivity and the results of our cross-validation.

Fig. 1 provides a flow diagram of this research.

Fig. 1

The method used for molecular simulation analysis.

Results

Evolutionary distances for s gene variants and results of docking of RBD with ACE2 protein

Table 2 shows the evolutionary distances between the wild type and the variants, as well as the binding affinities of RBD with ACE2. After multiple alignments of wild-type and variant S genes by using the ClustalW program, FastTree was used to determine the phylogenetic tree (data not shown) and evolutionary distances. The binding affinity of the SARS-CoV-2 spike protein RBD and ACE2 was determined by using the ZDOCK score as an indicator. ZDOCK scores were utilized in cluster analysis to determine the average ZDOCK score of the most stable cluster and ZDOCK score of the most stable complex for each variant. We used these results to develop a mathematical model to predict infectivities of SARS-CoV-2 variants.

Table 2

Evolutionary distances and binding affinities of SARS-CoV-2 spike proteins with ACE2.

Measure	Wild type	Alpha (B.1.1.7)	Beta (B.1.351)	Gamma (P.1)	Delta (B.1.617.2)	Omicron BA.1 (B.1.1.529/BA.1)	Omicron BA.2 (B.1.1.529/BA.2)
Evolutionary distance	0	0.00420	0.00446	0.01077	0.01209	0.01426	0.02600
ZDOCK scores of the most stable complexes	216.076	253.965	264.905	282.820	454.240	509.856	615.442
ZDOCK scores of the most stable clusters (mean ± SD)	219.633 ± 50.513 (n = 10)	247.478 (n = 1)	260.243 (n = 1)	270.947 (n = 1)	428.721 (n = 1)	417.824 ± 4.936 (n = 5)	595.340 ± 4.736 (n = 3)

Development of a mathematical model to estimate infectivities of SARS-CoV-2 wild type and variants

We developed a mathematical model to predict the infectivities of SARS-CoV-2 wild type and variants according to our previous research (Takaoka et al., 2019). According to Fig. 2 , the relationship between the evolutionary distance and infectivity of each variant can be represented by the following equation:where x is the evolutionary distance between the wild-type (wild) and a mutant SARS-CoV-2 spike protein. Constant values p and q were estimated by minimizing the sum of squared error between the calculated infectivity and the reported infectivity (ratio per wild type) (Allen et al., 2022; Araf et al., 2022; Chen and Wei, 2022; Davies et al., 2021; Faria et al., 2021; Radvak et al., 2021). We then derived the mathematical model for the binding affinity of RBD and ACE2 by using the results of docking analyses. The binding affinity for the mathematical model was derived from the following equation:where A is the average ZDOCK score of the most stable cluster and B is the ZDOCK score of the most stable complex. The relationship between the binding affinity and the reported infectivity (Fig. 3 ) can be represented by the following equation:

Fig. 2

Correlation between evolutionary distance and infectivity of SARS-CoV-2 variants. Collinearity between the evolutionary distance and infectivity was noted. The regression line was indicated by dotted line as an equation of “y = 201.06x + 0.5784″.

Fig. 3

Correlation between the biding affinity based on our docking simulation analyses and infectivity of SARS-CoV-2 variants. The binding affinity for our mathematical model which is based on the molecular simulation data can be used to predict the infectivity of SARS-CoV-2 variants. The regression curve was indicated by dotted line as an equation of “y = 1.9455ln(x) + 1.2029″.

Constant values s and t were estimated by minimizing the sum of squared error between the calculated infectivity and the reported infectivity (ratio per wild type) (Allen et al., 2022; Araf et al., 2022; Chen and Wei, 2022; Davies et al., 2021; Faria et al., 2021; Radvak et al., 2021). The infectivity of SARS-CoV-2 (P) is defined by the evolutionary distance and the docking result with ACE2:

Constant values u, v, and w were estimated by minimizing the sum of squared error between the calculated infectivity P and the reported infectivity of each SARS-CoV-2 variant V:where set M represents the SARS-CoV-2 variants whose infectivities are reported. These results suggest that our mathematical model can predict the infectivities of SARS-CoV-2 wild type and variants. Finally, the infectivity of SARS-CoV-2 variants is defined by the following equation:

In addition, the result of the cross-validation of our mathematical model showed that the predicted infectivity of each variant, which was calculated by our model using the data of the other five variants was strongly correlated with the reported infectivity (Fig. 4 ).

Fig. 4

Cross-validation of predicted infectivity of SARS-CoV-2 with reported infectivity. Predicted infectivity of SARS-CoV-2 significantly correlated with reported infectivity (Pearson correlation coefficients R = 0.9527, p = 0.0009).

Discussion

Reported infectivities of SARS-CoV-2 variants were correlated with both evolutionary distance of the S gene and the binding affinity for the mathematical model according to our in silico docking data (Figs. 2 and 3). In addition, our mathematical model including these two factors reproduced the reported infectivity (Fig. 4).

In previous research, only the RBD domain of the spike protein was used in the analysis of binding affinity with ACE2: no MD analysis of the structure optimization (Celik et al., 2022; Gan et al., 2021; Shah and Woo, 2021; Sitthiyotha and Chunsrivirot, 2021; Tragni et al., 2022) and/or local optimization of only the RBD domain of the spike protein (Celik et al., 2022; Gan et al., 2021; Kumar et al., 2022; Shah and Woo, 2021; Sitthiyotha and Chunsrivirot, 2021; Tragni et al., 2022; Xue et al., 2021). However, we chose the three-dimensional structure of the complete trimeric spike protein (PDB ID: 6ZGG) in homology modeling in this research. After MD analysis, we confirmed the trajectory of the root mean square deviation and the Ramachandran plot for proper quality of each variant spike protein structure. The RBD domain of the variant trimeric spike protein structures was then provided for the docking analyses with ACE2. Therefore, our experimental strategy achieved high accuracy because we used a molecular simulation method which possibly imitates the behavior of the biomolecules. In addition, our previous research suggested the value of this prediction that utilized our molecular simulation methods (Nakamura et al., 2017; Ohta et al., 2018; Qiang et al., 2017; Sakaeda et al., 2018; Sugawara et al., 2019; Yamamura et al., 2020).

Conclusion

We developed a mathematical model for prediction of SARS-CoV-2 variant infectivity. Additional analyses of binding affinities of other spike proteins of SARS-CoV-2 variants with ACE2 may provide greater accuracy, and such experiments are now in progress.

Funding

This work was supported by Grant-in-Aid for Scientific Research [grant numbers 21K12110 to Y.T. and 19K12202 to A.S].

Author contributions

Y.T. conceived and designed this research. Y.T. and S.A. preformed the analyses and acquired the data. Y.T., S.A., Y.M., and M.O. interpreted the data. Y.T., S.A., and K.M derived mathematical model. Y.T. and S.A. wrote the draft, and all authors reviewed and approved the manuscript.

Ethical approval statement

This research is not applicable because we performed computer analyses by using sequence data obtained from public database.

Author statement

Y.T. conceived and designed this research. Y.T. and S.A. preformed the analyses and acquired the data. Y.T., S.A., Y.M., and M.O. interpreted the data. Y.T., S.A., and K.M derived mathematical model. Y.T. and S.A. wrote the draft, and all authors reviewed and approved the manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.