Free Energies of the Disassembly of Viral Capsids from a Multiscale Molecular Simulation Approach

Coupling finite and boundary element methods to solve the Poisson-Boltzmann equation for electrostatics in molecular solvation

The Poisson-Boltzmann equation is widely used to model electrostatics in molecular systems. Available software packages solve it using finite difference, finite element, and boundary element methods, where the latter is attractive due to the accurate representation of the molecular surface and partial charges, and exact enforcement of the boundary conditions at infinity. However, the boundary element method is limited to linear equations and piecewise constant variations of the material properties. In this work, we present a scheme that couples finite and boundary elements for the linearised Poisson-Boltzmann equation, where the finite element method is applied in a confined solute region and the boundary element method in the external solvent region. As a proof-of-concept exercise, we use the simplest methods available: Johnson-Nédélec coupling with mass matrix and diagonal preconditioning, implemented using the Bempp-cl and FEniCSx libraries via their Python interfaces. We showcase our implementation by computing the polar component of the solvation free energy of a set of molecules using a constant and a Gaussian-varying permittivity. As validation, we compare against well-established finite difference solvers for an extensive binding energy data set, and with the finite difference code APBS (to 0.5%) for Gaussian permittivities. We also show scaling results from protein G B1 (955 atoms) up to immunoglobulin G (20,148 atoms). For small problems, the coupled method was efficient, outperforming a purely boundary integral approach. For Gaussian-varying permittivities, which are beyond the applicability of boundary elements alone, we were able to run medium to large-sized problems on a single workstation. The development of better preconditioning techniques and the use of distributed memory parallelism for larger systems remains an area for future work. We hope this work will serve as inspiration for future developments that consider space-varying field parameters, and mixed linear-nonlinear schemes for molecular electrostatics with implicit solvent models.

Studying the ssDNA loaded adeno-associated virus aggregation using coarse-grained molecular dynamics simulations

The aggregation of adeno-associated viral (AAV) capsids in an aqueous environment was investigated via coarse-grained molecular dynamics (CG-MD) simulations. The primary driving force and mechanism of the aggregation were investigated with or without single-strand DNA (ssDNA) loaded at various process temperatures. Capsid aggregation appeared to involve multiple residue interactions (i.e., hydrophobic, polar and charged residues) leading to complex protein aggregation. In addition, two aggregation mechanisms (i.e., the fivefold face-to-face contact and the edge-to-edge contact) were identified from this study. The ssDNA with its asymmetric structure could be the reason for destabilizing protein subunits and enhancing the interaction between the charged residues, and further result in the non-reversible face-to-face contact. At higher temperature, the capsid structure was found to be unstable with the significant size expansion of the loaded ssDNA which could be attributed to reduced number of intramolecular hydrogen bonds, the increased conformational deviations of protein subunits and the higher residue fluctuations. The CG-MD model was further validated with previous experimental and simulation data, including the full capsid size measurement and the capsid internal pressure. Thus, a good understanding of AAV capsid aggregation, instability and the role of ssDNA were revealed by applying the developed computational model.

Understanding Virus Structure and Dynamics through Molecular Simulations

Viral outbreaks remain a serious threat to human and animal populations and motivate the continued development of antiviral drugs and vaccines, which in turn benefits from a detailed understanding of both viral structure and dynamics. While great strides have been made in characterizing these systems experimentally, molecular simulations have proven to be an essential, complementary approach. In this work, we review the contributions of molecular simulations to the understanding of viral structure, functional dynamics, and processes related to the viral life cycle. Approaches ranging from coarse-grained to all-atom representations are discussed, including current efforts at modeling complete viral systems. Overall, this review demonstrates that computational virology plays an essential role in understanding these systems.

Calcium inhibits penetration of Alzheimer's Aβ1 -42 monomers into the membrane

