Macrocycle Cell Permeability Measured by Solvation Free Energies in Polar and Apolar Environments.

The relation of surface polarity and conformational preferences is decisive for cell permeability and thus bioavailability of macrocyclic drugs. Here, we employ grid inhomogeneous solvation theory (GIST) to calculate solvation free energies for a series of six macrocycles in water and chloroform as a measure of passive membrane permeability. We perform accelerated molecular dynamics simulations to capture a diverse structural ensemble in water and chloroform, allowing for a direct profiling of solvent-dependent conformational preferences. Subsequent GIST calculations facilitate a quantitative measure of solvent preference in the form of a transfer free energy, calculated from the ensemble-averaged solvation free energies in water and chloroform. Hence, the proposed method considers how the conformational diversity of macrocycles in polar and apolar solvents translates into transfer free energies. Following this strategy, we find a striking correlation of 0.92 between experimentally determined cell permeabilities and calculated transfer free energies. For the studied model systems, we find that the transfer free energy exceeds the purely water-based solvation free energies as a reliable estimate of cell permeability and that conformational sampling is imperative for a physically meaningful model. We thus recommend this purely physics-based approach as a computational tool to assess cell permeabilities of macrocyclic drug candidates.

Introduction

Macrocycles are a potent new class of molecules for drug discovery.[1,2] Approximately 75% of disease relevant proteins still cannot be targeted, neither with small molecules nor with biopharmaceuticals.[3] A major portion of these yet undruggable targets are intracellular protein–protein interfaces (PPIs), including several notorious cancer-associated targets.[4] Biologics, such as antibodies, are the prime class of pharmaceuticals to target extracellular PPIs with uncontested specificities and affinities.[5,6] However, with a few exceptions, they are generally not able to cross through the cell membrane.[7] Small-molecule drugs, on the other hand, are extremely well-studied, and clear models and guidelines to achieve oral bioavailability and membrane permeability are well-established.[8] However, they mostly require deep apolar binding pockets to achieve the desired affinities and physiological effects, which are lacking in typical PPIs with extensive flat surface areas.[3] Macrocycles bridge these two medication strategies in terms of physicochemical and pharmacological features.[9−12]

Macrocyclic compounds have repeatedly been established as drugs without fulfilling all or even any of Lipinski’s rule of 5 for bioavailability of small-molecule drugs. Nevertheless, it has been shown that they can be designed to achieve cell permeability and even oral bioavailability.[3,13−18] As they are substantially larger than typical small molecules, macrocycles are able to target the characteristic shallow and broad surfaces of protein–protein interaction sites.[16,19−21] Furthermore, their proteolytic stability and thus bioavailability are increased due to the cyclic scaffold.[16,22,23] Compared to their non-cyclic analogues, the cyclization additionally decreases the entropic loss upon binding, which can enhance their binding affinity to magnitudes that are usually only achievable by biopharmaceuticals.[24−26] However, despite the continuous advancement in experimental strategies, the synthesis of macrocyclic compounds is still challenging.[27−30] Reliable computational tools to identify and optimize promising scaffolds are thus paramount for the efficient design of macrocyclic drugs.[31−36]

Substantial scientific efforts in this field have led to a fast-growing number of theoretical methodologies for characterizing physicochemical properties of macrocycles.[9,37,38] A major aspect of these approaches is concerned with the development and testing of conformational sampling algorithms suitable for macrocyclic molecules.[35,36,39] The development of specialized strategies for conformer generation is imperative as the conformational restraints introduced by the ring closure entails unique structural characteristics to this compound class, which are generally not captured by conventional conformer generators.[34,40] Furthermore, cyclization can also induce a strain energy within the ring leading to high energetic barriers between relevant conformational states.[41] Numerous of the proposed sampling algorithms for macrocycles are force-field based.[36,38,42−44] While classical molecular dynamics (MD) simulations have often failed to overcome the energetic barriers between the diverse conformational states of macrocycles within a feasible simulation time, several enhanced sampling strategies have been shown to capture structurally accurate ensembles.[31] These more sophisticated sampling techniques, such as replica exchange MD,[45] multicanonical MD,[43] metadynamics,[46,47] and accelerated MD,[31] allow comprehensive and efficient profiling of the conformational space of macrocycles.

