Thermodynamic Integration Calculations Research Articles

Discovering potential WRN inhibitors from natural product database through computational methods

Microsatellite instability (MSI) is a relatively common feature associated with multiple cancers, and Werner syndrome (WRN) ATP-dependent helicase has been recognized as a novel target for treating MSI cancers, such as colorectal cancer. A small-molecule inhibitor targeting WRN would be a promising strategy for treating colorectal cancer with high MSI expression. In this study, we employed a computer-assisted drug discovery strategy to screen over 30,000 natural product molecules. By using a combination of docking, ligand efficiency, Molecular Mechanics/Generalized Born Surface Area (MM/GBSA), and thermodynamic integration (TI) calculations, we identified MOL008980, MOL010740, MOL011832, T4743, TN1166, and TNP-002173 as potential WRN inhibitors. Subsequent molecular dynamics simulation revealed that these screened natural products possessed better binding dynamic characteristics than ATP substrate and were capable of inhibiting the dynamic process of WRN, making them potential strong ATP competitive inhibitors. In conclusion, our computational approach revealed potential WRN inhibitors from a natural product database, providing a theoretical basis for future research.

Development and test of highly accurate endpoint free energy methods. 1: Evaluation of ABCG2 charge model on solvation free energy prediction and optimization of atom radii suitable for more accurate solvationfreeenergy prediction by the PBSA method.

Accurate estimation of solvation free energy (SFE) lays the foundation for accurate prediction of binding free energy. The Poisson-Boltzmann (PB) or generalized Born (GB) combined with surface area (SA) continuum solvation method (PBSA and GBSA) have been widely used in SFE calculations because they can achieve good balance between accuracy and efficiency. However, the accuracy of these methods can be affected by several factors such as the charge models, polar and nonpolar SFE calculation methods and the atom radii used in the calculation. In this work, the performance of the ABCG2 (AM1-BCC-GAFF2) charge model as well as other two charge models, that is, RESP (Restrained Electrostatic Potential) and AM1-BCC (Austin Model 1-bond charge corrections), on the SFE prediction of 544 small molecules in water by PBSA/GBSA was evaluated. In order to improve the performance of the PBSA prediction based on the ABCG2 charge, we further explored the influence of atom radii on the prediction accuracy and yielded a set of atom radius parameters for more accurate SFE prediction using PBSA based on the ABCG2/GAFF2 by reproducing the thermodynamic integration (TI) calculation results. The PB radius parameters of carbon, oxygen, sulfur, phosphorus, chloride, bromide and iodine, were adjusted. New atom types, on, oi, hn1, hn2, hn3, were introduced to further improve the fitting performance. Then, we tuned the parameters in the nonpolar SFE model using the experimental SFE data and the PB calculation results. By adopting the new radius parameters and new nonpolar SFE model, the root mean square error (RMSE) of the SFE calculation for the 544 molecules decreased from 2.38 to 1.05 kcal/mol. Finally, the new radius parameters were applied in the prediction of protein-ligand binding free energies using the MM-PBSA method. For the eight systems tested, we could observe higher correlation between the experiment data and calculation results and smaller prediction errors for the absolute binding free energies, demonstrating that our new radius parameters can improve the free energy calculation using the MM-PBSA method.

Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery.

Rigorous physics-based methods to calculate binding free energies of protein-ligand complexes have become a valued component of structure-based drug design. Relative and absolute binding free energy calculations have been deployed prospectively in support of solving diverse drug discovery challenges. Here we review recent applications of binding free energy calculations to fragment growing and linking, scaffold hopping, binding pose validation, virtual screening, covalent enzyme inhibition, and positional analogue scanning. Furthermore, we discuss the merits of using protein models and highlight recent efforts to replace costly binding free energy calculations with predictions from machine learning models trained on a limited number of free energy perturbation or thermodynamic integration calculations thereby allowing for extended chemical space exploration.

Molecular Insights into the Calcium Binding in Troponin C through a Molecular Dynamics Study

