Allosteric effects in cyclophilin mutants may be explained by changes in nano-microsecond time scale motions

This work investigates the connection between stochastic protein dynamics and function for the enzyme cyclophilin A (CypA) in wild-type form, and three variants that feature several mutations distal from the active site. Previous biophysical studies have suggested that conformational exchange between a ‘major’ active and a ‘minor’ inactive state on millisecond timescales plays a key role in catalysis for CypA. Here this hypothesis is addressed by a variety of molecular dynamics simulation techniques. Strikingly we show that exchange between major and minor active site conformations occurs at a rate that is 5 to 6 orders of magnitude faster than previously proposed. The minor active site conformation is found to be catalytically impaired, and decreased catalytic activity of the mutants is caused by changes in Phe113 motions on a ns-μs timescale. Therefore millisecond timescale motions may not be necessary to explain allosteric effects in cyclophilins. The relationship between molecular motion and catalysis in enzymes is debated. Here, simulations of cyclophilin A and three catalytically-impaired mutants reveal a nanosecond-scale interconversion between active and inactive conformations, orders of magnitude faster than previously suggested.


Abstract:
There is much debate about the mechanisms by which molecular motions influence catalysis in enzymes. This work investigates the connection between stochastic protein dynamics and function for the enzyme cyclophilin A (CypA) in wild-type (WT) form, and three variants that features several mutations that are distal from the active site. Previous biophysical studies have suggested that conformational exchange between a 'major' active and a 'minor' inactive state on millisecond timescales plays a key role in catalysis for CypA. Here this hypothesis was addressed by a variety of molecular dynamic (MD) simulation techniques. The simulations reproduce X-ray crystallography derived evidence for a shift in populations of major and minor active site conformations between the wild-type and mutant forms. Strikingly, exchange between these active site conformations occurs at a rate that is 5 to 6 orders of magnitude faster than previously proposed. Further analyses indicate that the minor active site conformation is catalytically impaired, and that decreased catalytic activity of the mutants may be explained by changes in Phe113 motions on a ns-s timescale. Therefore previously described millisecond time scale motions may not be necessary to explain allosteric effects in Cy pA mutants.
3 A major goal of modern molecular biophysics is to clarify the connection between protein motions and enzymatic catalysis . [1][2][3] A wide range of experimental methods, e.g. neutron scattering, X-ray crystallography, NMR, or vibrational spectroscopy have been used to characterize internal protein motions occurring from femtosecond to second timescales. 4,5 While there is broad consensus that protein motions are implicated in catalysis, there is much debate around the role of conformational changes occurring on a millisecond timescale, and several studies have linked changes in millisecond protein motions with changes in enzymatic function. [6][7][8][9] However, it remains unclear whether such motions have a causal link to catalysis, or are merely a manifestation of the inherent flexibility of proteins over a broad range of timescales.
There have been vigorous debates about the meaning of dynamics in the context of enzymatic catalysis. [10][11][12] In the framework of transition state theory, the reaction rate is given by equation 1: where T is the temperature and R the gas constant. The pre-exponential term A(T) includes contributions from non-statistical motions such as re-crossing or tunnelling, The exponential term involves the activation free energy of the chemical step ∆ ‡ ( ). If transitions between reactant states are fast compared to the time scale of the chemical reaction, ∆ ‡ ( ) is the free energy difference between the thermally equilibrated ensembles describing the reactant and transition states. 13,14 Non-statistical motions described by A(T) have typically been found to make a small contribution to rate constants with respect to the exponential term that involves equilibrium fluctuations of the protein and solvent degrees of freedom. 15 The current work is concerned with the connection between rates of thermally equilibrated motions, and catalysis in enzymes. Specifically, the focus is on clarifying the nature of protein motions implicated in catalysis for the well-studied enzyme cyclophilin A (CypA). CypA is a member of the cyclophilin family of peptidyl-prolyl isomerases which catalyzes the cis/trans isomerization of amide groups in proline residues. 16 CypA plays an essential role in protein folding and regulation, gene expression, cellular signaling and the immune system. Notably, CypA is involved in the infectious activity and the viral replication of HIV-1. 17 Accordingly, CypA has been the subject of structure-based drug design efforts for decades. [18][19][20] Because of its significance as a medical target, the catalytic mechanism of CypA has been the subject of extensive studies. 2,3,[21][22][23][24][25][26][27][28][29][30] Computational studies have shown that the speedup of the of cis/trans isomerization rate of the prolyl peptide bond is a result of preferential transition-state stabilization through selective hydrogen bonding interactions in the active site of CypA. 26,30 Figure 1(a) depicts key interactions between the substrate and active site residues, whereas   31 More recently, two further mutants were reported in an effort to rescue the lost enzymatic activity of ST. These mutations were S99T and C115S (now only referred to as STCS), or S99T, C115S, and I97V (now only referred to as STCSIV). The two newly introduced mutants recover the enzyme activity to some extent, which correlates with an increase in → values. 32 While this body of work suggested a link between millisecond time scale motions and catalysis in enzymes, there is currently no detailed mechanistic explanation for the decreased catalytic activity of the mutants. The present study uses a variety of extensive equilibrium and biased molecular dynamics (MD) simulations to clarify the link between catalytic activity and rates of molecular motions of CypA in wild-type and the three mutant variants. in the 'minor' conformation is  1 -60°.This will be referred to as the 'out' conformation. In contrast, the 'major' state  1 60°, takes an 'in' conformation. In Figure 2 suggest that in apo WT the Phe113 'in' and 'out' orientations are equally likely, which is consistent with the relatively similar occupancies of the two rotamers in the X-ray structure (occupancies = 0.63 and 0.37 respectively). 31 In apo ST there is a significant population shift towards the 'out' orientation ( 1 =-60º), and the 'in' orientation has a marginal population (ca. 1%), see Figure 2(d). This agrees with the X-ray structure of ST where only the Phe113 'out' rotamer is observed (occupancy = 1.0). This also agrees with J-coupling measurements that show the dominant Phe113  1 angle is ca. -60 in ST. 31 In the STCS and STCSIV mutants the 'in' rotamer is also destabilized with respect to wild-type but to a lesser extent (populations of ca. 16% and 17% respectively). Though only one 'out' rotamer was resolved in the X-ray  Remarkably these values are five orders of magnitude faster than the exchange rates that have been determined by NMR measurements for motions involving Phe113. and Hamelberg have previously shown that fixed-charge classical force fields reproduce the energetics of amide bond rotation reasonably well due to relatively small changes in intramolecular polarization during this process . 28 The calculated activation free energy for the uncatalyzed cistrans isomerization process in water is consistent with experimental data (20.1  0.1 kcal•mol -1 vs ca. 19.3 kcal•mol -1 for the related substrate Suc-AAPF-pNA at 283 K). 37,38 The free energy profile for the substrate bound to CypA WT and ST in the 'in' conformation shows that the enzyme catalyzes the isomerization reaction in both directions via a transition state with a positive  value (ca. 90-100º) equally well (Figure 4(a)). There is a more significant decrease in activation free energy for transcis (ca. -6 kcal•mol -1 ) than for cistrans with less than 1 kcal•mol -1 difference between WT and ST, because the cis form is more tightly bound to CypA than the trans form. According to Figure 4 39 the free enzyme rapidly equilibrates between catalytically active and inactive forms before substrate binding (Figure 6(a)). In the mutants form, the interconversion rates between catalytically active and inactive forms are still within the s -1 timescale, but the equilibrium is shifted towards the catalytically inactive form (Figure 6(b)), thus the mutants are less pre-organized than WT and the overall catalytic activity is decreased.

