Conformation of the Macrocyclic Drug Lorlatinib in Polar and Nonpolar Environments: A MD Simulation and NMR Study.

The replica exchange molecular dynamics (REMD) simulation is demonstrated to readily predict the conformations of the macrocyclic drug lorlatinib, as validated by solution NMR studies. In aqueous solution, lorlatinib adopts a conformer identical to its target bound structure. This conformer is stabilized by an extensive hydrogen bond network to the solvents. In chloroform, lorlatinib populates two conformers with the second one being less polar, which may contribute to lorlatinib's ability to cross cell membranes.

Introduction

Macrocycles are of particular interest in medicinal chemistry. They often display semirigid conformations, which not only allow them to match one of their major conformations to their target proteins, thereby achieving an increase in both binding affinity and selectivity as compared to a nonmacrocyclic analogue, but also facilitate their membrane penetration.[1−3] Despite having molecular weights that often exceed Lipinski’s rule of five,[4] macrocyclization may provide drug-like properties,[1,2,5,6] such as significantly improved metabolic stability and good oral bioavailability.[7,8] In the past decades, more than 100 macrocyclic drugs and clinical candidates have been marketed or studied in drug discovery programs,[1] most of which are natural products, such as erythromycin, cyclosporin A, rifampicin, spiramycin, and sirolimus.[1,9] Synthetic macrocycles are attracting an increasing interest due to their molecular diversity and complexity.[10] Several recently approved drugs including vaniprevir,[11] danoprevir,[12] and lorlatinib[13] show the vast potential of macrocycles in drug discovery.[6,14] In particular, macrocycles hold great promise for challenging targets such as protein–protein interactions (PPIs).[15,16] Recently macrocyclic drugs have been reported to target a variety of PPIs, such as MDM2/p53,[17] BCL-2,[18] and IL-17A/IL-17RA.[19]

The target bound conformer of macrocycles is expected to be present in their solution ensembles to a measurable extent.[20−25] Moreover, conformational dynamics involving bond rotations may result in dynamic shielding or exposure of their functional groups in response to changes in the polarity of the environment.[26−28] Such “chameleonic” behavior has been considered to be important for the penetration of biological membranes and the bioavailability of macrocycles as drugs.[26−29] Experimental methods, for example, NMR spectroscopy and X-ray crystallography, can be used to accurately determine the solution and solid state conformations of macrocycles. However, they not only require real samples but also are expensive and time-consuming.

Recent investigations[30−38] have revealed a large and growing interest in the theoretical prediction of the conformational landscape of macrocycles. Algorithms used include geometric methods (distance geometry,[30] inverse kinematics,[34] etc.), Monte-Carlo torsional sampling,[33] and molecular dynamics (MD)-based sampling methods (LowModeMD,[32] accelerated MD,[39] CoCo-MD,[40] etc.). In addition, a comparison of the macrocycle sampling performance of MacroModel,[33] LowModeMD,[32] and a distance-geometry-based algorithm in OMEGA[41] was presented in our previous work.[42] In general, conventional MD (cMD) simulations are considered to be able to capture structural interconversions at the atomic level.[43] However, cMD simulations of biological systems tend to get trapped in a large number of local minima. To overcome the low sampling efficiency, the replica exchange molecular dynamics simulation (REMD), which allows the exchange of noninteracting replicas of the system at different temperatures, has been applied to varepsilon-cyclodextrin and cyclic peptides.[44,45] Despite recent advances, only a few studies have begun to investigate how the properties of the environment (e.g., its polarity, hydrogen bond donor and acceptor capability) affect the conformational preferences of macrocycles and how these effects can be predicted.[27,42] As body fluids are aqueous in nature while the bilipid layer of cell membranes is nonpolar, advances that allow prediction of macrocycle conformations in different environments would constitute a major step toward prediction of their pharmaceutical properties.

