ORIGINAL RESEARCH article
Aromatic Rings Commonly Used in Medicinal Chemistry: Force Fields Comparison and Interactions With Water Toward the Design of New Chemical Entities
- 1Grupo de Bioinformática Estrutural, Centro de Biotecnologia, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil
- 2Swiss National Supercomputing Centre, Lugano, Switzerland
- 3Instituto de Informática, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil
- 4Instituto Aggeu Magalhães, Fundação Oswaldo Cruz, Recife, Brazil
The identification of lead compounds usually includes a step of chemical diversity generation. Its rationale may be supported by both qualitative (SAR) and quantitative (QSAR) approaches, offering models of the putative ligand-receptor interactions. In both scenarios, our understanding of which interactions functional groups can perform is mostly based on their chemical nature (such as electronegativity, volume, melting point, lipophilicity etc.) instead of their dynamics in aqueous, biological solutions (solvent accessibility, lifetime of hydrogen bonds, solvent structure etc.). As a consequence, it is challenging to predict from 2D structures which functional groups will be able to perform interactions with the target receptor, at which intensity and relative abundance in the biological environment, all of which will contribute to ligand potency and intrinsic activity. With this in mind, the aim of this work is to assess properties of aromatic rings, commonly used for drug design, in aqueous solution through molecular dynamics simulations in order to characterize their chemical features and infer their impact in complexation dynamics. For this, common aromatic and heteroaromatic rings were selected and received new atomic charge set based on the direction and module of the dipole moment from MP2/6-31G* calculations, while other topological terms were taken from GROMOS53A6 force field. Afterwards, liquid physicochemical properties were simulated for a calibration set composed by nearly 40 molecules and compared to their respective experimental data, in order to validate each topology. Based on the reliance of the employed strategy, we expanded the dataset to more than 100 aromatic rings. Properties in aqueous solution such as solvent accessible surface area, H-bonds availability, H-bonds residence time, and water structure around heteroatoms were calculated for each ring, creating a database of potential interactions, shedding light on features of drugs in biological solutions, on the structural basis for bioisosterism and on the enthalpic/entropic costs for ligand-receptor complexation dynamics.
The development of a drug is a multi step process, usually starting with the identification of hit compounds. The challenging task of optimizing these compounds into leads and finally into drugs is commonly facilitated by computer aided drug design (CADD) techniques (Anderson, 2003; Sliwoski et al., 2013; Bajorath, 2015). With the growing information on protein structure on the last years, structure based drug design (SBDD) has become a significant tool for hit discovery (Anderson, 2003; Lounnas et al., 2013; Lionta et al., 2014). When structural information of the receptor is absent, molecular fingerprints of approved drugs are also used to search for new ligands in a process also known as ligand based drug design (LBDD) (Lee et al., 2011). Nevertheless, there are still considerable challenges associated to the predictiveness of ligand potency and affinity via computational methods (Paul et al., 2010; Csermely et al., 2012).
In general, optimization of lead compounds is based in qualitative or quantitative structure-activity relationships (SAR or QSAR, respectively) (Shahlaei, 2013). These relationships are usually based in molecular descriptors to predict ligand pharmacodynamics and pharmacokinetics, such as logP to access lipophilicity, logS to access solubility or pKa to access the ionic state of a compound, along with other topological, geometrical and physicochemical descriptors (Danishuddin and Khan, 2016). While some correlations have reasonable power of predictiveness, many descriptors have no biological meaning and can mislead the optimization process. As highlighted by Hopkins et al. (2014), high-throughput screening methods have been linked to the rise of hits with inflated physicochemical properties during the optimization process (Keserü and Makara, 2009). Also, recent reviews have shown an increase of molar mass in the recent medicinal chemistry efforts (Leeson and Springthorpe, 2007) and many authors correlate this strategy with the likelihood of poor results of such compounds (Gleeson, 2008; Waring, 2009, 2010; Gleeson et al., 2011).
Many chemical moieties are regularly used in medicinal chemistry to produce chemical diversity (Bemis and Murcko, 1996; Welsch et al., 2010; Taylor et al., 2014), a practice well-known as fragment based drug design (FBDD), and its use for pharmacophore modeling and to prevent high toxicity is not recent (Gao et al., 2010). Particularly, aromatic rings are extensively used in drugs due to their well known synthetic and modification paths (Aldeghi et al., 2014). For example, at least, one aromatic ring can be found in 99% of a database containing more than 3,500 evaluated by the medicinal chemistry department of Pfizer, AstraZeneca (AZ) and GlaxoSmithKlin (GSK) (Roughley and Jordan, 2011). Still, little is known about their chemical features in biological solution, such as H-bonds availability, lifetime of H-bonds, solvent accessibility, and conformational ensemble. In this sense, molecular dynamics (MD) simulations can provide useful information with atomistic resolution and access the aforementioned features of chemical groups in water, providing fundamental data to drive medicinal chemistry approaches.
Still, dynamical properties of chemical moieties in biological solution are usually neglected in drug design and very difficult to access (Ferenczy and Keseru, 2010; Reynolds and Holloway, 2011; Hopkins et al., 2014). Even though MD simulations have been used in medicinal chemistry to generate different receptor conformers and to validate binding poses predicted by docking (Zhao and Caflisch, 2015; Ganesan et al., 2017), simulations of free ligand in solution is rarely used to access the conformational ensemble and energies associated with solvation due to the challenge on solving conformational flexibility and internal energies (Butler et al., 2009; Blundell et al., 2016). When solvated, the enthalpic and entropic costs of disrupting a H-bond or dismantling the entire solvation shell of a ligand can be the determinant step to provide the proper energy of binding (Biela et al., 2012; Blundell et al., 2013; Mondal et al., 2014). Yet, free-energy of binding is often predicted via geometrical or alchemical transformations (Zwanzig, 1954; Aqvist et al., 1994; Woo and Roux, 2005; Gumbart et al., 2013), alongside with recent developments in funnel metadynamics (Limongelli et al., 2013). More recently, thermodynamical features of ligands have been experimentally investigated in order to enhance binding and efficiency (Freire, 2009; Ferenczy and Keseru, 2010; Reynolds and Holloway, 2011). Ligand features such as H-bonds lifetime, effects of vicinity in H-bonds availability and strength, accessible surface area and water structure around binding sites can provide substantial information for designing new molecular entities (Blundell et al., 2016).
Different force fields have been used for drug design purposes, such as MMF94 (Halgren, 1996), OPLS-AA (Jorgensen et al., 1996), and GAFF (Wang et al., 2004). While these force fields parameterized their electrostatic terms using ab initio calculations, the GROMOS force fields (derived from the Groningen Molecular Simulation package) used free-energy of solvation as target (Daura et al., 1998; Oostenbrink et al., 2004) to empirically assign atomic partial charges. Thus, in this work, we have chosen the GROMOS force field to simulate the dynamical behavior of 103 aromatic rings (including a calibration subset of 42 molecules) mostly commonly used in drug design and their interactions with solvent in order to access thermodynamical properties in solution. These interactions, in turn, offer a reference for future rational drug design studies, as describe in details how several functional groups interact with their surroundings.
2.1. Selection of Rings
A series of 103 aromatic rings commonly used in drug design were selected for this study (Broughton and Watson, 2004; Jordan and Roughley, 2009; Welsch et al., 2010; Taylor et al., 2014, 2017). Among them, a calibration set of 42 molecules (Table 1), for which physical-chemical properties are known, were selected from the benchmark developed by Caleman et al. (2012). Briefly, both works of Taylor et al. (2014, 2017) employed a detailed search of substructure frequencies from FDA Orange Book and cross referenced with ChEMBL, DrugBank, Nature, Drug Reviews, the FDA Web site, and the Annual Reports in Medicinal Chemistry; the work of Broughton and Watson (2004) employed search of substructure frequencies in MDL Drug Data Report database by using a “Phase II” keyword; and the work of Welsch et al. (2010) have pinpointed privileged scaffolds from natural-products works throughout literature.
2.2. Topology Construction
Structures for these aromatic rings were built using Avogadro (Hanwell et al., 2012). Molecular mechanical (MM) topological parameters as bonds, angles, and Lennard-Jones parameters were taken from GROMOS53A6 (Oostenbrink et al., 2004). Due to the well–known good performance of MP2 methods for small aromatic rings (Li et al., 2015; Matczak and Wojtulewski, 2015), atomic partial charges were based on quantum mechanical (QM) calculations using MP2 theory (Møller and Plesset, 1934), 6-31G* (Petersson et al., 1988) basis set and implicit solvent Polarizable Continuum Model (PCM) (Mennucci and Tomasi, 1997) followed by a RESP fitting (Bayly et al., 1993). The so obtained partial charges were adjusted in the MM to reproduce the QM dipole moment of the ring. The angle θ formed between the QM and MM model dipole moment vectors was monitored through an in house script to make sure the angle had the lowest value possible, guaranteeing the conservation of the QM dipole moment direction. For our calibration set, the module of the MM partial charges were adjusted to better reproduce the physicochemical properties of the organic liquids. Following the philosophy of charge group assignment, groups were limited, at maximum, to the atoms at the ortho position on each ring. In more complex substitution patterns, a superimposition of two charge groups was required to correctly describe the chemical group. In such cases, the Coulombic terms of the overlapping atoms were adjusted to correctly describe the direction of the total dipole moment of the ring. For molecules containing linear constraints (benzonitrile), virtual sites were added in order to preserve the total moment of inertia and mass, thus preserving the linearity of these groups (Feenstra et al., 1999).
2.3. New Torsional Potentials
The quantum mechanical torsional profile of every dihedral angle was calculated using Gaussian (Frisch et al., 2016) (RRID:SCR_014897). Molecular structures were built using Avogadro (Hanwell et al., 2012) and their geometry were optimized using Hartree-Fock method (Fock, 1930; Hartree and Hartree, 1935) and basis set 3-21G* (Dobbs and Hehre, 1986). Afterwards, the Scan routine was used to calculate the total energy of the molecule conformation for each dihedral orientation, adopting a tight convergence criteria, with geometric optimization, MP2/6-31G* and steps of 30°. In order to calculate the torsional profile for molecular mechanics model, dihedral orientations were kept fixed during minimization using restraint forces for the same angles evaluated by quantum calculations. Both profiles were submitted to the Rotational Profiler server (Rusu et al., 2014) to obtain appropriate sets of classical mechanics parameters that provided a better fitting to the QM-obtained torsional profile.
2.4. General Simulation Settings
All simulations were carried out using the GROMACS 5.0.7 package (Abraham et al., 2015) (RRID:SCR_014565). In order to create parameters compatible with the GROMOS family, we have followed previous literature (Daura et al., 1998; Schuler et al., 2001; Oostenbrink et al., 2004) settings: twin-range scheme was used with short- and long-range cutoff distances of 0.8 and 1.4 nm, respectively. Also, the reaction-field method was applied to correct the effects of electrostatic interactions beyond the long-range cutoff distance (Barker and Watts, 1973; Tironi et al., 1995), using the dielectric constant as εRF for organic liquid simulations and εRF = 62 for simulations in water (Heinz et al., 2001; Oostenbrink et al., 2004). The LINCS algorithm (Hess et al., 1997; Hess, 2008) was used to constrain all covalent bonds, using a cubic interpolation, a Fourier grid of 0.12 nm and timestep of 2 fs. Configurations were saved at every 2 ps for analysis.
2.4.1. Organic Liquids Simulations
In order to build the organic liquid systems, cubic boxes of 2 × 2 × 2 nm were created, each with a single organic molecule. A total of 125 of these boxes were stacked, forming an unique box with conventional periodic boundary conditions treatment of 10 × 10 × 10 nm which was simulated under high pressure (100 bar) to induce liquid phase. The systems were then simulated and equilibrated at 1 bar. Afterwards, the boxes were staggered to obtain systems with 1000 molecules in liquid phase and simulated at 1 bar until the total energy drift converged to values below 0.5 J/(mol × ns × Degrees of Freedom). Such criterion is necessary to make sure that the fluctuating properties could be accurately calculated (Caleman et al., 2012). All simulations were carried out with Berendsen pressure and temperature coupling algorithm due to their efficiency in molecular relaxations (Berendsen et al., 1984), using τT = 0.2 ps and τP = 0.5 ps. When available, experimental values of isothermal compressibility and dielectric constant were used as an additional parameter for liquid simulations. Otherwise, the compressibility of the most chemically similar molecule was used. The experimental dielectric constants from each liquid were also used as parameters in the simulations (Oostenbrink et al., 2004).
In order to calculate the densities of liquids (ρ), simulations at constant pressure were carried out for 10 ns and ρ were calculated using block averages of 5 blocks. Enthalpy of vaporization (ΔHvap) were calculated by block averaging the same 10 ns of liquid simulation to obtain Epot(l) and another 100 ns of gas phase simulation using a stochastic dynamics integrator (SD) (Van Gunsteren and Berendsen, 1988) with a single molecule in vacuum, to obtain Epot(g) as the equation:
Aiming to calculate the dielectric constant (ε), the simulation of the liquid boxes from which ρ were obtained were extended up to 60 ns. Convergence calculations of ε were done using running averages and ε were evaluated only after convergence. In order to calculate thermal expansion coefficients (αP) and classic isobaric heat capacities (CPcla), three constant pressure simulations were carried out for 5 ns each, with temperatures T, T+10K, and T-10K, for each liquid. The calculations of αP and CPcla were done using the finite difference method (Kunz and van Gunsteren, 2009):
In order to calculate isothermal compressibilities (κT), three constant volume simulations were carried out for 5 ns each, with pressures 1, 0.9, and 1.1 bar. The calculations of κT was also done using the finite difference method:
2.4.2. Solvation Free Energy Simulations
Simulations in water were carried out to evaluate the solvation free energies (ΔGhyd) of 30 molecules at 1 bar and 298 K. Each aromatic ring (solute) was centered into a cubic box with appropriate dimensions to reproduce the density of SPC water models (0.997 g/cm3). In free-energy calculations using thermodynamic integration (TI) method, a coupling parameter λ is used to perturb solute-solvent interactions.
in which H is the Halmiltonian, λ = 0 refers to the state in which the solute fully interacts with the solvent and λ = 1 refers to the state in which the solute-solvent interactions do not exist. In our setup, Coulombic interactions were decoupled first, and the Lennard-Jones interactions after, using a soft-core potential to avoid issues related to strong Lennard-Jones interactions (Beutler et al., 1994). A soft-core power was set to 1 and αLJ set to 0.5, following recommendations of Shirts and Pande (2005). Both interactions were decoupled using λ values: 0, 0.02, 0.04, 0.07, 0.1, 0.15, 0.2, …, 0.8, 0.85, 0.9, 0.93, 0.96, 0.98, 1, totalizing 50 λ simulations.
Our simulation protocol consisted of an initial steepest-descent minimization, followed by a L-BFGS minimization until a maximum force of 10 kJ/(mol−1 nm−1) was reached. After, initial velocities were assigned and the systems were equilibrated for 100 ps using a NVT ensemble at each λ. The systems were subjected to another 100 ps of equilibration on a NPT ensemble, using the Parrinello-Rahman pressure coupling algorithm (Parrinello and Rahman, 1981), a τt = 5 ps time constant for coupling and a compressibility of 4.5 × 105 bar−1. Finally, production simulations were done using the Langevin integrator (Van Gunsteren and Berendsen, 1988) to sample the 〈∂H/∂λ〉λ until convergence. Therefore, simulations time varied between 1 and 5 ns. In addition, the last frame of the production phase of each λ was used as input for the next subsequent λ.
2.4.3. Simulation of Rings in Water
After an extensive comparison of simulated and experimental physicochemical properties of our calibration set and consequent validation, the same strategy of topological construction was applied to other 61 rings commonly used in drug design (Table 2) for which experimental properties are not available, totalizing 103 aromatic rings in this study. Hence, in order to evaluate chemical features and interactions of aromatic rings with their surroundings, a total set of 103 aromatic was simulated in water, including all 42 molecules present in the calibration set (Table 1). Each solute was placed in a cubic box with a distance of 1.0 nm to its edges. The boxes were then filled with SPC water model and minimized long enough eliminate any possible clashes until convergence at a maximum force of 0.1 kJ/mol × nm. After, the system was equilibrated in a NVT ensemble at 298.15 K using the Nosé-Hoover algorithm (Nosé, 1984) for temperature coupling. Production runs of 250 ns were carried out with temperature and pressure coupling handled by V-rescale (Bussi et al., 2007) and Parrinelo-Rahman (Parrinello and Rahman, 1981) algorithms, using τT = 0.1 ps and τP = 2.0 ps. The GROMACS tools hbond, rdf, and sorient were used to calculate H-bonds related properties and solvation structure around the heteroatom using a block-averaging approach over 5 box of 50 ns.
3.1. New Torsional Profiles
In order to accurately describe the torsional angles of the selected aromatic rings, a total of 15 new dihedral potentials were derived by fitting the MM profiles to the corresponding QM-calculated ones (Table S1). Fittings were conducted using the Rotational Profiler server (Rusu et al., 2014). For all cases, the use of new parameters yield almost identical values of minimum and barrier amplitudes to those calculated by QM (Figure 1). Dihedral distribution throughout simulations was also evaluated.
Figure 1. Evaluation of torsional parameters and dihedral distribution. QM and adjusted MM torsional profiles are shown in black and green, respectively. In red, the dihedral distribution during simulations.
3.2. Physical-Chemical Properties
In order to validate our strategy of topology building, boxes of organic liquids were simulated to obtain physical-chemical properties for each compound. Reference experimental values (Table S2) were used to calculate the absolute error of each property and to guide adjustments on the coulombic terms in order to mitigate deviations. We have calculated the θ angle between QM and MM dipole moments and the final version of our calibration set (Table 1) yielded an average θ angle of 2.5° ± 6.1°, suggesting that our MM models conserve the direction of the QM dipole moment, preserving the electrostatic potential of each molecule.
Following the GROMOS philosophy (Oostenbrink et al., 2004; Horta et al., 2016), density (ρ), enthalpy of vaporization (ΔHvap), and free energy of solvation (ΔGhyd) were used as targets for the parametrization, while isothermal coefficient (αP), isothermal compressibility (κT), dielectric constant (ε), and classic isobaric heat capacity (CPcla) were calculated as benchmarks for GROMOS performance and compared with the results obtained in Caleman et al. (2012) and Horta et al. (2016) (Table 3). Linear regression between experimental and simulated values were calculated in order to access the prediction power of the employed strategy (Figure 2). The equations further reported were calculated excluding outliers (values higher than 2 standard deviations).
Table 3. Average deviation between experimental and simulated physicochemical properties of aromatic rings evaluated in our calibration set. Simulated GAFF and OPLS-AA values were obtained from Caleman et al. (2012) and 2016H66 values from Horta et al. (2016). Density (ρ) in g/cm3, enthalpy of vaporization (ΔHvap) in kJ/mol, thermal expansion coefficient (αP) in 10−3/K, isothermal compressibility (κT) in 1/GPa, dielectric constant (ε), classic isobaric heat capacity (Cpcla) in J/mol × K, and free-energy of solvation (ΔGhyd) in kJ/mol.
Figure 2. Correlation between experimental and calculated physical-chemical properties of organic liquids for 42 aromatic compounds on the calibration set. Standard deviations are shown as bars, linear regressions are shown as green and empty dots represent outliers.
Regarding the targeted properties, our calibration set yielded the equations y = 0.9118x + 0.1001 for density, y = 1.0699x − 1.6491 for enthalpy of vaporization and y = 0.8676x + 0.8929 for free energy of solvation, with correlation coefficients of R = 0.92, R = 0.96, and R = 0.89, respectively. In terms of average deviation (AVED), our calibration set overestimates ρ in 0.008 g/cm3, ΔHvap in 1.51 kJ/mol and underestimates ΔGhyd in 3.35 kJ/mol. Without the outliers, the AVED for ΔGhyd improves to 2.83 kJ/mol.
Non-targeted properties were calculated to evaluate how they behaved in our simulations. Linear regressions yielded equations y = 0.93825x + 0.1406 for αP (R = 0.82), y = +0.90079x − 0.0140 for κT (R = 0.70), y = 0.2581x + 1.8961 for ε (R = 0.65), and y = 0.8989x + 100.5 for Cpcla (R = 0.77). In terms of AVED, αP is overestimated in 0.14 10−3/K and κT is overestimated in 0.0465 1/GPa. As expected (Caleman et al., 2012; Horta et al., 2016), ε is poorly described due to the lack of polarization effects, resulting in a underestimation of −4.52 in the dielectric constant. On other hand, Cpcla was overestimated by 88.2 J/mol × K, a behavior aligned with recent works in literature (Caleman et al., 2012; Horta et al., 2016). Individual AVED and absolute errors can be found in Tables S4, S5 in Supplementary Material, along with experimental properties in Table S3.
3.3. Interactions in Water
In order to quantitatively evaluate the behavior of heteroaromatic rings in water and their interactions with the aqueous surrounding, some properties were calculated throughout 250 ns of simulation. From these calculations, we were capable to assess the average H-bond (AverHB) of each heteroatom along with its residence time (τHB), lifetime (lifetimeHB), the free-energy of breakage of a H-bond (ΔGHB), and the percentage of simulation time that a given heteroatom was involved in, at least, one H-bond (Percent). We were also capable to obtain the optimal binding distance between an heteratom and water (OBDHB), along with the coordination number (CNHB) at the OBDHB and the average orientation of water molecules surrounding the heteroatom. These data are compiled in Tables 4, 5.
Table 4. Properties of heteroaromatic rings in water. Average H-bonds (AverHB), H-bond residence time (τHB) is ps, H-bond lifetime (lifetimeHB) in 1/ps, free-energy of H-bond breakage (ΔGHB) in kJ/mol, percentage of simulation with at least one formed H-bond (Percent.), coordination number of water (CN), optimal binding distance with water (OBDHB) in nm, and overall water orientation around the heteroatom (Orientation).
Table 5. Properties of heteroaromatic rings in water. Average H-bonds (AverHB), H-bond residence time (τHB) is ps, H-bond lifetime (lifetimeHB) in 1/ps, free-energy of H-bond breakage (ΔGHB) in kJ/mol, percentage of simulation with at least one formed H-bond (Percent.), coordination number of water (CN), optimal binding distance with water (OBDHB) in nm, and overall water orientation around the heteroatom (Orientation).
4.1. Topology Building Strategy
The accurate description of organic compounds' chemical diversity, mainly in the context of drugs and medicinal chemistry, is a challenging task in molecular mechanics since it must be described as broadly as possible by the force field fragments. However, the most common sets of MM parameters employed in biomolecules simulations are usually centered on the monomeric constituents of biopolymers and lipids, while parameters for synthetic compounds, as well as other common non-polymeric biological molecules (e.g., natural products), must be included from specific calculations or external sets of parameters.
In this sense, a proper description of torsional terms will impact directly the dynamical behavior of these small molecules, even considering that, when evaluating ligand-receptor complexes, the influence of these terms might be mitigated due to the ligand movement restriction inside the binding pocket. Still, accommodation of flexible docking derived poses, fine tunning of induced fit, and characterization of ligands conformational induction vs. selection (with potential inferences of the entropic costs of binding) require dihedrals potentials specifically adjusted to organic compounds. Hence, new parameters were generated in this work exclusively for 15 dihedrals in aromatic rings in our calibration set (Figure 1). In general, our results revealed that our MM parameters yielded a good description of the QM torsional profile, with the exceptions of  tiophenol,  phenoxybenzene,  phenylmethanol, and  trifluoromethylbenzene. For these molecules, the distribution profile was almost evenly spread, most likely due to the low energy barrier (below 2.5 kJ/mol), indicating that transient states are commonly achieved during our simulations in SPC water model. Simulations of these particular molecules in vacuum revealed little influence of water solvation in the dihedral profile (data not shown).
In another sense, the choice of an atomic charge set for ligands can drastically impact thermodynamical binding properties such as complexation free-energy and desolvation. Therefore, we employed in this work a dipole moment based strategy to describe the Coulombic contribution using physicochemical properties of organic liquids as target. The prediction power of our strategy was compared to recent comparisons of aromatic compounds in liquid phase (Caleman et al., 2012; Horta et al., 2016) and summarized in Table 3. In general, our calibration set yielded similar or lower average deviations than benchmarks made with OPLS-AA, GAFF, and 2016H66 sets for all physicochemical properties evaluated in this work. The main difference was in terms of Cpcla, for which GAFF and OPLS-AA overestimate nearly 40 J/mol × K more than our parameters. Still, all four parameters sets overestimates Cpcla. In addition, the GROMOS53A5 force field was designed to reproduce physicochemical properties, and later on adjusted to reproduce free energy of solvation and hydration (GROMOS53A6) (Oostenbrink et al., 2004). The average deviation on density, enthalpy of vaporization and free-energy of solvation of GROMOS53A5 were 0.0389 g/cm3, −0.4 and 3.8 kJ/mol, respectively. These values are very similar to our results, as shown in Table 3, reiterating the quality of our parameters.
It is important to mention that the employed benchmark set was built using the same Lennard-Jones parameters used in the benzene ring of phenylalanine in GROMOS53A6. While GROMOS53A6 produces a ΔGhyd = 0.0 kJ/mol for benzene (phenylalanine side-chain), our benzene parameters yield a ΔGhyd = −3.4 kJ/mol, a much closer value to the experimental data (ΔGhyd = −3.6 kJ/mol). Nevertheless, the AVED value reveals a underestimation for free energy of hydration in our parameter set. A possible reason is that chemical functions such as nitro, fluorine, chlorine, and aldehydic carbonyls are not commonly found in biomolecules and, therefore, the LJ parameters used in GROMOS53A6 may not be properly extrapolated to synthetic compounds. Moreover, we have tested ether oxygens LJ parameters reported in Horta et al. (2011) in our pure liquid simulations of furan and methoxybenzene, leading to approximately the same behavior in their respective physical-chemical properties (data not shown).
4.2. Properties in Solution: Influence of Nearby Substitutions in H-Bonds
In order to access quantitative informations regarding how aromatic rings interact with their surroundings, we performed molecular dynamics simulations for 103 aromatic rings most commonly used in drug design, including our 42 molecules calibration set. These information are condensed in the Tables 4, 5. Simulations were carried for 250 ns to properly sample multiple events of H-bond breakages and solvation shell rearrangements.
Our results reveal non-obvious information about the H-bond availability and strength, as in the case of pyridine/pyrimidine/pyrazine/pyridazine/triazine series (Figure 3). While exchanging a pyridine by a pyrimidine ring might lead to apparent gain of a H-bond acceptor, nitrogens of pyrimidine present a ΔGHB of nearly 1 kJ/mol lower than pyridine. Moreover, the Percent of time with at least one formed H-bond between water and pyridine nitrogen is higher than the ones in pyrimidine. When comparing pyridine with pyrazine (an addition of another N in para), H-bonds are very similar, so as the second and third solvation layers. Also, acceptance capacity in pyrimidine ring is very similar to triazine, where all three nitrogens are located in meta. Intriguingly, values for pyridine are very similar to the ones calculated for pyridazine, with a slight increase in OBDHB and a more compact second layer of solvation, as shown in Figure 3A. These results suggest that another nitrogen acceptor in meta decreases nitrogen acceptance capacity, while another nitrogen acceptor in ortho has low effect in H-bond capacity, but a considerable effect in the solvation layers structures. In this sense, these features can impact the binding inside receptors. Pyridazine, for example, has a larger OBDHB than pyridine, suggesting that these molecules can occupy the binding pocket in a different manner, impacting the entropic cost of binding.
Figure 3. (A) Methyl substituitions: 2-Me (green), 3-Me (yellow), (B) Nearby N substitution: Northo (green), Nmeta (yellow), 4-Me (purple) and 2,4,6-Me (pink). Npara (purple). Solvation properties of aromatic rings in pyridine family. Radial distribution functions (RDFs) and H-bonding strength of N1 (blue) are affected by substitutions in ortho, meta, and para.
Other cases have been equally surprising, like the quinoline/isoquinoline. The main difference between them is the location of the acceptor nitrogen (closer to C8 in the quinoline fused ring). Counterintuitively, the AverHB of isoquinoline is slightly lower than for quinoline, such as the τHB, and the ΔGHB is almost 1.25 kJ/mol lower. The same properties for pyridine ring are somewhat between these values of quinoline and isoquinoline. In addition, ΔGHB for quinazoline and quinoxaline rings are almost 3 kJ/mol lower than quinoline and isoquinoline. In this sense, quinazoline and quinoxaline would be better candidates in fragment-based drug design due to the lower energetic cost of desolvation, while maintaining the H-bond capacity inside the receptor. Another case in terms of aromatic nitrogen hydrogen bond acceptor is the 2,4,6-trimethylpyridine (Figure 3B). The presence of methyl groups in both ortho positions drastically reduces the availability of H-bonds, as shown in Figure 3, and diminish the residence time of the accepted H-bond. But the presence of only one methyl group in ortho appears to have a modest effect, slightly favoring the presence of H-bond in nitrogen of 2-methylpyridine. Moreover, the second and third solvation layers of 2- and 2,4,6-trimethylpyridine are dismantled, while the same behavior is not observed for 3- and 4-methylpyridine.
Other non-obvious events can be observed regarding H-bond donation in hydroxyls groups. In case of phenol, the necessary energy to break a donated H-bond (~10 kJ/mol) is almost the double to break an accepted one (~5.70 kJ/mol), in alignment with the QM data reported by Parthasarath et al. (2005) in HF, MP2, and DFT level. And while phenol and phenylmethanol might appear interchangeable during the lead optimization process, the ΔGHB of accepted and donated H-bonds in the hydroxyl group is almost 1 kJ/mol higher for phenylmethanol. While targeting thermodynamics of binding during drug design, these energy costs of desolvation can play a crucial role. As expected, benzenethiol was revealed to be a poor acceptor of hydrogen bonds in our simulations, but a reasonable H-bond donator. In terms of vicinity effects, methylation in ortho seems to have little effect on hydroxyl groups, since the properties evaluated for the series phenol/2-methylphenol/3-methylphenol/4-methylphenol have very similar behavior.
It is well– know that halogens are widely used for drug design, and the role of halogen bonds (X-bonds) and H-bonds role have been investigated thoroughly (Rendine et al., 2011; Ford and Ho, 2016; Lin and Mackerell, 2017). In general, the H-bonding strength decreases with the halogen radius (F > Cl > Br > I), while the halogen bond strength increases (Rendine et al., 2011). In this work, we investigated how fluorine and chlorine behave as H-bond acceptors in water. In the case of fluorobenzene, the ΔGHB = 1.54 ± 0.24 is in accordance with a weak H-bond (Domagała et al., 2017). The other fluorinated rings in the series (1,2-, 1,3-, 1,2,3,4-, and 1,2,3,5-tetrafluorobenzene [8-11]) have similar values, varying from 1.5 to 2.2 kJ/mol. Regarding the chlorinated rings series (chlorobenzene, 1,2-, 1,3-, 1,2,3,4-, and 1,2,3,5-tetrachlorobenzene [94–98]), ΔGHB ranged from 1.80 to 3.24 kJ/mol, contradicting the expected behavior. X-bonding are often poorly described in MM, since it treats atoms as a sphere with isoelectric surface and thus not describing the necessary positive potential required for such interaction. In fact, we have visually evaluated that waters surrounding fluorine and chlorine have their hydrogens oriented toward the halogens, confirming our measure of H-bonds and not X-bonds.
Regarding oxygen atoms within the aromatic ring, AverHB are generally lower than expected. It is well known that oxygens in heterocycles act as H-bond acceptor (Kaur and Khanna, 2011), but our model does not reproduce this tendency. It is important to notice that GROMOS53A6 does not have specific parameters for oxygens within aromatic rings, and LJ parameters from ethers were employed. Not surprisingly, the calculated properties for the oxygen atom in furan and benzofuran are very similar to methoxybenzene and phenoxybenzene. This result suggests that the description of the properties in aqueous solutions of aromatic rings containing oxygen might be improved by specific LJ parameters. Moreover, we have tested ether LJ parameters reported in Horta et al. (2011) for our simulations of furan and methoxybenzene in water, yielding lower AverHB and ΔGHB (data not shown). The new force field parameters developed in this work can be obtained upon request.
4.3. Impacts in Drug Design
Recently, several authors have questioned the LE approach as optimization tool and its actual power to lead to high affinity compounds (Abad-Zapatero, 2007; Morgan et al., 2011; Cavalluzzi et al., 2017). Another recent review (DeGoey et al., 2017) has pointed out the emergence of approved drugs that violate Lipinski's rules of 5 and correlated them to properties such as number of aromatic rings and rotatable bonds. Freire (2009) have proposed an experimental thermodynamic approach to guide the drug design process and these results led to believe that tweaking ligand enthalpy and entropy of binding is not only experimentally possible, but also possible to predict. Therefore, the GROMOS series of force fields present an extra advantage here due to their calibration to reproduce free-energy of solvation and other thermodynamical properties.
In this sense, we have parameterized and validated a calibration set of 42 aromatic rings commonly used in drug design using thermodynamical properties in condensed phase. After, we performed a study with a larger dataset of 103 heteroaromatic rings in order to understand how these molecules interact with water and to prospect and map potential interactions with target-receptors. The water molecules probe the occurrence of hydrogen bonds, and the absence of these interactions, as well as the distance from the first solvation sphere, may probe sites for hydrophobic interactions. With these information at hand, medicinal chemists and pharmacologists may employ quantitative estimations on how each functional group may or may not interact with its target protein, as well as identify the potential influence of close chemical modifications. These properties (and a handful of others) are compiled in Tables 4, 5, and can be used as reference during lead optimization process.
The strategy employed here could be used to amplify the spectrum of drug fragments with accurate description of chemical events simulated by molecular dynamics. In addition, it can improve the description of drug-receptor complexation dynamics of other molecules of interest, molecular recognition of drugs and signal transduction mediated by conformational changes of ligands. In fact, by assessing the strength and availability of interactions between aromatic rings and water solvent, the results presented here not only offer detailed quantitative information about potential interactions that each individual aromatic ring can make with its surrounding, but also shed light upon the energetics of biological events, such as dismantling solvation shells — an important step in the ligand binding process.
In this work, we have successfully produced topologies for a calibration set of 42 aromatic rings using as target physicochemical properties of respective organic liquids. Our strategy revealed a very competitive prediction power when compared alongside with other force fields, while presenting a simple approach to describe aromatic rings through molecular dynamics simulations that can be easily extrapolated to other rings. In addition to that, H-bond availability and solvent accessibility are difficult and non-obvious informations to predict from bidimensional data, but still essential for medicinal chemistry purposes. Here, we have simulated in aqueous solvent more than 100 aromatic rings commonly used in drug design in order to assess dynamical chemical properties, such as average H-bonds, their lifetime, residence time and free energy of breakage. Thus, we have described a low cost approach based on molecular dynamics simulations to access valuable information that could be useful both to predict the enthalpic cost of desolvation and for interpretation of pharmacological data by a medicinal chemist or pharmacologist. Our results provide a large database of quantitative information for a total of 103 aromatic rings most commonly used in drug design that can guide medicinal chemists in future drug design efforts.
MP carried out quantum calculations, molecular dynamics simulations, data analyses, and drafted the manuscript. VR contributed in the simulations protocols and manuscript draft. BG wrote in house scripts for dipole-based charge assignment and data analyses. MD contributed to manuscript draft. RL contributed to simulations protocols and manuscript draft. HV contributed to data analyses and manuscript draft.
The authors thank the funding agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Fundação de Amparo à Pesquisa do Rio Grande do Sul (FAPERGS). This work was partially supported by grants from FAPERGS/PRONUPEQ (16/2551-0000520-6).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The reviewer GT and handling Editor declared their shared affiliation.
Research developed with support of the Centro Nacional de Supercomputação (CESUP), from Universidade Federal do Rio Grande do Sul (UFRGS). We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2018.00395/full#supplementary-material
Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). Gromacs: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1-2, 19–25. doi: 10.1016/j.softx.2015.06.001
Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97, 10269–10280.
Beutler, T. C., Mark, A. E., van Schaik, R. C., Gerber, P. R., and van Gunsteren, W. F. (1994). Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 222, 529–539.
Biela, A., Khayat, M., Tan, H., Kong, J., Heine, A., Hangauer, D., et al. (2012). Impact of ligand and protein desolvation on ligand binding to the S1 pocket of thrombin. J. Mol. Biol. 418, 350–366. doi: 10.1016/j.jmb.2012.01.054
Blundell, C. D., Nowak, T., and Watson, M. J. (2016). Measurement, interpretation and use of free ligand solution conformations in drug discovery. Prog. Med. Chem. 55, 45–147. doi: 10.1016/bs.pmch.2015.10.003
Blundell, C. D., Packer, M. J., and Almond, A. (2013). Quantification of free ligand conformational preferences by NMR and their relationship to the bioactive conformation. Bioorg. Med. Chem. 21, 4976–4987. doi: 10.1016/j.bmc.2013.06.056
Caleman, C., van Maaren, P. J., Hong, M., Hub, J. S., Costa, L. T., and van der Spoel, D. (2012). Force field benchmark of organic liquids: Density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theor. Comput. 8, 61–74. doi: 10.1021/ct200731v
Cavalluzzi, M. M., Mangiatordi, G. F., Nicolotti, O., and Lentini, G. (2017). Ligand efficiency metrics in drug discovery: the pros and cons from a practical perspective. Expert Opin. Drug Dis. 12, 1087–1104. doi: 10.1080/17460441.2017
Csermely, P., Korcsmáros, T., Kiss, H. J., London, G., and Nussinov, R. (2012). Structure and dynamics of molecular networks: A novel paradigm of drug discovery. A comprehensive review. Pharmacol. Ther. 138, 333–408. doi: 10.1016/j.pharmthera.2013.01.016
DeGoey, D. A., Chen, H. J., Cox, P. B., and Wendt, M. D. (2017). Beyond the rule of 5: lessons learned from AbbVie's drugs and compound collection. J. Med. Chem. 61, 2636–2651. doi: 10.1021/acs.jmedchem.7b00717
Dobbs, K. D., and Hehre, W. J. (1986). Molecular orbital theory of the properties of inorganic and organometallic compounds 4. Extended basis sets for third-and fourth-row, main-group elements. J. Comput. Chem. 7, 359–378.
Gleeson, M. P., Hersey, A., Montanari, D., and Overington, J. (2011). Probing the links between in vitro potency, ADMET and physicochemical parameters. Nat. Rev. Drug Discov. 10, 197–208. doi: 10.1038/nrd3367
Hanwell, M. D., Curtis, D. E., Lonie, D. C., Vandermeersch, T., Zurek, E., and Hutchison, G. R. (2012). Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. Cheminformatics 4:17. doi: 10.1186/1758-2946-4-17
Heinz, T. N., van Gunsteren, W. F., and Hünenberger, P. H. (2001). Comparison of four methods to compute the dielectric permittivity of liquids from molecular dynamics simulations. J. Chem. Phys. 115, 1125–1136. doi: 10.1063/1.1379764
Hopkins, A. L., Keserü, G. M., Leeson, P. D., Rees, D. C., and Reynolds, C. H. (2014). The role of ligand efficiency metrics in drug discovery. Nat. Rev. Drug Discov. 13, 105–121. doi: 10.1038/nrd4163
Horta, B. A., Merz, P. T., Fuchs, P. F., Dolenc, J., Riniker, S., and Hünenberger, P. H. (2016). A GROMOS-compatible force field for small organic molecules in the condensed phase: the 2016H66 parameter set. J. Chem. Theor. Comput. 12, 3825–3850. doi: 10.1021/acs.jctc.6b00187
Horta, B. A., Fuchs, P. F., van Gunsteren, W. F., and Hünenberger, P. H. (2011). New interaction parameters for oxygen compounds in the GROMOS force field: Improved pure-liquid and solvation properties for alcohols, ethers, aldehydes, ketones, carboxylic acids, and esters. J. Chem. Theor. Comput. 7, 1016–1031. doi: 10.1021/ct1006407
Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996). Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225–11236.
Kunz, A. P., and van Gunsteren, W. F. (2009). Development of a nonlinear classical polarization model for liquid water and aqueous solutions: COS/D. J. Phys. Chem. A 113, 11570–11579. doi: 10.1021/jp903164s
Li, S., Smith, D. G., and Patkowski, K. (2015). An accurate benchmark description of the interactions between carbon dioxide and polyheterocyclic aromatic compounds containing nitrogen. Phys. Chem. Chem. Phys. 17, 16560–16574. doi: 10.1039/c5cp02365c
Lin, F. Y., and Mackerell, A. D. (2017). Do halogen-hydrogen bond donor interactions dominate the favorable contribution of halogens to ligand-protein binding? J. Phys. Chem. B 121, 6813–6821. doi: 10.1021/acs.jpcb.7b04198
Lionta, E., Spyrou, G., Vassilatis, D. K., and Cournia, Z. (2014). Structure-based virtual screening for drug discovery: principles, applications and recent advances. Curr. Top. Med. Chem. 14, 1923–1938. doi: 10.2174/1568026614666140929124445
Lounnas, V., Ritschel, T., Kelder, J., McGuire, R., Bywater, R. P., and Foloppe, N. (2013). Current progress in structure-based rational drug design marks a new mindset in drug discovery. Comput. Struct. Biotechnol. J. 5:e201302011. doi: 10.5936/csbj.201302011
Matczak, P., and Wojtulewski, S. (2015). Performance of Møller-Plesset second-order perturbation theory and density functional theory in predicting the interaction between stannylenes and aromatic molecules. J. Mol. Model. 21, 41. doi: 10.1007/s00894-015-2589-1
Oostenbrink, C., Villa, A., Mark, A. E., and van Gunsteren, W. F. (2004). A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 25, 1656–1676. doi: 10.1002/jcc.20090
Paul, S. M., Mytelka, D. S., Dunwiddie, C. T., Persinger, C. C., Munos, B. H., Lindborg, S. R., et al. (2010). How to improve R&D productivity: the pharmaceutical industry's grand challenge. Nat. Rev. Drug Discov. 9, 203. doi: 10.1038/nrd3078
Petersson, G. A., Bennett, A., Tensfeldt, T. G., Al-Laham, M. A., Shirley, W. A., and Mantzaris, J. (1988). A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements. J. Chem. Phys. 89, 2193–2218.
Rendine, S., Pieraccini, S., Forni, A., and Sironi, M. (2011). Halogen bonding in ligand—receptor systems in the framework of classical force fields. Phys. Chem. Chem. Phys. 13:19508. doi: 10.1039/c1cp22436k
Shirts, M. R., and Pande, V. S. (2005). Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration. J. Chem. Phys. 122:144107. doi: 10.1063/1.1873592
Taylor, R. D., MacCoss, M., and Lawson, A. D. (2017). Combining molecular scaffolds from FDA approved drugs: application to drug discovery. J. Med. Chem. 60, 1638–1647. doi: 10.1021/acs.jmedchem.6b01367
Waring, M. J. (2009). Defining optimum lipophilicity and molecular weight ranges for drug candidates-Molecular weight dependent lower log D limits based on permeability. Bioorg. Med. Chem. Lett. 19, 2844–2851. doi: 10.1016/j.bmcl.2009.03.109
Keywords: drug design, GROMOS, aromatic rings, functional groups, interactions
Citation: Polêto MD, Rusu VH, Grisci BI, Dorn M, Lins RD and Verli H (2018) Aromatic Rings Commonly Used in Medicinal Chemistry: Force Fields Comparison and Interactions With Water Toward the Design of New Chemical Entities. Front. Pharmacol. 9:395. doi: 10.3389/fphar.2018.00395
Received: 09 November 2017; Accepted: 05 April 2018;
Published: 24 April 2018.
Edited by:Adriano D. Andricopulo, University of São Paulo, Brazil
Reviewed by:Antonio Monari, Université de Lorraine, France
Gustavo Trossini, Universidade de São Paulo, Brazil
Copyright © 2018 Polêto, Rusu, Grisci, Dorn, Lins and Verli. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hugo Verli, firstname.lastname@example.org