Calcium-binding proteins play critical roles in various biological processes such as signal transduction, cell growth, and transcription factor regulation. Ion binding and target binding of Ca2+-binding proteins are highly related. Therefore, understanding the ion binding mechanism will benefit the relevant inhibitor design toward the Ca2+-binding proteins. The EF-hand is the typical ion binding motif in Ca2+-binding proteins. Previous studies indicate that the ion binding affinity of the EF-hand increases with the peptide length, but this mechanism has not been fully understood. Herein, using molecular dynamics simulations, thermodynamic integration calculations, and molecular mechanics Poisson-Boltzmann surface area analysis, we systematically investigated four Ca2+-binding peptides containing the EF-hand loop in site III of rabbit skeletal troponin C. These four peptides have 13, 21, 26, and 34 residues. Our simulations reproduced the observed trend that the ion binding affinity increases with the peptide length. Our results implied that the E-helix motif preceding the EF-hand loop, likely the Phe99 residue in particular, plays a significant role in this regulation. The E-helix has a significant impact on the backbone and side-chain conformations of the Asp103 residue, rigidifying important hydrogen bonds in the EF-hand and decreasing the solvent exposure of the Ca2+ ion, hence leading to more favorable Ca2+ binding in longer peptides. The present study provides molecular insights into the ion binding in the EF-hand and establishes an important step toward elucidating the responses of Ca2+-binding proteins toward the ion and target availability.

Grid inhomogeneous solvation theory for cross-solvation in rigid solvents.

Grid Inhomogeneous Solvation Theory (GIST) has proven useful to calculate localized thermodynamic properties of water around a solute. Numerous studies have leveraged this information to enhance structure-based binding predictions. We have recently extended GIST toward chloroform as a solvent to allow the prediction of passive membrane permeability. Here, we further generalize the GIST algorithm toward all solvents that can be modeled as rigid molecules. This restriction is inherent to the method and is already present in the inhomogeneous solvation theory. Here, we show that our approach can be applied to various solvent molecules by comparing the results of GIST simulations with thermodynamic integration (TI) calculations and experimental results. Additionally, we analyze and compare a matrix consisting of 100 entries of ten different solvent molecules solvated within each other. We find that the GIST results are highly correlated with TI calculations as well as experiments. For some solvents, we find Pearson correlations of up to 0.99 to the true entropy, while others are affected by the first-order approximation more strongly. The enthalpy-entropy splitting provided by GIST allows us to extend a recently published approach, which estimates higher order entropies by a linear scaling of the first-order entropy, to solvents other than water. Furthermore, we investigate the convergence of GIST in different solvents. We conclude that our extension to GIST reliably calculates localized thermodynamic properties for different solvents and thereby significantly extends the applicability of this widely used method.

Partitioning of noble gases (He, Ne, Ar, Kr, Xe) during Earth’s core segregation: A possible core reservoir for primordial noble gases

Noble gas isotopes in plume-influenced oceanic-island basalts (OIBs) clearly indicate the existence of distinct reservoirs somewhere in the deep Earth that are characterized by their more primitive noble gas isotopic signatures. As one of the largest differentiated layers on Earth, whether the Earth’s liquid outer core can preserve ‘primordial’ noble gases is fundamental for better understanding the Earth’s formation and the mantle evolution. Since Earth’s core is believed to have formed via liquid iron segregated from magma ocean during the planetary accretion stage, how much noble gases could have been incorporated into Earth’s core depends on the partitioning of noble gases between liquid metal and molten silicates. In this study, we perform thermodynamic integration (TI) calculations based on first-principles molecular dynamics (FPMD) simulation to investigate the partition coefficients of noble gases (NG: He, Ne, Ar, Kr, Xe) between liquid Fe and MgSiO3 melt. The simulations have been performed by adopting conditions that the pressures are from 10 ∼ 135 GPa and the temperatures from 2300 ∼ 5000 K to fully account for possible core-mantle equilibrium conditions. By analyzing the temperature and pressure dependence of the results, we suggest that the metal-silicate partitioning of noble gases (DNG) merely depends on the temperature, with extremely low DNG < 0.001 at temperatures < 2500 K (20 GPa), but increases to different degrees at high temperatures. The DNG tends to be more affected by temperature as the noble gas atomic mass increases, with the least temperature effect on the partition coefficient of He (DHe), which increases from ∼ 10−3 at 2300 K (10 GPa) to ∼ 10−1 (5000 K, 135 GPa) near the core-mantle boundary (CMB), and DNe = 10−6 ∼ 10−3, DAr = 10−6 ∼ 10−2, DKr = 10−7 ∼ 10−1 and DXe = 10−5 ∼ 100 from 2300 to 5000 K, respectively. Analysis of electronic interactions indicates that noble gases have positive chemical bonding with Fe, Mg, and Si atoms, and the bonding strengths increase with pressure as well as increasing noble gas atomic mass. Based on our calculated partition coefficients, we suggest that Earth’s core may have preserved significant amounts of noble gases considering the volatile-rich environment from which the core segregated.