Lorlatinib (PF-06463922, Pfizer) is a selective and reversible 12-membered macrocyclic drug approved for treatment of nonsmall cell lung cancer (NSCLC), possessing good effectiveness, safety, and bioavailablilty.[46,47] As compared to linear analogues, the macrocyclic framework reduces the number of rotatable bonds, leading to a more compact bioactive conformer when bound to its target proteins, ALK (anaplastic lymphoma kinase) and ROS1 (c-ros oncogene1 kinase).[48] X-ray structures revealed that lorlatinib adopts the same bioactive conformer in the two target protein crystals.[13,46,49] However, high membrane and blood brain barrier permeability[13] indicates the possibility that lorlatinib may display molecular chameleonicity. It is therefore an excellent model to investigate conformational sampling of macrocycles in solvents of different polarities.

In this study, we evaluated the ability of REMD simulations to predict the conformation ensembles of macrocycles in polar and nonpolar environments using lorlatinib as an example. Aqueous solution was used to mimic the body fluids and chloroform to mimic the hydrophobic interior of the cell membrane. Chloroform was used as it has a dielectric constant (ε = 4.8), close to that determined for a lipid bilayer (ε = 3.0),[50] and as it is the most common nonpolar solvent used for NMR spectroscopy. Our simulations revealed that the solvent has an essential impact on the conformational distribution of lorlatinib, exerted via intermolecular interactions. The simulation results were validated against solution NMR experiments, indicating that REMD simulation is capable of predicting the ensemble of lorlatinib in both environments.

Results

REMD Predictions

By periodical exchange attempts of neighboring replicas at different temperatures, REMD simulations are expected to efficiently explore larger conformational space than cMD, which might remain in local minima. Using REMD, lorlatinib was simulated in two solvents, viz., CHCl3 and a H2O/DMSO (6:4) mixture, respectively, at temperatures ranging from 300 to 600 K. The H2O/DMSO mixture was used due to lorlatinib’s poor aqueous solubility.

We calculated the one-dimensional free-energy profile as a function of heavy-atom RMSD (Figure S1), which revealed an excellent convergence of the simulations in both solvents, as presented by the free-energy profiles at three periods of time: 80, 90, and 100 ns. Only one conformation, P1, was found in aqueous solution, while two conformers were predicted to exist in CHCl3 (P1 and P2). P1 is similar to the crystal structure of lorlatinib complexed with ALK (PDB ID: 4CLI; Pc, RMSD 0.27 Å), and P2 has an RMSD of 1.2 Å to Pc (Figure S2). This suggests that lorlatinib binds to its target proteins (ALK and ROS1) in the conformation it adopts in aqueous solution.

As the RMSD of the overall structure might hide local motion of the structural motifs, lorlatinib was divided into two parts, Part 1 and Part 2 (Figure A). RMSDs of these two parts were calculated for all the snapshots taken from the simulation at 300 K. While the RMSD of Part 1 (RMSD_part1) mainly describes the ring flip motion originating from rotation about the bond between the pyridine and pyrazole rings, RMSD of Part 2 (RMSD_part2) is sensitive to the benzene ring flip motion originating from rotation about the single bonds in Part 2 and rotation about the amide bond (Figure B). Two-dimensional free-energy landscapes were thus calculated with the two RMSDs as reaction coordinates. The population of P1 in CHCl3 was calculated to be 80%, while that of P2 is 20% (Figure B,D). RMSDs of P1 to Pc are 0.25 Å both for Part 1 and Part 2, further supporting that P1 and Pc are identical conformations. P2 has RMSDs of ∼1.5 and ∼0.9 Å to Pc for Part 1 and Part 2, respectively, correlating with the major structural change as compared to Pc. Furthermore, P1 is ∼11.3 kJ/mol more stable than P2 in terms of the calculated free energy. The free-energy barrier of P1 changing to P2 is ∼22.6 kJ/mol, indicating that P1 and P2 interconvert in chloroform. According to the minimum free-energy pathway (Figure D), the conformational shift from P1 to P2 is initiated by an inward flip of both Part 1 and Part 2, followed by Part 1 flipping further inward; indicating the importance of flip motions of the pyridine ring. Only conformer P1 (100%) was predicted to exist in H2O/DMSO (6:4) (Figure C).

