Alchemical Free Energy Calculations Research Articles

Superimposing Ligands with a Ligand Overlay as an Alternate Topology Model for λ-Dynamics-Based Calculations.

Alchemical free energy (AFE) calculations can predict binding affinity changes as a function of structural modifications and have become powerful tools for lead optimization and drug discovery. Central to the setup and performance of AFE calculations is the manner of mapping alchemical transformations, known as the topology model. Single, dual, and hybrid topology models have been used with various AFE methods in the field. In recent works, λ-dynamics (λD) free energy calculations, specifically, have preferred the use of a hybrid multiple topology (HMT) for sampling multiple ligand perturbations. In this work, we evaluate a new topology method called ligand overlay (LO) for use with λD-based calculations, including the recently introduced λ-dynamics with a bias-updated Gibbs sampling (LaDyBUGS) approach. LO is a full multiple topology model that allows entire ligands to be sampled and restrained within a λ-dynamics framework. Relative binding free energies were computed with HMT or LO topology models with LaDyBUGS for 45 ligands across five protein benchmark systems. An overall Pearson R correlation of 0.98 and mean unsigned error of 0.32 kcal/mol were observed, suggesting that LO is a viable alternative topology model for λD-based calculations. We discuss the merits of using an HMT or LO model for future ligand studies with λD or LaDyBUGS calculations.

Comparison of Methodologies for Absolute Binding Free Energy Calculations of Ligands to Intrinsically Disordered Proteins.

Modulating the function of Intrinsically Disordered Proteins (IDPs) with small molecules is of considerable importance given the crucial roles of IDPs in the pathophysiology of numerous diseases. Reported binding affinities for ligands to diverse IDPs vary broadly, and little is known about the detailed molecular mechanisms that underpin ligand efficacy. Molecular simulations of IDP ligand binding mechanisms can help us understand the mode of action of small molecule inhibitors of IDP function, but it is still unclear how binding energies can be modeled rigorously for such a flexible class of proteins. Here, we compare alchemical absolute binding free energy calculations (ABFE) and Markov-State Modeling (MSM) protocols to model the binding of the small molecule 10058-F4 to a disordered peptide extracted from a segment of the oncoprotein c-Myc. The ABFE results produce binding energy estimates that are sensitive to the choice of reference structure. In contrast, the MSM results produce more reproducible binding energy estimates consistent with weak mM binding affinities and transient intermolecular contacts reported in the literature.

Comparative Assessment of Water Models in Protein-Glycan Interaction: Insights from Alchemical Free Energy Calculations and Molecular Dynamics Simulations.

Accurate computational simulations of protein-glycan dynamics are crucial for a comprehensive understanding of critical biological mechanisms, including host-pathogen interactions, immune system defenses, and intercellular communication. The accuracy of these simulations, including molecular dynamics (MD) simulation and alchemical free energy calculations, critically relies on the appropriate parameters, including the water model, because of the extensive hydrogen bonding with glycan hydroxyl groups. However, a systematic evaluation of water models' accuracy in simulating protein-glycan interaction at the molecular level is still lacking. In this study, we used full atomistic MD simulations and alchemical absolute binding free energy (ABFE) calculations to investigate the performance of five distinct water models in six protein-glycan complex systems. We evaluated water models' impact on structural dynamics and binding affinity through over 5.8 μs of simulation time per system. Our results reveal that most protein-glycan complexes are stable in the overall structural dynamics regardless of the water model used, while some show obvious fluctuations with specific water models. More importantly, we discover that the stability of the binding motif's conformation is dependent on the water model chosen when its residues form weak hydrogen bonds with the glycan. The water model also influences the conformational stability of the glycan in its bound state according to density functional theory (DFT) calculations. Using alchemical ABFE calculations, we find that the OPC water model exhibits exceptional consistency with experimental binding affinity data, whereas commonly used models such as TIP3P are less accurate. The findings demonstrate how different water models affect protein-glycan interactions and the accuracy of binding affinity calculations, which is crucial in developing therapeutic strategies targeting these interactions.

Automated Adaptive Absolute Binding Free Energy Calculations.

