Reduced Point Charge Models of Proteins: Effect of ProteinâWater

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Reduced Point Charge Models of Proteins - Effect of Protein-Water Interactions in Molecular Dynamics Simulations of Ubiquitin Systems Laurence Leherte, and Daniel P. Vercauteren J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b06355 • Publication Date (Web): 02 Oct 2017 Downloaded from http://pubs.acs.org on October 4, 2017

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.

The Journal of Physical Chemistry B 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 48

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

The Journal of Physical Chemistry

Reduced Point Charge Models of Proteins - Effect of Protein-Water Interactions in Molecular Dynamics Simulations of Ubiquitin Systems

Laurence Leherte*, Daniel P. Vercauteren Laboratoire de Physico-Chimie Informatique Unité de Chimie Physique Théorique et Structurale Department of Chemistry NAmur MEdicine & Drug Innovation Center (NAMEDIC) Namur Institute of Structured Matter (NISM) University of Namur, Rue de Bruxelles 61, B-5000 Namur (Belgium)

*Corresponding author. Email: [email protected], Phone: +32-81-72.45.60, Fax: +3281-72.54.66

1

ACS Paragon Plus Environment


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 48

Abstract

We investigate the influence of various solvent models on the structural stability and proteinwater interface of three Ubiquitin complexes (PDB access codes: 1Q0W, 2MBB, 2G3Q) modeled using the Amber99sb force field (FF) and two different point charge distributions. A previously developed reduced point charge model (RPCM), wherein each amino acid residue is described by a limited number of point charges, is tested and compared to its all-atom (AA) version. The complexes are solvated in TIP4P-Ew or TIP3P type water molecules, involving either the scaling of the Lennard-Jones protein-Owater interaction parameters, or the coarse-grain (CG) SIRAH water description. The best agreements between the RPCM and AA models were obtained for structural, protein-water, and ligand-Ubiquitin properties when using the TIP4P-Ew water FF with a scaling factor γ of 0.7. At the RPCM level, a decrease in γ, or the inclusion of SIRAH particles, allows to weaken the protein-water interactions. It results in a slight collapse of the protein structure and a less compact hydration shell, and thus, in a decrease in the number of protein-water and water-water H-bonds. The dynamics of the surface protein atoms and of the water shell molecules are also slightly refrained, which allow to generate stable RPCM trajectories.

2


Page 3 of 48

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


Introduction Polar solvent molecules play an important role in the conformation and dynamics of proteins.1 They are necessary to stabilize hydrophilic residues, as well as interactions with counter-ions.2 They also play an important role on protein flexibility, and thus on protein activity.3,4 Numerous theoretical studies are completing experimental studies in allowing a description of proteinsolvent interfaces at the molecular level. For example, Caliskan et al. reported strong couplings with solvents such as glycerol and trehalose whose structure rigidity/flexibility affects the energy barriers for the conformational transitions of the protein.5 As well, Hinsen and Keller studied the influence of water on the slow internal dynamics of proteins using normal mode analysis of solvated and non-solvated lysozyme.6 They found out that few slow modes are actually affected by the solvent and that the presence of water facilitates the transitions between conformational sub-states. Experimental works achieved by Young and colleagues7,8 allowed to actually identify two types of protein motions due to a coupling with the solvent, i.e., (i) large-scale conformational changes due to dielectric fluctuations of the bulk and controlled by the solvent viscosity, and (ii) internal protein motions due to fast fluctuations in the hydration shell. Viscosity was actually reported to control both the translation and rotational diffusion of a protein, as observed by Takemura and Kitao from Molecular Dynamics (MD) simulations of Ubiquitin in pure water modeled using various force fields (FF).9 Conversely, the influence of a biomolecule on the solvent behavior has also been studied, especially through MD simulations.10 Water mobility is hampered at the protein surface,2,10 and can be clearly reduced when molecules are tightly bound to the protein10,11. By analyzing perpendicular and parallel diffusion rates to the surface, Makarov et al. proposed that the restriction of the solvent mobility at the protein-water interface could be due to (i) the reduction 3



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 48

of the local space dimensionality, (ii) the protein surface roughness, and (iii) solvent structuring by means of hydrogen bonds (H-bonds).10 The tetrahedral order parameter of hydration water, as studied by Bhattacharjee and Biswas from 188 high-resolution PDB structures,12 indeed appears to be affected for molecules very close to the protein surface, due to the non-compliance with the H-bond optimal distance. Also, Kovacs et al. proposed, from MD simulations of proteins in water, methanol, and chloroform, that the dynamical behavior of the solvent can be strongly influenced by the interactions with the protein side chains while the solvent nature has no influence on its ordering.13 Rather, through their MD simulation results, Fogarty and Laage concluded that the reorientation dynamics of water molecules is only moderately perturbed in the hydration shell of, e.g., Ubiquitin.14 However, from MD simulations as well, Schröder et al. observed that Ubiquitin and other proteins slow down the dynamics of water molecules up to distances of 1.35 nm, i.e., well beyond the first two hydration shells located in the ranges 0 to 0.35 and 0.35 to 0.60 nm from the protein surface.15 A recent review covers both theoretical and experimental aspects of water dynamics in the vicinity of biomolecules.16 Theoretical studies such as MD simulations have allowed to distinguish between different interaction types occurring between a protein and water. For example, from the analysis of radial distribution functions (RDF) calculated between Ubiquitin surface atoms and the oxygen atoms of water (Owater), Dastidar et al. clearly identified the water molecules H-bonded to acceptor groups of the surface from molecules interacting with non-hydrogen atoms of the protein,17 as earlier reported by Bizzarri and Cannistraro. 18 In the present paper, we investigate the influence of various solvent models on the proteinwater interface of Ubiquitin complexes modeled using different point charge distributions, and their consequences on the MD trajectory stability. In particular, a previously developed reduced 4


Page 5 of 48

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


point charge model (RPCM) for proteins19 is tested for solvated Ubiquitin complexes and compared to its all-atom (AA) version. It was indeed previously shown that the nature of the solvent can strongly affect the dynamics and structure of Ubiquitin described with a RPCM.20 Also, while a RPCM is able to approximate the electrostatic potential of a rigid protein,21 it affects the structure and dynamics of the system and it leads to variations in the conformations sampled through MD simulations. Nevertheless, new conformations so obtained may appear to be stable when simulated at an AA level.22 Contrarily, in the present work, we have been interested in the determination of FF modifications that generate RPCM MD trajectories similar to AA ones. That aspect is important to be evaluated when, for example, one wishes to consider large protein systems with hybrid point charge distributions, involving a combination of AA and RPCM charges. Besides the already considered TIP4P-Ew water FF, a rigid four-site model used before,19,20,22 we also apply the less time-consuming TIP3P FF, a rigid three-site model. Additionally, we have used two approaches to modulate the protein-water interactions, i.e., either a scaling of the Lennard-Jones (LJ) protein-Owater interaction parameters according to the approach of Best et al.,23 or the combination of the TIP3P water model with the coarse-grain (CG) SIRAH water description.24

Materials and Computational Methods Three Ubiquitin complexes have been considered, i.e., yeast Vps27 UIM-1 (PDB code: 1Q0W)25, polymerase (Pol) iota UBM1 (PDB code: 2MBB)26, and Ede1 UBA (PDB code: 2G3Q)27. They involve 24, 38, and 43 amino acid residues and are characterized by one, two, and three α-helix motifs, respectively (Figure 1 and Table S1).

5



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

1Q0W

2MBB

Page 6 of 48

2G3Q

Figure 1. PDB structure of Vps27 UIM-1 (PDB code: 1Q0W), Pol ι UBM1 (PDB code: 2MBB), and Ede1 UBA (PDB code: 2G3Q) ligands complexed to Ubiquitin. Ligands and Ubiquitin are presented in black and light gray, respectively. Residues Leu8, Ile44, and Val70 of Ubiquitin are in red. Residues of the ligands that are reported to interact with Leu8, Ile44, and Val70 of Ubiquitin are in green (Table S2). In the PDB structures, all ligands interact with the hydrophobic Leu8, Ile44, and Val70 Ubiquitin residues (Figure 1).25-30 A detailed analysis of the contacts was achieved using PDBsum31,32 and is presented in Table S2. In particular, the His68 residue of Ubiquitin is part of the ligand-Ubiquitin interaction network. In complexes 1Q0W and 2G3Q, His68 is protonated and is involved in H-bonds and electrostatic interactions with the ligand,25,27 while in structure 2MBB, the unprotonated His68 interacts only weakly with UBM1, as reported by Burschowsky et al.30 Following different literature sources, the ligand of structure 1Q0W presents a lower affinity than UBM1 and UBA, with dissociation constants of 277, 90, and 83 µM, respectively.25,27,30

6


Page 7 of 48

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


a. Molecular Dynamics (MD) simulations The MD protocol used to simulate the protein systems is briefly given hereafter. MD trajectories of the systems were run using the GROMACS 4.5.5 program package33,34 with the Amber99sb FF35 under particle mesh Ewald periodic boundary conditions and a Coulomb cut-off distance of 1.0 nm. The van der Waals cut-off distance was set equal to 1.4 nm. Long-range dispersion corrections to energy and pressure were applied. The initial configurations were retrieved from the Protein Data Bank36 and histidine residues occurring in 1Q0W, 2MBB, and 2G3Q were set, according to the PDB file content, to a His+, Hisδ, and His+ protonation state, respectively. Two point charge distributions were considered, i.e., either the original Amber99sb atom point charges, or the corresponding RPCM parameterized from electrostatic forces of amino acids in various conformations.19 In that last case, the point charges are located on selected atoms, except for residues His+, Phe, and Trp where some of the side chain charges must be treated as virtual sites. In a RPCM, each amino acid backbone is described by two charges, located on the C and O atoms, while the number of side chain charges differs. Apolar amino acids have no side chain charges while, e.g., tyrosine is characterized by six charges located on the Cɛ, Hδ, O(-H), and H(-O) atoms. All end residues bear an extra charge on either N (of NH3) or O (of COO-). A full description of the models can be found in the Supplementary Information File SI4 of a previous paper.19 Coordinates and topology files of the initial MD conformations are provided in Annex S1. The systems were optimized using a steepest descent algorithm with an initial step size of 0.10 nm. Then, they were solvated using TIP4P-Ew37 or TIP3P38 water molecules so that protein atoms lie at least at 1.2 nm from the cubic box walls. On the whole, from 10,500 to 12,100 water molecules were considered (Table S1). To strongly reduce the calculation time, the CG SIRAH water FF was also tested.24,39,40 In the CG cases, the 7



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 8 of 48

