Recent advances in electron microscopy have enabled the imaging of single cells in 3D at nanometer length scale resolutions. An uncharted frontier for in silico biology is the ability to simulate cellular processes using these observed geometries. Enabling such simulations requires watertight meshing of electron micrograph images into 3D volume meshes, which can then form the basis of computer simulations of such processes using numerical techniques such as the finite element method. In this paper, we describe the use of our recently rewritten mesh processing software, GAMer 2, to bridge the gap between poorly conditioned meshes generated from segmented micrographs and boundary marked tetrahedral meshes which are compatible with simulation. We demonstrate the application of a workflow using GAMer 2 to a series of electron micrographs of neuronal dendrite morphology explored at three different length scales and show that the resulting meshes are suitable for finite element simulations. This work is an important step towards making physical simulations of biological processes in realistic geometries routine. Innovations in algorithms to reconstruct and simulate cellular length scale phenomena based on emerging structural data will enable realistic physical models and advance discovery at the interface of geometry and cellular processes. We posit that a new frontier at the intersection of computational technologies and single cell biology is now open.

We introduce CASC: a new, modern, and header-only C++ library which provides a data structure to represent arbitrary dimension abstract simplicial complexes (ASC) with user-defined classes stored directly on the simplices at each dimension. This is accomplished by using the latest C++ language features including variadic template parameters introduced in C++11 and automatic function return type deduction from C++14. Effectively CASC decouples the representation of the topology from the interactions of user data. We present the innovations and design principles of the data structure and related algorithms. This includes a metadata aware decimation algorithm which is general for collapsing simplices of any dimension. We also present an example application of this library to represent an orientable surface mesh.

BACE-1 is the β-secretase responsible for the initial amyloidogenesis in Alzheimer's disease, catalyzing hydrolytic cleavage of substrate in a pH-sensitive manner. The catalytic mechanism of BACE-1 requires water-mediated proton transfer from aspartyl dyad to the substrate, as well as structural flexibility in the flap region. Thus, the coupling of protonation and conformational equilibria is essential to a full in silico characterization of BACE-1. In this work, we perform constant pH replica exchange molecular dynamics simulations on both apo BACE-1 and five BACE-1-inhibitor complexes to examine the effect of pH on dynamics and inhibitor binding properties of BACE-1. In our simulations, we find that solution pH controls the conformational flexibility of apo BACE-1, whereas bound inhibitors largely limit the motions of the holo enzyme at all levels of pH. The microscopic pKa values of titratable residues in BACE-1 including its aspartyl dyad are computed and compared between apo and inhibitor-bound states. Changes in protonation between the apo and holo forms suggest a thermodynamic linkage between binding of inhibitors and protons localized at the dyad. Utilizing our recently developed computational protocol applying the binding polynomial formalism to the constant pH molecular dynamics (CpHMD) framework, we are able to obtain the pH-dependent binding free energy profiles for various BACE-1-inhibitor complexes. Our results highlight the importance of correctly addressing the binding-induced protonation changes in protein-ligand systems where binding accompanies a net proton transfer. This work comprises the first application of our CpHMD-based free energy computational method to protein-ligand complexes and illustrates the value of CpHMD as an all-purpose tool for obtaining pH-dependent dynamics and binding free energies of biological systems.

The nonmevalonate pathway is responsible for isoprenoid production in microbes, including H. pylori, M. tuberculosis and P. falciparum, but is nonexistent in humans, thus providing a desirable route for antibacterial and antimalarial drug discovery. We coordinate a structural study of IspH, a [4Fe-4S] protein responsible for converting HMBPP to IPP and DMAPP in the ultimate step in the nonmevalonate pathway. By performing accelerated molecular dynamics simulations on both substrate-free and HMBPP-bound [Fe4S4](2+) IspH, we elucidate how substrate binding alters the dynamics of the protein. Using principal component analysis, we note that while substrate-free IspH samples various open and closed conformations, the closed conformation observed experimentally for HMBPP-bound IspH is inaccessible in the absence of HMBPP. In contrast, simulations with HMBPP bound are restricted from accessing the open states sampled by the substrate-free simulations. Further investigation of the substrate-free simulations reveals large fluctuations in the HMBPP binding pocket, as well as allosteric pocket openings - both of which are achieved through the hinge motions of the individual domains in IspH. Coupling these findings with solvent mapping and various structural analyses reveals alternative druggable sites that may be exploited in future drug design efforts.

