Computational exploration of the binding mode of heme‐dependent stimulators into the active catalytic domain of soluble guanylate cyclase

Soluble guanylate cyclase (sGC), the main target of nitric oxide (NO), has been proven to have a significant role in coronary artery disease, pulmonary hypertension, erectile dysfunction, and myocardial infarction. One of its agonists, BAY 41‐2272 (Riociguat), has been recently approved for treatment of pulmonary arterial hypertension (PHA), while some others are in clinical phases of development. However, the location of the binding sites for the two known types of agonists, heme‐dependent stimulators and heme‐independent activators, is a matter of debate, particularly for the first group where both a location on the regulatory (H‐NOX) and on the catalytic domain have been suggested by different authors. Here, we address its potential location on the catalytic domain, the unique well characterized at the structural level, by an “in silico” approach. Homology models of the catalytic domain of sGC in “inactive” or “active” conformations were constructed using the structure of previously described crystals of the catalytic domains of “inactive” sGCs (2WZ1, 3ET6) and of “active” adenylate cyclase (1CJU). Each model was submitted to six independent molecular dynamics simulations of about 1 μs. Docking of YC‐1, a classic heme‐dependent stimulator, to all frames of representative trajectories of “inactive” and “active” conformations, followed by calculation of absolute binding free energies with the linear interaction energy (LIE) method, revealed a potential high‐affinity binding site on the “active” structure. The site, located between the pseudo‐symmetric and the catalytic site just over the loop β2–β3, does not overlap with the forskolin binding site on adenylate cyclases. Proteins 2016; 84:1534–1548. © 2016 Wiley Periodicals, Inc.


INTRODUCTION
Soluble guanylate cyclase (sGC), the main target of nitric oxide (NO), plays an important role in several key biological processes such as vasodilation, Ca 21 cycling, endothelium permeability, myocardial contraction, long-term depression (LTD) and inflammation. Recently, several drugs have been found to stimulate sGC. One of them, BAY 41-2272 or Riociguat, has been recently approved for treatment of pulmonary arterial hypertension (PHA) and other related pathologies. Other agonists, as BAY 58-2667 or Cinaciguat, are in clinical phases of development. sGC catalyzes the cyclization of alpha phosphate in GTP to form cyclic GMP (cGMP). Vertebrate sGC is heterodimeric and two different isoforms have been described: a 1 /b 1 , predominantly found in the cardiovascular system, and a 2 /b 1 , found in brain. It is constituted by a regulatory domain (an H-NOX domain, located at the N-terminal of both subunits) and a catalytic domain (at the C-terminals) separated by a central region, containing H-NOXA (H-NOX-associated domain, that adopts a Per/Arnt/Sim or PAS fold) and a coiled-coil (CC) domain. 1 NO binds to a heme prostetic group located at the H-NOX domain causing a conformational change that affects the whole regulatory domain 3 and a 100-to 200-fold activation of the protein. 2 How the NO interaction with the heme-group results in an increase of the catalytic activity at a domain relatively far away is not known. NO seems to release heme distortion in the unbound enzyme, 3 what evokes a conformational change that might travel across the entire enzyme to reach the catalytic domain and produce an increase in the catalytic activity (as occurs probably in the activation of membrane guanylate cyclases). In the case of sGC, the activation signal could also reach the catalytic domain by a direct interaction between regulatory and catalytic domains. 4 Recent experimental data supports such direct contact between the two domains. 5 Independently of the mechanism of transmission from the regulatory domain, it also remains unclear what changes take place when the activation signal arrives into the catalytic domain. Comparison between "active" and "inactive" structures of the catalytic domain of adenylate cyclase and the "inactive" conformation of algae sGC 6 suggests that during activation both component subunits become closer. Since the catalytic site is localized just in the interface of both subunits, a and b, this approximation would also implicate a tighter GTP binding.
A great effort has been performed during the last years to find sGC modulators. Thanks to this work, several new drugs have been described to activate this enzyme. Some of these compounds [BAY 58-2667 (Cinaciguat), S3448, HMR1766 (Ataciguat), BAY 60-2770, and others] do not require the presence of the heme group (hemeindependent activators) to produce their stimulatory effect. 7 They are supposed to substitute the heme group in the regulatory domain and to force an active conformation of the enzyme. Compounds like CFM-1571, BAY 41-2272, BAY 41-8543, BAY 63-2521 (Riociguat), and others, have been developed using the classic agonist YC-1 as lead and require a native heme to exert their function (heme-dependent stimulators). 8 These compounds have a limited effect when alone, but in synergy with NO or CO they exert a potent sGC activation. Their effect has been associated predominantly to an increase of the maximal catalytic rate for the GTP cyclization reaction. 8 Despite this group of drugs is known for more than a decade, 9 the precise site of interaction with sGC is controversial.
Evolutionary old sGC presents two catalytic sites with mutual allosteric interactions. In mammalian sGC, some critical residues for the catalysis have been lost in one of these sites, although significant sequence similarity is maintained. This site is frequently named the pseudosymmetric site. Several authors consider this pseudosymmetric pocket as the most probable site of interaction of heme-dependent stimulators. 10,11 This thesis is based on two different classes of data, crystallization of the homologous enzyme adenylate cyclase 12 in the presence of its classical activator, forskolin, a compound with an effect very similar to heme-dependent stimulators on sGC, and mutational studies which have revealed that some residues close or forming part of the pseudosymmetric site are relevant to the function of YC-1. 10, 11 On the contrary, two other group of evidences argue in favor of an interaction of YC-1 with the regulatory domain, labeling with a photolabile derivative of BAY 41-2272 13 and data obtained using resonance Raman. 14 The absence of a reliable structure of sGC has precluded until now an optimal computational approach to analyze interaction of heme-dependent stimulators. However, recent studies have added valuable information of the catalytic domain. In 2008, the first structure of an eukaryotic sGC was published (3ET6, catalytic domain of sGC from Chlamydomonas reinhardtii). 6 Some years later, the X-ray structures of a non-physiological homodimer (b 1 b 1 ; PDB ID: 2WZ1 15 ) and heterodimers (a 1 b 1 ; PDB ID: 3UVJ, 15 4NI2 16 ) of human sGC have also been resolved. In both cases, structures of the catalytic domain were very similar to that described previously for adenylate cyclase. Thus, the relative conformations between the two subunits composing the dimer resembled adenylate cyclase in the "inactive" conformation. 6 In this study, we constructed in silico structures of the catalytic domain of mammalian sGC either in "inactive" (mimicking the global configuration of the 3ET6 structure) and "active" (based on the "active" catalytic domain of adenylate cyclase: 1CJU) 12 conformations. The quality of both homology models was analyzed throughout long molecular dynamics (MD) (approx. 1 ls each). Results suggest an important role of the loop b 2 -b 3 in the precise location of GTP for the catalytic reaction. Two Arg, 592 in the aand 539 in the b-chain, that interact with residues located both on the catalytic and on the pseudo-symmetric sites, seem to drive critical changes in the shape of this loop. Docking of YC-1 to all frames of selected simulations of the "inactive" and "active" models defined a relatively small area located in the center of the cleft as the potential target for agonist binding to the catalytic domain. The best pose of YC-1 was determined by free energy, and is located in the interface between catalytic and pseudo-symmetric sites and just over the loop b 2 -b 3 . This site, then, would be different from that previously proposed in previous studies that predicted the binding site on the same pseudosymmetric site in similarity with forskolin. Binding of agonists to the site defined in this study could stabilize a proper conformation of the loop and directly regulate the position of GTP inside the catalytic site. Interestingly binding of YC-1 to this area instead of the pseudosymmetric site would make compatible the simultaneous regulation of sGC by YC-1 and ATP as it has been suggested recently. 17