initial protein structures were solvated so as protein atoms lie at least at a distance dm from the cubic box walls. A shell of thickness ds of TIP3P water molecules was defined around the protein, and the remaining space of the solvation box was filled with SIRAH water beads named WT4, i.e., each water is composed of four interaction sites and represents about 11 water molecules. Three sets of dm/ds distances were considered, i.e., 2.0/1.0, 1.4/0.8, and 1.2/1.2. The first set of distances was selected from a SIRAH tutorial available on the Web.41 The second set allowed to test a computationally cheaper system, and the third set was similar to the MD conditions used to simulate Ubiquitin complexes in pure TIP4P-Ew and TIP3P solvents. The three sets of distances led to a ratio of AA over SIRAH solvent particles of about 1.4, 2, and 4.9, respectively (Table S1). For example, with dm = 2.0 and ds = 1.0, the system 1Q0W described at the RPCM level is solvated by 2266 TIP3P water molecules and 1623 SIRAH water grains. The TIP3P/SIRAH MD trajectories carried out under the three solvation conditions are labeled SIRAH2.0/1.0, SIRAH1.4/0.8, and SIRAH1.2/1.2 further in the text. Na+ ions were taken into account to neutralize the electric charge of the proteins. The whole systems were again optimized, using a steepest descent algorithm with an initial step size of 0.10 nm, to eliminate large forces and then heated to 50 K through a 10 ps canonical (NVT) MD, with a time step of 2 fs and LINCS constraints acting on bonds involving H atoms. The trajectory was followed by two successive 20 ps heating stages, at 150 K and at the final temperature of 300 K, under the same conditions. Next, each system was equilibrated during 50 ps in the NPT ensemble, at P = 1 bar, to relax the solvent molecules, and for a further 60 ns MD equilibration run. The ‘V-Rescale’ and ‘Parrinello-Rahman’ algorithms were selected to constrain T and P, respectively. A final production run of 100 ns (50 106 steps) was performed for the evaluation of

8


Page 9 of 48

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


the structural and dynamical properties of each system. Trajectory data were saved every 20 ps for systems 1Q0W and 2MBB, and every 40 ps for 2G3Q.

b. Modification of the protein-water interactions As RPCM allows to reproduce rather well the long-range electrostatic interactions,21 it was avoided to alter the Coulomb interaction terms. Rather, a scaling of the shorter range proteinwater Lennard-Jones (LJ) terms was tested as a way to modulate the RPCM impact on protein structure and dynamics in comparison with the AA model. Let us mention that revisiting Amber99sb/TIP4P-Ew interaction parameters has already been proposed before by Nerenberg et al. to refine solvation free energies.42 In the present work, the protein-Owater interactions were modified according the scheme proposed by Best et al.,23 who justify their approach by the fact that the original Amber parametrization did not explicitly involve protein-water interactions. Precisely, the Best's approach consists in a strengthening of all protein-Owater LJ terms by setting γ > 1 in: ε i −Ow = γ ε i ε Ow

(1)

where i and Ow stand for protein (Ubiquitin complex) atom i and Owater, respectively. Contrarily to the work of Nerenberg et al.,42 the σi and σOw values are left unchanged in the work of Best et al.:23 σ i −Ow =

(σ i + σ Ow ) (2) 2

In practice, the modified εi-Ow parameters are listed in the [nonbonded-params] list of the GROMACS topology files. Best et al. found that a value of γ = 1.1 allows to better reproduce the solvation effect when coupled to the Amber03/TIP3P and Amber03w/TIP4P(2005) FFs.23 For 9



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 10 of 48

example, it relaxes the usually too collapsed unfolded states, thus increasing the gyration radius, and favors the stability of non-native states. In the present paper, six values of γ, ranging from 0.6 to 1.1, were selected in combination with the water FF TIP4P-Ew to solvate the RPCM of complex 1Q0W. The intermediate value of 0.75 was also taken into account as it appeared, during the MD analysis stage, that a value between 0.7 and 0.8 might also be considered to approximate an AA MD conformational ensemble. All calculations were performed on 2.2 and 2.66 GHz Intel processors. The gain in CPU time brought by a RPCM compared to an AA one is usually very low when the protein system is buried in numerous AA solvent molecules. As expected, the use of a CG water model like SIRAH allows to significantly reduce the CPU time. As an example, the 100 ns MD trajectories of complex 1Q0W simulated at the AA level with the TIP4Ew, TIP3P, and SIRAH2.0/1.0 solvation conditions required, on 8 nodes, 198h, 190h, and 143h, respectively. The simulation of the RPCM of 1Q0W with the hybrid SIRAH2.0/1.0 conditions required 110h, i.e., 23 % less time than needed for the corresponding AA simulation under the same SIRAH2.0/1.0 conditions.

Results and discussion a. Damping of the 1Q0W-water interactions The present Section reports results and discussions regarding the analysis of MD trajectories obtained for system 1Q0W described at the RPCM level, solvated with TIP4P-Ew water molecules at different γ values. It is intended to determine a γ value that allows the RPCM of 1Q0W to remain conformationally stable and to approximate the corresponding AA MD conformational ensemble. Indeed, during the RPCM MD simulations, the ligand was previously shown to progressively adopt a bent structure, due to a region of lower α-helix propensity values, 10


Page 11 of 48

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


located at the level of residues Ala266 to Leu269.19,43 The selected value of γ will then be used for the RPCM MD simulations of the two other Ubiquitin complexes, i.e., 2MBB and 2G3Q. Structure and interface properties of complex 1Q0W obtained from the analysis of the 100 ns RPCM MD trajectories are reported in Tables 1 and SI3.

Table 1. Mean and standard deviation values of the propertiesa obtained from the analysis of the 100 ns Amber99sb AA and RPCM MD trajectories of complex 1Q0Wb using the TIP4P-Ew water FF with various scaling γ values, at 300 K and 1 bar.

γ RMSD (nm) Cplx Lig Ubq No. of Cplx surface atoms RMSF of Cplx surf. atoms (nm) No. of clusters for Lig/Ubq Rg (nm) No. of intra. H-bonds in Lig/Ubq Shell water molecules Number dmin Cplx- Owater (nm) dmin Cplx- Hwater (nm) No. of Cplx-water H-bonds No. of water-water H-bonds D of water (10-5 cm2/s) Layer 0-0.226 nm Layer 0.226-0.35 nm Layer > 0.35 nm

0.7

RPCM 0.9

1.0

AA 1.0

0.468 ± 0.031 0.562 ± 0.032 0.295 ± 0.025 377 0.239 ± 0.131 6/1 1.349 ± 0.011 6 ± 2/18 ± 3

0.679 ± 0.052 1.012 ± 0.118 0.354 ± 0.034 408 0.464 ± 0.321 16/1 1.440 ± 0.046 3 ± 2/17 ± 3

1.248 ± 0.430 1.000 ± 0.092 0.615 ± 0.108 437 1.138 ± 0.291 37/21 1.752 ± 0.235 1 ± 1/12 ± 3

0.341 ± 0.023 0.434 ± 0.037 0.234 ± 0.021 364 0.404 ± 0.168 1/1 1.374 ± 0.012 16 ± 2/53 ± 4

281 0.138 0.118 294 ± 8 431 ± 20

351 0.148 0.118 299 ± 10 541 ± 33

407 0.152 0.120 308 ± 12 656 ± 52

261 0.150 0.138 283 ± 9 436 ± 20

1.15 ± 0.20 1.41 ± 0.10 2.49 ± 0.06

1.25 ± 0.21 1.43 ± 0.10 2.49 ± 0.06

1.28 ± 0.19 1.48 ± 0.10 2.48 ± 0.05

1.21 ± 0.19 1.46 ± 0.10 2.51 ± 0.06

Ligand-Ubiquitin Tanimoto index vs. AA TIP4P-Ew 0.71 0.70 0.00 1.00 distance map Min. Lig-Ubq distance (nm) 0.179 ± 0.013 0.197 ± 0.016 0.275 ± 0.244 0.177 ± 0.008 No. of Lig-Ubq H-bonds 2.5 ± 1.5 0.7 ± 0.9 0.8 ± 1.0 4.8 ± 1.4 a A full version of the Table, including results at all tested γ values, is given as Table S3. 11



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

b

Page 12 of 48

Lig = ligand; Ubq = Ubiquitin; Cplx = Ubiquitin complex.

Particularly, the number of water molecules surrounding the Ubiquitin complex was calculated by integrating up to the first main peak of the RDF g(r), r designating the distance between Owater and the closest protein atom, named surface atom, of the Ubiquitin complex. H-bonds were determined with cut-off values of 30° and 0.35 nm for the angle Hydrogen-Donor-Acceptor and the distance Donor-Acceptor, respectively. To evaluate the self-diffusion coefficient D of water molecules located in a given shell around the protein, the MD trajectories were cut into 100 subtrajectories of 1 ns each. To limit the influence of water migration between the shell and the bulk, a mean Mean Square Displacement (MSD) value was calculated by averaging over the first 0.1 ns of each sub-trajectory. Structural stability. The stability of structure 1Q0W was monitored with (i) Root Mean Square Deviation (RMSD) calculations of the Ubiquitin complex conformations from the PDB structure, and (ii) with a clustering approach using the so-called 'gromos' method.44 For the first approach, all atoms of the structures were considered in the calculation of RMSD. For the second method, cut-off values of 0.45 and 0.3 nm were selected to determine conformation sets of the ligand and Ubiquitin, respectively. The larger cut-off of 0.45 nm was applied to the ligand which presents a more flexible structure. Very high RMSD values result from a complete deconstruction of the complex starting at γ = 0.9, with a mean value of 0.679 nm (Table 1 and Figure S1). Indeed, final conformations show that, from γ = 0.9, the ligand begins to completely unfold and, at γ = 1.0, it separates from Ubiquitin (Figure 2) after an initial bending occurring during the MD equilibration stage at the level of residues Ala266 to Leu269 (Figure S2). 12