Figure 1

Simulated conformations and free-energy profiles of lorlatinib. (A) Lorlatinib and its two parts divided for RMSD analysis; (B) One-dimensional free-energy profile of the structural transition from P1 to P2 in CHCl3; (C and D) Two-dimensional free-energy landscapes of the conformation distribution of lorlatinib in H2O/DMSO and in CHCl3, respectively. The minimum free-energy pathway is shown as a gray line.

NMR

Lorlatinib was analyzed in D2O/DMSO-d6 (6:4) and CDCl3. Experimental population-averaged distances were determined by acquisition of nuclear Overhauser enhancement (NOE) build-ups at 800 MHz and by conversion of the initial build-up rates into interproton distances (Tables S1, S2, and S3). A theoretical ensemble covering the entire available conformational space was generated by unrestrained Monte Carlo conformational search using the GB/SA water and chloroform solvation models. Following redundant conformation elimination, conformations from all individual searches were combined and used as theoretical inputs for the NAMFIS analyses (Table S4). The NAMFIS algorithm has previously been successfully applied for the description of solution ensemble of various flexible macrocycles.[25,42,51,52] Solution ensembles were determined by varying the probability of each conformation and fitting the back-calculated distances for each computationally predicated conformation to the experimentally determined population-averaged distances derived from NMR studies for the two solvents (D2O/DMSO-d6 and CDCl3) (Tables S6 and S7).

These analyses (Table S7 and Figure S7) resulted in one conformer (1) in D2O/DMSO-d6 (Figure A) and two conformers in CDCl3 (1 and 2) (Figure B). In CDCl3, 1 is also the major conformer (76%), while 2 is the minor one with a population of 24% (Table S7). Conformation 1 showed higher overall similarity (RMSD 0.35 Å) than 2 (RMSD 1.02 Å) to the crystal structure Pc (Table S8). This experimental observation is in excellent agreement with the REMD prediction (Figure ). The overall difference between 1 and P1 (RMSD 0.15 Å), and 2 and P2 difference (RMSD 0.17 Å) is very small.

Figure 2

Solution conformations of lorlatinib in D2O/DMSO-d6 (blue) and CDCl3 (pink) aligned with the predicted conformations (gray). (A) 1 vs P1 in aqueous solution, (B) 1 vs P1 (left), and 2 vs P2 (right) in chloroform.

Discussion

The ability of REMD simulations to accurately predict the solution ensembles of lorlatinib in both polar and nonpolar environments was confirmed by NMR conformational studies. The different conformational spaces populated in different solvents reveal that the solvent plays an important role in determining lorlatinib’s conformations.

To understand the observed molecular chameleonicity on an atomic level, a detailed analysis of the REMD simulation trajectory was performed. In H2O/DMSO, lorlatinib was found to maintain an extensive hydrogen bond network throughout the simulation in which both lorlatinib and H2O could act as both hydrogen bond donor and acceptor (Figure A). This hydrogen bond network may reduce the flexibility of the compound, explaining the observation of a single conformer in aqueous solution. In contrast, there is no strong interaction between the solvent and lorlatinib in CHCl3 (Figure B), allowing bond rotations.

Figure 3

Intermolecular interaction between lorlatinib and the solvents. (A) There is a hydrogen bond network (red dashed lines) in H2O/DMSO, (B) while no strong intermolecular interaction in CHCl3. (C) Alignment of P1 and P2. The average distance between the HB donor and acceptor during the REMD simulations is provided along the dashed line (Å). The percentage time of HB during the REMD simulations is shown in the brackets.