Slow Growth Approach for Lithium Ion Deposition on Lithium Metal Anode Surfaces

Lithium metal anodes are indispensable for next-generation batteries including Lithium-Sulfur and Lithium-air batteries which show high promise for electric vehicles and other energy storage applications. Challenges such as extreme reactivity and lithium dendrite formation have kept lithium metal anodes away from practical applications. Lithium dendrite formation on anode surfaces during the charging process originates from the inhomogeneous deposition of metallic lithium during the electrochemical reduction of lithium cations. However, the underlying mechanisms of Li ion deposition from the electrolyte solution onto the anode surface are still poorly understood due to its inherent complexity. In this work, Density Functional Theory (DFT) simulations and thermodynamic integration calculations are conducted to study the electron and ion transfer between lithium metal slab and the electrolyte in absence of an external field. The slow growth approach allow the energy barriers on the diffusion and deposition pathway to be studied. This method illuminate the underlying mechanisms on lithium-ion deposition and provide a better understanding of the early stages of dendrite nucleation. We explore the effect of the solvent chemistry and structure, distance of the solvated complex from the surface, anion-cation separation, and concentration of Li-salts on the deposition of lithium ions from the electrolyte phase onto the surface. It is found that the structure and properties of the solvation shell around the lithium cation has a great influence on the ability of the cation to diffuse as well as on its surrounding electron environment.

Phase Boundary and Salt Partitioning in Coacervate Complexes Formed between Poly(acrylic acid) and Poly(N,N-dimethylaminoethyl methacrylate) from Detailed Atomistic Simulations Combined with Free Energy Perturbation and Thermodynamic Integration Calculations

With the help of detailed atomistic simulations, we have determined the phase boundary of a complex coacervate system resulting from the complexation of two oppositely and fully charged weak polyel...

Polar residues lining the binding cleft of a Serratia marcescens family 18 chitinase position the substrate for attack and stabilize associative interactions

Glycoside hydrolases (GHs) bind carbohydrates to convert polysaccharides into soluble sugars. GH binding sites are lined with aromatic and polar residues that facilitate function. The role of polar residues in GH function, relative to aromats, is often less straightforward. Using molecular dynamics simulation, thermodynamic integration and umbrella sampling molecular dynamics free energy calculations, and biochemical assays, we examined how Asp-316, Tyr-145, Arg-294, and Glu-221 aid in substrate binding, activity, and processive ability of a GH family 18 bacterial chitinase, Serratia marcescens chitinase B (ChiB). Despite what appears to be a ‘roof’ between Asp-316 and Trp-97 over the substrate in crystal structures, our results indicated ‘roof’ formation was not thermodynamically favourable. Asp-316 appeared to be supportive, providing interactions aiding in ligand stability at the −2 and −3 product sites. Mutation of Tyr-145 to alanine inactivated the enzyme; simulations revealed Tyr-145 stabilised the binding cleft strand containing the catalytic motif. Arg-294 appeared to play several roles, both stabilising Asp-215 and the +1 pyranose and forcing Trp-220 to stack with the +2 pyranose. Arg-294 may even modulate active site pKa. Glu-221 hydrogen bonded with the +3 pyranose and appeared to stabilise the adjacent carbohydrate-π stacking between Phe-190 and the +3 pyranose.

Predicting changes in protein stability caused by mutation using sequence-and structure-based methods in a CAGI5 blind challenge.