Page 13 of 48

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


RPCM γ = 0.6 (0.496)

RPCM γ = 0.7 (0.460)

RPCM γ = 0.75 (0.508)

RPCM γ = 0.8 (0.595)

RPCM γ = 0.9 (0.761)

RPCM γ = 1.0 (1.814)

RPCM γ = 1.1 (1.979)

AA γ = 1.0 (0.354)

Figure 2. Last frames of the 100 ns Amber99sb AA and RPCM MD trajectories of complex 1Q0W in TIP4P-Ew water with various γ values, at 300 K and 1 bar. Vps27 UIM-1 and Ubiquitin are presented in black and light gray, respectively. RMSD values (nm) from the PDB structure are given in parentheses. The RMSD profiles of the whole Ubiquitin complex structure shows a sharp increase, from 0.8 to 1.08 nm at about 60 ns in the production stage, due to a complete dissociation of the complex, while the individual RMSD profiles of both complex partners are quite stable (Figure S1). Additionally, a plot of the mean square displacement (MSD) as a function of time calculated from conformations aligned with the PDB complex structure shows that the dynamics of both partners is similar to the AA one when γ = 0.6 and 0.7 for Ubiquitin and the ligand, respectively (Figure S3). Contrarily, starting from γ = 0.9, the MSD rises rapidly. 13



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 48

The RMSD values and MSD profiles are consistent with the evolution in the numbers of clusters which deviate significantly from 1 (Tables 1 and S3). In parallel, minimal ligandUbiquitin distances stay well above the PDB value of 0.168 nm, especially at γ = 1.0 and 1.1 with mean values of 0.275 and 0.485 nm, respectively. Root Mean Square Fluctuations (RMSF) of the surface atoms of the Ubiquitin complex, defined here as the atoms that are located at a distance ≤ 0.25 nm from the water atoms, arbitrarily selected at mid-trajectory, i.e., at time t = 50 ns, are also expanding to reach mean values of 0.464, 1.138, and 1.381 nm, at γ = 0.9, 1.0, and 1.1, respectively. It is interesting to note that, at γ = 0.7, their number, 377, reaches a minimum, thus approaching the AA value of 364 (Tables 1 and S3). When ligand-Ubiquitin separation occurs, the gyration radius Rg of the complex obviously increases, with a mean value of 1.752 nm at γ = 1.0, as well as the number of water molecules surrounding the protein complex, due to an enlarged solvent accessible surface. In that case, 407 molecules are reported in the first solvation layer, against, e.g., 257 at γ = 0.6, where a smaller Rg value is obtained, i.e., 1.319 nm (Table S3). The use of RPCM tends to increase Rg as well as the number of conformation clusters for Ubiquitin compared to the AA charge distribution. Thus, a γ value below 1.0 appears to cancel, at least partially, that effect. A visual inspection of the Dictionary of Secondary Structure of Proteins (DSSP) plots confirms that values of γ from 0.6 to 0.8 are better suited to depict the secondary structure elements of the protein complex (Figure S4). Indeed, the main α-helix and four over five β-strands of Ubiquitin (residues 25 to 100 in the DSSP plots) still occur versus the TIP4P-Ew AA plot. Even if the RPCM of 1Q0W appears to unfold helices more easily than it does for β-strands, α-helix motifs of the ligand (residues 1 to 24 in the DSSP plots) are slightly longer at γ = 0.7 during the 100 ns MD trajectories. Plots of the 14


Page 15 of 48

helicity degree adopted by the amino acid residues of complex 1Q0W reveal that γ = 0.7 indeed allows to partly preserve the helical content of the ligand as depicted by higher peaks from residues 1 to 24 (Figure 3). At γ = 0.6, 1.0, and 1.1, no helical character is observed. As already reported, protein RPCM-based representations systematically yield a decrease in the number of H-bonds occurring both in the Ubiquitin structure and between the two partners.20 As an example, one observes a total of 24 intramolecular H-bonds with γ values in the range 0.6 to 0.75, while the AA corresponding value is 69. The continuous decrease of intramolecular Hbonds with γ above 0.8 is due to the complete deconstruction of the protein system.

100 0.7

0.75

0.8

0.9

AA

80

% of time

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


60 40 20 0 0

20

40

60

80

100

Residue

Figure 3. Helicity profiles obtained from the analysis of the 100 ns Amber99sb RPCM MD trajectories of complex 1Q0W using the TIP4P-Ew water FF with various γ values, at 300 K and 1 bar. Residues 1 to 24 and 25 to 100 belong to the ligand and Ubiquitin, respectively. The allatom (AA) TIP4P-Ew (γ = 1.0) profile is given for comparison. With respect to an AA representation, RPCM MD results are impacted by differences in the charge distribution and by changes in the conformations of the complex. Both effects are obviously strongly linked. To separate the influence of the two effects on the unfolding of the 15



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 16 of 48

ligand, four additional MD simulations were carried out. The three first were achieved with all the atoms of the ligand, except the H atoms, kept fixed. The mobile H atoms were allowed to relax due to the motion of the protein partners and the surrounding water molecules. The Ubiquitin structure remained flexible. All conditions involved γ = 1.0 with (i) AA point charges and TIP4P-Ew applied to the PDB structure of 1Q0W, (ii) RPCM and TIP4P-Ew applied to the PDB structure of 1Q0W, and (iii) RPCM of a bent ligand and TIP4P-Ew (Figure S2). A fourth simulation was carried out on the PDB structure with rigid partners modeled at the RPCM level. Results are presented in Table 2.

Table 2. Mean and standard deviation values of the properties obtained from the analysis of the 100 ns Amber99sb AA and RPCM MD trajectories of complex 1Q0Wa using the TIP4P-Ew water FF at γ = 1.0, at 300 K and 1 bar. All atoms of the ligand, except H atoms, are kept fixed during the simulations.

Starting conformation RMSD (nm) Cplx Lig Ubq RMSF of Cplx surf. atoms (nm) No. of clusters for Lig/Ubq Rg (nm) No. of intra. H-bonds in Lig/Ubq Shell water molecules Number dmin Cplx-Owater (nm) dmin Cplx-Hwater (nm) No. of Cplx-water H-bonds No. of water-water H-bonds D of water (10-5 cm2/s)

AA

RPCM

RPCM

PDB

PDB

bent

fully rigid RPCMb PDB

0.320 ± 0.020 0.409 ± 0.003 0.266 ± 0.033 0.232 ± 0.092 1/1 1.358 ± 0.008 14 ± 2/53 ± 4

0.531 ± 0.039 0.111 ± 0.006 0.538 ± 0.033 0.207 ± 0.140 1/6 1.407 ± 0.013 11 ± 2/20 ± 3

0.781 ± 0.024 0.866 ± 0.003 0.427 ± 0.035 0.264 ± 0.109 1/6 1.441 ± 0.015 1 ± 1/17 ± 3

0.093 ± 0.001 0.115 ± 0.007 0.126 ± 0.004 0.052 ± 0.030 1/1 1.359 ± 0.001 12 ± 2/33 ± 3

273 0.150 0.140 280 ± 9 421 ± 18

335 0.150 0.106 280 ± 8 550 ± 28

342 0.150 0.112 299 ± 9 563 ± 25

290 0.148 0.124 251 ± 7 465 ± 20 16


Page 17 of 48

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


Layer 0-0.226 nm Layer 0.226-0.35 nm Layer > 0.35 nm

1.12 ± 0.16 1.41 ± 0.09 2.49 ± 0.05

1.17 ± 0.21 1.34 ± 0.09 2.47 ± 0.05

Ligand-Ubiquitin Min. Lig-Ubq distance (nm) No. of Lig-Ubq H-bonds

0.174 ± 0.006 0.180 ± 0.010 6±1 3±1 a Lig = ligand; Ubq = Ubiquitin; Cplx = Ubiquitin complex. b

1.11 ± 0.19 1.36 ± 0.08 2.47 ± 0.05

1.13 ± 0.21 1.31 ± 0.09 2.47 ± 0.05

0.192 ± 0.013 1±1

0.203 ± 0.011 0.4 ± 0.5

All atoms of the ligand and Ubiquitin, except H atoms, are kept fixed.

When the system is fully flexible, Ubiquitin adopts conformations that are gathered in 21 clusters (Table 1), while, with a rigid ligand, only 6 clusters are observed. The ligand flexibility thus favors both Ubiquitin structural changes and complex dissociation. As expected, in the fully rigid RPCM-based simulation, there is only one cluster left for Ubiquitin (Table 2). Regarding the total number of intramolecular H-bonds, i.e., within the ligand and Ubiquitin, the rigidity of the ligand at the AA level implies only two H-bonds less than in the non-rigid case, 67 against 69, respectively (Tables 1 and 2). At the RPCM level, the ligand rigidity also allows to preserve a number of intramolecular H-bonds. In both RPCM-based MD simulations carried out with a rigid PDB ligand, either with a flexible or with a rigid Ubiquitin partner, one gets 11 and 12 ligand H-bonds, respectively. Such values, which are not strongly dependent on the Ubiquitin flexibility, are close to the AA ones, i.e., 14 and 16 obtained for the rigid and flexible MD simulations (Tables 2 and 1, respectively). Contrarily, the RPCM-based MD simulation of the fully rigid system shows that there is only an average of 33 intramolecular Ubiquitin H-bonds left, while the corresponding AA value is 53 regardless of the ligand flexibility degree (Tables S1 and 2). When the ligand flexibility is considered, only 4 to 7 H-bonds are observed from γ = 0.6 to 0.8 (Tables 1 and S3).

Conformational changes thus more strongly affect the MD

trajectories than the charge distribution changes. Particularly, a drastic conformational change, 17



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 18 of 48

