Dynamic design: manipulation of millisecond timescale motions on the energy landscape of cyclophilin A

Proteins need to interconvert between many conformations in order to function, many of which are formed transiently, and sparsely populated. Particularly when the lifetimes of these states approach the millisecond timescale, identifying the relevant structures and the mechanism by which they interconvert remains a tremendous challenge. Here we introduce a novel combination of accelerated MD (aMD) simulations and Markov state modelling (MSM) to explore these ‘excited’ conformational states. Applying this to the highly dynamic protein CypA, a protein involved in immune response and associated with HIV infection, we identify five principally populated conformational states and the atomistic mechanism by which they interconvert. A rational design strategy predicted that the mutant D66A should stabilise the minor conformations and substantially alter the dynamics, whereas the similar mutant H70A should leave the landscape broadly unchanged. These predictions are confirmed using CPMG and R1ρ solution state NMR measurements. By efficiently exploring functionally relevant, but sparsely populated conformations with millisecond lifetimes in silico, our aMD/MSM method has tremendous promise for the design of dynamic protein free energy landscapes for both protein engineering and drug discovery.


Introduction
Capturing the complete dynamics of proteins to understand their many functions has long been a goal of structural biology. Experimentally, trapping proteins into relevant conformational states and characterising these at atomic resolution provides tremendous insight. 1-3 Such static 'snapshots' are inherently incomplete and especially in cases where proteins are particularly dynamic, understanding function can require characterisations of conformational states that are formed only transiently, and populated as sparsely as a few per cent. [4][5][6] These conformers are typically challenging to characterise experimentally. Solution-state NMR has proven unique in its ability to detect and structurally characterise these functionally relevant states, whose lifetimes can be on the order of a few milliseconds and are otherwise 'invisible' to experimental measures. [6][7][8] Such states can play important roles in processes as diverse as protein folding, molecular recognition and catalysis, but remain challenging to characterise. [9][10][11][12][13][14][15][16] To engineer proteins with new functions, and for drug discovery, there is a clear need to be able to characterise, explore and manipulate the full set of conformational states accessed by proteins, and the mechanism by which they interconvert.
Molecular dynamics (MD) simulations provide an attractive means to accomplishing this goal in silico. To characterise low amplitude structural uctuations in local minima, simulations on timescales of a few hundred nanoseconds may be sufficient. 17,18 Substantially longer simulations are required in order to sample sparsely populated 'excited' conformational states with millisecond lifetimes. Such a 'brute force' calculation is possible at present only for relatively small proteins using dedicated hardware that is not widely available. 19 To make progress, two broad categories of enhanced sampling methods have been developed with the goal of increasing total effective time covered by a simulation. Predened collective variables can be exhaustively sampled (metadynamics, [20][21][22] umbrella sampling, [23][24][25] steered MD [26][27][28] ), or the potential energy function of the forceeld can be explicitly altered (scaled MD [29][30][31][32] accelerated MD (aMD) 33,34 ). The former class of methods require the degrees of freedom underpinning the motion of interest to be supplied in advance, which pre-supposes that the relevant motions are already known. The latter class of methods can overestimate the population of high energy regions of the free energy surface (FES) and so explore unphysical regions of conformational space.
Aiming to improve upon these methods, Markov State Modelling (MSM) based sampling approaches have emerged as an exciting method for sampling protein dynamics on millisecond timescales. 2,[35][36][37][38][39][40][41] These techniques may be used to sample conformational space using a series of relatively short MD trajectories. These trajectories are analysed with the goal of obtaining a discrete set of interconverting conformational 'states', where each has no knowledge of any previously visited state. This method allows a series of independently generated MD trajectories to be combined into a complete description of the energy landscape of the protein. MSM have been shown to reproduce experimental observables in model systems where the structures of the relevant states were known in advance, 35,42,43 and in cases where the timescales for conformational transitions were sufficiently short to be validated with 'brute force' MD methods. 19,37,44 In this work, we generate a protocol suitable for sampling transitions between structurally ambiguous states separated by millisecond kinetic barriers (Fig. 1b). This is achieved by harnessing the ability of aMD to efficiently explore conformational space, while relying on the MSM procedures to recover a more accurate thermodynamic and kinetic description of the system. Because this approach provides atomically detailed structural ensembles and detailed information about how they interconvert, it is well suited to both generate experimentally testable mechanistic hypotheses and enable the rational design of conformational ensembles.
To manipulate the energy surface, we implemented a design strategy based on hydrogen bonding patterns to stabilise specically sparsely populated conformers with millisecond lifetimes (Fig. 2). This procedure led us to a previously unidentied mutation D66A. Whereas WT has the 70s loop predominantly in the closed state, consistent with the majority of X-ray structures of CypA solved to date, 52,55,[71][72][73][74][75][76] D66A was predicted to predominantly occupy the open state of the 70s loop. Similarly, as a negative control, a second similar mutation H70A was predicted by aMD/MSM to closely resemble the WT (Fig. 3).
Using a combination of R 1r and CPMG NMR experiments, with X-ray crystallography we tested these predictions ( Fig. 4 and 5) and veried that exactly as predicted, the 70s loop is substantially disordered in D66A, and adopts a conformation very similar to the minor conformations populated in WT ensemble. Also consistent with the predictions of aMD/MSM, the free energy landscape of H70A was largely the same as for WT.
Overall, we demonstrate that the aMD/MSM method accurately predicts the conformational uctuations of the challenging and dynamic protein CypA including states with millisecond lifetimes. The resulting free energy surface is sufficiently accurate to enable the design of mutants that stabilise specic conformations. aMD/MSM provides a powerful tool for not only the atomistic characterisation, but also the manipulation of the free energy surface of dynamic proteins, including controlling states with millisecond lifetimes.

