The molecular models of CRF-BP, CRF_2αR, and CRF_2βR were constructed using MODELER (29, 30), as implemented in the Protein Modeling module of Discovery Studio v2.1 (Accelrys Inc., San Diego, CA, USA). Human CRF-BP, CRF_2αR, and CRF_2βR reference sequences were retrieved from the Uniprot database, with accession numbers P24387, Q13324-1, and Q13324-2, respectively (31). CRF was modeled using the crystal structure of human CRF inactive analog (PDB: 1GO9) containing a D-Phe residue at position 12 and alpha-aminoisobutyric acid in position 15 (32).

For CRF-BP, top scoring models produced by threading-based approaches identified by Muster and Phyre2 servers (33, 34) were retrieved, aligned, and used as starting templates to generate a human CRF-BP model. Fragments from gastric intrinsic factor receptor cubilin (PDB: 3KQ4) (35) and neuropilin (PDB: 2QQL) (36) were used to construct the model. Secondary structure elements restraints such as α-helices and β-sheets as predicted by PCI-SS server (37) were included, as well as experimentally determined disulfide bridges (38) during modeling (Figure S1A in Supplementary Material).

For CRF receptors modeling, we used the crystal structure of CRF₁R (PDB: 4K5Y) (27), the N-terminal ECD of human CRF_2αR in complex with UCNI (PDB: 3N96) (23), and murine CRF_2βR in complex with Astressin analog peptide (PDB: 2JND) (24) as templates. In addition, the crystal structure of the transmembrane bundle of glucagon receptor (PDB: 4L6R) was used as guide to model the N and C-terminal portions absent from the available CRF₁R crystal structure (26) (Figures S1B,C in Supplementary Material). For each protein model, a set of 100 models were constructed and the best model according to Modeler internal PDF score was subjected to a molecular minimization protocol using the CHARMM22 force field available within Discovery Studio (39, 40). The protocol consisted of 5,000 steps of steepest descent method, followed by 10,000 steps of conjugate gradient method to reach a final root-mean-square (RMS) gradient of 0.001 kcal/mol/Ų.

The overall quality of the final models was assessed by Ramachandran plot using the RAMPAGE server and quality model assessment with ProSA (protein structure analysis) server, respectively (41, 42). ProSA web was used to check and compare the obtained protein structural models with those experimentally determined by X-ray crystallography or NMR (42). The ProSA z-score indicates overall model quality and measures the deviation of the total energy of the structure with respect to an energy distribution derived from random conformations. The APBS software was used to calculate the spatial distribution of electrostatic potential on protein atoms using a two-dielectric implicit solvent model and the finite difference method to solve the Poisson–Boltzmann Equation (43). The dielectric constant used was 4 for proteins and 80 for the solvent.

The CRF₂Rs and the CRF_2αR/CRF-BP complex were inserted into a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine POPC lipid membrane considering the spatial arrangements of the protein with respect to the hydrocarbon core of the lipid bilayer, as obtained from the OPM database (44). For the CRF₂R systems, a 150 Å × 150 Å × 120 Å box consisting of the protein, lipids, classic TIP3P model for water molecules, and 150 mM KCl was generated using the membrane builder module of CHARMM-GUI (45, 46). In a similar fashion, the CRF_2αR/CRF-BP complex was embedded in a 140 Å × 140 Å × 160 Å box. MD simulations were carried out with the NAMD 2.9 simulation package (47), using the CHARMM36 force field parameters for proteins and lipids (48, 49). Periodic boundary conditions were imposed in all three directions and the Particle Mesh Ewald method was used to account for full long-range electrostatic interactions within the selected boundary condition within a relative tolerance of 1 × 10⁻⁶ (50). The final systems were composed of nearly 235,000 atoms for CRF₂Rs, and nearly 255,000 atoms for the CRF_2αR/CRF-BP complex. The simulations were started from different seeds, and three replicas of 100 ns for each CRF₂ receptors were performed, while a single 100 ns simulation was performed for the CRF_2αR/CRF-BP complex. A 12 Å cutoff was used to compute non-bonded interactions with a smooth switching function applied at a distance of 10 Å. To impose the thermal exchange with an external thermostat, the isobaric–isothermal ensemble (NPT) with constant number of particles N, pressure P, and temperature T was used. Constant temperature was maintained by coupling the system to a thermal bath whose temperature is maintained via Langevin dynamics with a friction coefficient of 1 ps⁻¹. Constant pressure was maintained using a Langevin piston at a nominal value of 1 atm (51). The SHAKE algorithm, with a tolerance of 1 × 10⁻⁸ Å, was applied to constrain the length of all covalent bonds involving hydrogen, thus allowing the use of a 2 fs integration time step along with the r-RESPA integrator, which allows a multiple time step scheme where bonded, short-range non-bonded, and long-range electrostatic terms are calculated every 2, 2, and 4 fs, respectively. By plotting Cα-root-mean-square deviation (RMSD) and RMS fluctuation (RMSF) along the MD simulation, we assessed the structural equilibration reached by our models. To further characterize the structure of the CRF₂Rs, three parameters were calculated: angle phi, defined between the hinge region (connecting the TM with the ECD) and the center of mass (COM) of the TM domain; angle theta, defined between the hinge region and the COM of the ECD; and the distance between the COM of both domains. Both angles helped define the orientation of the ECD with respect to the XY plane (parallel to the membrane plane) and the Z axis (perpendicular to the membrane) (52). Also, for the CRF_2αR/CRF-BP complex, the total internal energy of the complex was calculated, as well as the total interaction energy between the CRF_2αR and the CRF-BP in terms of its electrostatic and van der Waals components. These calculations were performed using the NAMD Energy analysis tool available in the Visual Molecular Dynamics v1.9.3 (VMD) software (53).