such as the bending of the ligand, leads to a single intramolecular H-bond left (Table 2). At γ = 1.0, the separation stage does not occur when the ligand is kept rigid, regardless of its conformation. It is consistent with the ability of RPCM to approximate longer range Coulomb interactions while a change in the point charge distribution considerably affects the intramolecular electrostatic interactions so as to lead to unfolding. Complex-water interface. In the present Section, structural and dynamical aspects of the complex-water interactions are analyzed in terms of H-bonds and RDFs. The mobility of the surface atoms of the Ubiquitin complex and of the water molecules in the vicinity of the complex is also investigated. The number of complex-water H-bonds in the bound system is relatively stable, with mean values ranging between 289 at γ = 0.75 and 301 at γ = 0.6, which are nevertheless above the AA value of 283 (Tables 1 and S3). The increase observed between γ = 0.8 and 1.1, from 291 to 321 complex-water H-bonds, is due to the separation of both partners, which leads to a larger solvent accessible surface, and where ligand-Ubiquitin H-bonds are progressively replaced by complex-water H-bonds. As already stated, while the number of complex-water H-bonds is clearly linked to the protein complex conformations, it does not appear to depend on the point charge model, with 280 H-bonds present both at the AA and RPCM levels for a rigid ligand (Table 2). The RPCM complex with the bent ligand is characterized by additional complex-water H-bonds compensating less numerous intramolecular and ligand-Ubiquitin H-bonds, with values of 299, 1, and 17, respectively. Finally, one obtains a better approximation of the number of water-water H-bonds occurring in a 0.35 nm shell surrounding the complex at γ = 0.7, with a value of 431 for 281 molecules to be compared to 436 for 261 molecules in the AA case (Table 1).

18


Page 19 of 48

As observed previously,20 only one hydration peak is visible in the RDFs (Figure 4, left), contrarily to AA results which show that the first hydration shell actually involves two types of water molecules, visualized as a small peak around 0.2 nm and a larger one at about 0.3 nm.11

6000

6000 0.6

0.7

0.75

0.8

0.9

1.0

1.1

AA

4000 3000 2000 1000

5000

g(complex-Hw)

5000

g(complex-Ow)

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


4000 3000 2000 1000

0

0 0

0.2

0.4 d (nm)

0.6

0

0.2

0.4 d (nm)

0.6

Figure 4. RDF of the Ubiquitin complex surface atoms - Owater and Hwater obtained from the analysis of the 100 ns Amber99sb RPCM MD trajectories of complex 1Q0W using the TIP4PEw water FF with various γ values, at 300 K and 1 bar. The all-atom (AA) TIP4P-Ew (γ = 1.0) profile is given for comparison. Nevertheless, plots of the RDF of the distances between the Ubiquitin complex surface atoms and Hwater let clearly appear a peak at distances below 0.2 nm, as reported in literature,11,45 (Figure 4, right). Such differences in the RDFs between RPCM and AA charge distributions are due to a change in the mean orientation of the water molecules facing the surface of the Ubiquitin complex. RDFs present a slight increase in the complex-Owater shortest distance dmin with γ, i.e., from 0.138 to 0.148 nm for the bound complex structures (Tables 1 and S3), even if steric LJ parameters were left unmodified (Equation 2), while no precise trend is seen in the

19



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 20 of 48

complex-Hwater RDFs. On the whole, dmin values always stay below the corresponding AA distance of 0.150 and 0.138 nm for complex-Owater and complex-Hwater, respectively. The influence of γ on the diffusion of the solvent molecules closest to the Ubiquitin complex surface, i.e., in a layer of 0.226 nm from the surface atoms, is rather weak. It nevertheless increases monotonically as γ increases, with D going from 1.04 to 1.42 10-5 cm2/s (Tables 1 and S3). Between 0.226 and 0.35 nm from the surface, such a trend is less marked, except for the dissociated complex at γ = 1.0 and 1.1. Beyond 0.35 nm, D is rather constant and close to the reported bulk value in literature, i.e., 2.4 10-5 cm2/s.37 Even if the dynamics of proximal water molecules is not significantly affected by γ, there is a correlation between their self-diffusion coefficient and the mobility of the complex surface atoms, the latter being quantified through their RMSF. Both D in layer 0 to 0.226 nm and RMSD increase with γ, i.e., when the complexOwater LJ interactions become stronger. However, changes in D stay within their standard deviation values, and RMSF is likely to be more sensitive to ligand-Ubiquitin separation than to γ, with values well above 1.0 nm starting at γ = 1.0 (Table 1 and S3). It is interesting to note that, in the particular case of γ = 0.9, the RMSF profile corresponding to the 124 first complex surface atoms, which belong to the ligand, show extremely large values at the peptide ends (Figure 5). They are thus the most sensitive ones to a change in the complex-Owater LJ interactions. The complete separation stage is obtained at larger γ values (Figure 2).

20


Page 21 of 48

0.6

0.7

0.75

0.8

0.9

1

1.1

2

RMSF (nm)

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


1

0 0

50

100

150

200 250 300 Atom number

350

400

450

Figure 5. RMSF of the Ubiquitin complex surface atoms obtained from the analysis of the 100 ns Amber99sb RPCM MD trajectories of complex 1Q0W using the TIP4P-Ew water FF with various γ values, at 300 K and 1 bar. The number of water molecules in a ~0.35 nm layer surrounding the protein complex is drastically increased at the RPCM level compared to the AA level. Even for the RPCM MD simulations involving rigid ligands, one observes 335 and 342 molecules with the PDB and bent ligand, respectively, while one gets only 273 water molecules from the AA trajectory (Table 2). For the fully rigid RPCM case, one also detects a relatively low number of surrounding water molecules, only 290, for a Rg value of 1.359 nm that is close to the AA value of 1.358 nm. The number of surrounding molecules is thus related to the size of the complex, but not to the minimal ligand-Ubiquitin distance which is larger for the fully rigid RPCM complex than it is for the AA complex, with values of 0.203 and 0.174 nm, respectively (Table 2). It is again noted that the minimal distance observed between the complex surface atoms and Hwater is reduced for the RPCMs, with a value of 0.124 nm for the fully rigid complex, and even lower values of 0.106 21



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 22 of 48

and 0.112 nm (Table 2), compared to the AA TIP4P-Ew case with a value of 0.138 nm (Table 1). However, complex-Owater distances, with a constant value of 0.150 nm (Table 2), are not affected by the use of a RPCM when the ligand is kept rigid. Slightly lower values, between 0.138 and 0.146 nm, are obtained when the ligand is flexible (Tables 1 and S3). Ligand-Ubiquitin contacts. The rigid bent ligand is located at a minimal distance from the receptor that is slightly longer than the one observed for the PDB structure, i.e., 0.192 and 0.167 nm, respectively. Less H-bonds are detected, i.e., 1 rather than 3 for the rigid ligand in its PDB conformation (Table 2). The number of ligand-Ubiquitin H-bonds is also inversely related to the ligand-Ubiquitin distances, with values of 3, 1, and 0.4, for distances of 0.180, 0.192, and 0.203 nm, respectively (Table 2). Distance maps between the atoms of Vps27 UIM-1 and Ubiquitin were built by considering the minimal distance observed during a MD simulation between the atoms of each ligand residue and the atoms of each Ubiquitin residue. Values are averages over the whole production trajectories (Figures 6 and S5). In the map generated by the analysis of the AA MD trajectory simulated in TIP4P-Ew water molecules, one clearly distinguishes three regions extended along the Vps27 UIM-1 chain. This extension is due to the spatial alignment of Vps27 UIM-1 with a number of β-strands of Ubiquitin. The first region corresponds to the contacts occurring between segment 259 to 272 of Vps27 UIM-1 and β-strand residues 4 to 10 of Ubiquitin, while the second and third regions are due to contacts with β-strands residues 40 to 45 and 48 to 49, and 66 to 72, respectively.

22


Page 23 of 48

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


1Q0W

(1.00)

(0.71)

(0.60)

(1.00)

(0.86)

(0.69)

(1.00)

(0.62)

(0.81)

AA TIP4P-Ew (γ = 1.0)

RPCM TIP4P-Ew (γ = 0.7)

RPCM SIRAH2.0/1.0

2MBB

2G3Q

Figure 6. Ligand-Ubiquitin minimal distance maps obtained from the analysis of the 100 ns Amber99sb RPCM MD trajectories of Ubiquitin complexes under various solvation conditions, at 300 K and 1 bar. The gray scale goes from 0 to 1.5 nm (black to white) using a distance increment of 0.30 nm. Tanimoto indices calculated from the map of the AA complex solvated in TIP4P-Ew water are reported in parentheses. 23



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 24 of 48

The shortest ligand-Ubiquitin distances, below 0.3 nm, occur in the two last areas. All trajectory analyses, but the ones at γ ≥ 1.0, show the three characteristic regions mentioned above. A similarity score between a reference map, i.e., the TIP4P-Ew AA distance map, and all other RPCM maps, was calculated using the well-known Tanimoto index applied to each of the N pixels of the maps. Pixels corresponding to distances between 0 and 0.15, 0.15 and 0.30 nm, ... receive a value of 10, 9, ..., respectively. All pixels corresponding to distances beyond the truncation value of 1.50 nm have a value of 0. A total of 11 distance ranges, or bins, are thus selected. The similarity measure between two maps is calculated by summing over the product of their pixel values p and q: N

Sij = ∑ pi qi

(3)

i =1

The Tanimoto index T is obtained from the following relationship:

T=

S aa

S ab + S bb − S ab

(4)

where a and b differentiate the two maps, respectively. The highest indices for the RPCM trajectories are obtained at γ = 0.7 and 0.9, with values of 0.71 and 0.70, respectively (Tables 1 and S3 and Figures 6 and S5). However, at γ = 0.9, both ligand extremities form fewer contacts with Ubiquitin than at γ = 0.7. It is consistent with the increased atom fluctuations at the same γ value, as illustrated in Figure 5. The Tanimoto indices appear to be rather insensitive to the number of bins dividing the truncation distance. With a larger number of bins, 40 rather than 11, the index differences vary from -0.01 to +0.04. Regarding the ligand-Ubiquitin H-bonds, analyzed using VMD,46 the AA number of 4.8 is never reached with the RPCM representations (Tables 1 and S3). For RPCM, the maximum 24


Page 25 of 48

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