As described above, a particularly intriguing feature of macrocycles is their ability to cross the cell membrane.[15,48] While high passive membrane permeability has been demonstrated for a multitude of macrocycles, not all macrocyclic scaffolds are inherently membrane-permeable.[17,49] In order to achieve permeability, macrocyclic compounds have to also balance an intricate interplay of physicochemical properties to ensure solubility in the polar extra- and intracellular environments as well as within the mostly apolar membrane.[50,51] The surprisingly high permeability of these large molecules is commonly explained based on solvent-dependent conformational rearrangements ( e.g., cyclosporin A).[38,43,52] The fundamental idea is that macrocycles with high passive membrane permeabilities are able to adapt to different solvent polarities via a population shift in their conformational ensembles: In aqueous solution, the most favorable conformational state is “open” with polar groups turned outward to interact with the polar solvent. Upon entering a less polar environment, the conformational ensemble shifts toward a “closed” conformational state. Here, polar groups are turned inward increasing the number of intramolecular hydrogen bonds, while apolar groups rearrange to maximize the apolar surface area. Several comprehensive studies, including extensive enhanced sampling and Markov state modeling, have fostered this hypothesis and also extended it to more generalized multistate models.[53−55]

To quantify conformational preferences in polar and apolar environments, most of these studies perform simulations in water and an apolar solvent, such as chloroform.[18,38,43] In previous studies, atomistic models for membranes led to promising results in estimating permeability of small molecules. However, this approach is still challenging and computationally costly.[56] Despite its striking simplicity, the approximation of a membrane by organic solvents has repeatedly demonstrated its suitability to estimate cell permeability in a multitude of approaches.[17,51,57,58]

In this study, we use accelerated MD (aMD) simulations[59,60] to capture the diverse conformational ensembles of a series of six macrocyclic compounds in water and chloroform (Figure ).[61] We have previously shown the reliability of aMD simulations in characterizing the structural ensemble and thermodynamic quantities of macrocycles consistent with experiments.[31] Here, we perform aMD simulations to profile and quantify solvent-induced shifts of ensemble populations. Furthermore, we track how structural rearrangements translate into changes in surface properties by calculating solvation free energies with grid inhomogeneous solvation theory (GIST).[62,63] The fundamental idea of this approach is to estimate thermodynamic solvation properties of the solvent in the vicinity of the solute by tracing the solvent distribution in MD simulations. We have demonstrated previously that GIST solvation free energies can be used to describe surface hydrophobicity.[64,65] Furthermore, we have recently reimplemented this approach on the GPU achieving superior computational performance. Although, compared to the calculation of polar surface area values, GIST analysis is still computationally noticeably more demanding and estimating surface polarity via solvation free energies offers several convincing advantages. First, GIST accounts for non-additive effects, for example, polar atoms, which are solvent-exposed but form intramolecular hydrogen bonds, show fewer interactions with the solvent and thus contribute less to surface polarity. Second, the latest reimplementation of the GIST algorithm introduces the possibility to calculate solvation free energies in chloroform. Hence, while previous GIST studies estimate differences in hydration free energy referenced to the solute in vacuum,[62,64,66−69] we now are able to compare differences between solvation in water and in chloroform.[70] In a preceding study, we demonstrated the accuracy of the approach in estimating partition coefficients, i.e., differences in solvation free energies, between water and chloroform for a set of rigid small molecules.[70]

Figure 1

Series of macrocyclic model systems.[61] The model compounds share a common ring scaffold but vary in their side chains. These modifications were specifically designed to achieve different conformational preferences and cell permeabilities. Reprinted with permission from ref (61). Copyright 2018 American Chemical Society.

For the present study, we extend this approach by further incorporating conformational aspects into the calculations of solvation free energy differences. By combining information on the state populations in chloroform and water with the respective solvation free energies, we derive an estimate of the transfer free energy of each compound. We evaluate the reliability of this approach against the experimentally measured cell permeabilities and analyze contributions of conformational sampling and solvation free energies from both solvents.

Theory and Methods

Grid Inhomogeneous Solvation Theory