Calcium ion regulation plays a crucial role in maintaining neuronal functions such as neurotransmitter release and synaptic plasticity. Copper (Cu²⁺ ) coordination to amyloid-β (Aβ) has accelerated Aβ_1-42 aggregation that can trigger calcium dysregulation by enhancing the influx of calcium ions by extensive perturbing integrity of the membranes. Aβ_1-42 aggregation, calcium dysregulation, and membrane damage are Alzheimer disease (AD) implications. To gain a detail of calcium ions' role in the full-length Aβ_1-42 and Aβ_1-42 -Cu²⁺ monomers contact, the cellular membrane before their aggregation to elucidate the neurotoxicity mechanism, we carried out 2.5 μs extensive molecular dynamics simulation (MD) to rigorous explorations of the intriguing feature of the Aβ_1-42 and Aβ_1-42 -Cu²⁺ interaction with the dimyristoylphosphatidylcholine (DMPC) bilayer in the presence of calcium ions. The outcome of the results compared to the same simulations without calcium ions. We surprisingly noted robust binding energies between the Aβ_1-42 and membrane observed in simulations containing without calcium ions and is two and a half fold lesser in the simulation with calcium ions. Therefore, in the case of the absence of calcium ions, N-terminal residues of Aβ_1-42 deeply penetrate from the surface to the center of the bilayer; in contrast to calcium ions presence, the N- and C-terminal residues are involved only in surface contacts through binding phosphate moieties. On the other hand, Aβ_1-42 -Cu²⁺ actively participated in surface bilayer contacts in the absence of calcium ions. These contacts are prevented by forming a calcium bridge between Aβ_1-42 -Cu²⁺ and the DMPC bilayer in the case of calcium ions presence. In a nutshell, Calcium ions do not allow Aβ_1-42 penetration into the membranes nor contact of Aβ_1-42 -Cu²⁺ with the membranes. These pieces of information imply that the calcium ions mediate the membrane perturbation via the monomer interactions but do not damage the membrane; they agree with the western blot experimental results of a higher concentration of calcium ions inhibit the membrane pore formation by Aβ peptides.

Electrostatic Interactions Are the Primary Determinant of the Binding Affinity of SARS-CoV-2 Spike RBD to ACE2: A Computational Case Study of Omicron Variants

To explore the mechanistic origin that determines the binding affinity of SARS-CoV-2 spike receptor binding domain (RBD) to human angiotensin converting enzyme 2 (ACE2), we constructed the homology models of RBD-ACE2 complexes of four Omicron subvariants (BA.1, BA.2, BA.3 and BA.4/5), and compared them with wild type complex (RBD_WT-ACE2) in terms of various structural dynamic properties by molecular dynamics (MD) simulations and binding free energy (BFE) calculations. The results of MD simulations suggest that the RBDs of all the Omicron subvariants (RBD_OMIs) feature increased global structural fluctuations when compared with RBD_WT. Detailed comparison of BFE components reveals that the enhanced electrostatic attractive interactions are the main determinant of the higher ACE2-binding affinity of RBD_OMIs than RBD_WT, while the weakened electrostatic attractive interactions determine RBD of BA.4/5 subvariant (RBD_BA.4/5) lowest ACE2-binding affinity among all Omicron subvariants. The per-residue BFE decompositions and the hydrogen bond (HB) networks analyses indicate that the enhanced electrostatic attractive interactions are mainly through gain/loss of the positively/negatively charged residues, and the formation or destruction of the interfacial HBs and salt bridges can also largely affect the ACE2-binding affinity of RBD. It is worth pointing out that since Q493R plays the most important positive contribution in enhancing binding affinity, the absence of this mutation in RBD_BA.4/5 results in a significantly weaker binding affinity to ACE2 than other Omicron subvariants. Our results provide insight into the role of electrostatic interactions in determining of the binding affinity of SARS-CoV-2 RBD to human ACE2.

Coarse-Graining ddCOSMO through an Interface between Tinker and the ddX Library

The domain decomposition conductor-like screening model is an efficient way to compute the solvation energy of solutes within a polarizable continuum medium in a linear scaling computational time. Despite its efficiency, the application to very large systems is still challenging. A possibility to further accelerate the algorithm is resorting to coarse-graining strategies. In this paper we present a preliminary interface between the molecular dynamics package Tinker and the ddX library. The interface was used to test a united atom coarse-graining strategy that allowed us to push ddCOSMO to its limits by computing solvation energies on systems with up to 7 million atoms. We first present benchmarks to find an optimal discretization, and then, we discuss the performance and results obtained with fine- and coarse-grained solvation energy calculations.

Quantitative electrostatic force tomography for virus capsids in interaction with an approaching nanoscale probe

Electrostatic interactions are crucial for the assembly, disassembly and stability of proteinaceous viral capsids. Moreover, at the molecular scale, elucidating the organization and structure of the capsid proteins in response to an approaching nanoprobe is a major challenge in biomacromolecular research. Here, we report on a generalized electrostatic model, based on the Poisson-Boltzmann equation, that quantifies the subnanometric electrostatic interactions between an AFM tip and a proteinaceous capsid from molecular snapshots. This allows us to describe the contributions of specific amino acids and atoms to the interaction force. We show validation results in terms of total electrostatic forces with previous semi-empirical generalized models at available length scales (d > 1 nm). Then, we studied the interaction of the Zika capsid with conical and spherical AFM tips in a tomography-type analysis to identify the most important residues and atoms, showing the localized nature of the interaction. This method can be employed for the interpretation of force microscopy experiments in fundamental virological characterization and in diverse nanomedicine applications, where specific regions of the protein cages are aimed to electrostatically interact with molecular sized functionalized inhibitors, or tailoring protein-cage functional properties for nucleic acid delivery.