number is 3 at γ = 0.75. At γ = 1.0, both ligand ends interact with Ubiquitin through H-bonds, while at γ = 1.1, the still occurring H-bonds involve electrically charged residues at one of the ligand ends, i.e., Asp258, Glu259, and Glu261, which interact with Lys6 and Lys11 of Ubiquitin (Table S4). In the complex, the H-bond between Glu273 and Ser65 is characterized by an extremely high occurrence frequency at γ = 0.75, i.e., 81.6 %. The AA simulation shows that His+68 of Ubiquitin is actually preferentially H-bonded to Glu273 of the ligand, as also observed in RPCM trajectories at γ = 0.75 and 0.9. The frequency is almost as high at TIP4P-Ew γ = 0.75 as it is at the AA level, with values of 47.4 and 49.8 %, respectively. At γ = 0.7, His+68 is Hbonded to Ser270, a contact also reported in the PDBsum results (Table S2). Interactions involving Arg42 and Arg72 of Ubiquitin appear as H-bond types only at the AA level, in pure TIP4P-Ew (Table S4) and TIP3P solvent (Table S5). The two water FFs involve the Arg residues to be H-bonded to glutamate residues of the ligand, more precisely to Glu260 and Glu259 in TIP4P-Ew water, and to Glu259 and Glu257 in TIP3P water, respectively. Finally, numerous and highly frequent H-bonds are observed at the AA TIP3P level, all involving glutamate residues of the ligand as well as Ser270 and Asp258 (Table S5). It is due to shorter ligand-Ubiquitin distances compared to the TIP4P-Ew solvation conditions, with values of 0.170 and 0.177 nm, respectively.

All the observations reported above allow to select a γ value of 0.7 as a possible optimal parameter to simulate RPCM protein systems in TIP4P-Ew water. The two major points that led to the choice of γ = 0.7 rather than γ = 0.75 are (i) the minimal distance map that is better approximated at γ = 0.7 than it is at γ = 0.75, with Tanimoto coefficients of 0.71 and 0.65, respectively, and (ii) the DSSP plots (Figure S4) which show that secondary structure elements 25



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 26 of 48

of Ubiquitin particularly, have a longer life time. Additionally, the mean RMSD values of the ligand and Ubiquitin at γ = 0.7, i.e., 0.562 and 0.295 nm, are slightly below the corresponding values obtained at γ = 0.75, i.e., 0.647 and 0.345 nm, respectively (Table S3), while the overall RMSD value is only slightly larger. With respect to the unscaled LJ protein-Owater energy terms (γ = 1.0), a decrease in γ weakens the Ubiquitin-water and ligand-water interactions. It induces a slight collapse of the protein complex structure, which is consistent with the conclusions obtained by Best et al.23 A weakening of complex-water interactions allows the water molecules of the solvation shell to be less packed around the protein system, as seen through the decrease in the number of such molecules. As a consequence, the number of H-bonds they form between themselves decreases too. The dynamics of the complex surface atoms is slightly refrained when γ decreases, while the dynamics of water molecules is only very slightly affected by γ and by the dynamics of the protein complex. D is reduced for the fully rigid RPCM of 1Q0W in comparison to the RPCM where only the ligand is fixed (Table 2). Conversely, D of the water molecules around the AA model is larger than the corresponding values for the AA where only the ligand is fixed (Tables 1, S3, and 2). The motion of the protein atoms thus affects, at least partly, water diffusion.

b. Hybrid all-atom/coarse-grained water description level applied to system 1Q0W The present Section reports results and discussions based on the analysis of MD trajectories obtained for the test system 1Q0W at the RPCM level, solvated with a combination of TIP3P and SIRAH water models. As described above, three box size and water shell conditions were tested on the RPCM description of system 1Q0W to determine the factors that favor stable conformational states. Results of the MD simulations compared to pure TIP3P MD are also 26


Page 27 of 48

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


reported for comparison in Table S6. A gain of about 40 % in the simulation time is observed between the most time-consuming TIP3P/SIRAH and the pure TIP3P MD simulations of the RPCM of complex 1Q0W, i.e., 110 h and 183 h, respectively (Table S6). Structural stability. Only the SIRAH2.0/1.0 solvation conditions retain some secondary structure elements of complex 1Q0W, especially four of the five β-strands of Ubiquitin. The αhelix component of Ubiquitin is almost totally absent compared to the pure AA TIP4P-Ew FF at γ = 0.7 (Figure S4). The RMSD of the whole complex is the lowest, with a mean value of 0.604 nm (Table S6), which is nevertheless larger than the corresponding TIP4P-Ew values generated with γ ≤ 0.8, i.e., from to 0.452 to 0.551 nm (Tables 1 and S3). However, the RMSD profiles obtained with the SIRAH2.0/1.0 conditions are not as stabilized as with TIP4P-Ew and suggest that longer simulations should be carried out (Figure S1). The apparent convergence of the SIRAH1.2/1.2 RMSD profile does not occur for the two other Ubiquitin systems described later. A consequence is that the gain in calculation time due to coarse-graining the solvent will be partly compensated by longer simulation times. Thus, among the three simulations conditions, the procedure described in reference41 can be better suited to simulate the structural properties of the RPCM of 1Q0W, but is however slow to converge. General conclusions resulting from the analysis of the SIRAH-based MD trajectories are given below. Additional details are reported in the Supplementary Information (Annex S2). Complex-water interface. The diffusion coefficient of the water molecules in a close vicinity to the RPCM of 1Q0W is reduced compared to the pure TIP3P case. For example, a value of 1.98 10-5 cm2/s is obtained with SIRAH2.0/1.0 in comparison with 3.4 10-5 cm2/s with the pure TIP3P solvent.

27



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 28 of 48

As already mentioned, SIRAH beads are initially placed at a distance of 1.0 nm from the structure of the complex. During the MD simulation, the CG beads diffuse and are finally distributed at distances beyond 0.21 nm from the complex surface (Figure S6). The RDFs g(surface-WT4) reveal a small peak at very short distances, centered at 0.35 nm, followed by a minimum at about 0.45 nm (Figure S6). It is due to a very limited number of closest SIRAH beads (Figure S7 top). Below 0.21 nm, only water molecules surround the protein complex, while a mixture of WT4 beads and water molecules are present at distances beyond 0.45 nm from the surface (Figure S7 bottom). Thus, even if the closest water molecules from the surface do not mix with SIRAH beads, their self-diffusion coefficient is reduced in comparison with the pure TIP3P solvation conditions. It is assumed that long-range SIRAH-TIP3P interactions are strong enough to affect the mobility of the water molecules, even those in direct interaction with the protein structure. From Figure S7 bottom, it is clear that the distribution of WT4 sites is not homogeneous around the protein system. A visual analysis shows that WT4 beads preferentially avoid large regions of space where the electrostatic potential is mainly negative, as also observed for the other Ubiquitin complexes (Figure S8). Ligand-Ubiquitin contacts. The best agreement between the RPCM TIP3P/SIRAH and the AA TIP4P-Ew minimal distance maps is observed with SIRAH2.0/1.0 solvation conditions, with a Tanimoto similarity coefficient of only 0.60 (Table S6 and Figures 6 and S5) compared to the values obtained with pure TIP4P-Ew water at γ = 0.7. Also, as seen before with TIP4P-Ew, the use of a RPCM leads to a decrease in the number of ligand-Ubiquitin H-bonds. Only 2 to 4 Hbonds are detected while 11 are observed at the AA level (Table S6).

28


Page 29 of 48

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


Adding SIRAH particles to a TIP3P water FF has thus an effect that is similar to a reduction of γ in Amber99sb-TIP4P-Ew LJ terms, i.e., it slows down the dynamics of the Ubiquitin complex atoms and reduces the gyration radius of the complex, thus preventing the ligand and Ubiquitin from separation. Contrarily to the effect of scaling Amber99sb-TIP4P-Ew LJ interactions with γ = 0.7, SIRAH significantly affects the diffusion of water molecules compared to pure TIP3P. Indeed, D of AA water molecules is reduced by 5 % (1.15 10-5 versus 1.21 10-5 cm2/s) (Tables 1 and S3), and by 42 % (1.98 10-5 versus 3.43 10-5 cm2/s) (Table S6), respectively.

c. MD simulations of the 1Q0W, 2MBB, and 2G3Q Ubiquitin complexes In this Section, we analyze and discuss the AA and RPCM MD trajectories of two additional Ubiquitin complexes named 2MBB and 2G3Q. RPCMs were simulated under the best solvation conditions discussed above for complex 1Q0W, i.e., either with pure TIP4P-Ew water molecules at γ = 0.7, or with the hybrid SIRAH2.0/1.0 water FF. Structural stability. Final AA MD conformations of the complexes are close to their corresponding PDB starting structure (Figure 1), as shown in Figure S9. Their RMSD values calculated from the corresponding PDB structure, as well as the mean RMSD calculated from the corresponding 100 ns MD trajectory, are indeed lower than 0.4 nm. RMSD profiles are given in Figure S1. As for complex 1Q0W simulated at the AA level, the ligand of 2MBB is rather flexible. That behavior is hampered for the ligand of complex 2G3Q, especially with the pure TIP3P water FF. At the RPCM level, the remark previously raised for system 1Q0W is also valid for 2MBB and 2G3Q, i.e., the SIRAH RMSD profiles present a drift towards larger values. The analysis of the AA MD trajectories shows that Ubiquitin is actually characterized by single-cluster conformations, regardless of the FF conditions, while the ligands adopt 1 to 3 main 29



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 30 of 48

