MMBioS Research Highlights

2016     2015     2014    2013


2016 Highlights


Modeling with BioNetGen Gives MMBioS Team Insight into How Immune System Decides to Attack — or Let Be

T-Cell Receptor SignalingFriend or Foe?

A mix of computer modeling and laboratory experiments has helped reveal how the body differentiates “friend from foe.” Using their BioNetGen computer tool for simulating biochemistry, MMBioS members and colleagues have painted a sharper picture of how T cells, the advance scouts of the immune system, decide when to protect bodily tissues from immune attack—and when to lead the attack. The finding may guide future efforts to control human diseases like diabetes and cancer.

When T cells—a type of white blood cell—encounter other cells in the body, “they have to decide what kind of a response to make,” says James Faeder, project co-leader of MMBioS’s Technology Research and Development Project 2 Team and associate professor of computational & systems biology, University of Pittsburgh. “Is that a threat, or is it something benign? Depending on their assessment of a possible threat they can either become activated, immune-boosting cells that kill pathogens, or they can tamp down those responses.”

Read the entire article

See more about MMBioS research in cell modeling.



Tools for determining the spatial relationships between different cell components


An important task for understanding how cells are organized is determining which components have spatial patterns that are related to each other. We can think of this as learning the “order of assembly” of cells – which parts affect the patterns of other parts. We added tools for analyzing this to CellOrganizer and used them to demonstrate that different organelles that show “punctate” patterns (such as peroxisomes and coated pits) can be distinguished from each other by measuring how much their spatial patterns depend upon microtubules1. We further extended this to perform rigorous statistical testing of the extent to which these organelles depend upon four different components: the cell membrane, nuclear membrane, endoplasmic reticulum membrane, and microtubules2.

[Image Left] The panels show maps of different factors that might influence punctate organelle distribution. The factors are calculated from images of cells labeled with probes for DNA, microtubules, and endoplasmic reticulum. The values of each factor at each position in a typical cell are shown.

1Johnson GR, Buck TE, Sullivan DP, Rohde GK, Murphy RF Joint modeling of cell and nuclear shape variationMolecular biology of the cell26(22): 4046-56. doi: 10.1091/mbc.E15-06-0370.

2Li Y, Majarian TD, Naik AW, Johnson GR, Murphy RF (2016) Point process models for localization and interdependence of punctate cellular structures Cytometry Part A89:in press. doi: 10.1002/cyto.a.22873.

See more about MMBioS research in image processing


Pipeline for creation of spatiotemporal maps of T cell signaling proteins

4d rtd