Lactoferrin Inhibition of the Complex Formation between ACE2 Receptor and SARS CoV-2 Recognition Binding Domain

The present investigation focuses on the analysis of the interactions among human lactoferrin (LF), SARS-CoV-2 receptor-binding domain (RBD) and human angiotensin-converting enzyme 2 (ACE2) receptor in order to assess possible mutual interactions that could provide a molecular basis of the reported preventative effect of lactoferrin against CoV-2 infection. In particular, kinetic and thermodynamic parameters for the pairwise interactions among the three proteins were measured via two independent techniques, biolayer interferometry and latex nanoparticle-enhanced turbidimetry. The results obtained clearly indicate that LF is able to bind the ACE2 receptor ectodomain with significantly high affinity, whereas no binding to the RBD was observed up to the maximum “physiological” lactoferrin concentration range. Lactoferrin, above 1 µM concentration, thus appears to directly interfere with RBD–ACE2 binding, bringing about a measurable, up to 300-fold increase of the K_D value relative to RBD–ACE2 complex formation.

Towards optimal boundary integral formulations of the Poisson-Boltzmann equation for molecular electrostatics

The Poisson-Boltzmann equation offers an efficient way to study electrostatics in molecular settings. Its numerical solution with the boundary element method is widely used, as the complicated molecular surface is accurately represented by the mesh, and the point charges are accounted for explicitly. In fact, there are several well-known boundary integral formulations available in the literature. This work presents a generalized expression of the boundary integral representation of the implicit solvent model, giving rise to new forms to compute the electrostatic potential. Moreover, it proposes a strategy to build efficient preconditioners for any of the resulting systems, improving the convergence of the linear solver. We perform systematic benchmarking of a set of formulations and preconditioners, focusing on the time to solution, matrix conditioning, and eigenvalue spectrum. We see that the eigenvalue clustering is a good indicator of the matrix conditioning, and show that they can be easily manipulated by scaling the preconditioner. Our results suggest that the optimal choice is problem-size dependent, where a simpler direct formulation is the fastest for small molecules, but more involved second-kind equations are better for larger problems. We also present a fast Calderón preconditioner for first-kind formulations, which shows promising behavior for future analysis. This work sets the basis towards choosing the most convenient boundary integral formulation of the Poisson-Boltzmann equation for a given problem.

Martini 3 Model of Cellulose Microfibrils: On the Route to Capture Large Conformational Changes of Polysaccharides

High resolution data from all-atom molecular simulations is used to parameterize a Martini 3 coarse-grained (CG) model of cellulose I allomorphs and cellulose type-II fibrils. In this case, elementary molecules are represented by four effective beads centred in the positions of O2, O3, C6, and O6 atoms in the D-glucose cellulose subunit. Non-bonded interactions between CG beads are tuned according to a low statistical criterion of structural deviation using the Martini 3 type of interactions and are capable of being indistinguishable for all studied cases. To maintain the crystalline structure of each single cellulose chain in the microfibrils, elastic potentials are employed to retain the ribbon-like structure in each chain. We find that our model is capable of describing different fibril-twist angles associated with each type of cellulose fibril in close agreement with atomistic simulation. Furthermore, our CG model poses a very small deviation from the native-like structure, making it appropriate to capture large conformational changes such as those that occur during the self-assembly process. We expect to provide a computational model suitable for several new applications such as cellulose self-assembly in different aqueous solutions and the thermal treatment of fibrils of great importance in bioindustrial applications.

Assessing the Stability of Biological Fibrils by Molecular-Scale Simulations

The nanomechanical characterization of several biological fibrils that are the result of protein aggregation via molecular dynamics simulation is nowadays feasible, and together with atomic force microscopy experiments has widened our understanding of the forces in the regime of pN-nN and system sizes of about hundreds of nanometers. Several methodologies have been developed to achieve this target, and they range from the atomistic representation via molecular force fields to coarse-grained strategies that provide comparable results with experiments in a systematic way. In this chapter, we discuss several methodologies for the calculation of mechanical parameters, such as the elastic constants of relevant biological systems. They are presented together with details about parameterization and current limitations. Then, we discuss some of the applications of such methodologies for the description of bacterial filament and β-amyloid systems. Finally, the latest lines of development are discussed.