Modeling of the catalytic domain
Models of the catalytic domain of sGC in active and inactive states were constructed using Modeller 9v8. 18 This domain is multimeric and the active site is located just in the interface between the two monomers. In consequence, the precise location of the residues of one chain respect to those of the other one was considered critical for the recreation of a particular state or conformation. Then, specific crystal structures were used to sculpture each chain (chains A and B from 3ET6 for the a chain and chains A and B from 2WZ1 for the b chain) and in addition two dimeric templates were used to model each state. Alignment was done using 3D-Coffee. 19 For the active state of the domain, we used the crystal structure of adenylate cyclase with PDB code 1CJU in which the substrate L22 0 ,3 0dideoxyadenosine-5 0 -triphosphate (L-ddATP) has been substituted by GTP (rigidly superposed using PyMOL 20 ). Two water molecules located close to the substrate and the two atoms of Mg 21 were maintained in the chimeric template. For the inactive state, the dimeric structure 3ET6 in which GTP-Mg 21 and the two water molecules had been added was used as template. These molecules were superposed as heteroatoms with the STAMP program 21 using the two beta chains, of 1CJU and 3ET6, as references. The best of five models, as assessed by DOPE (Discrete Optimized Protein Energy) method, was selected for ulterior analysis. Models were also analyzed by the Verify3D Structure Evaluation Server, Prosa and PROCHECK. Ionization was assessed by three different programs Propka, PHEMTO and MOLARIS. Final ionization for each amino acid was decided by visual inspection and comparison of the different results.

Molecular dynamics
Antechamber was used to prepare models for MD using Amber Force field with the modification ff99SBildn. Models were solvated with TIP3P water molecules (15 Å as minimum distance to the edge of the box) and Cl 2 and K 1 ions added (to neutralize the protein and to increase the salt concentration of the solvent to 0.29 Osm). Minimization (40 ps) and thermalization (500 ps) of the protein were performed with NAMD 22 at CESCA and BSC supercomputing facilities (Barcelona). Molecular dynamic runs for structure equilibration of about 1 ls were performed using ACEMD 23 at a temperature of 310 K, timestep of 4 fs due to the use of the hydrogen mass repartition scheme 24 implemented in ACEMD, a cut-off distance for non-bonded interactions of 9 Å and Particle-Mesh-Ewald for long-range electrostatic interactions. Six different trajectories were simulated for each system. All figures and analysis of the protein dynamics analysis were done using the VMD program. 25

