Conformational selection of vasopressin upon V1a receptor binding

Graphical abstract

VP ligand-receptor interaction has been the subject of longstanding research and drug design efforts [28,32,36,37]. Nevertheless, structural dynamics and details of its receptor interactions remain unclear, particularly due to the lack of crystal or electron microscopy structures, and therapeutic targeting of its interactions has had only limited success to date [36,37].
VP comprises a six amino acid containing, disulfide-cyclic subunit linked to an amidated three-residue C-terminal tail. Pro 7 functions as a hinge between these two subunits [38,39]. Arg 8 in the three-residue tail is a primary modulator of VP-receptor (VPR) selectivity (compared to OTR) [40,41], whose spatial position is crucial for binding [41]. Similarly, Pro 7 plays a role in V 1a R and V 2 R selectivity [41]. However, the structure-activity relationships (SARs) underpinning these features remain elusive.
The discrepancy between medical interest and lack of insights into structural dynamics is not least due to the limited experimental access to atomic-level details. In particular, nuclear magnetic resonance (NMR) spectroscopy -the major structure elucidation technique for non-crystallizable substrates in solution -is limited by the absence of long-range distance information (e.g., by nuclear Overhauser effects; NOE [42]). Consequently, combinations of NMR and computational techniques [38,43] have become increasing popular to reveal the complex dynamics of VP and related SARs. In this context, we capitalized in this work on an integrative approach combining in-silico modelling of the VP-ligand-receptor complexes and NMR-based structure validation. More specifically, we adapted virtual screening techniques [36,44,45] to assess the structural dynamics involved in VPR recognition. We used the VP-V 1a R interaction as a model system, given its particular importance for regulating cardiovascular functions.
With our approach, we demonstrate how conformational selection-type events upon encounter complex formation can be assessed. In particular, we can assess in-silico which structures within the heterogeneous peptide conformational space favor complex formation. Thereby, we introduce a new approach of 'virtual conformational space screening' at the example of a medically important target. This advancement of the virtual screening methodology enables several novel avenues for understanding peptide activities and ligand modifications in the context of SARs and drug development.
Our approach comprises two stages: First, we identified distinct sub-ensembles within the total VP structural ensemble by integrating molecular dynamics (MD) simulations and NMR spectroscopy. Then, we independently docked the identified sub-ensembles to an inactive V 1a R-homology model [46] and screened the resulting complex structures for variations in binding energies and dissociation constants to identify preferentially selected peptide conformations. V 1a R preferentially selected VP conformations featuring solvent-exposed tail units, as these enabled adoption of energetically favorable complex conformations. This contrasts with two other identified sub-ensembles that feature structural restraints due to intramolecular side-chain contacts or cis-configuration of Pro 7 and three-fold reduced binding energies.
Virtual screening techniques are well-established in the context of in-silico drug discovery [45,47,48] and were here expanded to peptide-receptor interactions. For virtual screening in drugreceptor interactions [44,49], a library of drugs is typically docked to a receptor and compounds displaying the desirable binding properties are selected for further analysis. We replaced the library of drugs with a library of VP conformational ensembles to assess the conformational selectivity involved in VP-V 1a R recognition.
Comparable approaches by Bonvin and co-workers were already successful in generating flexible protein-protein [50] and protein-peptide [51] complexes. This was achieved either by screening conformations observed in NMR structure predictions for varying binding energies [50], or by individually docking different peptides with varying degrees of conformational freedom to a target receptor [50].
Here, we first ran an MD simulation of VP in solution to obtain the library of conformations to be docked.