conformations due to less numerous intramolecular interactions (Table S7). TIP4P-Ew is the water FF that restricts the most the ligand flexibility, with 1 cluster only for the three Ubiquitin complexes (Table S7). Mean RMSD values are rather low, with values lower than 0.3 nm for the largest complex 2G3Q. At the RPCM level, Vps27 UIM-1, whose helix can be easily deformed, is not the only ligand to separate from its receptor. Indeed, a scaling factor γ = 1.1 applied to the TIP3P water FF reinforces the separation of the ligand of structures 2MBB and 2G3Q from Ubiquitin (Figure 7). Huge numbers of ligand clusters are thus obtained with TIP3P at γ = 1.1, i.e., 838 and 500, respectively (Table S8). With TIP4P-Ew, the ligand Vps27 UIM-1 of structure 1Q0W has the particularity to bend and then to separate from Ubiquitin, while the two other complex structures are stable. It is correlated to the larger dissociation constant reported for complex 1Q0W.25 Let us mention that we detected only the bending stage in our previous work because of shorter simulation times.19 Contrarily to TIP3P at γ = 1.1, both Ubiquitin and ligand conformational changes are clearly restricted with the solvation model TIP4P-Ew at γ = 0.7. It is confirmed by the lowest mean RMSD values, i.e., 0.468, 0.438, and 0.480 nm, for structures 1Q0W, 2MBB, and 2G3Q, respectively (Table S8). At the AA level, the mean gyration radius Rg values of structures 1Q0W, 2MBB, and 2G3Q, are all rather close to the PDB values of 1.347, 1.305, and 1.456 nm, respectively (Table S7). At the RPCM level, Rg can be larger (Table S8). For example, for 1Q0W in TIP3P water, one obtains a gyration radius Rg = 1.415 nm and 1.316 nm at the RPCM and AA levels, respectively (Tables S7 and S8).

30


Page 31 of 48

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


1Q0W

(1.814)

(0.460)

(1.326)

(1.599)

(0.707)

(0.684)

(0.482)

(0.677)

(1.338)

(0.757)

(0.746)

(0.476)

(0.657)

(2.898)

(0.701)

TIP4P-Ew (γ = 1.0)

TIP4P-Ew (γ = 0.7)

TIP3P (γ = 1.0)

TIP3P (γ = 1.1)

SIRAH2.0/1.0

2MBB

2G3Q

Figure 7. Last frame of the 100 ns Amber99sb RPCM MD trajectories of the Ubiquitin complexes simulated under various solvation conditions, at 300 K and 1 bar. Ligands and Ubiquitin are presented in black and light gray, respectively. RMSD values (nm) from the PDB structures are given in parentheses. The Rg of the AA Ubiquitin complexes solvated in TIP3P water is slightly affected by γ. Mean values reported in Table S7 show, as also reported by Best et al.,23 a slight increase of Rg with γ, which however stay below the corresponding standard deviation values for structures 2MBB and 2G3Q.

31



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 32 of 48

On the whole, the best solvation condition that lead to stable RPCM complexes is TIP4P-Ew at γ = 0.7 that reduces the influence of water on the structure and dynamics of the Ubiquitin complexes. At the AA level, the RMSF of the surface atoms of the Ubiquitin complexes is not significantly modified when changing the water FF. All values stay below averages of 0.3 nm, except for complex 1Q0W with mean values that are consistently higher than for 2MBB and 2G3Q, i.e., between 0.345 and 0.583 nm as compared to between 0.199 and 0.282 nm (Table S7). Contrarily, the RPCM MD simulations provide mean RMSF values that differ with the solvation model, with the highest values, above 1 nm, obtained when unfolding is observed, i.e., with TIP3P at γ = 1.1. The lowest values are obtained with TIP4P-Ew at γ = 0.7 and SIRAH2.0/1.0. As an example, for 1Q0W, values of 0.239 and 0.333 nm are respectively obtained, to be compared to 1.138 nm for TIP4P-Ew with γ = 1.0, 0.764 nm for TIP3P with γ = 1.0, and 1.312 nm for TIP3P with γ = 1.1 (Table S8). Despite the difficulty of MD simulations to stabilize helices in RPCMs, the DSSP analysis of TIP4P-Ew MD trajectories carried out with γ = 0.7 shows that helical segments of the ligands are preserved, especially for 2G3Q (Figure 8). Plots displayed in Figure 8 were obtained as time averages over the corresponding DSSP data (Figure S10).

32


1

1

0.8

0.8

0.6

0.6

0.6

0.4 0.2

0.2

0.2 0

0 0

10

30

40 50 60 Residue number

70

80

90

0

100

10 20 30 40 50 60 70 80 90 100 Residue number

0 0

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0 0

2G3Q

20

Secondary structure occurrence

2MBB

0.4

0.4

10

20

30

40 50 60 70 80 Residue number

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

10

20

30

40


80

90

100 110

0

10

20

30

40 50 60 70 80 Residue number

90 100 110

AA TIP4P-Ew (γ = 1.0)

20

30


70


80

80

90

100

0

10

20

30

40

90

100 110

1

0

0

10

0 0

90 100 110


1Q0W

1 0.8


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 33 of 48

0.8

0.6

0.4

0.2

0

0

10

20

30

40

50 60 70 80 Residue number

90 100 110 120

0

RPCM TIP4P-Ew (γ = 0.7)

10 20 30 40 50 60 70 80 90 100 110 120 Residue number

RPCM SIRAH2.0/1.0

Figure 8. Occurrence frequency of DSSP secondary structure elements obtained from the analysis of the 100 ns Amber99sb AA and RPCM MD trajectories of Ubiquitin complexes using various solvent conditions, at 300 K and 1 bar. Secondary structure elements are color-coded as: α-helix (blue), β-sheet (red), bend (green), turn (yellow). Residues 1 to 24, 1 to 38, and 1 to 43 belong to the ligand of structures 1Q0W, 2MBB, and 2G3Q, respectively. The size of the Ede1 UBA ligand helps in stabilizing its structure through numerous ligand intramolecular interactions, as seen in the distance maps (Figure S11), specific ligand-ligand intramolecular contacts appearing in the bottom left area of the maps. As an example, for complex 2G3Q, a total of 55 to 85 pairs of residues (excluding up to i-i+4 pairs, i.e., pairs of Hbonded residues in α-helices) separated by distances ≤ 0.60 nm are detected, while one gets 42 to 61 for 2MBB and 6 to 17 for 1Q0W (Table S9).

33



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 34 of 48

Complex-water interface. The observations reported for complex 1Q0W in previous Sections b. and c. are also valid for structures 2MBB and 2G3Q. In parallel to the slight increase of Rg with γ in the AA TIP3P MD trajectories, one observes an increase in the number of first shell solvent molecules located within a distance of ~0.35 nm from the surface atoms of the Ubiquitin complexes. For example, the number of TIP3P water molecules surrounding complex 1Q0W increases from 247 at γ = 1.0 to 258 molecules at γ = 1.1 (Table S7). Thus, one reaches values that are closer to the corresponding ones obtained under the TIP4P-Ew solvation, i.e., 261 for 1Q0W. For all solvation FFs, the RPCMs involve more compact first water shells, with larger numbers of molecules. For example, with TIP4P-Ew, one observes 396, 375, and 381 water molecules in the first hydration shell of the three RPCM complexes (Table S8), rather than 261, 294, and 288 molecules at the AA level (Table S7). Under TIP4P-Ew at γ = 0.7 and SIRAH2.0/1.0 solvation conditions, the number of solvent molecules surrounding, e.g., complex 1Q0W, is drastically reduced, with values of 265 and 292 molecules, respectively (Table S8). Those two solvation conditions thus allow to recover the hydration shell density of the AA calculations, as also observed for structures 2MBB and 2G3Q. Consistently with more compact first hydration shells, the use of RPCM increases the number of complex-water and water-water H-bonds in comparison with the AA cases. It also prevents a clear distinction of the two solvation peaks within ~0.35 nm of the protein surface,11 as already discussed above. Such a distinction can however be seen at the level of RDFs g(complex-Hwater), as illustrated for complex 1Q0W in Figure 4. Experimental and MD studies have established that water motions within the first hydration shell are 3 to 5 times slower than in bulk.2,10 Indeed, for AA TIP4P-Ew, the self-diffusion coefficient D values of 1.50, 1.50, and 1.42 10-5 cm2/s are obtained for complex 1Q0W, 2MBB, 34


Page 35 of 48

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


and 2G3Q, respectively (Table S7), while D of bulk water is equal to 2.4 10-5 10-5 cm2/s.37 Such a behavior was also reported by Abseher et al. for solvated SPC proteins,11 and Dastidar and Mukhopadhyay for solvated Ubiquitin with the TIP3P FF, whose bulk D values are 3.85 and 5.19 10-5 cm2/s, respectively.17 Increasing γ when using the TIP3P FF does not significantly affect D at the AA level (Table S7). Differences between D values are at most equal to 0.08 10-5 cm2/s, as observed for complex 2MBB in the shell ranging from 0 to 0.226 nm. Conversely, for RPCMs solvated in TIP3P or TIP4P-Ew water, D in layers 0 to 0.226 and 0.226 to 0.35 nm systematically increases with γ, while the opposite phenomenon occurs for water molecules beyond 0.35 nm (Table S8). However, variations often stay within the standard deviation values associated with D. It is interesting to note that D calculated from RPCM trajectories are larger than the corresponding values at the AA level for a same FF. The differences almost vanish for molecules located beyond 0.35 nm. For example, the water molecules surrounding the AA model of 1Q0W at distances beyond 0.35 nm from the protein surface are characterized by D values of 2.51, 5.67, and 5.65 10-5 cm2/s (Table S7), while the corresponding D for RPCMs adopts similar values, 2.48, 5.63, and 5.55 10-5 cm2/s, under TIP4P-Ew (γ = 1.0), TIP3P (γ = 1.0), and TIP3P (γ = 1.1), respectively (Table S8). On the whole, for RPCMs solvated using TIP4P-Ew at γ = 0.7 and SIRAH2.0/1.0, the number of water molecules located in the first hydration shell, i.e., from 0 to ~0.35 nm from the protein surface, is reduced compared to the pure TIP4P-Ew and TIP3P conditions, as well as their mobility D. It comes with a decrease in the mobility of the surface atoms of the complexes, RMSF, as reported earlier in the paper. Ligand-Ubiquitin contacts. All distance maps calculated from the AA MD trajectories contain the three typical extended contact areas of Ubiquitin (Figures 6 and S12), i.e., at the level 35



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

of strands β1 (residues 1 to 7), β3 (residues 40 to 45), and β5 (residues 66 to 72).

Page 36 of 48

A contact