Predicting the impact of mutations on proteins remains an important problem. As part of the CAGI5 frataxin challenge, we evaluate the accuracy with which Provean, FoldX, and ELASPIC can predict changes in the Gibbs free energy of a protein using a limited data set of eight mutations. We find that different methods have distinct strengths and limitations, with no method being strictly superior to other methods on all metrics. ELASPIC achieves the highest accuracy while also providing a web interface which simplifies the evaluation and analysis of mutations. FoldX is slightly less accurate than ELASPIC but is easier to run locally, as it does not depend on external tools or datasets. Provean achieves reasonable results while being computational less expensive than the other methods and not requiring a structure of the protein. In addition to methods submitted to the CAGI5 community experiment, and with the aim to inform about other methods with high accuracy, we also evaluate predictions made by Rosetta's ddg_monomer protocol, Rosetta's cartesian_ddg protocol, and thermodynamic integration calculations using Amber package. ELASPIC still achieves the highest accuracy, while Rosetta's catesian_ddg protocol appears to perform best in capturing the overall trend in the data.

Energetic Fingerprinting of Ligand Binding to Paralogous Proteins: The Case of the Apoptotic Pathway.

Networks of biological molecules are key to cellular function, governing processes ranging from signal cascade propagation to metabolic pathway regulation. Genetic duplication processes give rise to sets of regulatory proteins that have evolved from a common ancestor, leading to interactomes whose dysregulation is often associated with disease. A better understanding of the determinants of specificity at interfaces shared by functionally related proteins is crucial to the rational design of novel pharmacotherapeutic agents. To this end, a comprehensive data set of drug and drug-like binders was assembled for the Bcl-xL and Bcl-2 antiapoptotic proteins-archetypal examples of regulatory systems governed by evolutionarily conserved protein-protein interactions. These were first used to derive a two-dimensional quantitative structure-activity relationship (2D QSAR) model, predicting ligand specificity for these homologous proteins. The strengths and weaknesses of high-throughput 2D QSAR were then compared and contrasted to those of theoretically rigorous thermodynamic integration calculations performed on 14 complexes of Bcl-xL-specific, Bcl-2-specific, and potent dual binders bound to the Bcl-xL and Bcl-2 proteins. We demonstrate that free energy calculations provide an added layer of essential information, which traditional QSAR cannot capture. Moreover, we show that protein energetic responses to different ligands, expressed as per-residue energy values, can be used to fingerprint the protein-ligand interaction, extending the framework of four-dimensional molecular dynamics/quantitative structure-activity relationships (4D-MD/QSAR) toward the facilitation of future drug design strategies.

Effects of Single Nucleotide Polymorphisms on the Binding of Afatinib to EGFR: A Potential Patient Stratification Factor Revealed by Modeling Studies

The anticancer drug afatinib has been found to be more effective at inhibiting the oncogenic EGFR mutant exon 19 deletion (19del) over the oncogenic EGFR mutant L858R. The underlying mechanism has been hypothesized to result from differences in structural constraints introduced by the mutations and stabilizing interactions afforded by a buried water molecule in 19del (Kannan, S.; et al. Sci. Rep. 2017, 7, 1540). The COSMIC cancer database is mined for EGFR sequences to discover that several mutations in the form of Single nucleotide polymorphisms (SNPs) line this hydration cavity. In this work, the effects of these SNPs on the affinity of afatinib for EGFRWT and oncogenic mutants EGFRL858R and EGFR19del were studied using free energy perturbation and thermodynamic Integration calculations. The simulations reveal that several SNPs have significant effects on the affinity of afatinib for the mutant EGFRs carrying the SNPs and may thus have clinical implications relating to emergence of resistance to afatinib, thus potentially impacting the choice of EGFR inhibitors in the clinic.

Classical and quantum calculations of the temperature dependence of the free energy of argon

The free energy is central to statistical mechanics and thermodynamics, and its accurate calculation via. computational modelling is important for a large number of applications, especially when its experimental value is hard to obtain. Several established and general methods for calculating the Helmholtz free energy across different length scales, including continuum, atomistic and quantum mechanical, are compared and analyzed. A computational approach is then proposed to calculate the temperature dependences of internal energy and absolute Helmholtz free energy for solid and liquid phases with the coupling of thermodynamic integration (TI) and harmonic approximation calculations from both classical molecular dynamics (MD) and density functional theory (DFT). We use the Lennard-Jones system as an example (i.e. argon) for the demonstration of the approach. It is observed that the free energy transits smoothly from being describable by the harmonic approximation to including anharmonic effects at a transition temperature around 0.56 Tm; below this temperature, the quantum behavior of atoms is important. At higher temperatures (T > 0.56 Tm), the TI and harmonic approximation results for the Helmholtz free energy functions become increasingly divergent with the increase of temperature. This work demonstrates that a multiscale approach employing TI, MD, and DFT can provide accurate calculations of the temperature dependence of absolute Helmholtz free energy for both solid and liquid phases.