Protein–protein docking was performed using Hex v8.0 with default parameters (54). Briefly, for the generation of the top scoring solutions, we used an initial Steric Scan at N = 16, followed by a Final Search at N = 25, obtained by using just the steric contribution to the docking energy. We used the Shape only correlations, the 3D Fast Lite as FFT mode, with a grid dimension of 0.6 Å. These orientations are sorted by calculated energy, and a new set of trial orientations are generated for the top scoring. 10,000–20,000 orientations using the Scan Step and SubSteps were used to construct new distance samples in steps of ±(Scan Step 0.75 Å)/(Substeps 2) from the initial orientations, 1 Å resolution was used to scan the search space and a 0.5 Å resolution was used to perform the high-resolution scoring (55). A final minimization protocol for the top scoring solution complexes consisted of 20,000 steps of steepest descent method, followed by 10,000 steps of conjugate gradient method to reach a final RMS gradient of 0.001 kcal/mol/Ų to obtain the final models. Protein interactions such as disulfide bonds, hydrophobic interactions, ionic interactions, hydrogen bonds, aromatic–aromatic interactions, aromatic–sulfur interactions, and cation–π interactions within a protein or between proteins in a complex were calculated using the PPI server (56).

We have previously demonstrated that CRF-BP and CRF_2αR interact (20). In order to further characterize this interaction and predict which residues are forming the interacting interface, we reasoned that the prediction of structural models for CRF-BP and CRF₂R were necessary. There are no crystal structures or modeling for CRF-BP. As CRF-BP sequence is conserved among species but displays no significant sequence similarity to any other known protein experimentally resolved; a threading approach was used to predict a model of the structure of CRF-BP. Fragments from neuropilin (PDB: 2QQL) (36) and gastric intrinsic factor receptor cubilin (PDB: 3KQ4) (35) were used as starting templates. The predicted structural model for CRF-BP fold consisted in two modules. The first module containing residues 50–180 which displayed two short alpha helices, six antiparallel beta-sheets, and the first pair of disulfide bridges C60–C81 and C104–C141. The second module comprised residues 180–245 of the protein, with a pair of beta-sheets and one alpha helix that contained a second pair of disulfide bridges C183–C205 and C237–C264 (Figure 1A). An electrostatic potential surface (EPS) was obtained for CRF-BP. In the protein, two acidic patches were observed (Figure 1B left, red colored) and when the protein was turned in 180°, two basic patches were observed (Figure 1B right, blue colored). In order to validate the predicted model, a Ramachandran plot distribution and a ProSA protein quality analysis were performed. Ramachandran statistics showed that more than 95% of the residues of the predicted model were in the allowed geometric regions for amino acids (Table 1; Figure S2A in Supplementary Material). Although some amino acids were positioned in the non-allowed regions, they were residues participating in protein turns. This result indicates that the obtained fold is feasible. ProSA protein quality analysis casted out a Z-score = −3.46 (Table 1), value that falls in the range of native structures (Figure S2B in Supplementary Material), indicating that a good quality protein model was predicted. The computational engine used for the calculation of z-score and plots uses knowledge-based potentials of mean force to evaluate model accuracy (57). The potentials of mean force compiled from the PDB database provide a statistical average over the known structures. A z-score within the range characteristic for native proteins is indicative of a correct structure (42).

