Search for: All records

Mutations in the new variants of SARS-CoV-2 spike protein modulates the dynamics of the neutralizing antibodies. Capturing such modulations from MD simulations and graph network model identifies the role of mutations in facilitating immune evasion. 

Infection by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) involves the attachment of the receptor-binding domain (RBD) of its spike proteins to the ACE2 receptors on the peripheral membrane of host cells. Binding is initiated by a down-to-up conformational change in the spike protein, the change that presents the RBD to the receptor. To date, computational and experimental studies that search for therapeutics have concentrated, for good reason, on the RBD. However, the RBD region is highly prone to mutations, and is therefore a hotspot for drug resistance. In contrast, we here focus on the correlations between the RBD and residues distant to it in the spike protein. This allows for a deeper understanding of the underlying molecular recognition events and prediction of the highest-effect key mutations in distant, allosteric sites, with implications for therapeutics. Also, these sites can appear in emerging mutants with possibly higher transmissibility and virulence, and preidentifying them can give clues for designing pan-coronavirus vaccines against future outbreaks. Our model, based on time-lagged independent component analysis (tICA) and protein graph connectivity network, is able to identify multiple residues that exhibit long-distance coupling with the RBD opening. Residues involved in the most ubiquitous D614G mutation and the A570D mutation of the highly contagious UK SARS-CoV-2 variant are predicted ab initio from our model. Conversely, broad-spectrum therapeutics like drugs and monoclonal antibodies can target these key distant-but-conserved regions of the spike protein. 

Genetic information is encoded in the DNA double helix, which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleo-base pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen (HG) base pairing pattern in naked DNA duplexes, and estimated its relative stability and lifetime. In contrast with DNA, such structures were not observed in RNA duplexes. Understanding HG base pairing is important because the underlying "breathing" motion between the two conformations can significantly modulate protein binding. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force-bias method) and Markov state modeling to obtain accurate free energy, kinetics, and the intermediates in the transition pathway between Watson-Crick and HG base pairs for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extrahelical intermediates in the seemingly simple process of base flipping in B-DNA. Extending our calculation to A-RNA, for which HG base pairing is not observed experimentally, resulted in relatively unstable, single-hydrogen-bonded, distorted Hoogsteen-like bases. Unlike B-DNA, the transition pathway primarily involved base paired and intrahelical intermediates with transition timescales much longer than that of B-DNA. The seemingly obvious flip-over reaction coordinate (i.e., the glycosidic torsion angle) is unable to resolve the intermediates. Instead, a multidimensional picture involving backbone dihedral angles and distance between hydrogen bond donor and acceptor atoms is required to gain insight into the molecular mechanism. 

To directly simulate rare events using atomistic molecular dynamics is a significant challenge in computational biophysics. Well-established enhanced-sampling techniques do exist to obtain the thermodynamic functions for such systems. However, developing methods for obtaining the kinetics of long timescale processes from simulation at atomic detail is comparatively less developed an area. Milestoning and the weighted ensemble (WE) method are two different stratification strategies; both have shown promise for computing long timescales of complex biomolecular processes. Nevertheless, both require a significant investment of computational resources. We have combined WE and milestoning to calculate observables in orders-of-magnitude less central processing unit and wall-clock time. Our weighted ensemble milestoning method (WEM) uses WE simulation to converge the transition probability and first passage times between milestones, followed by the utilization of the theoretical framework of milestoning to extract thermodynamic and kinetic properties of the entire process. We tested our method for a simple one-dimensional double-well potential, for an eleven-dimensional potential energy surface with energy barrier, and on the biomolecular model system alanine dipeptide. We were able to recover the free energy profiles, time correlation functions, and mean first passage times for barrier crossing events at a significantly small computational cost. WEM promises to extend the applicability of molecular dynamics simulation to slow dynamics of large systems that are well beyond the scope of present day brute-force computations.