RNA Secondary Structures Regulate Adsorption of Fragments onto Flat Substrates

RNA is a functionally rich molecule with multilevel, hierarchical structures whose role in the adsorption to molecular substrates is only beginning to be elucidated. Here, we introduce a multiscale simulation approach that combines a tractable coarse-grained RNA structural model with an interaction potential of a structureless flat adsorbing substrate. Within this approach, we study the specific role of stem-hairpin and multibranch RNA secondary structure motifs on its adsorption phenomenology. Our findings identify a dual regime of adsorption for short RNA fragments with and without the secondary structure and underline the adsorption efficiency in both cases as a function of the surface interaction strength. The observed behavior results from an interplay between the number of contacts formed at the surface and the conformational entropy of the RNA molecule. The adsorption phenomenology of RNA seems to persist also for much longer RNAs as qualitatively observed by comparing the trends of our simulations with a theoretical approach based on an ideal semiflexible polymer chain.

Acidification induces condensation of the adenovirus core

The adenovirus (AdV) icosahedral capsid encloses a nucleoprotein core formed by the dsDNA genome bound to numerous copies of virus-encoded, positively charged proteins. For an efficient delivery of its genome, AdV must undergo a cascade of dismantling events from the plasma membrane to the nuclear pore. Throughout this uncoating process, the virion moves across potentially disruptive environments whose influence in particle stability is poorly understood. In this work we analyze the effect of acidic conditions on AdV particles by exploring their mechanical properties, genome accessibility and capsid disruption. Our results show that under short term acidification the AdV virion becomes softer and its genome less accessible to an intercalating dye, even in the presence of capsid openings. The AFM tip penetrates deeper in virions at neutral pH, and mechanical properties of genome-less particles are not altered upon acidification. Altogether, these results indicate that the main effect of acidification is the compaction of the nucleoproteic core, revealing a previously unknown role for chemical cues in AdV uncoating. STATEMENT OF SIGNIFICANCE: Studying the behavior of virus particles under changing environmental conditions is key to understand cell entry and propagation. One such change is the acidification undergone in certain cell compartments, which is thought to play a role in the programmed uncoating of virus genomes. Mild acidification in the early endosome has been proposed as a trigger signal for human AdV uncoating. However, the actual effect of low pH in AdV stability and entry is not well defined. Understanding the consequences of acidification in AdV structure and stability is also relevant to define storage conditions for therapeutic vectors, or design AdV variants resistant to intestinal conditions for oral administration of vaccines.

Structural 3D Domain Reconstruction of the RNA Genome from Viruses with Secondary Structure Models

Three-dimensional RNA domain reconstruction is important for the assembly, disassembly and delivery functionalities of a packed proteinaceus capsid. However, to date, the self-association of RNA molecules is still an open problem. Recent chemical probing reports provide, with high reliability, the secondary structure of diverse RNA ensembles, such as those of viral genomes. Here, we present a method for reconstructing the complete 3D structure of RNA genomes, which combines a coarse-grained model with a subdomain composition scheme to obtain the entire genome inside proteinaceus capsids based on secondary structures from experimental techniques. Despite the amount of sampling involved in the folded and also unfolded RNA molecules, advanced microscope techniques can provide points of anchoring, which enhance our model to include interactions between capsid pentamers and RNA subdomains. To test our method, we tackle the satellite tobacco mosaic virus (STMV) genome, which has been widely studied by both experimental and computational communities. We provide not only a methodology to structurally analyze the tertiary conformations of the RNA genome inside capsids, but a flexible platform that allows the easy implementation of features/descriptors coming from both theoretical and experimental approaches.

An Effective MM/GBSA Protocol for Absolute Binding Free Energy Calculations: A Case Study on SARS-CoV-2 Spike Protein and the Human ACE2 Receptor

