Subscriber access provided by CORNELL UNIVERSITY LIBRARY
Article
Probing Allosteric Inhibition Mechanisms of the Hsp70 Chaperone Proteins Using Molecular Dynamics Simulations and Analysis of the Residue Interaction Networks Gabrielle Stetz, and Gennady M. Verkhivker J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.5b00755 • Publication Date (Web): 22 Jul 2016 Downloaded from http://pubs.acs.org on July 27, 2016
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
Probing Allosteric Inhibition Mechanisms of the Hsp70 Chaperone Proteins Using Molecular Dynamics Simulations and Analysis of the Residue Interaction Networks Gabrielle Stetz1, Gennady M. Verkhivker1, 2
1
‡
Graduate Program in Computational and Data Sciences, Department of Computational Sciences,
Schmid College of Science and Technology, Chapman University, One University Drive, Orange, CA 92866, USA 2
Chapman University School of Pharmacy, Irvine, CA 92618, USA
‡corresponding author E-mail:
[email protected] ACS Paragon Plus Environment
1
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 114
Abstract Although molecular mechanisms of allosteric regulation in the Hsp70 chaperones have been extensively studied at both structural and functional levels, allosteric study,
the current understanding of
inhibition of chaperone activities by small molecules is still lacking. In the current using a battery of computational approaches we
probed
allosteric inhibition
mechanisms of E. coli Hsp70 (DnaK) and human Hsp70 proteins by small molecule inhibitors PET-16 and novolactone. were combined with
Molecular dynamics simulations and binding free energy analysis network-based modeling
of
residue interactions and allosteric
communications to systematically characterize and compare molecular signatures of the apo form, substrate-bound, and inhibitor-bound chaperone complexes. mechanism by which the allosteric inhibitors may
The results suggested a
leverage binding energy hotspots in the
interaction networks to stabilize a specific conformational state and impair the inter-domain allosteric control.
Using the network-based centrality analysis and community detection, we
demonstrated that substrate binding may strengthen the connectivity of local interaction communities, leading to a dense interaction network that can promote communication.
an efficient allosteric
In contrast, binding of PET-16 to DnaK may induce significant dynamic
changes and lead to a fractured interaction network and impaired allosteric communications in the DnaK complex. By using a mechanistic-based analysis of distance fluctuation maps and allosteric propensities of protein residues, we determined that the allosteric network in the PET16 complex
may be
small and localized due to the reduced communication and low
cooperativity of the substrate binding loops, which may promote the higher rates of substrate dissociation and the decreased substrate affinity. In comparison with the significant PET-16, binding of novolactone to HSPA1A may cause
effect of
only moderate network changes and
ACS Paragon Plus Environment
2
Page 3 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
preserve allosteric coupling between the allosteric pocket and the substrate binding region. The impact of novolactone on the conformational dynamics and allosteric communications in the HSPA1A complex was comparable to the substrate effect, which is consistent with the experimental evidence that PET-16, but not novolactone binding, can significantly decrease substrate affinity. We argue that the unique dynamic and network signatures of PET-16 and novolactone may be linked with the experimentally observed
functional effects
of these
inhibitors on allosteric regulation and substrate binding.
Introduction
The mammalian heat-inducible cytosolic 70-kilodalton (kDa) Heat shock protein Hsp70 along with its evolutionarily conserved bacterial orthologue, DnaK belong to a conserved and ubiquitous
chaperone family that is involved in regulation of diverse cellular processes,
including heat shock response and maintenance of protein stability.1-14 The functions of the Hsp70 chaperones in the protein folding networks and their remarkable adaptability to various substrates are regulated through an allosteric mechanism that involves cooperative involvement of the nucleotide, substrate, co-chaperones and nucleotide exchange factor proteins.15,16 functional cycle of the Hsp70 chaperones is
The
mainly driven by the nucleotide-dependent
allosteric transitions between ATP-bound and ADP-bound chaperone states. In the course of this cycle, binding of ATP in the nucleotide-binding domain (NBD) can allosterically induce the increased rates of substrate dissociation and lower the substrate binding affinity in the substratebinding domain (SBD).
Structural understanding of the Hsp70 mechanisms has been steadily
advancing in recent years,
evolving from the initial structural studies of the isolated domains
ACS Paragon Plus Environment
3
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 114
to high-resolution structures of the full-length two-domain constructs and complexes with the nucleotides and small molecule inhibitors. Crystal structures and nuclear magnetic resonance (NMR) studies of the NBD constructs17-20 and SBD complexes with the peptide substrates21-24 have provided first important insights into the structural basis of substrate specificity with Hsp70.
Amide hydrogen exchange and mass spectrometry (HX-MS) studies of the two-
domain
DnaK constructs25-28
along with
subsequent biophysical experiments29,30 have
shown that the NBD and SBD units are largely independent in a domain-undocked form of the ADP-bound and nucleotide-free states, whereas ATP binding could stabilize a domain-docked structure that facilitates the inter-domain signaling. The solution NMR structure of the fulllength DnaK has confirmed that the NBD and SBD are spatially separated and only weakly bound in the ADP-bound state.31 Crystal structures of the ATP-bound DnaK constructs32,33 have revealed the molecular details of the extensive NBD-SBD interface in the domaindocked
form.
Electron paramagnetic resonance spectroscopy34
and
single molecule
fluorescence experiments35,36 have confirmed the highly dynamic nature of both ADP-bound and ATP-bound DnaK structures, revealing significant movements of the SBD-α helical subdomain (often termed as a helical lid) and presence of multiple conformational states. Thermodynamics parameters of
functional DnaK states obtained
in isothermal titration
calorimetry (ITC) experiments37 have suggested a dynamic equilibrium between ATP-bound states and ADP-bound chaperone conformations that is accompanied by switching between low and high substrate affinity. NMR studies of allosteric regulation in DnaK have shown that
ATP binding could promote conformational flexibility of the substrate binding loops,
leading to the low affinity and fast kinetics of substrate release and association.38 This study
ACS Paragon Plus Environment
4
Page 5 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
has proposed that structural changes and allosteric communication between the substrate binding loops and the SBD-NBD interface may be central to the allosteric mechanism in DnaK.
Site-directed mutagenesis studies have
been widely employed to examine role of specific
functional residues in allosteric regulation of DnaK.39-44 A number of functionally important sites are located in the SBD regions, and N415G, P419A and D326V)
may have
mutations of these residues ( for instance K414I, significant detrimental effect on
allostery.41,44
Functional analysis of DnaK-SBD mutants has confirmed that the major structural elements responsible
for the thermodynamics and kinetics of substrate binding include the central
hydrophobic pocket, the helical lid and the “arch” region in the substrate binding site. 45,46 According to these studies, a partial removal of the helical lid by truncation of the helix B can reduce the substrate affinity by up to 10-fold, largely by adversely affecting the off-rates. Mutagenesis experiments with different lid constructs in the DnaK chaperone have
further
highlighted the role of lid helices in controlling the stability and function of the Dnak-substrate complexes.47-49
Recent experimental studies of allosteric regulation in the DnaK
have
shown that some mutant variants, including I438A, V440A, L484A, D481A, D481K, could interfere
with distinct steps of allosteric communication and produce direction-specific
communication pathways between the NBD and SBD regions.50 Collectively, these functional studies have suggested a dynamic mechanism of substrate binding and release in which the nucleotide binding signal may be transmitted through a network of allosteric
interactions
connecting the NBD-SBD interface with the substrate binding loops.
ACS Paragon Plus Environment
5
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 114
While DnaK chaperone is by far the most extensively studied member of the Hsp70 family, structural studies of other Hsp70 proteins, such as human protein Hsp70 (HSPA1A), have only recently gained a momentum due to their therapeutic importance and growing efforts to develop effective modulators of chaperone activities.
The crystal structure of the HSPA1A complex
with a peptide substrate has shown considerable similarities to the DnaK-substrate complex, yet revealing conformational differences in the loop Lα, β that bridges the SBD-α and the SBDβ subdomains, and changes in the inter-domain linker (sometimes referred as loop LL,1) .51 Deregulations of signal transduction pathways are often linked to the Hsp70-regulated processes, and selective inhibition of the Hsp70 chaperone oncogenic pathways,
could lead to the disruption of
while achieving tumor cell specificity.52,53
multiple
These revelations have
strengthened the arguments about viability of Hsp70 as an attractive drug target, fostering the development of Hsp70 inhibitors for the treatment of cancer, particularly to combat tumor cell resistance. Despite a great interest in the therapeutic potential of Hsp70 chaperones52,53 discovery of effective Hsp70 modulators has been challenging, and only a small number of inhibitors has been identified and characterized at structural and biophysical levels.54-58 The early screening studies have discovered a series of functionalized dihydropyrimidines that could stimulate the ATPase activity of Hsp70, whereas other structurally related molecules inhibit this function.59-61 High-throughput screening (HTS) assays have identified small molecules that can adjust the folding activity of the DnaK and its co-chaperones DnaJ and GrpE.62 Structure-based design inhibitors
of adenosine-derived inhibitors has led to the first structurally characterized Hsp70 that target the ATPase binding domain.63 The
identification and molecular
characterization of VER-155008 inhibitor64 in the complex with the NBD of human Hsp70 has shown that this ligand can arrest the NBD of human Hsp70 in a specific conformation and
ACS Paragon Plus Environment
6
Page 7 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
disrupt
the allosteric control.65 A small molecule, 2-phenylethynesulfonamide (PES), was
shown to bind to the C-terminal domain and inhibit the protein folding activity of Hsp70 and DnaK chaperones.65-67 Chemical biology approaches combined with structural modeling have identified a new site in the NBD of human Hsp70 chaperone68 which led to the development of several series
of allosteric
and reversible Hsp70 inhibitors.69,70
irreversible
Recently,
structural and biophysical studies have identified a new small molecule inhibitor (often referred as PET-16) that can inhibit Hsp70 and DnaK chaperones by targeting a novel C-terminal allosteric pocket.71,72
Cellular-based functional approaches have also identified and
characterized a natural product, novolactone that can affect the cellular chaperone functions through blockage of an ATP-induced substrate release and inhibition of refolding activities.73 Biochemical and structural analyses of novolactone have shown that this inhibitor is covalently bound with HSPA1A in an allosteric pocket,
which may block the inter-domain interface and
impede allosteric transitions between open and closed chaperone forms.73
High-throughput
screening of inhibitors targeting protein-protein interactions in DnaK complexes with nonenzyme partners
DnaJ and
GrpE
have identified molecules with distinct inhibitory
mechanisms.74 Computational studies have characterized conformational dynamics of the DnaK states and have offered
several atomistic models of allosteric
communications in these molecular
chaperones.75-83 Although structure and function of the Hsp70 chaperones have been actively investigated in recent years and have produced significant insights, mechanisms of allosteric regulation of chaperone activities by small molecules remained fairly elusive
due to limited
structural and dynamic information about inhibitor binding with the Hsp70 proteins.
Among
ACS Paragon Plus Environment
7
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
important questions is
to
clarify
Page 8 of 114
differences in the ligand-induced modulation of allosteric
regulation by peptide substrates and small molecule inhibitors. In the current study,
to probe allosteric inhibition of DnaK and HSPA1A chaperones, we
combined molecular dynamics (MD) simulations and binding free energy computations of the chaperone
structures with the modeling of allosteric communications and structure-based
network analysis. For network-based analysis, we employed a graph-based representation of protein structures84,85 that provides a convenient description of the residue interaction networks using topological parameters such as residue centrality.86-91 Using these approaches, we investigated and systematically compared conformational dynamics, allosteric propensities of protein residues
and organization of the residue interaction networks in
the DnaK
and
HSPA1A structures, including the apo forms and complexes with the peptide substrates and allosteric inhibitors PET-16 and novolactone (Figure 1).
To characterize
ligand-induced
changes in the allosteric interaction networks, we also determined structural organization and dynamic evolution of local stable communities in the DnaK and HSPA1A structures. show that
We
small molecule inhibitors may target central mediating sites in the residue
interaction networks and differentially modulate allosteric regulation by inhibiting specific allosteric interactions and impairing the global network connectivity. By using a mechanisticbased
analysis of distance fluctuation
maps and
residue-based modeling of allosteric
propensities, we highlight principal differences in the allosteric communication networks of the PET-16 and novolactone complexes. The results of this study unveil molecular details of allosteric inhibition in the Hsp70 proteins from both atomistic and network-centric perspectives, offering a plausible mechanism that could rationalize distinct functional effects of PET-16 and novolactone on the substrate binding and chaperone cycle.
ACS Paragon Plus Environment
8
Page 9 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
Materials and Methods Structure Preparation 500 ns MD simulations were
performed independently for the DnaK-SBD and HSPA1A-SBD
crystal structures including the unliganded structures, the complexes with the peptide substrate NRLLLTG, and the complexes with the allosteric inhibitors PET-16 and novolactone respectively (Figure 1). In addition, 5 independent 200 ns simulations were carried out for each of the studied systems to ensure computations of the system properties by averaging over multiple and
independent trajectories and also to minimize potential biases that may be
encountered in a single trajectory analysis.
We simulated
all reported apo structures of
DnaK-SBD, including the 1.4 Å crystal structure of the DnaK-SBD (residues 389-607, pub id 4F01),92 the NMR structure of the DnaK-SBD-β (residues 393-507, pub id 1DG4),24 and the three independently solved crystal structures of the unliganded DnaK-SBD (2.36 Å DnaKSBD-A, pub id 4R5J, 1.75 Å DnaK-SBD-B, pub id 4R5K, and 2.97 Å DnaK-SBD-C, pub id 4R5L).71 We report the results of MD simulations for the apo-DnaK structure (residues 381-605, pub id 4R5K). MD simulations of the DnaK-SBD complexes included DnaK-NRLLLTG X-ray structures (pub id 1DKZ, 1DKY),18
the 1.98 Å DnaK-SBD complex with the NRLLLTG
peptide ( residues 384-605, pub id 4R5I),71 and the DnaK-SBD complex with the PET-16 inhibitor (residues 378-595, pub id 4R5G).71 MD simulations were also conducted for the HSPA1A crystal structures in the unliganded form (residues 395-543, pub id 4WV5),73,the complex with the peptide substrate NRLLLTG (residues 386-613, pub id 4PO2),51 and the complex with novolactone (residues 395-543, pub id 4WV7).73 In simulations we employed several
SBD constructs with different SBD-α helical lid length,
in the presence and in the
absence of the intermolecular linker. The crystal structure of the DnaK-SBD complex with the
ACS Paragon Plus Environment
9
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
substrate (residues 384-605)
Page 10 of 114
provided coordinates for the complete SBD construct, with all
helices in the helical lid structurally resolved. In contrast, in the crystal structure of the DnaKSBD complex with PET-16 (residues 378-595), helices A and B of the helical lid do not form a continuous helix and the coordinates of the helices C, D, an E were missing. Structural changes, induced by lid truncations may have important biological consequences,
as the DnaK variant
which lacks the helical bundle (helices C-E of the lid) revealed a substantially less stable secondary structure for the helix B, causing the reduced substrate-activated ATPase rate (about two times lower than that of full-length DnaK).25 All crystal structures were obtained from the Protein Data Bank (RCSB PDB www.rcsb.org).93 Structure preparation process consisted of several stages and included checking and correcting for errors, adding missing residues and atoms, filling valences with hydrogens, predicting pK values for titratable amino acids, assigning predefined partial charges and radii to all atoms, and generating force field parameter and topology files for MD simulations. Here, we described the main components and details of structure analysis and preparation. First, protein residues in all studied crystal structures were systematically inspected for missing residues and
protons.
Hydrogen atoms and missing residues were initially added and assigned according to the WHATIF program web interface.94,95 Protonation states for titratable groups were initially determined using WHAT IF pKa calculations in aqueous solution at pH 7.0.96,97 The initial assignments of protonation states were checked for consistency and subjected to a second round of refinement using H++ approach and web server (http://biophysics.cs.vt.edu/H++).98,99 Using this approach, missing hydrogens were added to the studied protein structures based on predicted pKa of titratable groups. These calculations are based on the continuum solvent approach100
in the framework of the Poisson- Boltzmann (PB) model.97 The unresolved
ACS Paragon Plus Environment
10
Page 11 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
structural segments in the crystal structures of the DnaK complex with PET-16 (residues 596610) , the apo form of HSPA1A (residues 544-613) and the HSPA1A complex with novolactone (residues 544-613) were modeled using MODELLER program.101 modeling and a template-based loop
By combining homology
prediction approaches ModLoop102,103
ArchPRED104 we reconstructed and refined
and the
missing loops in the DnaK and HSPA1A
complexes. The chaperone structures were then optimized using the 3Drefine method that is based on an atomic-level energy minimization using a composite physics and knowledge-based force fields.105 The crystallographic coordinates of the DnaK-PET16 and HSPA1A-novolactone complexes were used to extract the inhibitor conformation, followed by optimization of the inhibitor geometry with density functional theory at the B3LYP/6-31G level using the Gaussian 09 package.106 In this stage of structure preparation, we followed the strategy of comprehensive and compatible parameterization for the protein and inhibitor atoms CHARMM General Force Field (CGenFF)-compatible parameters.107
utilizing a set of The atomistic charge
parameters for PET-16 and novolactone inhibitors were initially probed using ParamChem approach and webserver by assigning parameters based on analogy with a set of CGenFF parameters.108,109
The inhibitor parameters were subsequently generated an refined using the
restrained electrostatic
potential (RESP) charge fitting procedure110-112 on the electrostatic
potentials produced by single-point quantum mechanical calculations at the Hartree-Fock level with a 6-31G* basis set. The adopted protocol is consistent with the accepted best practices in the employment of CHARMM force field for simulations of protein-inhibitor complexes that include RESP fitting procedure or extracting charges directly from quantum mechanical calculations. A VMD plugin ffTK was used to complete computation of partial atomic
ACS Paragon Plus Environment
11
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 114
charges and parameters for small molecule inhibitors PET-16 and novolactone.113 This protocol is based on the CGenFF force field and parameterization strategy that allows for compatibility between
CHARMM22 force field parameters for proteins
topology and parameter
file information
for PET-16 and
Information). It is worth noting that a set of
and CGenFF parameters (see novolactone
in Supporting
conceptually similar tools such as
ANTECHAMBER114 and ACPYPE115 could simplify the automatic generation of topology and parameters and assignment of partial charges for small molecules from quantum mechanical computations that are compatible with the General Amber Force Field (GAFF).116 MD Simulations 500 ns MD simulations were carried out using NAMD 2.6 package117 with the CHARMM22 force field118,119 and the explicit TIP3P water
model.120
The employed MD protocol is
consistent with the overall setup that was described in details in our earlier studies.81,121,122 The initial structures were solvated in a TIP3 water box with the buffering distance of 10 Å from the protein surface to the solvation boundary using Solvate and Autoionize Plug-ins in VMD program.123 Sodium and chloride ions were added to neutralize the simulated systems and to ensure an ion concentration of 150 mM ionic strength. Solvated structures were minimized for 50,000 steps using fixed heavy atoms followed by relaxation for another 50,000 steps without applying any restraints. The following protocol was used before the production stage of MD simulations. All atoms of the complex were first restrained at their crystal structure positions with a force constant of 10 Kcal mol-1 Å-2. The system was subjected to the following simulation annealing to ensure the proper equilibration. The temperature was increased from 0K to 500K at a rate of 1K per 1ps and was kept at 500K for 500ps. The temperature was then decreased from 500 K to 300K at a rate of 1K per 1ps and was kept at 300K for additional 500ps.
ACS Paragon Plus Environment
12
Page 13 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
500 ns MD simulations were run in an NPT ensemble with a time step of 2 fs keeping the temperature at 300 K through a Langevin thermostat with a damping coefficient of 0.5 ps-1. The constant pressure was maintained at 1 atm with a Langevin-piston barostat.124 Long-range nonbonded van der Waals interactions were treated using an atom-based cutoff of 12Å with switching van der Waals potential beginning at 10Å. The smooth particle mesh Ewald (PME) method125,126 was employed to treat the long-range electrostatics. Numerical integration was performed using the leap-frog Verlet algorithm with 2fs time step. The SHAKE algorithm was employed to constrain the covalent bond lengths involving hydrogen atoms in the solute and solvent molecules.127 Principal component analysis (PCA) of the MD conformational ensembles was used to obtain the dynamic cross-correlation matrix between residues. PCA of simulation trajectories was based on the set of backbone heavy atoms (N, Cα, Cβ, C, O) to determine the essential dynamics of the chaperone structures. The calculations were performed using the CARMA program128 and PCA_NEST web service.129 Coarse-Grained Modeling of Functional Dynamics: The Gaussian Network Model Coarse-grained approaches and elastic network models (ENM) such as Gaussian network model130-132 (GNM) combined with the normal mode analysis133,134 can efficiently probe functional movements by reducing nodes where the
protein structure
to a network of uniformly connected
native interactions in the equilibrium structure determine conformational
dynamics of the system.135-137 In the GNM representation, the protein structure is modeled as a network of N nodes identified by the α-carbon atoms of proteins where the fluctuations of each node are isotropic and Gaussian. For nodes i and j , equilibrium position vectors, are R 0i and R 0j , equilibrium distance vector is R 0ij , instantaneous fluctuation vectors are ΔR i and ΔR j , and
ACS Paragon Plus Environment
13
Journal of Chemical Information and Modeling
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 114
instantaneous distance vector is R ij . The topology of the protein structure is described by N × N Kirchhoff matrix of inter-residue contacts, Γ, where the off-diagonal elements are −1 if the nodes are within a cutoff distance and zero otherwise. Bonded and nonbonded pairs of residues located within an interaction cutoff distance rc =7.0 Å are assumed to be connected by springs with a uniform spring constant γ, which is the single parameter (force constant) of the Hookean potential.130-132 The equilibrium dynamics of the structure results from the superposition of N-1 nonzero modes found by the eigenvalue decomposition of Γ. By assigning a spring constant, γ, to all contacts, the cross-correlations between the fluctuations ΔR i and ΔR j of residues (nodes) i and j are computed. The GNM normal modes are found by diagonalization of the Kirchhoff matrix, Γ = UΛUT .
ΔR i ΔR j
3kBT
[Γ-1 ]ij
3kBT
[UΛUT ]ij
3k BT
N 1
λ k 1
-1 k
[ukukT ]ij (1)
Here, U is a unitary matrix, UT U1 of the eigenvectors u k of Γ and Λ is the diagonal matrix of eigenvalues λ k .The elements of the kth eigenvector u k describe the displacements of the residues along the kth mode coordinate, and the kth eigenvalue, λ k , scales with the frequency of the kth mode, where 1 ≤ k ≤ N – 1. KB is the Boltzmann constant, T is the absolute temperature. The normal mode analysis generated eigenvalues λ k and eigenvectors u k that were obtained using the WEBnm@ approach.138 The normal mode analysis included computation of 10 low frequency modes using an approximation that was proven to be highly efficient and accurate in computations of low-frequency motions.139 The mean square fluctuations of a given residue can
ACS Paragon Plus Environment
14
Page 15 of 114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Information and Modeling
be then evaluated as a sum over the contributions of all modes. The fluctuation of the ith atomic degree of freedom along the eigenvector u k reflects the mobility of residue i in the kth mode.
(ΔR i )2
3kBT
N 1
(λ k 1
u ukT )ii (2)
-1 k k
The residue-based conformational mobility profiles in the essential space of low frequency modes were estimated using oGNM computation of structural dynamics based on the GNM approximation.132 Using this approach, we assessed cumulative contributions of the three slowest modes to the essential residue mobility and computed the mean square displacements of protein residues driven by these slow modes. Protein Stability Calculations A systematic alanine scanning of protein residues in the DnaK and HSPA1A structures was performed using FoldX approach.140,141 We utilized a graphical user interface for the FoldX calculations142 that was implemented as a plugin for the YASARA molecular graphics suite.143 If a free energy change between a mutant and the wild type (WT) proteins ΔΔG= ΔG (MT)-ΔG (WT) > 0, the mutation is destabilizing, while when ΔΔG