Docking and free energy calculations
The MD trajectories of MD-A 4 and MD-I 2 were sampled with frames extracted every one ns for the second half of the trajectory (that is, 500-1000 ns time frame). Protein structures corresponding to each frame were prepared for docking with the script prepare-receptor4.py of AutoDock-Tools. 26 GTP was included, but atoms of Mg 21 were excluded because Autodock Vina 27 does not recognize this atom type. The ligand (YC-1) was built with Maestro (Schr€ odinger, LLC) and prepared for docking with Auto-DockTools. Docking computations were performed using Autodock Vina, using a search box covering the whole catalytic domain (56.25 Å side), with the exhaustiveness parameter set to 20, and yielding 4 docked conformations per frame. The top-ranked docked conformations with an energy <28.6 kcal/mol were retained and clustered according to and RMSD criteria of 2.5 Å . The cluster representatives were then subject to binding affinities calculations using the linear interaction energy (LIE) method. 28 Here, binding affinities (DG calc bind ) are estimated on the basis of the difference (D) in average ligand-surrounding interaction energies hU l2s i extracted from the MD simulations of the ligand in the two relevant states considered for the binding process, that is in water and embedded in the active site of the solvated protein, respectively. The non-polar hU vdw l2s i and polar hU el l2s i interaction energies are treated separately, and their difference is scaled by different scaling factors, a and b: The constant a has been empirically set to 0.18 while b varies depending on the properties of the ligand. 28,29 In this study b was set to 0.37, which is the value determined for neutral ligands with <2 hydroxyl groups. 29 MD simulations for the LIE were performed using the program Q 30 with the OPLS-AA force field. 31 The parameters needed for the ligands that were not present in the original version of the force field were retrieved from the automatic parameterization performed with Macromodel (Schr€ odinger, LLC). Spherical boundary conditions were used, with a sphere of 25 Å radius centered on the center of mass of the ligand. The ionization states of titratable residues inside the simulation sphere were manually assessed with the following residues treated as ionized: Asp15, Asp59, Asp222, Glu137, Glu218, Glu 299, Arg103, Arg122, Arg284, Arg297, Lys135, Lys223. Ionizable residues close to the sphere boundary were modeled in their neutral form to account for dielectric screening. This sphere was solvated with TIP3P 32 water molecules and subject to polarization and radial constraints according to the surface constrained all-atom solvent (SCAAS) model 33 at the sphere surface, to mimic the properties of bulk water. Non-bonded interactions were calculated explicitly up to a 10 Å cutoff, except for the ligand atoms for which no cutoff was used. Beyond the cutoff, long-range electrostatics were treated with the local reaction field (LRF) multipole expansion method. 34 All solvent bonds and angles were treated with the SHAKE 35 algorithm and a 1 fs MD time step size was used. Protein atoms outside it were restrained to their initial positions and only interacted with the system through bonds, angles, and torsions. All simulations were carried out at a temperature of 298 K. Non-bonded pair lists were updated every 25 steps and the same interval was used for the sampling of the ligand-surrounding interaction energies. The equilibration phase lasted for 0.05 ns, during which all solute heavy atoms were restrained and gradually released as the temperature increased. The subsequent data collection phase lasted for 1 ns. The convergence was estimated by calculating the associated error as the difference between the average energies of the first and second half of the data collection phase.
In addition, the effect of the mutations described by Lamothe et al. 11 on the binding of YC-1 was evaluated as described below. Models of the mutated sGCs were created applying the VMD Mutator Plugin on the sGC structure resulting on the highest binding affinity (LIE calculations, previous paragraph). Mutated sGCs were solvated, minimized, thermalized and submitted to short MD of 10 ns (NAMD and Amber force-field). The MD trajectories were sampled with frames extracted every one ns for the second half of the trajectory (i.e., 5-10 ns time frame). In the first place, protein structures corresponding to each frame were prepared for docking with the script prepare-receptor4.py of AutoDockTools and docked with YC-1 using Autodock Vina. For each mutated sGC, the cluster with the best docking score was subjected to binding affinities calculations using the LIE method exactly as described previously. In the second place, the best YC-1 docking pose obtained for the native protein was superposed to the last frame of the MD trajectory of each mutated sGCs and subjected to the calculation of binding affinity using the LIE method.