To further validate the importance of water and the hydrogen bond network in constraining the conformation of lorlatinib, we performed another REMD simulation in pure water. The free-energy profile showed one conformation only (P1, Figure S3), which is the same as the P1 in CHCl3 and H2O/DMSO. Again, a strong HB network was formed between water and lorlatinib where both water and lorlatinib could act as the HB donor and acceptor, similar to what we observed in H2O/DMSO (Figure A). This confirms that the hydrogen bond donor and acceptor capabilities of the solvent may significantly affect the conformational preferences of macrocycles.

We found that the carbonyl group of the amide bond is fully exposed in aqueous solution (1/P1) but partially shielded by the pyridine moiety in conformation 2/P2 in chloroform (Figure C and Figure S4), indicating that 2/P2 might be more hydrophobic. Indeed, the calculated solvent accessible 3D polar surface area[26] of conformation 2/P2 is significantly smaller than that of 1/P1 (193.5:209.8 vs 212.5:216.5 Å2), while the 3D lipophilicity (MLP) of 2/P2 is higher than that of 1/P1 (1.28:1.34 vs 1.19:1.14). In addition, 2/P2 has a lower dipole moment (μ = 3.5/3.6 D) compared to 1/P1 (μ = 4.1/4.4 D) (Figure S4). The more hydrophobic profile of 2/P2 indicates that it should prefer low polarity environments. Therefore, lorlatinib adopts 1/P1 in body fluids, while conformation 2/P2 may help penetrating cell membranes. The low free-energy barrier from P1 to P2 allows the drug to switch its conformation while crossing membranes, thus favoring its cell permeability.

To evaluate the scope of using REMD for conformational sampling of macrocycles, a simulation of the ensemble of the macrocycle radicicol in CHCl3 was carried out. A previous NAMFIS analysis has found that radicicol populates the target (Hsp90) bound conformation to 21% in CDCl3.[21] The REMD simulation found six conformations in CHCl3, with populations ranging from 5 to 33% (Figure S5). The second-ranked conformer (28%) is almost identical to the crystal structure of radicicol complexed with Hsp90 (RMSD 0.18 Å, PDB ID: 1BGQ). Thus, the REMD prediction is in excellent agreement with the reported NAMFIS analysis.[21]

Conclusions

With lorlatinib as an example, we demonstrate that REMD simulations are capable of predicting the conformational ensembles of macrocycles in different solvents. The predicted structures and populations are in excellent agreement with the experimental NMR data. Lorlatinib adopts the same conformer in H2O/DMSO (6:4) as when bound to its protein target, while an additional conformer exists in CHCl3. The hydrogen bond network formed by lorlatinib and water is proposed to play an important role in constraining the macrocycle’s conformation in water. Population of a different and less polar conformation in chloroform, accessed by crossing a low energy barrier, most likely explains how lorlatinib penetrates cell membranes. Additional validation of the REMD simulation approach was obtained for the macrocycle radicicol. This simulation identified a major conformation very similar to radicicol’s target bound conformation, which was previously found by NAMFIS analysis of solution NMR data. Therefore, the REMD simulation shows a potential for prediction of the conformational landscape of macrocycles and for understanding of the membrane permeability of macrocyclic drugs.

Experimental Section

System Preparation and REMD Simulation

The initial structure of lorlatinib was obtained by ChemDraw, followed by an optimization using Gaussian 09 at the B3LYP/6-31G* level. The RESP partial atomic charge was obtained based on optimized structures by the antechamber program. The system was immersed in a cubic solvent box, including CHCl3 or H2O/DMSO (6:4) using the tleap module in the Amber16 package. The total number of atoms in each simulation system was kept at ∼1850. The general Amber force field (GAFF) was applied to obtain the force field parameters of the compound.

