{'fraction_non_alphanumeric': 0.02399398704902868, 'fraction_numerical': 0.02179694727104533, 'mean_word_length': 5.081926863572433, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Influence of a Ser111‐phosphorylation on Rab1b GTPase conformational dynamics studied by advanced sampling simulations

Rab GTPases constitute the largest branch of the Ras protein superfamily that regulate intra‐cellular membrane trafficking. Their signaling activity is mediated by the transition between an active GTP‐bound state and an inactive GDP‐bound state. In the inactive state the switch I and II segments adopt largely disordered flexible conformations, whereas in the active state these regions are in well‐defined conformations. The switch I and II states are central for recognition of Rab GTPases by interacting partners. Phosphorylation of the Rab1b‐GTPase at residue Ser111 (pS111) results in modulation of the signaling activity due to alterations of the protein interaction interface and also due to modulation of the conformational flexibility. We have studied the flexibility of native and pS111‐Rab1b in complex with GTP or GDP using extensive Molecular Dynamics (MD) simulations and an advanced sampling method called DIhedral Angle‐biasing potential Replica‐Exchange Molecular dynamics (DIA‐REMD). The DIA‐REMD method promotes backbone and side chain dihedral transitions along a series of replica simulations in selected protein segments and through exchanges also improves sampling in an unbiased reference simulation. Application to the Rab1b system results in significantly enhanced sampling of different switch I/II conformational states in the GDP‐bound Rab1b state. The pS111 modification is found to reduce the conformational flexibility even in the presence of GDP, which may influence signaling activities. The stabilizing effect can be attributed to the formation of additional surface salt bridges between Arg‐residues and pS111 not present in the native structure. The DIA‐REMD method could be a valuable approach for studying also other signaling proteins that contain flexible segments.

GTPases are associated with the switch I and II segments. Specific posttranslational modification (PTM) of amino-acid residues can further modulate the activity and regulation of small GTPases. [3][4][5] Such modifications include adenylylation, phosphorylation and phosphocholination of residues and reside either in the spatial vicinity or are directly located in the two conformational switch regions. [4][5][6][7] Growing evidence suggests that there are links between regulation of Rab proteins and Parkinson disease (PD)-related proteins. [8][9][10][11] PTEN-induced kinase 1 (PINK1), for instance, is a serine/threonine kinase that functions as a mitochondrial damage sensor, whose mutations cause autosomal recessive PD. 12 It has been shown that once activated, PINK1 regulates a number of Rab GTPases, namely Rab1b, 8a, 8b, and 13, through an indirect phosphorylation of a highly conserved residue, Serine 111. 10 This modification has been proven to diminish Rab8a's binding to its cognate GEF11 likely by interfering with the network of surface sidechain interactions. 13 The determination of the X-ray structure of Rab8a with phosphorylated Serine 111 (pS111) indicated only little difference compared to the native non-phosphorylated GTPase structures. 11 These include in the case of the pS111 variant small adjustments in the backbone structure of switch II with bound GTP and a better resolved residual switch II structure in the presence of GDP. 11 Also, Nuclear Magnetic Resonance (NMR) spectroscopy of the closely related Rab1b indicated no major structural alterations in solution due to the pS111 modification. However, these studies could not resolve the behavior of the switch I and switch II regions of the Rab1b GTPase. 11 Thermal melting experiments revealed an overall structural stabilization of the pS111-Rab1b both with bound GDP and in the presence of the GTP-analog GppNHp. 11 Hence, the influence of the pS111 modification on the conformational flexibility of Rab1b is still not understood. The increased stability of the pS111 variant indicates that the pS111 modification potentially affects the ensemble of conformational states of the switch regions not detected by the structural studies and that it can play a role for the interaction with other signaling proteins. To complement the biochemical and structural studies we employ comparative MD simulations to investigate the effects of this PTM on the conformational dynamics of Rab1b both in complex with GDP and GTP. MD simulations have played already a substantial role in gaining a deeper understanding of GTPases. For example, the catalytic cleavage of GTP at high resolution was studied by a combination of X-ray, Fourier-transform infrared spectroscopy (FTIR) and combined quantum mechanics and molecular mechanics (QM/MM) in atomic detail by Gerwert and coworkers. 14 Several computational studies have addressed the transition from active to inactive states and have illustrated the interplay of interactions between the protein's backbone atoms, the nucleotide's phosphate atoms and the Mg 2+ ion along the transition process. 15,16 Furthermore, accelerated MD simulations indicate that the mechanism by which the small GTPases bind to the nucleotide is based on conformational selection due to sampling of multiple conformations regardless of the nucleotide bound. Previous free energy simulations on the adenylylation of a Rab GTPase predicted a stabilization of the GTPase active form even in the presence of GDP 17 which was later also confirmed experimentally. 18 In the present study we investigate the influence of the pS111 modification on the Rab1b dynamics using extensive comparative MD simulations of the phosphorylated and non-phosphorylated S111 vari- The stabilization is mainly mediated by an interaction of the phosphate group with the R79 residue located in the switch II segment (and another Arg residue) and provides an explanation for the experimentally observed increased stability of the pS111 variant of Rab1b.
The study also demonstrates that the application of the DIA-REMD technique could be helpful to study other flexible protein or peptide segments that may play a role in signal transduction processes.

| Structural changes in S111phosphorylated Rab1b
In the next step, we studied the structural changes caused by the phosphorylation of Rab1b on S111 employing exactly the same DIA-  Figure S4).
Hence, the simulations indicate that S111-phosphorylation has a stabilizing influence on the switch I and II conformations compared to the wild type Rab1b. The cluster analysis based on the sampled switch I and II conformations resulted in fewer clusters (only five) and smaller conformational deviation from the start structure observed within those clusters ( Figure 3D).
Interestingly, the cluster analysis of the pS111-Rab1b:GDP complex revealed new hydrogen bonds involving the pS111 residuemainly formed between pS111 and the neighboring arginines (see  RMSDs based on a kernel density estimator. The plots were obtained from the DIA-REMD simulations of wild type and S111-phosphorylated Rab1b in the presence of GTP and GDP. In the wild type variants ( Figure 5A) and in the presence of GDP, the RMSD distribution covers a broad range, indicating a significant deviation from the conformation found in the active GTP-bound state. The time evolution of the sampled distribution shows reasonable convergence (see Supporting Information Figure S5). In contrast, for the GTP-bound Rab1b the distribution of RMSDs throughout the simulation was limited to a region near the GTP-bound start structure. Basically, the same distribution was obtained in the GTP-bound complex when S111 was phosphorylated, while in the GDP-bound pS111-Rab1b variant, we observed a contraction in comparison to the wild type. Overall, our data suggests that the phosphorylation of S111 may stabilize the switch regions and reduces the flexibility of the GDP-bound Rab1b.

| DISCUSSION
GTPases such as Rab1b alternate between active (GTP-bound) and inactive (GDP-bound) states which involves structural changes mainly in the switch I and II regions. Rab1b is a member of the Rab F I G U R E 2 Conformational changes in the switch regions during DIA-REMD simulations RMSD (non-hydrogen atoms) of the switch I, in A, B, and switch II segments, in C, D, with respect to the active Rab1b conformation during the DIA-REMD simulations (data from the reference replica). Amino-acid residues involved in the RMSD calculation of switch I and II are residues 31-43 and 68-79, respectively. Data point are taken in 20 ps intervals branch of GTPases and can undergo post-translational modifications modulating its function. Phosphorylation of S111 is frequently found in several Rab members (e.g. Rab1b and the closely related Rab8a) and in the case of Rab8a, experimental structures of both active and inactive forms with and without S111 phosphorylation are known. 11 Our simulation study demonstrates that even during relatively long Other enhanced sampling techniques have been used previously to sample conformational states of GTPases. They include enhanced sampling along specific reaction coordinates to drive the switch folding and unfolding transition using Umbrella sampling 17 or Metadynamics. 21 Alternatively, accelerated MD simulations using a biased/ deformed force field have also been used. 22 The disadvantage in the former approach is the necessity of defining a reaction coordinate prior the simulation and in the latter case it is the likely oversampling of irrelevant conformations, due to the altered force field. In contrast, in the present approach no definition of a reaction coordinate is necessary and only relevant conformations fully compatible with the original reference force field are consideredand thus, no re-weighting of sampled states is necessary. The technique could also be useful in issues, current force fields may over-stabilize folded structures relative to disordered conformational states. 23,24 Indeed, it has been demonstrated that current force fields yield varying results for disordered segments of proteins 24 and this could also explain the significant sampling of the active state in the presence of bound GDP. In the future, the DIA-REMD approach might be helpful to systematically test force field improvements for modeling disordered segments in proteins.
Our results suggest that the phosphorylation of S111 leads to a stabilization of the closed form in both switch regions of the GDPbound Rab1b complex, while it has no significant effect on the GTP-F I G U R E 4 Non-native Hydrogen bonds are formed in the pS111-Rab1b variant. Representative conformations of the most populated clusters obtained from the simulations of the pS111-Rab1b:GDP complex. In each case A-D, the percentage of the frames represented by the cluster is also indicated. The yellow dashed lines represent hydrogen bonds F I G U R E 5 Probability distribution of the switch regions' RMSD. Switch I and switch II amino acid residues deviations from the active Rab1b conformation for wild type (left) and S111-phosphorylated (right) variants. Probability density levels are given in 1/20 fractions of the probability maximum bound variant. The RMSD plots from the switch regions in Figure 5 show a limited range of sampling of open conformations in the presence of phosphorylated serine 111 in Rab1b bound to GDP compared to the wild-type variant. However, the X-ray crystal structures of the closely related inactive and active Rab8a states do not change significantly upon S111 phosphorylation. It is important to emphasize that our simulations do not exclude the possibility that the pS111-Rab1b switch I and II adopt disordered conformations also in solution, but indicate an ensemble of structures that is influenced by the pS111 modification. S111 resides in the vicinity of the switch II and once phosphorylated, it forms an intra-molecular bond with the side chain of R79 on the switch II and the R172 side chain. This arrangement was found in a large fraction of the entire simulation time (Figure 4) and can be interpreted as a likely reason for the reduced flexibility and stabilization of the structure in the presence of pS111. It is important to note, that such hydrogen bonding between pS111 and R79 was also found in the X-ray structure of the closely related pS111-Rab8a GTPase. 11 Furthermore, the pS111-Rab1b exhibits an enhanced thermal stability compared to the wild-type, 11 which can be explained qualitatively by the additional pS111-mediated hydrogen bonding. In addition to a possible stabilization of the switch region the appearance of a hydrogen bonded salt bridge between pS111 and R79 can also directly interfere with Rab-effector binding. 11,13 4 | MATERIALS AND METHODS

| Molecular dynamics simulations
All Rab1b simulations started from the crystal structure of Rab1b (PDB: 3nkv, the adenylylation modification at Y77 was removed). 7 The nitrogen atom between the βand γ-phosphates in GppNHp was exchanged with an oxygen atom to model the Rab1b:GTP structure.
The γ-phosphate was removed to generate a start model of the Rab1b:GDP complex. T72 and S111 residues were replaced with unprotonated phosphothreonine and phosphoserine to form the phosphorylated Rab1b structures. Amber ff14SB forcefield 25 was used for the proteins. Additional parameters for GDP, GTP, phosphothreonine and phosphoserine were taken from Amber parameter database at their fully unprotonated states. 26,27 An experimental study by Mann et al. 28 confirmed the unprotonated state of γ-GTP as the dominant protonation state. For Rab1b:GDP simulations, a Mg 2+ ion with two bound water molecules was placed next to GDP β-phosphate. Using sodium and chloride ions the salt concentrations were adjusted to 0.1 M and the complexes were solvated with TIP3P water model. 29 The solvated complexes were equilibrated by energy minimization (2500 steps of steepest descent) followed by heating to 300 K (100 ps) and 1 ns of density equilibration at constant pressure (1 bar) and a temperature of 300 K. During heating and equilibration, the protein heavy backbone and nucleotide atoms as well as magnesium ions were restraint with a harmonic potential at force constant of 5.0 kcal mol -1 Å -2 . All data gathering continuous MD simulations and replica-exchange (REMD) simulations were performed without any restraints. The pmemd.cuda module of the Amber 14 software package 30 with a time-step of 2 fs, periodic boundary conditions and the particle mesh Ewald (PME) method to account for long range interactions was employed. Trajectories were processed and analyzed using cpptraj program 31

| DIA-REMD technique
The DIA-REMD technique includes a biased potential that promotes transitions from low-energy backbone and sidechain dihedrals states.
The starting setup is composed of eight parallel MD simulationsreferred to as "replicas"-of the solvated Rab GTPase systems. Favorable dihedral angle combinations of the protein backbone (such as α-helices, β-sheets and left-handed α-helical regimes in the Ramachandran plot) as well as the first side chain dihedral angles of a preset peptide segment are penalized by adding a penalty potential along the replica simulations. An advantage of the approach is the possibility to limit the biasing to certain peptide segments in the protein structure (only the flexible switch I, residues 31 to 43, and switch II, residues 68 to 79, were included in the dihedral angle biasing). The first reference replica, however, runs under the control of the unmodified force field. Every 1000 steps an exchange was attempted between the neighboring replicas and it was allowed or rejected based on the Metropolis criterion. For the biasing potential on the backbone dihedral angles two-dimensional potentials that depend on the φ and ψ backbone dihedral angles were employed. The potentials have a maximum at the favorable states (φ c ,ψ c )) in the Ramachandran plot and fall off smoothly (parameters given in Table 1) . Let x i be the shortest distance between the position of residue i on the φ À ψ diagram (Ramachandran plot) and an energetically stable conformation, R(φ c ,ψ c ); The penalty force is a circular plateau with radius r 1 centered at R with the maximum potential (E max ), which continuously decreases left-handed α -helix (45, 45) 17. 5 10 down to zero at r 2 . Moreover, a one-dimensional potential was applied to promote side-chain rotations. Hence, the total potential was calculated by summing the 2D backbone and the 1D sidechain dihedral-angles penalty force.
8 > > > < > > > : : Test simulations showed that with eight replicas sufficient "mixing" between the windows was achieved. In order to assure high replica exchange rates, the heights of the potentials were adjusted on the fly based on an evaluation of the acceptance ratio during the last 100 exchange attempts; if any of the eight average rates fell below 20 %, the difference in the energy potential between replicas were reduced by 10 %. If all windows show 60 % of average successful exchange, then the difference is increased by 10 %. Typically, after less than 20 ns the biasing levels reached stable levels. Adjustment was therefore stopped after 50 ns and only data gathering beyond 50 ns was used for analysis.