Modeling of the catalytic domain
Our initial objective was to model the catalytic domain of mammalian sGC. This domain is formed by about 200 residues of the C-terminal of both chains (a 1 and b 1 ) sharing similar sequence (60% of homology), secondary (babbab) and tertiary structures (b-sheets forming the nucleus of each chain). Subunits present a head-to-tail association forming a groove at the interface. Catalytic and pseudo-symmetric sites are located inside this cleft. The goal was to use structural models to get a better understanding of the conformational changes that take place during enzyme activation and to analyze the potential binding of sGC stimulators to this domain. Thus, two different sets of models were constructed; on the one hand, an enzyme in a conformation as close as possible to the transition state and, on the other hand, a more relaxed or inactive structure. Then, not only good structural templates were needed for each protein chain (a 1 and b 1 ), but a suitable global structure of the dimer. Comparison of crystallized catalytic domains of guanylate and adenylate cyclases in "inactive" and "active" forms 6 suggest that during activation the a-chain twists respect to the b-chain what makes both subunits become closer at the catalytic cleft.
To construct the model, rat sGC sequence was chosen since a large quantity of experimental data was available. Rat a 1 and b 1 chains in this domain have 100% and 97% of identity with their human counterparts, respectively (see Fig. 1). At the moment of starting the project structural data of the whole mammalian heterodimer was not available. Only the homodimeric structures of human (2WZ1) 15 and C. reinhardtii (3ET6) 6 catalytic domains of sGC had been described. 2WZ1 is a nonphysiological b 1 /b 1 homodimer lacking a functional catalytic site (it does not contain the required catalytic residues present in the a-chain). Algae sGC contains two catalytic sites and is made of two identical subunits. Their sequences share a 53% and a 45% of identity (or a 66% and a 64% of similarity) with rat a 1 and b 1 chains, respectively. Conformation of the crystallized structure of this protein has been found to be very similar to "inactive" conformations of adenylate cyclase. 6 Then, the structure of 2WZ1 was used as template of the b 1 chain of our rat model and 3ET6 for the a 1 chain (Fig. 1). The final segment of this chain was not included in the final model because there was not apparent similarity with 3ET6 chains (Fig. 1). Alignment were performed taking into account secondary structure information by means of 3D-Coffee. 19 In addition, 3ET6 was used as template of the global configuration for the "inactive" model (see Material and Methods for a more detailed description). Since no "active" templates of sGC had been described, the 1CJU structure of adenylate cyclase 12 was used as template of the "active" conformation. Recently, structures of physiological heterodimers of the catalytic domain of human sGC have been published (PDB ID: 3UVJ, 15 4NI2 16 ). However, although they could have been considerably better templates for the individual sculpture of the a-chain, both have been reported as "inactive" conformations. Thus, the use of 1CJU as template to model the interaction between the two subunits would have also been required. 1CJU was also used as template for the inclusion of GTP [rigidly superposed onto L22 0 ,3 0 -dideoxyadenosine-5 0 -triphosphate (L-ddATP) in 1CJU], two atoms of Mg 21 and two water molecules within the catalytic site.
The best models of "inactive" (iGC 0 ) and "active" (aGC 0 ) structures created using Modeller, 18 as assessed by DOPE (Discrete Optimized Protein Energy) method, are shown in [ Fig. 2(A)]. The approximate location of the YC-1 Binding to Soluble Guanylate Cyclase pseudo-symmetric site and of the loop b 2 -b 3 , that will be referred to later, is also shown. As occurs between "inactive" 3ET6 and "active" 1CJU, 6 chain a in the "active" model or aGC 0 showed a shift with respect to the b-chain, showing some similarity to a very small turn of a screw. To quantify this displacement we defined an angle (h act , see Fig. 2, inset) between the a carbons of three amino acids: Met572 (a-chain), Met480 (b-chain), and Thr555 (b-chain). These residues are located in the nucleus of each subunit and form approximately a transverse plane to the vertical axis of the dimer (Note: in 1CJU, the first Met is substituted by an Iso). The value of h act is 7.9 lower in 1CJU than in 3ET6. iGC 0 (70.98) and aGC 0 (62.88) presented values of h act that were practically identical to their respective templates (Table I). Other measured parameters, like the distance between the catalytic residues of the a-chain and of the b-chain (d cat , distance dbetween the center of mass of a carbons of residues at both sides of the active center: Arg573, Asp485 and Asp529 in the a-chain and Lys593, Arg552, and Asn548 in the b-chain) and the minimal distance between a and b-chains at the cleft (d min , distance between the C carbon of Thr490 and the a carbon of Gly594; Fig. 2, inset) were also extremely similar (Table I). In [ Fig. 2(B)], the relative positions of the a carbons of the catalytic residues of both constructed models and templates are shown. The characteristic differences between "inactive" (3ET6) and "active" (1CJU) templates [e.g., a shorter separation of catalytic residues of 3D-Coffee alignment a 1 of and b 1 sequences of rat sGC (sequences of final models) with those of the PDB used as structural templates (3ET6, sGC of the eukaryotic algae C. reinhardtii; 1CJU, mammalian adenylate cyclase; 2WZ1, b 1 dimers of the human sGC). Colors represent the quality of the alignment. Blue rectangles enclosed: (A) the last part of the a 1 sequence that was not included in final models because of its low score and (B) residues implicated in GTP cyclization (M, in metal binding; P, in triphosphate-positioning; B, in nucleotide base discrimination; S, in riboseorienting). Amino acids that have been related with binding to YC-1 are labeled with asterisks.
a and b-chains and a displacement of Asp529 toward Asp485 in 1CJU; see yellow arrows in Fig. 2(B)] were maintained in the respective models, iGC 0 and aGC 0 . The location of GTP (transparent surface; no difference between iGC 0 and aGC 0 ) and Mg 21 (green spheres) after their rigid superposition onto L-ddATP/Mg 21 in 1CJU (CPK representation, ddATP; yellow spheres, Mg 21 ) is also represented [ Fig. 2(B)].