GIST calculates thermodynamic properties of a solvent around a solute. Calculation of the free energy around the solute can be split into an energetic and an entropic part. These two parts are then calculated individually (eq ). Here, we aim at a short summary of the fundamental theory of GIST. For a comprehensive overview of GIST[62,63,71,72] as well as of the underlying concepts from Lazaridis et al.,[73,74] we want to refer the reader to the original literature.

For the calculation of the energetic contribution, the force field of the simulation is used. At each grid voxel where a water molecule is found, the total energy of this molecule is calculated. This total energy consists of two parts (eq ), the solvent–solvent interaction energy and the solute–solvent interaction energy (ΔE()). To avoid double counting, the total solvent–solvent energy is divided by two, and we denote the result of this division as ΔE().[62]

Finally, the entropic contribution only considers the two-body term. This two-body term is approximated via a nearest neighbor method, which estimates the translational and the orientational contributions to the entropy together (ΔSsix). The nearest neighbor is calculated by an l2 norm of the distance in translational and orientational space, resulting in eq . Finally, the translational distance is calculated as a simple Euclidean distance and the orientational distance is calculated as a quaternion distance.

Transfer Free Energies from GIST Calculations in Water and Chloroform

The transfer free energy between two different solvents can be readily computed from GIST using the thermodynamic cycle.[70] The transition from chloroform into water can be partitioned into the transition from chloroform into vacuum and then from vacuum into water. The first transition is the solvation free energy of the compound in chloroform. The second transition is simply the negative hydration free energy of the compound. Both values can be calculated using GIST.

For the calculation of the solvation free energies, GIST analyses were performed on simulation trajectories of the various macrocycles. For this analysis, the reference density for water was set according to the values in the AMBER manual (0.0329 Å−3),[75] and the reference density for chloroform was set to the same value as found by Kraml et al. (0.00768 Å−3).[70]

The values for the solvation free energy were then derived from the GIST calculations. In a first step, the reference solvent–solvent energy was subtracted. For the TIP3P water model, the value present in the AMBER manual was used, and for chloroform, the value found by Kraml et al. was used.[70] In the second step, the two energy contributions were summed up, following eq . Then, the entropic contribution was subtracted, following eq , to yield the solvation free energy in the respective solvent.

For the calculation of the transfer free energy (ΔAtransfer), the two solvation free energies were subtracted from each other, following eq .

Studied Series of Macrocyclic Compounds

We test the reliability of our approach on six macrocyclic molecules introduced by Tyagi et al., which were inspired by natural products.[61] This series was specifically designed to investigate molecular determinants of passive membrane permeability, in particular, the role of shielding NH−π interactions. Within their comprehensive work, Tyagi et al. provide high-quality experimental data including cell permeabilities, log POW values, and a crystal structure. Despite their high similarities, the six macrocycles in this series, 1a to 1f, show clear distinction in their passive membrane permeability (Figure ). Macrocycle 1e clearly exhibits the lowest permeability followed by the macrocycles 1a, 1c, and 1d, which are nearly identical in their permeabilities. The highest permeabilities are observed for the macrocycles 1b and 1f. Hence, the permeability distribution for these six is not spread evenly but rather represents three to four clusters. Nevertheless, the high quality and consistency of the available experimental data render these molecules an ideal test set for the presented study.

Structure Preparation and Simulation Setup

We performed our calculation using a set of macrocyclic compounds published by Tyagi et al.[61] The single available crystal structure (compound f in the terminology of Tyagi et al.; CCDC identifier 1853494) was used directly, while all other studied compounds were modeled based on this structure with the molecular operating environment (MOE).[76] All compounds were parameterized using the AMBER ff14SB[77] and GAFF[78] force fields. Missing parameters were derived with the antechamber module of AmberTools 18.[79] We assigned partial charges using Gaussian16[80] and the RESP[81] procedure using the HF/6-31G* basis set. Subsequently, topology and coordinate files were generated with the tLEaP module of AmberTools 18.[75] All compounds were solvated in a cubic box of TIP3P water[82] and chloroform[83] with a minimum wall distance of 12 Å for both solvents. Before production simulations, all systems were equilibrated following an extensive heating and cooling protocol.[84] The conformational space was sampled with the accelerated molecular dynamics (aMD) framework, as implemented in AMBER 18, using the dual boost approach.[59,85−87] Boosting parameters were obtained from short cMD simulations following a procedure adapted from Pierce et al.[60] For each macrocycle, 10 aMD simulations with random starting velocities were performed in both solvents to maximize conformational sampling. The simulations were run for 100 ns each, resulting in an aggregate simulation time of 1 μs/system. We used a Langevin thermostat[88] with a collision frequency of 2 ps–1 to keep the system at 300 K, together with a Berendsen barostat[89] with a relaxation time of 2 ps to keep the system at atmospheric pressure. All bonds involving hydrogen were constrained with the SHAKE algorithm, which enabled the use of a 2 fs time step for the water simulations.[90] Chloroform simulations were run with a 1 fs time step. All simulations were carried out with the GPU implementation of the PMEMD module in AMBER.[91,92]