Atomic-level insights in millisecond timescale CypA loop motion mechanisms
All-atom accelerated molecular dynamics (aMD) simulations 33,34 were initiated from an X-ray structure of wild type CypA (PDB ID 1AK4 72 ). In aMD conformational sampling is accelerated by adding a biasing potential to the molecular mechanics potential energy function when it falls below a threshold (ESI methods †). Structures sampled from the aMD trajectory were used to seed sets of equilibrium MD trajectories, which were used for Markov state modelling (MSM) analysis. In a MSM, the protein dynamics are modelled as a set of Markov jump processes between discrete conformational states. These Markov processes can be thought of as jumps between nodes in a network, with the edges representing transition probabilities between these nodes, and where no knowledge of previously visited nodes is required to estimate the probability for a particular transition. The procedure enables reconstruction of a conformational ensemble from sets of independently run short MD trajectories that have each individually explored only a small portion of that conformational ensemble. Three rounds of MD simulations and MSM estimations were required to produce a 100 microstate MSM (Fig. 1c) that was numerically stable, corresponding to a cumulative MD sampling time of 187.5 ms. Details about model validation against experimental observables are given in the ESI (Fig. S1-S5 †).
To obtain physical insight into the calculated ensemble, the 100 MSM microstates were further lumped into a coarser 5macrostate model (see ESI methods †). This ve state conformational landscape reveals that both the 100s and 70s loops can adopt 'open' and 'closed' conformations. While the interconversions of both loops were largely independent, the interconversions of the 100s loop was one order of magnitude faster than the 70s loop ( Fig. 1d and S6 †), suggesting that the slower loop dynamics of the protein are largely associated with the opening and closing of the 70s loop. The most populated state for WT CypA, (100s open/70s closed, Fig. 1c and d orange, 41 AE 6%) closely resembles the conformation typically observed in X-ray structures of CypA ( Fig. S7 †), and uctuations within the state were small (Ca RMSF values of 2.7Å and 1.1Å for the 100s and 70s loops, respectively). From this state, the 100s loop can rapidly close to adopt the 100s closed/70s closed state ( Fig. 1c and d, red, RMSF 2.6Å and 1.1Å). While these two states share a common conformation of the 70s loop, the 100s loop is in a signicantly different position with values of Ca RMSD with respect to the crystal structure >4.5Å. No major conformational uctuations outside the 70s and 100s loop regions are observed between macrostates (Fig. S8 †).
These To validate these ndings we sought a mutation that would stabilise the group of minor, 70s open, conformations (teal/ blue/purple). The number of hydrogen bonds each residue makes to the 70s loop was assessed in both the open and closed forms (Fig. 2b). Notably, D66 was found to adopt a large number of H-bonds in the closed form, and substantially fewer in the open form (Fig. 2c), making this residue stand out as 'designable' and provide a target for altering the free energy landscape of CypA. The critical role of D66 in the proposed 70s loop motions mechanism is not supported by a recent NMR structural ensemble that has been proposed to describe the CypA 70s loop conformational changes measured by CPMG experiments. 63 Analysis of the structures in that NMR bundle shows little changes in hydrogen bonding interactions with the 70sloop (Fig. 2b). Further none of the structures in the NMR bundle contain a 70s-loop conformation similar to that commonly observed in high-resolution X-ray structures of CypA (Fig. S7 †). aMD/MSM simulations of D66A were performed as described for WT CypA (Fig. 1). The resulting ensemble predicted that the population of the previous ground state was substantially reduced (from ca. 40% to ca. 15%, Fig. 3a), and the intermediate (teal) state was increased in a compensatory fashion to become the new ground state (from ca. 5% to ca. 40%). Overall with respect to CypA the population of states in which the 70s-loop is closed/open transitioned from ca. 70 : 30% to 25 : 75%. Moreover, the mean rst passage times between the conformational states were predicted to decrease by 10-15 fold with respect to wild type (Fig. 3a and c). Finally, as a negative control, we performed aMD/MSM simulations on the mutation H70A. 70 The hydrogen bonding pattern for this residue was predicted to be broadly similar when the 70s loop is either open or closed (Fig. 2c). Consistent with this, aMD/MSM simulations determined that the conformational ensemble for H70A should resemble WT (Fig. 3b).