Using a combination of diffeomorphic methods and improved cell segmentation, we developed a CellOrganizer pipeline for use in DPB4 to construct models of the 4D distributions of actin and 8 of its regulators during the response of T cells to antigen presentation. This work, published in Science Signaling (Roybal et al) revealed that WAVE2 and cofilin play a critical role in co-stimulation through CD28, and this was confirmed by reconstitution experiments. Spatiotemporal maps of all nine proteins were made publically available ( along with movies of the concentration of three representative proteins in 3D over time as T cell costimulation proceeds either as 3D projections.

(Figure right shows the first frame) (View full movie) or (Central 2d slices).

Frame from movie showing 4D distributions of cofilin (red), MRLC (green) and WAVE2 (blue). At this timepoint, cofilin and WAVE2 have translocated to the synapse but MRLC is lagging.

Roybal KT, Buck TE, Ruan X, Cho BH, Clark DJ, Ambler R, et al. (2016) Computational spatiotemporal analysis identifies WAVE2 and cofilin as joint regulators of costimulation-mediated T cell actin dynamics. Science Signaling. 9(424):1-13. doi: 10.1126/scisignal.aad4149.

See more about MMBioS research in image processing



New Release of iGNM Database

Gaussian network model (GNM) is a simple yet powerful model for investigating the intrinsic dynamics of proteins and their complexes. Intrinsic dynamics refers to conformational changes intrinsically favored by 3D structure, which often underlie the adaptation of biomolecules to functional interactions. Which structural elements undergo large fluctuations, or which ones provide adequate flexibility to enable conformational changes (e.g. hinge-bending sites) that may be relevant to function? Furthermore, which structural elements are subject to strongly correlated (or anticorrelated) motions, toward gaining insights into allosterically coupled regions? The GNM addresses these questions. It further allows to dissect these properties into the contributions of individual modes, thus elucidating the cooperative couplings underlied by low frequency modes.

The updated iGNM 2.0 covers more than 95% of the structures in the Protein Data Bank (PDB). Advanced search and visualization capabilities, both 2D and 3D, permit users to retrieve information on inter-residue and inter-domain cross-correlations, cooperative modes of motion, the location of hinge sites and energy localization spots. A unique feature of iGNM 2.0 is the accessibility of data for biological assemblies (BAs), rather than the single chains or asymmetric units (Asym), thus providing insights into the dynamic properties of biologically functional entities.

Examination of structures of even 104 residues in iGNM 2.0 showed that the accuracy of the results did not decrease with increasing size (N).

Li,H., Chang,YY, Yang,LW and Bahar,I. (2016) iGNM 2.0: The Gaussian Network Model Database for Biomolecular Structural Dynamics. Nucleic Acids Res., 44, D415-D422.

See more about MMBioS research in molecular modeling

Improved Sampling of Cell-Scale Models using the Weighted Ensemble Strategy

The “weighted ensemble” (WE) strategy for orchestrating a large set of parallel simulations has been established as an effective tool for efficiently calculating kinetic and equilibrium observables in molecular systems – and now has been extended to spatially resolved cell-scale systems by MMBioS researchers. In a collaboration among several MMBioS groups, the WESTPA implementation of WE has been used to control MCell simulations, including models built using a BioNetGen-CellOrganizer pipeline for situating complex biochemistry within spatially realistic cell models. As sketched, the WE strategy enables computational effort to be focused on rare, difficult-to-sample events – essentially increasing the precision of measurements in the tails of distributions. In an application to an MCell model of the frog neuro-muscular junction (NMJ), WE simulation enabled calculation of vesicle release probabilities under particularly challenging low-calcium concentrations that could not be well sampled by conventional MCell simulation: see red-circled points in the figure. This in turn, enabled validation of the model by confirmation of an established NMJ empirical power-law relationship between calcium-concentration and release probability. WE yielded estimates of observables in less overall computing time than would be required in ordinary parallelization, thus exhibiting super-linear parallel performance.

Donovan RM, Tapia JJ, Sullivan DP, Faeder JR, Murphy RF, Dittrich M, Zuckerman DM (2016). Unbiased Rare Event Sampling in Spatial Stochastic Systems Biology Models Using a Weighted Ensemble of Trajectories PLoS Comput Biol. 12(2):e1004611

See more about MMBioS research in molecular modeling


Anatomy and Function of an Excitatory Network in the Visual Cortex 

Mouse visual cortex

MMBioS researcher Greg Hood’s collaboration with Wei-Chung Allen Lee of Harvard University and R. Clay Reid of the Allen Institute for Brain Science concerning the reconstruction of an excitatory nerve-cell network in the mouse brain cortex at a subcellular level using the AlignTK software has been published in Nature.

Lee WA, Bonin V, Reed M, Graham BJ, Hood G, Glattfelder K, Reid RC (2016) Anatomy and function of an excitatory network in the visual cortex Nature 532: 370–374.

See more about MMBioS research in image processing


Memory Capacity of Brain 10 Times More than Previously Thought

brain tissue in the hippocampus

Salk researchers and MMBioS team members Tom Bartol and Terry Sejnowski, along with their collaborators, have achieved critical insight into the size of neural connections, putting the memory capacity of the brain far higher than common estimates. The new work also answers a longstanding question as to how the brain is so energy efficient and could help engineers build computers that are incredibly powerful but also conserve energy.

Read the entire story.

In a computational reconstruction of brain tissue in the hippocampus, Salk scientists and UT-Austin scientists found the unusual occurrence of two synapses from the axon of one neuron (translucent black strip) forming onto two spines on the same dendrite of a second neuron (yellow). Separate terminals from one neuron’s axon are shown in synaptic contact with two spines (arrows) on the same dendrite of a second neuron in the hippocampus. The spine head volumes, synaptic contact areas (red), neck diameters (gray) and number of presynaptic vesicles (white spheres) of these two synapses are almost identical.

Credit: Salk Institute

TM Bartol, C Bromer, J Kinney, MA Chirillo, JN Bourne, KM Harris, TJ Sejnowski (2015) Nanoconnectomic Upper Bound on the Variability of Synaptic Plasticity. eLife 2015;4:e10778.

See more about MMBioS research in cell modeling.


Development and improvements to MCell

The MCell modeling and simulation platform provides the core capabilities for spatial simulation reaction-diffusion dynamics at complex biological interfaces. Highlights of progress include:

libMCell API We have developed a high-level set of routines, which enable the creation and simulation of MCell models via API calls alone, i.e., without invoking any parsed MDL description of the model. Current functionality includes creation of species, defining reactions, defining molecule releases, creating model geometry, and simulation and runtime control.

MCell testing framework This framework, called nutmeg (, has completely replaced our previous, python based unit test framework. During the year we made nutmeg much more robust and added many additional test cases.

Implementation of new simulation capabilities that extend the range of simulations that can be performed in order to address important biological questions, such as the how membrane dynamics and signal transduction interact or how multisite phosphorylation, binding, and cooperativity affect signaling. Highlights include:

  • Dynamic geometries Users can now create dynamical meshes using Blender and use these to drive an MCell simulation of reaction and diffusion with moving boundaries.
  • Modeling pipeline Users can build complex biochemical models using BioNetGen and sample complex cellular geometries from imaging data using CellOrganizer and combine those into a single model that can be simulated in CellBlender (in collaboration with TRD3).
  • Recovery of protein-protein interactions and other aspect of biochemical mechanisms from reaction network models using Atomizer. The automated conversion of these models into a rule-based format allows for comparative analysis of models and reuse of existing model components.
  • Visualization of regulatory structurein rule-based models New visualization methods have been developed that enable global visualization of models exhibiting combinatorial complexity.

D. P. Sullivan, J. J. Tapia, R. Arepally, J. Czech, R. F. Murphy, M. Dittrich, and J. R. Faeder, “Design Automation for Biological Models : A Pipeline that Incorporates Spatial and Molecular Complexity,” in 25th edition of the Great Lakes Symposium on VLSI, 2015, pp. 321–323.

J. A. P. Sekar, J.-J. Tapia, and J. R. Faeder, “Visualizing Regulation in Rule-based Models.” Submitted to Bioinformatics . arXiv:1509.00896 [q-bio.QM].

See more about MMBioS research in cell modeling.


Development of CellBlender

CellBlender is a graphical interface for model construction, simulation, and analysis of complex spatial models of reaction-diffusion systems  Highlights of progress  include

  • A major redesign of the interface to consolidate functionality based on a single toolbar.
  • Implementation of a robust parameter handling system
  • Functionality for running simulations from within the interface taking advantage of parallel computing resources allowing users to manage on-going simulation streams.
  • Direct access to run-time error logs within the Blender interface.
  • Model and geometry importers that allow import of compartmental reaction network models in either SBML or compartmental BNGL format.
  • Implementation of internal data model provides backward compatibility with future CellBlender versions.
  • Release of CellBlender version 1.0 and publication of an article featuring MCell and CellBlender in the Encyclopedia of Computational Neuroscience

T. Bartol, M. Dittrich, and J. Faeder, “MCell,” in Encyclopedia of Computational Neuroscience, D. Jaeger and R. Jung, Eds. Springer New York, 2014, pp. 1–5.

See more about MMBioS research in cell modeling.


2015 Highlights

A Mutation in Transmembrane Domain 7 (TM7) of Excitatory Amino Acid Transporters Disrupts the Substrate-dependent Gating of the Intrinsic Anion Conductance and Drives the Channel into a Constitutively Open State


In the mammalian central nervous system, excitatory amino acid transporters (EAATs) are responsible for the clearance of glutamate after synaptic release. This energetically demanding activity is crucial for precise neuronal communication and for maintaining extracellular glutamate concentrations below neurotoxic levels. In addition to their ability to recapture glutamate from the extracellular space, EAATs exhibit a sodium and glutamate gated anion conductance. Here we show that substitution of a conserved positively charged residue (Arg388, hEAAT1) in transmembrane domain 7 with a negatively charged amino acid eliminates the ability of glutamate to further activate the anion conductance.

When expressed in oocytes, R388D or R388E mutants show large anion currents that display no further increase in amplitude after application of saturating concentrations of Na(+) and glutamate.

They also show a substantially reduced transport activity. The mutant transporters appear to exist preferentially in a sodium and glutamate independent constitutive open channel state that rarely transitions to complete the transport cycle. In addition, the accessibility of cytoplasmic residues to membrane permeant modifying reagents supports the idea that this substrate independent open state correlates with an intermediate outward facing conformation of the transporter. Our data provide additional insights into the mechanism by which substrates gate the anion conductance in EAATs and suggest that in EAAT1, Arg388 is a critical element for the structural coupling between the substrate translocation and the gating mechanisms of the EAATassociated anion channel.

TorresSalazar, D.; Jiang, J.; Divito, C. B.; GarciaOlivares, J.; Amara, S. G. A Mutation in Transmembrane Domain 7 (TM7) of Excitatory Amino Acid Transporters Disrupts the Substrate dependent Gating of the Intrinsic Anion Conductance and Drives the Channel into a Constitutiv ely Open StateJ Biol Chem 2015, 290 (38), 2297722990.


Development of diffeomorphic and dynamic models of cell and nuclear shape


We added a major new capability to our open source CellOrganizer system, the ability to construct diffeomorphic models of cell and nuclear shape. The models were developed because of the needs of DBP4; to carry out initial benchmarking we applied them to publicly available datasets and published an extensive study (Johnson et al, Mol. Biol. Cell). We measured the statistical dependency of nuclear shape on cell shape (below), and showed that drug treatment and gene inhibition can alter that dependency.

Projections of cells (right) and nuclei (left) onto the 2D subspace of highest-variance components in 7-dim diffeomorphic shape space. Cell shapes (left) are colored by the p-value that the nuclear shape predicted from the cell shape is closer to the actual shape compared to that expected by randomly choosing a nuclear shape from all examples; and the nuclear shapes on the right are colored similarly for predicting cell shape. Most cell shapes are predicted well from nuclear shape and vice versa.

Johnson GRBuck TESullivan DPRohde GKMurphy RF(2015) Joint modeling of cell and nuclear shape variationMolecular biology of the cell26(22):4046-56. doi: 10.1091/mbc.E15-06-0370. 

Amphetamine activates Rho GTPase signaling to mediate dopamine transporter internalization and acute behavioral effects of amphetamine

Acute amphetamine (AMPH) exposure elevates extracellular dopamine through a variety of mechanisms that include inhibition of dopamine reuptake, depletion of vesicular stores, and facilitation of dopamine efflux across the plasma membrane. Recent work has shown that the DAT substrate AMPH, unlike cocaine and other nontransported blockers, can also stimulate endocytosis of the plasma membrane dopamine transporter (DAT). Here, we show that when AMPH enters the cytoplasm it rapidly stimulates DAT internalization through a dynamin-dependent, clathrin-independent process. This effect, which can be observed in transfected cells, cultured dopamine neurons, and midbrain slices, is mediated by activation of the small GTPase RhoA. Inhibition of RhoA activity with C3 exotoxin or a dominant-negative RhoA blocks AMPH-induced DAT internalization. These actions depend on AMPH entry into the cell and are blocked by the DAT inhibitor cocaine. AMPH also stimulates cAMP accumulation and PKA-dependent inactivation of RhoA, thus providing a mechanism whereby PKA- and RhoA-dependent signaling pathways can interact to regulate the timing and robustness of AMPH’s effects on DAT internalization. Consistent with this model, the activation of D1/D5 receptors that couple to PKA in dopamine neurons antagonizes RhoA activation, DAT internalization, and hyperlocomotion observed in mice after AMPH treatment. These observations support the existence of an unanticipated intracellular target that mediates the effects of AMPH on RhoA and cAMP signaling and suggest new pathways to target to disrupt AMPH action.

Wheeler, D. S.; Underhill, S. M.; Stolz, D. B.; Murdoch, G. H.; Thiels, E.; Romero, G.; Amara, S. G. Amphetamine activates Rho GTPase signaling to mediate dopamine transporter internalization and acute behavioral effects of amphetamine. Proc. Natl. Acad. Sci. U. S. A112 (51), E7138-E7147.


Energy landscape for LeuT


Development of a multiscale hybrid methodology for constructing LeuT energy landscape and delineating the molecular mechanism of substrate/neurotranmitter transport cycle by NSS family members The hybrid methodology, coMD, that we have recently developed [1] has been recently extended to construct the energy landscape near the functional states of LeuT (Fig 1) [2]. This is the first energy landscape constructed for this NSS family member. It is obtained by combining 10s of microseconds of Anton trajectories [3] and coMD sampling of transition regions. It permits to visualize the accessible states and transition paths.

Fig 1. Energy landscape for LeuT, sampled during its substrate transport cycle. OF and IF refer to outward-facing and inward-facing states, and each has multiple substates (o, c, etc.).

References: 1. Gur M, Madura J, Bahar I. (2013) Biophys J 105: 164352. 2. Zomot E, Gur M, Bahar I. (2015) J Biol Chem 290: 54455. 3. Gur M, Zomot E, Cheng M, Bahar I. (2015) J Chem Phys 143: 243134


Insights into the Cooperative Dynamics of AMPAR

Interconversion from AMPAR to NMDAR conformation by moving in a space of soft modes. The dark blue bars show the overlap (correlation cosine) between ANM modes of AMPAR and the structural difference vector between AMPAR and NMDAR; the red curve displays the cumulative overlap; the dotted green curve shows the expected cumulative overlap if the modes were equally contributing each. The inset shows the initial structures of AMPAR (left) and NMDAR (right).

Ionotropic glutamate receptors (iGluRs) mediate the majority of excitatory signal transmission in the brain. Their activation underlies learning and higher-order cognitive functions. iGluRs are organized in three domain layers: The trans-membrane domain (TMD) ion channel and two domains in the extracellular region (ECR) involved in ligand sensing - the ligand binding domain (LBD) and N-terminal domain (NTD). The NTD is a powerful allosteric module in NMDAR-type iGluRs but its function in the related AMPAR-type is less clear.

We performed structural dynamics analyses including all the three domains NTD-LBD-TMD, using the crystal structures of intact AMPAR (PDB: 3KG2) and NMDAR (PDB: 4PE5), based on elastic network models (Bahar et al. Annu. Rev. Biophys. 39, 23-42, 2010). We characterized the motions within and between domain layers, identifying movements that have been related to the gating properties and confirming a prominent motion underlying AMPAR desensitization using cysteine cross-linking in vivo. Our study disclosed conserved mechanisms of motions between the two receptors, despite their substantial differences in domain packing (RMSD ~ 18 �). The figure shows one of our computations for the correlation cosines between their structural difference and each of the softest 100 ANM modes accessible to the AMPAR in order to assess the ability of the two iGluR subtypes to interconvert conformations. It is noticeable that a single ANM mode 6 yields a correlation cosine of more than 0.50 - higher than that of a random mode by a factor of 40. These data show that AMPARs possess an intrinsic potential to 'easily' (via soft modes) convert to the NMDAR structure (that is, easy compression of the AMPAR towards the more tightly packed NMDAR). Moreover, we identified hotspot residues that are strategically localized to mediate allosteric signal transmission throughout the domain layers down to the channel gate. Overall, our data provided a first glimpse into the dynamic spectrum of AMPAR and NMDAR and delineated conserved mechanisms underlying allosteric communication in iGluRs.

Go to movies of global dynamics of AMPAR and NMDAR

Dutta A, Krieger J, Lee JY,  Garcia-Nafria J, Greger IH, Bahar I. (2015) Cooperative Dynamics of Intact AMPA and NMDA Glutamate Receptors: Similarities and Subfamily-Specific Differences. Structure 23, 1692-1704.

See more about MMBioS research in molecular modeling

Molecular Mechanism of Dopamine Transport by Human Dopamine Transporter 

Dopamine transporter (DAT) controls neurotransmitter dopamine (DA) homeostasis by reuptake of excess DA from the synapse into the presynaptic neuron, assisted by the co-transport of two sodium ions. Malfunction of human DAT (hDAT) has been implicated in many neurological disorders. The first DAT structure (dDAT, from Drosophila melanogaster) has been recently resolved, which permit us to conduct a structure-based computational study of the time-resolved mechanism of DA reuptake by hDAT at full atomic scale. Using homology modeling and full-atomic microseconds accelerated simulation, we investigated the complete DA translocation including uptake from the EC region to its intracellular release and highlighted the key interactions that mediate DA translocation through hDAT. Our major observations are: spontaneous closure of extracellular gates prompted by DA binding; stabilization of a holo-occluded intermediate distinguished by close association of TM1b and TM6a with TM10 and resulting cluster of hydrophobic residues that prevents the hydration of the binding site; subsequent exposure to intracellular water triggered by Na2 dislocation, accompanied by a redistribution of salt bridges at the cytosolic surface; concerted tilting of TM5 and TM7, critical to opening the exit pore for substrate/ion channeling, enabled by the stretching of the G258-G263 loop, disruption of N82-N353 hydrogen bond, which drives the release of a Na+ ion and the Cl- ion; and DA release induced only after protonation of D79. Following Figure illustrates the transport mechanism, deduced from our over 13 sets micro-seconds MD simulations.

Proposed dopamine transport mechanism

Movie: Dopamine extracellular binding, translocation, and intracellular release, captured by MD simlutions.

Cheng MH, Bahar I (2015) Molecular Mechanism of Dopamine Transport by Human Dopamine Transporter Structure 23: 2171-81.

See more about MMBioS research in molecular modeling


Short-Term Synaptic Facilitation Revealed 

Short-term facilitation at synapses occurs during high-frequency stimulation, is known to be dependent on presynaptic calcium ions, and persists for tens of milliseconds after a presynaptic action potential. In a recent computational study, we have used the frog neuromuscular junction as a model synapse for both experimental and computer simulation studies aimed at testing various mechanistic hypotheses proposed to underlie short-term synaptic facilitation.

Our study built on our recently developed excess-calcium-binding-site model of synaptic vesicle release at the frog neuromuscular junction. Using MCell simulations, we investigated several mechanisms of short-term facilitation at the frog neuromuscular junction. Our studies place constraints on previously proposed facilitation mechanisms and conclude that the presence of a second class of calcium sensor proteins distinct from synaptotagmin can explain known properties of facilitation observed at the frog neuromuscular junction. We were further able to identify a novel facilitation mechanism, which relied on the persistent binding of calcium-bound synaptotagmin molecules to lipids of the presynaptic membrane. In a real physiological context, both mechanisms identified in our study - as well as others - may act simultaneously to give rise to the experimentally observed facilitation.

Thus, by combining computer simulations and physiological recordings we have developed a stochastic computer model of synaptic transmission at the frog neuromuscular junction, which enables microscopic insight into the biochemical dynamics of this model synapse and sheds light on its facilitation mechanism.

Ma J, Kelly L, Ingram J, Price TJ, Meriney SD, Dittrich M (2015) New insights into short-term synaptic facilitation at the frog neuromuscular junction Journal of Neurophysiology, 113:71-87.


See more about MMBioS research in cell modeling.


2014 Highlights

Learning Sequence Determinants of Protein:protein Interaction Specificity with Sparse Graphical Models 

langmead2 550In studying the strength and specificity of interaction between members of two protein families, key questions center on which pairs of possible partners actually interact, how well they interact, and why they interact while others do not. Recently, we developed a method, DgSpi (Data-driven Graphical models of Specificity in Protein:protein Interactions), for learning and using graphical models that explicitly represent the amino acid basis for interaction specificity (why) and extend earlier classification-oriented approaches (which) to predict the free energy of binding (how well). Our predicted free energy values are highly predictive of the experimentally measured ones, reaching correlation coefficients of 0.69 in 10-fold cross-validation. Furthermore, the model serves as a compact representation of amino acid constraints underlying the interactions, enabling protein-level free energy predictions to be naturally understood in terms of residue-level constraints. Finally, DgSpi readily enables the design of new interacting partners, and we demonstrated that designed ligands are novel and diverse.


Kamisetty H, Ghosh B, Langmead CJ and Bailey-Kellogg C. Learning Sequence Determinants of Protein: Protein Interaction Specificity with Sparse Graphical Models. Proceedings of the 18th Annual International Conference Research in Computational Molecular Biology (RECOMB 2014) 129-143. PMCID: PMC4235964


See more about MMBioS research in molecular modeling


Advancing Parallel Bio-simulations 

A central challenge in computer simulation is the extraction of timescales longer than the time for which biological systems can be simulated. The weighted ensemble (WE) parallel simulation approach, previously developed for a restricted set of problems, was recently advanced in MMBios work. Suarez and coworkers showed that a new non-Markovian analysis is capable of eliminating bias in estimates of long-timescale behavior, such as the mean first-passage time (MFPT) shown in the figure for the dissociation of methane molecules in explicit solvent. The analysis permits extraction of both equilibrium and non-equilibrium quantities.

Suarez E, Lettieri S, Zwier MC, Stringer CA, Subramanian SR, Chong LT, Zuckerman DM (2014) Simultaneous Computation of Dynamical and Equilibrium Information Using a Weighted Ensemble of Trajectories J. Chem. Theory Comput. 10:2658−2667.  PMCID: PMC4168800


See more about MMBioS research in molecular modeling


Stochastic Modeling and Analysis of Apoptotic Events 

stochasticModelingBing 322

Developing pharmacological strategies for controlling ionizing radiation (IR)-induced cell death is important for both mitigating radiation damage and alleviating the side effects of anti-cancer radiotherapy manifested in surrounding tissue morbidity. Exposure to IR often triggers the onset of apoptotic events mediated by the tumor suppressor protein p53, which regulates IR-induced apoptosis via both transcription-dependent and -independent pathways. In a recent study, we constructed and calibrated a stochastic model of the p53 signaling network comprised of coupled modules of nuclear p53 activation, mitochondrial cytochrome c release and cytosolic caspase activation. The model takes into account cellular heterogeneity and is able to reproduce previously published experimental observations.

Through a detailed examination of the dynamics of p53 network in response to IR damage using stochastic simulation and statistical modeling checking techniques, we found that the strength of p53 transcriptional activity and its coupling (or timing with respect) to mitochondrial pore opening are major determinants of cell fate: for systems where apoptosis is elicited via a p53-transcription-independent mechanism, direct activation of Bax by p53 becomes critical to IR-induced-damage initiation. Our results also highlight the importance of a truncated Bid/caspase-3 feedback loop in driving apoptosis and determining treatment efficacy: if the feedback is sufficiently strong, inhibition of pro-apoptotic proteins like PUMA mitigates damage; but the effect weakens with delay in inhibitor administration. Further, we show that the combined inhibition of Bid and Bax elicits an anti-apoptotic response that is effective over a range of time delays. These results offer new insights into novel polypharmacological strategies for alleviating IR damage as well as controlling cell susceptibility to apoptosis under different disease states and environmental challenges.

Liu B, Bhatt D, Oltvai ZN, Greenberger JS, Bahar I (2014) Significance of p53 dynamics in regulating apoptosis in response to ionizing radiation, and polypharmacological strategies Sci. Rep. 4: 6245


See more about MMBioS research in cell modeling.

2013 Highlights


Gating events in LeuT

Figure1 SubstrateBinding LeuT 400

Top panel: Substrate Ala (purple vDW) binding to primary site S1 prompts the closure of EC gates. (A) Salt bridge formed intermittently between R30 (blue) and D404 (red); (B) Isomerization of F253 (cyan) brings its aromatic side chain into close proximity of Y108 (violet). White sticks show the side chain orientations in the OFo crystal structure. Results are from aMD simulation of Ala binding from the EC solution.

Bottom panel: The secondary substrate-binding site S2, the functional relevance of which has been debated, appears to stably bind an alanine only when the transporter assumes an intermediate conformer close to IF state. Results are from cMD simulations of IF state.


Unraveling the molecular mechanism of function of NSS family members has been a challenge due to the involvement of both local (EC or IC gate opening/closure) and global (between outward- and inward-facing) changes in structure. These events usually occur at different time scales (e.g. tens of nanoseconds for local, microseconds or slower for global). Their examination necessitates to adoption of multiscale methods. We implemented a multiscale approach that combines conventional molecular dynamics (MD), targeted MD and accelerated MD (aMD), to investigate substrate binding events from the extracellular region. Our study shed light at the same time at a controversial issue, the role of a secondary site for binding substrate, that have drawn much attention in recent years.

Mary H. Cheng and Bahar I (2013) Coupled Global and Local Changes Direct Substrate Translocation by Neurotransmitter-Sodium Symporter Ortholog LeuT Biophys J. 105:630-639. PMID: 23931311

Anton simulations confirm the global motions predicted by ANM for GltPh

DBP1 500GltPh is the only structurally characterized member of the family of excitatory amino acid transporters. Even though various conformations of GltPh have been resolved to date and several simulations have been performed using these structural data, their mechanism of function is yet to be understood. GltPh is a homotrimer (left panel). Its monomers are composed each of two domains: a ‘scaffold’ (transmembrane TM helices 1, 2, 4 and 5) that provides support for the transport core and forms the interface between the monomers, and a transport core (TM3, 6, 7, 8 and two helical hairpins, HP1 and HP2), which binds and translocates the substrate and co-transported sodium ions. Substrate transport is enabled by alternating access to EC and IC environment, which is accomplished by the global transitions of GltPh between its outward-facing (OF) and inward-facing (IF) states, during substrate uptake and release, respectively. Experimental and computational studies have pointed to the role of HP1 and HP2 in gating the substrate-binding site. The conformations of HP1 and HP2 (bottom left panel) are determinants of the functional substates sampled during two Anton trajectories of 6 µs each generated for the IF state of GltPh. The essential motions extracted from the Anton simulations for GltPh in the IF state show close overlap with low frequency modes of Anisotropic Network Model (ANM) (right panel), a simple physics-based model of beads and springs, which exclusively depends on inter-residue contact topology.

M Gur, Zomot E, and Bahar I (2013) Global Motions Exhibited by Proteins in Micro- to Milliseconds Simulations Concur with Anisotropic Network Model. J Chem Phys 139 (12):121912. PMCID: PMC3739829

On the Intracellular Gating of Glutamate Transporters 

Sodium-coupled neurotransmitter transporters play a key role in neuronal signaling by clearing excess transmitter from the synapse. Structural data on a trimeric archaeal aspartate transporter, GltPh, have provided valuable insights into structural features of human excitatory amino acid transporters. However, the time-resolved mechanisms of substrate binding and release, as well as that of coupling to sodium co-transport, remain largely unknown for this important family.

Our recent study highlights the role of the helical hairpin HP2 as an intracellular gate, in addition to its role as an extracellular gate. In a recent study, we have elucidated the pathway of aspartate release from the inward-facing structure to the intracellular medium by generating multiple microsecond trajectories (using the Anton supercomputer), which consistently show that the helical hairpin HP2, not HP1, serves as an intracellular gate (in addition to its extracellular gating role). In contrast to previous proposals, HP1 can neither initiate nor accommodate neurotransmitter release without prior opening of HP2 by at least 4.0 Å. Aspartate release invariably follows that of a sodium ion located near the HP2 gate entrance. Asp394 on TM8 and Arg276 on HP1 emerge as key residues that promote the reorientation and diffusion of substrate toward the cell interior.

Zomot E and Bahar I (2013) Intracellular Gating in an Inward-facing State of Aspartate Transporter GltPh Is Regulated by the Movements of the Helical Hairpin HP2  J. Biol. Chem. 288:8231-37


See more about MMBioS research in molecular modeling and in modeling neurotransmitter transport.


Junction Crossing Ahead 

azSynapses are the connections between neurons in the brain or between neurons and muscle fibers - so called neuromuscular junctions - and thus play a crucial role in human health. Despite decades of intense experimental study we still lack a detailed understanding of how synapses work.

Many neuro-degenerative diseases such as dementia, depression and anxiety can be traced to improperly functioning synapses. Using a combination of computer simulations and experimental study, we have recently developed a spatially-realistic model of an active zone in the neuromuscular junction of the frog [1]. Using MCell's Monte Carlo algorithms to simulate the model, we could predict the calcium binding stoichiometry and dynamics that underlie
neurotransmitter release. Our model revealed that 20-40 independent calcium binding sites on synaptic vesicles, only a fraction of which need to bind calcium to trigger fusion, are sufficient to predict physiological release.

Our excess-calcium-binding site model has many functional advantages, agrees with recent data on synaptotagmin copy number, and is the first to link detailed physiological observations with the molecular machinery of calcium-triggered exocytosis. Importantly, our model provides detailed microscopic insight into the underlying calcium dynamics during synapse activation and thus constitutes a stepping stone for more detailed investigations in the future.

Dittrich M, Pattillo JM, King JD, Cho S, Stiles JR , Meriney SD (2013) An Excess-Calcium-Binding-Site Model Predicts Neurotransmitter Release at the Neuromuscular Junction Biophysical Journal 104:2751-2763.


See more about MMBioS research in cell modeling.


Stochastic Simulation 

fig1 weightedensemble biggreySystems biology, the quantitative study of complex interacting biological systems, is becoming increasingly demanding of computational resources. As models become able to capture truly physiological behavior, some of the most difficult computations may also be the most important, such as a rare transition from normal to pathological behavior.

A novel approach to meet growing computational demands is described in a paper recently accepted to the Journal of Chemical Physics ("Efficient Stochastic Simulation of Chemical Kinetics Networks Using A Weighted Ensemble Of Trajectories"). First-author Donovan, in collaboration with MMBioS investigators Faeder and Zuckerman, applied a sophisticated parallel simulation algorithm originally developed for small-scale molecular systems to systems biology models.

The results were eye-catching: the new approach was more computationally efficient at modeling rare events than standard methods by orders of magnitude, even for a very complex model with thousands of reacting species. 

Donovan RM, Sedgewick AJ, Faeder JR, Zuckerman DM (2013) Efficient Stochastic Simulation of Chemical Kinetics Networks Using A Weighted Ensemble Of Trajectories J. Chem. Phys. 139:115105.


See more about MMBioS research in molecular modeling.


CoMD: a Hybrid Methodology for Allosteric Changes

MMBios BPJ cover


(Above) Three step (LID-NMP-LID) transition mechanism of adenylate kinase. Residues whose interactions undergo notable redistributions are shown in space-filling and labeled.

Efficient and accurate mapping of transition pathways is a challenging problem in allosteric proteins. We propose here a new methodology, called collective molecular dynamics (coMD) that takes advantage of the collective modes of motions encoded by the fold, while evaluating the interactions and energetics via a full-atomic molecular dynamics simulation protocol. The basic approach of the coMD is to deform the structure collectively along the modes predicted by the anisotropic network model (ANM)1 upon selecting them via a Monte Carlo/Metropolis algorithm from amongst the complete pool of all accessible modes. The method consists of a series of cycles, denoted by index k. Pairs of intermediate conformers are simultaneously generated starting from both ends at each cycle k.  Each cycle is composed of multiple collectivemoves, s = 0, 1, 2,…, stot(k) where all Ca atoms are simultaneously displaced, guided by the ANM. The Ca-coordinates of the stot(k)th conformer are used to guide the full-atomic targeted molecular dynamics simulations (TMD), which followed by as short energy minimization (EM).

Application to adenylate kinase (AK), an allosteric enzyme composed of three domains, CORE, LID and NMP, shows that both open-to-closed and closed-to-open transitions of AK are readily sampled by coMD, being dominated by large-scale motions of the LID. An energy-barrier crossing occurs during the NMP movements. The energy barrier originates from a switch between the salt bridges K136- D118 at LID-CORE interface and K57-E170 and D33-R156 at CORE-NMP and LID-NMP interfaces, respectively. Despite its simplicity and computing efficiency, coMD yields ensembles of transition pathways in close accord with detailed full atomic simulations, lending support to its utility as a multiscale hybrid method for efficiently exploring the allosteric transitions of multidomain or multimeric proteins.

See more about MMBioS research in molecular modeling.


Micro-to-milliseconds Simulations vs Elastic Network Models 

JCP cover 1c

(Above) Structures representative of substates S1-S8 (colored) superimposed onto the crystal structure (grey). Transitions observed during simulations are indicated by arrows.

The Anton supercomputing technology recently developed for efficient molecular dynamics (MD) simulations permits us to examine micro- to millisecond events at full atomic resolution for proteins in explicit water and lipid bilayer. It also permits us to investigate to what extent the collective motions predicted by network models (that have found broad use in molecular biophysics) agree with those exhibited by full-atomic long simulations. The present study focuses on Anton trajectories generated for two systems: the bovine pancreatic trypsin inhibitor (BPTI), and an archaeal aspartate transporter, GltPh.  The former, a thoroughly studied system, helps benchmark the method of comparative analysis, and the latter provides new insights into the mechanism of function of glutamate transporters. The principal components (PCs) of motion derived from the 1.013 ms long MD trajectory of  BPTI closely overlapped with those predicted by the anisotropic network model (ANM)1, which a simple physics-based model of beads and springs that exclusively depends on inter-residue contact topology. 

Notably, the ANM modes define the collective mechanisms, or the pathways on conformational energy landscape, that underlie the passage between the crystal structure and substates visited in simulations. In particular, the lowest frequency ANM modes facilitate the conversion between the most probable substates, lending support to the view that easy access to functional substates is a robust determinant of evolutionarily selected native contact topology.

Similar to BPTI, the PCs of the 12 µs long MD trajectory of the GltPh homotrimer3 exhibited good overlap with its ANM modes. Moreover, these ANM modes appear to naturally facilitate the transitions between the substates sampled in the simulations.

Gur M, Zomot E, Bahar I (2013) Global motions exhibited by proteins in micro- to milliseconds simulations concur with anisotropic network model predictions J Chem Phys 139(12): 121912.


See more about MMBioS research in molecular modeling.


ANMPathway for Exploring Transitions

Biomolecular conformational transitions are essential to biological functions. Most experimental methods report on the long-lived functional states of biomolecules, but information about the transition pathways between these stable states is generally scarce. Such transitions involve short-lived conformational states that are difficult to detect experimentally.

For this reason, computational methods are needed to produce plausible hypothetical transition pathways that can then be probed experimentally. Here we propose a simple and computationally efficient method, called ANMPathway, for constructing a physically reasonable pathway between two endpoints of a conformational transition. We adopt a coarse-grained representation of the protein and construct a two-state potential by combining two anisotropic network models  (ANMs)1 representative of the experimental structures resolved for the endpoints.

MMBios Plos Cover

(Above)  Conformational transition between states with two protomers facing inward and all protomers facing inward of the glutamate transporter, GltPh. (A). Left : Structure of the intermediate state (iOF) state where two protomers are in IF and the third (green) protomer is in the OF; right : IF state where all protomers are in IF. The central diagram (stick representation) is the Cα trace of the transition state produced by ANMPathway. (B) The energy of the system along the transition. Total number of images in the pathway is 79, RMSD between two consecutive images is ~0.1Å. (C) Cumulative squared cosines between the structural change to reach a few selected conformers/images along the transition pathway and the ANM modes accessible to the starting (iOF)  structure.  The force constants and cut-offs for both the end states were set to 0.1 kcal/mol and 15 Å respectively. No energy offsets were used for either of the end  structures.

The two-state potential has a cusp in the configuration space where the energies from both the ENMs are equal. We first search for the minimum energy structure along the cusp hyper-surface and then treat it as the transition state. The continuous pathway is subsequently constructed by following the steepest descent energy minimization trajectories starting from the transition state on each side of the cusp. Application to several systems of broad biological interest such as adenylate kinase, ATP-driven calcium pump SERCA, leucine transporter and glutamate transporter (See figure at left) shows that ANMPathway yields results in good agreement with those from other similar methods and with data obtained from all-atom molecular dynamics simulations, in support of the utility of this simple and efficient approach. Notably the method provides experimentally testable predictions, including the formation of non-native contacts during the transition which we were able to detect in two of the systems we studied.



See more about MMBioS research in molecular modeling.


Major New Release of CellOrganizer 2.0 Published 


(Above) Example cell image generated by CellOrganizer showing the nuclear membrane (red), cell boundary (blue) and individual lysosomes (yellow).

A major new release of the CellOrganizer system for creating image-derived models of cell shape and organization has just been published. The software is a major focus of research supported by the NIH through the National Center for Multiscale Modeling of Biological Systems (MMBioS). It creates statistical, generative models of cell and nuclear shape, microtubule patterns, and vesicular organelles that can be used as the basis for cell simulations (another focus of MMBioS). Generative models capture not just a description of a pattern but the ability to produce new examples of that pattern (analogously to the way in which humans capture models of letters or spoken words not just by describing them but by learning to produce new examples). Support for CellOrganizer has also come from NIH grants GM075205 and GM090033. Specific improvements or additions in the new release include:

Cell image synthesis

  • Improved vesicular organelles model
    • Eliminate object/object and object/boundary overlap during generation
  • Ability to combine models learned from images of different resolution, and synthesize images at desired resolution
  • Ability to synthesis random walks in shape space from diffeomorphic models of cell and nuclear shape
    • Includes directed random walks using Willmore energy and shape space density

Model training

  • New capabilities for cell and nuclear shape model learning
    • Build nuclear models from images without nuclear marker
    • Build joint diffeomorphic models of cell and nuclear shape
  • Per-cell representations for easy model building and comparison


  • Export to SBML-spatial and mesh formats for interfacing with tools such as CellBlender and VCell
  • Parallelization of model learning pipeline
  • Reporter tools
    • Assess learned models
    • Compare models
    • Compare per-cell parameters within or across models


See more about MMBioS research in image processing and analysis.


Release of AlignTK Image Alignment Software 

mvc small

A cutaway view of a 450x350x52μm volume reconstructed from ~1150 45nm sections of mouse visual cortex. (Missing sections are left black.) Closeups of the x-z and y-z faces show that neural features can easily be traced in a direction perpendicular to the original cutting planes. 

We have issued the first public release of the AlignTK toolkit for image alignment. This package is primarily targeted at large 3-D datasets that are acquired in sections, such as in serial-section electron microscopy. The sectioning process physically distorts each tissue section in a slightly different way, and imaging these sections typically introduces additional electron-beam and optical distortions. Aligning these sections to produce a coherent 3-D volume is challenging, and AlignTK provides a set of tools that can be used in a batch environment to process millions of camera images into a volumetric dataset in which axons and dendrites of individual neurons can be traced.

An earlier version of these tools was used at NRBSC and the Pittsburgh Supercomputing Center to align the mouse visual cortex dataset in the Nature paper by Bock et al. These tools have also been applied to optical images of stained mouse brain sections and time-series fluorescence images of the frog neuromuscular junction.

Bock, DD, Lee, WA, Kerlin, AM, Andermann, ML, Hood, G, Wetzel, AW, Yurgenson, S, Soucy, ER, Kim, HS, & Reid, RC (2011). Network anatomy and in vivo physiology of visual cortical neurons, Nature, 471:177-182.


See more about MMBioS research in image processing and analysis.