Vasopressin's conformational sub-ensembles
In our MD trajectories (all-atom simulations, pH 7.4, 37°C), VP sampled three well-distinguishable sub-ensembles including two in which Pro 7 adopted its natural trans state. One of these is characterized by a compact three-residue tail and the other by an extended tail (Fig. 1a), which is in agreement with reported conformational switches [52]. The third sub-ensemble was observed when Pro 7 adopts a cis state, in which the peptide tail also adopted an extended conformation. A compact sub-ensamble with Pro 7 in a cis state was not detected (Fig. 1b). We have denoted these three observed ensembles as C TRANS , E TRANS and E CIS indicating compacted (C) and extended (E) states of the three-residue tail as well as the cis and trans conformations of Pro 7 in the index.
In the case of the compact C TRANS tail conformation, Arg 8 folds transiently back towards Tyr 2 and is stabilized by side-chain interactions between both residues (see Fig. 1a). In particular, the Tyr 2 -OH g ---H g -Arg 8 hydrogen bonds appear to stabilize the compacted fold. This 'back-folded' structure is distinct from the extended solvent-exposed tail conformations E TRANS (Fig. 1a). The differential structural sampling of the two trans forms constitutes a bimodal conformational space, which is reflected in the intramolecular distances between atoms Arg 8 (H g ) and Tyr 2 (H g ). The conformational space is visualized by the distance distributions in Fig. 1b. Evidently, a bimodal distance distribution is observed. Arg 8 (H g )-Tyr 2 (H g ) distances r as short as 0.3 nm were detected for the compacted C TRANS fold. The distance distribution of this compact structural ensemble housing the back-folded state is centered Fig. 1. a) Representative snapshots of limiting cases for VP structures in the Pro 7trans extended (left), Pro 7 -trans compact (middle, most compacted) and Pro 7 -cis extended (right) state. The distance between residues Tyr 2 and Arg 8 is indicated by the dashed line. b) Distributions of Arg 8 (H g )-Tyr 2 (H g ) distances r found in MD simulations of VP in the Pro 7 -trans (left) Pro 7 -cis (right) state. While the former displayed a bimodal distribution with two distinct maxima, the latter had only a single maximum. The blue shade indicates distances extracted from the compact sub-ensemble. c) Trajectories of Arg 8 (H g )-Tyr 2 (H d ) distances r underlying the distributions in panel (a). The blue shaded area guides the eye to distances < 1 nm between residues Tyr 2 and Arg 8 . These were only observed for the Pro 7 -trans state. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) at~0.5 nm (Fig. 1b, left). By contrast, the extended mode that includes all other disordered states spans distances up to 2.5 nm whilst being centered at 1.9 nm (Fig. 1b, right).
When the conformation of Pro 7 switches from a trans to cis state, only extended peptide tails were sampled (Fig. 1b, right panel, monomodal histogram). Here, the Arg 8 (H g )-Tyr 2 (H g ) distances span a broad dominant distribution centered at $ 2.2 nm, consisting of species having separations ranging again from 1-2.8 nm. Fig. 1c displays the time trajectories underlying the distance distributions. It can be seen how VP in the Pro 7 -trans state switches between the two sub-ensembles, i.e., compact (indicated by the blue area) and extended. No such switches were observed for the Pro 7 -cis state.
The three different peptide tail ensembles were confirmed in three independent MD runs ( Fig. 1 shows one replica; see Supporting Information for the other two replicas; note that each run contains two independent data sets for the trans-and cis-conformations, respectively. Hence, six MD trajectories were computed in total). Note that we refrained from a quantitative assessment of the MD data due to the finite length of the trajectories. Instead, we restricted our interpretation on the observation of the three different sub-ensembles, but did not quantify their relative populations or the frequencies of conformational switches.

Virtual conformational space screening
Next, we randomly selected 50 structures from the three subensembles observed in the MD trajectories and docked these insilico to an established V 1a R model [46]. For computation, we used VINA (Vina Is Not Autodock [53,54]); further details can be found in the Experimental section. This 'virtual conformational space screening' approach is illustrated in Fig. 2a. From the all-atoms MD trajectories, snapshots of 50 different VP conformations were independently docked to the V 1a R model. The VP structures for each of the three sub-ensembles (C TRANS , E TRANS and E CIS ) that led to the lowest dissociation constants are depicted in Fig. 2b and the associated ligand-receptor complex structures in Fig. 2c.
It was interesting to observe how residues Tyr 2 and Phe 3 interacted with transmembrane helices (TMH) 3 and 4 in all cases to provide optimal complex stability (in agreement with published results [41]). In contrast, the position of the three-residue tail varied for the three docked cases C TRANS , E TRANS and E CIS . The preformed tail structures were conserved upon docking to the receptor which led to different complex conformations. Most prominently, tail residue Arg 8 interacted with TMH 6 in the E TRANS , but with TMH 5 for the C TRANS and E CIS states, while tail residue Gly 9 interacted with TMH 6 in the E TRANS and C TRANS states, but with TMH 4 in the E CIS state. These findings are in partial agreement with published structures of the VP-V 1a R complex. These displayed similar interaction patterns for the three-residue tail with TMH 4 [41,55]. Yet, the interaction with TMH 6 was not reported.
For statistical analysis, each docking experiment for the 50 selected VP conformations was repeated 25 times yielding 1250 VP-V 1a R complexes per sub-ensemble. Each computed complex structure was analyzed with respect to binding energies and dissociation constants K D . Fig. 2d displays box plots for the three subensembles C TRANS , E TRANS and E CIS , calculated in each case from the top five ligand-receptor complex structures with respect to binding energies. The same analysis was performed for three independent runs. In each run, the E TRANS state led to significantly lower dissociation constants compared to C TRANS and E CIS . The p-values for comparing the compacted and extended conformations the three runs were 1.0Á10 À5 , 1.6Á10 À3 and 1.9Á10 À1 . The full list of computed values is in the Supporting Information (Table S1). Our findings were further confirmed by comparison of the dissociation constants calculated for the complex structure with the best VINA score out of the 1250 structures (Fig. 2e). E TRANS led in all runs to two-to three-fold reduced dissociation constants. Note that several structures with comparably low dissociation constants have been computed (see Table S1). Since only results with the five highest VINA scores were used for computation of the box plots in Fig. 2d, the bottom edges converge to the values shown in Fig. 2e.