Molecular dynamics
Each model of sGC (aGC 0 and iGC 0 ) was submitted to 6 parallel MD (MD-A 1-6 and MD-I 1-6 , respectively) of about 1ls each (Table II). In addition, aGC 0 and iGC 0 models in which both GTP and Mg 21 had been removed were submitted to shorter simulations (0.2-0.3 ls; named MD-A w/o and MD-I w/o , respectively). Main structural changes (variations of RMSD of 2-3 Å ) took place during the equilibration period or during the first 1-5 ns of the MD. Thereafter, RMSD remained relatively stable throughout dynamics with standard deviations (SDs) below 0.5 Å , independently of the model ("active" or "inactive") and the presence or not of substrate [mean values are shown at Table II; see values at every frame in Figure 1 Table  II and Fig. 3(B)] showed some divergence between the different dynamics. In trajectories 2 and 4 of aGC 0 (MD-A 2 and MD-A 4 , see Table II) the angle become smaller than in the 1CJU template. Also, in some dynamics of the iGC 0 model (MD-I 4 and MD-I 6 ), h act was considerably reduced reaching values comparable to the "active" template. Similarly, in dynamics MD-A 2 and MD-A 4 , the value of d min reached values smaller than that found in 1CJU (Table II), and in MD-I 6 values close to that observed in "active" models. Then, in a great part of the simulations the conformations of the dimer lost some of the 1CJU-or 3ET6-like characteristics of the initial models. It allows the classification of the MD-A dynamics into two main groups: those in which the association between monomers become even tighter than in the original template (MD-A 2 and MD-A 4 ) and those that maintain or even loose that relationship (MD-A 1 , MD-A 3 , MD-A 5 , MD-A 6 ). A similar classification is difficult to do with MD-I dynamics because changes in h act and d min do not always go in the same direction.
Other of the differential features of the initial models, the smaller distance between catalytic residues of aand b-chains (d) in the "inactive" model, completely vanished in MD-I 1-6 trajectories. Interestingly, when iGC 0 is simulated in the absence of GTP-Mg 21 this difference persisted (MD-I w/o , Table II), suggesting that the presence of substrate limits the separation between the two chains.

Figure 2
Left, ribbon representations of "active" (red) and "inactive" (blue) models of rat sGC where subunits (transparent) have been superposed. GTP-Mg 21 is shown by a ball-and-stick representation. The location of the b 2b 3 loop and the pseudo-symmetric site is indicated by the corresponding labels, and the relative displacement of the a 1 respect to the b 1 subunit after "activation" (from iGC 0 to aGC 0 ) is represented by a green arrow. Right, comparison of the catalytic site in final models (iGC 0 and aGC 0 ) and "active" (1CJU, adenylate cyclase) and "inactive" (3ET6, 2WZ1) templates. GTP in the two models is displayed by a molecular surface representation and the compound used as template for GTP location in the catalytic site (ddATP, in PDB 1CJU) by a ball-and-stick representation.
Only the a carbons of the different catalytic residues are displayed. Numeration and names refer to those in the rat sGC sequence. Residues displacement during "activation" is represented by yellow arrows. Inset  a Interelationship between a and ß chains in constructed models of the catalytic domain of sGC (iGC 0 , "inactive" and aGC 0 , "active") and comparison with the crystal structures used as templates for the global conformation of the dimer in each case (3ET6, "inactive" algae guanylate cyclase, and 1CJU, "active" mammalian adenylate cyclase). Data for the recently described human structure 3UVJ is also included. As it is shown in Table II [see also Fig. 3(A) for specific trajectories], RMSD for GTP is considerably higher than that measured for the whole protein or the catalytic site. In addition, in average it is significantly bigger for trajectories of the inactive model (7.5 for MD-I 1-6 vs. 4.4 Å for MD-A 1-6 ). Fluctuation (SD of the RMSD), however, is smaller    in MD-I 1-6 than in MD-A 1-6 dynamics, suggesting that the position of GTP has suffered an abrupt change during the equilibration period of the iGC 0 model and that the new adopted location is relatively stable [see Fig. 3

(A) for comparison of MD-A 4 and MD-I 2 ]. Fluctuation of RMSD for GTP in certain MD-A dynamics was impressive [see for example the RMSD for GTP in MD-A 6 , Fig. 3(A)]
and higher in trajectories in which the conformation of the dimers was far away from the 1CJU structure (Table II). Visual inspection of the different dynamics showed that the g-phosphate group remained relatively still during the fluctuations of GTP (RMSD of 2.2 6 0.9 for MD-I 1-6 vs. 3.3 6 0.5 for MD-A 1-6 ). In addition, lateral movements are very limited. For that reason, fluctuation of GTP mainly consists of a twist of the less charged segments of GTP toward the pseudo-symmetric site and can be defined at an important extent by its vertical inclination angle. We defined it as the angle x GTP formed by Val528 (in the core of the b-chain) and g phosphorus and carbon C2 (base segment) of GTP (Fig. 4, inset). We used Val528, a residue located in a plane lower than the catalytic site, to avoid negative angles in some of the simulations. Then, in 1CJU and in the initial models, aGC 0 and iGC 0 , the angle had values of about 308 (378 for 1CJU and 378 and 358 for aGC 0 and iGC 0 , respectively). In Figure 4(A), different positions of GTP in a typical trajectory of the "inactive" model (MD-I 1 ) and in one of the trajectories of the "active" model that maintained a 1CJU-like structure (MD-A 4 ) are shown. In MD-I 1 , GTP had an inclination close to a right angle and completely out of the catalytic bed. Very similar values could be found for the other dynamics of this model [ Fig. 4(C)]. As it is shown for MD-A 2 , this position was reached during the equilibration period and did not change significantly throughout trajectory [ Fig. 4(B)]. In MD-A 4 , GTP maintained values of GTP similar to that measured for 1CJU through the whole trajectory [ Fig.  4(B)]. On the contrary, in dynamics MD-A 3 , MD-A 6 and perhaps MD-A 5 [ Fig. 4(C)], GTP oscillates continuously throughout simulation [see x GTP values for MD-A 6 in Fig. 4(B)].