analysis of the PDB structures using PDBsum31,32 shows that, among the 24 residues of Vps27 UIM-1, 11 of them are involved in 65 interactions with Ubiquitin, while ligand Pol ι UBM1 is characterized by 17 (over 38) residues involved in 110 interactions and one H-bond (Asp71Leu208) with Ubiquitin (Table S2). The case of 2G3Q is particular in that central residues Gly1353 to Leu1370 of ligand Ede1 UBA are not involved in close contacts with Ubiquitin. The folding of the Ede1 UBA ligand is such as a limited number of residues are actually in contact with Ubiquitin. On the whole, only 7 ligand residues over 43 are involved in 41 contacts and 2 H-bonds, the latter occurring between Asn1375 and Arg42. In each of the three complexes, Leu8 of β1, Arg42, Ile44, Gly47, Gln49 of β3 and β4, His68 of β5, and the hydrophobic patch Val70, Leu71, and Leu73 of β5 in Ubiquitin are involved. As for 1Q0W, TIP4P-Ew at γ = 0.7 allows to rather well reproduce the AA contact maps of structures 2MBB and 2G3Q (Figure 6), as well as SIRAH2.0/1.0 conditions which, as already mentioned, presents a slower convergence. At the AA MD level, the H-bond occurrence frequency analysis shows that His+68 of Ubiquitin is involved in H-bonds with all ligands. As already mentioned in Section a., His+68 of 1Q0W is H-bonded to Glu273 with an occurrence frequency of 49.8 % (Tables S10 and S11). Even at the RPCM level, His+68 stays H-bonded either to Ser270 with 17.0 %, or Glu273 with 37.1 %. Rather, in 2MBB, Hisδ68 H-bonding is detected as a rare H-bond type, especially with Gln77 with percentages of 1.3, 1.8, and 0.3 % under the AA TIP4P-Ew, RPCM TIP4P-Ew γ = 0.7, and RPCM SIRAH2.0/1.0 conditions, respectively. It is consistent with the findings of Burschowski et al. who also noticed that His68 of Ubiquitin interacts only weakly with UBMs.30 Contrarily to complex-water H-bonds, RPCM prevents intramolecular and ligand-protein Hbonds to be frequently formed. In 2MBB, at the AA TIP4P-Ew level, the H-bond occurrence 36


Page 37 of 48

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


frequency value of Asp71-Leu208, a contact also detected by PDBsum,31,32 is by far the most frequent with 86.3 % (Tables S2 and S11). At the RPCM level, it also appears under SIRAH2.0/1.0 solvation, with a strongly reduced occurrence frequency of 5.6 %. Under TIP4P-Ew at γ = 0.7, the occurrence frequency is even lower, i.e., 3.2 %. At the AA level, the second most frequent Hbond is formed between Ser91 and Arg274; it is also detected in the RPCM TIP4P-Ew MD trajectory at γ = 0.7, with a frequency of 8.2 %. In the case of 2G3Q, the Asn1375-Arg42 H-bond reported by PDBsum (Table S2) occurs in the AA MD conformational ensemble, but to the lesser extent of 14.4 % than the Glu1348-His+68 (37.9 %) and the Met1352-Gly47 H-bonds (43.0 %) (Table S11). No H-bond is characterized by occurrence frequencies above 50 %, contrarily to 1Q0W and 2MBB. Lys 6, His+68, and Arg74 of Ubiquitin are involved in H-bonds with various residues of the ligand, whatever the point charge distribution, especially with TIP4P-Ew water (Table S11). The analysis of the AA minimum distance maps shows that Leu8 and Ile44 of Ubiquitin are both involved in short contacts characterized by distances from 0.30 to 0.45 nm with the ligands (Figure S12 and Table S12). However, within the RPCM descriptions solvated using TIP4P-Ew at γ = 0.7, Leu8 and Ile44 are not involved in any close contacts with UIM-1 and UBA, respectively. For the RPCMs, there is no particular increase in the number and frequency of ligand-Ubiquitin H-bonds, except for 1Q0W where separation is now prevented, and for 2G3Q where the occurrence frequency of the Asn1375-Leu73 H-bond is significantly increased (38.5 %) (Table S11). At the AA level, similarity indices show that TIP3P with γ = 1.1 performs better than TIP3P in approaching the TIP4P-Ew result (Table S7). The Tanimoto index of the RPCM maps from their corresponding AA TIP4P-Ew version clearly shows that TIP4P-Ew at γ = 0.7 and SIRAH2.0/1.0 solvation allows to simultaneously generate values that are among the highest 37



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 38 of 48

for the three complexes (Table S8). Particularly, to approximate ligand-Ubiquitin contacts, TIP4P-Ew at γ = 0.7 is the best FF choice for structures 1Q0W and 2MBB, while SIRAH2.0/1.0 performs the best for 2G3Q for the current trajectory length.

Conclusions and perspectives In this work, we investigated the influence of various solvent models on the structure stability and the water-protein interface of Ubiquitin complexes modeled using different point charge distributions. In particular, a previously developed reduced point charge model (RPCM) has been tested and compared to its all-atom (AA) version for the complexes solvated by TIP4P-Ew and TIP3P water molecules involving either a scaling of the Lennard-Jones protein-Owater interaction parameters,23 or the coarse-grain (CG) SIRAH water description.24 Three ligands bound to Ubiquitin were studied, i.e., Vps27 UIM-1 (PDB code: 1Q0W)25, polymerase (Pol) iota UBM1 (PDB code: 2MBB)26, and Ede1 UBA (PDB code: 2G3Q)27. The ligand of complex 1Q0W is the most sensitive to the water force field (FF), due to its ability to bend due to a low α-propensity region occurring in its amino acid sequence. It was thus selected as the test case of the study. Results were further confirmed by Molecular Dynamics (MD) simulations carried out on complexes 2MBB and 2G3Q. The use of a RPCM rather than an AA charge distribution tends to increase the gyration radius of the complex, as well as the number of conformation clusters. It is due to a reduced number of stabilizing hydrogen bonds (H-bonds). MD simulations of rigid complexes show that the number of Ubiquitin complex-water H-bonds does not depend on the point charge model, while it is obviously linked to protein complex conformations and dissociation. Additional Ubiquitin complex-water H-bonds compensate less numerous intramolecular and ligand-Ubiquitin H38


Page 39 of 48

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


bonds. When the interactions with the solvent are increased using a scaling factor > 1, a separation between the ligand and Ubiquitin is observed. MD simulation results show that the best agreements between the RPCM and AA models are obtained either (i) with TIP4P-Ew at a scaling factor of γ of 0.7, or (ii) with the hybrid TIP3P/SIRAH FF with an initial Ubiquitin complex structure that is solvated in a large solvation box. In the last case however, convergence to a stable conformational ensemble is not satisfied. A decrease in γ, or the inclusion of SIRAH particles, allows to weaken the complex-solvent interactions, i.e., to reduce the influence of the solvent on the structure and dynamics of the protein complexes. It induces a slight collapse of the protein structure, thus involving a limited number of surface atoms in the complexes, partly preserved secondary structure elements, and a better approximation of the number of water-water H-bond occurring in a 0.35 nm shell surrounding the complex. A decrease in γ values also allows the water molecules of the solvation shell to be less packed around the protein system, as seen through the decrease in the number of such molecules and, thus, in the number of H-bonds they form between themselves. The dynamics of the surface atoms of the complexes and of the neighboring water molecules is slightly refrained when γ decreases. It also partly affects the diffusional behavior of the closest water molecules. Contrarily to the effect of scaling Amber99sb-TIP4P-Ew LJ interactions, SIRAH significantly affects the diffusion of water molecules compared to pure TIP3P. A scaling of the non-electrostatic terms in the FF thus partly compensates for the changes in electrostatic interactions induced by RPCM. A similar effect might possibly be obtained by modifying, e.g., the permittivity at a cut-off distance that remains to be determined. First attempts, not reported here, to modify the protein-water LJ interactions under TIP3P/SIRAH solvation conditions have shown that a scaling factor of 0.8 can be used to 39



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 40 of 48

stabilize MD trajectories of RPCMs, particularly to get convergence of the RMSD profiles of the protein complex and of each of its partners. However, such a procedure still does not seem to improve the secondary structure of the protein complexes. Additional investigations could be considered in that direction to understand how AA and CG water descriptions should be scaled. As SIRAH water beads are located farther than water molecules from the protein surface, it is expected that, at a given γ value, the scaling of protein-SIRAH interactions has a limited effect in comparison with the scaling of protein-TIP3P interactions. Different scaling factors could also be considered for atoms involved in H-bonds or not. Finally, a more comprehensive analysis could be achieved to clarify the apparent correlation between the dissociation of RPCM complexes and the experimental dissociation constants. In the paper, we focused on the approximation of protein AA MD conformational ensembles using RPCMs with modified protein-solvent interactions. The effect of protein-solvent modifications on the sampling efficiency of the RPCM-based MD simulations is a possible extension to our work. In that case, AA MD trajectories would no longer be considered as targets to the RPCM MD trajectories. In practice, the RPCM-based MD simulations should be continued at the AA level to determine whether the protein systems are stable, whether they return to their initial conformation, or whether they pursue their deconstruction.

Acknowledgments Frédéric Wautelet and Laurent Demelenne are gratefully acknowledged for program installation and maintenance. We also thank the reviewers for the extremely useful comments that helped to improve the paper. This research used resources of the ‘Plateforme Technologique 40


Page 41 of 48

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