Alchemical absolute binding free energy (ABFE) calculations have substantial potential in drug discovery, but are often prohibitively computationally expensive. To unlock their potential, efficient automated ABFE workflows are required to reduce both computational cost and human intervention. We present a fully automated ABFE workflow based on the automated selection of λ windows, the ensemble-based detection of equilibration, and the adaptive allocation of sampling time based on inter-replicate statistics. We find that the automated selection of intermediate states with consistent overlap is rapid, robust, and simple to implement. Robust detection of equilibration is achieved with a paired t-test between the free energy estimates at initial and final portions of a an ensemble of runs. We determine reasonable default parameters for all algorithms and show that the full workflow produces equivalent results to a nonadaptive scheme over a variety of test systems, while often accelerating equilibration. Our complete workflow is implemented in the open-source package A3FE (https://github.com/michellab/a3fe).

Exploring the molecular interactions of dorzolamide hydrochloride in aqueous solutions: a comprehensive study using thermophysical techniques and molecular simulations across various temperatures

This research begins with the initial assessment of key physical properties—density ( ρ), speed of sound ( u), and refractive index ( nD ) in binary aqueous solutions of dorzolamide hydrochloride. These properties were measured in solutions with varying molalities, ranging from 0.0 to 0.0276, under standard atmospheric conditions and at temperatures between 300.15 and 320.15 K. From these measurements, we were able to calculate a range of important parameters, including the apparent molar volume ( Vϕ), its limiting values ([Formula: see text]), limiting apparent molar expansivity ([Formula: see text]), Helper’s constant ([Formula: see text]), apparent molar isentropic compression ( κϕ), and its limiting values ([Formula: see text]), specific refraction ( RD), molar refraction ( RM), and hydration number ( nh). These calculated values provided insights into the strong interactions between the solute and solvent in the liquid system. They also shed light on the presence of various molecular interactions, such as inter- and intramolecular hydrogen bonds, dipole–dipole, and dipole-induced-dipole interactions, as well as the ability of the solute to influence the structure of the surrounding water molecules. Additionally, we derived the densities of dorzolamide hydrochloride solutions from molecular dynamics (MD) simulations, which showed a strong correlation with our experimental findings, underscoring the reliability of MD simulations in such studies. In the end, our combined study on dorzolamide hydrochloride in water demonstrates that it shapes structures because it has lower volumetric properties and stronger hydrogen bonds in MD simulations. More alchemical free energy calculations show that there are big changes in solvation energies. This is important for understanding how dorzolamide hydrochloride works in water.

High throughput mutational scanning of a protein via alchemistry on a high-performance computing resource.

Antibiotic resistance presents a significant challenge to public health, as bacteria can develop resistance to antibiotics through random mutations during their life cycles, making the drugs ineffective. Understanding how these mutations contribute to drug resistance at the molecular level is crucial for designing new treatment approaches. Recent advancements in molecular biology tools have made it possible to conduct comprehensive analyses of protein mutations. Computational methods for assessing molecular fitness, such as binding energies, are not as precise as experimental techniques like deep mutational scanning. Although full atomistic alchemical free energy calculations offer the necessary precision, they are seldom used to assess high throughput data as they require significantly more computational resources. We generated a computational library using deep mutational scanning for dihydrofolate reductase (DHFR), a protein commonly studied in antibiotic resistance research. Due to resource limitations, we analyzed 33 out of 159 positions, identifying 16 single amino acid replacements. Calculations were conducted for DHFR in its drug-free state and in the presence of two different inhibitors. We demonstrate the feasibility of such calculations, made possible due to the enhancements in computational resources and their optimized use.

Hit-to-Lead Optimization of Heterocyclic Carbonyloxycarboximidamides as Selective Antagonists at Human Adenosine A3 Receptor.

Antagonism of the human adenosine A3 receptor (hA3R) has potential therapeutic application. Alchemical relative binding free energy calculations of K18 and K32 suggested that the combination of a 3-(2,6-dichlorophenyl)-isoxazolyl group with 2-pyridinyl at the ends of a carbonyloxycarboximidamide group should improve hA3R affinity. Of the 25 new analogues synthesized, 37 and 74 showed improved hA3R affinity compared to K18 (and K32). This was further improved through the addition of a bromine group to the 2-pyridinyl at the 5-position, generating compound 39. Alchemical relative binding free energy calculations, mutagenesis studies and MD simulations supported the compounds' binding pattern while suggesting that the bromine of 39 inserts deep into the hA3R orthosteric pocket, so highlighting the importance of rigidification of the carbonyloxycarboximidamide moiety. MD simulations highlighted the importance of rigidification of the carbonyloxycarboximidamide, while suggesting that the bromine of 39 inserts deep into the hA3R orthosteric pocket, which was supported through mutagenesis studies 39 also selectively antagonized endogenously expressed hA3R in nonsmall cell lung carcinoma cells, while pharmacokinetic studies indicated low toxicity enabling in vivo evaluation. We therefore suggest that 39 has potential for further development as a high-affinity hA3R antagonist.

Replica Exchange of Expanded Ensembles: A Generalized Ensemble Approach with Enhanced Flexibility and Parallelizability.

Generalized ensemble methods such as Hamiltonian replica exchange (HREX) and expanded ensemble (EE) have been shown effective in free energy calculations for various contexts, given their ability to circumvent free energy barriers via nonphysical pathways defined by states with different modified Hamiltonians. However, both HREX and EE methods come with drawbacks, such as limited flexibility in parameter specification or the lack of parallelizability for more complicated applications. To address this challenge, we present the method of replica exchange of expanded ensembles (REXEE), which integrates the principles of HREX and EE methods by periodically exchanging coordinates of EE replicas sampling different yet overlapping sets of alchemical states. With the solvation free energy calculation of anthracene and binding free energy calculation of the CB7-10 binding complex, we show that the REXEE method achieves the same level of accuracy in free energy calculations as the HREX and EE methods, while offering enhanced flexibility and parallelizability. Additionally, we examined REXEE simulations with various setups to understand how different exchange frequencies and replica configurations influence the sampling efficiency in the fixed-weight phase and the weight convergence in the weight-updating phase. The REXEE approach can be further extended to support asynchronous parallelization schemes, allowing looser communications between larger numbers of loosely coupled processors such as cloud computing and therefore promising much more scalable and adaptive executions of alchemical free energy calculations. All algorithms for the REXEE method are available in the Python package ensemble_md, which offers an interface for REXEE simulation management without modifying the source code in GROMACS.

Ligand Binding Affinity Prediction for Membrane Proteins with Alchemical Free Energy Calculation Methods.

Alchemical relative binding free energy (ΔΔG) calculations have shown high accuracy in predicting ligand binding affinity and have been used as important tools in computer-aided drug discovery and design. However, there has been limited research on the application of ΔΔG methods to membrane proteins despite the fact that these proteins represent a significant proportion of drug targets, play crucial roles in biological processes, and are implicated in numerous diseases. In this study, to predict the binding affinity of ligands to G protein-coupled receptors (GPCRs), we employed two ΔΔG calculation methods: thermodynamic integration (TI) with AMBER and the alchemical transfer method (AToM) with OpenMM. We calculated ΔΔG values for 53 transformations involving four class A GPCRs and evaluated the performance of AMBER-TI and AToM-OpenMM. In addition, we conducted tests using different numbers of windows and varying simulation times to achieve reliable ΔΔG results and to optimize resource utilization. Overall, both AMBER-TI and AToM-OpenMM show good agreement with the experimental data. Our results validate the applicability of AMBER-TI and AToM-OpenMM for optimization of lead compounds targeting membrane proteins.

Correct Nucleotide Selection Is Confined at the Binding Site of Polymerase Enzymes.

DNA polymerases (Pols) add incoming nucleotides (deoxyribonucleoside triphosphate (dNTPs)) to growing DNA strands, a crucial step for DNA synthesis. The insertion of correct (vs incorrect) nucleotides relates to Pols' fidelity, which defines Pols' ability to faithfully replicate DNA strands in a template-dependent manner. We and others have demonstrated that reactant alignment and correct base pairing at the Pols catalytic site are crucial structural features to fidelity. Here, we first used equilibrium molecular simulations to demonstrate that the local dynamics at the protein-DNA interface in the proximity of the catalytic site is different when correct vs incorrect dNTPs are bound to polymerase β (Pol β). Formation and dynamic stability of specific interatomic interactions around the incoming nucleotide influence the overall binding site architecture. This explains why certain Pols' mutants can affect the local catalytic environment and influence the selection of correct vs incorrect nucleotides. In particular, this is here demonstrated by analyzing the interaction network formed by the residue R283, whose mutant R283A has an experimentally measured lower capacity of differentiating correct (G:dCTP) vs incorrect (G:dATP) base pairing in Pol β. We also used alchemical free-energy calculations to quantify the G:dCTP →G:dATP transformation in Pol β wild-type and mutant R283A. These results correlate well with the experimental trend, thus corroborating our mechanistic insights. Sequence and structural comparisons with other Pols from the same family suggest that these findings may also be valid in similar enzymes.

Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions.

In drug discovery, the in silico prediction of binding affinity is one of the major means to prioritize compounds for synthesis. Alchemical relative binding free energy (RBFE) calculations based on molecular dynamics (MD) simulations are nowadays a popular approach for the accurate affinity ranking of compounds. MD simulations rely on empirical force field parameters, which strongly influence the accuracy of the predicted affinities. Here, we evaluate the ability of six different small-molecule force fields to predict experimental protein-ligand binding affinities in RBFE calculations on a set of 598 ligands and 22 protein targets. The public force fields OpenFF Parsley and Sage, GAFF, and CGenFF show comparable accuracy, while OPLS3e is significantly more accurate. However, a consensus approach using Sage, GAFF, and CGenFF leads to accuracy comparable to OPLS3e. While Parsley and Sage are performing comparably based on aggregated statistics across the whole dataset, there are differences in terms of outliers. Analysis of the force field reveals that improved parameters lead to significant improvement in the accuracy of affinity predictions on subsets of the dataset involving those parameters. Lower accuracy can not only be attributed to the force field parameters but is also dependent on input preparation and sampling convergence of the calculations. Especially large perturbations and nonconverged simulations lead to less accurate predictions. The input structures, Gromacs force field files, as well as the analysis Python notebooks are available on GitHub.

Computational Investigation of BMAA and Its Carbamate Adducts as Potential GluR2 Modulators.

Beta-N-methylamino-l-alanine (BMAA) is a potential neurotoxic nonprotein amino acid, which can reach the human body through the food chain. When BMAA interacts with bicarbonate in the human body, carbamate adducts are produced, which share a high structural similarity with the neurotransmitter glutamate. It is believed that BMAA and its l-carbamate adducts bind in the glutamate binding site of ionotropic glutamate receptor 2 (GluR2). Chronic exposure to BMAA and its adducts could cause neurological illness such as neurodegenerative diseases. However, the mechanism of BMAA action and its carbamate adducts bound to GluR2 has not yet been elucidated. Here, we investigate the binding modes and the affinity of BMAA and its carbamate adducts to GluR2 in comparison to the natural agonist, glutamate, to understand whether these can act as GluR2 modulators. Initially, we perform molecular dynamics simulations of BMAA and its carbamate adducts bound to GluR2 to examine the stability of the ligands in the S1/S2 ligand-binding core of the receptor. In addition, we utilize alchemical free energy calculations to compute the difference in the free energy of binding of the beta-carbamate adduct of BMAA to GluR2 compared to that of glutamate. Our findings indicate that carbamate adducts of BMAA and glutamate remain stable in the binding site of the GluR2 compared to BMAA. Additionally, alchemical free energy results reveal that glutamate and the beta-carbamate adduct of BMAA have comparable binding affinity to the GluR2. These results provide a rationale that BMAA carbamate adducts may be, in fact, the modulators of GluR2 and not BMAA itself.

Peptide Drug Design Using Alchemical Free Energy Calculation: An Application and Validation on Agonists of Ghrelin Receptor.

With recent large-scale applications and validations, the relative binding free energy (RBFE) calculated using alchemical free energy methods has been proven to be an accurate measure to probe the binding of small-molecule drug candidates. On the other hand, given the flexibility of peptides, it is of great interest to find out whether sufficient sampling could be achieved within the typical time scale of such calculation, and a similar level of accuracy could be reached for peptide drugs. However, the systematic evaluation of such calculations on protein-peptide systems has been less reported. Most reported studies of peptides were restricted to a limited number of data points or lacking experimental support. To demonstrate the applicability of the alchemical free energy method for protein-peptide systems in a typical real-world drug discovery project, we report an application of the thermodynamic integration (TI) method to the RBFE calculation of ghrelin receptor and its peptide agonists. Along with the calculation, the synthesis and in vitro EC50 activity of relamorelin and 17 new peptide derivatives were also reported. A cost-effective criterion to determine the data collection time was proposed for peptides in the TI simulation. The average of three TI repeats yielded a mean absolute error of 0.98 kcal/mol and Pearson's correlation coefficient (R) of 0.77 against the experimental free energy derived from the in vitro EC50 activity, showing good repeatability of the proposed method and a slightly better agreement than the results obtained from the arbitrary time frames up to 20 ns. Although it is limited by having one target and a deduced binding pose, we hope that this study can add some insights into alchemical free energy calculation of protein-peptide systems, providing theoretical assistance to the development of peptide drugs.

Multistate Method to Efficiently Account for Tautomerism and Protonation in Alchemical Free-Energy Calculations.

The majority of drug-like molecules contain at least one ionizable group, and many common drug scaffolds are subject to tautomeric equilibria. Thus, these compounds are found in a mixture of protonation and/or tautomeric states at physiological pH. Intrinsically, standard classical molecular dynamics (MD) simulations cannot describe such equilibria between states, which negatively impacts the prediction of key molecular properties in silico. Following the formalism described by de Oliveira and co-workers (J. Chem. Theory Comput. 2019, 15, 424-435) to consider the influence of all states on the binding process based on alchemical free-energy calculations, we demonstrate in this work that the multistate method replica-exchange enveloping distribution sampling (RE-EDS) is well suited to describe molecules with multiple protonation and/or tautomeric states in a single simulation. We apply our methodology to a series of eight inhibitors of factor Xa with two protonation states and a series of eight inhibitors of glycogen synthase kinase 3β (GSK3β) with two tautomeric states. In particular, we show that given a sufficient phase-space overlap between the states, RE-EDS is computationally more efficient than standard pairwise free-energy methods.

Application of an alchemical free energy method for the prediction of thermostable DuraPETase variants

Non-equilibrium (NEQ) alchemical free energy calculations are an emerging tool for accurately predicting changes in protein folding free energy resulting from amino acid mutations. In this study, this method in combination with the Rosetta ddg monomer tool was applied to predict more thermostable variants of the polyethylene terephthalate (PET) degrading enzyme DuraPETase. The Rosetta ddg monomer tool efficiently enriched promising mutations prior to more accurate prediction by NEQ alchemical free energy calculations. The relative change in folding free energy of 96 single amino acid mutations was calculated by NEQ alchemical free energy calculation. Experimental validation of ten of the highest scoring variants identified two mutations (DuraPETaseS61M and DuraPETaseS223Y) that increased the melting temperature (Tm) of the enzyme by up to 1 °C. The calculated relative change in folding free energy showed an excellent correlation with experimentally determined Tm resulting in a Pearson’s correlation coefficient of r = − 0.84. Limitations in the prediction of strongly stabilizing mutations were, however, encountered and are discussed. Despite these challenges, this study demonstrates the practical applicability of NEQ alchemical free energy calculations in prospective enzyme engineering projects.Key points• Rosetta ddg monomer enriches stabilizing mutations in a library of DuraPETase variants• NEQ free energy calculations accurately predict changes in Tmof DuraPETase• The DuraPETase variants S223Y, S42M, and S61M have increased TmGraphical

Binding Sites of Bicarbonate in Phosphoenolpyruvate Carboxylase.

Phosphoenolpyruvate carboxylase (PEPC) is used in plant metabolism for fruit maturation or seed development as well as in the C4 and crassulacean acid metabolism (CAM) mechanisms in photosynthesis, where it is used for the capture of hydrated CO2 (bicarbonate). To find the yet unknown binding site of bicarbonate in this enzyme, we have first identified putative binding sites with nonequilibrium molecular dynamics simulations and then ranked these sites with alchemical free energy calculations with corrections of computational artifacts. Fourteen pockets where bicarbonate could bind were identified, with three having realistic binding free energies with differences with the experimental value below 1 kcal/mol. One of these pockets is found far from the active site at 14 Å and predicted to be an allosteric binding site. In the two other binding sites, bicarbonate is in direct interaction with the magnesium ion; neither sequence alignment nor the study of mutant K606N allowed to discriminate between these two pockets, and both are good candidates as the binding site of bicarbonate in phosphoenolpyruvate carboxylase.

Adenine Methylation Enhances the Conformational Flexibility of an RNA Hairpin Tetraloop.

The N6-methyladenosine modification is one of the most abundant post-transcriptional modifications in ribonucleic acid (RNA) molecules. Using molecular dynamics simulations and alchemical free-energy calculations, we studied the structural and energetic implications of incorporating this modification in an adenine mononucleotide and an RNA hairpin structure. At the mononucleotide level, we found that the syn configuration is more favorable than the anti configuration by 2.05 ± 0.15 kcal/mol. The unfavorable effect of methylation was due to the steric overlap between the methyl group and a nitrogen atom in the purine ring. We then probed the effect of methylation in an RNA hairpin structure containing an AUCG tetraloop, which is recognized by a "reader" protein (YTHDC1) to promote transcriptional silencing of long noncoding RNAs. While methylation had no significant conformational effect on the hairpin stem, the methylated tetraloop showed enhanced conformational flexibility compared to the unmethylated tetraloop. The increased flexibility was associated with the outward flipping of two bases (A6 and U7) which formed stacking interactions with each other and with the C8 and G9 bases in the tetraloop, leading to a conformation similar to that in the RNA/reader protein complex. Therefore, methylation-induced conformational flexibility likely facilitates RNA recognition by the reader protein.

Resolving coupled pH titrations using alchemical free energy calculations.

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured . Here, we derive a formalism based on double free energy differences ( ) for quantifying the individual site values of coupled residues. As values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue s in practice. We demonstrate that our approach and a previously proposed microscopic formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the values of both uncoupled and coupled residues.

Quantitative Assessment of Energetic Contributions of Residues in a SARS-CoV-2 Viral Enzyme/Nanobody Interface.

The highly conserved protease enzyme from SARS-CoV-2 (MPro) is crucial for viral replication and is an attractive target for the design of novel inhibitory compounds. MPro is known to be conformationally flexible and has been stabilized in an extended conformation in a complex with a novel nanobody (NB2B4), which inhibits the dimerization of the enzyme via binding to an allosteric site. However, the energetic contributions of the nanobody residues stabilizing the MPro/nanobody interface remain unresolved. We probed these residues using all-atom MD simulations in combination with alchemical free energy calculations by studying the physical residue-residue interactions and discovered the role of hydrophobic and electrostatic interactions in stabilizing the complex. Specifically, we found via mutational analysis that three interfacial nanobody residues (Y59, R106, and L109) contributed significantly, two residues (L107 and P110) contributed moderately, and two residues (H112 and T113) contributed minimally to the overall binding affinity of the nanobody. We also discovered that the nanobody affinity could be enhanced via a charge-reversal mutation (D62R) that alters the local interfacial electrostatic environment of this residue in the complex. These findings are potentially useful in designing novel synthetic nanobodies as allosteric inhibitors of MPro.

Influence of Iodine Merz–Singh–Kollman Radius on the Calculated Charges and Hydration Free Energies of Iodinated Molecules

Empirical force field methods typically rely on point charges to describe the electrostatic interactions, which is problematic when anisotropy needs to be considered, as in the case of the electrostatic potential of covalently bound halogens that possess a positive site, termed [Formula: see text]-hole, surrounded by a large negative belt. To address this, an off-center point charge (extra point, EP) is usually placed at a given distance from the halogen to emulate the [Formula: see text]-hole and commonly used implementations are based on the restrained electrostatic potential (RESP) procedure to fit atomic charges, being one of the most used charge models. In this context, no specific Merz–Singh–Kollman (MK) radius for iodine is available in the literature, which is an essential parameter in the RESP fitting procedure. In this work, we explored the impact of the iodine MK radius on the obtained RESP charges for a set of 12 iodinated molecules. We verified that the relative root mean square (RRMS) values obtained with and without an EP kept decreasing with increasing radii for most compounds, thus impairing optimization using such a procedure. Nevertheless, the use of an iodine MK radius lower than 2 Å is not advisable since the RRMS kept decreasing considerably until this value was reached. Moreover, the performance of three iodine MK radii was studied with the estimation of the free energy of hydration ([Formula: see text] values using alchemical free energy calculations, which are particularly sensitive to the charges used. Despite the usage of different radii not leading to remarkable differences, our results indicate that using a value of 2.70 Å leads to lower mean absolute errors (MAE) and root mean squared error (RMSE) values when comparing the calculated with the experimental [Formula: see text] values.