Figure 4
Left, different snap-shots of GTP (stick representation) in representative dynamics of the "active" (MD-A 4 ) and "inactive" (MD-I 1 ) models. Amino acids of the catalytic site, that have been superposed to be used as reference of GTP movements, are displayed by a molecular surface representation. Inset, visualization of the angle formed by the a carbon of Val528 (in the core of the b-chain), and g phosphorus and C2 carbon (base segment) of GTP; x GTP , that is used to define the displacement of GTP respect to the catalytic site. Right-top, value of x GTP for representative dynamics (MD-A 4 and MD-A 6 , "active"; MD-I 2 , "inactive") at every frame of the simulations. Right-bottom, mean values of x GTP and of extra-coil width (distance between Ca of Val524 and Iso527, see Figure 6A to visualize the extra-coil in the loop b 2 -b 3 ) for all dynamics (the red horizontal line corresponds to the values found in the "active" template 1CJU). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] To understand the different behavior of GTP between MD-A and MD-I simulations and between the different MD-A trajectories, we analyzed in detail the conformation of the protein in the catalytic site and in its surroundings. We observed that the smallest width of the catalytic site was between loop (Gly528) and the catalytic residue Asn548. Interestingly that width (d width ) was of 8.7 Å for MD-A and only of 7.6 for MD-I trajectories (Table II). In fact, in those MD-A trajectories where GTP showed a higher degree of oscillation the value of d width were lower than 8 Å (MD-A 3 and MD-A 6 ). However, some MD-I trajectories had values close to 10 Å and, in spite of that, GTP was expelled from the catalytic bed. The explanation seems to be related with the position of the loop b 2 -b 3 (where the catalytic residue Asp529 is located). We commented two differences in the distribution of the catalytic amino acids between 1CJU-and 3ET6-like structures. One was the mean distance between band a-chain catalytic residues, difference that vanished in the presence of GTP at the beginning of the dynamics, and the other was the displacement of Asp529 toward Asp485 in 1CJU-like structures [ Fig. 2(B)]. The degree of this displacement, that we monitored measuring the separation between Asp529 and Cys541 (d cys , see Table II) Two MD, MD-A 2 and MD-A 4 , maintained values of h act , d cat , and d min similar to those of the "active" template of adenylate cyclase. In the MD-A 4 trajectory, in addition, the position of the lateral chains of residues Arg552 and Asn548 in relation with GTP were extremely close to 1CJU. For that reason, this dynamics was chosen to study the interactions that maintain the configuration of the catalytic site in an "active" state and in particular the conformation of the loop b 2 -b 3 . In spite of the stability of the general structure of the catalytic site, moderate changes took place in the position of the lateral chains of the catalytic residues during the first 450 ns of the dynamics. Then the structure acquired a stable  Distances between the base or P a of GTP to the closest residues of the catalytic site (Arg552 and Cys541, respectively) were extremely constant [ Fig. 6(C)]. This period of the dynamics was also characterized by a steady curvature of the extra-coil in the loop b 2 -b 3 [ Fig. 6(D)]. Thr526, the residue located at the crest of the curvature, was very close to Arg529 [ Fig. 6(E)] and to the catalytic residue Glu473. In fact, the frequent hydrogen bonds formed between these amino acids seemed to be the main responsible for the stabilization of this structure [ Fig. 6(E)]. After 650 ns, these bonds started to be less common allowing interactions of Arg592 with the equivalent Arg in the other chain (Arg289). At 907 ns, Arg592 definitively got departed from Thr526 and bound to the previous residue in the loop b 2b 3 , Glu525 [ Fig. 6(E)], about 3.9 Å further on [ Fig. 6(E)]. Previously, Glu525 had been bound by a hydrogen bond to Met537. The consequence of these changes was a relaxation of the extra coil, the progressive extension of the loop toward GTP [ Fig. 6(B)] that finally put pressure on it. It increased the GTP folding, pulled it around 2-3 Å toward the opening of the catalytic cleft and removed its close contact with the catalytic residue Arg552 [ Fig. 6(C)]. In Figure  7, there is a simplified representation of the changes in the pattern of hydrogen bonds observed in MD-A 4 . In those dynamics of the aGC 0 model more distant to the 1CJU structure (MD-A 3 and MD-A 6 ), the interaction between Arg592 and Glu525 was relevant since the beginning (30 6 50% and 57 6 68% of hydrogen bond, respectively).  Then, Arg592 seemed to act in this system as a switch of a hypothetical change from an "active" to an "inactive" state.
Interestingly, Arg592 belongs to the loop b 4 -b 5 that forms part of the pseudo-symmetric site. Cys594 (a-chain), two residues apart from Arg592, is the counterpart of Cys541 (b-chain) in the functional catalytic site. Winger at al. 6 propose that in the sGC of C. reinhardtii communication between both catalytic sites is mediated by direct interactions between loops b 2 -b 3 of each monomer. However, in the MD of our models we only found direct hydrogen bonds between residues of these two loops on rare occasions. Interactions seemed to be indirect and mediated by Arg529 and Arg289 (Fig. 7, inset). Arg529 can make hydrogen bonds simultaneously with Thr526 (loop b 2 -b 3 , a-chain) and Glu473 (loop b 2 0 -b 3 0 , b-chain). Arg289 makes contact with Glu473 (b 2 0 -b 3 0 ) and Pro591 (loop b 4 -b 5 , a-chain).