NMR validation of the VP structure: The heterogeneous conformational space of the VP Pro 7 -trans form
To validate the MD-based results, we experimentally probed VP's conformational space by means of solution-state NMR spectroscopy. Evidence for the tri-modal sampling space of the threeresidue tail was observed in 1 H-1 H NOESY and 1 H-1 H TOCSY experiments. An overlay of TOCSY (blue) and NOESY (red-green) NMR spectra of VP in phosphate buffered solution at pH 6.5 is presented in Fig. 3a. The resonances of Gly 9 and Arg 8 are both split into two clear signals along the direct dimension, indicating the presence of the cis and trans Pro 7 isoforms (further confirmation by repeated annealing is presented in the Supporting Information). Additionally, three features are notable in the NOESY data: (i) Arg 8 (H a )-Gly 9 (H a ) cross-peaks were observed for the Pro 7 -trans state, but not for the cis state, (ii) Arg 8 (H a )-Tyr 2 (H d ) cross-peaks were observed ( Fig. 3b; note that these NOEs contacts stem from different atoms than the H-bond described in Fig. 1; as the Arg 8 (H g ) could not be identified by NMR.), and, most importantly, (iii) the Arg 8 (H b )-Arg 8 (H N ) and the Arg 8 (H c )-Arg 8 (H N ) cross-peaks had a negative amplitude for the Pro 7 -cis and a positive one for the Pro 7 -trans state (with respect to positively phased diagonal peaks; green signals represent fast motions and red restricted motions; Fig. 3c).
The presence of these features confirmed the conformational sampling of the MD simulations for the acyclic Pro 7 -Arg 8 -Gly 9 tail. In particular, the positive amplitudes for the Arg 8   By contrast, in the Pro 7 -cis state, only negative cross-peaks were observed indicating faster side-chain dynamics. This validated the heterogeneous conformational space of the VP three-residue tail, as only transient intramolecular attraction between the tail and the macrocycle residues could account for these observations [52]. Given the observation of Arg 8 (H b )-Tyr 2 (H d ) cross-peaks, attractive interaction between tail residue Arg 8 and the Tyr 2 aromatic side-chain appear intuitive as an underlying force, which aligned well with the structures computed for the C TRANS case.
We conclude that VP, when Pro 7 adopts the trans state, indeed populates a bimodal conformational space due to transient sidechain interactions between Tyr 2 and Arg 8 that slows down conformational sampling, i.e., the intramolecular dynamics. In the Pro 7 -cis state, no indications of a conformationally restrained sub-ensemble were observed, and hence no compacted state could be inferred, in line with the MD results. Besides, the observation of slow conformational dynamics for the trans-state indicates that the C TRANS forms dominate the conformational ensemble when Pro 7 adopts the trans-state. Indeed, a predominance of extended conformations likely would result in fast backbone and side-chain dynamics and a similar cross-peak sign as observed for the E CIS structural ensemble.
Note that a line shape analysis supports this conclusion. The trans-state leads to a broader linewidth compared to the cis-state ( Fig. 3c and S5-S7 in the Supporting Information, the trans-state G 9 peak even overlaps with the cis-state one). This indicates, in agreement with the above interpretation, an exchange-based line broadening as the trans-state exchanges between two conformations. However, the cis-state only samples extended forms such that no exchange broadening is observed.
It has been suggested that both the cis and trans conformations of VP are populated under physiological conditions, owing to the similarity in energy of their respective proline backbone bond configurations [39,56,57]. At the same time, the energy barrier for bond rotation is sufficiently high to permit observation of the separate resonances for both isomers on the NMR time scale [58]. Indeed, previous NMR investigations of VP structure observed spectroscopic fingerprints of two isomers, observing Pro 7 -cis popu-lations in the range of 5-9% [56,38] under near-physiological conditions. Here, under our conditions (pH 6.5, 25°C), we were able to detect the presence both states simultaneously. To this end, we used the almost pure cis-form as starting material and obtained6 0% trans populations (as judged from the NMR peak intensities) after annealing at 60°C for 8 h.

