The protein structures for studies were taken from RCSB PDB (Research
Collaboratory for Structural Bioinformatics Protein Data Bank) (Table 1). In study I, a crystal structure of PPARγ‐RXRα heterodimer was utilized as a template to construct the PPARα‐RXRα heterodimer by replacing the PPARγ with a PPARα structure using Sybyl 7.3 sofware (Tripos, 2006) with the coactivator peptides being removed. The structures for ligands GW409544 (GW, (Xu et al., 2001)) and 9-cis-retinoic acid (9cRA) were obtained from the above-mentioned crystal structures and the parameters and charges were assigned in line with the GROMOS 45a3 united atom force field (Schuler et al., 2001).
Table 1. Protein structures used in studies I–III.
Study Protein PDB code Reference
I PPARγ‐RXRα 1K74 (Xu et al., 2001)
I PPARα 1K7L (Xu et al., 2001)
II ERα 1X7R (Manas, Xu et al., 2004)
II ERβ 1X7B (Manas, Unwalla et al., 2004)
III AR with DHT 1T7T (Hur et al., 2004)
III AR with
bicalutamide 1Z95 (Bohl, Gao et al., 2005)
PDB= Protein data bank, PPAR= Peroxisome proliferator-activated receptors, RXR= Retinoid X receptor, ER= Estrogen receptor, AR= androgen receptor, DHT= dihydrotestosterone
In study II, the crystal structures of ERα and ERβ, both crystallized with an agonistic compound, genistein, were used in docking studies. The ligands were constructed using Accelrys Discovery Studio (Accelrys Software Inc., 2006).
In study III, the protein coordinates were taken from a crystal structure of AR in a complex with dihydrotestosterone (DHT). Accelrys Discovery Studio (Accelrys Software Inc., 2006) was used to construct the novel AR ligand structures. The
ligands were minimized using AMBER 9.0 (Case et al., 2006) program with Amber ff99SB force field (Hornak et al., 2006).
5.1.1 MD simulation methods (I, II, III)
In study I four simulation systems were constructed:
1) PPARα-RXRα heterodimer with no ligands,
2) PPARα-RXRα heterodimer with ligand 9cRA bound to RXRα 3) PPARα-RXRα heterodimer with ligand GW bound to PPARα 4) PPARα-RXRα heterodimer with both ligands
The systems were solvated with SPC water (Berendsen, Postma et al. 1981).
GROMOS 45a3 force field (Schuler et al., 2001) was used in the simulations with GROMACS 3.0 program (Lindahl et al., 2001), constraining all the bond lengths of the protein and ligands with LINCS (Hess et al., 1997) and bond lengths and angles of SPC waters (Berendsen et al., 1981) with SETTLE (Miyamoto and Kollman, 1992).
The time step was 2 fs. After equilibration of 300 ps, new starting velocities were assigned randomly to create 10 parallel MD simulations for each system. The 2 ns simulations were run at 300 K at a constant pressure, controlling temperature and pressure with Berendsen weak coupling algorithm (Berendsen et al., 1984) and storing snapshots every 0.1 ps.
In study II, Amber 9 program (Case et al., 2006) was used with ff99SB forcefield (Hornak et al., 2006). The simulations were performed at 300K. The simulation length was 600 ps after a 300 ps equilibrium run.
In study III the crystal structures of AR with DHT (PDB code 1T7T, (Bohl, Gao et al., 2005)) was used for starting coordinates for the protein. MD simulations (not reported) for ligand-protein complex were performed with Amber 9 and NAMD (Phillips et al., 2005) using ff99SB force field (Hornak et al., 2006).
5.1.2 Fluctuation analysis (I)
To measure movements of the PPARα-RXRα heterodimer in the simulations, root mean square fluctuations (RMSF, equation 4) around the average position of Cα atoms were measured from the last 1 ns of the MD simulations. In equation 4, ri
represents the position of atom i and
〈𝑟
𝑖〉
the average position of atom i over all snapshots.𝜌
𝑖𝑅𝑀𝑆𝐹= √〈(𝑟
𝑖− 〈𝑟
𝑖〉)
2〉
(4)To measure the RMSF, the snapshot structures were fitted over the rigid core of the receptors, each monomer separately. P-values of t-testing were used to visualize the differences between fluctuations in the simulations of PPARα-RXRα-complex with different combinations of ligands.
5.1.3 Clustering analysis (I)
The clustering method used is based on fitting snapshot structures of the
simulations and calculating root mean square deviation (RMSD, equation 5) of the fitted structures (Daura et al., 1999).
𝜌
𝑖𝑅𝑀𝑆𝐷(𝑡) = √
1𝑁
∑
𝑁𝑖=1(𝑟
𝑖(𝑡) − 𝑟
𝑖𝑟𝑒𝑓)
2(5)
In the equation 5, N is the number of atoms, ri(t) is the position of atom ri at time t and riref represents the position of atom i in the reference coordinates.
Structures with RMSD under a specified cut-off value (Table 2) are considered as neighbours. A structure with the largest number of neighbours is designated as the central member of the first cluster, and taken out of the pool along with all its neighbours. The same procedure is repeated until all structures are clustered. The cut-off values were selected visually to show maximum differences between simulation systems. Multiple cut-offs were used to avoid selection bias.
The clustering was performed separately for the whole LBD of PPARα and RXRα, for the helix 8-9 loop and for the co-activator binding site with multiple different cut-off values (Table 2).
Table 2. Clustering analysis.
residues
LBD 202–229, 268–468 227–241,264–458 0.075-0.085
helix 8-9 loop area 383–394 376–387 0.03-0.045
co-activator binding site 282–293, 302–311 274–285, 294–303 0.035-0.05 LBD= Ligand binding domain
5.1.4 Molecular docking methods (II, III)
In study II, the molecular docking calculations were performed with GLIDE program in the Schrödinger molecular modelling package (Friesner et al., 2006;
Friesner et al., 2004), using flexible docking and extra precision mode. The possibility of hydrogen bonding was inferred from distances between potential donor and acceptor pairs in the docked structures.
In study III, the new structures were sketched directly into the AR ligand-binding pocket with the Accelrys Discover Studio (Accelrys Software Inc., 2006). The crystal structures of AR with bicalutamide (PDB code 1Z95, (Bohl, Gao et al., 2005)) and a related compounds (Bohl, Miller et al., 2005) were used to find the correct positioning of the A ring, which is similar in all these ligands. The rest of the novel structure was designed with a drawing tool, guided by the interactions with the receptor and consultations with the synthesis group about possible structures.
5.1.5 Binding free energies (II)
The binding free energies to ERα and ERβ were estimated using the MM-PBSA (Kuhn and Kollman, 2000; Kollman, Peter A. et al., 2000). The structures were extracted from the first 150 ps of the MD trajectories. MM-PBSA calculations included the structural waters involved in the hydrogen bonding network of the hydroxyl group of the ligand, Arg394/Arg346 (ERα/ERβ) and Glu353/Glu305.
In addition, the MM-PBSA results were confirmed with thermodynamic integration (TI) method (Kollman, Peter, 1993). The end structure of a 600 ps simulation with ERα-(R)50 complex was used as a starting structure and atomic point charges of (R)50 were used for both R and S isomers of 50 in the calculation. The same procedure was repeated, using the end structure of a simulation with the ERα-(S)50 complex and ERα-(S)50 point charges. 100 ps equilibration simulations followed by 450 ps production simulations were used for each of the nine different lambda values.