Role of Urea-Aromatic Stacking Interactions in Stabilizing the Aromatic Residues of the Protein in Urea-Induced Denatured State.

A delicate balance of different types of intramolecular interactions makes the folded states of proteins marginally more stable than the unfolded states. Experiments use thermal, chemical, or mechanical stress to perturb the folding equilibrium for examining protein stability and the protein folding process. Elucidation of the mechanism by which chemical denaturants unfold proteins is crucial; this study explores the nature of urea-aromatic interactions relevant in urea-assisted protein denaturation. Free energy profiles corresponding to the unfolding of Trp-cage miniprotein in the presence and absence of urea at three different temperatures demonstrate the distortion of the hydrophobic core to be a crucial step. Exposure of the Trp6 residue to the solvent is found to be favored in the presence of urea. Previous experiments showed that urea has a high affinity for aromatic groups of proteins. We show here that this is due to the remarkable ability of urea to form stacking and NH-π interactions with aromatic groups of proteins. Urea-nucleobase stacking interactions have been shown to be crucial in urea-assisted RNA unfolding. Examination of these interactions using microsecond-long unrestrained simulations shows that urea-aromatic stacking interactions are stabilizing and long lasting. Further MD simulations, thermodynamic integration, and quantum mechanical calculations on aromatic model systems reveal that such interactions are possible for all the aromatic amino acid side-chains. Finally, we validate the ubiquitous nature of urea-aromatic stacking interactions by analyzing experimental structures of urea transporters and proteins crystallized in the presence of urea or urea derivatives.

Surfactant pKa Calculations using Molecular Dynamics Simulations

The protonation state for ionizable surfactants is crucial for self-assembly, phase behavior, and morphology. For example, oleic acid forms vesicles at intermediate pH conditions (∼7-9), micelles at high pH, and oil phases at low pH. The pKa for oleic acid shifts from ∼4.8 as monomers in water to ∼7.5 in model phospholipid vesicles. It remains challenging to accurately model surfactant self-assembly and protonation behavior with molecular dynamics simulations. We use constant pH molecular dynamics simulations and thermodynamic integration calculations for protonating a single fatty acid in different chemical environments.

Accuracy comparison of several common implicit solvent models and their implementations in the context of protein-ligand binding

In this study several commonly used implicit solvent models are compared with respect to their accuracy of estimating solvation energies of small molecules and proteins, as well as desolvation penalty in protein-ligand binding. The test set consists of 19 small proteins, 104 small molecules, and 15 protein-ligand complexes. We compared predicted hydration energies of small molecules with their experimental values; the results of the solvation and desolvation energy calculations for small molecules, proteins and protein-ligand complexes in water were also compared with Thermodynamic Integration calculations based on TIP3P water model and Amber12 force field. The following implicit solvent (water) models considered here are: PCM (Polarized Continuum Model implemented in DISOLV and MCBHSOLV programs), GB (Generalized Born method implemented in DISOLV program, S-GB, and GBNSR6 stand-alone version), COSMO (COnductor-like Screening Model implemented in the DISOLV program and the MOPAC package) and the Poisson-Boltzmann model (implemented in the APBS program). Different parameterizations of the molecules were examined: we compared MMFF94 force field, Amber12 force field and the quantum-chemical semi-empirical PM7 method implemented in the MOPAC package. For small molecules, all of the implicit solvent models tested here yield high correlation coefficients (0.87–0.93) between the calculated solvation energies and the experimental values of hydration energies. For small molecules high correlation (0.82–0.97) with the explicit solvent energies is seen as well. On the other hand, estimated protein solvation energies and protein-ligand binding desolvation energies show substantial discrepancy (up to 10kcal/mol) with the explicit solvent reference. The correlation of polar protein solvation energies and protein-ligand desolvation energies with the corresponding explicit solvent results is 0.65–0.99 and 0.76–0.96 respectively, though this difference in correlations is caused more by different parameterization and less by methods and indicates the need for further improvement of implicit solvent models parameterization. Within the same parameterization, various implicit methods give practically the same correlation with results obtained in explicit solvent model for ligands and proteins: e.g. correlation values of polar ligand solvation energies and the corresponding energies in the frame of explicit solvent were 0.953–0.966 for the APBS program, the GBNSR6 program and all models used in the DISOLV program. The DISOLV program proved to be on a par with the other used programs in the case of proteins and ligands solvation energy calculation. However, the solution of the Poisson-Boltzmann equation (APBS program) and Generalized Born method (implemented in the GBNSR6 program) proved to be the most accurate in calculating the desolvation energies of complexes.