Discussion
V 1a R is the main receptor in mediating these functions [59] and it is thus not surprising that significant efforts have been made to develop selective V 1a R agonists and antagonists for therapeutic applications as well as to study this signaling system [60,61].
Targeting V 1a R has been a long-standing research challenge [28,37,46,62,63] mainly due to the selectivity problems to its closely-related receptor subtypes, namely V 1b R, V 2 R and OTR. Our findings might therefore be of interest to the medicinal chemistry community for the development of more potent and selective V 1a R ligands. In particular, one could envisage to implement unnatural amino acids that shift the C TRANS /E TRANS population equilibrium to improve receptor affinity and potentially also receptor subtype selectivity. Our simulations indicate that the Tyr 2 -OH---NH-Arg 8 hydrogen bond particularly stabilizes the C TRANS state. Hence, modification of either the Tyr 2 or the Arg 8 residue might foster an overpopulation of the better binding E TRANS state, e.g., by introducing a methylated Tyr 2 modification that lacks the possibility to form the hydrogen bond in question. Alternatively, the binding affinity of VP to V 1a R could be modulated by directly manipulating the cis-trans equilibrium of the Pro 7 residue [64,65]. Such strategic modifications could deliver important potency and selectivity improvements for V 1a R binding, which is critical for therapeutic lead development.
Concerning the computed complex structures, we found that the dissociation constant of the VP-V 1a R complex depends on the spatial conformation of the VP three-residue tail. The lower dissociation constant for E TRANS indicates a preferential binding of the corresponding structural sub-ensemble to the receptor. From  Fig. 2c we deduce that the conformation with an extended peptide tail enables simultaneous interactions of VP residue Arg 8 with TMH 6, accompanied simultaneously by interactions aromatic sidechains of Tyr 2 and Phe 3 with TMHs 2 and 4, respectively. This complex configuration led to the highest binding energies in our simulations but was not observed in the E CIS nor C TRANS states, suggesting that this unique configuration drives the observed conformational selectivity. This is supported by other studies stating that the VP ligand-receptor complex is characterized by an intricate network of hydrogen bond interactions across several residues, rather than by a few well-defined points of contact [66]. The ligand-receptor interactions in the complex consist primarily of contacts between the residues of the cyclic part of VP and TMHs 3 and 4 and the extracellular loops of the receptor [55,66]. At the same time, the C-terminal tail modulates binding affinities via a few intermolecular contacts [66]. Indeed, SAR studies of VP revealed that modification or replacement of Pro 7 in VP can affect V 1a /V 2 receptor selectivity [67,68].
It should be noted that the herein computed differences in dissociation constants between the compacted and extended VP conformations likely do not suffice to lead to a 'pure' conformational selection-type binding event. Instead, a more complex binding process that features aspects of structural selectivity appears to be more probable. Moreover, due to the nature of virtual docking experiments with limited receptor flexibility, the computed dissociation constants correspond to conformations found upon ligandreceptor encounter. In other words, the E TRANS state is preferentially selected by V 1a R for encounter complex formation as the likelihood of dissociation is significantly reduced in comparison to the C TRANS and E CIS VP states . However, subsequent conformational sampling and structural adaption in the complex are not excluded by the presented methodology. In contrast, it is not unlikely that the selectivity towards a particular encounter complex structure can guide the VP-V 1a R system towards a particular final complex structure.
VP has received substantial attention from both, the structural biology [36,41,55] and medicinal chemistry communities [28,69,70], and V 1a R remains a drug target of interest. However, crystallization and electron microscopy efforts so far have failed to provide a high-resolution structure of this receptor. With longer MD trajectories, higher computing power, more refined V 1a R homology models, and high-sensitivity NMR becoming more accessible, deeper insights into the conformational characteristics of VP and its receptor interactions can be provided. We capitalized on these developments in this study and adapted a wellestablished in-silico methodology for screening receptor interactions of potential drug candidates to the screening of solution conformations for favorable binding properties. In classical virtual screening experiments, the rational selection of the ligand set often constitutes a bottleneck, but we could overcome this by confirming the selected conformations experimentally by NMR spectroscopy.

