Biophysical principles predict fitness of SARS-CoV-2 variants

Abstract

SARS-CoV-2 employs its spike protein's receptor binding domain (RBD) to enter host cells. The RBD is constantly subjected to immune responses, while requiring efficient binding to host cell receptors for successful infection. However, our understanding of how RBD's biophysical properties contribute to SARS-CoV-2's epidemiological fitness remains largely incomplete. Through a comprehensive approach, comprising large-scale sequence analysis of SARS-CoV-2 variants and the identification of a fitness function based on binding thermodynamics, we unravel the relationship between the biophysical properties of RBD variants and their contribution to viral fitness. We developed a biophysical model that uses statistical mechanics to map the molecular phenotype space, characterized by dissociation constants of RBD to ACE2, LY-CoV016, LY-CoV555, REGN10987, and S309, onto an epistatic fitness landscape. We validate our findings through experimentally measured and machine learning (ML) estimated binding affinities, coupled with infectivity data derived from population-level sequencing. Our analysis reveals that this model effectively predicts the fitness of novel RBD variants and can account for the epistatic interactions among mutations, including explaining the later reversal of Q493R. Our study sheds light on the impact of specific mutations on viral fitness and delivers a tool for predicting the future epidemiological trajectory of previously unseen or emerging low-frequency variants. These insights offer not only greater understanding of viral evolution but also potentially aid in guiding public health decisions in the battle against COVID-19 and future pandemics.

Keywords: SARS-CoV-2; antibody; fitness landscape; receptor binding domain; viral evolution.

Fig. 1.

Model illustration. Link between fitness and thermodynamics of protein folding and binding. High fitness variants may exhibit improved stability in the folded state or in the ACE2-bound state to facilitate cellular entry or have the capacity to destabilize bound-to-antibody states, thereby enabling evasion.

Fig. 2.

Biophysical model analysis. (A) Fit of the predictive model on the training set. Each dot represents a variant in the training dataset, plotted against the fitness derived from population data on the x-axis and the predicted fitness by the model on the y-axis ($R^2=0.97$). (B) Model’s performance on the testing set. The model’s predictions align well with the fitness observed in population, as reflected by $R^2=0.94$ suggesting that the model maintains a strong predictive power on unseen data. (C–G) The dependence of fitness function on the logarithm of each dissociation constant. Existing variants are depicted as yellow regions on the curve, and a histogram showing the distribution of the dissociation constants is provided beneath each plot. Red and blue dashed lines represent wt and BA.1 respectively.

Fig. 3.

Assessing predictive power of the model using experimental $K_D$ (A–C) and computed $K_D$ (D–F). (A) Fitness prediction for variants carrying the G446S mutation which was not included in the training set yields $R^2=0.92$. (B) $R^2$ derived from a model trained on variants excluding a specific mutation, then used to predict fitness of variants exhibiting that mutation. (C) Predictions of fitness compared with actual fitness trend for variants between Wuhan-Hu-1 and Omicron BA.1. Variants observed before May 2021 are used as training set for the model. The model uses experimental $K_D$. (D) Predicted fitness from biophysical model against actual fitness derived from population data. The model uses $K_D$ acquired from ML and is fit on Wuhan-Omicron set. Selected variants are highlighted. (E) Predicted fitness compared with actual fitness for variants observed in 2020, 2021, and 2022. (F) Predictions of fitness using ML derived $K_D$ over three-month rolling windows, compared with the actual fitness trends for variants between 2021 and 2023.

Fig. 4.

Epistasis analysis. (A) Predicted fitness values with the nonepistatic model of Obermeyer et al. and with our epistatic biophysical model plotted against the genome’s mutation count, for all mutation combinations with mutation T478K. We define “Max Fitness” as the maximum fitness prediction from our biophysical model. Max Fitness curve begins to plateau with a higher mutation count, demonstrating a diminishing returns effect in epistatic. (B) Pairwise (second order) interaction coefficients against the spatial distances between the corresponding residues, with mutations colored in accordance with the absolute value of their pairwise coefficient. (C) Coefficients of epistasis: Diagonal coefficients denote first-order interactions, whereas off-diagonal coefficients represent second-order interactions. Coefficients smaller than 0.01 have been masked for clarity.