Group VI Ca²⁺-independent phospholipase A₂ (iPLA₂) is a water-soluble enzyme that is active when associated with phospholipid membranes. Despite its clear pharmaceutical relevance, no X-ray or NMR structural information is currently available for the iPLA₂ or its membrane complex. In this paper, we combine homology modeling with coarse-grained (CG) and all-atom (AA) molecular dynamics (MD) simulations to build structural models of iPLA₂ in association with a phospholipid bilayer. CG-MD simulations of the membrane insertion process were employed to provide a starting point for an atomistic description. Six AA-MD simulations were then conducted for 60 ns, starting from different initial CG structures, to refine the membrane complex. The resulting structures are shown to be consistent with each other and with deuterium exchange mass spectrometry (DXMS) experiments, suggesting that our approach is suitable for the modeling of iPLA₂ at the membrane surface. The models show that an anchoring region (residues 710-724) forms an amphipathic helix that is stabilized by the membrane. In future studies, the proposed iPLA₂ models should provide a structural basis for understanding the mechanisms of lipid extraction and drug-inhibition. In addition, the dual-resolution approach discussed here should provide the means for the future exploration of the impact of lipid diversity and sequence mutations on the activity of iPLA₂ and related enzymes.

Academic researchers and many in industry often lack the financial resources available to scientists working in "big pharma." High costs include those associated with high-throughput screening and chemical synthesis. In order to address these challenges, many researchers have in part turned to alternate methodologies. Virtual screening, for example, often substitutes for high-throughput screening, and click chemistry ensures that chemical synthesis is fast, cheap, and comparatively easy. Though both in silico screening and click chemistry seek to make drug discovery more feasible, it is not yet routine to couple these two methodologies. We here present a novel computer algorithm, called AutoClickChem, capable of performing many click-chemistry reactions in silico. AutoClickChem can be used to produce large combinatorial libraries of compound models for use in virtual screens. As the compounds of these libraries are constructed according to the reactions of click chemistry, they can be easily synthesized for subsequent testing in biochemical assays. Additionally, in silico modeling of click-chemistry products may prove useful in rational drug design and drug optimization. AutoClickChem is based on the pymolecule toolbox, a framework that may facilitate the development of future python-based programs that require the manipulation of molecular models. Both the pymolecule toolbox and AutoClickChem are released under the GNU General Public License version 3 and are available for download from http://autoclickchem.ucsd.edu.

Troponin C (TnC) is implicated in the initiation of myocyte contraction via binding of cytosolic Ca²⁺ and subsequent recognition of the Troponin I switch peptide. Mutations of the cardiac TnC N-terminal regulatory domain have been shown to alter both calcium binding and myofilament force generation. We have performed molecular dynamics simulations of engineered TnC variants that increase or decrease Ca²⁺ sensitivity, in order to understand the structural basis of their impact on TnC function. We will use the distinction for mutants that are associated with increased Ca²⁺ affinity and for those mutants with reduced affinity. Our studies demonstrate that for GOF mutants V44Q and L48Q, the structure of the physiologically-active site II Ca²⁺ binding site in the Ca²⁺-free (apo) state closely resembled the Ca²⁺-bound (holo) state. In contrast, site II is very labile for LOF mutants E40A and V79Q in the apo form and bears little resemblance with the holo conformation. We hypothesize that these phenomena contribute to the increased association rate, k(on), for the GOF mutants relative to LOF. Furthermore, we observe significant positive and negative positional correlations between helices in the GOF holo mutants that are not found in the LOF mutants. We anticipate these correlations may contribute either directly to Ca²⁺ affinity or indirectly through TnI association. Our observations based on the structure and dynamics of mutant TnC provide rationale for binding trends observed in GOF and LOF mutants and will guide the development of inotropic drugs that target TnC.

