The 2019 novel coronavirus, SARS-CoV-2, is an emerging pathogen of critical significance to international public health. Knowledge of the interplay between molecular-scale virus-receptor interactions, single-cell viral replication, intracellular-scale viral transport, and emergent tissue-scale viral propagation is limited. Moreover, little is known about immune system-virus-tissue interactions and how these can result in low-level (asymptomatic) infections in some cases and acute respiratory distress syndrome (ARDS) in others, particularly with respect to presentation in different age groups or pre-existing inflammatory risk factors like diabetes. A critical question for treatment and protection is why it appears that the severity of infection may correlate with the initial level of virus exposure. Given the nonlinear interactions within and among each of these processes, multiscale simulation models can shed light on the emergent dynamics that lead to divergent outcomes, identify actionable "choke points" for pharmacologic interactions, screen potential therapies, and identify potential biomarkers that differentiate response dynamics. Given the complexity of the problem and the acute need for an actionable model to guide therapy discovery and optimization, we introduce a prototype of a multiscale model of SARS-CoV-2 dynamics in lung and intestinal tissue that will be iteratively refined. The first prototype model was built and shared internationally as open source code and interactive, cloud-hosted executables in under 12 hours. In a sustained community effort, this model will integrate data and expertise across virology, immunology, mathematical biology, quantitative systems physiology, cloud and high performance computing, and other domains to accelerate our response to this critical threat to international health.

The weighted ensemble (WE) strategy enables direct simulation of atomistic, fully-continuous protein–protein binding pathways in explicit solvent, yielding rigorous kinetics.

Despite more than three decades of effort with molecular dynamics simulations, long-timescale (ms and beyond) biologically relevant phenomena remain out of reach in most systems of interest. This is largely because important transitions, such as conformational changes and (un)binding events, tend to be rare for conventional simulations (< 10 µs). That is, conventional simulations will predominantly dwell in metastable states instead of making large transitions in complex biomolecular energy landscapes. In contrast, path sampling approaches focus computing effort specifically on transitions of interest. Such approaches have been in use for nearly 20 years in biomolecular systems and enabled the generation of pathways and calculation of rate constants for ms processes, including large protein conformational changes, protein folding, and protein (un)binding.

The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present five tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. Users are expected to already have significant experience with running standard molecular dynamics simulations using the underlying dynamics engine of interest (e.g. Amber, Gromacs, OpenMM). The tutorials range from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding.

An essential baseline for determining the extent to which electrostatic interactions enhance the kinetics of protein–protein association is the “basal” kon, which is the rate constant for association in the absence of electrostatic interactions. However, since such association events are beyond the milliseconds time scale, it has not been practical to compute the basal kon by directly simulating the association with flexible models. Here, we computed the basal kon for barnase and barstar, two of the most rapidly associating proteins, using highly efficient, flexible molecular simulations. These simulations involved (a) pseudoatomic protein models that reproduce the molecular shapes, electrostatic, and diffusion properties of all-atom models, and (b) application of the weighted ensemble path sampling strategy, which enhanced the efficiency of generating association events by >130-fold. We also examined the extent to which the computed basal kon is affected by inclusion of intermolecular hydrodynamic interactions in the simulations.

A grand challenge in the field of biophysics has been the complete characterization of proteinprotein binding processes at atomic resolution. This characterization requires the direct simulation of binding pathways starting from the initial unbound state and proceeding through states that are too transient to be captured by experiment. Here we applied the weighted ensemble path sampling strategy to enable atomistic simulation of protein-protein binding pathways. Our simulation generated 203 fully continuous binding pathways for the bacterial proteins, barnase and barstar, yielding a computed kon that is within error of experiment. Results reveal that the formation of the "encounter complex" intermediate is rate limiting with ~11% of all diffusional collisions being productive. Consistent with experiment, our simulations identify R59 as the most kinetically important barnase residue for the binding process. Furthermore, protein desolvation occurs late in the binding process during the rearrangement of the encounter complex to the native complex. Notably, the positions of interfacial crystallographic water molecules that bridge hydrogen bonds between barnase and barstar are occupied upon formation of the native complex in our simulations. Our simulations were completed within a month using 1600 CPU cores at a time, demonstrating that it is now practical to carry out atomistic simulations of protein-protein binding processes, particularly using the latest GPU-accelerated computing.

Many intrinsically disordered proteins, which are prevalent in nature, fold only upon binding their structured partner proteins. Such proteins have been hypothesized to have a kinetic advantage over their folded, preorganized analogues in binding their partner proteins. Here we determined the effects of ligand preorganization on the kon for a biomedically important system: an intrinsically disordered p53 peptide ligand and the MDM2 protein receptor. Based on direct simulations of binding pathways, computed kon values for fully disordered and preorganized p53 peptide analogues were within error of each other, indicating little if any kinetic advantage to being disordered or preorganized for binding the MDM2 protein. We also examined the effects of increasing the concentration of MDM2 on the extent to which its mechanism of binding to the p53 peptide is induced fit vs conformational selection. Results predict that the mechanism is solely induced fit if the unfolded state of the peptide is more stable than its folded state; otherwise, the mechanism shifts from being dominated by conformational selection at low MDM2 concentration to induced fit at high MDM2 concentration. Taken together, our results are relevant to any protein binding process that involves a disordered peptide of a similar length that forms a single α-helix upon binding a partner protein. Such disorder-to-helix transitions are common among protein interactions of disordered proteins and are therefore of fundamental biological interest.

Summary: Spatial heterogeneity can have dramatic effects on the biochemical networks that drive cell regulation and decision-making. For this reason, a number of methods have been developed to model spatial heterogeneity and incorporated into widely used modeling platforms. Unfortunately, the standard approaches for specifying and simulating chemical reaction networks become untenable when dealing with multi-state, multi-component systems that are characterized by combinatorial complexity. To address this issue, we developed MCell-R, a framework that extends the particle-based spatial Monte Carlo simulator, MCell, with the rule-based model specification and simulation capabilities provided by BioNetGen and NFsim. The BioNetGen syntax enables the specification of biomolecules as structured objects whose components can have different internal states that represent such features as covalent modification and conformation and which can bind components of other molecules to form molecular complexes. The network-free simulation algorithm used by NFsim enables efficient simulation of rule-based models even when the size of the network implied by the biochemical rules is too large to enumerate explicitly, which frequently occurs in detailed models of biochemical signaling. The result is a framework that can efficiently simulate systems characterized by combinatorial complexity at the level of spatially-resolved individual molecules over biologically relevant time and length scales.