Docking and free energy calculations
One of the critical questions about the biochemistry of sGC is the precise site of interaction of the different agonists recently discovered. No doubt exist that hemeindependent stimulator bind to the regulatory domain, probably by the substitution of heme itself. Some models of the potential binding site for these compounds have been reported. 10,11 For heme-dependent stimulators, on the contrary, two different hypotheses exist. Raman and labeling studies with a modified heme-dependent stimulator suggest binding of these drugs to the same domain. 13,14 However, a latter study found that the response to NO and to YC-1 was not altered in sGC in which the segment that contains labeled residues had been removed. 36 These authors propose another location on the N-terminal region as the binding site for this molecule. The other hypothesis is that YC-1 binds to the  catalytic domain. This is based on studies that mutate residues in sGC equivalent to those that had been found to be relevant for binding of forskolin to adenylate cyclase. Additional studies analyzing the effect of different mutations suggest that YC-1, and by analogy the other heme-dependent drugs, binds at the catalytic domain near the pseudo-symmetric site. 11 Depending on the mutated amino acid, authors find blockage of its direct stimulatory effect (in the absence of NO) or of its synergic effect with NO. All these studies measuring cGMP synthesis in the presence or absence of stimulators have the limitation that mutations can be affecting the transmission of the effect and not the binding itself. Up to our knowledge, the binding of YC-1 or other stimulators has not been directly measured in any study in which specific residues had been mutated.
A computational approach had been difficult until now because of the absence of a reliable structural model of mammalian sGC, although a recent study has used an homology model of the regulatory domain based on prokaryotic structures to perform docking studies and suggest a potential binding site near the heme-group. 37 Now, the structures of human regulatory and catalytic domains had been reported, but at the beginning of this project only sGC structures of the catalytic domain had been describd for eukaryotic cells. Then, we only studied here the potential binding of YC-1 to this domain. With this objective, we performed an exhaustive docking study of YC-1 to every frame of equilibrated MD-A 4 and MD-I 2 trajectories with Autodock Vina. 20 In total, 4 poses were retained for each of the 5,000 frames considered from each trajectory. Figure 8 shows that docking poses were mostly concentrated in one area of the domain, the cleft formed at the interface between both monomers. In MD-A 4 , docking poses occupied an area that starts just over GTP and extend until the pseudo-symmetric site, including the highest ranked poses. In addition, some poses were located on a cavity that begins at this pseudo-symmetric site and that in some structures almost went through the entire depth of the protein. In MD-I 2 , because of the position of GTP toward the pseudo-symmetric site, most of the poses found in MD-A 4 were impeded. Instead, most populated area lay under GTP, where the substrate in "active" models is located. These poses presented generally lower scores than those corresponding to the highest populated area from the MD-A 4 conformations. Docking of YC-1 to the initial models (aGC 0 and iGC 0 ) or 1CJU gave poses concentrated on a small groove between the a 1 helix and the loop formed at the end of the a 3 helix or just on the forskolin binding site, respectively (Fig. 3

of the Supplementary Material).
The best docking poses for each frame of the MD-A 4 trajectory were considered for further analysis, following a filtering according to a minimum scoring threshold (28.6 kcal/mol) and clustered using a 2.5 Å RMSD cutoff, resulting in 10 clusters. A total of 18 poses (i.e., between 1 and 6 from each cluster) were retained after visual inspection. Eight representative poses with a score above the initial threshold were additionally considered for the next stage, to account for different positions of the cavity. It followed a rescoring of all the 26 docking Docking of YC-1 to every frame (each 0.1 ns, since 500 ns) of dynamics MD-A 4 (A) and MD-I 2 (B). The best pose of YC-1 for each snapshot is shown by its center of mass (red points). Structures of GTP and a 1 subunit at 500 ns are displayed by surface representation. Black circle shows the area where the poses with the best predicted DG bind,LIE are located (see Table III and Figure 9). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] YC-1 Binding to Soluble Guanylate Cyclase poses using the LIE method. 28 At this stage, we discarded the docking results from the MD-I 2 trajectory for further analysis, since all of them laid on the GTP catalytic site and therefore could not explain the stimulatory effect of YC-1.
From the LIE analysis, two binding modes (clusters 1 7 and 5 2 ) clearly emerged as the most energetically favored (Table III). We will further refer to the best ranked 1 7 . It lies between loops b 2 -b 3 and b 2 0 -b 3 0 , a location just in the interface between catalytic and pseudo-symmetric sites. The polar group of YC-1 strongly interacts with first and second atoms of Mg 21 . The rest of the drug molecule lays horizontally over the b 2 -b 3 loop [ Fig. 9(A)]. Main interactions (apart from that with Mg 21 atoms) are with residues of the two loops previously mentioned, the helices a 2 and a 2 0 , and the strands b 1 and b 1 ' (Table IV). This proposed binding site for YC-1 does not correspond with that previously suggested on the pseudo-symmetric site. 11 Docking of other heme-dependent drugs (Bay 41-2272, Bay 41-8543, Bay 63-2521) and less studied stimulators (A350619, A778935, Benzydamine, CFM-1571) to the protein structure used for cluster 1 7 resulted in poses very similar to that observed for YC-1 [see Fig. 9(B) for the docking result with best score for Bay 63-2521, and Table I  Of the four mutations described to modify YC-1 response, 11 three of the substituted residues are surrounding the YC-1 binding site proposed in this study (Fig. 5 of Supplementary Material). Models of the catalytic domain with the different mutations analyzed experimentally (a 1 C594D/b 1 , a 1 C594Y/b 1 , a 1 C594D-E525K/b 1 , a 1 /b 1D 477A, and a 1 / b 1 M537N) were constructed, as described in Materials and Methods, and submitted to short MDs (10 ns). Five different structures were extracted from the MDs for each mutated model and used to perform docking with YC-1. Although the short time of the MD precluded big structural changes of the mutated proteins, none of the docking poses obtained laid close to the pose with the highest affinity for native sGC (pose in cluster 1_7). Most of them were located on the pseudo-symmetric site and their affinities (calculated by LIE) were quite low in comparison with the best poses previously achieved for native sGC (Table V). On the other hand, we directly superposed the pose of YC-1 in cluster 1_7 of the native protein onto the different mutated sGC to evaluate the resulting affinity by the LIE method (after the corresponding periods of equilibration, as described in Materials and Methods for LIE calculations) ( Table V). Substitutions a 1 C594D/ b 1 and a 1 C594Y/b 1 reduced the binding energy by 1.6 and 3.2 kcal/mol, respectively (what would make 5 and 24 times bigger the dissociation constant). Furthermore, the double substitution a 1 C594D-E525K/b 1 , that suppress the response to NO, YC-1 and the synergy between both in the study of   Similarly, the described subtitutions of residues of the b chain, a 1 /b 1D 477A and a 1 /b 1 M537N, also resulted in significantly lower affinities for YC-1 in this location. On the other hand, E525 and M537 have a key role on the conformational changes of the catalytic domain described above (Fig. 7). Our hypothesis is that these structural changes might be implicated on the modulation of enzyme activity (both on direct stimulatory responses or on synergic effects) and thus any substitution of these residues, apart from the potential consequences on drug binding, would have an effect on the transmission of the YC-1 response.
Although there are no direct evidences demonstrating nucleotide binding to the pseudo-symmetric site, it is known that mammalian sGC can bind two molecules of substrate. On the other hand, sGC can use (at least in vitro) ATP as substrate and ATP is a known inhibitor of cGMP synthesis. The nucleotide regulatory site must be found on the catalytic domain because ATP also inhibits enzymatic activity on dimers containing only this part of the protein. 17 Interestingly, YC-1 effect is not influenced by the presence of ATP in full-length sGC, suggesting that YC-1 does not share the same binding site of ATP. Direct stimulatory effect of YC-1 (activity measured in the presence of Mn 21 ) was absent on the catalytic domain alone. These data suggest that YC-1 require a native regulatory domain to stimulate sGC in the absence of NO. This could mean that the binding site of YC-1 is in that domain or that the heme-domain has a significant influence in the transmission of its activation signal. It is difficult to extrapolate more conclusive data from this result because mutations of residues in the catalytic domain (such as the above commented mutation C594Y) can suppress the response to NO or YC-1 alone, but not the effect when they are together. 11 From the mutational data, it seems that there are two different signaling effects that crosses the sGC and converge finally on the catalytic site. One could come from the heme-group after binding of NO and the other from the pseudo-symmetric site after binding of a hypothetical allosteric modulator (e.g., ATP). YC-1 could be located on the middle of the transmission pathways of both activatory signals, then acting both as a direct stimulator and as a synergic agent.
Our study cannot rule out binding of YC-1 to the regulatory domain, but propose a potential binding site in the catalytic domain with a considerable binding energy that could explain some of the mutational studies published previously. On the contrary to previous proposals, the binding site of heme-dependent stimulators at the catalytic domain would not be located on the pseudo-symmetric site, but in the interface of catalytic and pseudo-symmetric sites. This  location would be compatible with the simultaneous binding of ATP or GTP to the pseudo-symmetric site to exert allosteric modulation.