Experimental validation of the calculated ensembles
To establish the validity of mechanistic hypotheses originating from the aMD/MSM models, we sought to compare our results to experimental data. Initially, we examined binding of an inhibitor against both WT and D66A. The K d for D66A was 12fold weaker than for WT (Fig. 4a), indicating a substantial change in the underlying protein conformations. Crystallisation of D66A proved challenging as we might expect for a protein with increased disorder, although both WT and D66A were successfully crystallised in the presence of an inhibitor 51 and their structures were solved (Table S2 †). Comparison of electron density proles revealed that while the 70s-loop adopts as expected a closed conformation in the WT complex (Fig. 4b  orange), electron density in the 70s loop for the D66A complex was largely absent (Fig. 4b blue), consistent with the notion that this loop occupies a predominantly disordered conformation in D66A.
To probe the conformational dynamics more directly, isotopically enriched samples of WT and D66A CypA were prepared for NMR analysis (Fig. S9 and S10 †). 1 H-15 N HSQC spectra of both proteins revealed that while both were predominantly well structured ( Fig. 4c and S11 †), there are substantial differences in the ground-state ensembles, exactly as predicted from the simulations. On a per-residue basis, the largest deviations in observed chemical shi positions between D66A and WT were found in the 70s loop region (Fig. 4e). The changes in chemical shi were quantitatively estimated by calculating the chemical shis of individual conformers in the aMD/MSM ensembles using the program shix2. 77 The calculated CSP values are consistent with the experimental observation of signicant changes in the 70s loop region, though the magnitude of the CSPs is underestimated (Fig. 4f).
To determine if there are any changes in fast timescale motions of the 70s loop, the random coil index (RCI) was calculated for WT and D66A from the observed chemical shis (Fig. 4d). These conrmed that the 70s loop is substantially more disordered in D66A, as expected from the simulations. This nding was further corroborated by measurement of the H-N steady-state heteronuclear NOE, where substantially lower values were recorded in the 70s loop for the D66A versus the WT. The heteronuclear NOEs were calculated from the MD ensembles, and also found to be in reasonable qualitative agreement with the experimental measurements ( Fig. 4g and h).
To further characterise the fast timescale dynamics, R 2 and R 1 values were additionally measured at two magnetic eld strengths, and a model-free analysis was performed. Such an analysis requires successively more complex motional models to be tted to individual residues, in order to achieve an overall self-consistent picture of the dynamics. The procedure was unable to provide a self-consistent picture of relaxation for either CypA or D66A, with large variations in parameters between adjacent residues and 'fast' local uctuations This journal is © The Royal Society of Chemistry 2020 Chem. Sci., 2020, 11, 2670-2680 | 2673 measured to be comparable to the global tumbling rate (Fig. S12 †). This is likely because of the substantial contributions of chemical exchange to the individual measurements (vide infra). Of the measurements obtained, the heteronuclear NOE nevertheless strikes a reasonable compromise between being relatively resistant to the effects of chemical exchange, whilst still being sensitive to relatively rapid local uctuations within the molecule (Fig. 4g and h).
While the fast ps/ns timescale dynamics indicate substantial differences in the 70s loop, the major goal was to determine if the structure of D66A resembles the excited state of WT CypA, thus validating the design strategy. A combination of R 1r and CPMG experiments were employed to study these conformational uctuations on the millisecond timescale (Fig. 5a, b and S13 †). The CPMG data for WT strongly resembles that obtained previously. 69,70 The combined data were analysed globally within a modied version of the soware CATIA. 78 Such an analysis can ascertain a minimally complex binding model, and provide experimental measures of the kinetic and thermodynamic parameters that govern the exchange. Moreover, the analysis also provides chemical shis (Fig. 5c and d), which characterise the structural changes that accompany the millisecond conformational rearrangements. For WT, a two-state model was found to explain well the combined dataset to within experimental error. The exchange rate k ex was determined to be approximately 2200 AE 90 s À1 and a population of excited state was determined to be approximately 2.0 AE 0.1% (Table S3 †).
Remarkably, the chemical shi differences directly measured from WT and D66A HSQC spectra were in good agreement with the chemical shi differences between ground and excited states of WT obtained by CPMG/R 1r measurements (Fig. 5e). This correlation strongly supports the in silico design strategy demonstrating that the ground state of D66A is similar to the 'excited' minor conformational state adopted by the WT, and that the design process effectively swapped the major and minor conformers.
Having established that the ground state of D66A resembles the 'excited' state of WT, we acquired CPMG and R 1r data to further characterise its millisecond dynamics. The raw exchange rates R ex for D66A in the 70s loop R ex were found substantially smaller than those detected for the WT protein, and residues that had large R ex values in WT were reduced to zero in D66A (e.g. A66, G74 Fig. 5) indicating the success of the design protocol. Consistent with this, the pattern of specic residues involved with the conformational dynamics was very different (Fig. 5 and Table S3 †) to WT. As with WT, the data could be well explained by tting to a 2-state model, where the exchange rate was similar to WT (k ex ca. 2000 AE 350 s À1 ), but the population of the excited state was much smaller, very close to our detection threshold (p B ca. 0.5 AE 0.1%). Analysis of the excited state chemical shis shows that this new state is not moving towards a particularly disordered state, nor the ground state of CypA. Instead, since a number of residues in the 80s loop were found to have relatively large R ex values in D66A this suggests that the excited state involves rearrangements of the 80s loop, outwith the scope of the present simulations that were designed to resolve millisecond timescale motions of the 70s and 100s loops. Overall the key nding is that the NMR data shows that the 'ground' conformational state of D66A has the 70s loop in a less structured conformation, consistent with what is expected for the 'open' conformations identied by MD simulations. The correlation of chemical shis suggests that this conformation is similar to the 'excited' state adopted by the WT protein.