Preorganization can explain observed decreases in catalytic activity of the mutants
In the case of the ST mutant and WT forms, Fraser et al. have reported bi-directional onenzyme isomerization rates ( → + → ) by NMR spectroscopy, and found a ratio of 68±13 between WT and ST. 31 According to the model proposed in Figure 6 and by combining the MSM-derived populations and the US-derived activation free energies, a ratio of 12 < 46 < 176 can be derived from the simulations (see SI for details). The uncertainty from the simulations is larger than that of the measurements because small variations in activation free energies contribute large change in catalytic rates. Thus the model described in Figure 6  A defining feature of this model is that the 1 rotamers of a number of active-site side-chains such as Gln63, Ile/Val97, Phe113, Cys/Ser115 flip in WT and mutants on ns -s timescales.
Back-calculation of C  -C  order parameters shows that this effect is captured by a decrease in S 2 values upon increasing the averaging window from 10 ns to 100 ns (Supplementary Figure S19).
Motions on these timescales are too rapid to be detected by CPMG or CEST experiments that have been used extensively to study s-ms processes in cyclophilin A. 3,25,31,40,41 Likewise relaxation experiments cannot detect motions on this timescale as they are limited to processes occurring faster than the tumbling time  c of cyclophilin A (ca. 10 ns). 42 Residual Dipolar Couplings (RDCs) can, however, provide information about dynamic orientation of inter-nuclear vectors on the supra- c time scale. 43 Such experiments have been reported for backbone and methyl-RDCs in ubiquitin. 43,44 Therefore the model predictions can be experimentally tested with combined nuclear spin relaxation and RDC based model-free analyses coupled with a labelling scheme that resolves 1 side-chain motions. 43,45