Clustering and GIST Calculations

To determine representative structures and their population in each solvent, we performed a hierarchical cluster analysis based on the polar surface area (PSA) using cpptraj from AmberTools 18.[93] We estimate the PSA as the solvent-accessible surface area of polar heavy atoms, oxygen and nitrogen, for the concatenated trajectories in water and chloroform for each model system. We then applied the implemented average linkage clustering algorithm with a cutoff distance of 12 Å2 on the resulting data. This clustering strategy resulted in three to five representative structures per macrocycle. To reweight the respective population of each conformational cluster, we adapt a strategy proposed by Miao et al. in which we approximate the Boltzmann factor with a Maclaurin expansion up to the 20th order.[94,95]

To perform subsequent GIST analysis, we solvated each representative structure in a cubic box of TIP3P waters and chloroform molecules, ensuring a minimum distance of 20 Å for each solute atom to the box faces. We equilibrate the solvent as described above. Subsequently, we performed restrained simulations with a weight of 100 kcal/(mol·Å2) to fix the solute heavy atoms. Furthermore, all bonds involving hydrogen atoms were restrained following the SHAKE algorithm.[90] As described above, we employed the Langevin thermostat with a collision frequency of 2 ps–1 to ensure a temperature of 300 K and used the CUDA implementation of the particle mesh Ewald MD (pmemd.cuda) module of AMBER 18 to perform simulation in the NVT ensemble.[91,92] For each macrocycle/solvent combination, we performed 100 ns of restrained MD simulations where frames were collected every 2 ps, resulting in 50,000 frames per trajectory.

Conformational Space Characterization

To characterize the captured conformational space of each compound in water and chloroform, we performed principal component analysis (PCA) based on the dihedrals along the ring atoms using cpptraj.[93] To represent the common and deviating structural preferences in each solvent, we calculated the principal components based on the concatenated trajectories. In order to compare the structural data in water and chloroform, we project the respective data on the two first eigenvectors (with the highest eigenvalues) of the matrix with the combined features. To retrieve the unbiased population from the aMD simulations, we applied a Boltzmann reweighting scheme described above[94] (Figure S1). For the two macrocycles with the highest and lowest permeability, i.e., 1f and 1e, we performed a density-based clustering of the PC space to visualize structural differences corresponding to the free energy minima (Figures S3 to S5). Furthermore, we profiled the conformational space using a more global measure based on ratios of the principal moments of inertia (Figure S6).

As an additional measure of solvent-dependent structural rearrangements, we calculated PSAs, i.e., solvent-accessible surface areas of oxygen and nitrogen atoms. In order to remove the bias introduced by the aMD approach, we again reweighted the individual distributions as described above (Figure S7).

Results