CRF-BP is a protein that binds CRF and UCNI with high affinity, and these interactions have been well characterized (58, 59). To further validate the CRF-BP structural model, protein–protein docking experiments between CRF-BP with CRF and UCNI were performed and the predicted residues involved in this interaction were obtained and compared with previously published data. The obtained binding modes predicted that the C-terminal domain of CRF and UCNI may bind over the positively charged surface at the N-terminal domain of CRF-BP (Figure 2A). The CRF-BP/CRF interaction comprises mainly CRF-BP residues R55, R56, C60, L61, D62, M63, L64, T71, F72, and T73, and CRF residuesV18, M21, A22, E25, Q26, A28, and Q29 (Figure 2B). In addition, the CRF-BP/UCNI interaction comprises mainly CRF-BP E51, R55, R56, C60, L61, D62, M63, L64, S65, I86, and W116, and UCNI L18, L21, A22, S26, E29, E32, Q33, N34, I36, and D39 (Figure 2C). This obtained binding mode was in agreement with site-directed mutagenesis data from CRF-BP (12). Furthermore, it has been previously described that CRF-BP binds CRF as a dimer (59). Therefore, we also performed protein–protein docking experiments for two CRF-BP alone and with CRF in order to further validate our CRF-BP structural model. The results showed that the CRF-BP model was permissive for a symmetrical homodimerization arrangement (Figure 3A) and for interacting with CRF as a dimer (Figure 3B). Thus, all the aforementioned validation approach suggest that the predicted fold is feasible, of good quality, and in agreement to previously published data (12, 60). In addition, the obtained binding model for (CRF-BP/CRF)₂ predicted that residues 19–38 of CRF are sandwiched by the CRF-BPs. In addition to already described interaction between CRF and a CRF-BP monomer, CRF may also interact with two patches within the C-terminal domain of CRF-BP (Figure 3C). CRF-BP residues 235–238 and 257–264 contacts CRF, with prediction of the residue E238 from CRF-BP displaying H-bond interactions with R23 and Q26 of CRF. CRF-BP D262 main chain carbonyl group also contact R23 of CRF. An additional H-bond interaction is predicted to occur between the side chain of T258 from CRF-BP and the main chain NH group of N34 from CRF.

There are still no full-length CRF receptor crystal structures available (21, 22). For CRF_2αR and CRF_2βR, only the ECD structures obtained by NMR are available, human CRF_2αR-ECD in complex with UCNI, and murine CRF_2βR-ECD in complex with an Astressin analog peptide (PDBs: 3N96 and 2JND, respectively) (23, 24). Thus, the recently solved crystal structure of CRF₁R, in addition to the CRF_2αR and CRF_2βR ECDs available structures, were used as templates to predict a model for the structure of CRF₂Rs (Figure 4A). The crystal structure of the transmembrane bundle from glucagon receptor (PDB: 4L6R) was used to model the extended helix 1 in the N-terminal region (TM1stalk region), the intracellular loop 2 (IC2), and the helix 8 in the C-terminal region, which are absent from the CRF₁R crystal structure (21, 61). Considering the already available structural information and guided by similar works in class B GPCRs, the N-terminal domain of CRF was located in a position able to interact with the J-domain on the CRF₂R TM bundle (62, 63). The EPS obtained for CRF₂Rs showed that the electrostatic potential is similar for both receptors, with minor differences in the N-terminal region. The CRF_2βR (Figure 4C) showed a more extensive basic patch than CRF_2αR (Figure 4B).

Molecular dynamics were performed in order to test the stability of the CRF₂R models including the ECD and TM regions. The receptors were embedded in a pre-equilibrated POPC lipid bilayer and solvated using the Membrane Builder in the CHARMM-GUI web server (Figure S3 in Supplementary Material). Each system was subjected to a 100 ns of MD simulations with three replicas. The RMSD and RMSF were computed over the course of the simulation for the Cα atoms of the proteins to measure structural stability and qualitatively characterize the dynamics of the proteins (Figure 5). The CRF_2αR and CRF_2βR MD trajectory analyses showed no significant changes in the RMSD values for the ECD (red lines) and TM (blue lines) regions (Figures 5A,D). However, the RMSD values calculated for the full-length receptors showed a significant change for CRF_2αR and CRF_2βR (Figures 5A,D, black lines). The changes obtained for CRF_2αR can be attributed to translational and rotational movements of the ECD relative to the TM domain. The values obtained are coincident with three main different conformational states for CRF_2αR: open-like, semi-closed, and closed-like (Figure 5B). The closed-like state was the one obtained at the end of the simulation, indicating that the receptor has a higher tendency for that conformation. In the case of CRF_2βR, the values obtained are coincident with only one open-like conformational state, which varies the angle of extension (Figure 5E). RMSF describes the average fluctuation of each Cα atom of the amino acid residues in the proteins over the simulation time (Figures 5C,F). The general fluctuations of specific regions of the proteins are similar for both CRF₂R. Within the TM region, the peaks of higher movement are coincident with the intra and extracellular loops. Moreover, in the ECD region, the peaks with higher movement are coincident with the loop that connects the α-helix with the β sheet bundle, and with the loop that connects the ECD with the TM region (stalk region).