Conclusions
We developed a 'virtual conformational space screening' that provided new insights into the interactions of VP with V 1a R. We could establish that the binding event has different energetic preferences based on the sampled and experimentally verified conformations of the VP three-residue tail. In particular, the extended conformations of the Pro 7 -trans configuration (E TRANS ) led to a reduced dissociation constant and favorable binding energies.
Our results highlight the key role of the three-residue 'tail' in determining VP conformational selectivity, providing missing atomic-level detail to the current understanding of VP-V 1a R interaction -an aspect that opens novel research avenues for under-standing the functionality of the evolutionary selected conformational properties of VP, as well as guidance for ligand design strategies to provide more potent and selective VP analogues. Our VCSS is fast, easy to implement, applicable to other peptide-receptor systems, and data interpretation is straightforward. We expect this methodological advance to be interesting for a spectrum of chemists and biologists as it enables the elucidation of peptide-receptor recognition events at the atomic level.

NMR spectroscopy
For NMR spectroscopy VP was dissolved at 3 mg/mL in PBS buffer at pH 6.0 (90% H 2 O and 10% D 2 O; 25 mM Na 2 HPO 4 , 25 mM KH 2 PO 4 , 25 mM NaCl). NMR spectra were recorded on a 600 MHz Bruker NEO spectrometer equipped with a Prodigy TCI probe head. All spectra were acquired at 25°C. TOCSY and NOESY data were recorded using the Bruker 'dipsi2gpph19 0 and 'noesyfpgpphrs19 0 pulse sequences for TopSpin 4. For annealing, the samples were heated to a maximum of 60°C for 8 h. The high temperatures were necessary as the cis-trans conversion features high energy barriers and is rather slow. All TOCSY spectra were recorded with a spectral width of 8196.7 Hz in both dimensions and 32 scans. Mixing time was 150 ms. QUADRATURE detection was done using States-TPPI. All NOESY spectra were measured with width of 9615.3 Hz in the F2 dimension and 7202.1 Hz in the F1 dimension. We recorded 32 scans. The mixing time was 300 ms. QUADRATURE detection was employed again using States-TPPI sampling schemes.
Spectral processing was achieved using TopSpin, NMRPipe [71], and Sparky [72]. All data were zero filled to twice the original number of data points and apodized using a 60°shifted sine bell function prior to Fourier transformation. This was followed by a polynomial baseline correction in the frequency space.
Line shape analysis were conducted using home-written scripts for the MATLAB software package employing the 'fitnlorentian.m' function.

MD simulations
We performed all-atoms MD simulations of VP, beginning from either Pro 7 -cis-or -trans conformations. The published Pro 7 -trans VP structure was used as starting model [73]. To obtain the Pro 7cis form the N-C' bond was switched accordingly. This was followed in both cases by energy minimization and simulated annealing in explicit solvent (1% NaCl in water at pH 7.4). Finally, trajectories were recorded at 37 C at 2 fs time steps.
MD-simulations were performed using the YASARA softwarepackage [74,75]. The AMBER03 force field was employed with periodic boundary conditions [76]. Non-bonded interactions were cut off at 1.05 nm Long-range Coulombic interactions were treated by a smoothed particle-mesh Ewald method [77,78]. Non-canonic amino acids were built using YASARA and semiquantum-mechanically parameterized (YAPAC-AM1). MD trajectories of > 120 ns length were accumulated for the two systems. In total, six trajectories were computed, i.e., for both states (cis and trans) three trajectories were computed. Intermolecular forces were recalculated at every second simulation sub-step. Temperature rescaling was employed with a set-temperature of 37°C. The box dimensions (cubic of 37 Å side length) were controlled to yield a solvent pressure of 1 bar. Snapshots of the simulations were taken every 10,000 fs.

Receptor docking
Insights into ligand-receptor interactions were generated using VINA [53] as implemented in YASARA. We used a recently reported homology model of V 1a R [46]. 25 structures were generated for each docked structure. Side-chains flexibility was enabled upon docking the ligand structures. Dissociation constants were calculated from the binding energies DG binding following the Autodock method and conventions: The binding energy DG binding was calculated using the 'BindEnergy' macro implemented in YASARA, i.e., by determining the energy at infinite distance (between the object and the rest of the soup, i.e., the unbound state) and then subtracting the energy of the soup (the bound state). The more positive the binding energy, the more favorable the interaction. It should be noted that such determined energies are representing approximations that should only be compared among each other. They should not be understood as independent absolute measures of binding affinities. For each of the MD runs, the three sub-ensembles were identified, and the respective docking experiments were performed. This led to a total of nine docking experiments. The initial VP structure for docking to V 1a R were chosen randomly from the respective conformational sub-ensembles using the MATLAB random number generator. Box plots (Fig. 2) were calculated using the 'boxplot' function implemented in the MATLAB software package.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.