Discussion
This work highlights the potential of detailed molecular simulation studies to guide the interpretation of biophysical measurements for the elucidation of protein allosteric mechanisms. 46 Previous work has suggested that exchange on millisecond timescales between conformational states in CypA are linked to its catalytic cycle, 27 leading to a proposal for a slow exchange between a 'major' and a 'minor' state of a set of side chain rotamers linking distal residue Ser99 to active-site residues. 27,31 The present results do not support or reject this hypothesis because the MD simulations used here do not resolve motional processes occurring on timescales slower than microseconds. However a major finding of this study is that transitions between 'in' and 'out' rotamers of Phe 113 in WT and mutants occur on a time scale of ns-s, thus five to six orders of magnitude faster than suggested by earlier NMR relaxation dispersion measurements. 31 Nevertheless the simulations reproduce well the population shifts in Phe113 rotamers observed in room-temperature X-ray crystallography experiments. This suggests that the X-ray structures may have resolved motional processes occurring on a distinct timesc ales from the processes resolved by CPMG experiments. Indeed in reported CPMG experiments the millisecond motions of Phe113 are coupled to a large network of ca. 30 residues, 31 whereas the 1 rotameric flip observed in the simulations, is a largely local motion.
Nevertheless, the simulations suggest that a local 'in' to 'out' rotation of Phe113 is sufficient to abrogate catalysis in cyclophilin A, and variations of exchange parameters on the ns-s timescale between these two conformational states appear sufficient to explain the decreased catalytic activity of the ST, STCS, STCSIV mutants with respect to WT. Therefore it is advisable to carry out additional experiments to confirm the existence of Phe113 1 rotations on the ns-s timescale before causally linking catalysis to millisecond time scale motions. On the computational side, efforts should focus on advancing MD methodologies such that millisecond timescale processes observed in experiments can be resolved in atomistic details.
The contribution of protein flexibility on the ps -ns and s-ms timescales to enzymatic catalysis has been the focus of several computational and experimental studies. 3,8,10,13,15,25,27,31,32,47 Our work suggests that more efforts should be directed at resolving conformational processes on the ns -s timescale. This has important conceptual implications for enzyme design and optimization strategies.

Systems preparation
Models for apo/substrate bound human CypA of the WT and ST were prepared for MD simulations from PDB structures 3K0N (R=1.39 Å) and 3K0O (R= 1.55 Å) respectively. For apo STCS and STCSIV two structures were prepared from PDB structures 6BTA (R=1.5 Å) and 5WC7 (R=1.43Å) and also by mutating residues in WT using Schrödinger's Maestro. 47 For WT the major conformation of 3K0N (altloc A, occupancy 0.63) was retained. For STCS and STCSIV the residues with higher occupancy were chosen for initial structures . Supplementary Tables S1-S2 summarise all simulations conducted in this study. The proteins were solvated in a rhombic dodecahedron box of TIP3P water molecules with edges extending 1 nm away from the proteins and chloride counter-ions were added to neutralise the overall net-charge. The Charmm22* forcefield 48 was used to describe protein atoms in the apo simulations because previous work from Papaleo et al. has shown that this forcefield reproduces more accurately conformational changes in CypA. 33 Steepest descent minimized was used for 50,000 steps followed by equilibration for 100 ps in an NVT ensemble, and 100 ps NPT ensemble, with heavy protein atoms restraint using a harmonic force constant of 1000 kJmol -1 nm -2 .
Models of CypA WT and other mutants in complex with the Ace-AAPF-Nme substrate were prepared.The amber99sb forcefield was used for the complex simulations because Doshi and co-workers have reported optimised  angle parameters for amides to simulate cis/trans isomerisation reactions. 49 The crystal structure of the CypA-cis AAPF peptide complex (PDB ID: 1RMH) 50 Tables S3-S4 and Figure   S7). For 'out' US calculation the distances between the proline ring of substrate and the phenyl rings of Phe113 and Phe60 were restrained using a flat-bottom harmonic restraint with force constants of 200 kJmol -1 rad -2 and 300 kJmol -1 rad -2 , respectively. Simulations were performed serially initially for 7 ns, with the starting conformation for a given target angle  k taken from the preceding run performed at the neighbouring  k+Δω value. Each US was then extended to 20 ns. A total of 22 (substrate in solution) or 24 (substrate bound to protein) umbrellas were used. In order to estimate uncertainties of free energy profiles six repeats of the entire procedure were performed for 'in' and 'out' US. All simulations were carried out using a PLUMED2 patched version of Gromacs 5.0 with simulation parameters identical to the previously described apo MD simulation protocols unless otherwise mentioned. The weighted histogram analysis method (WHAM) was used to produce a free energy profile from the pool of US simulations. 65

Other trajectory analyses
Average proton-proton distances were derived as = 〈 −6 〉 https://github.com/michellab/CypAcatalysis_input and are also provided as Supplementary dataset S1. Other data that support the findings of this study are available from the corresponding author upon reasonable request.