All-atom molecular dynamics simulations for Targeting Human Prostaglandin Reductase 1 with Licochalcone A: Insights from Molecular Dynamics and Covalent Docking Studies

The dataset comprises simulations of the PTGR1 protein under four different conditions: in its apo (unbound) form, bound to the cofactor NADH, and in complex with both covalently and non-covalently bound licochalcone A. Each simulation was conducted using the ff19SB force field and the OPC water model, with water molecules excluded from the trajectories.

For the apo form, NADH-bound, and non-covalently bound licochalcone A conditions, each trajectory consists of 5000 snapshots, representing a total of 500 nanoseconds of simulation. However, the trajectory for the no covalently bound licochalcone A condition includes only 1000 frames, corresponding to 150 nanoseconds. This discrepancy in frame count and simulation length across conditions is important to consider when comparing dynamics and structural behavior within the dataset.

PTGR1-NADPH.tar.xz - PTGR1 dimer in complex with NADPHPTGR1-apo.tar.xz - PTGR1 apo dimerPTGR1-monomer_NADPH.tar.xz - PTGR1 monomer in complex with NADPHPTGR1-LicA_covalent.tar.xz - PTGR1 dimer with covalently bound licochalcone APTGR1-LicA_NO_covalent.tar.xz - PTGR1 dimer with covalently bound licochalcone A

All folders contain:

*.parm7 - dry topology in amber format

*.nc - dry trajectories in netcdf format

Plain molecular dynamics simulations. Two structures of human PTGR1 have been deposited in the PDB: one bound to NADPH and the raloxifene inhibitor (PDB ID 2Y05, 2.2 Å resolution), and another in apo form (PDB ID 1ZSV, 2.3 Å resolution). In the first structure, it is reported as a monomer, whereas in the second as a dimer. However, the protomers exhibit highly similar conformations in both structures (backbone RMSD ~0.5 Å). Experimental evidence, akin to PTGR1 orthologs and many other MDR enzymes, indicates that the functional form of human PTGR1 is a homodimer (Mesa et al., 2015). Thus, to construct the dimeric form of the coenzyme complex, we duplicated the NADPH-bound protomer from 2Y05 and aligned the two subunits with the dimer from 1ZSV, deleting the raloxifene molecule. In the resulting structure, no steric clashes between the protomers were observed. This structure was used as the starting point for the simulations. Additionally, the apo dimer was generated by removing the NADPH from both subunits, and the monomeric form in complex with NADPH was derived from the initial structure.  

The pmemd.cuda module of AMBER 22 was used to perform the MD simulations, employing the force field FF19SB and the OPC water model (Case et al., 2022; Izadi, Anandakrishnan, & Onufriev, 2014; Salomon-Ferrer, Götz, Poole Duncan and Le Grand, & Walker, 2013; Tian et al., 2020). NADPH parameters were taken from (Cummins, Ramnarayan, Singh, & Gready, 1991). The system was protonated at pH 7.4 with PDBfixer (Eastman et al., 2017a) and placed in a truncated octahedral box, initially spanning 12 Å further from the solute in each direction using the AMBER tLeap module. The overall charge of the system was neutralized by the addition of four sodium ions. ParmEd (Eastman et al., 2017b) was used to implement the hydrogen mass repartitioning scheme (Hopkins, Le Grand, Walker, & Roitberg, 2015). Local clashes and solvent orientation were corrected using the steepest descent algorithm for 5,000 cycles. During the initial NVT equilibration, the velocities gradually increased through five steps of 200 ps each. The temperature progression started at 150 K and was raised to 200 K, 250 K, 300 K, and finally, 310 K. Position restraints were applied to heavy atoms of the protein, with the restraining forces progressively decreasing at each step. The spring constants were set at 4, 5, 3, and 1 kcal/mol Å2, respectively, to allow for the gradual relaxation of the protein. The system was further equilibrated for 1 ns in the NPT ensemble with no restraints. For treating long-range electrostatic interactions, periodic boundary conditions and Ewald sums were used with a 9 Å cutoff for direct interactions (Darden, York, & Pedersen, 1993; Simmonett & Brooks, 2021). The same cutoff was used for Lennard-Jones interactions. The Langevin thermostat (Sindhikara, Kim, Voter, & Roitberg, 2009) with a collision frequency of 4 ps-1 and the Monte Carlo barostat(Åqvist, Wennerström, Nervall, Bjelic, & Brandsdal, 2004) with a pressure relaxation time of 2 ps were used to control temperatures and pressures, respectively. The SHAKE algorithm was used to fix any bond involving hydrogen atoms (Ryckaert, Ciccotti, & Berendsen, 1977), and a 4-fs time step integration was used. This protocol was taken from (Cofas-Vargas et al., 2022; Medrano‐Cerano et al., 2024) Unless otherwise stated, no other constraints were used. Five replicas of 500 ns each per system were produced.

The topology and parameter files for a LicA molecule and for this inhibitor covalently bound to the sulfur atom of a cysteine residue were generated with Antechamber suite (J. Wang, Wang, Kollman, & Case, 2006), using the general Amber force field (GAFF2) for organic molecules (He, Man, Yang, Lee, & Wang, 2020). Atomic charges were derived using the AM1-BB method (Jakalian, Jack, & Bayly, 2002). The parameters are documented in Supplementary Tables SI-1 and SI-2. Trajectories for PTGR1 covalently and noncovalently bound to LicA were run using the same conditions as described above. All molecular structure representations were created using UCSF ChimeraX v1.8 (Meng et al., 2023; Pettersen et al., 2021).

Solvent-site identification and guided docking. Determination of solvent sites (SS) for ethanol and water molecules was conducted by employing the MDmix method. After removing both NADPH molecules from the enzyme dimer, the system was protonated at pH 7.4 with PDBfixer (Eastman et al., 2017a) and placed in a truncated octahedral box of water/ethanol 80/20% v/v, extending12 Å beyond the solute in each direction using the AMBER tLeap module. Five 20 ns replicas were run, using the same conditions described above, but applying Cartesian restrictions of 0.01 kcal/mol A2 over all heavy atoms. After the alignment of trajectories, density maps for probe atoms were generated by constructing a static mesh with cubic grids (0.5 Å edge length) over the entire simulation box. The occurrence of probe atoms within each grid were tracked across the trajectories. These density distributions were then converted into binding free energy using the Boltzmann relationship, comparing observed probe atom distributions against the expected bulk solvent distribution at 1.0 M. Solvent sites were then filtered by applying an energy threshold of 1 kcal/mol, as previously described (Alvarez-Garcia & Barril, 2014; Avila-Barrientos et al., 2022).

LicA docking. For covalent docking, LicA, bound through its Cb atom to the sulfur atom of C239, was docked onto the NADPH-binding site of human PTGR1 employing the covalent docking module of AutoDock4 v4.2.6 (Bianco, Forli, Goodsell, & Olson, 2016; Morris et al., 2009). The flexible side-chain methodology was used. In a subsequent non-covalent docking, solvent sites previously identified for ethanol and water were used as pharmacophoric element for rDock (Ruiz-Carmona et al., 2014). This docking involved defining the receptor system and generating a binding cavity using the NADPH as a reference molecule. During the non-covalent docking, a penalty score proportional to the square of the distance from each ligand conformation to a solvent site (SS) was applied when the separation exceeded 2 Å. The docking run included 100 simulations, generating a set of potential binding modes for LicA within the NADPH site.