The structure of lorlatinib was first minimized with 5000 steps of the steepest descent method, followed by 5000 steps of the conjugated gradient method. To equilibrate the systems, a short MD simulation of 2 ns was performed prior to the REMD simulation. Finally, 100 ns REMD simulations were performed in two different solvents, viz., CHCl3 and H2O/DMSO. Conformations at all temperatures were sampled at every 2 ps. Here, all REMD simulations were performed with 22 replicas at temperatures ranging from 300 to 600 K, which was predicted by the webserver (http://folding.bmc.uu.se/remd/). The conformational change of lorlatinib was monitored through the RMSD of lorlatinib using the cpptraj module in the Amber16 package.

NMR Spectroscopy

NMR spectra recorded at 25 °C on an 800 MHz Bruker Avance III HD NMR spectrometer equipped with a TCI cryogenic probe, for D2O/DMSO-d6 (6:4) and CDCl3 solutions. Assignments were deduced using 1D (1H and 13C) and 2D (COSY, TOCSY, HSQC, HMBC, and NOESY) NMR spectra. NOESY buildups were acquired with seven mixing times of 100, 200, 300, 400, 500, 600, and 700 ms, with 16 scans, 8096 points in the direct whereas 512 points in the indirect dimension and d1 relaxation delay set to 2.5 s, without solvent suppression. Interproton distances (Tables S2 and S3) were calculated according to the initial rate approximation from the linear part of the buildups (r2 > 0.95) using orthoaromatic protons (2.51 Å) and methine to methyl proton distances (2.43 Å) as internal distance references. NOE peak intensities were calculated using normalization of both cross peaks and diagonal peaks according to ([cross peak1 × cross peak2]/[diagonal peak1 × diagonal peak2])0.5. Initial NOE buildup rates were converted into distances using the equation rij = rref(σref/σij) where rij is the distance between protons i and j in angstrom, and σij is the normalized intensity obtained from NOESY experiments. Further information is provided in the Supporting Information.

NAMFIS Analyses

The theoretically available conformational space conformation ensembles of lorlatinib were identified by unrestrained Monte Carlo conformational analysis using five different (OPLS, OPLS-2005, OPLS3e, AMBER*, and MMFF) force fields, each with the GB/SA solvation models chloroform and water.[1] Conformational searches were done using the Monte Carlo algorithm with intermediate torsion sampling, 50,000 Monte Carlo steps, and an RMSD cut-off set to 2.0 Å. Molecular mechanics energy minimization was performed for each conformation using the Polak-Ribière-type conjugate gradient (PRCG), as implemented in the Macromodel BatchMin V12.1 of the Schrödinger Package, with a maximum of 5000 iterative steps. All conformations within 42 kJ/mol from the global minimum were saved (Table S4) and combined with elimination of redundant conformations by comparison of heavy atom coordinates applying an RMSD cutoff set to 1.0 Å, giving the final ensemble used for NAMFIS analysis.

Solution ensembles were determined by fitting the experimentally measured distances to those back-calculated for computationally predicted conformations following previously described protocols.[25,53] CH2 signals were treated according to the equation d = (((d1–6) + (d2–6))/2)−1/6, and methyl signals according to d = (((d1–6) + (d2–6) + (d3–6))/3)−1/6. The NAMFIS ensemble analyses were validated using standard methods, that is, through evaluation of the reliability of the conformational restraints by the addition of 10% random noise to the experimental data, by the random removal of individual restraints, and by comparison of the experimentally observed and back-calculated distances. The result of the NAMFIS analysis is given in Table S7 and Figure S7 for D2O/DMSO-d6 (6:4) and CDCl3. Further information about the NAMFIS analysis is provided in the Supporting Information.

Solvent accessible 3D polar surface area was calculated as previously described using a solvent probe radius of 1.4 Å and with inclusion of partial charges in addition to polar atoms (O, N, and attached H).[26] The molecular lipophilicity potential (MLP)-based virtual log P, a quantitative 3D description of lipophilicity, was calculated as implemented in the VEGA ZZ suite.[54,55]

X-ray Crystal Structures

The published crystal structures of lorlatinib in PDB with different codes have similar structures, and one structure with PDB code: 4CLI was selected and used for this analysis.