The binding free energy calculation of protein-ligand complexes is necessary for research into virus-host interactions and the relevant applications in drug discovery. However, many current computational methods of such calculations are either inefficient or inaccurate in practice. Utilizing implicit solvent models in the molecular mechanics generalized Born surface area (MM/GBSA) framework allows for efficient calculations without significant loss of accuracy. Here, GBNSR6, a new flavor of the generalized Born model, is employed in the MM/GBSA framework for measuring the binding affinity between SARS-CoV-2 spike protein and the human ACE2 receptor. A computational protocol is developed based on the widely studied Ras-Raf complex, which has similar binding free energy to SARS-CoV-2/ACE2. Two options for representing the dielectric boundary of the complexes are evaluated: one based on the standard Bondi radii and the other based on a newly developed set of atomic radii (OPT1), optimized specifically for protein-ligand binding. Predictions based on the two radii sets provide upper and lower bounds on the experimental references: -14.7(ΔGbindBondi)<-10.6(ΔGbindExp.)<-4.1(ΔGbindOPT1) kcal/mol. The consensus estimates of the two bounds show quantitative agreement with the experiment values. This work also presents a novel truncation method and computational strategies for efficient entropy calculations with normal mode analysis. Interestingly, it is observed that a significant decrease in the number of snapshots does not affect the accuracy of entropy calculation, while it does lower computation time appreciably. The proposed MM/GBSA protocol can be used to study the binding mechanism of new variants of SARS-CoV-2, as well as other relevant structures.

Efficient mesh refinement for the Poisson-Boltzmann equation with boundary elements

The Poisson-Boltzmann equation is a widely used model to study electrostatics in molecular solvation. Its numerical solution using a boundary integral formulation requires a mesh on the molecular surface only, yielding accurate representations of the solute, which is usually a complicated geometry. Here, we utilize adjoint-based analyses to form two goal-oriented error estimates that allow us to determine the contribution of each discretization element (panel) to the numerical error in the solvation free energy. This information is useful to identify high-error panels to then refine them adaptively to find optimal surface meshes. We present results for spheres and real molecular geometries, and see that elements with large error tend to be in regions where there is a high electrostatic potential. We also find that even though both estimates predict different total errors, they have similar performance as part of an adaptive mesh refinement scheme. Our test cases suggest that the adaptive mesh refinement scheme is very effective, as we are able to reduce the error one order of magnitude by increasing the mesh size less than 20% and come out to be more efficient than uniform refinement when computing error estimations. This result sets the basis toward efficient automatic mesh refinement schemes that produce optimal meshes for solvation energy calculations.

Characterization of Structural and Energetic Differences between Conformations of the SARS-CoV-2 Spike Protein

The novel coronavirus disease 2019 (COVID-19) pandemic has disrupted modern societies and their economies. The resurgence in COVID-19 cases as part of the second wave is observed across Europe and the Americas. The scientific response has enabled a complete structural characterization of the Severe Acute Respiratory Syndrome-novel Coronavirus 2 (SARS-CoV-2). Among the most relevant proteins required by the novel coronavirus to facilitate the cell entry mechanism is the spike protein. This protein possesses a receptor-binding domain (RBD) that binds the cellular angiotensin-converting enzyme 2 (ACE2) and then triggers the fusion of viral and host cell membranes. In this regard, a comprehensive characterization of the structural stability of the spike protein is a crucial step to find new therapeutics to interrupt the process of recognition. On the other hand, it has been suggested that the participation of more than one RBD is a possible mechanism to enhance cell entry. Here, we discuss the protein structural stability based on the computational determination of the dynamic contact map and the energetic difference of the spike protein conformations via the mapping of the hydration free energy by the Poisson-Boltzmann method. We expect our result to foster the discussion of the number of RBD involved during recognition and the repurposing of new drugs to disable the recognition by discovering new hotspots for drug targets apart from the flexible loop in the RBD that binds the ACE2.

Quantitative determination of mechanical stability in the novel coronavirus spike protein

We report on the novel observation about the gain in nanomechanical stability of the SARS-CoV-2 (CoV2) spike (S) protein in comparison with SARS-CoV from 2002 (CoV1). Our findings have several biological implications in the subfamily of coronaviruses, as they suggest that the receptor binding domain (RBD) (∼200 amino acids) plays a fundamental role as a damping element of the massive viral particle's motion prior to cell-recognition, while also facilitating viral attachment, fusion and entry. The mechanical stability via pulling of the RBD is 250 pN and 200 pN for CoV2 and CoV1 respectively, and the additional stability observed for CoV2 (∼50 pN) might play a role in the increasing spread of COVID-19.

Nanomechanical mapping of soft materials with the atomic force microscope: methods, theory and applications

Fast, high-resolution, non-destructive and quantitative characterization methods are needed to develop materials with tailored properties at the nanoscale or to understand the relationship between mechanical properties and cell physiology. This review introduces the state-of-the-art force microscope-based methods to map at high-spatial resolution the elastic and viscoelastic properties of soft materials. The experimental methods are explained in terms of the theories that enable the transformation of observables into material properties. Several applications in materials science, molecular biology and mechanobiology illustrate the scope, impact and potential of nanomechanical mapping methods.