@article {Chowdhury2021.05.17.444532, title = {Empowering drug off-target discovery with metabolic and structural analysis}, year = {Submitted}, publisher = {Cold Spring Harbor Laboratory}, abstract = {Understanding intracellular antibiotic targeting and the associated mechanisms leading to bacterial growth inhibition has been a difficult problem. Here, we discovered the additional intracellular targets of the novel {\textquotedblleft}evolution-drug{\textquotedblright} lead CD15-3 designed to delay the emergence of antibiotic resistance by inhibiting bacterial DHFR and its Trimethoprim resistant variants. Overexpression of DHFR only partially rescued inhibition of E. coli growth by CD15.3 suggesting that CD15.3 also inhibits a non-DHFR target in the cell. We utilized untargeted global metabolomics and the metabolic network analysis along with structural similarity search of the putative targets to identify the additional target of CD15-3. We validated in vivo and in vitro that besides DHFR CD15-3 inhibits HPPK (folK), an essential protein upstream of DHFR in bacterial folate metabolism. This bivalent cellular targeting makes CD15-3 a promising lead to develop a {\textquotedblleft}monotherapy analogue{\textquotedblright} of combination drugs.Competing Interest StatementThe authors have declared no competing interest.}, doi = {10.1101/2021.05.17.444532}, url = {https://www.biorxiv.org/content/early/2021/05/18/2021.05.17.444532}, author = {Chowdhury, Sourav and Zielinski, Daniel C. and Dalldorf, Christopher and Rodrigues, Joao V and Palsson, Bernhard O. and Shakhnovich, EugeneI.} } @article {serebryany2020native, title = {A Native Chemical Chaperone in the Human Eye Lens}, journal = {arXiv:2012.09805 [q-bio]}, year = {Submitted}, month = {dec}, abstract = {Cataract is one of the most prevalent protein aggregation disorders and still the biggest cause of vision loss worldwide. The human lens, in its core region, lacks turnover of any cells or cellular components; it has therefore evolved remarkable mechanisms for resisting protein aggregation for a lifetime. We now report that one such mechanism relies on an unusually abundant metabolite, myo-inositol, to suppress light-scattering aggregation of lens proteins. We quantified aggregation suppression by in vitro turbidimetry and characterized both macroscopic and microscopic mechanisms of myo-inositol action using negative-stain electron microscopy, differential scanning fluorometry, and a thermal scanning Raman spectroscopy apparatus. Given recent metabolomic evidence that it is dramatically depleted in human cataractous lenses compared to age-matched controls, we suggest that maintaining or restoring healthy levels of myo-inositol in the lens may be a simple, safe, and widely available strategy for reducing the global burden of cataract.}, keywords = {Quantitative Biology - Biomolecules}, url = {https://arxiv.org/abs/2012.09805}, author = {Serebryany, Eugene and Chowdhury, Sourav and Watson, Nicki E. and McClelland, Arthur and Shakhnovich, EugeneI.} } @article {wang2019improved, title = {Improved fragment-based movement with LRFragLib for all-atom Ab initio protein folding}, journal = {arXiv preprint arXiv:1906.05785}, year = {Submitted}, abstract = {Fragment-based assembly has been widely used in Ab initio protein folding simulation which can effectively reduce the conformational space and thus accelerate sampling. The efficiency of fragment-based movement as well as the quality of fragment library determine whether the folding process can lead the free energy landscape to the global minimum and help the protein to reach near-native folded state. We designed an improved fragment-based movement, "fragmove", which substituted multiple backbone dihedral angles in every simulation step. This movement strategy was derived from the fragment library generated by LRFragLib, an effective fragment detection algorithm using logistic regression model. We show in replica exchange Monte Carlo (REMC) simulation that "fragmove", when compared with a set of existing movements in REMC, shows significant improved ability at increasing secondary and tertiary predicted model accuracy by 11.24\% and 17.98\%, respectively and reaching energy minima decreased by 5.72\%. Our results demonstrate that this improved movement is more powerful to guide proteins faster to low energy regions of conformational space and promote folding efficiency and predicted model accuracy.}, url = {https://arxiv.org/abs/1906.05785}, author = {Wang, Tong and Gong, Haipeng and Shakhnovich, EugeneI.} } @article {razban2021avoidancea, title = {Avoidance of Protein Unfolding Constrains Protein Stability in Long-Term Evolution}, journal = {Biophysical Journal}, volume = {120}, year = {2021}, pages = {2413{\textendash}2424}, abstract = {Every amino acid residue can influence a protein{\textquoteright}s overall stability, making stability highly susceptible to change throughout evolution. We consider the distribution of protein stabilities evolutionarily permittable under two previously reported protein fitness functions: flux dynamics and misfolding avoidance. We develop an evolutionary dynamics theory and find that it agrees better with an extensive protein stability data set for dihydrofolate reductase orthologs under the misfolding avoidance fitness function rather than the flux dynamics fitness function. Further investigation with ribonuclease H data demonstrates that not any misfolded state is avoided; rather, it is only the unfolded state. At the end, we discuss how our work pertains to the universal protein abundance-evolutionary rate correlation seen across organisms{\textquoteright} proteomes. We derive a closed-form expression relating protein abundance to evolutionary rate that captures Escherichia coli, Saccharomyces cerevisiae, and Homo sapiens experimental trends without fitted parameters.}, issn = {0006-3495}, doi = {10.1016/j.bpj.2021.03.042}, url = {https://www.sciencedirect.com/science/article/pii/S0006349521003404}, author = {Rostam M. Razban and Pouria Dasmeh and Adrian W.R. Serohijos and Shakhnovich, EugeneI.} } @article {zhang2021development, title = {Development of Antibacterial Compounds That Constrain Evolutionary Pathways to Resistance}, journal = {eLife}, volume = {10}, year = {2021}, pages = {e64518}, publisher = {eLife Sciences Publications, Ltd}, abstract = {Antibiotic resistance is a worldwide challenge. A potential approach to block resistance is to simultaneously inhibit WT and known escape variants of the target bacterial protein. Here, we applied an integrated computational and experimental approach to discover compounds that inhibit both WT and trimethoprim (TMP) resistant mutants of E. coli dihydrofolate reductase (DHFR). We identified a novel compound (CD15-3) that inhibits WT DHFR and its TMP resistant variants L28R, P21L and A26T with IC50 50{\textendash} 75 \textmu M against WT and TMP-resistant strains. Resistance to CD15-3 was dramatically delayed compared to TMP in in vitro evolution. Whole genome sequencing of CD15-3-resistant strains showed no mutations in the target folA locus. Rather, gene duplication of several efflux pumps gave rise to weak (about twofold increase in IC50) resistance against CD15-3. Altogether, our results demonstrate the promise of strategy to develop evolution drugs - compounds which constrain evolutionary escape routes in pathogens.}, keywords = {antibiotic resistance, computational modelling, dihydrofolate reductase, evolution drug}, issn = {2050-084X}, doi = {10.7554/eLife.64518}, url = {https://elifesciences.org/articles/64518}, author = {Zhang, Yanmin and Chowdhury, Sourav and Rodrigues, Jo{\~a}o V and Eugene Shakhnovich}, editor = {Landry, Christian R and Perry, George H} } @article {zhao2021switching, title = {Switching an Active Site Helix in Dihydrofolate Reductase Reveals Limits to Subdomain Modularity}, journal = {Biophysical Journal}, year = {2021}, publisher = {Elsevier}, abstract = {To what degree are individual structural elements within proteins modular such that similar structures from unrelated proteins can be interchanged? We study subdomain modularity by creating 20 chimeras of an enzyme,\ Escherichia colidihydrofolate reductase (DHFR), in which a catalytically important, 10-residue\ α-helical sequence is replaced by\ α-helical sequences from a diverse set of proteins. The chimeras stably fold but have a range of diminished thermal stabilities and catalytic activities. Evolutionary coupling analysis indicates that the residues of this\ α-helix are under selection pressure to maintain catalytic activity in DHFR. Reversion to phenylalanine at key position 31 was found to partially restore catalytic activity, which could be explained by evolutionary coupling values. We performed molecular dynamics simulations using replica exchange with solute tempering. Chimeras with low catalytic activity exhibit nonhelical conformations that block the binding site and disrupt the positioning of the catalytically essential residue D27. Simulation observables and in\ vitro measurements of thermal stability and substrate-binding affinity are strongly correlated. Several\ E.\ coli\ strains with chromosomally integrated chimeric DHFRs can grow, with growth rates that follow predictions from a kinetic flux model that depends on the intracellular abundance and catalytic activity of DHFR. Our findings show that although\ α-helices are not universally substitutable, the molecular and fitness effects of modular segments can be predicted by the biophysical compatibility of the replacement segment.}, issn = {0006-3495}, doi = {10.1016/j.bpj.2021.09.032}, url = {https://www.cell.com/biophysj/fulltext/S0006-3495(21)00794-3}, author = {Zhao, Victor Y. and Rodrigues, Jo{\~a}o V. and Elena R. Lozovsky and Daniel L. Hartl and Shakhnovich, EugeneI.} } @article {graff2021accelerating, title = {Accelerating High-Throughput Virtual Screening through Molecular Pool-Based Active Learning}, journal = {Chemical Science}, volume = {12}, year = {2021}, pages = {7866{\textendash}7881}, publisher = {The Royal Society of Chemistry}, abstract = {Structure-based virtual screening is an important tool in early stage drug discovery that scores the interactions between a target protein and candidate ligands. As virtual libraries continue to grow (in excess of 108 molecules), so too do the resources necessary to conduct exhaustive virtual screening campaigns on these libraries. However, Bayesian optimization techniques, previously employed in other scientific discovery problems, can aid in their exploration: a surrogate structure{\textendash} property relationship model trained on the predicted affinities of a subset of the library can be applied to the remaining library members, allowing the least promising compounds to be excluded from evaluation. In this study, we explore the application of these techniques to computational docking datasets and assess the impact of surrogate model architecture, acquisition function, and acquisition batch size on optimization performance. We observe significant reductions in computational costs; for example, using a directed-message passing neural network we can identify 94.8\% or 89.3\% of the top-50 000 ligands in a 100M member library after testing only 2.4\% of candidate ligands using an upper confidence bound or greedy acquisition strategy, respectively. Such model-guided searches mitigate the increasing computational costs of screening increasingly large virtual libraries and can accelerate high-throughput virtual screening campaigns with applications beyond docking.}, issn = {2041-6539}, doi = {10.1039/D0SC06805E}, url = {https://pubs.rsc.org/en/content/articlelanding/2021/SC/D0SC06805E}, author = {Graff, David E. and Shakhnovich, EugeneI. and Coley, Connor W.} } @article {bhattacharyya2021metabolic, title = {Metabolic Response to Point Mutations Reveals Principles of Modulation of in Vivo Enzyme Activity and Phenotype}, journal = {Molecular Systems Biology}, volume = {17}, year = {2021}, pages = {e10200}, publisher = {John Wiley \& Sons, Ltd}, abstract = {The relationship between sequence variation and phenotype is poorly understood. Here, we use metabolomic analysis to elucidate the molecular mechanism underlying the filamentous phenotype of\ E.\ coli\ strains that carry destabilizing mutations in dihydrofolate reductase (DHFR). We find that partial loss of DHFR activity causes reversible filamentation despite SOS response indicative of DNA damage, in contrast to thymineless death (TLD) achieved by complete inhibition of DHFR activity by high concentrations of antibiotic trimethoprim. This phenotype is triggered by a disproportionate drop in intracellular dTTP, which could not be explained by drop in dTMP based on the Michaelis{\textendash}Menten-like\ in vitro\ activity curve of thymidylate kinase (Tmk), a downstream enzyme that phosphorylates dTMP to dTDP. Instead, we show that a highly cooperative (Hill coefficient 2.5)\ in vivo\ activity of Tmk is the cause of suboptimal dTTP levels. dTMP supplementation rescues filamentation and restores\ in vivo\ Tmk kinetics to Michaelis{\textendash}Menten. Overall, this study highlights the important role of cellular environment in sculpting enzymatic kinetics with system-level implications for bacterial phenotype.}, keywords = {filamentation, fractal kinetics, in vivo enzyme activity, metabolomics, thymine limitation}, issn = {1744-4292}, doi = {10.15252/msb.202110200}, author = {Bhattacharyya, Sanchari and Shimon Bershtein and Adkar, Bharat V. and Jaie C. Woodard and Shakhnovich, EugeneI.} } @article {ranganathan2021effect, title = {Effect of RNA on Morphology and Dynamics of Membraneless Organelles}, journal = {The Journal of Physical Chemistry B}, volume = {125}, year = {2021}, pages = {5035{\textendash}5044}, publisher = {American Chemical Society}, abstract = {Membraneless organelles (MLOs) are spatiotemporally regulated structures that concentrate multivalent proteins or RNA, often in response to stress. The proteins enriched within MLOs are often classified as high-valency {\textquoteleft}{\textquoteleft}scaffolds{\textquoteright}{\textquoteright} or low-valency {\textquoteleft}{\textquoteleft}clients{\textquoteright}{\textquoteright}, with the former being associated with a phase-separation promoting role. In this study, we employ a minimal model for P-body components, with a defined protein{\textendash} protein interaction network, to study their phase separation at biologically realistic low protein concentrations. Without RNA, multivalent proteins can assemble into solid-like clusters only in the regime of high concentration and stable interactions. RNA molecules promote cluster formation in an RNA-length-dependent manner, even in the regime of weak interactions and low protein volume fraction. Our simulations reveal that long RNA chains act as superscaffolds that stabilize large RNA{\textendash} protein clusters by recruiting low-valency proteins within them while also ensuring functional {\textquoteleft}{\textquoteleft}liquid-like{\textquoteright}{\textquoteright} turnover of components. Our results suggest that RNA-mediated phase separation could be a plausible mechanism for spatiotemporally regulated phase separation in the cell.}, issn = {1520-6106}, doi = {10.1021/acs.jpcb.1c02286}, url = {https://pubs.acs.org/doi/abs/10.1021/acs.jpcb.1c02286}, author = {Ranganathan, Srivastav and Eugene Shakhnovich} } @article {bitran2020validation, title = {Validation of DBFOLD: An Efficient Algorithm for Computing Folding Pathways of Complex Proteins}, journal = {PLOS Computational Biology}, volume = {16}, year = {2020}, pages = {e1008323}, publisher = {Public Library of Science}, abstract = {Atomistic simulations can provide valuable, experimentally-verifiable insights into protein folding mechanisms, but existing ab initio simulation methods are restricted to only the smallest proteins due to severe computational speed limits. The folding of larger proteins has been studied using native-centric potential functions, but such models omit the potentially crucial role of non-native interactions. Here, we present an algorithm, entitled DBFOLD, which can predict folding pathways for a wide range of proteins while accounting for the effects of non-native contacts. In addition, DBFOLD can predict the relative rates of different transitions within a protein{\textquoteright}s folding pathway. To accomplish this, rather than directly simulating folding, our method combines equilibrium Monte-Carlo simulations, which deploy enhanced sampling, with unfolding simulations at high temperatures. We show that under certain conditions, trajectories from these two types of simulations can be jointly analyzed to compute unknown folding rates from detailed balance. This requires inferring free energies from the equilibrium simulations, and extrapolating transition rates from the unfolding simulations to lower, physiologically-reasonable temperatures at which the native state is marginally stable. As a proof of principle, we show that our method can accurately predict folding pathways and Monte-Carlo rates for the well-characterized Streptococcal protein G. We then show that our method significantly reduces the amount of computation time required to compute the folding pathways of large, misfolding-prone proteins that lie beyond the reach of existing direct simulation. Our algorithm, which is available online, can generate detailed atomistic models of protein folding mechanisms while shedding light on the role of non-native intermediates which may crucially affect organismal fitness and are frequently implicated in disease.}, keywords = {Biochemical simulations, Computing methods, Free energy, Melting, Monte Carlo Method, protein folding, Protein structure, Simulation and modeling}, issn = {1553-7358}, doi = {10.1371/journal.pcbi.1008323}, url = {https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008323}, author = {Bitran, Amir and Jacobs, William M. and Eugene Shakhnovich} } @article {zhao2020effect, title = {Effect of Protein Structure on Evolution of Cotranslational Folding}, journal = {Biophysical Journal}, volume = {119}, year = {2020}, pages = {1123{\textendash}1134}, publisher = {Elsevier}, issn = {0006-3495}, doi = {10.1016/j.bpj.2020.06.037}, url = {https://www.cell.com/biophysj/fulltext/S0006-3495(20)30606-8}, author = {Zhao, Victor and Jacobs, William M. and Shakhnovich, EugeneI.} } @article {ranganathan2020dynamic, title = {Dynamic metastable long-living droplets formed by sticker-spacer proteins}, journal = {Elife}, volume = {9}, year = {2020}, pages = {e56159}, publisher = {eLife Sciences Publications Limited}, abstract = {Multivalent biopolymers phase separate into membrane-less organelles (MLOs) which exhibit liquid-like behavior. Here, we explore formation of prototypical MOs from multivalent proteins on various time and length scales and show that the kinetically arrested metastable multi-droplet state is a dynamic outcome of the interplay between two competing processes: a diffusion-limited encounter between proteins, and the exhaustion of available valencies within smaller clusters. Clusters with satisfied valencies cannot coalesce readily, resulting in metastable, long-living droplets. In the regime of dense clusters akin to phase-separation, we observe co-existing assemblies, in contrast to the single, large equilibrium-like cluster. A system-spanning network encompassing all multivalent proteins was only observed at high concentrations and large interaction valencies. In the regime favoring large clusters, we observe a slow-down in the dynamics of the condensed phase, potentially resulting in loss of function. Therefore, metastability could be a hallmark of dynamic functional droplets formed by sticker-spacer proteins.}, url = {https://elifesciences.org/articles/56159}, author = {Ranganathan, Srivastav and Shakhnovich, EugeneI.} } @article {razban2020effects, title = {Effects of single mutations on protein stability are Gaussian distributed}, journal = {Biophysical Journal}, year = {2020}, publisher = {Elsevier}, abstract = {The distribution of protein stability effects is known to be well approximated by a Gaussian distribution from previous empirical fits. Starting from first-principles statistical mechanics, we more rigorously motivate this empirical observation by deriving per-residue-position protein stability effects to be Gaussian. Our derivation requires the number of amino acids to be large, which is satisfied by the standard set of 20 amino acids found in nature. No assumption is needed on the number of residues in close proximity in space, in contrast to previous applications of the central limit theorem to protein energetics. We support our derivation results with computational and experimental data on mutant protein stabilities across all types of protein residues.}, url = {https://www.sciencedirect.com/science/article/pii/S0006349520303714}, author = {Razban, Rostam M and Shakhnovich, EugeneI.} } @article {Bitran201913207, title = {Cotranslational folding allows misfolding-prone proteins to circumvent deep kinetic traps}, journal = {Proceedings of the National Academy of Sciences}, year = {2020}, publisher = {National Academy of Sciences}, abstract = {Many proteins must adopt a specific structure to perform their functions, and failure to do so has been linked to disease. Although small proteins often fold rapidly and spontaneously to their native conformations, larger proteins are less likely to fold correctly due to the myriad incorrect arrangements they can adopt. Here, we provide mechanistic insights into how this problem can be alleviated if proteins start folding while they are being translated by the ribosome. This process of cotranslational folding biases certain proteins away from misfolded states that tend to hinder spontaneous refolding. Signatures of unusually slow translation suggest that some of these proteins have evolved to fold cotranslationally.Many large proteins suffer from slow or inefficient folding in vitro. It has long been known that this problem can be alleviated in vivo if proteins start folding cotranslationally. However, the molecular mechanisms underlying this improvement have not been well established. To address this question, we use an all-atom simulation-based algorithm to compute the folding properties of various large protein domains as a function of nascent chain length. We find that for certain proteins, there exists a narrow window of lengths that confers both thermodynamic stability and fast folding kinetics. Beyond these lengths, folding is drastically slowed by nonnative interactions involving C-terminal residues. Thus, cotranslational folding is predicted to be beneficial because it allows proteins to take advantage of this optimal window of lengths and thus avoid kinetic traps. Interestingly, many of these proteins{\textquoteright} sequences contain conserved rare codons that may slow down synthesis at this optimal window, suggesting that synthesis rates may be evolutionarily tuned to optimize folding. Using kinetic modeling, we show that under certain conditions, such a slowdown indeed improves cotranslational folding efficiency by giving these nascent chains more time to fold. In contrast, other proteins are predicted not to benefit from cotranslational folding due to a lack of significant nonnative interactions, and indeed these proteins{\textquoteright} sequences lack conserved C-terminal rare codons. Together, these results shed light on the factors that promote proper protein folding in the cell and how biomolecular self-assembly may be optimized evolutionarily.}, issn = {0027-8424}, doi = {10.1073/pnas.1913207117}, url = {https://www.pnas.org/content/early/2020/01/06/1913207117}, author = {Bitran, Amir and Jacobs, William M. and Zhai, Xiadi and Eugene Shakhnovich} } @article {chen2019semi, title = {Semi-rational design and molecular dynamics simulations study of the thermostability enhancement of cellobiose 2-epimerases}, journal = {International Journal of Biological Macromolecules}, year = {2020}, publisher = {Elsevier}, abstract = {Directed evolution using random mutation in vast sequence space leads to the low probability of obtaining target proteins. Emerging engineering strategies with computational tools are developed for more trustable outcomes. We used some semi-rational design methods to modify an industrial enzyme, namely cellobiose 2-epimerase (CE). A mutant was selected for its better thermostability and isomerization activity. The tradeoffs between thermostability, epimerization activity and isomerization activity of the CE mutants were different. To investigate the computational prediction performance of protein stability upon point mutations, molecular dynamics (MD) simulation analyses were conducted. The root mean square deviation (RMSD) and hydrogen bond analyses reproduced the correct trends in stability changes of the wild-type and mutated CEs with relatively high accuracy (correlation coefficients r\ ~\ 0.5{\textendash}0.8). The simulation temperature and time are important factors that influence the prediction performance. Our result shows that thermostability predictors calculated from MD simulation do better in predicting the thermostability changes of the mutated enzymes than the predictors using static-state information of the enzymes.}, url = {https://www.sciencedirect.com/science/article/pii/S0141813019371636}, author = {Chen, Qiuming and Xiao, Yaqin and Shakhnovich, EugeneI. and Zhang, Wenli and Wanmeng Mu} } @article {zhou2019common, title = {Common Activation Mechanism of Class A GPCRs}, journal = {eLife}, volume = {8}, year = {2019}, pages = {e50279}, publisher = {eLife Sciences Publications, Ltd}, abstract = {Class A G-protein-coupled receptors (GPCRs) influence virtually every aspect of human physiology. Understanding receptor activation mechanism is critical for discovering novel therapeutics since about one-third of all marketed drugs target members of this family. GPCR activation is an allosteric process that couples agonist binding to G-protein recruitment, with the hallmark outward movement of transmembrane helix 6 (TM6). However, what leads to TM6 movement and the key residue level changes of this movement remain less well understood. Here, we report a framework to quantify conformational changes. By analyzing the conformational changes in 234 structures from 45 class A GPCRs, we discovered a common GPCR activation pathway comprising of 34 residue pairs and 35 residues. The pathway unifies previous findings into a common activation mechanism and strings together the scattered key motifs such as CWxP, DRY, Na+ pocket, NPxxY and PIF, thereby directly linking the bottom of ligand-binding pocket with G-protein coupling region. Site-directed mutagenesis experiments support this proposition and reveal that rational mutations of residues in this pathway can be used to obtain receptors that are constitutively active or inactive. The common activation pathway provides the mechanistic interpretation of constitutively activating, inactivating and disease mutations. As a module responsible for activation, the common pathway allows for decoupling of the evolution of the ligand binding site and G-protein-binding region. Such an architecture might have facilitated GPCRs to emerge as a highly successful family of proteins for signal transduction in nature.}, keywords = {activation mechanism, allostery, Drug Discovery, genetic diseases, GPCR, Signal Transduction}, issn = {2050-084X}, doi = {10.7554/eLife.50279}, url = {https://elifesciences.org/articles/50279}, author = {Zhou, Qingtong and Yang, Dehua and Wu, Meng and Guo, Yu and Guo, Wanjing and Zhong, Li and Cai, Xiaoqing and Dai, Antao and Jang, Wonjo and Shakhnovich, EugeneI. and Liu, Zhi-Jie and Stevens, Raymond C and Lambert, Nevin A and Babu, M Madan and Wang, Ming-Wei and Zhao, Suwen}, editor = {Shan, Yibing and Boudker, Olga} } @article {razban2019protein, title = {Protein melting temperature cannot fully assess whether protein folding free energy underlies the universal abundance{\textendash}evolutionary rate correlation seen in proteins}, journal = {Molecular Biology and Evolution}, year = {2019}, abstract = {The protein misfolding avoidance hypothesis explains the universal negative correlation between protein abundance and sequence evolutionary rate across the proteome by identifying protein folding free energy (ΔG) as the confounding variable. Abundant proteins resist toxic misfolding events by being more stable, and more stable proteins evolve slower because their mutations are more destabilizing. Direct supporting evidence consists only of computer simulations. A study taking advantage of a recent experimental breakthrough in measuring protein stability proteome-wide through melting temperature (Tm) (Leuenberger et al. 2017), found weak misfolding avoidance hypothesis support for the\ Escherichia coli\ proteome, and no support for the\ Saccharomyces cerevisiae,\ Homo sapiens, and\ Thermus thermophilus\ proteomes (Plata and Vitkup 2018). I find that the nontrivial relationship between\ Tm\ and ΔG\ and inaccuracy in\ Tm\ measurements by Leuenberger et al. 2017 can be responsible for not observing strong positive abundance{\textendash}Tm\ and strong negative\ Tm{\textendash}evolutionary rate correlations.}, url = {https://academic.oup.com/mbe/article/36/9/1955/5489915}, author = {Razban, Rostam M} } @article {rodrigues2019adaptation, title = {Adaptation to mutational inactivation of an essential gene converges to an accessible suboptimal fitness peak}, journal = {eLife}, volume = {8}, year = {2019}, publisher = {eLife Sciences Publications, Ltd}, abstract = {The mechanisms of adaptation to inactivation of essential genes remain unknown. Here we inactivate\ E. coli\ dihydrofolate reductase (DHFR) by introducing D27G,N,F chromosomal mutations in a key catalytic residue with subsequent adaptation by an automated serial transfer protocol. The partial reversal G27- \> C occurred in three evolutionary trajectories. Conversely, in one trajectory for D27G and in all trajectories for D27F,N strains adapted to grow at very low metabolic supplement (folAmix) concentrations but did not escape entirely from supplement auxotrophy. Major global shifts in metabolome and proteome occurred upon DHFR inactivation, which were partially reversed in adapted strains. Loss-of-function mutations in two genes,\ thyA\ and\ deoB, ensured adaptation to low folAmix by rerouting the 2-Deoxy-D-ribose-phosphate metabolism from glycolysis towards synthesis of dTMP. Multiple evolutionary pathways of adaptation converged to a suboptimal solution due to the high accessibility to loss-of-function mutations that block the path to the highest, yet least accessible, fitness peak.}, url = {https://elifesciences.org/articles/50509}, author = {Rodrigues, Jo{\~a}o V and Shakhnovich, EugeneI.} } @article {Rodrigues_2019-2, title = {Chimeric dihydrofolate reductases display properties of modularity and biophysical diversity}, journal = {Protein Science}, year = {2019}, month = {May}, publisher = {Wiley}, abstract = {While reverse genetics and functional genomics have long affirmed the role of individual mutations in determining protein function, there have been fewer studies addressing how large-scale changes in protein sequences, such as in entire modular segments, influence protein function and evolution. Given how recombination can reassort protein sequences, these types of changes may play an underappreciated role in how novel protein functions evolve in nature. Such studies could aid our understanding of whether certain organismal phenotypes related to protein function{\textemdash}such as growth in the presence or absence of an antibiotic{\textemdash}are robust with respect to the identity of certain modular segments. In this study, we combine molecular genetics with biochemical and biophysical methods to gain a better understanding of protein modularity in dihydrofolate reductase (DHFR), an enzyme target of antibiotics also widely used as a model for protein evolution. We replace an integral α-helical segment of\ Escherichia coliDHFR with segments from a number of different organisms (many nonmicrobial) and examine how these chimeric enzymes affect organismal phenotypes (e.g., resistance to an antibiotic) as well as biophysical properties of the enzyme (e.g., thermostability). We find that organismal phenotypes and enzyme properties are highly sensitive to the identity of DHFR modules, and that this chimeric approach can create enzymes with diverse biophysical characteristics.}, issn = {1469-896X}, doi = {10.1002/pro.3646}, url = {https://onlinelibrary.wiley.com/doi/full/10.1002/pro.3646}, author = {Rodrigues, Jo{\~a}o V. and Ogbunugafor, C. Brandon and Daniel L. Hartl and Shakhnovich, EugeneI.} } @article {J_S_Loureiro_2019, title = {The Early Phase of β2m Aggregation: An Integrative Computational Study Framed on the D76N Mutant and the ΔN6 Variant}, journal = {Biomolecules}, volume = {9}, year = {2019}, month = {Aug}, pages = {366}, publisher = {MDPI AG}, abstract = {

Human β2-microglobulin (b2m) protein is classically associated with dialysis-related amyloidosis (DRA). Recently, the single point mutant D76N was identified as the causative agent of a hereditary systemic amyloidosis affecting visceral organs. To get insight into the early stage of the β2m aggregation mechanism, we used molecular simulations to perform an in depth comparative analysis of the dimerization phase of the D76N mutant and the ΔN6 variant, a cleaved form lacking the first six N-terminal residues, which is a major component of ex vivo amyloid plaques from DRA patients. We also provide first glimpses into the tetramerization phase of D76N at physiological pH. Results from extensive protein{\textendash}protein docking simulations predict an essential role of the C- and N-terminal regions (both variants), as well as of the BC-loop (ΔN6 variant), DE-loop (both variants) and EF-loop (D76N mutant) in dimerization. The terminal regions are more relevant under acidic conditions while the BC-, DE- and EF-loops gain importance at physiological pH. Our results recapitulate experimental evidence according to which Tyr10 (A-strand), Phe30 and His31 (BC-loop), Trp60 and Phe62 (DE-loop) and Arg97 (C-terminus) act as dimerization hot-spots, and further predict the occurrence of novel residues with the ability to nucleate dimerization, namely Lys-75 (EF-loop) and Trp-95 (C-terminus). We propose that D76N tetramerization is mainly driven by the self-association of dimers via the N-terminus and DE-loop, and identify Arg3 (N-terminus), Tyr10, Phe56 (D-strand) and Trp60 as potential tetramerization hot-spots.

}, issn = {2218-273X}, doi = {10.3390/biom9080366}, url = {http://dx.doi.org/10.3390/biom9080366}, author = {J. S. Loureiro, Rui and Vila-Vi{\c c}osa, Diogo and Machuqueiro, Miguel and Shakhnovich, EugeneI. and F. N. Fa{\'\i}sca, Patr{\'\i}cia} } @article {Xiao_2019, title = {Simulation-guided enzyme discovery: A new microbial source of cellobiose 2-epimerase}, journal = {International Journal of Biological Macromolecules}, volume = {139}, year = {2019}, month = {Oct}, pages = {1002{\textendash}1008}, publisher = {Elsevier BV}, abstract = {Cellobiose 2-epimerase (CE) is a promising industrial enzyme that can be utilized in the dairy industry. More thermostable CEs from different microorganisms are still needed for a higher lactulose productivity. This study demonstrated the feasibility to use molecular dynamics (MD) simulation as the preliminary computational filter for thermostable enzymes screening. Sequence information of eleven uncharacterized CEs were chosen to be analyzed by MD simulations. The CE from\ Dictyoglomus thermophilum\ (Dith-CE) was determined experimentally to be one of the most thermostable CEs with the highest epimerization (160 {\textpm} 6.5 U mg-1) and isomerization activities (3.52 {\textpm} 0.23 U mg-1) among all the reported CEs. This enzyme shows the highest isomerization activity at 85 {\textdegree}C and pH 7.0. The kinetic parameters (kcat\ and\ Km) of isomerization activity of this CE are 3.98 {\textpm} 0.3 s-1\ and 235.2 {\textpm} 11.2 mM, respectively. These results suggest that the CE from Dith-CE is a promising lactulose-producing enzyme.}, issn = {0141-8130}, doi = {10.1016/j.ijbiomac.2019.08.075}, url = {http://dx.doi.org/10.1016/j.ijbiomac.2019.08.075}, author = {Xiao, Yaqin and Chen, Qiuming and Shakhnovich, EugeneI. and Zhang, Wenli and Wanmeng Mu} } @article {Adkar201821447, title = {Substrate inhibition imposes fitness penalty at high protein stability}, journal = {Proceedings of the National Academy of Sciences}, year = {2019}, publisher = {National Academy of Sciences}, abstract = {Why highly stable proteins are not commonly observed in mesophilic organisms is an important evolutionary question. It has been suggested that high stability has a fitness cost in the form of loss of activity. However, this hypothesis has not yet been experimentally demonstrated. Here, we use an essential bacterial enzyme adenylate kinase (Adk) to explore this hypothesis and show that, as Adk{\textquoteright}s stability increases, one of its own substrates inhibits its activity. Bacterial strains carrying such stable Adks show substantial fitness defects, which can be mapped to the loss of Adk activity due to substrate inhibition. Overall, our study adds {\textquotedblleft}substrate inhibition{\textquotedblright} to the {\textquotedblleft}toolbox{\textquotedblright} that is used to rationalize the stability distribution of proteins.Proteins are only moderately stable. It has long been debated whether this narrow range of stabilities is solely a result of neutral drift toward lower stability or purifying selection against excess stability{\textemdash}for which no experimental evidence was found so far{\textemdash}is also at work. Here, we show that mutations outside the active site in the essential Escherichia coli enzyme adenylate kinase (Adk) result in a stability-dependent increase in substrate inhibition by AMP, thereby impairing overall enzyme activity at high stability. Such inhibition caused substantial fitness defects not only in the presence of excess substrate but also under physiological conditions. In the latter case, substrate inhibition caused differential accumulation of AMP in the stationary phase for the inhibition-prone mutants. Furthermore, we show that changes in flux through Adk could accurately describe the variation in fitness effects. Taken together, these data suggest that selection against substrate inhibition and hence excess stability may be an important factor determining stability observed for modern-day Adk.}, issn = {0027-8424}, doi = {10.1073/pnas.1821447116}, url = {https://www.pnas.org/content/early/2019/05/14/1821447116}, author = {Adkar, Bharat V. and Bhattacharyya, Sanchari and Amy I. Gilson and Zhang, Wenli and Shakhnovich, EugeneI.} } @article {PhysRevE.99.022424, title = {Mutation rate variability as a driving force in adaptive evolution}, journal = {Phys. Rev. E}, volume = {99}, year = {2019}, month = {Feb}, pages = {022424}, publisher = {American Physical Society}, abstract = {Mutation rate is a key determinant of the pace as well as outcome of evolution, and variability in this rate has been shown in different scenarios to play a key role in evolutionary adaptation and resistance evolution under stress caused by selective pressure. Here we investigate the dynamics of resistance fixation in a bacterial population with variable mutation rates, and we show that evolutionary outcomes are most sensitive to mutation rate variations when the population is subject to environmental and demographic conditions that suppress the evolutionary advantage of high-fitness subpopulations. By directly mapping a biophysical fitness function to the system-level dynamics of the population, we show that both low and very high, but not intermediate, levels of stress in the form of an antibiotic result in a disproportionate effect of hypermutation on resistance fixation. We demonstrate how this behavior is directly tied to the extent of genetic hitchhiking in the system, the propagation of high-mutation rate cells through association with high-fitness mutations. Our results indicate a substantial role for mutation rate flexibility in the evolution of antibiotic resistance under conditions that present a weak advantage over wildtype to resistant cells.}, doi = {10.1103/PhysRevE.99.022424}, url = {https://link.aps.org/doi/10.1103/PhysRevE.99.022424}, author = {Engelhardt, Dalit and Shakhnovich, EugeneI.} } @article {doi:10.1021/acs.jpcb.8b05842, title = {Accurate Protein-Folding Transition-Path Statistics from a Simple Free-Energy Landscape}, journal = {The Journal of Physical Chemistry B}, volume = {122}, year = {2018}, note = {PMID: 30091592}, pages = {11126-11136}, abstract = {A central goal of protein-folding theory is to predict the stochastic dynamics of transition paths{\textemdash}the rare trajectories that transit between the folded and unfolded ensembles{\textemdash}using only thermodynamic information, such as a low-dimensional equilibrium free-energy landscape. However, commonly used one-dimensional landscapes typically fall short of this aim, because an empirical coordinate-dependent diffusion coefficient has to be fit to transition-path trajectory data in order to reproduce the transition-path dynamics. We show that an alternative, first-principles free-energy landscape predicts transition-path statistics that agree well with simulations and single-molecule experiments without requiring dynamical data as an input. This {\textquotedblleft}topological configuration{\textquotedblright} model assumes that distinct, native-like substructures assemble on a time scale that is slower than native-contact formation but faster than the folding of the entire protein. Using only equilibrium simulation data to determine the free energies of these coarse-grained intermediate states, we predict a broad distribution of transition-path transit times that agrees well with the transition-path durations observed in simulations. We further show that both the distribution of finite-time displacements on a one-dimensional order parameter and the ensemble of transition-path trajectories generated by the model are consistent with the simulated transition paths. These results indicate that a landscape based on transient folding intermediates, which are often hidden by one-dimensional projections, can form the basis of a predictive model of protein-folding transition-path dynamics.}, doi = {10.1021/acs.jpcb.8b05842}, url = {https://doi.org/10.1021/acs.jpcb.8b05842}, author = {Jacobs, William M. and Shakhnovich, EugeneI.} } @article {cite-key, title = {Growth tradeoffs produce complex microbial communities on a single limiting resource}, journal = {Nature Communications}, volume = {9}, year = {2018}, pages = {3214}, abstract = {The relationship between the dynamics of a community and its constituent pairwise interactions is a fundamental problem in ecology. Higher-order ecological effects beyond pairwise interactions may be key to complex ecosystems, but mechanisms to produce these effects remain poorly understood. Here we model microbial growth and competition to show that higher-order effects can arise from variation in multiple microbial growth traits, such as lag times and growth rates, on a single limiting resource with no other interactions. These effects produce a range of ecological phenomena: an unlimited number of strains can exhibit multistability and neutral coexistence, potentially with a single keystone strain; strains that coexist in pairs do not coexist all together; and a strain that wins all pairwise competitions can go extinct in a mixed competition. Since variation in multiple growth traits is ubiquitous in microbial populations, our results indicate these higher-order effects may also be widespread, especially in laboratory ecology and evolution experiments.}, isbn = {2041-1723}, doi = {10.1038/s41467-018-05703-6}, url = {https://doi.org/10.1038/s41467-018-05703-6}, author = {Manhart, Michael and Shakhnovich, EugeneI.} } @article {serebryany2018dynamic, title = {Dynamic disulfide exchange in a crystallin protein in the human eye lens promotes cataract-associated aggregation}, journal = {Journal of Biological Chemistry}, volume = {293}, year = {2018}, pages = {17997{\textendash}18009}, publisher = {ASBMB}, abstract = {Increased light scattering in the eye lens due to aggregation of the long-lived lens proteins, crystallins, is the cause of cataract disease. Several mutations in the gene encoding human γD-crystallin (HγD) cause misfolding and aggregation. Cataract-associated substitutions at Trp42\ cause the protein to aggregate\ in vitro\ from a partially unfolded intermediate locked by an internal disulfide bridge, and proteomic evidence suggests a similar aggregation precursor is involved in age-onset cataract. Surprisingly, WT HγD can promote aggregation of the W42Q variant while itself remaining soluble. Here, a search for a biochemical mechanism for this interaction has revealed a previously unknown oxidoreductase activity in HγD. Using\ in vitrooxidation, mutational analysis, cysteine labeling, and MS, we have assigned this activity to a redox-active internal disulfide bond that is dynamically exchanged among HγD molecules. The W42Q variant acts as a disulfide sink, reducing oxidized WT and forming a distinct internal disulfide that kinetically traps the aggregation-prone intermediate. Our findings suggest a redox {\textquotedblleft}hot potato{\textquotedblright} competition among WT and mutant or modified polypeptides wherein variants with the lowest kinetic stability are trapped in aggregation-prone intermediate states upon accepting disulfides from more stable variants. Such reactions may occur in other long-lived proteins that function in oxidizing environments. In these cases, aggregation may be forestalled by inhibiting disulfide flow toward mutant or damaged polypeptides.}, url = {http://www.jbc.org/content/293/46/17997.abstract}, author = {Serebryany, Eugene and Yu, Shuhuai and Trauger, Sunia A and Budnik, Bogdan and Shakhnovich, EugeneI.} } @article {bhattacharyya2018accessibility, title = {Accessibility of the Shine-Dalgarno sequence dictates N-terminal codon bias in E. coli}, journal = {Molecular cell}, volume = {70}, year = {2018}, pages = {894{\textendash}905}, publisher = {Elsevier}, abstract = {Despite considerable efforts, no physical mechanism has been shown to explain\ N-terminalcodon bias\ in\ prokaryotic\ genomes. Using a systematic study of\ synonymous substitutions\ in two endogenous\ E.\ coli\ genes, we show that interactions between the\ coding region\ and the upstream\ Shine-Dalgarno\ (SD) sequence modulate the efficiency of translation initiation, affecting both intracellular mRNA and protein levels due to the inherent coupling of\ transcription\ and translation in\ E.\ coli. We further demonstrate that far-downstream mutations can also modulate mRNA levels by occluding the SD sequence through the formation of non-equilibrium secondary structures. By contrast, a non-endogenous\ RNA polymerase\ that decouples transcription and translation largely alleviates the effects of synonymous substitutions on mRNA levels. Finally, a complementary statistical analysis of the\ E.\ coli\ genome specifically implicates avoidance of intra-molecular base pairing with the SD sequence. Our results provide general physical insights into the coding-level features that optimize\ protein expression\ in prokaryotes.}, url = {https://www.sciencedirect.com/science/article/pii/S1097276518303587}, author = {Bhattacharyya, Sanchari and Jacobs, William M and Adkar, Bharat V and Jin Yan and Zhang, Wenli and Shakhnovich, EugeneI.} } @article {rotem2018evolution, title = {Evolution on the biophysical fitness landscape of an RNA virus}, journal = {Molecular biology and evolution}, year = {2018}, abstract = {Viral evolutionary pathways are determined by the fitness landscape, which maps viral genotype to fitness. However, a quantitative description of the landscape and the evolutionary forces on it remain elusive. Here, we apply a biophysical fitness model based on capsid folding stability and antibody binding affinity to predict the evolutionary pathway of norovirus escaping a neutralizing antibody. The model is validated by experimental evolution in bulk culture and in a drop-based microfluidics that propagates millions of independent small viral subpopulations. We demonstrate that along the axis of binding affinity, selection for escape variants and drift due to random mutations have the same direction, an atypical case in evolution. However, along folding stability, selection and drift are opposing forces whose balance is tuned by viral population size. Our results demonstrate that predictable epistatic tradeoffs between molecular traits of viral proteins shape viral evolution.}, url = {https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msy131/5046249?searchresult=1}, author = {Rotem, Assaf and Adrian W.R. Serohijos and Chang, Connie B and Wolfe, Joshua T and Fischer, Audrey E and Mehoke, Thomas S and Zhang, Huidan and Tao, Ye and Ung, W Lloyd and Jeong-Mo Choi and others} } @article {razban2018proteomevis, title = {ProteomeVis: a web app for exploration of protein properties from structure to sequence evolution across organisms{\textquoteright} proteomes}, journal = {Bioinformatics}, volume = {1}, year = {2018}, pages = {9}, publisher = {Oxford University Press}, abstract = {

Motivation:\ Protein evolution spans time scales and its effects span the length of an organism. A web app named ProteomeVis is developed to provide a comprehensive view of protein evolution in the\ Saccharomyces cerevisiae\ and\ Escherichia coli\ proteomes. ProteomeVis interactively creates protein chain graphs, where edges between nodes represent structure and sequence similarities within user-defined ranges, to study the long time scale effects of protein structure evolution. The short time scale effects of protein sequence evolution are studied by sequence evolutionary rate (ER) correlation analyses with protein properties that span from the molecular to the organismal level.

Results:\ We demonstrate the utility and versatility of ProteomeVis by investigating the distribution of edges per node in organismal protein chain universe graphs (oPCUGs) and putative ER determi- nants.\ S.cerevisiae\ and\ E.coli\ oPCUGs are scale-free with scaling constants of 1.79 and 1.56, respectively. Both scaling constants can be explained by a previously reported theoretical model describing protein structure evolution. Protein abundance most strongly correlates with ER among properties in ProteomeVis, with Spearman correlations of {\textendash}0.49 (P-value\ \<\ 10-10) and {\textendash}0.46 (P-value\<10-10) for\ S.cerevisiae\ and\ E.coli, respectively. This result is consistent with previous reports that found protein expression to be the most important ER determinant.

}, url = {https://academic.oup.com/bioinformatics/advance-article-pdf/doi/10.1093/bioinformatics/bty370/25018003/bty370.pdf}, author = {Razban, Rostam M and Gilson, Amy I and Durfee, Niamh and Strobelt, Hendrik and Dinkla, Kasper and Jeong-Mo Choi and Hanspeter Pfister and Shakhnovich, EugeneI.} } @article {klein2018stability, title = {Stability of the Influenza Virus Hemagglutinin Protein Correlates with Evolutionary Dynamics}, journal = {mSphere}, volume = {3}, year = {2018}, pages = {e00554{\textendash}17}, publisher = {Am Soc Microbiol}, abstract = {Protein thermodynamics are an integral determinant of viral fitness and one of the major drivers of protein evolution. Mutations in the influenza A virus (IAV) hemagglutinin (HA) protein can eliminate neutralizing antibody binding to mediate escape from preexisting antiviral immunity. Prior research on the IAV nucleoprotein suggests that protein stability may constrain seasonal IAV evolution; however, the role of stability in shaping the evolutionary dynamics of the HA protein has not been explored. We used the full coding sequence of 9,797 H1N1pdm09 HA sequences and 16,716 human seasonal H3N2 HA sequences to computationally estimate relative changes in the thermal stability of the HA protein between 2009 and 2016. Phylogenetic methods were used to characterize how stability differences impacted the evolutionary dynamics of the virus. We found that pandemic H1N1 IAV strains split into two lineages that had different relative HA protein stabilities and that later variants were descended from the higher-stability lineage. Analysis of the mutations associated with the selective sweep of the higher-stability lineage found that they were characterized by the early appearance of highly stabilizing mutations, the earliest of which was not located in a known antigenic site. Experimental evidence further suggested that H1N1 HA stability may be correlated with\ in vitro\ virus production and infection. A similar analysis of H3N2 strains found that surviving lineages were also largely descended from viruses predicted to encode more-stable HA proteins. Our results suggest that HA protein stability likely plays a significant role in the persistence of different IAV lineages.}, url = {http://msphere.asm.org/content/3/1/e00554-17}, author = {Eili Y. Klein and Blumenkrantz, Deena and Serohijos, Adrian and Eugene Shakhnovich and Jeong-Mo Choi and Rodrigues, Jo{\~a}o V and Smith, Brendan D and Lane, Andrew P and Feldman, Andrew and Andrew Pekosz} } @article {manhart2018trade, title = {Trade-offs between microbial growth phases lead to frequency-dependent and non-transitive selection}, journal = {Proc. R. Soc. B}, volume = {285}, year = {2018}, pages = {20172459}, publisher = {The Royal Society}, abstract = {Mutations in a microbial population can increase the frequency of a genotype not only by increasing its exponential growth rate, but also by decreasing its lag time or adjusting the yield (resource efficiency). The contribution of multiple life-history traits to selection is a critical question for evolutionary biology as we seek to predict the evolutionary fates of mutations. Here we use a model of microbial growth to show that there are two distinct components of selection corresponding to the growth and lag phases, while the yield modulates their relative importance. The model predicts rich population dynamics when there are trade-offs between phases: multiple strains can coexist or exhibit bistability due to frequency-dependent selection, and strains can engage in rock{\textendash}paper{\textendash}scissors interactions due to non-transitive selection. We characterize the environmental conditions and patterns of traits necessary to realize these phenomena, which we show to be readily accessible to experiments. Our results provide a theoretical framework for analysing high-throughput measurements of microbial growth traits, especially interpreting the pleiotropy and correlations between traits across mutants. This work also highlights the need for more comprehensive measurements of selection in simple microbial systems, where the concept of an ordinary fitness landscape breaks down.}, url = {http://rspb.royalsocietypublishing.org/content/285/1872/20172459}, author = {Manhart, Michael and Adkar, Bharat V and Shakhnovich, EugeneI.} } @article {valleau2017absence, title = {Absence of Selection for Quantum Coherence in the Fenna{\textendash}Matthews{\textendash}Olson Complex: A Combined Evolutionary and Excitonic Study}, journal = {ACS central science}, volume = {3}, year = {2017}, pages = {1086{\textendash}1095}, publisher = {ACS Publications}, abstract = {We present a study on the evolution of the Fenna{\textendash}Matthews{\textendash}Olson bacterial photosynthetic pigment{\textendash}protein complex. This protein complex functions as an antenna. It transports absorbed photons{\textemdash}excitons{\textemdash}to a reaction center where photosynthetic reactions initiate. The efficiency of exciton transport is therefore fundamental for the photosynthetic bacterium{\textquoteright}s survival. We have reconstructed an ancestor of the complex to establish whether coherence in the exciton transport was selected for or optimized over time. We have also investigated the role of optimizing free energy variation upon folding in evolution. We studied whether mutations which connect the ancestor to current day species were stabilizing or destabilizing from a thermodynamic viewpoint. From this study, we established that most of these mutations were thermodynamically neutral. Furthermore, we did not see a large change in exciton transport efficiency or coherence, and thus our results predict that exciton coherence was not specifically selected for.}, url = {https://pubs.acs.org/doi/full/10.1021/acscentsci.7b00269}, author = {Valleau, St{\'e}phanie and Studer, Romain A and H{\"a}se, Florian and Kreisbeck, Christoph and Saer, Rafael G and Blankenship, Robert E and Shakhnovich, EugeneI. and Aspuru-Guzik, Al{\'a}n} } @article {jacobs2017evidence, title = {Evidence of evolutionary selection for cotranslational folding}, journal = {Proceedings of the National Academy of Sciences}, volume = {114}, year = {2017}, pages = {11434{\textendash}11439}, publisher = {National Acad Sciences}, abstract = {Recent experiments and simulations have demonstrated that proteins can fold on the ribosome. However, the extent and generality of fitness effects resulting from cotranslational folding remain open questions. Here we report a genome-wide analysis that uncovers evidence of evolutionary selection for cotranslational folding. We describe a robust statistical approach to identify loci within genes that are both significantly enriched in slowly translated codons and evolutionarily conserved. Surprisingly, we find that domain boundaries can explain only a small fraction of these conserved loci. Instead, we propose that regions enriched in slowly translated codons are associated with cotranslational folding intermediates, which may be smaller than a single domain. We show that the intermediates predicted by a native-centric model of cotranslational folding account for the majority of these loci across more than 500\ Escherichia coli\ proteins. By making a direct connection to protein folding, this analysis provides strong evidence that many synonymous substitutions have been selected to optimize translation rates at specific locations within genes. More generally, our results indicate that kinetics, and not just thermodynamics, can significantly alter the efficiency of self-assembly in a biological context.}, url = {http://www.pnas.org/content/114/43/11434.short}, author = {Jacobs, William M and Shakhnovich, EugeneI.} } @article {williams2017statistical, title = {Statistical physics of the symmetric group}, journal = {Physical Review E}, volume = {95}, year = {2017}, pages = {042126}, publisher = {APS}, abstract = {Ordered chains (such as chains of amino acids) are ubiquitous in biological cells, and these chains perform specific functions contingent on the sequence of their components. Using the existence and general properties of such sequences as a theoretical motivation, we study the statistical physics of systems whose state space is defined by the possible permutations of an ordered list, i.e., the symmetric group, and whose energy is a function of how certain permutations deviate from some chosen correct ordering. Such a nonfactorizable state space is quite different from the state spaces typically considered in statistical physics systems and consequently has novel behavior in systems with interacting and even noninteracting Hamiltonians. Various parameter choices of a mean-field model reveal the system to contain five different physical regimes defined by two transition temperatures, a triple point, and a quadruple point. Finally, we conclude by discussing how the general analysis can be extended to state spaces with more complex combinatorial properties and to other standard questions of statistical mechanics models.}, url = {https://journals.aps.org/pre/abstract/10.1103/PhysRevE.95.042126}, author = {Williams, Mobolaji} } @article {cheron2017effect, title = {Effect of sampling on BACE-1 ligands binding free energy predictions via MM-PBSA calculations}, journal = {Journal of Computational Chemistry}, year = {2017}, publisher = {Wiley Online Library}, abstract = {The BACE-1 enzyme is a prime target to find a cure to Alzheimer{\textquoteright}s disease. In this article, we used the MM-PBSA approach to compute the binding free energies of 46 reported ligands to this enzyme. After showing that the most probable protonation state of the catalytic dyad is mono-protonated (on ASP32), we performed a thorough analysis of the parameters influencing the sampling of the conformational space (in total, more than 35 μs of simulations were performed). We show that ten simulations of 2 ns gives better results than one of 50 ns. We also investigated the influence of the protein force field, the water model, the periodic boundary conditions artifacts (box size), as well as the ionic strength. Amber03 with TIP3P, a minimal distance of 1.0 nm between the protein and the box edges and a ionic strength of\ I = 0.2 M provides the optimal correlation with experiments. Overall, when using these parameters, a Pearson correlation coefficient of\ R = 0.84 (R2 = 0.71) is obtained for the 46 ligands, spanning eight orders of magnitude of\ Kd\ (from 0.017 nm to 2000 μM, i.e., from -14.7 to -3.7 kcal/mol), with a ligand size from 22 to 136 atoms (from 138 to 937 g/mol). After a two-parameter fit of the binding affinities for 12 of the ligands, an error of RMSD = 1.7 kcal/mol was obtained for the remaining ligands.}, url = {http://onlinelibrary.wiley.com/doi/10.1002/jcc.24839/full}, author = {Ch{\'e}ron, Nicolas and Shakhnovich, EugeneI.} } @article {zhou2017exploring, title = {Exploring the Mutational Robustness of Nucleic Acids by Searching Genotype Neighborhoods in Sequence Space}, journal = {The Journal of Physical Chemistry Letters}, volume = {8}, year = {2017}, note = {PMID: 28045264}, pages = {407-414}, abstract = {To assess the mutational robustness of nucleic acids, many genome- and protein-level studies have been performed, where nucleic acids are treated as genetic information carriers and transferrers. However, the molecular mechanisms through which mutations alter the structural, dynamic, and functional properties of nucleic acids are poorly understood. Here we performed a SELEX in silico study to investigate the fitness distribution of the\ l-Arm-binding aptamer genotype neighborhoods. Two novel functional genotype neighborhoods were isolated and experimentally verified to have comparable fitness as the wild-type. The experimental aptamer fitness landscape suggests the mutational robustness is strongly influenced by the local base environment and ligand-binding mode, whereas bases distant from the binding pocket provide potential evolutionary pathways to approach the global fitness maximum. Our work provides an example of successful application of SELEX in silico to optimize an aptamer and demonstrates the strong sensitivity of mutational robustness to the site of genetic variation.}, doi = {10.1021/acs.jpclett.6b02769}, url = {http://pubs.acs.org/doi/abs/10.1021/acs.jpclett.6b02769}, author = {Zhou, Qingtong and Sun, Xianbao and Xia, Xiaole and Fan, Zhou and Luo, Zhaofeng and Zhao, Suwen and Eugene Shakhnovich and Liang, Haojun} } @article {choi2017graph, title = {Graph{\textquoteright}s Topology and Free Energy of a Spin Model on the Graph}, journal = {Physical Review Letters}, volume = {118}, year = {2017}, pages = {088302}, publisher = {APS}, abstract = {In this Letter we investigate a direct relationship between a graph{\textquoteright}s topology and the free energy of a spin system on the graph. We develop a method of separating topological and energetic contributions to the free energy, and find that considering the topology is sufficient to qualitatively compare the free energies of different graph systems at high temperature, even when the energetics are not fully known. This method was applied to the metal lattice system with defects, and we found that it partially explains why point defects are more stable than high-dimensional defects. Given the energetics, we can even quantitatively compare free energies of different graph structures via a closed form of linear graph contributions. The closed form is applied to predict the sequence-space free energy of lattice proteins, which is a key factor determining the designability of a protein structure.}, url = {https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.118.088302}, author = {Jeong-Mo Choi and Gilson, Amy I and Shakhnovich, EugeneI.} } @article {debroise2017hybrid, title = {A Hybrid Knowledge-Based and Empirical Scoring Function for Protein{\textendash}Ligand Interaction: SMoG2016}, journal = {Journal of Chemical Information and Modeling}, volume = {57}, year = {2017}, pages = {584{\textendash}593}, publisher = {ACS Publications}, abstract = {We present the third generation of our scoring function for the prediction of protein{\textendash}ligand binding free energy. This function is now a hybrid between a knowledge-based potential and an empirical function. We constructed a diversified set of \~{}1000 complexes from the PDBBinding-CN database for the training of the function, and we show that this number of complexes generates enough data to build the potential. The occurrence of 420 different types of atomic pairwise interactions is computed in up to five different ranges of distances to derive the knowledge-based part. All of the parameters were optimized, and we were able to considerably improve the accuracy of the scoring function with a Pearson correlation coefficient against experimental binding free energies of up to 0.57, which ranks our new scoring function as one of the best currently available and the second-best in terms of standard deviation (SD = 1.68 kcal/mol). The function was then further improved by inclusion of different terms taking into account repulsion and loss of entropy upon binding, and we show that it is capable of recovering native binding poses up to 80\% of the time. All of the programs, tools, and protein sets are released in the Supporting Information or as open-source programs.}, url = {http://pubs.acs.org/doi/abs/10.1021/acs.jcim.6b00610}, author = {Debroise, Th{\'e}au and Shakhnovich, EugeneI. and Ch{\'e}ron, Nicolas} } @article {srinivasan2017rational, title = {Rational design of novel allosteric dihydrofolate reductase inhibitors showing antibacterial-effects on drug-resistant E. coli escape-variants}, journal = {ACS Chemical Biology}, year = {2017}, publisher = {ACS Publications}, abstract = {In drug discovery, systematic variations of substituents on a common scaffold and bioisosteric replacements are often used to generate diversity and obtain molecules with better biological effects. However, this could saturate the small-molecule diversity pool resulting in drug resistance. On the other hand, conventional drug discovery relies on targeting known pockets on protein surfaces leading to drug resistance by mutations of critical pocket residues. Here, we present a two-pronged strategy of designing novel drugs that target unique pockets on a protein{\textquoteright}s surface to overcome the above problems. Dihydrofolate reductase, DHFR, is a critical enzyme involved in thymidine and purine nucleotide biosynthesis. Several classes of compounds that are structural analogues of the substrate dihydrofolate have been explored for their antifolate activity. Here, we describe 10 novel small-molecule inhibitors of\ Escherichia coli\ DHFR,\ EcDHFR, belonging to the stilbenoid, deoxybenzoin, and chalcone family of compounds discovered by a combination of pocket-based virtual ligand screening and systematic scaffold hopping. These inhibitors show a unique uncompetitive or noncompetitive inhibition mechanism, distinct from those reported for all known inhibitors of DHFR, indicative of binding to a unique pocket distinct from either substrate or cofactor-binding pockets. Furthermore, we demonstrate that rescue mutants of\ EcDHFR, with reduced affinity to all known classes of DHFR inhibitors, are inhibited at the same concentration as the wild-type. These compounds also exhibit antibacterial activity against\ E. coli\ harboring the drug-resistant variant of DHFR. This discovery is the first report on a novel class of inhibitors targeting a unique pocket on\ EcDHFR.}, url = {http://pubs.acs.org/doi/abs/10.1021/acschembio.7b00175}, author = {Srinivasan, Bharath and Rodrigues, Jo{\~a}o V and Tonddast-Navaei, Sam and Eugene Shakhnovich and Skolnick, Jeffrey} } @article {gilson2017role, title = {The role of evolutionary selection in the dynamics of protein structure evolution}, journal = {Biophysical Journal}, volume = {112}, year = {2017}, pages = {1350{\textendash}1365}, publisher = {Elsevier}, abstract = {Homology modeling is a powerful tool for predicting a protein{\textquoteright}s structure. This approach is successful because proteins whose sequences are only 30\% identical still adopt the same structure, while structure similarity rapidly deteriorates beyond the 30\% threshold. By studying the divergence of protein structure as sequence evolves in real proteins and in evolutionary simulations, we show that this nonlinear sequence-structure relationship emerges as a result of selection for protein folding stability in divergent evolution. Fitness constraints prevent the emergence of unstable protein evolutionary intermediates, thereby enforcing evolutionary paths that preserve protein structure despite broad sequence divergence. However, on longer timescales, evolution is punctuated by rare events where the fitness barriers obstructing structure evolution are overcome and discovery of new structures occurs. We outline biophysical and evolutionary rationale for broad variation in protein family sizes, prevalence of compact structures among ancient proteins, and more rapid structure evolution of proteins with lower packing density.}, url = {http://www.sciencedirect.com/science/article/pii/S0006349517302436}, author = {Gilson, Amy I and Marshall-Christensen, Ahmee and Jeong-Mo Choi and Shakhnovich, EugeneI.} } @article {1006961, title = {Optimization of lag phase shapes the evolution of a bacterial enzyme}, journal = {Nature Ecology \& Evolution}, volume = {1}, year = {2017}, month = {2017/04/28/onlin}, pages = {0149 - }, publisher = {Macmillan Publishers Limited, part of Springer Nature.}, abstract = {

Mutations provide the variation that drives evolution, yet their effects on fitness remain poorly understood. Here we explore how mutations in the essential enzyme adenylate kinase (Adk) of\ Escherichia coli\ affect multiple phases of population growth. We introduce a biophysical fitness landscape for these phases, showing how they depend on molecular and cellular properties of Adk. We find that Adk catalytic capacity in the cell (the product of activity and abundance) is the major determinant of mutational fitness effects. We show that bacterial lag times are at a well-defined optimum with respect to Adk{\textquoteright}s catalytic capacity, while exponential growth rates are only weakly affected by variation in Adk. Direct pairwise competitions between strains show how environmental conditions modulate the outcome of a competition where growth rates and lag times have a tradeoff, shedding light on the multidimensional nature of fitness and its importance in the evolutionary optimization of enzymes.

}, url = {http://dx.doi.org/10.1038/s41559-017-0149}, author = {Adkar, Bharat V. and Manhart, Michael and Bhattacharyya, Sanchari and Tian, Jian and Musharbash, Michael and Shakhnovich, EugeneI.} } @article {bershtein2017bridging, title = {Bridging the physical scales in evolutionary biology: from protein sequence space to fitness of organisms and populations}, journal = {Current Opinion in Structural Biology}, volume = {42}, year = {2017}, month = {31 October 2016}, pages = {31{\textendash}40}, publisher = {Elsevier}, abstract = {

Bridging the gap between the molecular properties of proteins and organismal/population fitness is essential for understanding evolutionary processes. This task requires the integration of the several physical scales of biological organization, each defined by a distinct set of mechanisms and constraints, into a single unifying model. The\ molecular\ scale is dominated by the constraints imposed by the physico-chemical properties of proteins and their substrates, which give rise to trade-offs and epistatic (non-additive) effects of mutations. At the\ systems\ scale, biological networks modulate protein expression and can either buffer or enhance the fitness effects of mutations. The\ population\ scale is influenced by the mutational input, selection regimes, and stochastic changes affecting the size and structure of populations, which eventually determine the evolutionary fate of mutations. Here, we summarize the recent advances in theory, computer simulations, and experiments that advance our understanding of the links between various physical scales in biology.

}, url = {http://www.sciencedirect.com/science/article/pii/S0959440X16301762}, author = {Shimon Bershtein and Adrian W.R. Serohijos and Shakhnovich, EugeneI.} } @article {zhou2016virtual, title = {Virtual Screening of Human O-GlcNAc Transferase Inhibitors}, journal = {Chinese Journal of Chemical Physics}, volume = {29}, year = {2016}, pages = {374{\textendash}380}, publisher = {CHINESE PHYSICAL SOC PO BOX 603, BEIJING 100080, PEOPLES R CHINA}, abstract = {O-GlcNAc transferase (OGT) is one of essential mammalian enzymes, which catalyze the transfer of N-acetylglucosamine from UDP-N-acetylglucosamine (UDP-GlcNAc) to hydroxyl groups of serines and threonines (Ser/Thr) in proteins. Dysregulations of cellular O-GlcNAc have been implicated in diabetes, neurodegenerative disease, and cancer, which brings great interest in developing potent and specific small-molecular OGT inhibitors. In this work, we performed virtual screening on OGT catalytic site to identify potential inhibitors. 7134792 drug-like compounds from ZINC (a free database of commercially available compounds for virtual screening) and 4287550 compounds generated by FOG (fragment optimized growth program) were screened and the top 116 compounds ranked by docking score were analyzed. By comparing the screening results, we found FOG program can generate more compounds with better docking scores than ZINC. The top ZINC compounds ranked by docking score were grouped into two classes, which held the binding positions of UDP and GlcNAc of UDP-GlcNAc. Combined with individual fragments in binding pocket, de novo compounds were designed and proved to have better docking score. The screened and designed compounds may become a starting point for developing new drugs.}, url = {https://cps.scitation.org/doi/abs/10.1063/1674-0068/29/cjcp1510211?journalCode=cjp}, author = {Zhou, Qing-tong and Liang, Hao-jun and Shakhnovich, EugeneI.} } @article {jacquin2016benchmarking, title = {Benchmarking inverse statistical approaches for protein structure and design with exactly solvable models}, journal = {PLoS Comput Biol}, volume = {12}, year = {2016}, month = {13 May, 2016}, pages = {e1004889}, publisher = {Public Library of Science}, abstract = {

Inverse statistical approaches to determine protein structure and function from Multiple Sequence Alignments (MSA) are emerging as powerful tools in computational biology. However the underlying assumptions of the relationship between the inferred effective Potts Hamiltonian and real protein structure and energetics remain untested so far. Here we use lattice protein model (LP) to benchmark those inverse statistical approaches. We build MSA of highly stable sequences in target LP structures, and infer the effective pairwise Potts Hamiltonians from those MSA. We find that inferred Potts Hamiltonians reproduce many important aspects of {\textquoteleft}true{\textquoteright} LP structures and energetics. Careful analysis reveals that effective pairwise couplings in inferred Potts Hamiltonians depend not only on the energetics of the native structure but also on competing folds; in particular, the coupling values reflect both positive design (stabilization of native conformation) and negative design (destabilization of competing folds). In addition to providing detailed structural information, the inferred Potts models used as protein Hamiltonian for design of new sequences are able to generate with high probability completely new sequences with the desired folds, which is not possible using independent-site models. Those are remarkable results as the effective LP Hamiltonians used to generate MSA are not simple pairwise models due to the competition between the folds. Our findings elucidate the reasons for the success of inverse approaches to the modelling of proteins from sequence data, and their limitations.

}, url = {http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004889}, author = {Jacquin, Hugo and Amy I. Gilson and Shakhnovich, EugeneI. and Cocco, Simona and Monasson, R{\'e}mi} } @article {rodrigues2016biophysical, title = {Biophysical principles predict fitness landscapes of drug resistance}, journal = {Proceedings of the National Academy of Sciences}, volume = {113}, year = {2016}, month = {27 January 2016}, pages = {E1470{\textendash}E1478}, publisher = {National Acad Sciences}, abstract = {

Fitness landscapes of drug resistance constitute powerful tools to elucidate mutational pathways of antibiotic escape. Here, we developed a predictive biophysics-based fitness landscape of trimethoprim (TMP) resistance for\ Escherichia coli\ dihydrofolate reductase (DHFR). We investigated the activity, binding, folding stability, and intracellular abundance for a complete set of combinatorial DHFR mutants made out of three key resistance mutations and extended this analysis to DHFR originated from\ Chlamydia muridarum\ and\ Listeria grayi. We found that the acquisition of TMP resistance via decreased drug affinity is limited by a trade-off in catalytic efficiency. Protein stability is concurrently affected by the resistant mutants, which precludes a precise description of fitness from a single molecular trait. Application of the kinetic flux theory provided an accurate model to predict resistance phenotypes (IC50) quantitatively from a unique combination of the in vitro protein molecular properties. Further, we found that a controlled modulation of the GroEL/ES chaperonins and Lon protease levels affects the intracellular steady-state concentration of DHFR in a mutation-specific manner, whereas IC50\ is changed proportionally, as indeed predicted by the model. This unveils a molecular rationale for the pleiotropic role of the protein quality control machinery on the evolution of antibiotic resistance, which, as we illustrate here, may drastically confound the evolutionary outcome. These results provide a comprehensive quantitative genotype{\textendash}phenotype map for the essential enzyme that serves as an important target of antibiotic and anticancer therapies.

}, url = {http://www.pnas.org/content/113/11/E1470.abstract}, author = {Rodrigues, Jo{\~a}o V and Shimon Bershtein and Li, Annabel and Lozovsky, Elena R and Hartl, Daniel L and Shakhnovich, EugeneI.} } @article {cheron2016evolutionary, title = {Evolutionary dynamics of viral escape under antibodies stress: A biophysical model}, journal = {Protein Science}, volume = {25}, year = {2016}, pages = {1332-1340}, publisher = {Wiley Online Library}, abstract = {

Viruses constantly face the selection pressure of antibodies, either from innate immune response of the host or from administered antibodies for treatment. We explore the interplay between the biophysical properties of viral proteins and the population and demographic variables in the viral escape. The demographic and population genetics aspect of the viral escape have been explored before; however one important assumption was the a priori distribution of fitness effects (DFE). Here, we relax this assumption by instead considering a realistic biophysics-based genotype-phenotype relationship for RNA viruses escaping antibodies stress. In this model the DFE is itself an evolvable property that depends on the genetic background (epistasis) and the distribution of biophysical effects of mutations, which is informed by biochemical experiments and theoretical calculations in protein engineering. We quantitatively explore in silico the viability of viral populations under antibodies pressure and derive the phase diagram that defines the fate of the virus population (extinction or escape from stress) in a range of viral mutation rates and antibodies concentrations. We find that viruses are most resistant to stress at an optimal mutation rate (OMR) determined by the competition between supply of beneficial mutation to facilitate escape from stressors and lethal mutagenesis caused by excess of destabilizing mutations. We then show the quantitative dependence of the OMR on genome length and viral burst size. We also recapitulate the experimental observation that viruses with longer genomes have smaller mutation rate per nucleotide.

}, url = {http://onlinelibrary.wiley.com/doi/10.1002/pro.2915/abstract}, author = {Ch{\'e}ron, Nicolas and Adrian W.R. Serohijos and Jeong-Mo Choi and Shakhnovich, EugeneI.} } @article {serebryany2016aggregation, title = {An internal disulfide locks a misfolded aggregation-prone intermediate in cataract-linked mutants of human γD-crystallin}, journal = {Journal of Biological Chemistry}, volume = {291}, year = {2016}, pages = {19172{\textendash}19183}, publisher = {ASBMB}, abstract = {

Considerable mechanistic insight has been gained into amyloid aggregation; however, a large number of non-amyloid protein aggregates are considered {\textquotedblleft}amorphous,{\textquotedblright} and in most cases, little is known about their mechanisms. Amorphous aggregation of γ-crystallins in the eye lens causes cataract, a widespread disease of aging. We combined simulations and experiments to study the mechanism of aggregation of two γD-crystallin mutants, W42R and W42Q: the former a congenital cataract mutation, and the latter a mimic of age-related oxidative damage. We found that formation of an internal disulfide was necessary and sufficient for aggregation under physiological conditions. Two-chain all-atom simulations predicted that one non-native disulfide in particular, between Cys32\ and Cys41, was likely to stabilize an unfolding intermediate prone to intermolecular interactions. Mass spectrometry and mutagenesis experiments confirmed the presence of this bond in the aggregates and its necessity for oxidative aggregation under physiological conditions\ in vitro. Mining the simulation data linked formation of this disulfide to extrusion of the N-terminal β-hairpin and rearrangement of the native β-sheet topology. Specific binding between the extruded hairpin and a distal β-sheet, in an intermolecular chain reaction similar to domain swapping, is the most probable mechanism of aggregate propagation.

}, url = {http://www.jbc.org/content/291/36/19172.short}, author = {Serebryany, Eugene and Woodard, Jaie C and Adkar, Bharat V and Shabab, Mohammed and King, Jonathan A and Shakhnovich, EugeneI.} } @article {cheron2016opengrowth, title = {OpenGrowth: An Automated and Rational Algorithm for Finding New Protein Ligands}, journal = {Journal of Medicinal Chemistry}, volume = {59}, year = {2016}, note = {

PMID: 26356253

}, pages = {4171-4188}, abstract = {

We present a new open-source software, called OpenGrowth, which aims to create de novo ligands by connecting small organic fragments in the active site of proteins. Molecule growth is biased to produce structures that statistically resemble drugs in an input training database. Consequently, the produced molecules have superior synthetic accessibility and pharmacokinetic properties compared with randomly grown molecules. The growth process can take into account the flexibility of the target protein and can be started from a seed to mimic R-group strategy or fragment-based drug discovery. Primary applications of the software on the HIV-1 protease allowed us to quickly identify new inhibitors with a predicted\ Kd\ as low as 18 nM. We also present a graphical user interface that allows a user to select easily the fragments to include in the growth process. OpenGrowth is released under the GNU GPL license and is available free of charge on the authors{\textquoteright} website and at\ http://opengrowth.sourceforge.net/.

}, doi = {10.1021/acs.jmedchem.5b00886}, url = {http://dx.doi.org/10.1021/acs.jmedchem.5b00886}, author = {Ch{\'e}ron, Nicolas and Jasty, Naveen and Shakhnovich, EugeneI.} } @article {woodard2016domainswap, title = {A Simple Model of Protein Domain Swapping in Crowded Cellular Environments}, journal = {Biophysical Journal}, volume = {110}, year = {2016}, month = { 7 June 2016}, pages = {2367{\textendash}2376}, publisher = {Elsevier}, abstract = {

Domain swapping in proteins is an important mechanism of functional and structural innovation. However, despite its ubiquity and importance, the physical mechanisms that lead to domain swapping are poorly understood. Here, we present a simple two-dimensional coarse-grained model of protein domain swapping in the cytoplasm. In our model, two-domain proteins partially unfold and diffuse in continuous space. Monte Carlo multiprotein simulations of the model reveal that domain swapping occurs at intermediate temperatures, whereas folded dimers and folded monomers prevail at low temperatures, and partially unfolded monomers predominate at high temperatures. We use a simplified amino acid alphabet consisting of four residue types, and find that the oligomeric state at a given temperature depends on the sequence of the protein. We also show that hinge strain between domains can promote domain swapping, consistent with experimental observations for real proteins. Domain swapping depends nonmonotonically on the protein concentration, with domain-swapped dimers occurring at intermediate concentrations and nonspecific interactions between partially unfolded proteins occurring at high concentrations. For folded proteins, we recover the result obtained in three-dimensional lattice simulations, i.e., that functional dimerization is most prevalent at\ intermediate temperatures and nonspecific interactions increase at low temperatures.

}, url = {http://www.sciencedirect.com/science/article/pii/S0006349516302363}, author = {Woodard, Jaie C and Dunatunga, Sachith and Shakhnovich, EugeneI.} } @article {jacobs2016paths, title = {Structure-based prediction of protein-folding transition paths}, journal = {Biophysical Journal}, volume = {111}, year = {2016}, month = {6 Sep, 2016}, pages = {925{\textendash}936}, publisher = {Elsevier}, abstract = {

We propose a general theory to describe the distribution of protein-folding transition paths. We show that transition paths follow a predictable sequence of high-free-energy transient states that are separated by free-energy barriers. Each transient state corresponds to the assembly of one or more discrete, cooperative units, which are determined directly from the native structure. We show that the transition state on a folding pathway is reached when a small number of critical contacts are formed between a specific set of substructures, after which folding proceeds downhill in free energy. This approach suggests a natural resolution for distinguishing parallel folding pathways and provides a simple means to predict the rate-limiting step in a folding reaction. Our theory identifies a common folding mechanism for proteins with diverse native structures and establishes general principles for the self-assembly of polymers with specific interactions.

}, url = {http://www.sciencedirect.com/science/article/pii/S0006349516304775}, author = {Jacobs, William M. and Shakhnovich, EugeneI.} } @article {938156, title = {Transient protein-protein interactions perturb E. coli metabolome and cause gene dosage toxicity}, journal = {Elife}, volume = {5}, year = {2016}, month = {10 Dec, 2016}, pages = {e20309}, abstract = {

Gene dosage toxicity (GDT) is an important factor that determines optimal levels of protein abundances, yet its molecular underpinnings remain unknown. Here, we demonstrate that overexpression of DHFR in\ E. coli\ causes a toxic metabolic imbalance triggered by interactions with several functionally related enzymes. Though deleterious in the overexpression regime, surprisingly, these interactions are beneficial at physiological concentrations, implying their functional significance\ in vivo. Moreover, we found that overexpression of orthologous DHFR proteins had minimal effect on all levels of cellular organization - molecular, systems, and phenotypic, in sharp contrast to\ E. coliDHFR. Dramatic difference of GDT between {\textquoteright}E. coli{\textquoteright}s self{\textquoteright} and {\textquoteright}foreign{\textquoteright} proteins suggests the crucial role of evolutionary selection in shaping protein-protein interaction (PPI) networks at the whole proteome level. This study shows how protein overexpression perturbs a dynamic metabolon of weak yet potentially functional PPI, with consequences for the metabolic state of cells and their fitness.

}, url = {https://elifesciences.org/content/5/e20309/ }, author = {Bhattacharyya, Sanchari and Shimon Bershtein and Jin Yan and Tijda Argun and Amy I. Gilson and Sunia A. Trauger and Shakhnovich, EugeneI.} } @article {784406, title = {Isolation and Analysis of Rare Norovirus Recombinants from Coinfected Mice Using Drop-Based Microfluidics}, journal = {J. Virol.}, volume = {89}, year = {2015}, pages = {7722-34}, abstract = {Human noroviruses (HuNoVs) are positive-sense RNA viruses that can cause severe, highly infectious gastroenteritis. HuNoV outbreaks are frequently associated with recombination between circulating strains. Strain genotyping and phylogenetic analyses show that noroviruses often recombine in a highly conserved region near the junction of the viral polyprotein (ORF1) and capsid (ORF2) genes and occasionally within the RNA-dependent RNA polymerase (RdRP) gene. Although genotyping methods are useful for tracking changes in circulating viral populations, they only report the dominant recombinant strains and do not elucidate the frequency or range of recombination events. Furthermore, the relatively low frequency of recombination in RNA viruses has limited studies to cell culture or\ in vitro\ systems that do not reflect the complexities and selective pressures present in an infected organism. Using two murine norovirus (MNV) strains to model co-infection, we developed a microfluidic platform to amplify, detect, and recover individual recombinants following\ in vitro\ and\ in vivo\ co-infection. One-step RT-PCR was performed in picoliter drops with primers that identified the wild-type and recombinant progenies and scanned for recombination breakpoints at approximately 1-kb intervals. We detected recombination between MNV strains at multiple loci spanning the viral protease, RdRP, and capsid ORFs and isolated individual recombinant RNA genomes that were present at a frequency of 1/300,000 or greater. This study is the first to examine norovirus recombination following co-infection of an animal and suggests that the exchange of RNA among viral genomes in the infected host occurs in multiple locations and is an important driver of genetic diversity.}, url = {http://jvi.asm.org/content/early/2015/05/11/JVI.01137-15.short}, author = {Zhang, H. and Cockrell, S. K. and Kolawole, A. O. and Rotem, A. and Serohijos, A. W. and Chang, C. B. and Tao, Y. and Mehoke, T. S. and Han, Y. and Lin, J. S. and Giacobbi, N. S. and Feldman, A. B. and Shakhnovich, E. and D. A. Weitz and Wobus, C. E. and Pipas, J. M.} } @article {784421, title = {Mechanical Model of Globular Transition in Polymers}, journal = {Chem Plus Chem}, volume = {80}, year = {2015}, pages = {37-41}, abstract = {In complex, multicomponent systems, polymers often undergo phase transitions between distinct conformations. This paper reports a millimeter-scale granular model of coil-to-globule transitions: one {\textquotedblleft}polymer{\textquotedblright} chain{\textemdash}a cylinders-on-a-string {\textquotedblleft}pearl necklace{\textquotedblright}{\textemdash}and many spheres, all shaken on a horizontal surface. It is possible to describe the behavior of this granular system by using formalisms generally used in statistical physics of polymers. Two sets of experiments allowed the observation of first- and second-order coil-to-globule transitions. The model shows that the competition between long- and short-range interactions leads to a first-order transition. Well-designed granular system represents another kind of approach to the study of polymer phase transitions and might inspire future designs of polymer-like mesoscale systems.}, url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/cplu.201402203}, author = {Ticard, S. and Shepherd, R. F. and Stan, C. A. and Snyder, P. W. and Cademartiri, R. and Zhu, D. and Aranson, I. S. and Shakhnovich, E. I. and G. M. Whitesides} } @article {784416, title = {Searching the Sequence Space for Potent Aptamers Using SELEX in Silico}, journal = {J. Chem. Theory Comput.}, volume = {11}, year = {2015}, pages = {5939-46}, abstract = {To isolate functional nucleic acids that bind to defined targets with high affinity and specificity, which are known as aptamers, the systematic evolution of ligands by exponential enrichment (SELEX) methodology has emerged as the preferred approach. Here, we propose a computational approach, SELEX in silico, that allows the sequence space to be more thoroughly explored regarding binding of a certain target. Our approach consists of two steps: (i) secondary structure-based sequence screening, which aims to collect the sequences that can form a desired RNA motif as an enhanced initial library, followed by (ii) sequence enrichment regarding target binding by molecular dynamics simulation-based virtual screening. Our SELEX in silico method provided a practical computational solution to three key problems in aptamer sequence searching: design of nucleic acid libraries, knowledge of sequence enrichment, and identification of potent aptamers. Six potent theophylline-binding aptamers, which were isolated by SELEX in silico from a sequence space containing 413\ sequences, were experimentally verified to bind theophylline with high affinity:\ Kd\ ranging from 0.16 to 0.52 μM, compared with the dissociation constant of the original aptamer-theophylline, 0.32 μM. These results demonstrate the significant potential of SELEX in silico as a new method for aptamer discovery and optimization.}, url = {http://dx.doi.org/10.1021/acs.jctc.5b00707}, author = {Zhou, Q. and Xia, X. and Luo, Z. and Liang, H. and Shakhnovich, E.} } @article {784411, title = {Systems-Level Response to Point Mutations in a Core Metabolic Enzyme Modulates Genotype-Phenotype Relationship}, journal = {Cell Reports}, volume = {11}, year = {2015}, pages = {645-6}, abstract = {Linking the molecular effects of mutations to fitness is\ central to understanding evolutionary dynamics. Here, we establish a quantitative relation between the\ global\ effect of mutations on the\ E.\ coliproteome\ and bacterial fitness. We created\ E.\ coli\ strains with specific destabilizing mutations in the\ chromosomalfolA\ gene encoding\ dihydrofolate reductase(DHFR) and quantified the ensuing changes in the abundances of 2,000+\ E.\ coli\ proteins in mutant strains using\ tandem mass tags\ with subsequent LC-MS/MS. mRNA abundances in the same\ E.\ coli\ strains were also quantified. The\ proteomic\ effects of mutations in\ DHFRare quantitatively linked to\ phenotype: the SDs of the distributions of logarithms of relative (to WT) protein abundances anticorrelate with bacterial growth rates. Proteomes hierarchically cluster first by media conditions, and within each condition, by the severity of the perturbation to DHFR function. These results highlight the importance of a systems-level layer in the genotype-phenotype relationship.}, url = {http://dx.doi.org/10.1016/j.celrep.2015.03.051 }, author = {Bershtein, S. and Choi, J.-M. and Bhattacharyya, S. and Budnik, B. and Shakhnovich, E.} } @article {732451, title = {Protein Homeostasis Imposes a Barrier on Functional Integration of Horizontally Transferred Genes in Bacteria.}, journal = {PLoS Genet}, volume = {11}, year = {2015}, month = {2015 Oct}, pages = {e1005612}, abstract = {Horizontal gene transfer (HGT) plays a central role in bacterial evolution, yet the molecular and cellular constraints on functional integration of the foreign genes are poorly understood. Here we performed inter-species replacement of the chromosomal folA gene, encoding an essential metabolic enzyme dihydrofolate reductase (DHFR), with orthologs from 35 other mesophilic bacteria. The orthologous inter-species replacements caused a marked drop (in the range 10-90\%) in bacterial growth rate despite the fact that most orthologous DHFRs are as stable as E.coli DHFR at 37{\textdegree}C and are more catalytically active than E. coli DHFR. Although phylogenetic distance between E. coli and orthologous DHFRs as well as their individual molecular properties correlate poorly with growth rates, the product of the intracellular DHFR abundance and catalytic activity (kcat/KM), correlates strongly with growth rates, indicating that the drop in DHFR abundance constitutes the major fitness barrier to HGT. Serial propagation of the orthologous strains for ~600 generations dramatically improved growth rates by largely alleviating the fitness barriers. Whole genome sequencing and global proteome quantification revealed that the evolved strains with the largest fitness improvements have accumulated mutations that inactivated the ATP-dependent Lon protease, causing an increase in the intracellular DHFR abundance. In one case DHFR abundance increased further due to mutations accumulated in folA promoter, but only after the lon inactivating mutations were fixed in the population. Thus, by apparently distinguishing between self and non-self proteins, protein homeostasis imposes an immediate and global barrier to the functional integration of foreign genes by decreasing the intracellular abundance of their products. Once this barrier is alleviated, more fine-tuned evolution occurs to adjust the function/expression of the transferred proteins to the constraints imposed by the intracellular environment of the host organism.}, issn = {1553-7404}, doi = {10.1371/journal.pgen.1005612}, url = {http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1005612}, author = {Shimon Bershtein and Adrian W.R. Serohijos and Bhattacharyya, Sanchari and Manhart, Michael and Jeong-Mo Choi and Wanmeng Mu and Jingwen Zhou and Shakhnovich, EugeneI.} } @article {732456, title = {Thermal stabilization of dihydrofolate reductase using monte carlo unfolding simulations and its functional consequences.}, journal = {PLoS Comput. Biol.}, volume = {11}, year = {2015}, month = {2015 Apr}, pages = {e1004207}, abstract = {

Design of proteins with desired thermal properties is important for scientific and biotechnological applications. Here we developed a theoretical approach to predict the effect of mutations on protein stability from non-equilibrium unfolding simulations. We establish a relative measure based on apparent simulated melting temperatures that is independent of simulation length and, under certain assumptions, proportional to equilibrium stability, and we justify this theoretical development with extensive simulations and experimental data. Using our new method based on all-atom Monte-Carlo unfolding simulations, we carried out a saturating mutagenesis of Dihydrofolate Reductase (DHFR), a key target of antibiotics and chemotherapeutic drugs. The method predicted more than 500 stabilizing mutations, several of which were selected for detailed computational and experimental analysis. We find a highly significant correlation of r=0.65-0.68 between predicted and experimentally determined melting temperatures and unfolding denaturant concentrations for WT DHFR and 42 mutants. The correlation between energy of the native state and experimental denaturation temperature was much weaker, indicating the important role of entropy in protein stability. The most stabilizing point mutation was D27F, which is located in the active site of the protein, rendering it inactive. However for the rest of mutations outside of the active site we observed a weak yet statistically significant positive correlation between thermal stability and catalytic activity indicating the lack of a stability-activity tradeoff for DHFR. By combining stabilizing mutations predicted by our method, we created a highly stable catalytically active E. coli DHFR mutant with measured denaturation temperature 7.2{\textdegree}C higher than WT. Prediction results for DHFR and several other proteins indicate that computational approaches based on unfolding simulations are useful as a general technique to discover stabilizing mutations.

}, issn = {1553-7358}, doi = {10.1371/journal.pcbi.1004207}, url = {https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004207}, author = {Tian, Jian and Woodard, Jaie C and Whitney, Anna and Shakhnovich, EugeneI.} } @article {639946, title = {Repurposing of rutin for the inhibition of norovirus replication}, journal = {Archives of Virology}, volume = {160}, year = {2015}, month = {2015}, pages = {2353-2358}, publisher = {Springer Vienna}, abstract = {

Drug repurposing is a strategy employed to circumvent some of the bottlenecks involved in drug development, such as the cost and time needed for developing new molecular entities. Noroviruses cause recurrent epidemics and sporadic outbreaks of gastroenteritis associated with significant mortality and economic costs, but no treatment has been approved to date. Herein, a library of molecules previously used in humans was screened to find compounds with anti-noroviral activity. Antiviral testing for four selected compounds against murine norovirus infection revealed that rutin has anti-murine norovirus activity in cell-based assays.

}, isbn = {0304-8608}, url = {http://dx.doi.org/10.1007/s00705-015-2495-y}, author = {Ch{\'e}ron, Nicolas and Yu, Chenchen and Kolawole, AbimbolaO. and Shakhnovich, EugeneI. and Wobus, ChristianeE.} } @article {601386, title = {Is Catalytic Activity of Chaperones a Selectable Trait for the Emergence of Heat Shock Response?}, journal = {Biophysical Journal}, volume = {108}, year = {2015}, note = {

doi: 10.1016/j.bpj.2014.11.3468

}, month = {2015}, pages = {438 - 448}, publisher = {Elsevier}, abstract = {

Although heat shock response is ubiquitous in bacterial cells, the underlying physical chemistry behind heat shock response remains poorly understood. To study the response of cell populations to heat shock we employ a physics-based ab initio model of living cells where protein biophysics (i.e., folding and protein-protein interactions in crowded cellular environments) and important aspects of proteins homeostasis are coupled with realistic population dynamics simulations. By postulating a genotype-phenotype relationship we define a cell division rate in terms of functional concentrations of proteins and protein complexes, whose Boltzmann stabilities of folding and strengths of their functional interactions are exactly evaluated from their sequence information. We compare and contrast evolutionary dynamics for two models of chaperon action. In the active model, foldase chaperones function as nonequilibrium machines to accelerate the rate of protein folding. In the passive model, holdase chaperones form reversible complexes with proteins in their misfolded conformations to maintain their solubility. We find that only cells expressing foldase chaperones are capable of genuine heat shock response to the increase in the amount of unfolded proteins at elevated temperatures. In response to heat shock, cells? limited resources are redistributed differently for active and passive models. For the active model, foldase chaperones are overexpressed at the expense of downregulation of high abundance proteins, whereas for the passive model; cells react to heat shock by downregulating their high abundance proteins, as their low abundance proteins are upregulated.

}, url = {http://dx.doi.org/10.1016/j.bpj.2014.11.3468}, author = {{\c C}etinba{\c s}, Murat and Shakhnovich, Eugene~I.} } @article {523021, title = {Minimalistic Predictor of Protein Binding Energy: Contribution of Solvation Factor to Protein Binding}, journal = {Biophysical Journal}, volume = {108}, year = {2015}, month = {2015/2/17/}, pages = {795 - 798}, abstract = {

It has long been known that solvation plays an important role in protein-protein interactions. Here, we use a minimalistic solvation-based model for predicting protein binding energy to estimate quantitatively the contribution of the solvation factor in protein binding. The factor is described by a simple linear combination of buried surface areas according to amino-acid types. Even without structural optimization, our minimalistic model demonstrates a predictive power comparable to more complex methods, making the proposed approach the basis for high throughput applications. Application of the model to a proteomic database shows that receptor-substrate complexes involved in signaling have lower affinities than enzyme-inhibitor and antibody-antigen complexes, and they differ by chemical compositions on interfaces. Also, we found that protein complexes with components that come from the same genes generally have lower affinities than complexes formed by proteins from different genes, but in this case the difference originates from different interface areas. The model was implemented in the software PYTHON, and the source code can be found on the Shakhnovich group webpage: http://faculty.chemistry.harvard.edu/shakhnovich/software.

}, isbn = {0006-3495}, url = {http://www.sciencedirect.com/science/article/pii/S0006349515000582}, author = {Jeong-Mo Choi and Serohijos, Adrian~W.R. and Murphy, Sean and Lucarelli, Dennis and Lofranco, Leo~L. and Feldman, Andrew and Shakhnovich, Eugene~I.} } @article {732461, title = {The influence of selection for protein stability on dN/dS estimations.}, journal = {Genome Biol Evol}, volume = {6}, year = {2014}, month = {2014 Oct}, pages = {2956-67}, abstract = {Understanding the relative contributions of various evolutionary processes-purifying selection, neutral drift, and adaptation-is fundamental to evolutionary biology. A common metric to distinguish these processes is the ratio of nonsynonymous to synonymous substitutions (i.e., dN/dS) interpreted from the neutral theory as a null model. However, from biophysical considerations, mutations have non-negligible effects on the biophysical properties of proteins such as folding stability. In this work, we investigated how stability affects the rate of protein evolution in phylogenetic trees by using simulations that combine explicit protein sequences with associated stability changes. We first simulated myoglobin evolution in phylogenetic trees with a biophysically realistic approach that accounts for 3D structural information and estimates of changes in stability upon mutation. We then compared evolutionary rates inferred directly from simulation to those estimated using maximum-likelihood (ML) methods. We found that the dN/dS estimated by ML methods (ωML) is highly predictive of the per gene dN/dS inferred from the simulated phylogenetic trees. This agreement is strong in the regime of high stability where protein evolution is neutral. At low folding stabilities and under mutation-selection balance, we observe deviations from neutrality (per gene dN/dS > 1 and dN/dS < 1). We showed that although per gene dN/dS is robust to these deviations, ML tests for positive selection detect statistically significant per site dN/dS > 1. Altogether, we show how protein biophysics affects the dN/dS estimations and its subsequent interpretation. These results are important for improving the current approaches for detecting positive selection.}, keywords = {Evolution, Molecular, Phylogeny, Protein Stability, Proteins, Selection, Genetic}, issn = {1759-6653}, doi = {10.1093/gbe/evu223}, author = {Pouria Dasmeh and Adrian W.R. Serohijos and Kasper P. Kepp and Eugene I. Shakhnovich} } @article {421591, title = {Evolution of Specificity in Protein-Protein Interactions}, journal = {Biophysical Journal}, volume = {107}, year = {2014}, month = {2014}, pages = {1686 - 1696}, abstract = {Hub proteins are proteins that maintain promiscuous molecular recognition. Because they are reported to play essential roles in cellular control, there has been a special interest in the study of their structural and functional properties, yet the mechanisms by which they evolve to maintain functional interactions are poorly understood. By combining biophysical simulations of coarse-grained proteins and analysis of proteins-complex crystallographic structures, we seek to elucidate those mechanisms. We focus on two types of hub proteins: Multi hubs, which interact with their partners through different interfaces, and Singlish hubs, which do so through a single interface. We show that loss of structural stability is required for the evolution of protein-protein-interaction (PPI) networks, and it is more profound in Singlish hub systems. In addition, different ratios of hydrophobic to electrostatic interfacial amino acids are shown to support distinct network topologies (i.e., Singlish and Multi systems), and therefore underlie a fundamental design principle of PPI in a crowded environment. We argue that the physical nature of hydrophobic and electrostatic interactions, in particular, their favoring of either same-type interactions (hydrophobic-hydrophobic), or opposite-type interactions (negatively-positively charged) plays a key role in maintaining the network topology while allowing the protein amino acid sequence to evolve.}, url = {http://dx.doi.org/10.1016/j.bpj.2014.08.004}, author = {Peleg, O. and Choi, J.-M. and Eugene I. Shakhnovich} } @article {407866, title = {Contribution of selection for protein folding stability in shaping the patterns of polymorphisms in coding region}, journal = {Mol. Biol \& Evol.}, volume = {31}, year = {2014}, month = {2014}, pages = {165 - 176}, abstract = {The patterns of polymorphisms in genomes are imprints of the evolutionary forces at play in nature. In particular, polymorphisms have been extensively used to infer the fitness effects of mutations and their dynamics of fixation. However, the role and contribution of molecular biophysics to these observations remain unclear. Here, we couple robust findings from protein biophysics, enzymatic flux theory, the selection against the cytotoxic effects of protein misfolding, and explicit population dynamics simulations in the polyclonal regime. First, we recapitulate results on the dynamics of clonal interference and on the shape of the DFE, thus providing them with a molecular and mechanistic foundation. Second, we predict that if evolution is indeed under the dynamic equilibrium of mutation-selection balance, the fraction of stabilizing and destabilizing mutations is almost equal among single-nucleotide polymorphisms segregating at high allele frequencies. This prediction is proven true for polymorphisms in the human coding region. Overall, our results show how selection for protein folding stability predominantly shapes the patterns of polymorphisms in coding regions.}, url = {http://dx.doi.org/10.1093/molbev/mst189}, author = {Serohijos, A.W.R. and Eugene I. Shakhnovich} } @article {407871, title = {Merging molecular mechanism and evolution: theory and computation at the interface of biophysics and evolutionary population genetics}, journal = {Curr. Opin. Struct. Biol.}, volume = {26}, year = {2014}, month = {2014}, pages = {84 - 91}, abstract = {The variation among sequences and structures in nature is both determined by physical laws and by evolutionary history. However, these two factors are traditionally investigated by disciplines with different emphasis and philosophy-molecular biophysics on one hand and evolutionary population genetics in another. Here, we review recent theoretical and computational approaches that address the crucial need to integrate these two disciplines. We first articulate the elements of these approaches. Then, we survey their contribution to our mechanistic understanding of molecular evolution, the polymorphisms in coding region, the distribution of fitness effects (DFE) of mutations, the observed folding stability of proteins in nature, and the distribution of protein folds in genomes.}, url = {http://dx.doi.org/10.1016/j.sbi.2014.05.005}, author = {Serohijos, A.W.R. and Eugene I. Shakhnovich} } @article {191071, title = {Influenza A H1N1 Pandemic Strain Evolution - Divergence and the Potential for Antigenic Drift Variants}, journal = {PLOS One}, volume = {9}, year = {2014}, pages = {e93632}, abstract = {

The emergence of a novel A(H1N1) strain in 2009 was the first influenza pandemic of the genomic age, and unprecedented surveillance of the virus provides the opportunity to better understand the evolution of influenza. We examined changes in the nucleotide coding regions and the amino acid sequences of the hemagglutinin (HA), neuraminidase (NA), and nucleoprotein (NP) segments of the A(H1N1)pdm09 strain using publicly available data. We calculated the nucleotide and amino acid hamming distance from the vaccine strain A/California/07/2009 for each sequence. We also estimated Pepitope{\textendash}a measure of antigenic diversity based on changes in the epitope regions{\textendash}for each isolate. Finally, we compared our results to A(H3N2) strains collected over the same period. Our analysis found that the mean hamming distance for the HA protein of the A(H1N1)pdm09 strain increased from 3.6 (standard deviation [SD]: 1.3) in 2009 to 11.7 (SD: 1.0) in 2013, while the mean hamming distance in the coding region increased from 7.4 (SD: 2.2) in 2009 to 28.3 (SD: 2.1) in 2013. These trends are broadly similar to the rate of mutation in H3N2 over the same time period. However, in contrast to H3N2 strains, the rate of mutation accumulation has slowed in recent years. Our results are notable because, over the course of the study, mutation rates in H3N2 similar to that seen with A(H1N1)pdm09 led to the emergence of two antigenic drift variants. However, while there has been an H1N1 epidemic in North America this season, evidence to date indicates the vaccine is still effective, suggesting the epidemic is not due to the emergence of an antigenic drift variant. Our results suggest that more research is needed to understand how viral mutations are related to vaccine effectiveness so that future vaccine choices and development can be more predictive.

}, url = {http://dx.doi.org/10.1371/journal.pone.0093632}, author = {Eili Y. Klein and Adrian W.R. Serohijos and Jeong-Mo Choi and Eugene I. Shakhnovich and Andrew Pekosz} } @article {138206, title = {Catalysis of Protein Folding by Chaperones Accelerates Evolutionary Dynamics in Adapting Cell Populations}, journal = {PLOS Computational Biology}, volume = {9}, year = {2013}, month = {2013}, pages = {e1003269}, abstract = {Although molecular chaperones are essential components of protein homeostatic machinery, their mechanism of action and impact on adaptation and evolutionary dynamics remain controversial. Here we developed a physics-based ab initio multi-scale model of a living cell for population dynamics simulations to elucidate the effect of chaperones on adaptive evolution. The 6-loci genomes of model cells encode model proteins, whose folding and interactions in cellular milieu can be evaluated exactly from their genome sequences. A genotype-phenotype relationship that is based on a simple yet non-trivially postulated protein-protein interaction (PPI) network determines the cell division rate. Model proteins can exist in native and molten globule states and participate in functional and all possible promiscuous non-functional PPIs. We find that an active chaperone mechanism, whereby chaperones directly catalyze protein folding, has a significant impact on the cellular fitness and the rate of evolutionary dynamics, while passive chaperones, which just maintain misfolded proteins in soluble complexes have a negligible effect on the fitness. We find that by partially releasing the constraint on protein stability, active chaperones promote a deeper exploration of sequence space to strengthen functional PPIs, and diminish the non-functional PPIs. A key experimentally testable prediction emerging from our analysis is that down-regulation of chaperones that catalyze protein folding significantly slows down the adaptation dynamics.}, url = {http://dx.doi.org/10.1371/journal.pcbi.1003269}, author = {{\c C}etinba{\c s}, Murat and Eugene I. Shakhnovich} } @article {94321, title = {

Positively Selected Sites in Cetacean Myoglobins Contribute to Protein Stability

}, journal = {PLOS Computational Biology}, volume = {9}, year = {2013}, pages = {e1002929}, abstract = {Since divergence ~50 Ma ago from their terrestrial ancestors, cetaceans underwent a series of adaptations such as a ~10{\textendash}20 fold increase in myoglobin (Mb) concentration in skeletal muscle, critical for increasing oxygen storage capacity and prolonging dive time. Whereas the O2-binding affinity of Mbs is not significantly different among mammals (with typical oxygenation constants of ~0.8{\textendash}1.2 {\textmu}M-1), folding stabilities of cetacean Mbs are ~2{\textendash}4 kcal/mol higher than for terrestrial Mbs. Using ancestral sequence reconstruction, maximum likelihood and Bayesian tests to describe the evolution of cetacean Mbs, and experimentally calibrated computation of stability effects of mutations, we observe accelerated evolution in cetaceans and identify seven positively selected sites in Mb. Overall, these sites contribute to Mb stabilization with a conditional probability of 0.8. We observe a correlation between Mb folding stability and protein abundance, suggesting that a selection pressure for stability acts proportionally to higher expression. We also identify a major divergence event leading to the common ancestor of whales, during which major stabilization occurred. Most of the positively selected sites that occur later act against other destabilizing mutations to maintain stability across the clade, except for the shallow divers, where late stability relaxation occurs, probably due to the shorter aerobic dive limits of these species. The three main positively selected sites 66, 5, and 35 undergo changes that favor hydrophobic folding, structural integrity, and intra-helical hydrogen bonds.}, url = {http://dx.doi.org/10.1371/journal.pcbi.1002929}, author = {Pouria Dasmeh and Adrian W.R. Serohijos and Kasper P. Kepp and Eugene I. Shakhnovich} } @article {85501, title = {Protein Quality Control Acts on Folding Intermediate to Shape the effects of Mutations on Organismal Fitness}, journal = {Molecular Cell}, volume = {49}, year = {2013}, pages = {133 - 144}, abstract = {

What are the molecular properties of proteins that fall on the radar of protein quality control (PQC)? Here we mutate the\ E.\ coli{\textquoteright}s gene encoding dihydrofolate reductase (DHFR) and replace it with bacterial orthologous genes to determine how components of PQC modulate fitness effects of these genetic changes. We find that chaperonins GroEL/ES and protease Lon compete for binding to molten globule intermediate of DHFR, resulting in a peculiar symmetry in their action: overexpression of GroEL/ES and deletion of Lon both restore growth of deleterious DHFR mutants and most of the slow-growing orthologous DHFR strains. Kinetic steady-state modeling predicts and experimentation verifies that mutations affect fitness by shifting the flux balance in cellular milieu between protein production, folding, and degradation orchestrated by PQC through the interaction with folding intermediates.

}, url = {http://dx.doi.org/10.1016/j.molcel.2012.11.004}, author = {Shimon Bershtein and Wanmeng Mu and Adrian W.R. Serohijos and Jingwen Zhou and Eugene I. Shakhnovich} } @article {85441, title = {Highly Abundant Proteins Favor More Stable 3D Structures in Yeast}, journal = {Biophysical Journal}, volume = {104}, year = {2013}, pages = {L1 - L3}, abstract = {To understand the variation of protein sequences in nature, we need to reckon with evolutionary constraints that are biophysical, cellular, and ecological. Here, we show that under the global selection against protein misfolding, there exists a scaling among protein folding stability, protein cellular abundance, and effective population size. The specific scaling implies that the several-orders-of-magnitude range of protein abundances in the cell should leave imprints on extant protein structures, a prediction that is supported by our structural analysis of the yeast proteome.}, url = {http://dx.doi.org/10.1016/j.bpj.2012.11.3838}, author = {Adrian W.R. Serohijos and S. Y. Ryan Lee and Eugene I. Shakhnovich} } @article {86671, title = {The interface of protein structure, protein biophysics, and molecular evolution}, journal = {Protein Science}, volume = {21}, year = {2012}, month = {2012}, pages = {769 - 785}, abstract = {Abstract The interface of protein structural biology, protein biophysics, molecular evolution, and molecular population genetics forms the foundations for a mechanistic understanding of many aspects of protein biochemistry. Current efforts in interdisciplinary protein modeling are in their infancy and the state-of-the art of such models is described. Beyond the relationship between amino acid substitution and static protein structure, protein function, and corresponding organismal fitness, other considerations are also discussed. More complex mutational processes such as insertion and deletion and domain rearrangements and even circular permutations should be evaluated. The role of intrinsically disordered proteins is still controversial, but may be increasingly important to consider. Protein geometry and protein dynamics as a deviation from static considerations of protein structure are also important. Protein expression level is known to be a major determinant of evolutionary rate and several considerations including selection at the mRNA level and the role of interaction specificity are discussed. Lastly, the relationship between modeling and needed high-throughput experimental data as well as experimental examination of protein evolution using ancestral sequence resurrection and in vitro biochemistry are presented, towards an aim of ultimately generating better models for biological inference and prediction.}, url = {http://dx.doi.org/10.1002/pro.2071}, author = {Liberles, D. A. and Teichmann, S. A. and Bahar, I. and Bastolla, U. and Bloom, J. and Bornberg-Bauer, E. and Colwell, L.J. and de Koning, A.P. and Dokholyan, N.V. and Echave, J. and Elofsson, A. and Gerloff, D.L. and Goldstein, R.A. and Grahnen, J.A. and Holder, M.T. and Lakner, C. and Lartillot, N. and Lovell, S.C. and Naylor, G. and Perica, T. and Pollock, D.D. and Pupko, T. and Regan, L. and Roger, A. and Rubinstein, N. and Shakhnovich, E.I. and Sj{\"o}lander, K. and Sunyaev, S. and Teufel, A.I. and Thorne, J.L. and Thornton, J.W. and Weinreich, D.M. and Whelan, S.} } @article {86661, title = {Mutation Induced Extinction in Finite Populations: Lethal Mutagenesis and Lethal Isolation}, journal = {PloS Computational Biology}, volume = {8}, year = {2012}, month = {2012}, pages = {e1002609}, abstract = {Reproduction is inherently risky, in part because genomic replication can introduce new mutations that are usually deleterious toward fitness. This risk is especially severe for organisms whose genomes replicate {\textquotedblleft}semi-conservatively,{\textquotedblright} e.g. viruses and bacteria, where no master copy of the genome is preserved. Lethal mutagenesis refers to extinction of populations due to an unbearably high mutation rate (U), and is important both theoretically and clinically, where drugs can extinguish pathogens by increasing their mutation rate. Previous theoretical models of lethal mutagenesis assume infinite population size (N). However, in addition to high U, small N can accelerate extinction by strengthening genetic drift and relaxing selection. Here, we examine how the time until extinction depends jointly on N and U. We first analytically compute the mean time until extinction (τ) in a simplistic model where all mutations are either lethal or neutral. The solution motivates the definition of two distinct regimes: a survival phase and an extinction phase, which differ dramatically in both how τ scales with N and in the coefficient of variation in time until extinction. Next, we perform stochastic population-genetics simulations on a realistic fitness landscape that both (i) features an epistatic distribution of fitness effects that agrees with experimental data on viruses and (ii) is based on the biophysics of protein folding. More specifically, we assume that mutations inflict fitness penalties proportional to the extent that they unfold proteins. We find that decreasing N can cause phase transition-like behavior from survival to extinction, which motivates the concept of {\textquotedblleft}lethal isolation.{\textquotedblright} Furthermore, we find that lethal mutagenesis and lethal isolation interact synergistically, which may have clinical implications for treating infections. Broadly, we conclude that stably folded proteins are only possible in ecological settings that support sufficiently large populations.}, url = {http://dx.doi.org/10.1371/journal.pcbi.1002609}, author = {C. Scott Wylie and Eugene I. Shakhnovich} } @article {86676, title = {Is there an en route folding intermediate for cold shock proteins?}, journal = {Protein Science}, volume = {21}, year = {2012}, month = {2012}, pages = {677 - 685}, abstract = {Cold shock proteins (Csps) play an important role in cold shock response of a diverse number of organisms ranging from bacteria to humans. Numerous studies of the Csp from various species showed that a two-state folding mechanism is conserved and the transition state (TS) appears to be very compact. However, the atomic details of the folding mechanism of Csp remain unclear. This study presents the folding mechanism of Csp in atomic detail using an all-atom Go model-based simulations. Our simulations predict that there may exist an en route intermediate, in which β strands 1-2-3 are well ordered and the contacts between β1 and β4 are almost developed. Such an intermediate might be too unstable to be detected in the previous fluorescence energy transfer experiments. The transition state ensemble has been determined from the P(fold) analysis and the TS appears even more compact than the intermediate state.}, url = {http://dx.doi.org/10.1002/pro.2053}, author = {Huang, Lei and Eugene I. Shakhnovich} } @article {86656, title = {Robustness of atomistic Go models in predicting native-like folding intermediates}, journal = {Journal of Chemical Physics}, volume = {137}, year = {2012}, month = {2012}, pages = {085102}, abstract = {Go models are exceedingly popular tools in computer simulations of protein folding. These models are native-centric, i.e., they are directly constructed from the protein{\textquoteright}s native structure. Therefore, it is important to understand up to which extent the atomistic details of the native structure dictate the folding behavior exhibited by Go models. Here we address this challenge by performing exhaustive discrete molecular dynamics simulations of a Go potential combined with a full atomistic protein representation. In particular, we investigate the robustness of this particular type of Go models in predicting the existence of intermediate states in protein folding. We focus on the N47G mutational form of the Spc-SH3 folding domain (x-ray structure) and compare its folding pathway with that of alternative native structures produced in silico. Our methodological strategy comprises equilibrium folding simulations, structural clustering, and principal component analysis.}, url = {http://dx.doi.org/10.1063/1.4747492}, author = {Estacio, S. G. and Fernandes, C. S. and Krobath, H. and Faisca, P. F. N. and Shakhnovich, E.I.} } @article {86136, title = {Identification of a Conserved Aggregation-Prone Intermediate State in the Folding Pathways of Spc-SH3 Amyloidogenic Variants}, journal = {Journal of Molecular Biology}, volume = {422}, year = {2012}, month = {2012/10/5/}, pages = {705 - 722}, abstract = {We compared the folding pathways of selected mutational variants of the α-spectrin SH3 domain (Spc-SH3) by using a continuum model that combines a full atomistic protein representation with the Go potential. Experimental data show that the N47G mutant shows very little tendency to aggregate while the N47A and triple mutant D48G(2Y) are both amyloidogenic, with the latter being clearly more aggregation prone. We identified a strikingly similar native-like folding intermediate across the three mutants, in which strand β1 is totally unstructured and more than half of the major hydrophobic core residues are highly solvent exposed. Results from extensive docking simulations show that the ability of the intermediates to dimerize is largely driven by strand β1 and is consistent with the in vitro aggregation behavior reported for the corresponding mutants. They further suggest that residues 44 and 53, which are key players in the nucleation{\textendash}condensation mechanism of folding, are also important triggers of the aggregation process.}, keywords = {aggregation, amyloids, intermediate states, molecular dynamics, protein folding}, isbn = {0022-2836}, url = {http://dx.doi.org/10.1016/j.jmb.2012.06.020}, author = {Krobath, H. and Est{\'a}cio, S.G. and Fa{\'\i}sca, P.F.N. and Shakhnovich, E.I.} } @article {86131, title = {Collapse of Unfolded Proteins in a Mixture of Denaturants}, journal = {J. Am. Chem. Soc.}, volume = {134}, year = {2012}, note = {doi: 10.1021/ja3031505}, month = {2012}, pages = {18266 - 18274}, publisher = {American Chemical Society}, abstract = {Both urea and guanidinium chloride (GdmCl) are frequently used as protein denaturants. Given that proteins generally adopt extended or unfolded conformations in either aqueous urea or GdmCl, one might expect that the unfolded protein chains will remain or become further extended due to the addition of another denaturant. However, a collapse of denatured proteins is revealed using atomistic molecular dynamics simulations when a mixture of denaturants is used. Both hen egg-white lysozyme and protein L are found to undergo collapse in the denaturant mixture. The collapse of the protein conformational ensembles is accompanied by a decreased solubility and increased non-native self-interactions of hydrophobic residues in the urea/GdmCl mixture. The increase of non-native interactions rather than the native contacts indicates that the proteins experience a simple collapse transition from the fully denatured states. During the protein collapse, the relatively stronger denaturant GdmCl displays a higher tendency to be absorbed onto the protein surface due to their stronger electrostatic interactions with proteins. At the same time, urea molecules also accumulate near the protein surface, resulting in an enhanced ?local crowding? for the protein near its first solvation shell. This rearrangement of denaturants near the protein surface and crowded local environment induce the protein collapse, mainly by burying their hydrophobic residues. These findings from molecular simulations are then further explained by a simple analytical model based on statistical mechanics.}, isbn = {0002-7863}, url = {http://dx.doi.org/10.1021/ja3031505}, author = {Xia, Zhen and Das, Payel and Eugene I. Shakhnovich and Zhou, Ruhong} } @article {85586, title = {Protein biophysics explains why highly abundant proteins evolve slowly}, journal = {Cell Reports}, volume = {2}, year = {2012}, pages = {249 - 256}, abstract = {The consistent observation across all kingdoms of life that highly abundant proteins evolve slowly demonstrates that cellular abundance is a key determinant of protein evolutionary rate. However, other empirical findings, such as the broad distribution of evolutionary rates, suggest that additional variables determine the rate of protein evolution. Here, we report that under the global selection against the cytotoxic effects of misfolded proteins, folding stability (ΔG), simultaneous with abundance, is a causal variable of evolutionary rate. Using both theoretical analysis and multiscale simulations, we demonstrate that the anticorrelation between the premutation ΔG and the arising mutational effect (ΔΔG), purely biophysical in origin, is a necessary requirement for abundance{\textendash}evolutionary rate covariation. Additionally, we predict and demonstrate in bacteria that the strength of abundance{\textendash}evolutionary rate correlation depends on the divergence time separating reference genomes. Altogether, these results highlight the intrinsic role of protein biophysics in the emerging universal patterns of molecular evolution.}, url = {http://dx.doi.org/10.1016/j.celrep.2012.06.022}, author = {Adrian W.R. Serohijos and Zilvinas Rimas and Eugene I. Shakhnovich} } @article {85511, title = {Soluble oligomerization provides a beneficial fitness effect on destabilizing mutations}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {109}, year = {2012}, pages = {4857 - 4862}, abstract = {Mutations create the genetic diversity on which selective pressures can act, yet also create structural instability in proteins. How, then, is it possible for organisms to ameliorate mutation-induced perturbations of protein stability while maintaining biological fitness and gaining a selective advantage? Here we used site-specific chromosomal mutagenesis to introduce a selected set of mostly destabilizing mutations into folA{\textemdash}an essential chromosomal gene of Escherichia coli encoding dihydrofolate reductase (DHFR){\textemdash}to determine how changes in protein stability, activity, and abundance affect fitness. In total, 27 E. coli strains carrying mutant DHFR were created. We found no significant correlation between protein stability and its catalytic activity nor between catalytic activity and fitness in a limited range of variation of catalytic activity observed in mutants. The stability of these mutants is strongly correlated with their intracellular abundance, suggesting that protein homeostatic machinery plays an active role in maintaining intracellular concentrations of proteins. Fitness also shows a significant correlation with intracellular abundance of soluble DHFR in cells growing at 30 {\textdegree}C. At 42 {\textdegree}C, the picture was mixed, yet remarkable: A few strains carrying mutant DHFR proteins aggregated, rendering them nonviable, but, intriguingly, the majority exhibited fitness higher than wild type. We found that mutational destabilization of DHFR proteins in E. coli is counterbalanced at 42 {\textdegree}C by their soluble oligomerization, thereby restoring structural stability and protecting against aggregation.}, url = {http://dx.doi.org/10.1073/pnas.1118157109}, author = {Shimon Bershtein and Wanmeng Mu and Eugene I. Shakhnovich} } @article {86691, title = {Diversity Against Adversity: How Adaptive Immune System Evolves Potent Antibodies}, journal = {Journal of Statistical Physics}, volume = {144}, year = {2011}, month = {2011}, pages = {241 - 267}, abstract = {Adaptive immunity is an amazing mechanism, whereby new protein functions{\textemdash}affinity of antibodies (Immunoglobulins) to new antigens{\textemdash}evolve through mutation and selection in a matter of a few days. Despite numerous experimental studies, the fundamental physical principles underlying immune response are still poorly understood. In considerable departure from past approaches, here, we propose a microscopic multiscale model of adaptive immune response, which consists of three essential players: The host cells, viruses, and B-cells in Germinal Centers (GC). Each moiety carries a genome, which encodes proteins whose stability and interactions are determined from their sequences using laws of Statistical Mechanics, providing an exact relationship between genomic sequences and strength of interactions between pathogens and antibodies and antibodies and host proteins (autoimmunity). We find that evolution of potent antibodies (the process known as Affinity Maturation (AM)) is a delicate balancing act, which has to reconcile the conflicting requirements of protein stability, lack of autoimmunity, and high affinity of antibodies to incoming antigens. This becomes possible only when antibody producing B cells elevate their mutation rates (process known as Somatic Hypermutation (SHM)) to fall into a certain range{\textemdash}not too low to find potency increasing mutations but not too high to destroy stable Immunoglobulins and/or already achieved affinity. Potent antibodies develop through clonal expansion of initial B cells expressing marginally potent antibodies followed by their subsequent affinity maturation through mutation and selection. As a result, in each GC the population of mature potent Immunoglobulins is monoclonal being ancestors of a single cell from initial (germline) pool. We developed a simple analytical theory, which provides further rationale to our findings. The model and theory reveal the molecular factors that determine the efficiency of affinity maturation, thereby providing insight into the variability of the immune response to cytopathic viruses (the direct response by germline antibodies) and poorly cytopathic viruses (a crucial role of SHM in the response).}, url = {http://dx.doi.org/10.1007/s10955-011-0241-8}, author = {Muyoung Heo and Zeldovich, Konstantin B. and Eugene I. Shakhnovich} } @article {86701, title = {The ensemble folding kinetics of the FBP28 WW domain revealed by an all-atom Monte Carlo simulation in a knowledge-based potential}, journal = {Proteins - Structure, function and bioinformatics}, volume = {79}, year = {2011}, month = {2011}, pages = {1704 - 1714}, abstract = {In this work, we apply a detailed all-atom model with a transferable knowledge-based potential to study the folding kinetics of Formin-Binding protein, FBP28, which is a canonical three-stranded β-sheet WW domain. Replica exchange Monte Carlo simulations starting from random coils find native-like (Cα RMSD of 2.68 {\r A}) lowest energy structure. We also study the folding kinetics of FBP28 WW domain by performing a large number of ab initio Monte Carlo folding simulations. Using these trajectories, we examine the order of formation of two β-hairpins, the folding mechanism of each individual β-hairpin, and transition state ensemble (TSE) of FBP28 WW domain and compare our results with experimental data and previous computational studies. To obtain detailed structural information on the folding dynamics viewed as an ensemble process, we perform a clustering analysis procedure based on graph theory. Further, a rigorous P(fold) analysis is used to obtain representative samples of the TSEs showing good quantitative agreement between experimental and simulated Φ values. Our analysis shows that the turn structure between first and second β strands is a partially stable structural motif that gets formed before entering the TSE in FBP28 WW domain and there exist two major pathways for the folding of FBP28 WW domain, which differ in the order and mechanism of hairpin formation.}, url = {http://dx.doi.org/10.1002/prot.22993}, author = {Xu, Jiabin and Huang, Lei and Eugene I. Shakhnovich} } @article {86696, title = {Multi-Scale Sequence Correlations Increase Proteome Structural Disorder and Promiscuity}, journal = {Journal of Molecular Biology}, volume = {409}, year = {2011}, month = {2011}, pages = {439 - 449}, abstract = {Numerous experiments demonstrate a high level of promiscuity and structural disorder in organismal proteomes. Here, we ask the question what makes a protein promiscuous, that is, prone to nonspecific interactions, and structurally disordered. We predict that multi-scale correlations of amino acid positions within protein sequences statistically enhance the propensity for promiscuous intra- and inter-protein binding. We show that sequence correlations between amino acids of the same type are statistically enhanced in structurally disordered proteins and in hubs of organismal proteomes. We also show that structurally disordered proteins possess a significantly higher degree of sequence order than structurally ordered proteins. We develop an analytical theory for this effect and predict the robustness of our conclusions with respect to the amino acid composition and the form of the microscopic potential between the interacting sequences. Our findings have implications for understanding molecular mechanisms of protein aggregation diseases induced by the extension of sequence repeats.}, url = {http://dx.doi.org/10.1016/j.jmb.2011.03.056}, author = {Afek, Ariel and Eugene I. Shakhnovich and Lukatsky, David B.} } @article {86686, title = {Sequence correlations shape protein promiscuity}, journal = {Journal of Chemical Physics}, volume = {135}, year = {2011}, month = {2011}, pages = {065104}, abstract = {We predict analytically that diagonal correlations of amino acid positions within protein sequences statistically enhance protein propensity for nonspecific binding. We use the term {\textquotedblleft}promiscuity{\textquotedblright} to describe such nonspecific binding. Diagonal correlations represent statistically significant repeats of sequence patterns where amino acids of the same type are clustered together. The predicted effect is qualitatively robust with respect to the form of the microscopic interaction potentials and the average amino acid composition. Our analytical results provide an explanation for the enhanced diagonal correlations observed in hubs of eukaryotic organismal proteomes [J. Mol. Biol. 409, 439 (2011), 10.1016/j.jmb.2011.03.056]. We suggest experiments that will allow direct testing of the predicted effect.}, url = {http://dx.doi.org/10.1063/1.3624332}, author = {Lukatsky, David B. and Afek, Ariel and Eugene I. Shakhnovich} } @article {86681, title = {Thermodynamics and Kinetics of the Hairpin Ribozyme from Atomistic Folding/Unfolding Simulations}, journal = {Journal of Molecular Biology}, volume = {411}, year = {2011}, month = {2011}, pages = {1128 - 1144}, abstract = {We report a set of atomistic folding/unfolding simulations for the hairpin ribozyme using a Monte Carlo algorithm. The hairpin ribozyme folds in solution and catalyzes self-cleavage or ligation via a specific two-domain structure. The minimal active ribozyme has been studied extensively, showing stabilization of the active structure by cations and dynamic motion of the active structure. Here, we introduce a simple model of tertiary-structure formation that leads to a phase diagram for the RNA as a function of temperature and tertiary-structure strength. We then employ this model to capture many folding/unfolding events and to examine the transition-state ensemble (TSE) of the RNA during folding to its active "docked" conformation. The TSE is compact but with few tertiary interactions formed, in agreement with single-molecule dynamics experiments. To compare with experimental kinetic parameters, we introduce a novel method to benchmark Monte Carlo kinetic parameters to docking/undocking rates collected over many single molecular trajectories. We find that topology alone, as encoded in a biased potential that discriminates between secondary and tertiary interactions, is sufficient to predict the thermodynamic behavior and kinetic folding pathway of the hairpin ribozyme. This method should be useful in predicting folding transition states for many natural or man-made RNA tertiary structures.}, url = {10.1016/j.jmb.2011.06.042}, author = {Nivon, Lucas G. and Eugene I. Shakhnovich} } @article {86141, title = {Protein folding: To knot or not to knot?}, journal = {Nature Materials}, volume = {10}, year = {2011}, pages = {84 - 86}, abstract = {A knot-containing protein is found to fold reversibly at biologically relevant timescales despite not having naturally evolved for this ability.}, url = {http://dx.doi.org/10.1038/nmat2953}, author = {Eugene I. Shakhnovich} } @article {85596, title = {Topology of protein interaction network shapes protein abundances and strengths of their functional and nonspecific interactions}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {108}, year = {2011}, pages = {4258 - 4263}, abstract = {How do living cells achieve sufficient abundances of functional protein complexes while minimizing promiscuous nonfunctional interactions? Here we study this problem using a first-principle model of the cell whose phenotypic traits are directly determined from its genome through biophysical properties of protein structures and binding interactions in a crowded cellular environment. The model cell includes three independent prototypical pathways, whose topologies of protein{\textendash}protein interaction (PPI) subnetworks are different, but whose contributions to the cell fitness are equal. Model cells evolve through genotypic mutations and phenotypic protein copy number variations. We found a strong relationship between evolved physical{\textendash}chemical properties of protein interactions and their abundances due to a {\textquotedblleft}frustration{\textquotedblright} effect: Strengthening of functional interactions brings about hydrophobic interfaces, which make proteins prone to promiscuous binding. The balancing act is achieved by lowering concentrations of hub proteins while raising solubilities and abundances of functional monomers. On the basis of these principles we generated and analyzed a possible realization of the proteome-wide PPI network in yeast. In this simulation we found that high-throughput affinity capture{\textendash}mass spectroscopy experiments can detect functional interactions with high fidelity only for high-abundance proteins while missing most interactions for low-abundance proteins.}, url = {http://dx.doi.org/10.1073/pnas.1009392108}, author = {Muyoung Heo and Sergei Maslov and Eugene Shakhnovich} } @article {85591, title = {A biophysical protein folding model accounts for most mutational fitness effects in viruses}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {108}, year = {2011}, pages = {9916 - 9921}, abstract = {Fitness effects of mutations fall on a continuum ranging from lethal to deleterious to beneficial. The distribution of fitness effects (DFE) among random mutations is an essential component of every evolutionary model and a mathematical portrait of robustness. Recent experiments on five viral species all revealed a characteristic bimodal-shaped DFE featuring peaks at neutrality and lethality. However, the phenotypic causes underlying observed fitness effects are still unknown and presumably, are thought to vary unpredictably from one mutation to another. By combining population genetics simulations with a simple biophysical protein folding model, we show that protein thermodynamic stability accounts for a large fraction of observed mutational effects. We assume that moderately destabilizing mutations inflict a fitness penalty proportional to the reduction in folded protein, which depends continuously on folding free energy (ΔG). Most mutations in our model affect fitness by altering ΔG, whereas based on simple estimates, \~{}10\% abolish activity and are unconditionally lethal. Mutations pushing ΔG > 0 are also considered lethal. Contrary to neutral network theory, we find that, in mutation/selection/drift steady state, high mutation rates (m) lead to less stable proteins and a more dispersed DFE (i.e., less mutational robustness). Small population size (N) also decreases stability and robustness. In our model, a continuum of nonlethal mutations reduces fitness by \~{}2\% on average, whereas \~{}10{\textendash}35\% of mutations are lethal depending on N and m. Compensatory mutations are common in small populations with high mutation rates. More broadly, we conclude that interplay between biophysical and population genetic forces shapes the DFE.}, url = {http://dx.doi.org/10.1073/pnas.1017572108}, author = {C. Scott Wylie and Eugene I. Shakhnovich} } @article {86716, title = {Interplay between Pleiotropy and Secondary Selection Determines Rise and Fall of Mutators in Stress Response}, journal = {PLoS Computational Biology}, volume = {6}, year = {2010}, month = {2010}, pages = {e1000710}, abstract = {Mutators are clones whose mutation rate is about two to three orders of magnitude higher than the rate of wild-type clones and their roles in adaptive evolution of asexual populations have been controversial. Here we address this problem by using an ab initio microscopic model of living cells, which combines population genetics with a physically realistic presentation of protein stability and protein-protein interactions. The genome of model organisms encodes replication controlling genes (RCGs) and genes modeling the mismatch repair (MMR) complexes. The genotype-phenotype relationship posits that the replication rate of an organism is proportional to protein copy numbers of RCGs in their functional form and there is a production cost penalty for protein overexpression. The mutation rate depends linearly on the concentration of homodimers of MMR proteins. By simulating multiple runs of evolution of populations under various environmental stresses{\textemdash}stationary phase, starvation or temperature-jump{\textemdash}we find that adaptation most often occurs through transient fixation of a mutator phenotype, regardless of the nature of stress. By contrast, the fixation mechanism does depend on the nature of stress. In temperature jump stress, mutators take over the population due to loss of stability of MMR complexes. In contrast, in starvation and stationary phase stresses, a small number of mutators are supplied to the population via epigenetic stochastic noise in production of MMR proteins (a pleiotropic effect), and their net supply is higher due to reduced genetic drift in slowly growing populations under stressful environments. Subsequently, mutators in stationary phase or starvation hitchhike to fixation with a beneficial mutation in the RCGs, (second order selection) and finally a mutation stabilizing the MMR complex arrives, returning the population to a non-mutator phenotype. Our results provide microscopic insights into the rise and fall of mutators in adapting finite asexual populations.}, url = {http://dx.doi.org/journal.pcbi.1000710}, author = {Muyoung Heo and Eugene I. Shakhnovich} } @article {86711, title = {Optimality of mutation and selection in germinal centers}, journal = {PLoS Computational Biology}, volume = {6}, year = {2010}, month = {2010}, pages = {e1000800}, abstract = {The population dynamics theory of B cells in a typical germinal center could play an important role in revealing how affinity maturation is achieved. However, the existing models encountered some conflicts with experiments. To resolve these conflicts, we present a coarse-grained model to calculate the B cell population development in affinity maturation, which allows a comprehensive analysis of its parameter space to look for optimal values of mutation rate, selection strength, and initial antibody-antigen binding level that maximize the affinity improvement. With these optimized parameters, the model is compatible with the experimental observations such as the ~100-fold affinity improvements, the number of mutations, the hypermutation rate, and the {\textquotedblleft}all or none{\textquotedblright} phenomenon. Moreover, we study the reasons behind the optimal parameters. The optimal mutation rate, in agreement with the hypermutation rate in vivo, results from a tradeoff between accumulating enough beneficial mutations and avoiding too many deleterious or lethal mutations. The optimal selection strength evolves as a balance between the need for affinity improvement and the requirement to pass the population bottleneck. These findings point to the conclusion that germinal centers have been optimized by evolution to generate strong affinity antibodies effectively and rapidly. In addition, we study the enhancement of affinity improvement due to B cell migration between germinal centers. These results could enhance our understanding of the functions of germinal centers.}, url = {http://dx.doi.org/10.1371/journal.pcbi.1000800}, author = {Jingshan Zhang and Eugene I. Shakhnovich} } @article {86721, title = {Stochastic Switching in Gene Networks Can Occur by a Single-Molecule Event or Many Molecular Steps}, journal = {Journal of Molecular Biology}, volume = {396}, year = {2010}, month = {2010}, pages = {230 - 244}, abstract = {Due to regulatory feedback, biological networks can exist stably in multiple states, leading to heterogeneous phenotypes among genetically identical cells. Random fluctuations in protein numbers, tuned by specific molecular mechanisms, have been hypothesized to drive transitions between these different states. We develop a minimal theoretical framework to analyze the limits of switching in terms of simple experimental parameters. Our model identifies and distinguishes between two distinct molecular mechanisms for generating stochastic switches. In one class of switches, the stochasticity of a single-molecule event, a specific and rare molecular reaction, directly controls the macroscopic change in a cell{\textquoteright}s state. In the second class, no individual molecular event is significant, and stochasticity arises from the propagation of biochemical noise through many molecular pathways and steps. As an example, we explore switches based on protein-DNA binding fluctuations and predict relations between transcription factor kinetics, absolute switching rate, robustness, and efficiency that differentiate between switching by single-molecule events or many molecular steps. Finally, we apply our methods to recent experimental data on switching in Escherichia coli lactose metabolism, providing quantitative interpretations of a single-molecule switching mechanism.}, url = {http://dx.doi.org/10.1016/j.jmb.2009.11.035}, author = {Choi, P.J. and Xie, X.S. and Shakhnovich, E.I.} } @article {86706, title = {Non-native interactions plan an effective role in protein folding dynamics}, journal = {Protein Science}, volume = {19}, year = {2010}, month = {2010}, pages = {2196 - 2209}, abstract = {Systematic Monte Carlo simulations of simple lattice models show that the final stage of protein folding is an ordered process where native contacts get locked (i.e., the residues come into contact and remain in contact for the duration of the folding process) in a well-defined order. The detailed study of the folding dynamics of protein-like sequences designed as to exhibit different contact energy distributions, as well as different degrees of sequence optimization (i.e., participation of non-native interactions in the folding process), reveals significant differences in the corresponding locking scenarios{\textemdash}the collection of native contacts and their average locking times, which are largely ascribable to the dynamics of non-native contacts. Furthermore, strong evidence for a positive role played by non-native contacts at an early folding stage was also found. Interestingly, for topologically simple target structures, a positive interplay between native and non-native contacts is observed also toward the end of the folding process, suggesting that non-native contacts may indeed affect the overall folding process. For target models exhibiting clear two-state kinetics, the relation between the nucleation mechanism of folding and the locking scenario is investigated. Our results suggest that the stabilization of the folding transition state can be achieved through the establishment of a very small network of native contacts that are the first to lock during the folding process.}, url = {http://dx.doi.org/10.1002/pro.498}, author = {Faisca, P. F. N. and Nunes, A. and Travasso, R.D.M. and Shakhnovich, E.I.} } @article {86146, title = {Thermal Adaptation of Viruses and Bacteria}, journal = {Biophysical Journal}, volume = {98}, year = {2010}, month = {2010/4/7/}, pages = {1109 - 1118}, abstract = {A previously established multiscale population genetics model posits that fitness can be inferred from the physical properties of proteins under the physiological assumption that a loss of stability by any protein confers the lethal phenotype to an~organism. Here, we develop this model further by positing that replication rate (fitness) of a bacterial or viral strain directly depends on the copy number of folded proteins, which determine its replication rate. Using this model, and both numerical and analytical approaches, we studied the adaptation process of bacteria and viruses at varied environmental temperatures. We found that a broad distribution of protein stabilities observed in the model and in experiment is the key determinant of thermal response for viruses and bacteria. Our results explain most of the earlier experimental observations: the striking asymmetry of thermal response curves; the absence of evolutionary tradeoff, which was expected but not found in experiments; correlation between denaturation temperature for several protein families and the optimal growth temperature of their carrier organisms; and proximity of bacterial or viral optimal growth temperatures to their evolutionary temperatures. Our theory quantitatively and with high accuracy described thermal response curves for 35 bacterial species using, for each species, only two adjustable parameters{\textemdash}the number of rate-determining genes and the energy barrier for metabolic reactions.}, isbn = {0006-3495}, url = {http://dx.doi.org/10.1016/j.bpj.2009.11.048}, author = {Chen, Peiqiu and Eugene I. Shakhnovich} } @article {85516, title = {De novo design: balancing novelty and confined chemical space}, journal = {Expert Opin Drug Discov.}, volume = {5}, year = {2010}, pages = {789 - 812}, abstract = {MPORTANCE OF THE FIELD: De novo drug design serves as a tool for the discovery of new ligands for macromolecular targets as well as optimization of known ligands. Recently developed tools aim to address the multi-objective nature of drug design in an unprecedented manner. AREAS COVERED IN THIS REVIEW: This article discusses recent advances in de novo drug design programs and accessory programs used to evaluate compounds post-generation. WHAT THE READER WILL GAIN: The reader is introduced to the challenges inherent in de novo drug design and will become familiar with current trends in de novo design. Furthermore, the reader will be better prepared to assess the value of a tool, and be equipped to design more elegant tools in the future. TAKE HOME MESSAGE: De novo drug design can assist in the efficient discovery of new compounds with a high affinity for a given target. The inclusion of existing chemoinformatic methods with current structure-based de novo design tools provides a means of enhancing the therapeutic value of these generated compounds}, url = {http://dx.doi.org/10.1517/17460441.2010.497534}, author = {Peter S. Kutchukian and Eugene I. Shakhnovich} } @article {101376, title = {SDR: a database of predicted specificity-determining residues in proteins}, journal = {Nucl. Acids Res.}, volume = {37}, year = {2009}, month = {2009}, pages = {D191 - D194}, abstract = {The specificity-determining residue database (SDR database) presents residue positions where mutations are predicted to have changed protein function in large protein families. Because the database pre-calculates predictions on existing protein sequence alignments, users can quickly find the predictions by selecting the appropriate protein family or searching by protein sequence. Predictions can be used to guide mutagenesis or to gain a better understanding of specificity changes in a protein family. The database is available on the web at http://paradox.harvard.edu/sdr.}, url = {http://dx.doi.org/10.1093/nar/gkn716}, author = {Donald, Jason E. and Eugene I. Shakhnovich} } @article {101276, title = {Slowly Replicating Lytic Viruses: Pseudolysogenic Persistence and Within-Host Competition}, journal = {Phys. Rev. Lett.}, volume = {102}, year = {2009}, month = {2009}, pages = {178103}, abstract = {We study the population dynamics of lytic viruses which replicate slowly in dividing host cells within an organism or cell culture, and find a range of viral replication rates that allows viruses to persist, avoiding extinction of host cells or dilution of viruses at too rapid or too slow viral replication. For the within-host competition between viral strains with different replication rates, a strain with a "stable{\textquotedblright} replication rate in the persistence range could outcompete another strain. However, when strains with higher and lower than the stable value replication rates are both present, competition between strains does not result in the dominance of one strain, but in their coexistence.}, url = {http://dx.doi.org/10.1103/PhysRevLett.102.178103}, author = {Jingshan Zhang and Eugene I. Shakhnovich} } @article {101166, title = {All-Atom Model for Stabilization of α-Helical Structure in Peptides by Hydrocarbon Staples}, journal = {J. Am. Chem. Soc.}, volume = {131}, year = {2009}, month = {2009}, pages = {4622 - 4627}, abstract = {Recent work has shown that the incorporation of an all-hydrocarbon ?staple? into peptides can greatly increase their α-helix propensity, leading to an improvement in pharmaceutical properties such as proteolytic stability, receptor affinity, and cell permeability. Stapled peptides thus show promise as a new class of drugs capable of accessing intractable targets such as those that engage in intracellular protein?protein interactions. The extent of α-helix stabilization provided by stapling has proven to be substantially context dependent, requiring cumbersome screening to identify the optimal site for staple incorporation. In certain cases, a staple encompassing one turn of the helix (attached at residues i and i+4) furnishes greater helix stabilization than one encompassing two turns (i,i+7 staple), which runs counter to expectation based on polymer theory. These findings highlight the need for a more thorough understanding of the forces that underlie helix stabilization by hydrocarbon staples. Here we report all-atom Monte Carlo folding simulations comparing unmodified peptides derived from RNase A and BID BH3 with various i,i+4 and i,i+7 stapled versions thereof. The results of these simulations were found to be in quantitative agreement with experimentally determined helix propensities. We also discovered that staples can stabilize quasi-stable decoy conformations, and that the removal of these states plays a major role in determining the helix stability of stapled peptides. Finally, we critically investigate why our method works, exposing the underlying physical forces that stabilize stapled peptides.}, url = {http://dx.doi.org/10.1021/ja805037p}, author = {Peter S. Kutchukian and Yang, Jae Shick and Verdine, Gregory L. and Eugene I. Shakhnovich} } @article {101126, title = {DNA unzipping phase diagram calculated via replica theory}, journal = {Phys. Rev. E}, volume = {79}, year = {2009}, month = {2009}, pages = {051923}, abstract = {We show how single-molecule unzipping experiments can provide strong evidence that the zero-force melting transition of long molecules of natural dsDNA should be classified as a phase transition of the higher-order type (continuous). Toward this end, we study a statistical-mechanics model for the fluctuating structure of a long molecule of dsDNA, and compute the equilibrium phase diagram for the experiment in which the molecule is unzipped under applied force. We consider a perfect-matching dsDNA model, in which the loops are volume-excluding chains with arbitrary loop exponent c. We include stacking interactions, hydrogen bonds, and main-chain entropy. We include sequence heterogeneity at the level of random sequences; in particular, there is no correlation in the base-pairing (bp) energy from one sequence position to the next. We present heuristic arguments to demonstrate that the low-temperature macrostate does not exhibit degenerate ergodicity breaking. We use this claim to understand the results of our replica-theoretic calculation of the equilibrium properties of the system. As a function of temperature, we obtain the minimal force at which the molecule separates completely. This critical-force curve is a line in the temperature-force phase diagram that marks the regions where the molecule exists primarily as a double helix versus the region where the molecule exists as two separate strands. We compare our random-sequence model to magnetic tweezer experiments performed on the 48 502 bp genome of bacteriophage λ. We find good agreement with the experimental data, which is restricted to temperatures between 24 and 50 {\textdegree}C. At higher temperatures, the critical-force curve of our random-sequence model is very different for that of the homogeneous-sequence version of our model. For both sequence models, the critical force falls to zero at the melting temperature Tc like |T-Tc|α. For the homogeneous-sequence model, α=1/2 almost exactly, while for the random-sequence model, α≈0.9. Importantly, the shape of the critical-force curve is connected, via our theory, to the manner in which the helix fraction falls to zero at Tc. The helix fraction is the property that is used to classify the melting transition as a type of phase transition. In our calculation, the shape of the critical-force curve holds strong evidence that the zero-force melting transition of long natural dsDNA should be classified as a higher-order (continuous) phase transition. Specifically, the order is 3rd or greater.}, url = {http://dx.doi.org/10.1103/PhysRevE.79.051923}, author = {Roland, C. Brian and Hatch, Kristi Adamson and Prentiss, Mara and Eugene I. Shakhnovich} } @article {101056, title = {Emergence of species in evolutionary "simulated annealing{\textquotedblright}}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {106}, year = {2009}, month = {2009}, pages = {1869 - 1874}, abstract = {Which factors govern the evolution of mutation rates and emergence of species? Here, we address this question by using a first principles model of life where population dynamics of asexual organisms is coupled to molecular properties and interactions of proteins encoded in their genomes. Simulating evolution of populations, we found that fitness increases in punctuated steps via epistatic events, leading to formation of stable and functionally interacting proteins. At low mutation rates, species form populations of organisms tightly localized in sequence space, whereas at higher mutation rates, species are lost without an apparent loss of fitness. However, when mutation rate was a selectable trait, the population initially maintained high mutation rate until a high fitness level was reached, after which organisms with low mutation rates are gradually selected, with the population eventually reaching mutation rates comparable with those of modern DNA-based organisms. This study shows that the fitness landscape of a biophysically realistic system is extremely complex, with huge number of local peaks rendering adaptation dynamics to be a glass-like process. On a more practical level, our results provide a rationale to experimental observations of the effect of mutation rate on fitness of populations of asexual organisms.}, url = {http://dx.doi.org/10.1073/pnas.0809852106}, author = {Muyoung Heo and Kang, Louis and Eugene I. Shakhnovich} } @article {101006, title = {Thymic Selection of T-Cell Receptors as an Extreme Value Problem}, journal = {Phys. Rev. Lett.}, volume = {103}, year = {2009}, month = {2009}, pages = {068103}, abstract = {T lymphocytes (T cells) orchestrate adaptive immune responses upon activation. T-cell activation requires sufficiently strong binding of T-cell receptors on their surface to short peptides (p) derived from foreign proteins, which are bound to major histocompatibility gene products (displayed on antigen-presenting cells). A diverse and self-tolerant T-cell repertoire is selected in the thymus. We map thymic selection processes to an extreme value problem and provide an analytic expression for the amino acid compositions of selected T-cell receptors (which enable its recognition functions).}, url = {http://dx.doi.org/10.1103/PhysRevLett.103.068103}, author = {Ko{\v s}mrlj, Andrej and Chakraborty, Arup K. and Kardar, Mehran and Eugene I. Shakhnovich} } @article {100381, title = {Lethal Mutagenesis in Viruses and Bacteria}, journal = {Genetics}, volume = {183}, year = {2009}, month = {2009}, pages = {639 - 650}, abstract = {In this work we study how mutations that change physical properties of cell proteins (stability) affect population survival and growth. We present a model in which the genotype is presented as a set folding free energies of cell proteins. Mutations occur upon replication, so stabilities of some proteins in daughter cells differ from those in the parent cell by amounts deduced from the distribution of mutational effects on protein stability. The genotype{\textendash}phenotype relationship posits that the cell{\textquoteright}s fitness (replication rate) is proportional to the concentration of its folded proteins and that unstable essential proteins result in lethality. Simulations reveal that lethal mutagenesis occurs at a mutation rate close to seven mutations in each replication of the genome for RNA viruses and at about half that rate for DNA-based organisms, in accord with earlier predictions from analytical theory and experimental results. This number appears somewhat dependent on the number of genes in the organisms and the organism{\textquoteright}s natural death rate. Further, our model reproduces the distribution of stabilities of natural proteins, in excellent agreement with experiments. We find that species with high mutation rates tend to have less stable proteins compared to species with low mutation rates.}, url = {http://dx.doi.org/10.1534/genetics.109.106492}, author = {Chen, Peiqiu and Eugene I. Shakhnovich} } @article {85521, title = {FOG: Fragment Optimized Growth Algorithm for the de Novo Generation of Molecules Occupying Druglike Chemical Space}, journal = {J. Chem. Inf. Model.}, volume = {49}, year = {2009}, pages = {1630 - 1642}, abstract = {An essential feature of all practical de novo molecule generating programs is the ability to focus the potential combinatorial explosion of grown molecules on a desired chemical space. It is a daunting task to balance the generation of new molecules with limitations on growth that produce desired features such as stability in water, synthetic accessibility, or druglikeness. We have developed an algorithm, Fragment Optimized Growth (FOG), which statistically biases the growth of molecules with desired features. At the heart of the algorithm is a Markov Chain which adds fragments to the nascent molecule in a biased manner, depending on the frequency of specific fragment-fragment connections in the database of chemicals it was trained on. We show that in addition to generating synthetically feasible molecules, it can be trained to grow new molecules that resemble desired classes of molecules such as drugs, natural products, and diversity-oriented synthetic products. In order to classify our grown molecules, we developed the Topology Classifier (TopClass) algorithm that is capable of classifying compounds, for example as drugs or nondrugs. The classification accuracies obtained with TopClass compare favorably with the literature. Furthermore, in contrast to {\textquotedblleft}black-box{\textquotedblright} approaches such as Neural Networks, TopClass brings to light characteristics of drugs that distinguish them from nondrugs.}, url = {http://dx.doi.org/10.1021/ci9000458}, author = {Peter S. Kutchukian and David Lou and Eugene I. Shakhnovich} } @article {101586, title = {Understanding ensemble protein folding at atomic detail}, journal = {J. Phys.: Condens. Matter}, volume = {20}, year = {2008}, pages = {283101}, abstract = {Although far from routine, simulating the folding of specific short protein chains on the computer, at a detailed atomic level, is starting to become a reality. This remarkable progress, which has been made over the last decade or so, allows a fundamental aspect of the protein folding process to be addressed, namely its statistical nature. In order to make quantitative comparisons with experimental kinetic data a complete ensemble view of folding must be achieved, with key observables averaged over the large number of microscopically different folding trajectories available to a protein chain. Here we review recent advances in atomic-level protein folding simulations and the new insight provided by them into the protein folding process. An important element in understanding ensemble folding kinetics are methods for analyzing many separate folding trajectories, and we discuss techniques developed to condense the large amount of information contained in an ensemble of trajectories into a manageable picture of the folding process.}, url = {http://dx.doi.org/10.1088/0953-8984/20/28/283101}, author = {Wallin, Stefan and Eugene I. Shakhnovich} } @article {101501, title = {Sensitivity-dependent model of protein{\textendash}protein interaction networks}, journal = {Phys. Biol.}, volume = {5}, year = {2008}, month = {2008}, abstract = {The scale free structure p(k) ~ k-γ of protein{\textendash}protein interaction networks can be reproduced by a static physical model in simulation. We inspect the model theoretically, and find the key reason for the model generating apparent scale free degree distributions. This explanation provides a generic mechanism of {\textquoteright}scale free{\textquoteright} networks. Moreover, we predict the dependence of γ on experimental protein concentrations or other sensitivity factors in detecting interactions, and find experimental evidence to support the prediction.}, url = {http://dx.doi.org/10.1088/1478-3975/5/3/036011}, author = {Jingshan Zhang and Eugene I. Shakhnovich} } @article {101386, title = {Improvisation in evolution of genes and genomes: whose structure is it anyway?}, journal = {Current Opinion in Structural Biology}, volume = {18}, year = {2008}, month = {2008}, pages = {375 - 381}, abstract = {Significant progress has been made in recent years in a variety of seemingly unrelated fields such as sequencing, protein structure prediction, and high-throughput transcriptomics and metabolomics. At the same time, new microscopic models have been developed that made it possible to analyze the evolution of genes and genomes from first principles. The results from these efforts enable, for the first time, a comprehensive insight into the evolution of complex systems and organisms on all scales {\textemdash} from sequences to organisms and populations. Every newly sequenced genome uncovers new genes, families, and folds. Where do these new genes come from? How do gene duplication and subsequent divergence of sequence and structure affect the fitness of the organism? What role does regulation play in the evolution of proteins and folds? Emerging synergism between data and modeling provides first robust answers to these questions.}, url = {http://dx.doi.org/10.1016/j.sbi.2008.02.007}, author = {Shakhnovich, Boris E. and Eugene I. Shakhnovich} } @article {101331, title = {Native atomic burials, supplemented by physically motivated hydrogen bond constraints, contain sufficient information to determine the tertiary structure of small globular proteins}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {70}, year = {2008}, month = {2008}, pages = {971 - 983}, abstract = {We investigate the possibility that atomic burials, as measured by their distances from the structural geometrical center, contain sufficient information to determine the tertiary structure of globular proteins. We report Monte Carlo simulated annealing results of all-atom hard-sphere models in continuous space for four small proteins: the all-β WW-domain 1E0L, the α/β protein-G 1IGD, the all-α engrailed homeo-domain 1ENH, and the α + β engineered monomeric form of the Cro protein 1ORC. We used as energy function the sum over all atoms, labeled by i, of |Ri - R i*|, where Ri is the atomic distance from the center of coordinates, or central distance, and R i* is the "ideal{\textquotedblright} central distance obtained from the native structure. Hydrogen bonds were taken into consideration by the assignment of two ideal distances for backbone atoms forming hydrogen bonds in the native structure depending on the formation of a geometrically defined bond, independently of bond partner. Lowest energy final conformations turned out to be very similar to the native structure for the four proteins under investigation and a strong correlation was observed between energy and distance root mean square deviation (DRMS) from the native in the case of all-β 1E0L and α/β 1IGD. For all α 1ENH and α + β 1ORC the overall correlation between energy and DRMS among final conformations was not as high because some trajectories resulted in high DRMS but low energy final conformations in which α-helices adopted a non-native mutual orientation. Comparison between central distances and actual accessible surface areas corroborated the implicit assumption of correlation between these two quantities. The Z-score obtained with this native-centric potential in the discrimination of native 1ORC from a set of random compact structures confirmed that it contains a much smaller amount of native information when compared to a traditional contact Go potential but indicated that simple sequence-dependent burial potentials still need some improvement in order to attain a similar discriminability. Taken together, our results suggest that central distances, in conjunction to physically motivated hydrogen bond constraints, contain sufficient information to determine the native conformation of these small proteins and that a solution to the folding problem for globular proteins could arise from sufficiently accurate burial predictions from sequence followed by minimization of a burial-dependent energy function. Proteins 2008. {\textcopyright} 2007 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/prot.21571}, author = {Pereira de Ara{\'u}jo, Ant{\^o}nio F. and Gomes, Antonio L. C. and Bursztyn, Alexandre A. and Eugene I. Shakhnovich} } @article {101291, title = {Constraints imposed by non-functional protein{\textendash}protein interactions on gene expression and proteome size}, journal = {Mol Syst Biol}, volume = {4}, year = {2008}, month = {2008}, pages = {210}, abstract = {Crowded intracellular environments present a challenge for proteins to form functional specific complexes while reducing non-functional interactions with promiscuous non-functional partners. Here we show how the need to minimize the waste of resources to non-functional interactions limits the proteome diversity and the average concentration of co-expressed and co-localized proteins. Using the results of high-throughput Yeast 2-Hybrid experiments, we estimate the characteristic strength of non-functional protein{\textendash}protein interactions. By combining these data with the strengths of specific interactions, we assess the fraction of time proteins spend tied up in non-functional interactions as a function of their overall concentration. This allows us to sketch the phase diagram for baker{\textquoteright}s yeast cells using the experimentally measured concentrations and subcellular localization of their proteins. The positions of yeast compartments on the phase diagram are consistent with our hypothesis that the yeast proteome has evolved to operate closely to the upper limit of its size, whereas keeping individual protein concentrations sufficiently low to reduce non-functional interactions. These findings have implication for conceptual understanding of intracellular compartmentalization, multicellularity and differentiation.}, url = {http://dx.doi.org/10.1038/msb.2008.48}, author = {Jingshan Zhang and Sergei Maslov and Eugene I. Shakhnovich} } @article {101176, title = {Statistically enhanced promiscuity of structurally correlated patterns}, journal = {Phys. Rev. E}, volume = {77}, year = {2008}, month = {2008}, pages = {020901}, abstract = {We predict that patterns with correlated surface density of atoms have statistically higher promiscuity (ability to bind stronger to an arbitrary pattern) as compared with noncorrelated patterns with the same average surface density. We suggest that this constitutes a generic design principle for highly connected proteins (hubs) in protein interaction networks. We develop an analytical theory for this effect. We show that our key predictions are generic and independent, qualitatively, on the specific form of the interatomic interaction potential, provided it has a finite range.}, url = {http://dx.doi.org/10.1103/PhysRevE.77.020901}, author = {Lukatsky, D. B. and Shakhnovich, E.I.} } @article {101086, title = {Understanding Protein Evolution: From Protein Physics to Darwinian Selection}, journal = {Annual Review of Physical Chemistry}, volume = {59}, year = {2008}, month = {2008}, pages = {105 - 127}, abstract = {Efforts in whole-genome sequencing and structural proteomics start to provide a global view of the protein universe, the set of existing protein structures and sequences. However, approaches based on the selection of individual sequences have not been entirely successful at the quantitative description of the distribution of structures and sequences in the protein universe because evolutionary pressure acts on the entire organism, rather than on a particular molecule. In parallel to this line of study, studies in population genetics and phenomenological molecular evolution established a mathematical framework to describe the changes in genome sequences in populations of organisms over time. Here, we review both microscopic (physics-based) and macroscopic (organism-level) models of protein-sequence evolution and demonstrate that bridging the two scales provides the most complete description of the protein universe starting from clearly defined, testable, and physiologically relevant assumptions.}, url = {http://dx.doi.org/10.1146/annurev.physchem.58.032806.104449}, author = {Zeldovich, Konstantin B. and Eugene I. Shakhnovich} } @article {101026, title = {Identifying critical residues in protein folding: Insights from ϕ-value and Pfold analysis}, journal = {J. Chem. Phys.}, volume = {129}, year = {2008}, month = {2008}, pages = {095108}, abstract = {We apply a simulational proxy of the ϕ-value analysis and perform extensive mutagenesis experiments to identify the nucleating residues in the folding "reactions{\textquotedblright} of two small lattice Go polymers with different native geometries. Our findings show that for the more complex native fold (i.e., the one that is rich in nonlocal, long-range bonds), mutation of the residues that form the folding nucleus leads to a considerably larger increase in the folding time than the corresponding mutations in the geometry that is predominantly local. These results are compared to data obtained from an accurate analysis based on the reaction coordinate folding probability Pfold and on structural clustering methods. Our study reveals a complex picture of the transition state ensemble. For both protein models, the transition state ensemble is rather heterogeneous and splits up into structurally different populations. For the more complex geometry the identified subpopulations are actually structurally disjoint. For the less complex native geometry we found a broad transition state with microscopic heterogeneity. These findings suggest that the existence of multiple transition state structures may be linked to the geometric complexity of the native fold. For both geometries, the identification of the folding nucleus via the Pfold analysis agrees with the identification of the folding nucleus carried out with the ϕ-value analysis. For the most complex geometry, however, the applied methodologies give more consistent results than for the more local geometry. The study of the transition state structure reveals that the nucleus residues are not necessarily fully native in the transition state. Indeed, it is only for the more complex geometry that two of the five critical residues show a considerably high probability of having all its native bonds formed in the transition state. Therefore, one concludes that, in general, the ϕ-value correlates with the acceleration/deceleration of folding induced by mutation, rather than with the degree of nativeness of the transition state, and that the "traditional{\textquotedblright} interpretation of ϕ-values may provide a more realistic picture of the structure of the transition state only for more complex native geometries.}, url = {http://dx.doi.org/doi:10.1063/1.2973624}, author = {Fa{\'\i}sca, P.F.N. and Travasso, R.D.M. and Ball, R. C. and Shakhnovich, E.I.} } @article {100981, title = {Universality and diversity of folding mechanics for three-helix bundle proteins}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {105}, year = {2008}, month = {2008}, pages = {895 - 900}, abstract = {In this study we evaluate, at full atomic detail, the folding processes of two small helical proteins, the B domain of protein A and the Villin headpiece. Folding kinetics are studied by performing a large number of ab initio Monte Carlo folding simulations using a single transferable all-atom potential. Using these trajectories, we examine the relaxation behavior, secondary structure formation, and transition-state ensembles (TSEs) of the two proteins and compare our results with experimental data and previous computational studies. To obtain a detailed structural information on the folding dynamics viewed as an ensemble process, we perform a clustering analysis procedure based on graph theory. Moreover, rigorous p fold analysis is used to obtain representative samples of the TSEs and a good quantitative agreement between experimental and simulated Φ values is obtained for protein A. Φ values for Villin also are obtained and left as predictions to be tested by future experiments. Our analysis shows that the two-helix hairpin is a common partially stable structural motif that gets formed before entering the TSE in the studied proteins. These results together with our earlier study of Engrailed Homeodomain and recent experimental studies provide a comprehensive, atomic-level picture of folding mechanics of three-helix bundle proteins.}, url = {http://dx.doi.org/10.1073/pnas.0707284105}, author = {Yang, Jae Shick and Wallin, Stefan and Eugene I. Shakhnovich} } @article {101471, title = {Energetics of protein{\textendash}DNA interactions}, journal = {Nucl. Acids Res.}, volume = {35}, year = {2007}, month = {2007}, pages = {1039 - 1047}, abstract = {Protein{\textendash}DNA interactions are vital for many processes in living cells, especially transcriptional regulation and DNA modification. To further our understanding of these important processes on the microscopic level, it is necessary that theoretical models describe the macromolecular interaction energetics accurately. While several methods have been proposed, there has not been a careful comparison of how well the different methods are able to predict biologically important quantities such as the correct DNA binding sequence, total binding free energy and free energy changes caused by DNA mutation. In addition to carrying out the comparison, we present two important theoretical models developed initially in protein folding that have not yet been tried on protein{\textendash}DNA interactions. In the process, we find that the results of these knowledge-based potentials show a strong dependence on the interaction distance and the derivation method. Finally, we present a knowledge-based potential that gives comparable or superior results to the best of the other methods, including the molecular mechanics force field AMBER99.}, url = {http://dx.doi.org/10.1093/nar/gkl1103}, author = {Donald, Jason E. and Chen, William W. and Eugene I. Shakhnovich} } @article {101421, title = {Quantifying fitness distributions and phenotypic relationships in recombinant yeast populations}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {104}, year = {2007}, month = {2007}, pages = {10553 - 10558}, abstract = {Studies of the role of sex in evolution typically involve a longitudinal comparison of a single ancestor to several intermediate descendants and to one terminally evolved descendant after many generations of adaptation under a given selective regime. Here we take a complementary, statistical approach to sex in evolution, by describing the distribution of phenotypic similarity in a population of yeast F1 meiotic recombinants. By applying graph theory to fitness measurements of thousands of Saccharomyces cerevisiae recombinants treated with 10 mechanistically distinct, growth-inhibitory small-molecule perturbagens (SMPs), we show that the network of phenotypic similarity among F1 recombinants exhibits a scale-free degree distribution. F1 recombinants are often phenotypically unique and sometimes exceptional, and their fitness strengths are unevenly distributed across the 10 compound treatments. By contrast, highly phenotypically similar F1 recombinants constitute failing hubs that display below-average fitness across all compound treatments and are candidate substrates for purifying selection. Comparison of the F1 generation with the parental strains reveals that (i) there is a specialist more fit in any given single condition than any of the parents but (ii) only rarely are there generalists that exhibit greater fitness than both parental strains across a majority of conditions. This analysis allows us to evaluate and to gain better theoretical understanding of the costs and benefits of sex in the F1 generation.}, url = {http://dx.doi.org/10.1073/pnas.0704037104}, author = {Perlstein, Ethan O. and Deeds, Eric J. and Ashenberg, Orr and Eugene I. Shakhnovich and Schreiber, Stuart L.} } @article {101396, title = {A First-Principles Model of Early Evolution: Emergence of Gene Families, Species, and Preferred Protein Folds}, journal = {PLoS Comput Biol}, volume = {3}, year = {2007}, month = {2007}, abstract = {Author Summary

Here, we address the question of how Darwinian evolution of organisms determines molecular evolution of their proteins and genomes. We developed a microscopic ab initio model of early biological evolution where the fitness (essentially lifetime) of an organism is explicitly related to the evolving sequences of its proteins. The main assumption of the model is that the death rate of an organism is determined by the stability of the least stable of their proteins. A lattice model is used to calculate stability of all proteins in a genome from their amino acid sequence. The simulation of the model starts from 100 identical organisms, each carrying the same random gene, and proceeds via random mutations, gene duplication, organism births via replication, and organism deaths. We find that exponential population growth is possible only after the discovery of a very small number of specific advantageous protein structures. The number of genes in the evolving organisms depends on the mutation rate, demonstrating the intricate relationship between the genome sizes and protein stability requirements. Further, the model explains the observed power-law distributions of protein family and superfamily sizes, as well as the scale-free character of protein structural similarity graphs. Together, these results and their analysis suggest a plausible comprehensive scenario of emergence of the protein universe in early biological evolution.

}, url = {http://dx.doi.org/10.1371/journal.pcbi.0030139}, author = {Zeldovich, Konstantin B. and Chen, Peiqiu and Shakhnovich, Boris E. and Eugene I. Shakhnovich} } @article {101321, title = {Structural Similarity Enhances Interaction Propensity of Proteins}, journal = {Journal of Molecular Biology}, volume = {365}, year = {2007}, month = {2007}, pages = {1596 - 1606}, abstract = {We study statistical properties of interacting protein-like surfaces and predict two strong, related effects: (i) statistically enhanced self-attraction of proteins; (ii) statistically enhanced attraction of proteins with similar structures. The effects originate in the fact that the probability to find a pattern self-match between two identical, even randomly organized interacting protein surfaces is always higher compared with the probability for a pattern match between two different, promiscuous protein surfaces. This theoretical finding explains statistical prevalence of homodimers in protein{\textendash}protein interaction networks reported earlier. Further, our findings are confirmed by the analysis of curated database of protein complexes that showed highly statistically significant overrepresentation of dimers formed by structurally similar proteins with highly divergent sequences ("superfamily heterodimers{\textquotedblright}). We suggest that promiscuous homodimeric interactions pose strong competitive interactions for heterodimers evolved from homodimers. Such evolutionary bottleneck is overcome using the negative design evolutionary pressure applied against promiscuous homodimer formation. This is achieved through the formation of highly specific contacts formed by charged residues as demonstrated both in model and real superfamily heterodimers.}, url = {http://dx.doi.org/10.1016/j.jmb.2006.11.020}, author = {Lukatsky, D. B. and Shakhnovich, B.E. and Mintseris, J. and Shakhnovich, E.I.} } @article {101311, title = {Protein stability imposes limits on organism complexity and speed of molecular evolution}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {104}, year = {2007}, month = {2007}, pages = {16152 - 16157}, abstract = {Classical population genetics a priori assigns fitness to alleles without considering molecular or functional properties of proteins that these alleles encode. Here we study population dynamics in a model where fitness can be inferred from physical properties of proteins under a physiological assumption that loss of stability of any protein encoded by an essential gene confers a lethal phenotype. Accumulation of mutations in organisms containing Γ genes can then be represented as diffusion within the Γ-dimensional hypercube with adsorbing boundaries determined, in each dimension, by loss of a protein{\textquoteright}s stability and, at higher stability, by lack of protein sequences. Solving the diffusion equation whose parameters are derived from the data on point mutations in proteins, we determine a universal distribution of protein stabilities, in agreement with existing data. The theory provides a fundamental relation between mutation rate, maximal genome size, and thermodynamic response of proteins to point mutations. It establishes a universal speed limit on rate of molecular evolution by predicting that populations go extinct (via lethal mutagenesis) when mutation rate exceeds approximately six mutations per essential part of genome per replication for mesophilic organisms and one to two mutations per genome per replication for thermophilic ones. Several RNA viruses function close to the evolutionary speed limit, whereas error correction mechanisms used by DNA viruses and nonmutant strains of bacteria featuring various genome lengths and mutation rates have brought these organisms universally ≈1,000-fold below the natural speed limit.}, url = {http://dx.doi.org/10.1073/pnas.0705366104}, author = {Zeldovich, Konstantin B. and Chen, Peiqiu and Eugene I. Shakhnovich} } @article {101216, title = {Positive and Negative Design in Stability and Thermal Adaptation of Natural Proteins}, journal = {PLoS Comput Biol}, volume = {3}, year = {2007}, month = {2007}, abstract = {Author Summary

What mechanisms does Nature use in her quest for thermophilic proteins? It is known that stability of a protein is mainly determined by the energy gap, or the difference in energy, between native state and a set of incorrectly folded (misfolded) conformations. Here we show that Nature makes thermophilic proteins by widening this gap from both ends. The energy of the native state of a protein is decreased by selecting strongly attractive amino acids at positions that are in contact in the native state (positive design). Simultaneously, energies of the misfolded conformations are increased by selection of strongly repulsive amino acids at positions that are distant in native structure; however, these amino acids will interact repulsively in the misfolded conformations (negative design). These fundamental principles of protein design are manifested in the "from both ends of the hydrophobicity scale{\textquotedblright} trend observed in thermophilic adaptation, whereby proteomes of thermophilic proteins are enriched in extreme amino acids{\textemdash}hydrophobic and charged{\textemdash}at the expense of polar ones. Hydrophobic amino acids contribute mostly to the positive design, while charged amino acids that repel each other in non-native conformations of proteins contribute to negative design. Our results provide guidance in rational design of proteins with selected thermal properties.

}, url = {http://dx.doi.org/10.1371/journal.pcbi.0030052}, author = {Berezovsky, Igor N and Zeldovich, Konstantin B. and Eugene I. Shakhnovich} } @article {101261, title = {Robust protein{\textendash}protein interactions in crowded cellular environments}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {104}, year = {2007}, month = {2007}, pages = {14952 - 14957}, abstract = {The capacity of proteins to interact specifically with one another underlies our conceptual understanding of how living systems function. Systems-level study of specificity in protein{\textendash}protein interactions is complicated by the fact that the cellular environment is crowded and heterogeneous; interaction pairs may exist at low relative concentrations and thus be presented with many more opportunities for promiscuous interactions compared with specific interaction possibilities. Here we address these questions by using a simple computational model that includes specifically designed interacting model proteins immersed in a mixture containing hundreds of different unrelated ones; all of them undergo simulated diffusion and interaction. We find that specific complexes are quite robust to interference from promiscuous interaction partners only in the range of temperatures T design > T > T rand. At T > T design, specific complexes become unstable, whereas at T < T rand, formation of specific complexes is suppressed by promiscuous interactions. Specific interactions can form only if T design > T rand. This condition requires an energy gap between binding energy in a specific complex and set of binding energies between randomly associating proteins, providing a general physical constraint on evolutionary selection or design of specific interacting protein interfaces. This work has implications for our understanding of how the protein repertoire functions and evolves within the context of cellular systems.}, url = {http://dx.doi.org/10.1073/pnas.0702766104}, author = {Deeds, Eric J. and Ashenberg, Orr and Gerardin, Jaline and Eugene I. Shakhnovich} } @article {101076, title = {All-Atom Ab Initio Folding of a Diverse Set of Proteins}, journal = {Structure}, volume = {15}, year = {2007}, month = {2007}, pages = {53 - 63}, abstract = {Natural proteins fold to a unique, thermodynamically dominant state. Modeling of the folding process and prediction of the native fold of proteins are two major unsolved problems in biophysics. Here, we show successful all-atom ab~initio folding of a representative diverse set of proteins by using a minimalist transferable-energy model that consists of two-body atom-atom interactions, hydrogen bonding, and~a local sequence-energy term that models sequence-specific chain stiffness. Starting from a random coil, the native-like structure was observed during replica exchange Monte Carlo (REMC) simulation for most proteins regardless of their structural classes; the lowest energy structure was close to native{\textemdash}in the range of 2{\textendash}6 {\r A} root-mean-square deviation (rmsd). Our results demonstrate that the successful folding of a protein chain to its native state is governed by only a few crucial energetic terms.}, url = {http://dx.doi.org/10.1016/j.str.2006.11.010}, author = {Yang, Jae Shick and Chen, William W. and Skolnick, Jeffrey and Eugene I. Shakhnovich} } @article {101051, title = {Description of atomic burials in compact globular proteins by Fermi-Dirac probability distributions}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {66}, year = {2007}, month = {2007}, pages = {304 - 320}, abstract = {We perform a statistical analysis of atomic distributions as a function of the distance R from the molecular geometrical center in a nonredundant set of compact globular proteins. The number of atoms increases quadratically for small R, indicating a constant average density inside the core, reaches a maximum at a size-dependent distance Rmax, and falls rapidly for larger R. The empirical curves turn out to be consistent with the volume increase of spherical concentric solid shells and a Fermi-Dirac distribution in which the distance R plays the role of an effective atomic energy ϵ(R) = R. The effective chemical potential μ governing the distribution increases with the number of residues, reflecting the size of the protein globule, while the temperature parameter β decreases. Interestingly, βμ is not as strongly dependent on protein size and appears to be tuned to maintain approximately half of the atoms in the high density interior and the other half in the exterior region of rapidly decreasing density. A normalized size-independent distribution was obtained for the atomic probability as a function of the reduced distance, r = R/Rg, where Rg is the radius of gyration. The global normalized Fermi distribution, F(r), can be reasonably decomposed in Fermi-like subdistributions for different atomic types τ, Fτ(r), with ΣτFτ(r) = F(r), which depend on two additional parameters μτ and hτ. The chemical potential μτ affects a scaling prefactor and depends on the overall frequency of the corresponding atomic type, while the maximum position of the subdistribution is determined by hτ, which appears in a type-dependent atomic effective energy, ετ(r) = hτr, and is strongly correlated to available hydrophobicity scales. Better adjustments are obtained when the effective energy is not assumed to be necessarily linear, or ετ*(r) = hτ*rα,, in which case a correlation with hydrophobicity scales is found for the product ατhτ*. These results indicate that compact globular proteins are consistent with a thermodynamic system governed by hydrophobic-like energy functions, with reduced distances from the geometrical center, reflecting atomic burials, and provide a conceptual framework for the eventual prediction from sequence of a few parameters from which whole atomic probability distributions and potentials of mean force can be reconstructed. Proteins 2007. {\textcopyright} 2006 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/prot.21137}, author = {Gomes, Antonio L. C. and de Rezende, J{\'u}lia R. and de Ara{\'u}jo, Ant{\^o}nio F. Pereira and Eugene I. Shakhnovich} } @article {101046, title = {The Folding Mechanics of a Knotted Protein}, journal = {Journal of Molecular Biology}, volume = {368}, year = {2007}, month = {2007}, pages = {884 - 893}, abstract = {An increasing number of proteins are being discovered with a remarkable and somewhat surprising feature, a knot in their native structures. How the polypeptide chain is able to "knot{\textquotedblright} itself during the folding process to form these highly intricate protein topologies is not known. Here we perform a computational study on the 160-amino-acid homodimeric protein YibK, which, like other proteins in the SpoU family of MTases, contains a deep trefoil knot in its C-terminal region. In this study, we use a coarse-grained Cα-chain representation and Langevin dynamics to study folding kinetics. We find that specific, attractive nonnative interactions are critical for knot formation. In the absence of these interactions, i.e., in an energetics driven entirely by native interactions, knot formation is exceedingly unlikely. Further, we find, in concert with recent experimental data on YibK, two parallel folding pathways that we attribute to an early and a late formation of the trefoil knot, respectively. For both pathways, knot formation occurs before dimerization. A bioinformatics analysis of the SpoU family of proteins reveals further that the critical nonnative interactions may originate from evolutionary conserved hydrophobic segments around the knotted region.}, url = {http://dx.doi.org/10.1016/j.jmb.2007.02.035}, author = {Wallin, Stefan and Zeldovich, Konstantin B. and Eugene I. Shakhnovich} } @article {101021, title = {Protein and DNA Sequence Determinants of Thermophilic Adaptation}, journal = {PLoS Comput Biol}, volume = {3}, year = {2007}, month = {2007}, abstract = {There have been considerable attempts in the past to relate phenotypic trait{\textemdash}habitat temperature of organisms{\textemdash}to their genotypes, most importantly compositions of their genomes and proteomes. However, despite accumulation of anecdotal evidence, an exact and conclusive relationship between the former and the latter has been elusive. We present an exhaustive study of the relationship between amino acid composition of proteomes, nucleotide composition of DNA, and optimal growth temperature (OGT) of prokaryotes. Based on 204 complete proteomes of archaea and bacteria spanning the temperature range from -10 {\textdegree}C to 110 {\textdegree}C, we performed an exhaustive enumeration of all possible sets of amino acids and found a set of amino acids whose total fraction in a proteome is correlated, to a remarkable extent, with the OGT. The universal set is Ile, Val, Tyr, Trp, Arg, Glu, Leu (IVYWREL), and the correlation coefficient is as high as 0.93. We also found that the G + C content in 204 complete genomes does not exhibit a significant correlation with OGT (R = -0.10). On the other hand, the fraction of A + G in coding DNA is correlated with temperature, to a considerable extent, due to codon patterns of IVYWREL amino acids. Further, we found strong and independent correlation between OGT and the frequency with which pairs of A and G nucleotides appear as nearest neighbors in genome sequences. This adaptation is achieved via codon bias. These findings present a direct link between principles of proteins structure and stability and evolutionary mechanisms of thermophylic adaptation. On the nucleotide level, the analysis provides an example of how nature utilizes codon bias for evolutionary adaptation to extreme conditions. Together these results provide a complete picture of how compositions of proteomes and genomes in prokaryotes adjust to the extreme conditions of the environment.}, url = {http://dx.doi.org/10.1371/journal.pcbi.0030005}, author = {Zeldovich, Konstantin B. and Berezovsky, Igor N and Eugene I. Shakhnovich} } @article {101001, title = {A knowledge-based move set for protein folding}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {66}, year = {2007}, month = {2007}, pages = {682 - 688}, abstract = {The free energy landscape of protein folding is rugged, occasionally characterized by compact, intermediate states of low free energy. In computational folding, this landscape leads to trapped, compact states with incorrect secondary structure. We devised a residue-specific, protein backbone move set for efficient sampling of protein-like conformations in computational folding simulations. The move set is based on the selection of a small set of backbone dihedral angles, derived from clustering dihedral angles sampled from experimental structures. We show in both simulated annealing and replica exchange Monte Carlo (REMC) simulations that the knowledge-based move set, when compared with a conventional move set, shows statistically significant improved ability at overcoming kinetic barriers, reaching deeper energy minima, and achieving correspondingly lower RMSDs to native structures. The new move set is also more efficient, being able to reach low energy states considerably faster. Use of this move set in determining the energy minimum state and for calculating thermodynamic quantities is discussed. Proteins 2007. {\textcopyright} 2006 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/prot.21237}, author = {Chen, William W. and Yang, Jae Shick and Eugene I. Shakhnovich} } @article {101581, title = {Physics and evolution of protein-protein interactions}, journal = {The FASEB Journal}, volume = {20}, year = {2006}, pages = {A1473}, abstract = {We will discuss recent developments in bioinformatics analysis and theoretical studies of protein-protein interactions in living cells {\textendash} an important aspect of systems biology. First, we consider the network of protein-protein interactions and demonstrate that two published independent measurements of these interactions produce graphs that are only weakly correlated with one another despite their strikingly similar topology. We then propose a physical model based on the fundamental principle that (de)solvation is a major physical factor in protein-protein interactions. This model reproduces not only the scale-free nature of such graphs but also a number of higher-order correlations in these networks. A key support of the model is provided by the discovery of a significant correlation between number of interactions made by a protein and the fraction of hydrophobic residues on its surface. Next, we discuss a number of fundamental models for specific protein-protein interactions that provide deep insight into how specific protein multimers could have evolved. In particular we show that homodimers are most likely to have been precursors of modern protein complexes (homodimers are prevalent in modern cells {\textendash} phenomenon of "molecular narcissism"). Subsequent evolution created homodimers in the process of {\textquoteleft}{\textquoteright}negative design{\textquoteright}{\textquoteright} against non-specific and homodimeric associations.}, url = {http://www.fasebj.org/cgi/content/meeting_abstract/20/5/A1473-c}, author = {Eugene I. Shakhnovich} } @article {101541, title = {Genetic instability and the quasispecies model}, journal = {Journal of Theoretical Biology}, volume = {241}, year = {2006}, month = {2006}, pages = {216 - 222}, abstract = {Genetic instability is a defining characteristic of cancers. Microsatellite instability (MIN) leads to by elevated point mutation rates, whereas chromosomal instability (CIN) refers to increased rates of losing or gaining whole chromosomes or parts of chromosomes during cell division. CIN and MIN are, in general, mutually exclusive. The quasispecies model is a very successful theoretical framework for the study of evolution at high mutation rates. It predicts the existence of an experimentally verified error catastrophe. This catastrophe occurs when the mutation rates exceed a threshold value, the error threshold, above which replicative infidelity is incompatible with cell survival. We analyse the semiconservative quasispecies model of both MIN and CIN tumors. We consider the role of post-methylation DNA repair in tumor cells and demonstrate that DNA repair is fundamental to the nature of the error catastrophe in both types of tumors. We find that CIN introduces a plateau in the maximum viable mutation rate for a repair-free model, which does not exist in the case of MIN. This provides a plausible explanation for the mutual exclusivity of CIN and MIN.}, url = {http://dx.doi.org/10.1016/j.jtbi.2005.11.018}, author = {Brumer, Yisroel and Michor, Franziska and Eugene I. Shakhnovich} } @article {101516, title = {PDB-UF: database of predicted enzymatic functions for unannotated protein structures from structural genomics}, journal = {BMC Bioinformatics}, volume = {7}, year = {2006}, month = {2006}, abstract = {The number of protein structures from structural genomics centers dramatically increases in the Protein Data Bank (PDB). Many of these structures are functionally unannotated because they have no sequence similarity to proteins of known function. However, it is possible to successfully infer function using only structural similarity.}, url = {http://dx.doi.org/10.1186/1471-2105-7-53}, author = {Grotthuss, Marcin von and Plewczynski, Dariusz and Ginalski, Krzysztof and Rychlewski, Leszek and Eugene I. Shakhnovich} } @article {101526, title = {Semiconservative quasispecies equations for polysomic genomes: The haploid case}, journal = {Journal of Theoretical Biology}, volume = {241}, year = {2006}, month = {2006}, pages = {791 - 805}, abstract = {This paper develops the semiconservative quasispecies equations for genomes consisting of an arbitrary number of chromosomes. We assume that the chromosomes are distinguishable, so that we are effectively considering haploid genomes. We derive the quasispecies equations under the assumption of arbitrary lesion repair efficiency, and consider the cases of both random and immortal strand chromosome segregation. We solve the model in the limit of infinite sequence length for the case of the static single fitness peak landscape, where the master genome has a first-order growth rate constant of k \> 1 , and all other genomes have a first-order growth rate constant of 1. If we assume that each chromosome can tolerate an arbitrary number of lesions, so that only one master copy of the strands is necessary for a functional chromosome, then for random chromosome segregation we obtain an equilibrium mean fitness of κ {\textasciimacron} ( t = $\infty$ ) = k 2 e - ( 1 / N ) μ λ / 2 + e - ( 1 / N ) μ ( 1 - λ / 2 ) 2 N - 1 , below the error catastrophe, while for immortal strand co-segregation we obtain κ {\textasciimacron} ( t = $\infty$ ) = k [ e - μ ( 1 - λ / 2 ) + e - μ λ / 2 - 1 ] (N denotes the number of chromosomes, λ denotes the lesion repair efficiency, and μ = ε L , where ε is the per base-pair mismatch probability, and L is the total genome length). It follows that immortal strand co-segregation leads to significantly better preservation of the master genome than random segregation when lesion repair is imperfect. Based on this result, we conjecture that certain classes of tumor cells exhibit immortal strand co-segregation.}, url = {http://dx.doi.org/10.1016/j.jtbi.2006.01.016}, author = {Tannenbaum, Emmanuel and Sherley, James L. and Eugene I. Shakhnovich} } @article {101491, title = {Physical Origins of Protein Superfamilies}, journal = {Journal of Molecular Biology}, volume = {357}, year = {2006}, month = {2006}, pages = {1335 - 1343}, abstract = {In this work, we discovered a fundamental connection between selection for protein stability and emergence of preferred structures of proteins. Using a standard exact three-dimensional lattice model we evolve sequences starting from random ones and determine the exact native structure after each mutation. Acceptance of mutations is biased to select for stable proteins. We found that certain structures, "wonderfolds{\textquotedblright}, are independently discovered numerous times as native states of stable proteins in many unrelated runs of selection. The strong dependence of lattice fold usage on the structural determinant of designability quantitatively reproduces uneven fold usage in natural proteins. Diversity of sequences that fold into wonderfold structures gives rise to superfamilies, i.e. sets of dissimilar sequences that fold into the same or very similar structures. The present work establishes a model of pre-biotic structure selection, which identifies dominant structural patterns emerging upon optimization of proteins for survival in a hot environment. Convergently discovered pre-biotic initial superfamilies with wonderfold structures could have served as a seed for subsequent biological evolution involving gene duplications and divergence.}, url = {http://dx.doi.org/10.1016/j.jmb.2006.01.081}, author = {Zeldovich, Konstantin B. and Berezovsky, Igor N and Eugene I. Shakhnovich} } @article {101316, title = {Common Motifs and Topological Effects in the Protein Folding Transition State}, journal = {Journal of Molecular Biology}, volume = {359}, year = {2006}, month = {2006}, pages = {1075 - 1085}, abstract = {Through extensive experiment, simulation, and analysis of protein S6 (1RIS), we find that variations in nucleation and folding pathway between circular permutations are determined principally by the restraints of topology and specific nucleation, and affected by changes in chain entropy. Simulations also relate topological features to experimentally measured stabilities. Despite many sizable changes in ϕ values and the structure of the transition state ensemble that result from permutation, we observe a common theme: the critical nucleus in each of the mutants share a subset of residues that can be mapped to the critical nucleus residues of the wild-type. Circular permutations create new N and C termini, which are the location of the largest disruption of the folding nucleus, leading to a decrease in both ϕ values and the role in nucleation. Mutant nuclei are built around the wild-type nucleus but are biased towards different parts of the S6 structure depending on the topological and entropic changes induced by the location of the new N and C termini.}, url = {http://dx.doi.org/10.1016/j.jmb.2006.04.015}, author = {Hubner, Isaac A. and Lindberg, Magnus and Haglund, Ellinor and Oliveberg, Mikael and Eugene I. Shakhnovich} } @article {101326, title = {On the origin and highly likely completeness of single-domain protein structures}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {103}, year = {2006}, month = {2006}, pages = {2605 - 2610}, abstract = {The size and origin of the protein fold universe is of fundamental and practical importance. Analyzing randomly generated, compact sticky homopolypeptide conformations constructed in generic simplified and all-atom protein models, all have similar folds in the library of solved structures, the Protein Data Bank, and conversely, all compact, single-domain protein structures in the Protein Data Bank have structural analogues in the compact model set. Thus, both sets are highly likely complete, with the protein fold universe arising from compact conformations of hydrogen-bonded, secondary structures. Because side chains are represented by their Cβ atoms, these results also suggest that the observed protein folds are insensitive to the details of side-chain packing. Sequence specificity enters both in fine-tuning the structure and thermodynamically stabilizing a given fold with respect to the set of alternatives. Scanning the models against a three-dimensional active-site library, close geometric matches are frequently found. Thus, the presence of active-site-like geometries also seems to be a consequence of the packing of compact, secondary structural elements. These results have significant implications for the evolution of protein structure and function.}, url = {http://dx.doi.org/10.1073/pnas.0509379103}, author = {Zhang, Yang and Hubner, Isaac A. and Arakaki, Adrian K. and Eugene Shakhnovich and Skolnick, Jeffrey} } @article {101281, title = {Identification of the minimal protein-folding nucleus through loop-entropy perturbations}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {103}, year = {2006}, month = {2006}, pages = {4083 - 4088}, abstract = {To explore the plasticity and structural constraints of the protein-folding nucleus we have constructed through circular permutation four topological variants of the ribosomal protein S6. In effect, these topological variants represent entropy mutants with maintained spatial contacts. The proteins were characterized at two complementary levels of detail: by φ-value analysis estimating the extent of contact formation in the transition-state ensemble and by Hammond analysis measuring the site-specific growth of the folding nucleus. The results show that, although the loop-entropy alterations markedly influence the appearance and structural location of the folding nucleus, it retains a common motif of one helix docking against two strands. This nucleation motif is built around a shared subset of side chains in the center of the hydrophobic core but extends in different directions of the S6 structure following the permutant-specific differences in local loop entropies. The adjustment of the critical folding nucleus to alterations in loop entropies is reflected by a direct correlation between the φ-value change and the accompanying change in local sequence separation.}, url = {http://dx.doi.org/10.1073/pnas.0508863103}, author = {Lindberg, Magnus O. and Haglund, Ellinor and Hubner, Isaac A. and Eugene I. Shakhnovich and Oliveberg, Mikael} } @article {101301, title = {A simple physical model for scaling in protein-protein interaction networks}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {103}, year = {2006}, month = {2006}, pages = {311 - 316}, abstract = {It has recently been demonstrated that many biological networks exhibit a "scale-free{\textquotedblright} topology, for which the probability of observing a node with a certain number of edges (k) follows a power law: i.e., p(k) \~{} k -γ. This observation has been reproduced by evolutionary models. Here we consider the network of protein-protein interactions (PPIs) and demonstrate that two published independent measurements of these interactions produce graphs that are only weakly correlated with one another despite their strikingly similar topology. We then propose a physical model based on the fundamental principle that (de)solvation is a major physical factor in PPIs. This model reproduces not only the scale-free nature of such graphs but also a number of higher-order correlations in these networks. A key support of the model is provided by the discovery of a significant correlation between the number of interactions made by a protein and the fraction of hydrophobic residues on its surface. The model presented in this paper represents a physical model for experimentally determined PPIs that comprehensively reproduces the topological features of interaction networks. These results have profound implications for understanding not only PPIs but also other types of scale-free networks.}, url = {http://dx.doi.org/10.1073/pnas.0509715102}, author = {Deeds, Eric J. and Ashenberg, Orr and Eugene I. Shakhnovich} } @article {101181, title = {Folding and binding: The conformational repertoire of proteins: folding, aggregation and structural recognition}, journal = {Current Opinion in Structural Biology}, volume = {16}, year = {2006}, month = {2006}, pages = {68 - 70}, url = {http://dx.doi.org/10.1016/j.sbi.2006.01.015}, author = {Oliveberg, Mikael and Eugene I. Shakhnovich} } @article {101146, title = {Understanding ensemble protein folding at atomic detail}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {103}, year = {2006}, month = {2006}, pages = {17747 - 17752}, abstract = {It has long been known that a protein{\textquoteright}s amino acid sequence dictates its native structure. However, despite significant recent advances, an ensemble description of how a protein achieves its native conformation from random coil under physiologically relevant conditions remains incomplete. Here we present a detailed all-atom model with a transferable potential that is capable of ab initio folding of entire protein domains using only sequence information. The computational efficiency of this model allows us to perform thousands of microsecond-time scale-folding simulations of the engrailed homeodomain and to observe thousands of complete independent folding events. We apply a graph-theoretic analysis to this massive data set to elucidate which intermediates and intermediary states are common to many trajectories and thus important for the folding process. This method provides an atomically detailed and complete picture of a folding pathway at the ensemble level. The approach that we describe is quite general and could be used to study the folding of proteins on time scales orders of magnitude longer than currently possible.}, url = {http://dx.doi.org/10.1073/pnas.0605580103}, author = {Hubner, Isaac A. and Deeds, Eric J. and Eugene I. Shakhnovich} } @article {101091, title = {Statistically Enhanced Self-Attraction of Random Patterns}, journal = {Phys. Rev. Lett.}, volume = {97}, year = {2006}, month = {2006}, pages = {178101}, abstract = {In this work we develop a theory of interaction of randomly patterned surfaces as a generic prototype model of protein-protein interactions. The theory predicts that pairs of randomly superimposed identical (homodimeric) random patterns have always twice as large magnitude of the energy fluctuations with respect to their mutual orientation, as compared with pairs of different (heterodimeric) random patterns. The amplitude of the energy fluctuations is proportional to the square of the average pattern density, to the square of the amplitude of the potential and its characteristic length, and scales linearly with the area of surfaces. The greater dispersion of interaction energies in the ensemble of homodimers implies that strongly attractive complexes of random surfaces are much more likely to be homodimers, rather than heterodimers. Our findings suggest a plausible physical reason for the anomalously high fraction of homodimers observed in real protein interaction networks.}, url = {http://dx.doi.org/10.1103/PhysRevLett.97.178101}, author = {Lukatsky, D. B. and Zeldovich, K. B. and Shakhnovich, E.I.} } @article {101536, title = {Geometric and physical considerations for realistic protein models}, journal = {Phys. Rev. E}, volume = {72}, year = {2005}, month = {2005}, pages = {022901}, abstract = {Protein structure is generally conceptualized as the global arrangement or of smaller, local motifs of helices, sheets, and loops. These regular, recurring secondary structural elements have well understood and standardized definitions in terms of amino acid backbone geometry and the manner in which hydrogen bonding requirements are satisfied. Recently, "tube{\textquotedblright} models have been proposed to explain protein secondary structure in terms of the geometrically optimal packing of a featureless cylinder. However, atomically detailed simulations demonstrate that such packing considerations alone are insufficient for defining secondary structure; both excluded volume and hydrogen bonding must be explicitly modeled for helix formation. These results have fundamental implications for the construction and interpretation of realistic and meaningful biomacromolecular models.}, url = {http://dx.doi.org/10.1103/PhysRevE.72.022901}, author = {Hubner, Isaac A. and Eugene I. Shakhnovich} } @article {101521, title = {High-resolution protein folding with a transferable potential}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {102}, year = {2005}, month = {2005}, pages = {18914 - 18919}, abstract = {A generalized computational method for folding proteins with a fully transferable potential and geometrically realistic all-atom model is presented and tested on seven helix bundle proteins. The protocol, which includes graph-theoretical analysis of the ensemble of resulting folded conformations, was systematically applied and consistently produced structure predictions of ≈3 {\r A} without any knowledge of the native state. To measure and understand the significance of the results, extensive control simulations were conducted. Graph theoretic analysis provides a means for systematically identifying the native fold and provides physical insight, conceptually linking the results to modern theoretical views of protein folding. In addition to presenting a method for prediction of structure and folding mechanism, our model suggests that an accurate all-atom amino acid representation coupled with a physically reasonable atomic interaction potential and hydrogen bonding are essential features for a realistic protein model.}, url = {http://dx.doi.org/10.1073/pnas.0502181102}, author = {Hubner, Isaac A. and Deeds, Eric J. and Eugene I. Shakhnovich} } @article {101486, title = {Evolutionary dynamics of adult stem cells: Comparison of random and immortal-strand segregation mechanisms}, journal = {Phys. Rev. E}, volume = {71}, year = {2005}, month = {2005}, pages = {041914}, abstract = {This paper develops a point-mutation model describing the evolutionary dynamics of a population of adult stem cells. Such a model may prove useful for quantitative studies of tissue aging and the emergence of cancer. We consider two modes of chromosome segregation: (1) random segregation, where the daughter chromosomes of a given parent chromosome segregate randomly into the stem cell and its differentiating sister cell and (2) "immortal DNA strand{\textquotedblright} co-segregation, for which the stem cell retains the daughter chromosomes with the oldest parent strands. Immortal strand co-segregation is a mechanism, originally proposed by [ Cairns Nature (London) 255 197 (1975)], by which stem cells preserve the integrity of their genomes. For random segregation, we develop an ordered strand pair formulation of the dynamics, analogous to the ordered strand pair formalism developed for quasispecies dynamics involving semiconservative replication with imperfect lesion repair (in this context, lesion repair is taken to mean repair of postreplication base-pair mismatches). Interestingly, a similar formulation is possible with immortal strand co-segregation, despite the fact that this segregation mechanism is age dependent. From our model we are able to mathematically show that, when lesion repair is imperfect, then immortal strand co-segregation leads to better preservation of the stem cell lineage than random chromosome segregation. Furthermore, our model allows us to estimate the optimal lesion repair efficiency for preserving an adult stem cell population for a given period of time. For human stem cells, we obtain that mispaired bases still present after replication and cell division should be left untouched, to avoid potentially fixing a mutation in both DNA strands.}, url = {http://dx.doi.org/10.1103/PhysRevE.71.041914}, author = {Tannenbaum, Emmanuel and Sherley, James L. and Eugene I. Shakhnovich} } @article {101476, title = {Protein structure and evolutionary history determine sequence space topology}, journal = {Genome Res.}, volume = {15}, year = {2005}, month = {2005}, pages = {385 - 392}, abstract = {Understanding the observed variability in the number of homologs of a gene is a very important unsolved problem that has broad implications for research into coevolution of structure and function, gene duplication, pseudogene formation, and possibly for emerging diseases. Here, we attempt to define and elucidate some possible causes behind the observed irregularity in sequence space. We present evidence that sequence variability and functional diversity of a gene or fold family is influenced by quantifiable characteristics of the protein structure. These characteristics reflect the structural potential for sequence plasticity, i.e., the ability to accept mutation without losing thermodynamic stability. We identify a structural feature of a protein domain{\textemdash}contact density{\textemdash}that serves as a determinant of entropy in sequence space, i.e., the ability of a protein to accept mutations without destroying the fold (also known as fold designability). We show that (log) of average gene family size exhibits statistical correlation (R2 > 0.9.) with contact density of its three-dimensional structure. We present evidence that the size of individual gene families are influenced not only by the designability of the structure, but also by evolutionary history, e.g., the amount of time the gene family was in existence. We further show that our observed statistical correlation between gene family size and contact density of the structure is valid on many levels of evolutionary divergence, i.e., not only for closely related sequence, but also for less-related fold and superfamily levels of homology.}, url = {http://dx.doi.org/10.1101/gr.3133605}, author = {Shakhnovich, Boris E. and Deeds, Eric and DeLisi, Charles and Eugene Shakhnovich} } @article {101511, title = {Selective advantage for conservative viruses}, journal = {Phys. Rev. E}, volume = {71}, year = {2005}, month = {2005}, pages = {031903}, abstract = {In this article we study the full semiconservative treatment of a model for the coevolution of a virus and an adaptive immune system. Regions of viability are calculated for both conservatively and semiconservatively replicating viruses interacting with a realistic semiconservatively replicating immune system. The conservative virus is found to have a selective advantage in the form of an ability to survive in regions with a wider range of mutation rates than its semiconservative counterpart, as well as an increased replication rate where both species can survive. This may help explain the existence of a rich range of viruses with conservatively replicating genomes, a trait that is found nowhere else in nature.}, url = {http://dx.doi.org/10.1103/PhysRevE.71.031903}, author = {Brumer, Yisroel and Eugene I. Shakhnovich} } @article {101446, title = {Reconstruction of the src-SH3 Protein Domain Transition State Ensemble using Multiscale Molecular Dynamics Simulations}, journal = {Journal of Molecular Biology}, volume = {350}, year = {2005}, month = {2005}, pages = {1035 - 1050}, abstract = {We use an integrated computational approach to reconstruct accurately the transition state ensemble (TSE) for folding of the src-SH3 protein domain. We first identify putative TSE conformations from free energy surfaces generated by importance sampling molecular dynamics for a fully atomic, solvated model of the src-SH3 protein domain. These putative TSE conformations are then subjected to a folding analysis using a coarse-grained representation of the protein and rapid discrete molecular dynamics simulations. Those conformations that fold to the native conformation with a probability (Pfold) of approximately 0.5, constitute the true transition state. Approximately 20\% of the putative TSE structures were found to have a Pfold near 0.5, indicating that, although correct TSE conformations are populated at the free energy barrier, there is a critical need to refine this ensemble. Our simulations indicate that the true TSE conformations are compact, with a well-defined central β sheet, in good agreement with previous experimental and theoretical studies. A structured central β sheet was found to be present in a number of pre-TSE conformations, however, indicating that this element, although required in the transition state, does not define it uniquely. An additional tight cluster of contacts between highly conserved residues belonging to the diverging turn and second β-sheet of the protein emerged as being critical elements of the folding nucleus. A number of commonly used order parameters to identify the transition state for folding were investigated, with the number of native Cβ contacts displaying the most satisfactory correlation with Pfold values.}, url = {http://dx.doi.org/10.1016/j.jmb.2005.05.017}, author = {Ding, Feng and Guo, Weihua and Dokholyan, Nikolay V. and Eugene I. Shakhnovich and Shea, Joan-Emma} } @article {101266, title = {CoC: a database of universally conserved residues in protein folds}, journal = {Bioinformatics}, volume = {21}, year = {2005}, month = {2005}, pages = {2539 - 2540}, abstract = {Summary: The conservatism of conservatism (CoC) database presents statistically analyzed information about the conservation of residue positions in folds across protein families.}, url = {http://dx.doi.org/10.1093/bioinformatics/bti360}, author = {Donald, Jason E. and Hubner, Isaac A. and Rotemberg, Veronica M. and Eugene I. Shakhnovich and Mirny, Leonid A.} } @article {101286, title = {The Emergence of Scaling in Sequence-Based Physical Models of Protein Evolution}, journal = {Biophysical Journal}, volume = {88}, year = {2005}, month = {2005}, pages = {3905 - 3911}, abstract = {It has recently been discovered that many biological systems, when represented as graphs, exhibit a scale-free topology. One such system is the set of structural relationships among protein domains. The scale-free nature of this and other systems has previously been explained using network growth models that, although motivated by biological processes, do not explicitly consider the underlying physics or biology. In this work we explore a sequence-based model for the evolution protein structures and demonstrate that this model is able to recapitulate the scale-free nature observed in graphs of real protein structures. We find that this model also reproduces other statistical feature of the protein domain graph. This represents, to our knowledge, the first such microscopic, physics-based evolutionary model for a scale-free network of biological importance and as such has strong implications for our understanding of the evolution of protein structures and of other biological networks.}, url = {http://dx.doi.org/10.1529/biophysj.104.051433}, author = {Deeds, Eric J. and Eugene I. Shakhnovich} } @article {101256, title = {Semiconservative replication, genetic repair, and many-gened genomes: Extending the quasispecies paradigm to living systems}, journal = {Physics of Life Reviews}, volume = {2}, year = {2005}, month = {2005}, pages = {290 - 317}, abstract = {Quasispecies theory has emerged as an important tool for modeling the evolutionary dynamics of biological systems. We review recent advances in the field, with an emphasis on the quasispecies dynamics of semiconservatively replicating genomes. Applications to cancer and adult stem cell growth are discussed. Additional topics, such as genetic repair and many-gene genomes, are covered as well.}, url = {http://dx.doi.org/10.1016/j.plrev.2005.08.001}, author = {Tannenbaum, Emmanuel and Eugene I. Shakhnovich} } @article {101141, title = {Determining functional specificity from protein sequences}, journal = {Bioinformatics}, volume = {21}, year = {2005}, month = {2005}, pages = {2629 - 2635}, abstract = {Motivation: Given a large family of homologous protein sequences, many methods can divide the family into smaller groups that correspond to the different functions carried out by proteins within the family. One important problem, however, has been the absence of a general method for selecting an appropriate level of granularity, or size of the groups. Results: We propose a consistent way of choosing the granularity that is independent of the sequence similarity and sequence clustering method used. We study three large, well-investigated protein families: basic leucine zippers, nuclear receptors and proteins with three consecutive C2H2 zinc fingers. Our method is tested against known functional information, the experimentally determined binding specificities, using a simple scoring method. The significance of the groups is also measured by randomizing the data. Finally, we compare our algorithm against a popular method of grouping proteins, the TRIBE-MCL method. In the end, we determine that dividing the families at the proposed level of granularity creates very significant and useful groups of proteins that correspond to the different DNA-binding motifs. We expect that such groupings will be useful in studying not only DNA binding but also other protein interactions.}, url = {http://dx.doi.org/10.1093/bioinformatics/bti396}, author = {Donald, Jason E. and Eugene I. Shakhnovich} } @article {101111, title = {Prokaryotic phylogenies inferred from protein structural domains}, journal = {Genome Res.}, volume = {15}, year = {2005}, month = {2005}, pages = {393 - 402}, abstract = {The determination of the phylogenetic relationships among microorganisms has long relied primarily on gene sequence information. Given that prokaryotic organisms often lack morphological characteristics amenable to phylogenetic analysis, prokaryotic phylogenies, in particular, are often based on sequence data. In this work, we explore a new source of phylogenetic information, the distribution of protein structural domains within fully sequenced prokaryotic genomes. The evolution of the structural domains we use has been studied extensively, allowing us to base our phylogenetic methods on testable theoretical models of structural evolution. We find that the methods that produce reasonable phylogenetic relationships are indeed the methods that are most consistent with theoretical evolutionary models. This work represents, to our knowledge, the first such theoretically motivated phylogeny, as well as the first application of structural information to phylogeny on this scale. Our results have strong implications for the phylogenetic relationships among prokaryotic organisms and for the understanding of protein evolution as a whole.}, url = {http://dx.doi.org/10.1101/gr.3033805}, author = {Deeds, Eric J. and Hennessey, Hooman and Eugene I. Shakhnovich} } @article {101071, title = {Entropic Stabilization of Proteins and Its Proteomic Consequences}, journal = {PLoS Comput Biol}, volume = {1}, year = {2005}, month = {2005}, abstract = {Evolutionary traces of thermophilic adaptation are manifest, on the whole-genome level, in compositional biases toward certain types of amino acids. However, it is sometimes difficult to discern their causes without a clear understanding of underlying physical mechanisms of thermal stabilization of proteins. For example, it is well-known that hyperthermophiles feature a greater proportion of charged residues, but, surprisingly, the excess of positively charged residues is almost entirely due to lysines but not arginines in the majority of hyperthermophilic genomes. All-atom simulations show that lysines have a much greater number of accessible rotamers than arginines of similar degree of burial in folded states of proteins. This finding suggests that lysines would preferentially entropically stabilize the native state. Indeed, we show in computational experiments that arginine-to-lysine amino acid substitutions result in noticeable stabilization of proteins. We then hypothesize that if evolution uses this physical mechanism as a complement to electrostatic stabilization in its strategies of thermophilic adaptation, then hyperthermostable organisms would have much greater content of lysines in their proteomes than comparably sized and similarly charged arginines. Consistent with that, high-throughput comparative analysis of complete proteomes shows extremely strong bias toward arginine-to-lysine replacement in hyperthermophilic organisms and overall much greater content of lysines than arginines in hyperthermophiles. This finding cannot be explained by genomic GC compositional biases or by the universal trend of amino acid gain and loss in protein evolution. We discovered here a novel entropic mechanism of protein thermostability due to residual dynamics of rotamer isomerization in native state and demonstrated its immediate proteomic implications. Our study provides an example of how analysis of a fundamental physical mechanism of thermostability helps to resolve a puzzle in comparative genomics as to why amino acid compositions of hyperthermophilic proteomes are significantly biased toward lysines but not similarly charged arginines.}, url = {http://dx.doi.org/10.1371/journal.pcbi.0010047}, author = {Berezovsky, Igor N and Chen, William W. and Choi, Paul J and Eugene I. Shakhnovich} } @article {101061, title = {Lessons from the design of a novel atomic potential for protein folding}, journal = {Protein Science}, volume = {14}, year = {2005}, month = {2005}, pages = {1741 - 1752}, abstract = {We investigate all-atom potentials of mean force for estimating free energies in protein folding and fold recognition. We search through the space potentials and design novel atomic potentials with a random mixing approximation and a contact-correlated Gaussian approximation of decoy states. We show that the two derived potentials are highly correlated, supporting the use of the random energy model as an accurate statistical description of protein conformational states. The novel atomic potentials perform well in a Z-score and fold decoy recognition test. Furthermore, the designed atomic potential performs slightly and significantly better than atomic potentials derived under a quasi-chemical assumption. While accounting for connectivity correlations between atom types does not improve the performance of the designed potential, we show these correlations lead to ambiguities in the distribution of energetic contributions for atoms on the same residue. Within the confines of the model then, many potentials may exist which stabilize all native folds in subtly different ways. Comparison of different protein conformations under the various atomic potentials reveals both a remarkable degree of correspondence in the estimated free energies and a remarkable degree of correspondence in the identity of the contacts types that make the dominant contributions to the estimated free energies. This consistency may be interpreted as a sign that the design procedure is extracting physically meaningful quantities.}, url = {http://dx.doi.org/10.1110/ps.051440705}, author = {Chen, William W. and Eugene I. Shakhnovich} } @article {101066, title = {Physics and evolution of thermophilic adaptation}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {102}, year = {2005}, month = {2005}, pages = {12742 - 12747}, abstract = {Analysis of structures and sequences of several hyperthermostable proteins from various sources reveals two major physical mechanisms of their thermostabilization. The first mechanism is "structure-based,{\textquotedblright} whereby some hyperthermostable proteins are significantly more compact than their mesophilic homologues, while no particular interaction type appears to cause stabilization; rather, a sheer number of interactions is responsible for thermostability. Other hyperthermostable proteins employ an alternative, "sequence-based{\textquotedblright} mechanism of their thermal stabilization. They do not show pronounced structural differences from mesophilic homologues. Rather, a small number of apparently strong interactions is responsible for high thermal stability of these proteins. High-throughput comparative analysis of structures and complete genomes of several hyperthermophilic archaea and bacteria revealed that organisms develop diverse strategies of thermophilic adaptation by using, to a varying degree, two fundamental physical mechanisms of thermostability. The choice of a particular strategy depends on the evolutionary history of an organism. Proteins from organisms that originated in an extreme environment, such as hyperthermophilic archaea (Pyrococcus furiosus), are significantly more compact and more hydrophobic than their mesophilic counterparts. Alternatively, organisms that evolved as mesophiles but later recolonized a hot environment (Thermotoga maritima) relied in their evolutionary strategy of thermophilic adaptation on "sequence-based{\textquotedblright} mechanism of thermostability. We propose an evolutionary explanation of these differences based on physical concepts of protein designability.}, url = {http://dx.doi.org/10.1073/pnas.0503890102}, author = {Berezovsky, Igor N and Eugene I. Shakhnovich} } @article {101096, title = {Predicting specificity-determining residues in two large eukaryotic transcription factor families}, journal = {Nucl. Acids Res.}, volume = {33}, year = {2005}, month = {2005}, pages = {4455 - 4465}, abstract = {Certain amino acid residues in a protein, when mutated, change the protein{\textquoteright}s function. We present an improved method of finding these specificity-determining positions that uses all the protein sequence data available for a family of homologous proteins. We study in detail two families of eukaryotic transcription factors, basic leucine zippers and nuclear receptors, because of the large amount of sequences and experimental data available. These protein families also have a clear definition of functional specificity: DNA-binding specificity. We compare our results to three other methods, including the evolutionary trace algorithm and a method that depends on orthology relationships. All of the predictions are compared to the available mutational and crystallographic data. We find that our method provides superior predictions of the known specificity-determining residues and also predicts residue positions within these families that deserve further study for their roles in functional specificity.}, url = {http://dx.doi.org/10.1093/nar/gki755}, author = {Donald, Jason E. and Eugene I. Shakhnovich} } @article {101016, title = {Nucleation and the Transition State of the SH3 Domain}, journal = {Journal of Molecular Biology}, volume = {349}, year = {2005}, month = {2005}, pages = {424 - 434}, abstract = {We present a verified computational model of the SH3 domain transition state (TS) ensemble. This model was built for three separate SH3 domains using experimental ϕ-values as structural constraints in all-atom protein folding simulations. While averaging over all conformations incorrectly considers non-TS conformations as transition states, quantifying structures as pre-TS, TS, and post-TS by measurement of their transmission coefficient ("probability to fold{\textquotedblright}, or pfold) allows for rigorous conclusions regarding the structure of the folding nucleus and a full mechanistic analysis of the folding process. Through analysis of the TS, we observe a highly polarized nucleus in which many residues are solvent-exposed. Mechanistic analysis suggests the hydrophobic core forms largely after an early nucleation step. SH3 presents an ideal system for studying the nucleation-condensation mechanism and highlights the synergistic relationship between experiment and simulation in the study of protein folding.}, url = {http://dx.doi.org/10.1016/j.jmb.2005.03.050}, author = {Hubner, Isaac A. and Edmonds, Katherine A. and Eugene I. Shakhnovich} } @article {101466, title = {Simulation, experiment, and evolution: Understanding nucleation in protein S6 folding}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {101}, year = {2004}, month = {2004}, pages = {8354 - 8359}, abstract = {In this study, we explore nucleation and the transition state ensemble of the ribosomal protein S6 using a Monte Carlo (MC) Go model in conjunction with restraints from experiment. The results are analyzed in the context of extensive experimental and evolutionary data. The roles of individual residues in the folding nucleus are identified, and the order of events in the S6 folding mechanism is explored in detail. Interpretation of our results agrees with, and extends the utility of, experiments that shift φ-values by modulating denaturant concentration and presents strong evidence for the realism of the mechanistic details in our MC Go model and the structural interpretation of experimental φ-values. We also observe plasticity in the contacts of the hydrophobic core that support the specific nucleus. For S6, which binds to RNA and protein after folding, this plasticity may result from the conformational flexibility required to achieve biological function. These results present a theoretical and conceptual picture that is relevant in understanding the mechanism of nucleation in protein folding.}, url = {http://dx.doi.org/10.1073/pnas.0401672101}, author = {Hubner, Isaac A. and Oliveberg, Mikael and Eugene I. Shakhnovich} } @article {101496, title = {Solution of the quasispecies model for an arbitrary gene network}, journal = {Phys. Rev. E}, volume = {70}, year = {2004}, month = {2004}, pages = {021903}, abstract = {In this paper, we study the equilibrium behavior of Eigen{\textquoteright}s quasispecies equations for an arbitrary gene network. We consider a genome consisting of N genes, so that the full genome sequence σ may be written as σ=σ1σ2...σN, where σi are sequences of individual genes. We assume a single fitness peak model for each gene, so that gene i has some "master{\textquotedblright} sequence σi,0 for which it is functioning. The fitness landscape is then determined by which genes in the genome are functioning and which are not. The equilibrium behavior of this model may be solved in the limit of infinite sequence length. The central result is that, instead of a single error catastrophe, the model exhibits a series of localization to delocalization transitions, which we term an "error cascade.{\textquotedblright} As the mutation rate is increased, the selective advantage for maintaining functional copies of certain genes in the network disappears, and the population distribution delocalizes over the corresponding sequence spaces. The network goes through a series of such transitions, as more and more genes become inactivated, until eventually delocalization occurs over the entire genome space, resulting in a final error catastrophe. This model provides a criterion for determining the conditions under which certain genes in a genome will lose functionality due to genetic drift. It also provides insight into the response of gene networks to mutagens. In particular, it suggests an approach for determining the relative importance of various genes to the fitness of an organism, in a more accurate manner than the standard "deletion set{\textquotedblright} method. The results in this paper also have implications for mutational robustness and what C.O. Wilke termed "survival of the flattest.{\textquotedblright}}, url = {http://dx.doi.org/10.1103/PhysRevE.70.021903}, author = {Tannenbaum, Emmanuel and Eugene I. Shakhnovich} } @article {101441, title = {Imperfect DNA lesion repair in the semiconservative quasispecies model: Derivation of the Hamming class equations and solution of the single-fitness peak landscape}, journal = {Phys. Rev. E}, volume = {70}, year = {2004}, month = {2004}, pages = {061915}, abstract = {This paper develops a Hamming class formalism for the semiconservative quasispecies equations with imperfect lesion repair, first presented and analytically solved in Y. Brumer and E.I. Shakhnovich (q-bio.GN/0403018, 2004). Starting from the quasispecies dynamics over the space of genomes, we derive an equivalent dynamics over the space of ordered sequence pairs. From this set of equations, we are able to derive the infinite sequence length form of the dynamics for a class of fitness landscapes defined by a master genome. We use these equations to solve for a generalized single-fitness-peak landscape, where the master genome can sustain a maximum number of lesions and remain viable. We determine the mean equilibrium fitness and error threshold for this class of landscapes, and show that when lesion repair is imperfect, semiconservative replication displays characteristics from both conservative replication and semiconservative replication with perfect lesion repair. The work presented here provides a formulation of the model which greatly facilitates the analysis of a relatively broad class of fitness landscapes, and thus serves as a convenient springboard into biological applications of imperfect lesion repair.}, url = {http://dx.doi.org/10.1103/PhysRevE.70.061915}, author = {Tannenbaum, Emmanuel and Sherley, James L. and Eugene I. Shakhnovich} } @article {101451, title = {Solvation versus freezing in a heteropolymer globule}, journal = {Phys. Rev. E}, volume = {70}, year = {2004}, month = {2004}, pages = {021802}, abstract = {We address the response of a random heteropolymer to preferential solvation of certain monomer types at the globule-solvent interface. For each set of monomers that can comprise the molecule{\textquoteright}s surface, we represent the ensemble of allowed configurations by a Gaussian distribution of energy levels, whose mean and variance depend on the set{\textquoteright}s composition. Within such a random energy model, mean surface composition is proportional to solvation strength under most conditions. The breadth of this linear response regime arises from the approximate statistical independence of surface and volume energies. Fluctuations play a crucial role in determining the excess of solvophilic monomers at the surface, and for a diverse set of monomer types can be overcome only by very strong solvent preference.}, url = {http://dx.doi.org/10.1103/PhysRevE.70.021802}, author = {Geissler, Phillip L. and Eugene I. Shakhnovich and Grosberg, Alexander Yu.} } @article {101406, title = {Host-parasite coevolution and optimal mutation rates for semiconservative quasispecies}, journal = {Phys. Rev. E}, volume = {69}, year = {2004}, month = {2004}, pages = {061909}, abstract = {In this paper, we extend a model of host-parasite coevolution to incorporate the semiconservative nature of DNA replication for both the host and the parasite. We find that the optimal mutation rate for the semiconservative and conservative hosts converge for realistic genome lengths, thus maintaining the admirable agreement between theory and experiment found previously for the conservative model and justifying the conservative approximation in some cases. We demonstrate that, while the optimal mutation rate for a conservative and semiconservative parasite interacting with a given immune system is similar to that of a conservative parasite, the properties away from this optimum differ significantly. We suspect that this difference, coupled with the requirement that a parasite optimize survival in a range of viable hosts, may help explain why semiconservative viruses are known to have significantly lower mutation rates than their conservative counterparts.}, url = {http://dx.doi.org/10.1103/PhysRevE.69.061909}, author = {Brumer, Yisroel and Eugene I. Shakhnovich} } @article {101306, title = {Error and repair catastrophes: A two-dimensional phase diagram in the quasispecies model}, journal = {Phys. Rev. E}, volume = {69}, year = {2004}, month = {2004}, pages = {011902}, abstract = {This paper develops a two-gene, single fitness peak model for determining the equilibrium distribution of genotypes in a unicellular population which is capable of genetic damage repair. The first gene, denoted by σvia, yields a viable organism with first-order growth rate constant k>1 if it is equal to some target "master{\textquotedblright} sequence σvia,0. The second gene, denoted by σrep, yields an organism capable of genetic repair if it is equal to some target "master{\textquotedblright} sequence σrep,0. This model is analytically solvable in the limit of infinite sequence length, and gives an equilibrium distribution which depends on μ=Lε, the product of sequence length and per base pair replication error probability, and εr, the probability of repair failure per base pair. The equilibrium distribution is shown to exist in one of the three possible "phases.{\textquotedblright} In the first phase, the population is localized about the viability and repairing master sequences. As εr exceeds the fraction of deleterious mutations, the population undergoes a "repair{\textquotedblright} catastrophe, in which the equilibrium distribution is still localized about the viability master sequence, but is spread ergodically over the sequence subspace defined by the repair gene. Below the repair catastrophe, the distribution undergoes the error catastrophe when μ exceeds lnk/εr, while above the repair catastrophe, the distribution undergoes the error catastrophe when μ exceeds lnk/fdel, where fdel denotes the fraction of deleterious mutations.}, url = {http://dx.doi.org/10.1103/PhysRevE.69.011902}, author = {Tannenbaum, Emmanuel and Eugene I. Shakhnovich} } @article {101271, title = {Native Atom Types for Knowledge-Based Potentials: Application to Binding Energy Prediction}, journal = {J. Med. Chem.}, volume = {47}, year = {2004}, month = {2004}, pages = {4538 - 4558}, abstract = {Knowledge-based potentials have been found useful in a variety of biophysical studies of macromolecules. Recently, it has also been shown in self-consistent studies that it is possible to extract quantities consistent with pair potentials from model structural databases. In this study, we attempt to extend the results obtained from these self-consistent studies toward the extraction of realistic pair potentials from the Protein Data Bank (PDB). The new method utilizes a clustering approach to define atom types within the PDB consistent with the optimal effective pairwise potential. The method has been integrated into the SMoG drug design package, resulting in an improved approach for the rapid and accurate estimation of binding affinities from structural information. Using this approach, it is possible to generate simple knowledge-based potentials that correlate (R = 0.61) with experimental binding affinities in a database of 118 diverse complexes. Furthermore, predictions performed on a random 1/3 of the database consistently show an average unsigned error of 1.5 log Ki units. It is also possible to generate specialized knowledge-based potentials, targeted to specific protein families. This approach is capable of generating potentials that correlate strongly with experimental binding affinities within these families (R = 0.8?0.9). Predictions on 1/3 of these family databases yield average unsigned errors ranging from 1.1 to 1.3 log Ki units. In summary, we describe a physically motivated approach to optimizing knowledge-based potentials for binding energy prediction that can be integrated into a variety of stages within a lead discovery protocol.}, url = {http://dx.doi.org/10.1021/jm0498046}, author = {Dominy, Brian N. and Eugene I. Shakhnovich} } @article {101251, title = {Semiconservative replication in the quasispecies model}, journal = {Phys. Rev. E}, volume = {69}, year = {2004}, month = {2004}, pages = {061916}, abstract = {This paper extends Eigen{\textquoteright}s quasispecies equations to account for the semiconservative nature of DNA replication. We solve the equations in the limit of infinite sequence length for the simplest case of a static, sharply peaked fitness landscape. We show that the error catastrophe occurs when μ, the product of sequence length and per base pair mismatch probability, exceeds 2 ln[2/(1+1/k)], where k>1 is the first-order growth rate constant of the viable "master{\textquotedblright} sequence (with all other sequences having a first-order growth rate constant of 1). This is in contrast to the result of ln k for conservative replication. In particular, as k{\textrightarrow}$\infty$, the error catastrophe is never reached for conservative replication, while for semiconservative replication the critical μ approaches 2 ln 2. Semiconservative replication is therefore considerably less robust than conservative replication to the effect of replication errors. We also show that the mean equilibrium fitness of a semiconservatively replicating system is given by k(2e-μ/2-1) below the error catastrophe, in contrast to the standard result of ke-μ for conservative replication (derived by Kimura and Maruyama in 1966). From this result it is readily shown that semiconservative replication is necessary to account for the observation that, at sufficiently high mutagen concentrations, faster replicating cells will die more quickly than more slowly replicating cells. Thus, in contrast to Eigen{\textquoteright}s original model, the semiconservative quasispecies equations are able to provide a mathematical basis for explaining the efficacy of mutagens as chemotherapeutic agents.}, url = {http://dx.doi.org/10.1103/PhysRevE.69.061916}, author = {Tannenbaum, Emmanuel and Deeds, Eric J. and Eugene I. Shakhnovich} } @article {101186, title = {Commitment and Nucleation in the Protein G Transition State}, journal = {Journal of Molecular Biology}, volume = {336}, year = {2004}, month = {2004}, pages = {745 - 761}, abstract = {An accurate characterization of the transition state ensemble (TSE) is central to furthering our understanding of the protein folding reaction. We have extensively tested a recently reported method for studying a protein{\textquoteright}s TSE, utilizing φ-value data from protein engineering experiments and computational studies as restraints in all-atom Monte Carlo (MC) simulations. The validity of interpreting experimental φ-values as the fraction of native contacts made by a residue in the TSE was explored, revealing that this definition is unable to uniquely specify a TSE. The identification of protein G{\textquoteright}s second hairpin, in both pre and post-transition conformations demonstrates that high experimental φ-values do not guarantee a residue{\textquoteright}s importance in the TSE. An analysis of simulations based on structures restrained by experimental φ-values is necessary to yield this result, which is not obvious from a simplistic interpretation of individual φ-values. The TSE that we obtain corresponds to a single, specific nucleation event, characterized by six residues common to all three observed, convergent folding pathways. The same specific nucleus was independently identified from computational and experimental data, and "Conservation of Conservation{\textquotedblright} analysis in the protein G fold. When associated strictly with complete nucleus formation and concomitant chain collapse, folding is a well-defined two state event. Once the nucleus has formed, the folding reaction enters a slow relaxation process associated with side-chain packing and small, local backbone rearrangements. A detailed analysis of φ-values and their relationship to the transition state ensemble allows us to construct a unified theoretical model of protein G folding.}, url = {http://dx.doi.org/10.1016/j.jmb.2003.12.032}, author = {Hubner, Isaac A. and Shimada, Jun and Eugene I. Shakhnovich} } @article {101206, title = {Proteomic Traces of Speciation}, journal = {Journal of Molecular Biology}, volume = {336}, year = {2004}, month = {2004}, pages = {695 - 706}, abstract = {Recent work has shown that the network of structural similarity between protein domains exhibits a power-law distribution of edges per node. The scale-free nature of this graph, termed the protein domain universe graph or PDUG, may be reproduced via a divergent model of structural evolution. The performance of this model, however, does not preclude the existence of a successful convergent model. To further resolve the issue of protein structural evolution, we explore the predictions of both convergent and divergent models directly. We show that when nodes from the PDUG are partitioned into subgraphs on the basis of their occurrence in the proteomes of particular organisms, these subgraphs exhibit a scale-free nature as well. We explore a simple convergent model of structural evolution and find that the implications of this model are inconsistent with features of these organismal subgraphs. Importantly, we find that biased convergent models are inconsistent with our data. We find that when speciation mechanisms are added to a simple divergent model, subgraphs similar to the organismal subgraphs are produced, demonstrating that dynamic models can easily explain the distributions of structural similarity that exist within proteomes. We show that speciation events must be included in a divergent model of structural evolution to account for the non-random overlap of structural proteomes. These findings have implications for the long-standing debate over convergent and divergent models of protein structural evolution, and for the study of the evolution of organisms as a whole.}, url = {http://dx.doi.org/10.1016/j.jmb.2003.12.066}, author = {Deeds, Eric J. and Shakhnovich, Boris and Eugene I. Shakhnovich} } @article {101121, title = {The evolution dynamics of model proteins}, journal = {J. Chem. Phys.}, volume = {121}, year = {2004}, month = {2004}, pages = {2381}, abstract = {Explicit simulations of protein evolution, where protein chains are described at a molecular, although simplified, level provide important information to understand the similarities found to exist between known proteins. The results of such simulations suggest that a number of evolutionary-related quantities, such as the distribution of sequence similarity for structurally similar proteins, are controlled by evolutionary kinetics and do not reflect an equilibrium state. An important result for phylogeny is that a subset of the residues of each protein evolve on a much larger time scale than the other residues. {\textcopyright} 2004 American Institute of Physics.}, url = {http://dx.doi.org/doi:10.1063/1.1768513}, author = {Tiana, Guido and Dokholyan, Nikolay V. and Broglia, Ricardo A. and Eugene I. Shakhnovich} } @article {101116, title = {Importance of DNA repair in tumor suppression}, journal = {Phys. Rev. E}, volume = {70}, year = {2004}, month = {2004}, pages = {061912}, abstract = {The transition from a normal to cancerous cell requires a number of highly specific mutations that affect cell cycle regulation, apoptosis, differentiation, and many other cell functions. One hallmark of cancerous genomes is genomic instability, with mutation rates far greater than those of normal cells. In microsatellite instability (MIN tumors), these are often caused by damage to mismatch repair genes, allowing further mutation of the genome and tumor progression. These mutation rates may lie near the error catastrophe found in the quasispecies model of adaptive RNA genomes, suggesting that further increasing mutation rates will destroy cancerous genomes. However, recent results have demonstrated that DNA genomes exhibit an error threshold at mutation rates far lower than their conservative counterparts. Furthermore, while the maximum viable mutation rate in conservative systems increases indefinitely with increasing master sequence fitness, the semiconservative threshold plateaus at a relatively low value. This implies a paradox, wherein inaccessible mutation rates are found in viable tumor cells. In this paper, we address this paradox, demonstrating an isomorphism between the conservatively replicating (RNA) quasispecies model and the semiconservative (DNA) model with post-methylation DNA repair mechanisms impaired. Thus, as DNA repair becomes inactivated, the maximum viable mutation rate increases smoothly to that of a conservatively replicating system on a transformed landscape, with an upper bound that is dependent on replication rates. On a specific single fitness peak landscape, the repair-free semiconservative system is shown to mimic a conservative system exactly. We postulate that inactivation of post-methylation repair mechanisms is fundamental to the progression of a tumor cell and hence these mechanisms act as a method for the prevention and destruction of cancerous genomes.}, url = {http://dx.doi.org/10.1103/PhysRevE.70.061912}, author = {Brumer, Yisroel and Eugene I. Shakhnovich} } @article {101011, title = {Imprint of evolution on protein structures}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {101}, year = {2004}, month = {2004}, pages = {2846 - 2851}, abstract = {We attempt to understand the evolutionary origin of protein folds by simulating their divergent evolution with a three-dimensional lattice model. Starting from an initial seed lattice structure, evolution of model proteins progresses by sequence duplication and subsequent point mutations. A new gene{\textquoteright}s ability to fold into a stable and unique structure is tested each time through direct kinetic folding simulations. Where possible, the algorithm accepts the new sequence and structure and thus a "new protein structure{\textquotedblright} is born. During the course of each run, this model evolutionary algorithm provides several thousand new proteins with diverse structures. Analysis of evolved structures shows that later evolved structures are more designable than seed structures as judged by recently developed structural determinant of protein designability, as well as direct estimate of designability for selected structures by thermodynamic sampling of their sequence space. We test the significance of this trend predicted on lattice models on real proteins and show that protein domains that are found in eukaryotic organisms only feature statistically significant higher designability than their prokaryotic counterparts. These results present a fundamental view on protein evolution highlighting the relative roles of structural selection and evolutionary dynamics on genesis of modern proteins.}, url = {http://dx.doi.org/10.1073/pnas.0306638101}, author = {Tiana, Guido and Shakhnovich, Boris E. and Dokholyan, Nikolay V. and Eugene I. Shakhnovich} } @article {100976, title = {All-atom Monte Carlo Simulation of GCAA RNA Folding}, journal = {Journal of Molecular Biology}, volume = {344}, year = {2004}, month = {2004}, pages = {29 - 45}, abstract = {We report a detailed all-atom simulation of the folding of the GCAA RNA tetraloop. The GCAA tetraloop motif is a very common and thermodynamically stable secondary structure in natural RNAs. We use our simulation methods to study the folding behavior of a 12-base GCAA tetraloop structure with a four-base helix adjacent to the tetraloop proper. We implement an all-atom Monte Carlo (MC) simulation of RNA structural dynamics using a Go potential. Molecular dynamics (MD) simulation of RNA and protein has realistic energetics and sterics, but is extremely expensive in terms of computational time. By coarsely treating non-covalent energetics, but retaining all-atom sterics and entropic effects, all-atom MC techniques are a useful method for the study of protein and now RNA. We observe a sharp folding transition for this structure, and in simulations at room temperature the state histogram shows three distinct minima: an unfolded state (U), a more narrow intermediated state (I), and a narrow folded state (F). The intermediate consists primarily of structures with the GCAA loop and some helix hydrogen bonds formed. Repeated kinetic folding simulations reveal that the number of helix base-pairs forms a simple 1D reaction coordinate for the I{\textrightarrow}N transition.}, url = {http://dx.doi.org/10.1016/j.jmb.2004.09.041}, author = {Niv{\'o}n, Lucas G. and Eugene I. Shakhnovich} } @article {100966, title = {The role of cotranslation in protein folding: a lattice model study}, journal = {Polymer}, volume = {45}, year = {2004}, month = {2004}, pages = {557 - 571}, abstract = {Computational studies of protein folding have implicitly assumed that folding occurs from a denatured state comprised of the entire protein. Cotranslational folding accounts for the linear production and release of a protein from the ribosome, allowing part of the protein to explore its conformation space before other parts have been synthesized. This gradual {\textquoteleft}extrusion{\textquoteright} from the ribosome can yield different folding kinetics than direct folding from the denatured state, for a lattice folding model. First, in model proteins containing chiefly short-ranged (local in sequence) contacts, cotranslational folding is shown to be significantly faster than direct folding from the denatured state. Secondly, for model proteins with two competing native states, cotranslational folding tilts the apparent equilibrium toward the state with a more local-contact dominant topology.}, url = {http://dx.doi.org/10.1016/j.polymer.2003.10.090}, author = {Morrissey, M.P. and Ahmed, Z. and Shakhnovich, E.I.} } @article {87106, title = {Structural Mining: Self-Consistent Design on Flexible Protein-Peptide Docking and Transferable Binding Affinity Potential}, journal = {J. Am. Chem. Soc.}, volume = {126}, year = {2004}, month = {2004}, pages = {8515 - 8528}, abstract = {A flexible protein?peptide docking method has been designed to consider not only ligand flexibility but also the flexibility of the protein. The method is based on a Monte Carlo annealing process. Simulations with a distance root-mean-square (dRMS) virtual energy function revealed that the flexibility of protein side chains was as important as ligand flexibility for successful protein?peptide docking. On the basis of mean field theory, a transferable potential was designed to evaluate distance-dependent protein?ligand interactions and atomic solvation energies. The potential parameters were developed using a self-consistent process based on only 10 known complex structures. The effectiveness of each intermediate potential was judged on the basis of a Z score, approximating the gap between the energy of the native complex and the average energy of a decoy set. The Z score was determined using experimentally determined native structures and decoys generated by docking with the intermediate potentials. Using 6600 generated decoys and the Z score optimization criterion proposed in this work, the developed potential yielded an acceptable correlation of R2 = 0.77, with binding free energies determined for known MHC I complexes (Class I Major Histocompatibility protein HLA-A*0201) which were not present in the training set. Test docking on 25 complexes further revealed a significant correlation between energy and dRMS, important for identifying native-like conformations. The near-native structures always belonged to one of the conformational classes with lower predicted binding energy. The lowest energy docked conformations are generally associated with near-native conformations, less than 3.0 {\r A} dRMS (and in many cases less than 1.0 {\r A}) from the experimentally determined structures.}, url = {http://dx.doi.org/10.1021/ja032018q}, author = {Liu, Zhijie and Dominy, Brian N. and Eugene I. Shakhnovich} } @article {101506, title = {Fold recognition with minimal gaps}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {51}, year = {2003}, month = {2003}, pages = {531 - 543}, abstract = {Here we present a simplified form of threading that uses only a 20 {\texttimes} 20 two-body residue-based potential and restricted number of gaps. Despite its simplicity and transparency the Monte Carlo-based threading algorithm performs very well in a rigorous test of fold recognition. The results suggest that by simplifying and constraining the decoy space, one can achieve better fold recognition. Fold recognition results are compared with and supplemented by a PSI-BLAST search. The statistical significance of threading results is rigorously evaluated from statistics of extremes by comparison with optimal alignments of a large set of randomly shuffled sequences. The statistical theory, based on the Random Energy Model, yields a cumulative statistical parameter, ϵ, that attests to the likelihood of correct fold recognition. A large ϵ indicates a significant energy gap between the optimal alignment and decoy alignments and, consequently, a high probability that the fold is correctly recognized. For a particular number of gaps, the ϵ parameter reaches its maximal value, and the fold is recognized. As the number of gaps further increases, the likelihood of correct fold recognition drops off. This is because the decoy space is small when gaps are restricted to a small number, but the native alignment is still well approximated, whereas unrestricted increase of the number of gaps leads to rapid growth of the number of decoys and their statistical dominance over the correct alignment. It is shown that best results are obtained when a combination of one-, two-, and three-gap threading is used. To this end, use of the ϵ parameter is crucial for rigorous comparison of results across the different decoy spaces belonging to a different number of gaps. Proteins 2003;51:531{\textendash}543. {\textcopyright} 2003 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/prot.10402}, author = {Chen, William and Mirny, Leonid and Eugene I. Shakhnovich} } @article {101461, title = {ELISA: Structure-Function Inferences based on statistically significant and evolutionarily inspired observations}, journal = {BMC Bioinformatics}, volume = {4}, year = {2003}, month = {2003}, abstract = {The problem of functional annotation based on homology modeling is primary to current bioinformatics research. Researchers have noted regularities in sequence, structure and even chromosome organization that allow valid functional cross-annotation. However, these methods provide a lot of false negatives due to limited specificity inherent in the system. We want to create an evolutionarily inspired organization of data that would approach the issue of structure-function correlation from a new, probabilistic perspective. Such organization has possible applications in phylogeny, modeling of functional evolution and structural determination. ELISA (Evolutionary Lineage Inferred from Structural Analysis, http://romi.bu.edu/elisa) is an online database that combines functional annotation with structure and sequence homology modeling to place proteins into sequence-structure-function}, url = {http://dx.doi.org/10.1186/1471-2105-4-34}, author = {Shakhnovich, Boris E. and Harvey, John M. and Comeau, Steve and Lorenz, David and DeLisi, Charles and Eugene Shakhnovich} } @article {101426, title = {Protein Evolution within a Structural Space}, journal = {Biophysical Journal}, volume = {85}, year = {2003}, month = {2003}, pages = {2962 - 2972}, abstract = {Understanding of the evolutionary origins of protein structures represents a key component of the understanding of molecular evolution as a whole. Here we seek to elucidate how the features of an underlying protein structural "space{\textquotedblright} might impact protein structural evolution. We approach this question using lattice polymers as a completely characterized model of this space. We develop a measure of structural comparison of lattice structures that is analogous to the one used to understand structural similarities between real proteins. We use this measure of structural relatedness to create a graph of lattice structures and compare this graph (in which nodes are lattice structures and edges are defined using structural similarity) to the graph obtained for real protein structures. We find that the graph obtained from all compact lattice structures exhibits a distribution of structural neighbors per node consistent with a random graph. We also find that subgraphs of 3500 nodes chosen either at random or according to physical constraints also represent random graphs. We develop a divergent evolution model based on the lattice space which produces graphs that, within certain parameter regimes, recapitulate the scale-free behavior observed in similar graphs of real protein structures.}, url = {http://dx.doi.org/10.1016/S0006-3495(03)74716-X}, author = {Deeds, Eric J. and Dokholyan, Nikolay V. and Eugene I. Shakhnovich} } @article {101381, title = {Equilibrium Distribution of Mutators in the Single Fitness Peak Model}, journal = {Phys. Rev. Lett.}, volume = {91}, year = {2003}, month = {2003}, pages = {138105}, abstract = {This Letter develops an analytically tractable model for determining the equilibrium distribution of mismatch repair deficient strains in unicellular populations. The approach is based on the single fitness peak model, which has been used in Eigen{\textquoteright}s quasispecies equations in order to understand various aspects of evolutionary dynamics. As with the quasispecies model, our model for mutator-nonmutator equilibrium undergoes a phase transition in the limit of infinite sequence length. This "repair catastrophe{\textquotedblright} occurs at a critical repair error probability of ϵr=Lvia/L, where Lvia denotes the length of the genome controlling viability, while L denotes the overall length of the genome. The repair catastrophe therefore occurs when the repair error probability exceeds the fraction of deleterious mutations. Our model also gives a quantitative estimate for the equilibrium fraction of mutators in Escherichia coli.}, url = {http://dx.doi.org/10.1103/PhysRevLett.91.138105}, author = {Tannenbaum, Emmanuel and Deeds, Eric J. and Eugene I. Shakhnovich} } @article {101371, title = {Structural Determinant of Protein Designability}, journal = {Phys. Rev. Lett.}, volume = {90}, year = {2003}, month = {2003}, pages = {218101}, abstract = {Here we present an approximate analytical theory for the relationship between a protein structure{\textquoteright}s contact matrix and the shape of its energy spectrum in amino acid sequence space. We demonstrate a dependence of the number of sequences of low energy in a structure on the eigenvalues of the structure{\textquoteright}s contact matrix, and then use a Monte Carlo simulation to test the applicability of this analytical result to cubic lattice proteins. We find that the lattice structures with the most low-energy sequences are the same as those predicted by the theory. We argue that, given sufficiently strict requirements for foldability, these structures are the most designable, and we propose a simple means to test whether the results in this paper hold true for real proteins.}, url = {http://dx.doi.org/10.1103/PhysRevLett.90.218101}, author = {England, Jeremy L. and Eugene I. Shakhnovich} } @article {101356, title = {Using protein design for homology detection and active site searches}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {100}, year = {2003}, month = {2003}, pages = {11361 - 11366}, abstract = {We describe a method of designing artificial sequences that resemble naturally occurring sequences in terms of their compatibility with a template structure and its functional constraints. The design procedure is a Monte Carlo simulation of amino acid substitution process. The selective fixation of substitutions is dictated by a simple scoring function derived from the template structure and a multiple alignment of its homologs. Designed sequences represent an enlargement of sequence space around native sequences. We show that the use of designed sequences improves the performance of profile-based homology detection. The difference in position-specific conservation between designed sequences and native sequences is helpful for prediction of functionally important residues. Our sequence selection criteria in evolutionary simulations introduce amino acid substitution rate variation among sites in a natural way, providing a better model to test phylogenetic methods.}, url = {http://dx.doi.org/10.1073/pnas.2034878100}, author = {Pei, Jimin and Dokholyan, Nikolay V. and Eugene I. Shakhnovich and Grishin, Nick V.} } @article {101241, title = {Amino acids determining enzyme-substrate specificity in prokaryotic and eukaryotic protein kinases}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {100}, year = {2003}, month = {2003}, pages = {4463 - 4468}, abstract = {The binding between a PK and its target is highly specific, despite the fact that many different PKs exhibit significant sequence and structure homology. There must be, then, specificity-determining residues (SDRs) that enable different PKs to recognize their unique substrate. Here we use and further develop a computational procedure to discover putative SDRs (PSDRs) in protein families, whereby a family of homologous proteins is split into orthologous proteins, which are assumed to have the same specificity, and paralogous proteins, which have different specificities. We reason that PSDRs must be similar among orthologs, whereas they must necessarily be different among paralogs. Our statistical procedure and evolutionary model identifies such residues by discriminating a functional signal from a phylogenetic one. As case studies we investigate the prokaryotic two-component system and the eukaryotic AGC (i.e., cAMP-dependent PK, cGMP-dependent PK, and PKC) PKs. Without using experimental data, we predict PSDRs in prokaryotic and eukaryotic PKs, and suggest precise mutations that may convert the specificity of one PK to another. We compare our predictions with current experimental results and obtain considerable agreement with them. Our analysis unifies much of existing data on PK specificity. Finally, we find PSDRs that are outside the active site. Based on our results, as well as structural and biochemical characterizations of eukaryotic PKs, we propose the testable hypothesis of "specificity via differential activation{\textquotedblright} as a way for the cell to control kinase specificity.}, url = {http://dx.doi.org/10.1073/pnas.0737647100}, author = {Li, Lewyn and Eugene I. Shakhnovich and Mirny, Leonid A.} } @article {101226, title = {Side-chain dynamics and protein folding}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {52}, year = {2003}, month = {2003}, pages = {303 - 321}, abstract = {The processes by which protein side chains reach equilibrium during a folding reaction are investigated using both lattice and all-atom simulations. We find that rates of side-chain relaxation exhibit a distribution over the protein structure, with the fastest relaxing side chains located in positions kinetically important for folding. Traversal of the major folding transition state corresponds to the freezing of a small number of side chains, belonging to the folding nucleus, whereas the rest of the protein proceeds toward equilibrium via backbone fluctuations around the native fold. The postnucleation processes by which side chains relax are characterized by very slow dynamics and many barrier crossings, and thus resemble the behavior of a glass. Proteins 2003;52:303{\textendash}321. {\textcopyright} 2003 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/prot.10426}, author = {Kussell, Edo and Shimada, Jun and Eugene I. Shakhnovich} } @article {101201, title = {Functional Fingerprints of Folds: Evidence for Correlated Structure{\textendash}Function Evolution}, journal = {Journal of Molecular Biology}, volume = {326}, year = {2003}, month = {2003}, pages = {1 - 9}, abstract = {Using structural similarity clustering of protein domains: protein domain universe graph (PDUG), and a hierarchical functional annotation: gene ontology (GO) as two evolutionary lenses, we find that each structural cluster (domain fold) exhibits a distribution of functions that is unique to it. These functional distributions are functional fingerprints that are specific to characteristic structural clusters and vary from cluster to cluster. Furthermore, as structural similarity threshold for domain clustering in the PDUG is relaxed we observe an influx of earlier-diverged domains into clusters. These domains join clusters without destroying the functional fingerprint. These results can be understood in light of a divergent evolution scenario that posits correlated divergence of structural and functional traits in protein domains from one or few progenitors.}, url = {http://dx.doi.org/10.1016/S0022-2836(02)01362-1}, author = {Shakhnovich, Boris E. and Dokholyan, Nikolay V. and DeLisi, Charles and Eugene I. Shakhnovich} } @inbook {101106, title = {Identifying Importance of Amino Acids for Protein Folding from Crystal Structures}, booktitle = {Methods in Enzymology}, volume = {374}, year = {2003}, month = {2003}, pages = {616 - 638}, abstract = {This chapter presents an overview of computational techniques for reconstructing the folding mechanisms of proteins from their crystal structures. The chapter describes analytical and computational tools for determining and characterizing protein-folding kinetics from crystal structures. It also discusses new protein model and shows that the thermodynamics of Src SH3 from molecular dynamics simulations is consistent with that observed experimentally. The revolution in protein crystallography has resulted in the identification of a large number of protein structures. The chapter discusses several studies that are based on the Go model. In one such study, it identifies the most evasive protein-folding transition state ensemble for Src SH3 protein and finds that it is consistent with experimental observations. These studies suggest that protein-folding kinetics can be determined to a reasonably detailed level from the knowledge of crystal structure. The chapter also dissects the transition state ensemble by studying the wiring properties of protein graphs. The structural properties of protein graphs are related to protein topology and thus may explain the kinetics of the folding process. These studies unveil the expanding possibilities for studying protein-folding kinetics from their crystal structures.}, url = {http://dx.doi.org/10.1016/S0076-6879(03)74025-7}, author = {Dokholyan, Nikolay V. and Borreguero, Jose M. and Buldyrev, Sergey V. and Ding, Feng and Stanley, H.Eugene and Eugene I. Shakhnovich} } @article {87101, title = {Natural selection of more designable folds: A mechanism for thermophilic adaptation}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {100}, year = {2003}, month = {2003}, pages = {8727 - 8731}, abstract = {An open question of great interest in biophysics is whether variations in structure cause protein folds to differ in the number of amino acid sequences that can fold to them stably, i.e., in their designability. Recently, we have shown that a novel quantitative measure of a fold{\textquoteright}s tertiary topology, called its contact trace, strongly correlates with the fold{\textquoteright}s designability. Here, we investigate the relationship between a fold{\textquoteright}s contact trace and its relative frequency of usage in mesophilic vs. thermophilic eubacteria. We observe that thermophilic organisms exhibit a bias toward using folds of higher contact trace when compared with mesophiles. We establish this difference both for the distributions of folds at the whole-proteome level and also through more focused structural comparisons of orthologous proteins. Our findings suggest that thermophilic adaptation in bacterial genomes occurs in part through natural selection of more designable folds, pointing to designability as a key component of protein fitness.}, url = {http://dx.doi.org/10.1073/pnas.1530713100}, author = {England, Jeremy L. and Shakhnovich, Boris E. and Eugene I. Shakhnovich} } @article {101531, title = {Reversible stretching of random heteropolymers}, journal = {Phys. Rev. E}, volume = {65}, year = {2002}, month = {2002}, pages = {056110}, abstract = {We analyze the equilibrium response of random heteropolymers to mechanical deformation. In contrast to homopolymer response, the stress-induced transformation of a heteropolymer from globule to coil need not be sharp. For chain lengths relevant to biological macromolecules, intermediate necklacelike structures dominate over a range of applied force. Stability of these conformations is primarily a consequence of solvation: In a typical necklace, relatively solvophilic regions of the chain are extended, while solvophobic regions remain compact. In the long-chain limit, homopolymeric behavior is recovered. Our results suggest that only select polypeptide sequences should unfold reproducibly at a specific force, explaining recent experimental observations.}, url = {http://dx.doi.org/10.1103/PhysRevE.65.056110}, author = {Geissler, Phillip L. and Eugene I. Shakhnovich} } @article {101401, title = {Direct Molecular Dynamics Observation of Protein Folding Transition State Ensemble}, journal = {Biophysical Journal}, volume = {83}, year = {2002}, month = {2002}, pages = {3525 - 3532}, abstract = {The concept of the protein transition state ensemble (TSE), a collection of the conformations that have 50\% probability to convert rapidly to the folded state and 50\% chance to rapidly unfold, constitutes the basis of the modern interpretation of protein engineering experiments. It has been conjectured that conformations constituting the TSE in many proteins are the expanded and distorted forms of the native state built around a specific folding nucleus. This view has been supported by a number of on-lattice and off-lattice simulations. Here we report a direct observation and characterization of the TSE by molecular dynamic folding simulations of the C-Src SH3 domain, a small protein that has been extensively studied experimentally. Our analysis reveals a set of key interactions between residues, conserved by evolution, that must be formed to enter the kinetic basin of attraction of the native state.}, url = {http://dx.doi.org/10.1016/S0006-3495(02)75352-6}, author = {Ding, Feng and Dokholyan, Nikolay V. and Buldyrev, Sergey V. and Stanley, H. Eugene and Eugene I. Shakhnovich} } @article {101341, title = {From Knowledge-Based Potentials to Combinatorial Lead Design in Silico}, journal = {Acc. Chem. Res.}, volume = {35}, year = {2002}, month = {2002}, pages = {261 - 269}, abstract = {Computational methods are becoming increasingly used in the drug discovery process. In this Account, we review a novel computational method for lead discovery. This method, called CombiSMoG for ?combinatorial small molecule growth?, is based on two components:? a fast and accurate knowledge-based scoring function used to predict binding affinities of protein?ligand complexes, and a Monte Carlo combinatorial growth algorithm that generates large numbers of low-free-energy ligands in the binding site of a protein. We illustrate the advantages of the method by describing its application in the design of picomolar inhibitors for human carbonic anhydrase.}, url = {http://dx.doi.org/10.1021/ar970146b}, author = {Bartosz A. Grzybowski and Alexey V. Ishchenko and Shimada, Jun and Eugene I. Shakhnovich} } @article {101351, title = {Glassy Dynamics of Side-Chain Ordering in a Simple Model of Protein Folding}, journal = {Phys. Rev. Lett.}, volume = {89}, year = {2002}, month = {2002}, pages = {168101}, abstract = {We introduce a modified version of protein lattice models in which monomers have several spin states, representing side-chain rotamers. Completion of folding corresponds to reaching the native backbone configuration with complete ordering of side chains. We find that as temperature is lowered, side-chain ordering becomes much slower than backbone folding. The presence of side chains leads to nonexponential kinetics and a broad distribution of relaxation times.}, url = {http://dx.doi.org/10.1103/PhysRevLett.89.168101}, author = {Kussell, Edo and Eugene I. Shakhnovich} } @article {101361, title = {Understanding conserved amino acids in proteins}, journal = {Physica A: Statistical Mechanics and its Applications}, volume = {314}, year = {2002}, month = {2002}, pages = {600 - 606}, abstract = {It has been conjectured that evolution exerted pressure to preserve amino acids bearing thermodynamic, kinetic, and functional roles. In this letter we show that the physical requirement to maintain protein stability gives rise to a sequence conservatism pattern that is in remarkable agreement with that found in natural proteins. Based on the physical properties of amino acids, we propose a model of evolution that explains conserved amino acids across protein families sharing the same fold.}, url = {http://dx.doi.org/10.1016/S0378-4371(02)01079-8}, author = {Dokholyan, Nikolay V. and Mirny, Leonid A. and Eugene I. Shakhnovich} } @article {101231, title = {Molecular Dynamics Simulation of the SH3 Domain Aggregation Suggests a Generic Amyloidogenesis Mechanism}, journal = {Journal of Molecular Biology}, volume = {324}, year = {2002}, month = {2002}, pages = {851 - 857}, abstract = {We use molecular dynamics simulation to study the aggregation of Src SH3 domain proteins. For the case of two proteins, we observe two possible aggregation conformations: the closed form dimer and the open aggregation state. The closed dimer is formed by "domain swapping{\textquotedblright}{\textemdash}the two proteins exchange their RT-loops. All the hydrophobic residues are buried inside the dimer so proteins cannot further aggregate into elongated amyloid fibrils. We find that the open structure{\textemdash}stabilized by backbone hydrogen bond interactions{\textemdash}packs the RT-loops together by swapping the two strands of the RT-loop. The packed RT-loops form a β-sheet structure and expose the backbone to promote further aggregation. We also simulate more than two proteins, and find that the aggregate adopts a fibrillar double β-sheet structure, which is formed by packing the RT-loops from different proteins. Our simulations are consistent with a possible generic amyloidogenesis scenario.}, url = {http://dx.doi.org/10.1016/S0022-2836(02)01112-9}, author = {Ding, Feng and Dokholyan, Nikolay V. and Buldyrev, Sergey V. and Stanley, H.Eugene and Eugene I. Shakhnovich} } @article {101221, title = {Phase diagram for unzipping DNA with long-range interactions}, journal = {Phys. Rev. E}, volume = {66}, year = {2002}, month = {2002}, pages = {032901}, abstract = {We present a critique and extension of the mean-field approach to the mechanical pulling transition in bound polymer systems. Our model is motivated by the theoretically and experimentally important examples of adsorbed polymers and double-stranded DNA. We focus on the case in which quenched disorder in the sequence of monomers is unimportant for the statistical mechanics, but we include excluded volume interactions between monomers. Our phase diagram for the critical pulling force shows an interesting reentrant phase at low temperatures which should be observable in both disordered and homogenous polymer systems. We also consider the case of nonequilibrium pulling, in which the external force probes the molecule{\textquoteright}s local, rather than global structure. The dynamics of the pulling transition in such experiments could illuminate the polymer{\textquoteright}s loop structure, which depends on the nature of excluded volume interactions.}, url = {http://dx.doi.org/10.1103/PhysRevE.66.032901}, author = {Mukamel, Eran A. and Eugene I. Shakhnovich} } @article {101236, title = {SMall Molecule Growth 2001 (SMoG2001): An Improved Knowledge-Based Scoring Function for Protein-Ligand Interactions}, journal = {J. Med. Chem.}, volume = {45}, year = {2002}, month = {2002}, pages = {2770 - 2780}, abstract = {Computational lead design procedures require fast and accurate scoring functions to rank millions of generated virtual ligands for protein targets. In this article, we present an improved version of the SMoG scoring function, called SMoG2001. This function is based on a knowledge-based approachthat is, the free energy parameters are derived from the observed frequencies of atom?atom contacts in the database of three-dimensional structures of protein?ligand complexes via a procedure based on statistical mechanics. We obtained the statistics from the set of 725 complexes. SMoG2001 reproduces the experimental binding constants of the majority of 119 complexes of the testing set with good accuracy. On similar testing sets, SMoG2001 performs better than two other widely used scoring functions, PMF and SCORE1(LUDI), and comparably to DrugScore. SMoG2001 poorly predicts the affinities of ligands interacting via quantum mechanical forces with metal ions and ligands that are large and flexible. We attribute significant improvement in accuracy over previous versions of the SMoG scoring function to a better description of the reference statethat is, the state of no interactions.}, url = {http://dx.doi.org/10.1021/jm0105833}, author = {Alexey V. Ishchenko and Eugene I. Shakhnovich} } @article {101191, title = {The ensemble folding kinetics of protein G from an all-atom Monte Carlo simulation}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {99}, year = {2002}, month = {2002}, pages = {11175 - 11180}, abstract = {Protein G is folded with an all-atom Monte Carlo simulation by using a Go potential. When folding is monitored by using burial of the lone tryptophan in protein G as the reaction coordinate, the ensemble kinetics is single exponential. Other experimental observations, such as the burst phase and mutational data, are also reproduced. However, more detailed analysis reveals that folding occurs over three distinct, three-state pathways. We show that, because of this tryptophan{\textquoteright}s asymmetric location in the tertiary fold, its burial (i) does not detect certain intermediates and (ii) may not correspond to the folding event. This finding demonstrates that ensemble averaging can disguise the presence of multiple pathways and intermediates when a non-ideal reaction coordinate is used. Finally, all observed folding pathways eventually converge to a common rate-limiting step, which is the formation of a specific nucleus involving hydrophobic core residues. These residues are conserved in the ubiquitin superfamily and in a phage display experiment, suggesting that fold topology is a strong determinant of the transition state.}, url = {http://dx.doi.org/10.1073/pnas.162268099}, author = {Shimada, Jun and Eugene I. Shakhnovich} } @article {101211, title = {Mechanical Response of Random Heteropolymers}, journal = {Macromolecules}, volume = {35}, year = {2002}, month = {2002}, pages = {4429 - 4436}, abstract = {We present an analytical theory for heteropolymer deformation, as exemplified experimentally by stretching of single protein molecules. Using replica mean-field theory, we determine phase diagrams for stress-induced unfolding of typical random sequences. This transition is sharp in the limit of infinitely long chain molecules. However, for chain lengths relevant to biological macromolecules, partially unfolded conformations prevail over an intermediate range of stress. These necklacelike structures, comprised of alternating compact and extended subunits, are stabilized by quenched variations in the composition of finite chain segments. The most stable arrangements of these subunits are largely determined by preferential extension of segments rich in solvophilic monomers. This predicted significance of necklace structures explains recent observations in protein stretching experiments. We examine the statistical features of select sequences that give rise to mechanical strength and may thus have guided the evolution of proteins that carry out mechanical functions in living cells.}, url = {http://dx.doi.org/10.1021/ma012008e}, author = {Geissler, Phillip L. and Eugene I. Shakhnovich} } @article {101161, title = {Field theory and segmental alignment analysis for a solution of sequence disordered liquid crystalline polymers}, journal = {J. Chem. Phys.}, volume = {116}, year = {2002}, month = {2002}, pages = {3134}, abstract = {A field theory for the many-chain mesogen/flexible sequence model is developed. The theory is solved using a creation{\textendash}annihilation formalism based on a reference state that recovers the Gupta{\textendash}Edwards solution for the many-chain homopolymer system. Orientational averages are performed, and a nonlinear equation for the orientational tensor-dependent free energy is derived. The effects of flexibility, stiffness and intersegment alignment on orientational ordering, the nematic/isotropic density threshold, and the segmental orientational ordering at the nematic/isotropic transition is determined. {\textcopyright} 2002 American Institute of Physics.}, url = {http://dx.doi.org/doi:10.1063/1.1430743}, author = {Gutman, Lorin and Eugene Shakhnovich} } @article {101131, title = {A structure-based method for derivation of all-atom potentials for protein folding}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {99}, year = {2002}, month = {2002}, pages = {5343 - 5348}, abstract = {A method for deriving all-atom protein folding potentials is presented and tested on a three-helix bundle protein, as well as on hairpin and helical sequences. The potentials obtained are composed of a contact term between pairs of atoms, and a local density term for each atom, mimicking solvent exposure preferences. Using this potential in an all-atom protein folding simulation, we repeatedly folded the three-helix bundle, with the lowest energy conformations having a Cα distance rms from the native structure of less than 2 {\r A}. Similar results were obtained for the hairpin and helices by using different potentials. We derived potentials for several different proteins and found a high correlation between the derived parameters, suggesting that a potential of this form eventually could be found that folds multiple, unrelated proteins at the atomic level of detail.}, url = {http://dx.doi.org/10.1073/pnas.072665799}, author = {Kussell, Edo and Shimada, Jun and Eugene I. Shakhnovich} } @article {100991, title = {Expanding protein universe and its origin from the biological Big Bang}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {99}, year = {2002}, month = {2002}, pages = {14132 - 14136}, abstract = {The bottom-up approach to understanding the evolution of organisms is by studying molecular evolution. With the large number of protein structures identified in the past decades, we have discovered peculiar patterns that nature imprints on protein structural space in the course of evolution. In particular, we have discovered that the universe of protein structures is organized hierarchically into a scale-free network. By understanding the cause of these patterns, we attempt to glance at the very origin of life.}, url = {http://dx.doi.org/10.1073/pnas.202497999}, author = {Dokholyan, Nikolay V. and Shakhnovich, Boris and Eugene I. Shakhnovich} } @article {100961, title = {Thermodynamics and Folding Kinetics Analysis of the SH3 Domain form Discrete Molecular Dynamics}, journal = {Journal of Molecular Biology}, volume = {318}, year = {2002}, month = {2002}, pages = {863 - 876}, abstract = {We perform a detailed analysis of the thermodynamics and folding kinetics of the SH3 domain fold with discrete molecular dynamic simulations. We propose a protein model that reproduces some of the experimentally observed thermodynamic and folding kinetic properties of proteins. Specifically, we use our model to study the transition state ensemble of the SH3 fold family of proteins, a set of unstable conformations that fold to the protein native state with probability 1/2. We analyze the participation of each secondary structure element formed at the transition state ensemble. We also identify the folding nucleus of the SH3 fold and test extensively its importance for folding kinetics. We predict that a set of amino acid contacts between the RT-loop and the distal hairpin are the critical folding nucleus of the SH3 fold and propose a hypothesis that explains this result.}, url = {http://dx.doi.org/10.1016/S0022-2836(02)00136-5}, author = {Borreguero, Jose M. and Dokholyan, Nikolay V. and Buldyrev, Sergey V. and Eugene I. Shakhnovich and Eugene Stanley, H.} } @article {100996, title = {Topological determinants of protein folding}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {99}, year = {2002}, month = {2002}, pages = {8637 - 8641}, abstract = {The folding of many small proteins is kinetically a two-state process that represents overcoming the major free-energy barrier. A kinetic characteristic of a conformation, its probability to descend to the native state domain in the amount of time that represents a small fraction of total folding time, has been introduced to determine to which side of the free-energy barrier a conformation belongs. However, which features make a protein conformation on the folding pathway become committed to rapidly descending to the native state has been a mystery. Using two small, well characterized proteins, CI2 and C-Src SH3, we show how topological properties of protein conformations determine their kinetic ability to fold. We use a macroscopic measure of the protein contact network topology, the average graph connectivity, by constructing graphs that are based on the geometry of protein conformations. We find that the average connectivity is higher for conformations with a high folding probability than for those with a high probability to unfold. Other macroscopic measures of protein structural and energetic properties such as radius of gyration, rms distance, solvent-accessible surface area, contact order, and potential energy fail to serve as predictors of the probability of a given conformation to fold.}, url = {http://dx.doi.org/10.1073/pnas.122076099}, author = {Dokholyan, Nikolay V. and Li, Lewyn and Ding, Feng and Eugene I. Shakhnovich} } @article {85526, title = {Combinatorial computational method gives new picomolar ligands for a known enzyme}, journal = { Proc. Natl. Acad. Sci. USA}, volume = {99}, year = {2002}, pages = {1270 - 1273}, abstract = {Combinatorial small molecule growth algorithm was used to design inhibitors for human carbonic anhydrase II. Two enantiomeric candidate molecules were predicted to bind with high potency (with R isomer binding stronger than S), but in two distinct conformations. The experiments verified that computational predictions concerning the binding affinities and the binding modes were correct for both isomers. The designed R isomer is the best-known inhibitor (Kd ~ 30 pM) of human carbonic anhydrase II.}, url = {http://dx.doi.org/10.1073/pnas.032673399}, author = {Bartosz A. Grzybowski and Alexey V. Ishchenko and Chu-Young Kim and George Topalov and Robert Chapman and David W. Christianson and George M. Whitesides and Eugene I. Shakhnovich} } @article {101431, title = {Characterization of the folding kinetics of a three-helix bundle protein via a minimalist Langevin model}, journal = {Journal of Molecular Biology}, volume = {310}, year = {2001}, month = {2001}, pages = {673 - 685}, abstract = {We use a simple off-lattice Langevin model of protein folding to characterize the folding and unfolding of a fast-folding, 46 residue three-helix bundle. Under conditions at which the C-terminal helix is 30 \% stable, we observe a clear three-state folding mechanism. In the on-pathway intermediate state, the middle and C-terminal helices are folded and in contact with each other, while the N-terminal region remains disordered. Nevertheless, under these conditions this intermediate is thermodynamically unstable relative to its unfolded state. The first and highest folding barrier corresponds to the organization of the hinge between the middle and C-terminal helices. A subsequent major barrier corresponds to the organization of the hinge between the middle and N-terminal helices. Hyperstabilizing the hinge regions leads to twice the folding rate that is obtained from hyperstabilizing the helices, even though much fewer contacts are involved in hinge hyperstabilization than in helix hyperstabilization. Unfolding follows single-exponential kinetics, even at temperatures only slightly above the folding transition temperature.}, url = {http://dx.doi.org/10.1006/jmbi.2001.4792}, author = {Berriz, Gabriel F. and Eugene I. Shakhnovich} } @article {101436, title = {Constructing, verifying, and dissecting the folding transition state of chymotrypsin inhibitor 2 with all-atom simulations}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {98}, year = {2001}, month = {2001}, pages = {13014 - 13018}, abstract = {Experimentally, protein engineering and φ-value analysis is the method of choice to characterize the structure in folding transition state ensemble (TSE) of any protein. Combining experimental φ values and computer simulations has led to a deeper understanding of how proteins fold. In this report, we construct the TSE of chymotrypsin inhibitor 2 from published φ values. Importantly, we verify, by means of multiple independent simulations, that the conformations in the TSE have a probability of ≈0.5 to reach the native state rapidly, so the TSE consists of true transition states. This finding validates the use of transition state theory underlying all φ-value analyses. Also, we present a method to dissect and study the TSE by generating conformations that have a disrupted α-helix (α-disrupted states) or disordered β-strands 3 and 4 (β-disrupted states). Surprisingly, the α-disrupted states have a stronger tendency to fold than the β-disrupted states, despite the higher φ values for the α-helix in the TSE. We give a plausible explanation for this result and discuss its implications on protein folding and design. Our study shows that, by using both experiments and computer simulations, we can gain many insights into protein folding.}, url = {http://dx.doi.org/10.1073/pnas.241378398}, author = {Li, Lewyn and Eugene I. Shakhnovich} } @article {101416, title = {Excluded volume in protein side-chain packing}, journal = {Journal of Molecular Biology}, volume = {311}, year = {2001}, month = {2001}, pages = {183 - 193}, abstract = {The excluded volume occupied by protein side-chains and the requirement of high packing density in the protein interior should severely limit the number of side-chain conformations compatible with a given native backbone. To examine the relationship between side-chain geometry and side-chain packing, we use an all-atom Monte Carlo simulation to sample the large space of side-chain conformations. We study three models of excluded volume and use umbrella sampling to effectively explore the entire space. We find that while excluded volume constraints reduce the size of conformational space by many orders of magnitude, the number of allowed conformations is still large. An average repacked conformation has 20 \% of its χ angles in a non-native state, a marked reduction from the expected 67 \% in the absence of excluded volume. Interestingly, well-packed conformations with up to 50 \% non-native χ angles exist. The repacked conformations have native packing density as measured by a standard Voronoi procedure. Entropy is distributed non-uniformly over positions, and we partially explain the observed distribution using rotamer probabilities derived from the Protein Data Bank database. In several cases, native rotamers that occur infrequently in the database are seen with high probability in our simulation, indicating that sequence-specific excluded volume interactions can stabilize rotamers that are rare for a given backbone. In spite of our finding that 65 \% of the native rotamers and 85 \% of χ1 angles can be predicted correctly on the basis of excluded volume only, 95 \% of positions can accommodate more than one rotamer in simulation. We estimate that, in order to quench the side-chain entropy observed in the presence of excluded volume interactions, other interactions (hydrophobic, polar, electrostatic) must provide an additional stabilization of at least 0.6 kT per residue in order to single out the native state.}, url = {http://dx.doi.org/10.1006/jmbi.2001.4846}, author = {Kussell, Edo and Shimada, Jun and Eugene I. Shakhnovich} } @article {101411, title = {Different circular permutations produced different folding nuclei in proteins: a computational study}, journal = {Journal of Molecular Biology}, volume = {306}, year = {2001}, month = {2001}, pages = {121 - 132}, abstract = {There have been many studies about the effect of circular permutation on the transition state/folding nucleus of proteins, with sometimes conflicting conclusions from different proteins and permutations. To clarify this important issue, we have studied two circular permutations of a lattice protein model with side-chains. Both permuted sequences have essentially the same native state as the original (wild-type) sequence. Circular permutant 1 cuts at the folding nucleus of the wild-type sequence. As a result, the permutant has a drastically different nucleus and folds more slowly than wild-type. In contrast, circular permutant 2 involves an incision at a site unstructured in the wild-type transition state, and the wild-type nucleus is largely retained in the permutant. In addition, permutant 2 displays both two-state and multi-state folding, with a native-like intermediate state occasionally populated. Neither the wild-type nor permutant 1 has a similar intermediate, and both fold in an apparently two-state manner. Surprisingly, permutant 2 folds at a rate identical with that of the wild-type. The intermediate in permutant 2 is stabilised by native and non-native interactions, and cannot be classified simply as on or off-pathway. So we advise caution in attributing experimental data to on or off-pathway intermediates. Finally, our work illuminates the results on α-spectrin SH3, chymotrypsin inhibitor 2 and β-lactoglobulin, and supports a key assumption in the experimental efforts to locate potential nucleation sites of real proteins via circular permutations.}, url = {http://dx.doi.org/10.1006/jmbi.2000.4375}, author = {Li, Lewyn and Eugene I. Shakhnovich} } @article {101391, title = {Understanding hierarchical protein evolution from first principles}, journal = {Journal of Molecular Biology}, volume = {312}, year = {2001}, month = {2001}, pages = {289 - 307}, abstract = {We propose a model that explains the hierarchical organization of proteins in fold families. The model, which is based on the evolutionary selection of proteins by their native state stability, reproduces patterns of amino acids conserved across protein families. Due to its dynamic nature, the model sheds light on the evolutionary time-scales. By studying the relaxation of the correlation function between consecutive mutations at a given position in proteins, we observe separation of the evolutionary time-scales: at short time intervals families of proteins with similar sequences and structures are formed, while at long time intervals the families of structurally similar proteins that have low sequence similarity are formed. We discuss the evolutionary implications of our model. We provide a "profile{\textquotedblright} solution to our model and find agreement between predicted patterns of conserved amino acids and those actually observed in nature.}, url = {http://dx.doi.org/10.1006/jmbi.2001.4949}, author = {Dokholyan, Nikolay V. and Eugene I. Shakhnovich} } @article {101346, title = {Mode-coupling theory for heteropolymers}, journal = {Phys. Rev. E}, volume = {63}, year = {2001}, month = {2001}, pages = {041501}, abstract = {We study the Langevin dynamics of a heteropolymer by means of a mode-coupling approximation scheme, giving rise to a set of coupled integro-differential equations relating the response and correlation functions. The analysis shows that there is a regime at low temperature characterized by out-of-equilibrium dynamics, with violation of time-translational invariance and of the fluctuation-dissipation theorem. The onset of aging dynamics at low temperatures gives insight into the nature of the slow dynamics of a disordered polymer. We also introduce a renormalization-group treatment of our mode-coupling equations, which supports our analysis, and might be applicable to other systems.}, url = {http://dx.doi.org/10.1103/PhysRevE.63.041501}, author = {Pitard, E. and Shakhnovich, E.I.} } @article {101156, title = {Energy Profile of the Space of Model Protein Sequences}, journal = {Journal of Biological Physics}, volume = {27}, year = {2001}, month = {2001}, pages = {147 - 159}, abstract = {A numerical study of the energy landscape of the space of model proteinsequences is carried out. As a consequence of the heterogeneity of thecontact energies among amino acids, the energy landscape displays a veryrough profile, a behaviour typical of frustrated systems. This givesraise to a hierarchical clustering of low-energy sequences and can have evolutionary consequences.}, url = {http://dx.doi.org/10.1023/A:1013151530254}, author = {Tiana, G. and Broglia, R. A. and Shakhnovich, E.I.} } @article {101136, title = {Exotic phase transitions in disordered globular networks}, journal = {J. Chem. Phys.}, volume = {114}, year = {2001}, month = {2001}, pages = {10968}, abstract = {We present a comprehensive numerical study of phase diagram and conformational transitions in disordered globular networks with physical/chemical crosslinks. Physically crosslinked random heteropolymer (RHP) networks at positive intersegment interaction dissimilarity, i.e., χf>0, exhibit three globular phases: frozen globular with microdomain structure, random globular, and frozen-random-globular phase, while chemically crosslinked RHPs at positive χf exhibit one random-globular phase and three frozen-globular phases; in contrast with RHPs with physical crosslinks, phase changes in chemically crosslinked RHP networks do not require reentrant passages via a random-globular phase. Based on physical grounds we explain the origins of the phase changes predicted and discuss the order of the phase transitions predicted. {\textcopyright} 2001 American Institute of Physics.}, url = {http://dx.doi.org/doi:10.1063/1.1372760}, author = {Gutman, Lorin and Eugene Shakhnovich} } @article {101081, title = {The folding thermodynamics and kinetics of crambin using an all-atom monte carlo simulation}, journal = {Journal of Molecular Biology}, volume = {308}, year = {2001}, month = {2001}, pages = {79 - 95}, abstract = {We present a novel Monte Carlo simulation of protein folding, in which all heavy atoms are represented as interacting hard spheres. This model includes all degrees of freedom relevant to folding, all side-chain and backbone torsions, and uses a Go potential. In this study, we focus on the 46 residue α/β protein crambin and two of its structural components, the helix and helix hairpin. For a wide range of temperatures, we recorded multiple folding events of these three structures from random coils to native conformations that differ by less than 1 {\r A} Cα dRMS from their crystal structure coordinates. The thermodynamics and kinetic mechanism of the helix-coil transition obtained from our simulation shows excellent agreement with currently available experimental and molecular dynamics data. Based on insights obtained from folding its smaller structural components, a possible folding mechanism for crambin is proposed. We observed that the folding occurs via a cooperative, first order-like process, and that many folding pathways to the native state exist. One particular sequence of events constitutes a "fast-folding{\textquotedblright} pathway where kinetic traps are avoided. At very low temperatures, a kinetic trap arising from the incorrect packing of side-chains was observed. These results demonstrate that folding to the native state can be observed in a reasonable amount of time on desktop computers even when an all-atom representation is used, provided the energetics sufficiently stabilize the native state.}, url = {http://dx.doi.org/10.1006/jmbi.2001.4586}, author = {Shimada, Jun and Kussell, Edo L. and Eugene I. Shakhnovich} } @article {101101, title = {Protein Folding Theory: From Lattice to All-Atom Models}, journal = {Annual Review of Biophysics and Biomolecular Structure}, volume = {30}, year = {2001}, month = {2001}, pages = {361 - 396}, abstract = {This review focuses on recent advances in understanding protein folding kinetics in the context of nucleation theory. We present basic concepts such as nucleation, folding nucleus, and transition state ensemble and then discuss recent advances and challenges in theoretical understanding of several key aspects of protein folding kinetics. We cover recent topology-based approaches as well as evolutionary studies and molecular dynamics approaches to determine protein folding nucleus and analyze other aspects of folding kinetics. Finally, we briefly discuss successful all-atom Monte-Carlo simulations of protein folding and conclude with a brief outlook for the future.}, url = {http://dx.doi.org/10.1146/annurev.biophys.30.1.361}, author = {Mirny, Leonid and Eugene Shakhnovich} } @article {100971, title = {Dynamic charge-density correlation function in weakly charged polyampholyte globules}, journal = {Phys. Rev. E}, volume = {64}, year = {2001}, month = {2001}, pages = {041802}, abstract = {We study solutions of statistically neutral polyampholyte chains containing a large fraction of neutral monomers. It is known that such solutions phase separate at very low concentrations, even if the quality of the solvent with respect to the neutral monomers is good. The precipitate is semidilute if the chains are weakly charged. This paper considers θ solvents and good solvents, and we calculate the dynamic charge density correlation function g(k,t) in the precipitate, using the quadratic approximation to the Martin-Siggia-Rose generating functional. It is convenient to express the results in terms of dimensionless space and time variables. Let ξ be the blob size, and let τ be the characteristic time scale at the blob level. Define the dimensionless wave vector q=ξk, and the dimensionless time s=t/τ. In the regime q<1, corresponding to length scales larger than the blob size, and 10.1, where entanglements are unimportant.}, url = {http://dx.doi.org/10.1103/PhysRevE.64.041802}, author = {Jan Angerman, Hindrik and Eugene Shakhnovich} } @article {87451, title = {Evolutionary conservation of the folding nucleus}, journal = {Journal of Molecular Biology}, volume = {308}, year = {2001}, month = {2001}, pages = {123 - 129}, abstract = {Here, we present statistical analysis of conservation profiles in families of homologous sequences for nine proteins whose folding nucleus was determined by protein engineering methods. We show that in all but one protein (AcP) folding nucleus residues are significantly more conserved than the rest of the protein. Two aspects of our study are especially important: (i) grouping of amino acid residues into classes according to their physical-chemical properties and (ii) proper normalization of amino acid probabilities that reflects the fact that evolutionary pressure to conserve some amino acid types may itself affect concentration of various amino acid types in protein families. Neglect of any of those two factors may make physical and biological "signals" from conservation profiles disappear.}, url = {http://dx.doi.org/10.1006/jmbi.2001.4602}, author = {Mirny, Leonid and Eugene Shakhnovich} } @article {101611, title = {Hiking in the energy landscape in sequence space: A bumpy road to good folders}, journal = {Proteins - Structure, function and bioinformatics}, volume = {39}, year = {2000}, pages = {244 - 251}, abstract = {With the help of a simple 20-letter lattice model of heteropolymers, we investigated the energy landscape in the space of designed good-folder sequences. Low-energy sequences form clusters, interconnected via neutral networks, in the space of sequences. Residues that play a key role in the foldability of the chain and in the stability of the native state are highly conserved, even among the chains belonging to different clusters. If, according to the interaction matrix, some strong attractive interactions are almost degenerate (i.e., they can be realized by more than one type of amino acid contacts), sequence clusters group into a few superclusters. Sequences belonging to different superclusters are dissimilar, displaying very small (≈10\%) similarity, and residues in key sites are, as a rule, not conserved. Similar behavior is observed in the analysis of real protein sequences. }, url = {http://dx.doi.org/10.1002/(SICI)1097-0134(20000515)39:3<244::AID-PROT70>3.0.CO;2-$\#$}, author = {Tiana, G. and Broglia, R. A. and Eugene I. Shakhnovich} } @article {101481, title = {Analysis of knowledge-based protein-ligand potentials using a self-consistent method}, journal = {Protein Science}, volume = {9}, year = {2000}, month = {2000}, pages = {765 - 775}, abstract = {We propose a self-consistent approach to analyze knowledge-based atom-atom potentials used to calculate protein-ligand binding energies. Ligands complexed to actual protein structures were first built using the SMoG growth procedure (DeWitte \& Shakhnovich, 1996) with a chosen input potential. These model protein-ligand complexes were used to construct databases from which knowledge-based protein-ligand potentials were derived. We then tested several different modifications to such potentials and evaluated their performance on their ability to reconstruct the input potential using the statistical information available from a database composed of model complexes. Our data indicate that the most significant improvement resulted from properly accounting for the following key issues when estimating the reference state: (1) the presence of significant nonenergetic effects that influence the contact frequencies and (2) the presence of correlations in contact patterns due to chemical structure. The most successful procedure was applied to derive an atom-atom potential for real protein-ligand complexes. Despite the simplicity of the model (pairwise contact potential with a single interaction distance), the derived binding free energies showed a statistically significant correlation (\~{}0.65) with experimental binding scores for a diverse set of complexes.}, url = {http://dx.doi.org/10.1110/ps.9.4.765}, author = {Shimada, J. and Ishchenko, A.v. and Shakhnovich, E.I.} } @article {101456, title = {Kinetics, thermodynamics and evolution of non-native interactions in a protein folding nucleus}, journal = {Nat Struct Mol Biol}, volume = {7}, year = {2000}, month = {2000}, pages = {336 - 342}, abstract = {A lattice model with side chains was used to investigate protein folding with computer simulations. In this model, we rigorously demonstrate the existence of a specific folding nucleus. This nucleus contains specific interactions not present in the native state that, when weakened, slow folding but do not change protein stability. Such a decoupling of folding kinetics from thermodynamics has been observed experimentally for real proteins. From our results, we conclude that specific non-native interactions in the transition state would give rise to -values that are negative or larger than unity. Furthermore, we demonstrate that residue Ile 34 in src SH3, which has been shown to be kinetically, but not thermodynamically, important, is universally conserved in proteins with the SH3 fold. This is a clear example of evolution optimizing the folding rate of a protein independent of its stability and function.}, url = {http://dx.doi.org/10.1038/74111}, author = {Li, Lewyn and Mirny, Leonid A. and Eugene I. Shakhnovich} } @article {101366, title = {Identifying the protein folding nucleus using molecular dynamics}, journal = {Journal of Molecular Biology}, volume = {296}, year = {2000}, month = {2000}, pages = {1183 - 1188}, abstract = {Molecular dynamics simulations of folding in an off-lattice protein model reveal a nucleation scenario, in which a few well-defined contacts are formed with high probability in the transition state ensemble of conformations. Their appearance determines folding cooperativity and drives the model protein into its folded conformation. Amino acid residues participating in those contacts may serve as "accelerator pedals{\textquotedblright} used by molecular evolution to control protein folding rate.}, url = {http://dx.doi.org/10.1006/jmbi.1999.3534}, author = {Dokholyan, Nikolay V. and Buldyrev, Sergey V. and Stanley, H.Eugene and Eugene I. Shakhnovich} } @article {101336, title = {Modeling the Kinetics of Acylation of Insulin using a Recursive Method for Solving the Systems of Coupled Differential Equations}, journal = {Biophysical Journal}, volume = {78}, year = {2000}, month = {2000}, pages = {652 - 661}, abstract = {This paper describes a theoretical method for solving systems of coupled differential equations that describe the kinetics of complicated reaction networks in which a molecule having multiple reaction sites reacts irreversibly with multiple equivalents of a ligand (reagent). The members of the network differ in the number of equivalents of reagent that have reacted, and in the patterns of sites of reaction. A recursive algorithm generates series, asymptotic, and average solutions describing this kinetic scheme. This method was validated by successfully simulating the experimental data for the kinetics of acylation of insulin.}, url = {http://dx.doi.org/10.1016/S0006-3495(00)76623-9}, author = {Bartosz A. Grzybowski and Anderson, Janelle R. and Colton, Ian and Brittain, Scott T. and Eugene I. Shakhnovich and George M. Whitesides} } @article {101296, title = {Orientational ordering in sequence-disordered liquid crystalline polymers}, journal = {Chemical Physics Letters}, volume = {325}, year = {2000}, month = {2000}, pages = {323 - 329}, abstract = {We study the nematic/isotropic (N/I) phase diagram and orientational ordering in heteropolymers consisting of stiff and flexible segments by field theory. It is shown that a finite and small variance in sequence disorder is sufficient to destroy orientational ordering; in the lyotropic case, the presence of sequence disorder triples the (N/I) coexistence density difference compared with that of stiff homopolymers at similar conditions. The typical N/I domain scale predicted using nucleation theory and field theory is \~{}μm in agreement with recent optical microscopy experimental results.}, url = {http://dx.doi.org/10.1016/S0009-2614(00)00693-X}, author = {Gutman, Lorin and Eugene Shakhnovich} } @article {101196, title = {Development of a Knowledge-Based Potential for Crystals of Small Organic Molecules: Calculation of Energy Surfaces for C=0{\textperiodcentered}{\textperiodcentered}{\textperiodcentered}H-N Hydrogen Bonds}, journal = {J. Phys. Chem. B}, volume = {104}, year = {2000}, month = {2000}, pages = {7293 - 7298}, abstract = {This paper describes the derivation of a Knowledge-Based Potential for intermolecular interactions from the statistical information stored in the Cambridge Structural Database. We develop a statistical mechanical method that relates the occurrences of intermolecular contacts in the database to their energies. Our approach allows us to quantify (in the form of energy) the geometrical preferences of interactions. We use our method to construct energy maps for a hydrogen bond between carbonyl oxygen and amino hydrogen. Our results demonstrate high orientational selectivity of this type of hydrogen bonding.}, url = {http://dx.doi.org/10.1021/jp000644t}, author = {Bartosz A. Grzybowski and Alexey V. Ishchenko and Robert S. DeWitte and George M. Whitesides and Eugene I. Shakhnovich} } @article {101171, title = {What can Disulfide Bonds Tell Us about Protein Energetics, Function and Folding: Simulations and Bioninformatics Analysis}, journal = {Journal of Molecular Biology}, volume = {300}, year = {2000}, month = {2000}, pages = {975 - 985}, url = {http://dx.doi.org/10.1006/jmbi.2000.3893}, author = {Abkevich, V.I. and Shakhnovich, E.I.} } @article {101031, title = {Comparison of two optimization methods to derive energy parameters for protein folding: Perceptron and Z score}, journal = {Proteins: Structure, Function, and Bioinformatics}, volume = {41}, year = {2000}, month = {2000}, pages = {192 - 201}, abstract = {Two methods were proposed recently to derive energy parameters from known native protein conformations and corresponding sets of decoys. One is based on finding, by means of a perceptron learning scheme, energy parameters such that the native conformations have lower energies than the decoys. The second method maximizes the difference between the native energy and the average energy of the decoys, measured in terms of the width of the decoys{\textquoteright} energy distribution (Z-score). Whereas the perceptron method is sensitive mainly to "outlier{\textquotedblright} (i.e., extremal) decoys, the Z-score optimization is governed by the high density regions in decoy-space. We compare the two methods by deriving contact energies for two very different sets of decoys: the first obtained for model lattice proteins and the second by threading. We find that the potentials derived by the two methods are of similar quality and fairly closely related. This finding indicates that standard, naturally occurring sets of decoys are distributed in a way that yields robust energy parameters (that are quite insensitive to the particular method used to derive them). The main practical implication of this finding is that it is not necessary to fine-tune the potential search method to the particular set of decoys used. Proteins 2000;41:192{\textendash}201. {\textcopyright} 2000 Wiley-Liss, Inc.}, url = {http://dx.doi.org/10.1002/1097-0134(20001101)41:2<192::AID-PROT40>3.0.CO;2-3}, author = {Vendruscolo, Michele and Mirny, Leonid A. and Eugene I. Shakhnovich and Domany, Eytan} } @article {101041, title = {Differential stabilization of two hydrophobic cores in the transition state of the villin 14T folding reaction}, journal = {Journal of Molecular Biology}, volume = {304}, year = {2000}, month = {2000}, pages = {99 - 115}, abstract = {We report the distribution of hydrophobic core contacts during the folding reaction transition state for villin 14T, a small 126-residue protein domain. The solution structure of villin 14T contains a central β-sheet with two flanking hydrophobic cores; transition states for this protein topology have not been previously studied. Villin 14T has no disulfide bonds or cis-proline residues in its native state; it folds reversibly, and in an apparently two-state manner under some conditions. To map the hydrophobic core contacts in the transition state, 27 point mutations were generated at positions spread throughout the two hydrophobic cores. After each point mutation, comparison of the change in folding kinetics with the equilibrium destabilization indicates whether the site of mutation is stabilized in the transition state. The results show that the folding nucleus, or the sub-region with the strongest transition state contacts, is located in one of the two hydrophobic cores (the predominantly aliphatic core). The other hydrophobic core, which is mostly aromatic, makes much weaker contacts in the transition state. This work is the first transition state mapping for a protein with multiple major hydrophobic cores in a single folding unit; the hydrophobic cores cannot be separated into individual folding subdomains. The stabilization of only one hydrophobic core in the transition state illustrates that hydrophobic core formation is not intrinsically capable of nucleating folding, but must also involve the right specific interactions or topological factors in order to be kinetically important.}, url = {http://dx.doi.org/10.1006/jmbi.2000.4190}, author = {Choe, Sung E. and Li, Lewyn and Matsudaira, Paul T. and Wagner, Gerhard and Eugene I. Shakhnovich} } @article {101036, title = {Proteins with selected sequences: A heteropolymeric study}, journal = {Phys. Rev. E}, volume = {62}, year = {2000}, month = {2000}, pages = {7100 - 7110}, abstract = {Protein sequences are expected not to be random but selected in order to form a stable native structure that is kinetically accessible. Therefore our model contains a selective temperature in sequence space (see [S. Ramanathan and E. Shakhnovich, Phys. Rev. E 50, 1303 (1994)] ) to optimize the sequence for the target conformation statistically. Replica calculations, which go beyond quadratic approximations in the field-theoretical Hamiltonian, are presented. A phase diagram indicating the temperatures and selective temperatures at which transitions to a frozen globule, i.e., the native state, occur is obtained. It is shown that going beyond the quadratic approximation in the field Hamiltonian is very important, since it results in a significant change of the phase diagram. Moreover, we suggest that a one-step replica permutation symmetry scheme is sufficient to solve the model. In addition to this we present a result for the sequence correlation function along the chain in the case of a short-ranged potential between the monomers. A correlation function between monomers that form a contact in the native state is given depending on the temperature and the interaction parameter.}, url = {http://dx.doi.org/10.1103/PhysRevE.62.7100}, author = {Wilder, J. and Shakhnovich, E.I.} } @article {100986, title = {Statistical significance of protein structure prediction by threading}, journal = {Proc. Natl. Acad. Sci. USA}, volume = {97}, year = {2000}, month = {2000}, pages = {9978 - 9983}, abstract = {In this study, we estimate the statistical significance of structure prediction by threading. We introduce a single parameter ɛ that serves as a universal measure determining the probability that the best alignment is indeed a native-like analog. Parameter ɛ takes into account both length and composition of the query sequence and the number of decoys in threading simulation. It can be computed directly from the query sequence and potential of interactions, eliminating the need for sequence reshuffling and realignment. Although our theoretical analysis is general, here we compare its predictions with the results of gapless threading. Finally we estimate the number of decoys from which the native structure can be found by existing potentials of interactions. We discuss how this analysis can be extended to determine the optimal gap penalties for any sequence-structure alignment (threading) method, thus optimizing it to maximum possible performance.}, url = {http://dx.doi.org/10.1073/pnas.160271197}, author = {Mirny, Leonid A. and Finkelstein, Alexey V. and Eugene I. Shakhnovich} } @article {85536, title = {SMoG: de Novo Design Method Based on Simple, Fast, and Accurate Free Energy Estimates. 2. Case Studies in Molecular Design}, journal = {J. Am. Chem. Soc.}, volume = {119}, year = {1997}, pages = {4608 - 4617}, abstract = {In this paper, we summarize three ligand design studies performed using the program SMoG, which was developed in our lab. The aim of this presentation is to communicate through examples the potential of this method: the richness of the molecules that can be developed and the ease with which they are found. In particular, we present suggestions for ligands to Src SH3 domain (specificity pocket and LP site) and CD4.}, url = {http://dx.doi.org/10.1021/ja963689+}, author = {Robert S. DeWitte and Alexey V. Ishchenko and Eugene I. Shakhnovich} } @article {85531, title = {SMoG: de Novo Design Method Based on Simple, Fast, and Accurate Free Energy Estimates. 1. Methodology and Supporting Evidence}, journal = {J. Am. Chem. Soc.}, volume = {118}, year = {1996}, pages = {11733 - 11744}, abstract = {In this paper, we present SMoG (Small Molecule Growth), a novel, straightforward method for de novo lead design and the evidence for its effectiveness. It is based on a simple model for ligand-protein interactions and a scoring that is directly related to the free energy through a knowledge-based potential. A large number of structures are examined by an efficient metropolis Monte Carlo molecular growth algorithm that generates molecules through the adjoining of functional groups directly in the binding region. Thus SMoG is a method that is able to rank a large number of potential compounds according to binding free energy in a short time. In this sense, SMoG represents a step toward an ideal computational tool for ligand design.}, url = {http://dx.doi.org/10.1021/ja960751u}, author = {Robert S. DeWitte and Eugene I. Shakhnovich} }