The importance of protonation and tautomerization in relative binding affinity prediction: a comparison of AMBER TI and Schrödinger FEP.

In drug discovery, protonation states and tautomerization are easily overlooked. Through a Merck-Rutgers collaboration, this paper re-examined the initial settings and preparations for the Thermodynamic Integration (TI) calculation in AMBER Free-Energy Workflows, demonstrating the value of careful consideration of ligand protonation and tautomer state. Finally, promising results comparing AMBER TI and Schrödinger FEP+ are shown that should encourage others to explore the value of TI in routine Structure-based Drug Design.

Rapid approximate calculation of water binding free energies in the whole hydration domain of (bio)macromolecules.

The evaluation of water binding free energies around solute molecules is important for the thermodynamic characterization of hydration or association processes. Here, a rapid approximate method to estimate water binding free energies around (bio)macromolecules from a single molecular dynamics simulation is presented. The basic idea is that endpoint free-energy calculation methods are applied and the endpoint quantities are monitored on a three-dimensional grid around the solute. Thus, a gridded map of water binding free energies around the solute is obtained, that is, from a single short simulation, a map of favorable and unfavorable water binding sites can be constructed. Among the employed free-energy calculation methods, approaches involving endpoint information pertaining to actual thermodynamic integration calculations or endpoint information as exploited in the linear interaction energy method were examined. The accuracy of the approximate approaches was evaluated on the hydration of a cage-like molecule representing either a nonpolar, polar, or charged water binding site and on α- and β-cyclodextrin molecules. Among the tested approaches, the linear interaction energy method is considered the most viable approach. Applying the linear interaction energy method on the grid around the solute, a semi-quantitative thermodynamic characterization of hydration around the whole solute is obtained. Disadvantages are the approximate nature of the method and a limited flexibility of the solute. © 2016 Wiley Periodicals, Inc.

Hydration free energies calculated using the AMBER ff03 charge model for natural and unnatural amino acids and multiple water models

In this work, we assess calculated hydration free energies of natural and unnatural amino acids compared to experimental hydration free energies of corresponding side-chain analogs using the ff03 charge model. Fine-grid, explicit water thermodynamic integration calculations using two widely used explicit water models (TIP3P and TIP4P-Ew) were performed on 19 natural amino acids and compared with experimental hydration free energies of corresponding side-chain analogs to establish expected accuracy levels for this charge model. Next, parameters previously derived for 17 unnatural amino acids and several new parameters optimized in this work were assessed using the same methodology. We found that the ff03 charge model is correlated with experimental hydration free energies but underestimates the solubilities of several polar natural and unnatural amino acids. Comparisons are presented with other forcefields and water models recently presented in the literature for both natural and modified amino acids.

Evaluation of Host-Guest Binding Thermodynamics of Model Cavities with Grid Cell Theory.

A previously developed cell theory model of liquid water was used to evaluate the excess thermodynamic properties of confined clusters of water molecules. The results are in good agreement with reference thermodynamic integration calculations, suggesting that the model is adequate to probe the thermodynamic properties of water at interfaces or in cavities. Next, the grid cell theory (GCT) method was applied to elucidate the thermodynamic signature of nonpolar association for a range of idealized host-guest systems. Polarity and geometry of the host cavities were systematically varied, and enthalpic and entropic solvent components were spatially resolved for detailed graphical analyses. Perturbations in the thermodynamic properties of water molecules upon guest binding are restricted to the immediate vicinity of the guest in solvent-exposed cavities, whereas longer-ranged perturbations are observed in buried cavities. Depending on the polarity and geometry of the host, water displacement by a nonpolar guest makes a small or large enthalpic or entropic contribution to the free energy of binding. Thus, no assumptions about the thermodynamic signature of the hydrophobic effect can be made in general. Overall the results warrant further applications of GCT to more complex systems such as protein-ligand complexes.