Structural Basis for Species Specific Inhibition of 17β-Hydroxysteroid Dehydrogenase Type 1 (17β-HSD1): Computational Study and Biological Validation

17β-Hydroxysteroid dehydrogenase type 1 (17β-HSD1) catalyzes the reduction of estrone to estradiol, which is the most potent estrogen in humans. Inhibition of 17β-HSD1 and thereby reducing the intracellular estradiol concentration is thus a promising approach for the treatment of estrogen dependent diseases. In the past, several steroidal and non-steroidal inhibitors of 17β-HSD1 have been described but so far there is no cocrystal structure of the latter in complex with 17β-HSD1. However, a distinct knowledge of active site topologies and protein-ligand interactions is a prerequisite for structure-based drug design and optimization. An elegant strategy to enhance this knowledge is to compare inhibition values obtained for one compound toward ortholog proteins from various species, which are highly conserved in sequence and differ only in few residues. In this study the inhibitory potencies of selected members of different non-steroidal inhibitor classes toward marmoset 17β-HSD1 were determined and the data were compared with the values obtained for the human enzyme. A species specific inhibition profile was observed in the class of the (hydroxyphenyl)naphthols. Using a combination of computational methods, including homology modelling, molecular docking, MD simulation, and binding energy calculation, a reasonable model of the three-dimensional structure of marmoset 17β-HSD1 was developed and inhibition data were rationalized on the structural basis. In marmoset 17β-HSD1, residues 190 to 196 form a small α-helix, which induces conformational changes compared to the human enzyme. The docking poses suggest these conformational changes as determinants for species specificity and energy decomposition analysis highlighted the outstanding role of Asn152 as interaction partner for inhibitor binding. In summary, this strategy of comparing the biological activities of inhibitors toward highly conserved ortholog proteins might be an alternative to laborious x-ray or site-directed mutagenesis experiments in certain cases. Additionally, it facilitates inhibitor design and optimization by offering new information on protein-ligand interactions.


Introduction
Human 17b-hydroxysteroid dehydrogenase type 1 (17b-HSD1) catalyzes the NAD(P)H dependent reduction of the weak estrogen estrone (E1) to the biologically most active estrogen estradiol (E2; Fig. 1) [1]. This reaction, which represents the last step in E2 biosynthesis, takes place in target cells where the estrogens exert their effects via the estrogen receptors a and b. Besides their physiological effects, estrogens are involved in the development and the progression of estrogen dependent diseases (EDDs) like breast cancer, endometriosis and endometrial hyperplasia [2][3][4].
In the past few years, aromatase inhibitors have been intensively investigated for the treatment of EDDs [5][6][7] but they lead to unwanted side effects due to their strong reduction of estrogen levels in the whole body. Therefore reducing local E2 levels by inhibition of 17b-HSD1 is a promising therapeutic approach for the treatment of EDDs. An analogous intracrine concept has already been proved successful for the treatment of androgen dependent diseases such as benign prostatic hyperplasia and alopecia by using 5a-reductase inhibitors [8][9][10][11]. 17b-HSD2 catalyzes the reverse reaction (oxidation of E2 to E1; Fig. 1) and inhibition of this enzyme must be avoided for the therapeutic concept to work. However, specific inhibition of 17b-HSD2 in bone cells may provide a novel approach to prevent osteoporosis [12].
17b-HSD1 is a cytosolic enzyme that belongs to the superfamiliy of short-chain dehydrogenases/reductases (SDRs) [13]. It consists of 327 amino acid residues (34.9 kDa) and the active form exists as homodimer [14]. 17b-HSD1 comprises a Rossmann fold, associated with cofactor binding, and a steroid-binding cleft [15]. The latter is described as a hydrophobic tunnel with polar residues at each end: His221/Glu282 on the C-terminal side, and Ser142/ Tyr155, belonging to the catalytic tetrad, which is present in the majority of characterized SDRs [16], on the other side [17]. To date 22 crystal structures of 17b-HSD1 are available as apoform, binary or ternary complexes [18][19][20]. All crystal structures show an overall identical tertiary structure, while major differences have been identified only for the highly flexible bFaG'-loop. It is not resolved in ten crystal structures, while the remaining twelve showed high b-factor values for this area, which is an additional hint for the flexibility of the bFaG'-loop. In some crystal structures a short a-helix was observed in the loop region but its occurrence seems not to be dependent on the presence of steroidal ligands, cofactor or inhibitor. However, the position and length of the ahelix changes: in the apoform (PDB entry 1bhs) the helix is limited to the beginning of the loop while in presence of steroidal ligands and/or cofactor it is shifted to the end (PDB entries 1dht, 1equ, and 1iol). Further, dependent on the presence of cofactor and ligands, the bFaG'-loop can occupy three possible orientations: an opened, semi-opened, and closed enzyme conformation [21].
However, a distinct knowledge of active site topologies and protein-ligand interactions is a prerequisite for structure-based drug design and optimization. To further increase this knowledge, inhibition values obtained for one compound toward proteins, differing only in few residues might be advantageous. For this purpose, wild type proteins and their mutants carrying a set of point mutations can be used. As alternative, proteins from various species, which are highly conserved in sequence and differ only in few residues, might be considered. This latter approach was applied in the present study employing human and marmoset monkey (callithrix jacchus) 17b-HSD1. Selected human 17b-HSD1 inhibitors, representative of our structurally diverse inhibitor classes, were tested toward the marmoset 17b-HSD1. The resulting inhibitory potencies were compared with those obtained for the human enzyme and remarkable differences were only observed in the class of the (hydroxyphenyl)naphthols. In order to rationalize the speciesspecific inhibition profiles at the structural basis, a homology model of the marmoset enzyme was built using a human 17b-HSD1 x-ray structure as template. Further, the docking poses of selected compounds into both the human crystal structure and the modelled marmoset 17b-HSD1 were considered. Notably, the marmoset homology model and the docking poses of the inhibitors presented herein were validated by their ability to explain inhibition data. Subsequently, the complexes of two representative inhibitors, docked into the marmoset model and the human crystal structure, respectively, were subjected to MD simulations to investigate their conformational equilibrium. In addition, binding energy calculations as well as energy decomposition analysis were performed, with the aim to investigate the influence of the marmoset amino acid variations on the inhibitory potencies. The current work provides new insights into the marmoset 17b-HSD1 active site topology, reveals probable inhibitor binding modes in human and marmoset 17b-HSD1, and identifies amino acids responsible for species specificity in 17b-HSD1 inhibition.