Periplasmic binding proteins (PBPs) are a large family of molecular transporters that play a key role in nutrient uptake and chemotaxis in Gram-negative bacteria. All PBPs have characteristic two-domain architecture with a central interdomain ligand-binding cleft. Upon binding to their respective ligands, PBPs undergo a large conformational change that effectively closes the binding cleft. This conformational change is traditionally viewed as a ligand induced-fit process; however, the intrinsic dynamics of the protein may also be crucial for ligand recognition. Recent NMR paramagnetic relaxation enhancement (PRE) experiments have shown that the maltose binding protein (MBP) - a prototypical member of the PBP superfamily - exists in a rapidly exchanging (ns to µs regime) mixture comprising an open state (approx 95%), and a minor partially closed state (approx 5%). Here we describe accelerated MD simulations that provide a detailed picture of the transition between the open and partially closed states, and confirm the existence of a dynamical equilibrium between these two states in apo MBP. We find that a flexible part of the protein called the balancing interface motif (residues 175-184) is displaced during the transformation. Continuum electrostatic calculations indicate that the repacking of non-polar residues near the hinge region plays an important role in driving the conformational change. Oscillations between open and partially closed states create variations in the shape and size of the binding site. The study provides a detailed description of the conformational space available to ligand-free MBP, and has implications for understanding ligand recognition and allostery in related proteins.

Chagas' disease, caused by the protozoan parasite Trypanosoma cruzi (T. cruzi), is a life-threatening illness affecting 11-18 million people. Currently available treatments are limited, with unacceptable efficacy and safety profiles. Recent studies have revealed an essential T. cruzi proline racemase enzyme (TcPR) as an attractive candidate for improved chemotherapeutic intervention. Conformational changes associated with substrate binding to TcPR are believed to expose critical residues that elicit a host mitogenic B-cell response, a process contributing to parasite persistence and immune system evasion. Characterization of the conformational states of TcPR requires access to long-time-scale motions that are currently inaccessible by standard molecular dynamics simulations. Here we describe advanced accelerated molecular dynamics that extend the effective simulation time and capture large-scale motions of functional relevance. Conservation and fragment mapping analyses identified potential conformational epitopes located in the vicinity of newly identified transient binding pockets. The newly identified open TcPR conformations revealed by this study along with knowledge of the closed to open interconversion mechanism advances our understanding of TcPR function. The results and the strategy adopted in this work constitute an important step toward the rationalization of the molecular basis behind the mitogenic B-cell response of TcPR and provide new insights for future structure-based drug discovery.

Paromomycin is an aminoglycosidic antibiotic that targets the RNA of the bacterial small ribosomal subunit. It binds in the A-site, which is one of the three tRNA binding sites, and affects translational fidelity by stabilizing two adenines (A1492 and A1493) in the flipped-out state. Experiments have shown that various mutations in the A-site result in bacterial resistance to aminoglycosides. In this study, we performed multiple molecular dynamics simulations of the mutated A-site RNA fragment in explicit solvent to analyze changes in the physicochemical features of the A-site that were introduced by substitutions of specific bases. The simulations were conducted for free RNA and in complex with paromomycin. We found that the specific mutations affect the shape and dynamics of the binding cleft as well as significantly alter its electrostatic properties. The most pronounced changes were observed in the U1406C∶U1495A mutant, where important hydrogen bonds between the RNA and paromomycin were disrupted. The present study aims to clarify the underlying physicochemical mechanisms of bacterial resistance to aminoglycosides due to target mutations.