Discussion
The work demonstrates that our new hybrid method aMD/MSM can efficiently provide a comprehensive description of protein free energy surfaces, including a characterisation of sparsely populated conformational states with millisecond lifetimes. The simulations are hypothesis generating and are sufficiently accurate to be amenable to rational design strategies whose goal is to adjust the relative populations of the different conformational states (Fig. 3).
It is interesting to compare different ensembles of WT CypA available in the literature. A notable ensemble of WT CypA 63 was calculated recently from solution-state NMR data. Our aMD/ MSM ensemble is also in good agreement with NOE, Jcoupling and RDC solution-state NMR data (Fig. S3 †). However while both ensembles suggest large amplitude motions of the 70s loop, the NMR ensemble does not support the critical role of D66 in a conformational change between open and closed loop conformations, and is also inconsistent with the closed 70s loop conformation observed by X-ray diffraction experiments.
To further test our method, we constructed an ensemble of CypA purely from our aMD trajectories. Unlike the aMD/MSM ensemble, the resulting aMD ensemble predicts that CypA adopts predominantly a 100s open/70s open conformation (Fig. S14 †), a prediction inconsistent with experimental data. Taken together, the aMD/MSM ensemble is the only method that characterises the conformational exchange, revealing the presence of both open and closed conformations of 70s loop in the WT with sufficient precision to provide the means to redesign the conformational landscape.
The resulting conformational states of CypA also provide insight into function. A double mutant (D66N/R69H) of the rhesus macaque homologue TRIMCyp adopts a 70s open conformation, which has been linked to potent HIV-2 restriction. 79 These observations lead to the interesting speculation that the sparsely populated open 70s conformation in human CypA plays an important and unrealised role in the pathway of HIV-2 infectivity, suggesting stabilisation of this conformation could lead to novel therapeutics.
Approaches have been developed based on using antibodies to stabilise transiently populated conformation states to facilitate ligand discovery efforts. [80][81][82] Our results open an exciting opportunity for determining this in silico as part of rational drug design endeavours. In the specic case of CypA, current treatment methods have insufficient isoform selectivity. 53,54 We envision that aMD/MSM can be used to selectively target transient states that are differently populated among different protein isoforms, hence paving the way for the rational design of isoform selective inhibitors.
There has been substantial discussion on the relationship between millisecond timescale motions in CypA and its catalytic cycle. 83 Initial measurements suggested that the two could be linked, although this has been disputed by later more detailed measurements. 64,68 Our results indicate that the opening and closing of the 70s loop is largely independent of the rest of the protein, supporting the notion that its uctuations are disconnected from the active site.
There remain technical challenges when comparing aMD/ MSM simulations to experimental data. A key issue is the different number of states used to explain the MD and NMR data. In the case of MD one starts with hundreds of thousands of snapshots that lead to a minimal 5-macrostates description aer application of standard MSM kinetic clustering methodologies. In the case of NMR the existence of states in intermediate exchange is inferred by tting phenomenological models to experimental data, leading almost always to a small number of states to avoid overtting the data. What constitutes a conformational state is also fundamentally different between both methodologies. To resolve different states in NMR a difference in chemical shis is needed, but this need not be the case in MD that uses a kinetic clustering approach.
To address this discrepancy it is necessary to calculate chemical shis from the structural ensembles and combine them appropriately to encapsulate the effects of intermediate exchange in order to accurately capture the effects of millisecond dynamics in NMR data. In this case, we found that there is a signicant degree of overlap of predicted chemical shis for residues in the 70s loop among MD macrostates, and the major separation is identied along the opening of the 70s loop (Fig. S15 †). We therefore suggest that the best mapping between CypA and D66A from HSQC measurements (x-axis) and from CPMG/R 1r measurements on the WT (y-axis). The agreement between the two is remarkably good. (f) Comparison of magnitude of Du between CypA and D66A from HSQC measurements (x-axis) and from CPMG/R 1r measurements on D66A (y-axis). The agreement is poor revealing that the 'excited state' in D66A is not the ground state of WT CypA.
MD and NMR states is along this variable between closed (orange + red) and open (teal + magenta + blue) states with a lowly populated intermediate state not evident in NMR experiments. We note that chemical shi calculation algorithms such as shiX2 are trained on static structures and not MD trajectories, and so their accuracy in such applications has not been extensively validated. Nevertheless, the broad agreement between the predicted and measured CSPs between WT/D66A and H70A/WT are encouraging. However agreement between chemical shis values calculated for one protein calculated in isolation is less impressive, perhaps suggesting that systematic deciencies in the chemical shi calculations are removed when one takes the difference between two simulated states. Moreover, using the expectation values for the CSPs from shiX2 with the rate constants predicted by the MSM leads to a large overestimation of R ex values, on the order of hundreds per second. This suggests that a combination of both the populations and the kinetics together with the estimated chemical shis are together not in perfect agreement. This can partially explained by noting that while the CPMG and R 1r data were acquired at 10 C, the simulations were performed at 25 C. These technical challenges aside, it is clear that the main predictions of aMD/MSM are fundamentally correct, and the key conformational ensembles predicted are sufficiently accurate to enable rational design efforts and to generate and test structural hypotheses.
Overall we demonstrate that aMD/MSM can provide atomiclevel insights into the thermodynamics and kinetics of independent conformational sub-ensembles that characterise the free energy surface of CypA. The picture is chemically intuitive, and sufficiently accurate for rational design efforts aimed at stabilising particular conformational states over others, even when the lifetimes of the individual states reaches the millisecond timescale. We anticipate that this protocol can greatly facilitate future protein structure-dynamics-function relationship studies and will unlock new opportunities in bioengineering and drug discovery.

Conflicts of interest
There are no conicts to declare.