Both CRF_2αR and CRF_2βR showed great variation in phi and theta angles (Figure 6). Average theta angles for CRF_2αR are near to the 40°–50° range, while for CRF_2βR the average value is close to 20° ± 10° (Figure 6B). According to the data described for the GCGR, an angle of less than 20° corresponds to the closed state of the receptor while values close to 40° are associated to an opened state (52). Upon analyzing the phi angle, CRF_2αR reaches an average value of 65°± 20°, while CRF_2β R stays at higher values at 85° ± 10° (Figure 6D). Analog to angle theta, low values of angle phi (~20°) have been associated with the closed conformation, which would mean that our simulations are either in the open conformation or in a semi-closed conformation. This is further supported by the distance between ECD and TM COMs, where both receptors reach similar values in the range of 55–60 Å (Figure 6C and Figure S4 in Supplementary Material), and this distance was associated with an open conformation. The results suggest that CRF_2βR is more stable and displays mainly only one conformational state and that CRF_2αR is less stable, reflected by more fluctuations within the Cα, and it has three main conformational states. In addition, in both receptors, the ECD region is the one with more fluctuations.

We previously demonstrated that CRF-BP physically interacts with the ECD of the α but not with the β CRF₂R isoform (20). In order to test our models, a protein–protein docking was performed to predict the potential binding mode for CRF-BP with the ECD region of CRF₂Rs. For the CRF-BP/CRF_2αR interaction, characterization (20) immunofluorescence co-localization analyses using the Santa Cruz Biotechnology N-20 anti CRF₂R antibody (residues W27-Q46 for CRF_2αR and I53-Q73 for CRF_2βR) were performed. We reasoned that, if the residues recognized by the antibody are available to bind the antibody they should not be participating on the CRF-BP/CRF_2αR interaction, thus, these residues were excluded from the search space during the protein–protein docking assay. The best solution obtained for CRF-BP/CRF_2αR docking casted out a total energy value of −838.1 kJ/mol and for CRF-BP/CRF_2βR −685.9 kJ/mol (Table 2, Figure S5 in Supplementary Material). Both protein have similar number of residues (382 and 372) and molecular mass (44.53 and 43.74 kDa); therefore, the estimated binding energy values suggest that CRF-BP could bind to CRF_2αR and form a more stable complex compared to CRF_2βR.

Molecular dynamics simulations were performed to test the stability of the predicted CRF-BP-CRF_2αR model and qualitatively characterize the dynamics of the complex (Figure 7). The CRF-BP-CRF_2αR complex MD trajectory analyses showed no significant changes in the RMSD values for the CRF-BP N-terminal domain (blue lines). However, the RMSD values calculated for the full-length complex, CRF-BP and CRF-BP C-terminal domain display a significant variation (Figure 7A green, black, and red lines, respectively). The changes can be attributed to translational and rotational movements of the ECD relative to the TM domain and particularly the C-terminal domain of CRF-BP, in agreement to the higher RMSF of the C-terminal domain of CRF-BP (Figure 7B). The energy of interaction of the CRF-BP-CRF_2αR complex during the dynamics indicates that the complex is stable and that the main contribution comes from electrostatics rather than van der Waals interactions (Figure 7C). This phenomenon was observed along with loss of connections between C-terminal domain of CRF-BP and CRF_2αR, as shown by enhanced flexibility with respect to starting conformation. The interaction between CRF-BP and CRF_2αR, the amino acids present in the interacting interface previously determined using the Protein interaction server (56), were also measured during the CRF-BP-CRF_2αR complex dynamics. For CRF-BP/CRF_2αR interaction interface, hydrophobic interactions (Figure 8A), ionic interactions within 6 Å (Figure 8B), cation–π interaction within 6 Å (Figure 8C), aromatic–aromatic interactions within 4.5–7 Å (Figure 8D), protein–protein side chain hydrogen bonds (Figure 8E), and protein–protein main chain hydrogen bonds were characterized (Figure 8F). Interactions remain in a similar range through most part of the dynamics upon loss of contacts in the last part of the simulation time, in agreement with a higher RMSF of this zone in CRF-BP.

In order to explore the relative contribution of these residues in the binding affinity for CRF-BP/CRF_2αR, different point mutations of CRF_2αR, and mutant combinations, were generated and the binding energy values were determined and compared with the WT (Table 3). The deletion of the first 12 amino acids (12aa) has a high contribution to the loss of the binding energy values predicted, and the loss was even higher in the 12aa/Y95Q/S96G/Q97E mutants. These results suggest that the first 12aa, that conform the α-helix, are important contributors and have synergistically effects with Y95, S96, and Q97 for CRF-BP/CRF_2αR binding affinity.