Comparison of 17b-HSD type 1 and type 2 sequences
To identify regions that are conserved through several species, the sequences of rodent, cynomolgus, marmoset and human 17b-HSDs 1 were aligned ( Fig. 2A). The N-terminal region (residues , which constitutes the common Rossmann fold as well as the catalytic tetrad, is highly conserved for the analyzed species. Remarkable differences were observed in the F/G segment (residues 191-230), which is lining the steroid-binding site, and the C-terminal region (residues 231-285). Sequence alignment of human and marmoset 17b-HSD1 revealed that they share 80% sequence identity and 85% similarity ( Fig. 2A). Focusing on the steroid-binding site (residues 94 to 196 and 214 to 284), the identity increases to 87%, with five major amino acid variations observed in the marmoset enzyme: A191P, E194Q, S222N, V225I and E282N. In contrast, comparing cynomolgus and human 17b-HSD1, which show even 91% identity, only one of the aforementioned amino acid variations can be found (E282H; Fig. 2A). With the exception of marmoset 17b-HSD1, Ala191 and Glu194 are conserved through the analyzed species ( Fig. 2A) indicating the significance of the observed variations between human and marmoset. 17b-HSDs 1 of mouse and rat are 83% similar to the human enzyme in the first 287 amino acids. In both analyzed rodent enzymes His221, which is involved in steroid binding [17], is mutated into a glycine (H221G). Moreover, several other amino acids of the substrate-binding site are replaced by more bulky residues (L96F, N152H, M193Y/H, S222Y). Interestingly, rat and mouse 17b-HSDs 1 are significantly less sensitive to inhibition by steroidal inhibitors compared to the human ortholog [39] and different classes of non-steroidal potent human 17b-HSD1 inhibitors turned out to be only weak inhibitors of E2 formation in rat liver preparations [40]. This might be partially explained by the absence of His221 as interaction partner as well as by the reduced volume of the active site (M193Y/H, S222Y).
Sequence identities between human 17b-HSD2 and 17b-HSD2 of the selected species range from 61% (mouse) to 93% (cynomolgus) while the F/G segment and the C-terminal region show the most pronounced variability (Fig. 2B). The comparison of 17b-HSD1 with the correspondent type 2 enzymes of the analyzed species revealed, that the type 2 enzymes have about 80 additional N-terminal residues relative to the type 1 enzymes. Sequence comparison showed that they share a very low overall sequence identity (#25%) with major differences in the F/G segment and the C-terminal domain, which constitute the active site. However, some amino acid motifs, characteristic for SDR enzymes [18], are highly conserved: the T-G-xxx-G-x-G motif, the Y-xxx-K sequence and the N-A-G motif (Fig. 2B). For marmoset 17b-HSD2 only a fragment of the primary sequence is available, which is 143 amino acids in length and constitutes the Rossman fold whereas the F/G segment and the C-terminal region are missing (Fig. 2B). This segment is 29% identical to marmoset 17b-HSD1.
As human and marmoset 17b-HSD1 differ in few residues of the active site, comparison of the inhibitory potencies of selected inhibitors, observed toward both enzymes, is suitable to further increase the knowledge of active site topologies and protein-ligand interactions. In particular, the amino acid variations in the flexible bFaG'-loop, which is suggested to play a crucial role in ligand binding, will be used to elucidate the function of the loop in more detail. Inhibition of marmoset 17b-HSD1 and 17b-HSD2 Marmoset placental tissue was used as enzyme source and the proteins were partially purified following a described procedure [41]. Tritiated E1 was incubated with 17b-HSD1, cofactor and inhibitor. The separation of substrate and product was performed by HPLC. In an assay similar to the 17b-HSD1 test, marmoset placental microsomes containing 17b-HSD2 were incubated with tritiated E2 in the presence of NAD + and inhibitor. Labelled product was quantified after HPLC separation.
The inhibition values of compounds 1-20 ( Fig. 3) are shown in Table 1. The bis(hydroxyphenyl) substituted arenes 1-11 showed comparable or even higher inhibitory potencies toward marmoset 17b-HSD1 compared to the human enzyme. For example, compound 8, which is a potent inhibitor of human 17b-HSD1 (IC 50 = 151 nM), showed a stronger inhibitory potency toward marmoset 17b-HSD1 (IC 50 = 4 nM). Comparable inhibition data for human and marmoset 17b-HSD1 were also measured for the bicyclic substituted hydroxyphenylmethanones 12 and 13. This indicates that the binding of the two aforementioned inhibitor classes is reinforced or at least not affected by the amino acid variations in marmoset 17b-HSD1. In the class of the (hydroxyphenyl)naphthols, the human 17b-HSD1 inhibitor 14 (IC 50 = 116 nM) showed also a good inhibitory potency toward marmoset 17b-HSD1 (IC 50 .50 nM). The introduction of space filling substituents in position 1 of the naphthol core of 14 is beneficial for the inhibition of the human enzyme (16, IC 50 = 26 nM), but for marmoset 17b-HSD1 no increase in inhibitory potency was observed (16, IC 50 .50 nM) compared to the unsubstituted 14. Interestingly, further enlargement of the substituents in 1-position of 14 led to highly active human 17b-HSD1 inhibitors (19, IC 50 = 15 nM) while a reduced potency toward marmoset 17b-HSD1 (19, IC 50 .50 nM) was found.
The inhibitory potencies of compounds 1-20 toward marmoset 17b-HSD2 were also determined to prove whether marmoset monkey is a suitable species for in vivo evaluation of 17b-HSD1 inhibitors (Table 1). Remarkably, lower selectivity of compounds toward non-target marmoset 17b-HSD2 was observed, when comparing to human 17b-HSD2.

Homology modelling and MD simulations of marmoset 17b-HSD1
In order to obtain a more precise picture of the threedimensional structure of the marmoset 17b-HSD1, a homology model was generated. A set of 100 models was built with MODELLER 9v7 [42] using the ternary complex E1-NADPHhuman 17b-HSD1 as template. This complex was obtained by docking E1 into human 17b-HSD1 employing the Protein Data Bank (PDB) entry 1fdt with conformation B for residues 187-200 (in the following determined as 1fdtB). This 3D structure was chosen as it represents a ternary complex with NADP + and E2 in the closed enzyme conformation, as the bFaG'-loop is resolved, and as this protein structure was already successfully used in previous docking studies with the investigated compound classes [30,36]. The best model was chosen according to internal DOPEscore [43] and PROCHECK [44] tests. Notably, it presents a short a-helix (residues 190-196) in the region between the bF- sheet (residues 178-186) and the aG'-helix (residues 209-227), with Pro191 in the first turn of the helix (N1 position; Fig. 4), thus differing from the secondary structure of the template.
The selected model was further refined both as holoenzyme with NADPH and as ternary complex with NADPH and E1 by MD simulations. The trajectories of the two MD simulations were stable with a a-carbon root-mean-square deviation (RMSD) to the starting structure below 3.0 Å (Fig. 5). No major structural differences were observed when comparing the simulated holoenzyme and the ternary complex. Therefore, the following results, based on analysis of the ternary complex, also apply to the holoform.
After an initial inward rotation in the MD simulation, the short a-helix (residues 190-196) remained stable for the last 5 ns of the  MD simulation, with an average backbone RMSD to the final conformation of 1.03 Å . Remarkably, a similar behaviour cannot be expected for the region between the bF-sheet and the aG'-helix in the human enzyme. The differences in the crystal structures suggest a high flexibility for the bFaG'-loop: it is not resolved in ten crystal structures and the remaining twelve show high b-factor values for this area. Some crystal structures present a short a-helix in the loop region whose position and length varies dependent on the presence of steroidal ligands, cofactor or inhibitor. In the apoform (PDB entry 1bhs) the helix is limited to the beginning of the loop, whereas in presence of steroidal ligands and/or cofactor it is shifted to the end (PDB entries 1dht, 1equ, and 1iol). Further, in the human enzyme, the bFaG'-loop axis occupies different orientations dependent on the presence of cofactor and ligands. In the marmoset however, both the holoform and the ternary complex show a helix starting already at the beginning of the bFaG'-loop with its axis in only one conformation.
The presence of the newly formed a-helix (residues 190-196) induced a different orientation of the side chain of Met193. Compared to the template structure, Met193 protrudes deeper into the substrate-binding site and stabilizes E1 by hydrophobic interactions (Fig. 4). Furthermore, during MD simulation a kink in the loop between the bD-sheet and the aE-helix was observed. Thereby the side chain of Leu96 was brought closer to E1 allowing Van der Waals contacts (Fig. 4). Summarizing, in the final part of the MD simulation E1 was stabilized by both lipophilic interactions and hydrogen bonds: Leu96, Leu149, Met193, and Phe259 constrained the steroidal scaffold while Ser142/Tyr155 and His221/Asn282 interacted with the carbonyl oxygen in 17-position and the 3 OH-group, respectively. The latter residue took over the H-bond acceptor abilities of Glu282, which is involved in forming an H-bond with the 3 OH-group of E1 in the human enzyme.
Employing CASTp [45], the active site volumes of human and marmoset 17b-HSD1 were calculated. The above-described conformational changes of Leu96 and Met193 as well as the S222N and V225I mutations resulted in a reduced volume of the marmoset 17b-HSD1 active site (478 Å 3 ) compared to that of the human ortholog (627 Å 3 ).
The stereochemical quality of the holoenzyme and the ternary complex models obtained from MD simulations was checked with PROCHECK [44]. The majority of the residues of the investigated structures were found to occupy the most favoured regions of the Ramachandran plots, while the other residues occupied the additional allowed regions. In detail, in the ternary complex 78.1% of the residues were placed in the most favoured region, 20.2% in the additional allowed region, 1.2% in the generously allowed region, and only 0.4% in the disallowed region (for the holoform: 81.8%, 17.4%, 0.8% and 0%).

Molecular docking
As the potent inhibitors of human 17b-HSD1 12 and 19 are structurally diverse and exhibit different potencies toward . Three-dimensional structure superimposition of marmoset and human 17b-HSD1. Low energy structure of marmoset 17b-HSD1 (green) in complex with E1 and NADPH (orange) from MD simulation superimposed onto the x-ray structure of human 17b-HSD1 (grey, PDB code: 1fdtB). When two residues are indicated, the first corresponds to marmoset 17b-HSD1 and the second to human 17b-HSD1. doi:10.1371/journal.pone.0022990.g004 marmoset 17b-HSD1 they were chosen as representatives to rationalize the observed species-specific inhibition profiles. Employing AutoDock 4 [46], both compounds were docked into the active sites of marmoset and human 17b-HSD1 using the equilibrated, ternary complex after 5537 ps and the x-ray structure 1fdtB (see experimental part), respectively. The inhibitor poses used for further investigation were selected considering binding energy and statistical representativity (cluster population; Table S1) and are shown in Figure 6.
In the human enzyme both inhibitors are placed in the substrate-binding site and they occupy an apolar subpocket consisting of the following amino acids: Gly94, Leu95, Leu96, Asn152, Tyr155, and Phe192 (Fig. 6A). While the carbonyl group of compound 12 mimics the D-ring keto function of E1, forming H-bonds with Ser142 and Tyr155, its para-OH group resembles the 3-OH of E1, which interacts with His221 via an H-bond. The meta-hydroxyphenyl moiety is projecting into the subpocket, where it forms an additional H-bond with Asn152 and is stabilized by pp-interactions with Tyr155 and Phe192.
The (hydroxyphenyl)naphthol-core of compound 19 occupies the substrate-binding site and is stabilized by three H-bonds: the 2-OH group interacts with Ser142 as well as with Tyr155 and the OH group in meta position of the phenyl ring in 6-position interacts with His221. In this case, the sulfonamide substituted phenyl ring in 1-position of the naphthol core protrudes into the subpocket, where it is stabilized by H-bonds with Asn152 and with the -NHof the backbone of Leu95 (Fig. 6A).
Also in the marmoset enzyme both compounds occupy the substrate-binding site, but only 12 protrudes into the apolar subpocket (Fig. 6B). Due to the altered side chain conformation of Leu96 in the marmoset enzyme, compound 12 is slightly displaced toward the C-terminus (Fig. 6C). Regarding the interaction pattern, only minor changes were observed: the carbonyl group forms only an H-bond with Ser142 but for the para-OH group a second H-bond with Asn282 was observed (Fig. 6C).
Interestingly, compound 19 resulted in a completely different binding mode when docked into the homology model of marmoset 17b-HSD1 with respect to its position in the human crystal structure 1fdtB (Fig. 6D). The OH-group in meta position of the phenyl ring makes an H-bond with the backbone carbonyl oxygen of Cys185 and the 2-OH function forms H-bonds with His221 and Asn282 in a bifurcated fashion. The sulfonamide substituted phenyl ring is located in the C-terminal gate and might be stabilized by p-p-interactions with Phe259 and an H-bond with the backbone -NH-of Leu262.

Validation of the docking complexes by means of MD simulations and free energy calculations (MM/PBSA)
With the aim to validate the docking results and to unravel possible induced-fit mechanisms, different MD simulations were run in explicit aqueous solution. Distance restraints were applied to inhibitors only in the first ns of the MD simulations with the aim of maintaining their proper orientation. For the rest of the MD simulation no restraints were used and the whole complexes were left free to move. This was done in order to avoid trapping the inhibitor in an unstable conformation, which could bias the results.
The RMSD values of the heavy atoms of the inhibitors and of the Caatoms of the enzymes were analyzed as a function of time to assess the degree of conformational drift, as shown in Figure 7.
In the simulation of 12 bound to human 17b-HSD1 the CaRMSD of the protein as well as the RMSD of the heavy atoms of 12 showed a stable plateau (,2.2 Å ) from 1.0 to 2.5 ns (Fig. 7A).
After 2.5 ns the CaRMSD of the protein increased and a minor fluctuation of the heavy atom RMSD of 12 was observed. The latter finding could be related to a slight shift of the inhibitor toward the C-terminal end of the enzyme. Notably, the hydrogen bonds between 12 and 17b-HSD1, which were observed in the initial structure, were conserved during the 4 ns simulation suggesting that both protein and 12 fluctuations do not impact the inhibitor binding.
Both human 17b-HSD1 and compound 19 were stable in the simulation of their complex (Fig. 7B). During the simulation (after 1.5 ns) the hydrogen bond of the meta-OH group with His221 is replaced by an H-bond with Glu282. However, after 2.1 ns, this hydrogen bond is interrupted for 0.2 ns allowing the hydroxyphenyl ring to rotate freely around the axis of the bond to the naphthol core thereby inducing a minor fluctuation of the heavy atom RMSD of 19. After this short interruption the H-bond interaction with Glu282 is re-established. On the other side, the hydrogen bonds between the 2-OH group and Ser142/Tyr155 as well as between the sulfonamide moiety and Asn152/Leu95 endured constant.
During the MD simulation of 12 in complex with marmoset 17b-HSD1, the meta-hydroxyphenyl moiety of 12 moved further out of the subpocket. This was reflected by the minor fluctuation of the RMSD of the heavy atoms of 12 after 1.4 ns (Fig. 7C). While this motion caused the break of the H-bond with Asn152, it also placed the meta-OH group in an appropriate distance to Asn222, thus allowing a new H-bond formation.
The analysis of the MD simulation of the marmoset 17b-HSD1-19 complex revealed an overall stable CaRMSD of the protein. However, after 1.8 ns, a ,0.5 Å fluctuation of the heavy atom RMSD of 19 was observed (Fig. 7D) corresponding to the rotation around the axis of the bond between the hydroxyphenyl ring and the naphthol core. Thereby the H-bond with the backbone amino group of Leu262 was lost.
Furthermore, for each of the four MD trajectories absolute free energy (DG) and relative binding affinity (DG bind ) were calculated applying MM/PBSA methods and NMODE analysis ( Table 2). This was done for the following stable sectors: 1000 to 2500 ps for 12 in complex with human 17b-HSD1 (Fig. 7A), 2300 to 4160 ps for 19 in complex with human 17b-HSD1 (Fig. 7B), 1500 to 4160 ps for 12 in complex with marmoset 17b-HSD1 (Fig. 7C) and 1800 to 4160 ps for 19 in complex with marmoset 17b-HSD1 (Fig. 7D). All four complexes showed favourable DG values, ranging from 27.3 kcal mol 21 to 23.2 kcal mol 21 . The free energies observed for compound 12 in complex with human (DG = 24.7 kcal mol 21 ) and marmoset 17b-HSD1 (DG = 25.8 kcal mol 21 ) are in the same range. This is in accordance with the inhibitory activities of 12, which are comparable for both species (see Table 1). Regarding compound 19, the complex with human 17b-HSD1 shows a more favourable free energy (DG = 27.3 kcal mol 21 ) than the one with marmoset 17b-HSD1 (DG = 23.2 kcal mol 21 ). This is mainly due to poor entropic contributions in the latter case. Remarkably, this finding is in concert with the experimentally determined inhibition data: compound 19 is a highly potent human 17b-HSD1 inhibitor (IC 50 = 15 nM) with reduced activity toward the marmoset enzyme (IC 50 .50 nM).

Analysis of the binding interactions using MM/GBSA methods
Focusing on the specific interactions, which mediate the binding of 12 and 19 to human and marmoset 17b-HSD1, we have analyzed the interaction energies of both inhibitors with the residues of the binding sites, employing a pairwise per-residue energy decomposition analysis.
Inspection of the interaction energies with human 17b-HSD1 (Table 3)  The interaction energies of 12 and 19 with marmoset 17b-HSD1 are listed in Table 4. Remarkably, for inhibitor 12, the energy contribution of the H-bond with the marmoset Asn152 (21.6 kcal mol 21 ) is 2.6 fold reduced compared to the human enzyme, while for compound 19 it is almost lost (20.2 kcalmol 21 ). Both for 12 and 19 reduced energies were also observed for interactions with Gly94, Leu95, and Leu96. Interestingly, inhibitor 12 showed an interaction with Asn222 in the marmoset enzyme (23.3 kcal mol 21 ), which was not observed for Ser222 in  Table 3. Interaction energies between the inhibitors 12 and 19 and the proximal (4.0 Å ) binding site residues of human 17b-HSD1.

Discussion
When the three-dimensional structures of marmoset and human 17b-HSD1 are compared, one of the most striking features is the small a-helix including the residues 190 to 196. It is formed in the segment between the bF-sheet and the aG'-helix starting from the interface residue Thr190, which is half in and half out of the helix (N-cap position). In contrast to the human enzyme, where this region is highly flexible, as suggested by the different crystal structures, the a-helix stayed stable during the MD simulation in both the holoform and the ternary complex. The observed conformational stability might be explained by the presence of a proline in position 191 instead of an alanine. Proline is a favourable candidate for N1 position because of its own conformational properties: with only one rotatable angle it loses less entropy than other amino acids in forming an a-helix and thereby it should have some stabilizing influence [47]. Furthermore, in an analysis of sequence-structural characteristics in protein crystal structures, proline was found to be a favoured residue at N1 position. Especially the residue pair involving threonine at N-cap and proline at N1 position, which is observed for marmoset 17b-HSD1, has a high prevalence [48].
In order to analyse the influence of the conformational changes in marmoset 17b-HSD1 on ligand binding, docking studies with subsequent MD simulations, free energy calculations, and energy decomposition analyses were carried out. While the conformational differences between the marmoset and the human enzyme did not affect the binding mode of 12 remarkably, the suggested binding mode of 19 differed strongly in 17b-HSD1 of both species. One possible explanation for that might be the lower sterical demand of 12 compared to 19. However, the energy contribution of the interaction between 12 and the marmoset Asn152 is reduced, whereas it was outstanding in complex with the human enzyme. Due to the minimal shift of 12 in the marmoset binding pocket, the geometric parameters for the H-bond with Asn152 are no longer optimal. Interestingly, in the marmoset enzyme an additional interaction of 12 with Asn222 is observed, which seems to compensate the deficit in interaction energy due to the absent interaction with Asn152 resulting in comparable binding energies for 12 in complex with human and marmoset 17b-HSD1. The latter finding is in accordance with the inhibition data observed for compound 12 and validates the marmoset 17b-HSD1 model.
Considering compound 19, no particular interactions with the subpocket residues of the marmoset enzyme exist. Although weak interactions between 19 and the C-terminal region of marmoset 17b-HSD1 are observed, the binding free energy is less favourable compared to that calculated for the human 17b-HSD1-19 complex. As the C-terminal part of the enzyme has already been discussed as a potential product exit gate of the enzyme [21], inhibitor 19 might be solvent exposed. This is consistent with the unfavourable entropy term of this complex resulting in the least favourable free energy.
Obviously, the presence of a proline in the flexible loop region and the thereby induced conformational changes in marmoset 17b-HSD1 are decisive for the species specific inhibition of 19. On one hand interactions with subpocket residues like Asn152, recently discussed as relevant interaction partner [49], are prevented and on the other hand the inhibitor is forced in an unfavourable solvent exposed conformation.
The bis(hydroxyphenyl) substituted arenes (compounds 1-11) show similar or increased inhibitory potencies toward marmoset 17b-HSD1 when comparing to human 17b-HSD1. Recently performed docking experiments proposed a steroidal binding mode when the human crystal structure 1fdtB was used [33]. The high inhibitory potencies toward marmoset 17b-HSD1 are in concert with the modelled structure of marmoset 17b-HSD1 as steroid-like binding is not affected by the proposed conformational changes. Obviously, they even stabilize the bis(hydroxyphenyl) substituted arenes in the marmoset 17b-HSD1 binding pocket as indicated by the observed inhibitory potencies.
Differing inhibitory potencies toward human 17b-HSD1 and 17b-HSD2 may arise from sequence variations in the regions 94-196 and 214-284 (numbering according to 17b-HSD1), which might lead to differences in the active sites of the two human subtypes. A lower selectivity of compounds toward non-target marmoset 17b-HSD2 was observed, when comparing to human 17b-HSD2. Obviously, the differences in the active sites of marmoset 17b-HSD1 and 17b-HSD2 are less pronounced compared to the human orthologs. However, as the available marmoset 17b-HSD2 sequence is missing the F/G segment and the C-terminal part this hypothesis cannot be proved. Table 4. Interaction energies between the inhibitors 12 and 19 and the proximal (4.0 Å ) binding site residues of marmoset monkey 17b-HSD1. The validity of the presented homology model is further substantiated by its ability to explain the reduced inhibitory potency of C-15 substituted estrone derivatives toward marmoset 17b-HSD1 [39]. The substituents in 15-position of the steroid were designed to occupy the hole between the flexible bFaG'-loop and the aG'-helix in the human enzyme [50]. Together with the helix formation and the conformational changes in the bD/aEsegment, the S222N mutation limits the size of the hole in marmoset 17b-HSD1 and thereby might reduce the inhibitory potency toward the marmoset enzyme.

Conclusion
An elegant strategy to gain more knowledge of active site topologies and, in particular, of protein-ligand interactions is to compare inhibition values obtained for one compound toward ortholog proteins from various species, which are highly conserved in sequence and differ only in few residues. Thereby, such an approach can be a valid alternative to site-directed mutagenesis. As human and marmoset 17b-HSD1 enzymes meet these criteria, selected human 17b-HSD1 inhibitors were assessed for their inhibitory potencies toward marmoset 17b-HSD1. While a species specific inhibition profile was observed in the class of the (hydroxyphenyl)naphthols, representatives of the other evaluated compound classes showed similar or even higher inhibition compared to those observed for the human enzyme. Using a combination of computational methods, including homology modeling, molecular docking, MD simulation, and binding energy calculation, a reasonable model of the three-dimensional structure of marmoset 17b-HSD1 was developed and inhibition data were rationalized on the structural basis. In the marmoset 17b-HSD1, residues 190 to 196 form a small a-helix, which is obviously stabilized by the presence of a proline in N-cap position (residue 191) and induces conformational changes that affect ligand binding. Furthermore energy decomposition analysis highlighted the important role of Asn152 as interaction partner for inhibitor binding.
This work could not only offer a better understanding of the active site topologies and of the protein-ligand interactions, but also provides novel structural clues that will help to design and optimize potent human 17b-HSD1 inhibitors with improved inhibitory potency toward marmoset 17b-HSD1. This is an important step to turn compounds, which show a promising pharmaceutical profile, into candidates for in vivo evaluation. Thus, our combined computational approach could also be considered as a valuable tool to achieve this goal.

MD Simulations
MD simulations were performed using the AMBER 9.0 suite program [55]. The partial atomic charges for E1 and the inhibitors were derived from the molecular electrostatic potential (MEP) previously calculated using GAMESS [56], according to the RESP methodology [57]. For the protein, partial atomic charges were read from the AMBER 9.0 libraries. The AMBER99SB force field [58] was employed to define atom types and potentials for the protein, while the general AMBER force field (gaff) [59] was used to define all needed atom types and parameters for E1 and the inhibitors. For NADPH (charge 24), the parameters previously reported by Ulf Ryde were applied (http://www.teokem.lu.se/,ulf/).
The input files for the MD simulation were prepared with the xLEaP module of AMBER. Each system was solvated with an octahedral box of TIP3P water molecules of 10 Å radius and neutralized by the addition of Na + ions. Finally, for each complex the topology and the coordinate files were written and used in the MD simulations.
Before starting the production-run phase, the following equilibration protocol was applied to all systems. At the beginning the system was energy-minimized in two stages: firstly, the solvent was relaxed while all the solute atoms were harmonically restrained to their original positions with a force constant of 100 kcal mol 21 Å 22 for 1000 steps; and secondly, the whole molecular system was minimized for 2500 steps by conjugate gradient. Subsequently, the system was heated during 60 ps from 0 to 300 K at constant volume conditions (NTV, PBC conditions), and then equilibrated keeping both temperature and pressure constant (NTP, PBC conditions, 300 K, 1 atm) during 100 ps. Electrostatic interactions were computed using the Particle Mesh Ewald method [60], and the SHAKE [61] algorithm was employed to keep all bonds involving hydrogen atoms rigid. NADPH, E1 and the inhibitors were constrained during the equilibration with a force constant of 20 kcal mol 21 Å 22 . After equilibration, a MD production stage (NTP, PBC conditions, 300 K, 1 atm) was performed. The total simulation length differed for the various complexes ranging from 4 to 6 ns. Distance restraints were applied to substrate/inhibitor with the aim of maintaining their proper orientation at the beginning (first ns) of production stage.
For the ternary complex of marmoset 17b-HSD1 with E1 and NADPH, two additional distance restraints were used: between the keto-oxygen of E1 and the side chain oxygen of Ser142 (d = 2.40-3.00 Å ; force constant: 10 kcal mol 21 Å 22 ) and between the oxygen of the OH-group in 17-position of E1 and the NE2 nitrogen of the His221 side chain (d = 2.70-3.40 Å ; force constant: 10 kcal mol 21 Å 22 ).
For the ternary complex of human 17b-HSD1 with 12 and NADPH, three additional distance restraints were used: between the keto-oxygen of 12 and the side chain oxygen of Tyr155 For the ternary complex of marmoset 17b-HSD1 with 12 and NADPH three additional distance restraints were used: between the keto-oxygen of 12 and the side chain oxygen of Ser142 (d = 2.50-3.10 Å ; force constant: 10 kcal mol 21 Å 22 ), between the oxygen of the OH-group in meta-position of 1 and the OD1 oxygen of the Asn152 side chain (d = 2.40-3.00 Å ; force constant: 10 kcal mol 21 Å 22 ), and between the oxygen of the para OHgroup of 12 and the OD1 oxygen of the Asn282 side chain (d = 2.50-3.10 Å ; force constant: 10 kcal mol 21 Å 22 ).
For the ternary complex of human 17b-HSD1 with 19 and NADPH three additional distance restraints were used: between the oxygen of the OH-group in 2-position of the naphthol core and the side chain oxygen of Tyr155 (d = 2.80-3.40 Å ; force constant: 10 kcal mol 21 Å 22 ), between the oxygen of the OHgroup in meta-position of the phenyl ring and the NE2 nitrogen of the His221 side chain (d = 3.00-4.20 Å ; force constant: 10 kcal mol 21 Å 22 ), and between the nitrogen of the sulfonamide moiety of 19 and the OD1 oxygen of the Asn152 side chain (d = 2.70-3.40 Å ; force constant: 10 kcal mol 21 Å 22 ).
For the ternary complex of marmoset 17b-HSD1 with 19 and NADPH two additional distance restraints were used: between the oxygen of the OH-group in 2-position of the naphthol core and the NE2 nitrogen of the His221 side chain (d = 2.70-3.40 Å ; force constant: 10 kcal mol 21 Å 22 ) and between the oxygen of the OHgroup in meta-position of the phenyl ring of 19 and the backbone carbonyl oxygen of Cys185 (d = 2.40-3.00 Å ; force constant: 10 kcal mol 21 Å 22 ).
Trajectories were analyzed using the AMBER ptraj module, the MMTSB toolset [62] and the molecular visualization program VMD (Visual Molecular Dynamics) [63]. The resulting low energy structures were extracted for the homology model (apoform and ternary complex) and subjected to a subsequent minimization of 1000 steps (500 steps of steepest descent followed by 500 steps of conjugate gradient), using the sander module of AMBER. The modified generalized Born solvation model (IGB = 2) [64] was used. Active site volumes of low energy structures were calculated using the CASTp [45].

Molecular docking
The three-dimensional structures used for docking studies were either retrieved from the PDB (1fdt, conformation B for residues 187-200) or from the homology modelling with subsequent MD simulation (equilibrated ternary complex after 5537 ps). The cocrystallized E2 and water molecules were removed from the PDB file. Hydrogen atoms and neutral end groups were added, NADP + was turned into NADPH and correct atom types were set. Ionization states and hydrogen positions were assigned using the Protonate 3D utility of MOE2009.10 (Chemical Computing Group Inc., Montreal, Canada). Ligand structures were built in MOE and RESP charges were assigned as described above. The 17b-HSD1 three-dimensional structures and ligand structures were prepared for docking studies through the graphical user interface AutoDockTools4 [46]. For the ligands, non-polar hydrogen atoms were deleted, rotatable bonds were defined and RESP charges were kept. For the protein, non-polar hydrogen atoms were deleted and charges were added to the structure. Autodock4.2 [46] was used to dock the ligands in the steroidal binding site of the processed protein structures. A box, centered on the steroid-binding site, was set to define the docking area. Grid points of 90690690 with 0.250 Å spacing were calculated around the docking area for all the ligand atom types using AutoGrid4.2. For each inhibitor, 50 separate docking calculations were performed. Each docking calculation consisted of 25610 5 energy evaluations using the Lamarckian genetic algorithm local search (GALS) method. Each docking run was performed with a population size of 250. A mutation rate of 0.02 and a crossover rate of 0.8 were used to generate new docking trials for subsequent generations. The docking results from each of the 50 calculations were clustered on the basis of root-mean-square deviation (RMSD = 2.0 Å ) between the Cartesian coordinates of the ligand atoms and were ranked on the basis of the free binding energy.
Free energy calculations using the MM/PBSA method The calculation of binding free energy was evaluated using the MM/PBSA (Molecular Mechanics/Poisson Boltzmann Surface Area) method as implemented in AMBER11 [65]. The electronic and Van der Waals energies were calculated using the sander module in AMBER11. The solvation free energy contributions may be further decomposed into an electrostatic and hydrophobic contribution. The electrostatic portion is calculated using the linearized PB equation. The hydrophobic contribution is approximated by the LCPO method [66] implemented within sander. The changes in entropy upon ligand association DS are estimated by normal mode analysis. For stable plateaus of the MD trajectories, snapshots were collected every 20 th frame (every 20 ps) and used to calculate relative binding affinity (DG bind ) and absolute free energy (DG).
Energy decomposition using the MM/GBSA method A free energy decomposition of the protein ligand complexes was performed on a pairwise per-residue basis using the MM/ GBSA (Molecular Mechanics/Generalized Born Surface Area) method as implemented in AMBER11. The GBSA implicitsolvent solvation model was used in order to avoid the retarding effect of the PBSA method.
Inhibition of 17b-HSD1. Inhibitory activities were evaluated by an established method with minor modifications [41]. Briefly, the enzyme preparation was incubated with NADH [500 mM] in the presence of potential inhibitors at 37uC in a phosphate buffer (50 mM) supplemented with 20% of glycerol and EDTA (1 mM). Inhibitor stock solutions were prepared in DMSO. The final concentration of DMSO was adjusted to 1% in all samples. The enzymatic reaction was started by addition of a mixture of unlabelled-and [2, 4, 6, 7-3 H]-E1 (final concentration: 500 nM, 0.15 mCi). After 10 min, the reaction was stopped by the addition of HgCl 2 (10 mM) and the mixture was extracted with diethylether. After evaporation, the steroids were dissolved in acetonitrile. E1 and E2 were separated using acetonitrile/water (45:55) as mobile phase in a C18 reverse phase chromatography column (Nucleodur C18 Gravity, 3 mm, Macherey-Nagel, Düren) connected to an HPLC-system (Agilent 1200 Series, Agilent Technologies, Waldbronn). Detection and quantification of the steroids were performed using a radioflow detector (Agilent 1200 Series, Agilent Technologies, Waldbronn). The conversion rate was calculated after analysis of the resulting chromatograms according to the following equation: %conversion~% E2 %E2z%E1 |100. Each value was calculated from at least three independent experiments.
Inhibition of 17b-HSD2. The 17b-HSD2 inhibition assay was performed similarly to the 17b-HSD1 procedure. The microsomal fraction was incubated with NAD + [1500 mM], test compound and a mixture of unlabelled-and [2, 4, 6, 7-3 H]-E2 (final concentration: 500 nM, 0.11 mCi) for 20 min at 37uC. Further treatment of the samples and HPLC separation was carried out as mentioned above. The conversion rate was calculated after analysis of the resulting chromatograms according to the following equation: |100. Figure S1 Energy profile drawn for the marmoset 17b-HSD1 model using PROSA. Energy profiles of marmoset 17b-HSD1 in complex with NADPH (orange) and of marmoset 17b-HSD1 in complex with NADPH and E1 (green). (TIF) Figure S2 Verify3D results for the marmoset 17b-HSD1 model. Verify-3D results are shown for the secondary complex of marmoset 17b-HSD1 (orange) with NADPH and for ternary complex (green) with NADPH and E1; residues with positive score are reasonably folded.

Supporting Information
(TIF) Table S1 Cluster analysis of molecular docking results.
All energies are expressed in kcal mol 21 . The lowest energy conformation of each cluster, which is marked in bold was used for further investigation. (DOC)