We characterize the captured conformational space of six macrocyclic compounds using PCA (Figure S2). Projecting the structural data captured in each solvent onto the combined PCA space allows a direct comparison of structural preferences. In Figure , we depict the conformational ensembles of the macrocycles 1e and 1f, which show the largest difference in cell permeability. For both systems, we clearly observe a solvent-dependent shift of ensemble populations. While the covered conformational space is similar in chloroform and water, the maxima in population are clearly shifted. These maxima in the probability density directly translate into minima in free energy to which we will further refer as favorable conformational states. For macrocycle 1e, we find four highly populated areas in chloroform of which the most favorable conformational state is located around [2, −0.5] in the PCA space shown in Figure A. On the other hand, in water, the same macrocycle populates seven to eight distinguishable areas and the most favorable conformational states are found in areas around [−1.8, −1] and [−1, −0.5] (Figure . The conformational space of macrocycle 1f in chloroform is clearly less restricted than that of 1e, showing a large number of highly populated areas in the space spanned by PC1 and PC2 (Figure C). The ensemble of the same macrocycle in water shows significantly fewer favorable conformational states and is overall shifted toward areas with more positive PC2 values (Figure D).

Figure 2

Solvent-dependent ensemble shift of macrocycles 1e and 1f. The structural ensembles of macrocycles 1e (A, B; top) and 1f (C, D; bottom) in chloroform (green) and water (blue) are projected onto the first two PC eigenvectors and color-coded according to their reweighted populations.

Furthermore, we profile distinctions between the ensembles captured in the polar and apolar environments in terms of intramolecular hydrogen bonds (IMHBs) and surface properties. It has been shown before that the structural differences of several peptidic macrocycles in varying solvents relate to changes in the pattern of IMHBs and in the polar surface area (PSA).[43] In Figure , we show the ensemble distributions of both descriptors for the most permeable macrocycle in our series, 1f, compared to the least permeable compound 1e. We clearly find that, for both systems, the number of IMHBs is higher in chloroform, while the polar surface area shifts toward smaller values in the apolar environment. However, for the least permeable macrocycle 1e, the average number of IMHBs observed during the simulation is significantly smaller than for 1f. Furthermore, we find that the distribution of PSA in chloroform is clearly skewed toward higher values for 1e, which is not observed for 1f. This observation suggests that 1e is not able to bury or shield its polar moieties as well as 1f.

Figure 3

Ensemble distributions of conformational descriptors in water and chloroform. (A, C) Histograms of the number of intramolecular hydrogen bonds (IMHB), and (B, D) the polar surface areas (PSAs) of macrocycles 1e and 1f are depicted in blue for simulations in water and in green for simulations in chloroform.

To achieve a more detailed analysis of the surface properties, we calculate solvation free energies using GIST. To do so, we determine representative structures using the PSA-based clustering strategy described in the Theory and Methods section. From the GIST solvation free energies in chloroform and water, we can then calculate transfer free energies ΔAtransfer for each representative structure. In Figure , we show the cluster populations and representative structures for the macrocycles with the highest and lowest membrane permeability, 1f and 1e, respectively. As described above, we find a more widespread PSA distribution in both solvents for macrocycle 1e than for 1f (Figure B,D). This trend translates into a higher number of clusters and a slightly broader distribution in terms of cluster populations for macrocycle 1e (Figure A). Furthermore, the cluster with the highest probability in water and chloroform is the same one for 1e. This is an additional result of the substantial overlap of 1e’s PSA distributions in water and chloroform. For macrocycle 1f, on the other hand, the ensemble probabilities shift from cluster 1, being most favorable in chloroform, toward cluster 2 when simulated in water (Figure C). Comparing representatives for the highest populated clusters of 1f in water and chloroform, we find that the conformation favored in chloroform shows a higher number of intramolecular hydrogen bonds (cyan dotted lines) (Figure D).

Figure 4

Conformational preferences in water and chloroform. (A, C) Cluster populations of macrocycles 1e (top row) and 1f (bottom row) are depicted in blue for simulations in water and in green for simulations in chloroform. (B, D) Representative structures from the respective clustering of macrocycles 1e and 1f.

To obtain a quantitative measure of the transfer free energy for each compound, we use the solvent-dependent populations of all representative structures and weight the transfer free energies accordingly, resulting in an ensemble average. In Figure A, we compare the resulting transfer free energies with the experimentally determined cell permeabilities and report a striking value of 0.92 for the Pearson correlation coefficient r. We find that the predicted ordering is perfectly in line with the experiment, as indicated by the Spearman correlation coefficient ρ of 1.0. If we follow the exact same sampling strategy but solely consider the solvation free energy in water, we obtain r = 0.71 and ρ = 0.49 (Figure B). Comparing the solvation free energy of chloroform alone with experimental cell permeabilities, we find similar agreement as quantified by r = −0.71 and ρ = −0.77 (Figure S10).

Figure 5

Impact of conformational sampling and contributions from both solvents on prediction accuracy. The top row depicts the cell permeabilities compared to (A) transfer free energies and (B) water solvation free energies considering the conformational ensemble of each macrocycle. The bottom row shows the results of the same calculations when only a single conformation is included (C - transfer free energy, D - water solvation free energy).

Furthermore, we obtain substantially lower agreement with the experiment if we do not consider the conformational variability of each macrocycle but only use the single starting structure of each compound for the GIST calculations (Figure C,D). By calculating transfer free energies from a single structure of each compound, we find r = 0.50 and ρ = 0.49 (Figure C). Calculating water solvation free energies of the single conformations lowers the correlation further, resulting in r = 0.29 and ρ = 0.43 (Figure D). Considering only contributions from chloroform solvation free energies of the single conformations also results in moderate agreement with the experiment with r = −0.58 and ρ = −0.43.

Discussion

Macrocyclic drugs could become a powerful alternative to biopharmaceutical drugs since they show superior membrane permeability while offering similar binding affinities.[25] However, their membrane permeability strongly depends on the interplay between surface polarity and conformational preferences.[51,96] In this study, we incorporate both of these aspects in our calculation of transfer free energies using ensemble-averaged solvation free energies.

We perform accelerated molecular dynamics simulations to capture the conformational diversity of six macrocycles in water and chloroform (Figure S1). We find that, for all model systems, the conformational space in water and chloroform overlaps to a varying extent. However, the change in solvent polarity consistently alters the population of the captured conformational states. Hence, we find a solvent-dependent shift of ensemble probabilities. This observation is consistent with current literature describing similar behavior for a broad range of peptidic and non-peptidic macrocycles.[38,43,57]

Comparing the conformational landscapes for the two macrocycles 1e and 1f, which vary most in their cell permeabilities, implies that the core scaffold of 1e has to undergo major structural rearrangements upon membrane penetration, which could translate into a significant cost in terms of free energy. The ring atoms in 1f (highest permeability) on the other hand are likely to pre-organize in conformations favorable in both solvents, which facilitates membrane crossing. This finding is perfectly in line with a model proposed by Witek et al. on the mechanism of cell permeability of macrocycles.[38,97] The authors propose that the traditional hypothesis of a macrocycle switching from one (open) conformation that is favorable in water to another (closed) conformation favorable in chloroform is too simplistic. Based on exhaustive sampling, they suggest a more generalized model, which includes one or more congruent conformational states. These congruent conformational states are significantly populated in both solvents and allow the compound to pass through the membrane.

However, the cost of conformational rearrangements is only one contribution to the transfer free energy, which determines the water membrane partition coefficient. While the PCA projections capture differences based on dihedral conformations, these analyses are not designed to reflect distinctions of surface properties. Since surface polarity is decisive for the solvation free energies involved in the membrane crossing, we calculate the PSA for the ensemble of each compound in both solvents (Figure S7). We find that, for all compounds, the PSA distribution is shifted toward smaller values in chloroform. This observation is consistent with several studies from Kihlberg and co-workers on the conformational ensembles of macrocycles in different environments.[96] Based on results from NMR, X-ray crystallography, and computational sampling, these studies report an increase of intramolecular hydrogen bonds and a decrease of the polar surface area for macrocyclic systems in an apolar environment compared to water.[51,57]

To assess structural aspects associated with the observed variations in surface polarity, we retrieve representative conformations using a PSA-based clustering (Figure A,C). For macrocycle 1f, we find the expected trend. In the PSA-based clustering of macrocycle 1e, however, the predominant conformational state of the 1e ensemble in chloroform does not show the typical features promoting cell permeability. This qualitative observation could be another contribution to the lower cell permeability of macrocycle 1e compared to 1f.

The results from the PSA-based clustering might appear counterintuitive compared to the trend captured with PCA. Hence, we emphasize that these two analyses focus on different aspects of conformational dynamics. The dihedral PCA depicted in Figure captures the structural preferences at the core of the macrocycles, i.e., the ring atoms. However, this representation does not consider motions of the macrocycles’ side chains. In Figures S2 to S5, we highlight that distinctly separated minima in the dihedral PCA space can be very similar in terms of heavy-atom RMSD and surface polarity. In contrast, the PSA-based clustering shown in Figure identifies structures based on their difference in surface polarity. Hence, by definition, we retrieve representative structures with clearly distinct surface polarity (Figures S8 and S9), which is dictated by the flexible side chains. Thus, the analyses presented in Figures and 4 elucidate complementary aspects of the underlying mechanism and molecular determinants of membrane permeability.

For a thorough quantification of the properties on the surface, we then compute the transfer free energies for each conformation using GIST. By calculating the respective ensemble averages, we find a striking agreement of GIST transfer free energies with cell permeability, resulting in a Pearson correlation coefficient of r of 0.92 and a Spearman correlation coefficient ρ of 1.0 (Figure A). This result becomes even more intriguing considering that the experimental log POW, i.e., the octanol–water partition coefficient, for the same compounds shows an r of 0.89 and ρ of 0.77 with cell permeabilities (Figure S11). The proposed workflow hence is comparable with experimental log POW values in accuracy and could thus tremendously reduce the costs at the stage of macrocycle design and optimization.

Furthermore, we profile the impact of individual aspects of the applied methodology. As described in the Introduction, we have previously shown that water solvation free energy relates to a molecule’s hydrophobicity. Therefore, we tested for the benefits of incorporating contributions from chloroform solvation. In Figure B, we use the captured conformational ensembles and the respective populations in both solvents to calculate the ensemble average of the water solvation free energy of each macrocycle. Thus, we apply the same workflow for Figure A but only consider solvation free energies from water (Figure B) or chloroform (Figure S10). The results are clearly inferior to the transfer free energies that consider the contributions from both water and chloroform solvation.

This finding is in line with our latest work, in which we introduce transfer free energies to estimate partition coefficients between water and chloroform. In this preceding study, we deliberately only considered rigid small molecules to avoid inaccuracies resulting from insufficient sampling. Here, we describe transfer free energies as an ensemble property. Consequently, we also benchmark whether conformational sampling has a beneficial impact. In Figure C,D, we again calculate transfer free energies and water solvation free energies, yet we only consider a single conformation of each macrocycle. We consistently find that ensemble averages lead to substantially higher agreement with the experiment than single conformations.

We want to note that, in recent years, several workflows have been established that achieve exhaustive and reliable sampling of macrocycle conformational ensembles.[32,35,36,38−40] Here, we perform aMD simulations as their computational demand is comparable to classic MD simulations while the conformational sampling is significantly enhanced.[98,99] Furthermore, we have previously benchmarked the high reliability of this technique in providing conformational state populations.[31] However, the presented GIST calculations are independent of the aMD approach and can be applied to ensembles from any sampling technique. We thus deem the proposed workflow as a highly generalizable and reliable tool to estimate cell permeabilities of macrocyclic drug candidates.

Conclusions

We present an approach to estimate cell permeabilities of macrocycles where we use aMD simulations and GIST to derive ensemble-averaged transfer free energies. We profile the conformational space in water and chloroform of a set of six macrocycles, which were designed to exhibit varying cell permeabilities and structural preferences.[61]

By applying our proposed approach, we find a remarkable correlation to the experimental membrane permeability with an r of 0.92 and ρ of 1.0, which is comparable with experimental log POW values. Furthermore, we find that exhaustive conformational sampling is an indispensable step to retrieve physically meaningful predictions for the physicochemical properties of macrocycles. Additionally, we show that incorporation of the contributions from chloroform considerably increases the reliability of our descriptor compared to the solely water-based approach.

Thus, we highlight the significance of reliable conformational sampling in macrocycle design. Additionally, we demonstrate the benefits of our recent re-implementation of the GIST algorithm, which now also allows us to employ chloroform as a solvent. We thus demonstrate a powerful workflow that provides extremely reliable estimates of macrocycle membrane permeabilities to enhance their design and optimization process.

The latest GPU implementation of GIST (GIGIST) is available free of charge from the github page of the Liedl Lab (https://github.com/liedllab/gigist.git).