de Calcul Intensif (PTCI)’ (http://www.ptci.unamur.be) located at the University of Namur, Belgium, which is supported by the F.R.S.-FNRS convention 2.5020.11. The PTCI is member of the ‘Consortium des Équipements de Calcul Intensif (CÉCI)’ (http://www.ceci-hpc.be). The authors also thank the Interuniversity Attraction Poles Programmes n° 7/05: ‘Functional supramolecular systems’ initiated by the Belgian Science Policy Office.

Supporting Information. The following files are available free of charge. Gromacs coordinate and topology files (ZIP); Detailed results calculated from the MD simulations of the Ubiquitin complexes (PDF) References 1

Raschke, T. M. Water Structure and Interactions with Protein Surfaces. Curr. Opin.

Struct. Biol. 2006, 16, 152-159. 2

Bellissent-Funel, M.-C.; Hassanali, A.; Havenith, M.; Henchman, R.; Pohl, P.; Sterpone,

F.; van der Spoel, D.; Xu, Y.; Garcia, A.E. Water Determines the Structure and Dynamics of Proteins, Chem. Rev. 2016, 116, 7673-7697. 3

Soares, C. M.; Teixeira, V. H.; Baptista, A. M. Protein Structure and Dynamics in

Nonaqueous Solvents: Insights from Molecular Dynamics Simulation Studies. Biophys. J. 2003, 84, 1628-1641. 4

Martini, S.; Bonechi, C.; Foletti, A.; Rossi, C. Water-Protein Interactions: The Secret of

Protein Dynamics. Scientific World J. 2013, 6, 138916.

41



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

5

Page 42 of 48

Caliskan, G.; Mechtani, D.; Roh, J. H.; Kisliuk, A.; Sokolov, A. P.; Azzam, S.; Cicerone,

M. T.; Lin-Gibson, S.; Peral, I. Protein and Solvent Dynamics: How Strong are they Coupled ? J. Chem. Phys. 2004, 121, 1978-1983. 6

Hinsen, K.; Kneller, G. R. Solvent Effects in the Slow Dynamics of Proteins. Proteins

2008, 70, 1235-1242. 7

Fenimore, P. W.; Frauenfelder, H.; McMahon, B. H.; Young, R. D. Bulk-Solvent and

Hydration-Shell Fluctuations, Similar to α- and β-Fluctuations in Glasses, Control Protein Motions and Functions. Proc. Natl. Acad. Sci. USA 2004, 101, 14408-14413. 8

Frauenfelder, H.; Chen, G.; Berendzen, J.; Fenimore, P. W.; Jansson, H.; McMahon, B.

H.; Stroe, I. R.; Swenson, J.; Young, R. D. A Unified Model of Protein Dynamics. Proc. Natl. Acad. Sci. USA 2009, 106, 5129-5134. 9

Takemura, K.; Kitao, A. Effects of Water Model and Simulation Box Size on Protein

Diffusional Motions. J. Phys. Chem. B 2007, 111, 11870-11872. 10

Makarov, V. A.; Feig, M.; Andrews, B. K.; Pettitt, B. M. Diffusion of Solvent around

Biomolecular Solutes: A Molecular Dynamics Simulation Study. Biophys. J. 1998, 75, 150-158. 11

Abseher, R.; Schreiber, H.; Steinhauser, O. The Influence of a Protein on Water

Dynamics in its Vicinity Investigated by Molecular Dynamics Simulation. Proteins 1996, 25, 366-378. 12

Bhattacharjee, N.; Biswas, P. Structure of Hydration Water in Proteins: A Comparison of

Molecular Dynamics Simulations and Database Analysis. Biophys. Chem. 2011, 158, 73-80.

42


Page 43 of 48

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


13

Kovacs, H.; Mark, A. E.; van Gunsteren, W. F. Solvent Structure at a Hydrophobic

Protein Surface. Proteins 1997, 27, 395-404. 14

Fogarty, A. C.; Laage, D. Water Dynamics in Protein Hydration Shells: The Molecular

Origins of the Dynamical Perturbation. J. Phys. Chem. B 2014, 118, 7715-7729. 15

Schröder, C.; Rudas, T.; Boresch, S.; Steinhauser, O. Simulation Studies of the Protein-

Water Interface. I. Properties at the Molecular Resolution. J. Chem. Phys. 2006, 124, 234907. 16

Laage, D.; Elsaesser, T.; Hynes, J. T. Water Dynamics in the Hydration Shells of

Biomolecules. Chem. Rev. 2017, 117, 10694-10725. 17

Dastidar, S. G.; Mukhopadhyay, C. Structure, Dynamics, and Energetics of Water at the

Surface of a Small Globular Protein: A Molecular Dynamics Simulation. Phys. Rev. E 2003, 68, 021921. 18

Bizzarri, A. R.; Cannistraro, S. Molecular Dynamics of Water at the Protein-Solvent

Interface. J. Phys. Chem. B 2002, 106, 6617-6633. 19

Leherte, L. Reduced Point Charge Models of Proteins: Assessment Based on Molecular

Dynamics Simulations. Mol. Simul. 2016, 42, 289-304. 20

Leherte, L.; Vercauteren, D. P. Comparison of Reduced Point Charge Models of Proteins:

Molecular Dynamics Simulations of Ubiquitin. Sci. China Chem. 2014, 57, 1340-1354. 21

Leherte, L.; Vercauteren, D. P. Coarse Point Charge Models for Proteins from Smoothed

Molecular Electrostatic Potentials. J. Chem. Theory Comput. 2009, 5, 3279-3298.

43



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 44 of 48

22 Leherte, L.; Vercauteren, D. P. Evaluation of Reduced Point Charge Models of Proteins through Molecular Dynamics Simulations: Application to the Vps27 UIM-1–Ubiquitin Complex. J. Mol. Graphics Model. 2014, 47, 44-61. 23

Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein-Water Interactions Improve

Properties of Disordered Proteins and Non-Specific Protein Association. J. Chem. Theory Comput. 2014, 10, 5113-5124. 24

Darré, L.; Machado, M. R.; Brandner, A. F.; González, H. C.; Ferreira, S.; Pantano, S.

SIRAH: A Structurally Unbiased Coarse-Grained Force Field for Proteins with Aqueous Solvation and Long-Range Electrostatics. J. Chem. Theory Comput. 2015, 11, 723-739. 25

Swanson, K. A.; Kang, R. S.; Stamenova, S. D.; Hicke, L.; Radhakrishnan, I. Solution

Structure of Vps27 UIM-Ubiquitin Complex Important for Endosomal Sorting and Receptor Downregulation. EMBO J. 2003, 22, 4597-4606. 26

Wang, S.; Zhou, P. Sparsely-Sampled, High-Resolution 4-D Omit Spectra for Detection

and Assignment of Intermolecular NOEs of Protein Complexes. J. Biomol. NMR 2014, 59, 5156. 27

Swanson, K. A.; Hicke, L.; Radhakrishnan, I. Structural Basis for Monoubiquitin

Recognition by the Ede1 UBA Domain. J. Mol. Biol. 2006, 358, 713-724. 28

Hicke, L.; Schubert, H. L.; Hill, C. P. Ubiquitin-Binding Domains, Nature Rev. 2005, 6,

610-621. 29

Hurley, J. H.; Lee, S.; Prag, G. Ubiquitin-Binding Domains. Biochem. J. 2006, 399, 361-

372. 44


Page 45 of 48

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


30

Burschowsky, D.; Rudolf, F.; Rabut, G.; Herrmann, T.; Matthias, P.; Wider, G. Structural

Analysis of the Conserved Ubiquitin-Binding Motifs (UBMs) of the Translesion Polymerase Iota in Complex with Ubiquitin. J. Biol. Chem. 2011, 286, 1364-1373. 31

de Beer, T. A. P.; Berka, K.; Thornton, J. M.; Laskowski, R. A. PDBsum Additions.

Nucleic Acids Res. 2014, 42, D292-D296. 32

PDBsum, Pictorial Database of 3D Structures in the Protein Data Bank.

http://www.ebi.ac.uk/pdbsum/, last accessed 27 February 2017. 33

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for

Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435-447. 34

Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.;

Smith, J. C.; Kasson, P. M.; van der Spoel, D., et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845-854. 35

Showalter, S. A.; Brüschweiler, R. Validation of Molecular Dynamics Simulations of

Biomolecules Using NMR Spin Relaxation as Benchmarks: Application to the AMBER99SB Force Field. J. Chem. Theory Comput. 2007, 3, 961-975. 36

Bernstein, F. C.; Koetzle, T. F.; Williams, G. J. B.; Meyer, E. F. Jr.; Brice, M. D.;

Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. The Protein Data Bank: A ComputerBased Archival File for Macromolecular Structures. J. Mol. Biol. 1977, 112, 535-542.

45



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

37

Page 46 of 48

Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-

Gordon, T. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665-9678. 38

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.

Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. 39

Darré, L.; Tek, A.; Baaden, M.; Pantano, S. Mixing Atomistic and Coarse Grain

Solvation Models for MD Simulations: Let WT4 Handle the Bulk. J. Chem. Theory Comput. 2012, 8, 3880-3894. 40

Gonzales, H. C.; Darré, L.; Pantano, S. Transferable Mixing of Atomistic and Coarse-

Grained Water Models. J. Phys. Chem. B 2013, 117, 14438-14448. 41

Machado, M. Tutorial 2. SIRAH Forcefield in GROMACS, Hybrid Solvation: Plugging

SIRAH

Solvent

to

Your

Atomistic

System.

seem.eu/images/trainingMaterial/LifeSciences/Tutorial_2_sirah4gmx.pdf,

https://training.vilast

accessed

24

February 2017. 42

Nerenberg, P. S.; Jo, B.; So, C.; Tripathy, A.; Head-Gordon, T. Optimizing Solute-Water

van der Waals Interactions to Reproduce Solvation Free Energies. J. Phys. Chem. B 2012, 116, 4524-4534. 43

Leherte, L. Applications of Leveling Methods to Properties of Small Molecules and

Protein Systems, In: Theoretical and Quantum Chemistry at the Dawn's End of the 21st Century, Eds. R. Carbó-Dorca, R.; Chakraborty, T. Apple Academic Press, in press. 46


Page 47 of 48

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


44

Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E.

Peptide Folding: When Simulation Meets Experiments. Angew. Chem. Int. Ed. 1999, 38, 236240. 45

Virtanen, J. J.; Makowski, L.; Sosnick, T. R.; Freed, K. F. Modeling the Hydration Layer

around Proteins: HyPred. Biophys. J. 2010, 99, 1611-1619. 46

Humphrey, W.; Dalk, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol.

Graphics 1996, 14, 33-38.

47



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 48 of 48

TOC graphic

Reduced point charge model γ = 0.7 Ubiquitin complex

γ = 0.8

MD with scaled (γ) protein-Ow LJ terms γ = 0.9

γ = 1.0

48


Reduced Point Charge Models of Proteins: Effect of ProteinâWater

Recommend Documents

Reduced Point Charge Models of Proteins: Effect of ProteinâWater