GroEL is an ATP dependent molecular chaperone that promotes the folding of a large number of substrate proteins in E. coli. Large-scale conformational transitions occurring during the reaction cycle have been characterized from extensive crystallographic studies. However, the link between the observed conformations and the mechanisms involved in the allosteric response to ATP and the nucleotide-driven reaction cycle are not completely established. Here we describe extensive (in total long) unbiased molecular dynamics (MD) simulations that probe the response of GroEL subunits to ATP binding. We observe nucleotide dependent conformational transitions, and show with multiple 100 ns long simulations that the ligand-induced shift in the conformational populations are intrinsically coded in the structure-dynamics relationship of the protein subunit. Thus, these simulations reveal a stabilization of the equatorial domain upon nucleotide binding and a concomitant "opening" of the subunit, which reaches a conformation close to that observed in the crystal structure of the subunits within the ADP-bound oligomer. Moreover, we identify changes in a set of unique intrasubunit interactions potentially important for the conformational transition.

The t-tubules of mammalian ventricular myocytes are invaginations of the cell membrane that occur at each Z-line. These invaginations branch within the cell to form a complex network that allows rapid propagation of the electrical signal, and hence synchronous rise of intracellular calcium (Ca(2+)). To investigate how the t-tubule microanatomy and the distribution of membrane Ca(2+) flux affect cardiac excitation-contraction coupling we developed a 3-D continuum model of Ca(2+) signaling, buffering and diffusion in rat ventricular myocytes. The transverse-axial t-tubule geometry was derived from light microscopy structural data. To solve the nonlinear reaction-diffusion system we extended SMOL software tool (http://mccammon.ucsd.edu/smol/). The analysis suggests that the quantitative understanding of the Ca(2+) signaling requires more accurate knowledge of the t-tubule ultra-structure and Ca(2+) flux distribution along the sarcolemma. The results reveal the important role for mobile and stationary Ca(2+) buffers, including the Ca(2+) indicator dye. In agreement with experiment, in the presence of fluorescence dye and inhibited sarcoplasmic reticulum, the lack of detectible differences in the depolarization-evoked Ca(2+) transients was found when the Ca(2+) flux was heterogeneously distributed along the sarcolemma. In the absence of fluorescence dye, strongly non-uniform Ca(2+) signals are predicted. Even at modest elevation of Ca(2+), reached during Ca(2+) influx, large and steep Ca(2+) gradients are found in the narrow sub-sarcolemmal space. The model predicts that the branched t-tubule structure and changes in the normal Ca(2+) flux density along the cell membrane support initiation and propagation of Ca(2+) waves in rat myocytes.

Conventional drug design embraces the "one gene, one drug, one disease" philosophy. Polypharmacology, which focuses on multi-target drugs, has emerged as a new paradigm in drug discovery. The rational design of drugs that act via polypharmacological mechanisms can produce compounds that exhibit increased therapeutic potency and against which resistance is less likely to develop. Additionally, identifying multiple protein targets is also critical for side-effect prediction. One third of potential therapeutic compounds fail in clinical trials or are later removed from the market due to unacceptable side effects often caused by off-target binding. In the current work, we introduce a multidimensional strategy for the identification of secondary targets of known small-molecule inhibitors in the absence of global structural and sequence homology with the primary target protein. To demonstrate the utility of the strategy, we identify several targets of 4,5-dihydroxy-3-(1-naphthyldiazenyl)-2,7-naphthalenedisulfonic acid, a known micromolar inhibitor of Trypanosoma brucei RNA editing ligase 1. As it is capable of identifying potential secondary targets, the strategy described here may play a useful role in future efforts to reduce drug side effects and/or to increase polypharmacology.

Ras proteins regulate signaling cascades crucial for cell proliferation and differentiation by switching between GTP- and GDP-bound conformations. Distinct Ras isoforms have unique physiological functions with individual isoforms associated with different cancers and developmental diseases. Given the small structural differences among isoforms and mutants, it is currently unclear how these functional differences and aberrant properties arise. Here we investigate whether the subtle differences among isoforms and mutants are associated with detectable dynamical differences. Extensive molecular dynamics simulations reveal that wild-type K-Ras and mutant H-Ras A59G are intrinsically more dynamic than wild-type H-Ras. The crucial switch 1 and switch 2 regions along with loop 3, helix 3, and loop 7 contribute to this enhanced flexibility. Removing the gamma-phosphate of the bound GTP from the structure of A59G led to a spontaneous GTP-to-GDP conformational transition in a 20-ns unbiased simulation. The switch 1 and 2 regions exhibit enhanced flexibility and correlated motion when compared to non-transitioning wild-type H-Ras over a similar timeframe. Correlated motions between loop 3 and helix 5 of wild-type H-Ras are absent in the mutant A59G reflecting the enhanced dynamics of the loop 3 region. Taken together with earlier findings, these results suggest the existence of a lower energetic barrier between GTP and GDP states of the mutant. Molecular dynamics simulations combined with principal component analysis of available Ras crystallographic structures can be used to discriminate ligand- and sequence-based dynamic perturbations with potential functional implications. Furthermore, the identification of specific conformations associated with distinct Ras isoforms and mutants provides useful information for efforts that attempt to selectively interfere with the aberrant functions of these species.

Ras mediates signaling pathways controlling cell proliferation and development by cycling between GTP- and GDP-bound active and inactive conformational states. Understanding the complete reaction path of this conformational change and its intermediary structures is critical to understanding Ras signaling. We characterize nucleotide-dependent conformational transition using multiple-barrier-crossing accelerated molecular dynamics (aMD) simulations. These transitions, achieved for the first time for wild-type Ras, are impossible to observe with classical molecular dynamics (cMD) simulations due to the large energetic barrier between end states. Mapping the reaction path onto a conformer plot describing the distribution of the crystallographic structures enabled identification of highly populated intermediate structures. These structures have unique switch orientations (residues 25-40 and 57-75) intermediate between GTP and GDP states, or distinct loop3 (46-49), loop7 (105-110), and alpha5 C-terminus (159-166) conformations distal from the nucleotide-binding site. In addition, these barrier-crossing trajectories predict novel nucleotide-dependent correlated motions, including correlations of alpha2 (residues 66-74) with alpha3-loop7 (93-110), loop2 (26-37) with loop10 (145-151), and loop3 (46-49) with alpha5 (152-167). The interconversion between newly identified Ras conformations revealed by this study advances our mechanistic understanding of Ras function. In addition, the pattern of correlated motions provides new evidence for a dynamic linkage between the nucleotide-binding site and the membrane interacting C-terminus critical for the signaling function of Ras. Furthermore, normal mode analysis indicates that the dominant collective motion that occurs during nucleotide-dependent conformational exchange, and captured in aMD (but absent in cMD) simulations, is a low-frequency motion intrinsic to the structure.

We used molecular dynamics (MD) simulations to explore the transport of single cations through the channel of the muscle nicotinic acetylcholine receptor (nAChR). Four MD simulations of 16 ns were performed at physiological and hyperpolarized membrane potentials, with and without restraints of the structure, but all without bound agonist. With the structure unrestrained and a potential of -100 mV, one cation traversed the channel during a transient period of channel hydration; at -200 mV, the channel was continuously hydrated and two cations traversed the channel. With the structure restrained, however, cations did not traverse the channel at either membrane potential, even though the channel was continuously hydrated. The overall results show that cation selective transport through the nAChR channel is governed by electrostatic interactions to achieve charge selectivity, but ion translocation relies on channel hydration, facilitated by a trans-membrane field, coupled with dynamic fluctuations of the channel structure.

UNLABELLED: An automated procedure for the analysis of homologous protein structures has been developed. The method facilitates the characterization of internal conformational differences and inter-conformer relationships and provides a framework for the analysis of protein structural evolution. The method is implemented in bio3d, an R package for the exploratory analysis of structure and sequence data. AVAILABILITY: The bio3d package is distributed with full source code as a platform-independent R package under a GPL2 license from: http://mccammon.ucsd.edu/~bgrant/bio3d/

The initial coupling between ligand binding and channel gating in the human alpha7 nicotinic acetylcholine receptor (nAChR) has been investigated with targeted molecular dynamics (TMD) simulation. During the simulation, eight residues at the tip of the C-loop in two alternating subunits were forced to move toward a ligand-bound conformation as captured in the crystallographic structure of acetylcholine binding protein (AChBP) in complex with carbamoylcholine. Comparison of apo- and ligand-bound AChBP structures shows only minor rearrangements distal from the ligand-binding site. In contrast, comparison of apo and TMD simulation structures of the nAChR reveals significant changes toward the bottom of the ligand-binding domain. These structural rearrangements are subsequently translated to the pore domain, leading to a partly open channel within 4 ns of TMD simulation. Furthermore, we confirmed that two highly conserved residue pairs, one located near the ligand-binding pocket (Lys145 and Tyr188), and the other located toward the bottom of the ligand-binding domain (Arg206 and Glu45), are likely to play important roles in coupling agonist binding to channel gating. Overall, our simulations suggest that gating movements of the alpha7 receptor may involve relatively small structural changes within the ligand-binding domain, implying that the gating transition is energy-efficient and can be easily modulated by agonist binding/unbinding.

The assembly of the ribosome has recently become an interesting target for antibiotics in several bacteria. In this work, we extended an analytical procedure to determine native state fluctuations and contact breaking to investigate the protein stability dependence in the 30S small ribosomal subunit of Thermus thermophilus. We determined the causal influence of the presence and absence of proteins in the 30S complex on the binding free energies of other proteins. The predicted dependencies are in overall agreement with the experimentally determined assembly map for another organism, Escherichia coli. We found that the causal influences result from two distinct mechanisms: one is pure internal energy change, the other originates from the entropy change. We discuss the implications on how to target the ribosomal assembly most effectively by suggesting six proteins as targets for mutations or other hindering of their binding. Our results show that by blocking one out of this set of proteins, the association of other proteins is eventually reduced, thus reducing the translation efficiency even more. We could additionally determine the binding dependency of THX--a peptide not present in the ribosome of E. coli--and suggest its assembly path.

Acetylcholinesterase (AChE) rapidly hydrolyzes acetylcholine in the neuromuscular junctions and other cholinergic synapses to terminate the neuronal signal. In physiological conditions, AChE exists as tetramers associated with the proline-rich attachment domain (PRAD) of either collagen-like Q subunit (ColQ) or proline-rich membrane-anchoring protein. Crystallographic studies have revealed that different tetramer forms may be present, and it is not clear whether one or both are relevant under physiological conditions. Recently, the crystal structure of the tryptophan amphiphilic tetramerization (WAT) domain of AChE associated with PRAD ([WAT]4PRAD), which mimics the interface between ColQ and AChE tetramer, became available. In this study we built a complete tetrameric mouse [AChE(T)]4-ColQ atomic structure model, based on the crystal structure of the [WAT]4PRAD complex. The structure was optimized using energy minimization. Block normal mode analysis was done to investigate the low-frequency motions of the complex and to correlate the structure model with the two known crystal structures of AChE tetramer. Significant low-frequency motions among the catalytic domains of the four AChE subunits were observed, while the [WAT]4PRAD part held the complex together. Normal mode involvement analysis revealed that the two lowest frequency modes were primarily involved in the conformational changes leading to the two crystal structures. The first 30 normal modes can account for more than 75% of the conformational changes in both cases. The evidence further supports the idea of a flexible tetramer model for AChE. This model can be used to study the implications of the association of AChE with ColQ.

We discuss data production rates and their impact on the performance of scientific applications using parallel computers. On one hand, too high rates of data production can be overwhelming, exceeding logistical capacities for transfer, storage and analysis. On the other hand, the rate limiting step in a computationally-based study should be the human-guided analysis, not the calculation. We present performance data for a biomolecular simulation of the enzyme, acetylcholinesterase, which uses the parallel molecular dynamics program EulerGROMOS. The actual production rates are compared against a typical time frame for results analysis where we show that the rate limiting step is the simulation, and that to overcome this will require improved output rates.