Using Membrane Partitioning Simulations To Predict Permeability of

Subscriber access provided by Columbia University Libraries

Computational Chemistry

Using membrane partitioning simulations to predict permeability of forty-nine drug-like molecules Callum J. Dickson, Viktor Hornak, Dallas Bednarczyk, and Jose S. Duca J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00744 • Publication Date (Web): 12 Dec 2018 Downloaded from http://pubs.acs.org on December 12, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Using Membrane Partitioning Simulations to Predict Permeability of Forty-Nine Drug-Like Molecules Callum J. Dickson,1* Viktor Hornak,1 Dallas Bednarczyk,2 Jose S. Duca1

1Computer-Aided

Drug Discovery, Global Discovery Chemistry, Novartis Institutes for BioMedical

Research, 181 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States 2PK

Sciences, Novartis Institutes for BioMedical Research, 250 Massachusetts Avenue, Cambridge,

Massachusetts 02139, United States

ABSTRACT A simple descriptor calculated from molecular dynamics simulations of the membrane partitioning event is found to correlate well with experimental measurements of passive membrane permeation from the high throughput MDCK-LE assay using a dataset of 49 drug-like molecules. This descriptor approximates the energy cost of translocation across the hydrophobic membrane core (flip-flop), which for many molecules limits permeability. Performance is found to be superior in comparison to calculated properties such as clogP, clogD or polar surface area. Furthermore, the atomistic simulations provide a structural understanding of the partitioned drug-membrane complex, facilitating medicinal chemistry optimization of membrane permeability.

KEYWORDS molecular dynamics, partitioning, small molecule membrane permeation, MDCK-LE, bioavailability

ABBREVIATIONS h-bonds: hydrogen bonds; QSAR: quantitative structure-activity relationship; MD: molecular dynamics; MDCK-LE: Madin-Darby Canine Kidney Low Efflux; clogP: calculated logarithm of octanol/water partition coefficient of unionized compound; clogD: calculated logarithm of octanol/water partition coefficient of ionized compound; TPSA: topological polar surface area; MW: molecular weight; HBD: hydrogen bond donor; HBA: hydrogen bond acceptor *Email: [email protected] *Telephone: +1-617-871-4968

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 18

INTRODUCTION Permeability is a desirable characteristic of most small molecule drugs because they frequently need to cross biological barriers to elicit pharmacological activity. Biological barriers may include the intestine for absorption, the blood-brain barrier for central nervous system targets, or cell membranes for intracellular targets. Permeability across an epithelial monolayer is an assay frequently employed to measure the passive transmembrane diffusion of a molecule. In particular, the MDCK-Low Efflux (MDCK-LE) cell line is homogeneous endothilial cell line that grows rapidly, and forms a tight reproducible monolayer in less than four days. The cell line is suitably robust to be amenable for high-throughput screening in 96-well format yielding permeability data familiar to most medicinal chemists. Furthermore, the MDCK-LE cell line is an MDCK cell line subclone that has been selected for low functional expression of P-glycoprotein and shows no functional activity of other common efflux transporters such as Bcrp or Mrp2. The lack of functional xenobiotic efflux transport minimizes underestimates of passive permeation due to efflux and maximizes the passive component of membrane permeation. However, a molecule must first be synthesized to make such passive permeability measurements. By contrast, computational methods offer the ability to influence the design of a molecule prior to synthesis and thus direct molecular design to enrich the synthesis of molecules with higher passive permeability. Computational methods to understand and predict the passive membrane permeability of small molecules typically involve either building knowledge-based models (quantitative structure-activity relationship (QSAR) models), or physics-based approaches (e.g. employing molecular dynamics (MD) simulations). The derivation of QSAR models involves collecting experimental permeation values for a range of molecules, calculating molecular descriptors and then fitting a statistical model capable of predicting permeation values for a new set of compounds, given their calculated descriptors. Provided a high quality training set is available, such a method may be rapidly deployed. However, it suffers from the drawback that predictions on novel molecules far removed from those in the training set may be poor.1 A structurally detailed yet computationally more costly view of membrane permeation may be gained with MD simulations. The potential of mean force (PMF) technique allows calculation of the free energy profile to cross the lipid membrane, from which an estimate of the permeation coefficient may be derived via the inhomogeneous solubility-diffusion model or kinetic framework.2-9 Such an approach often returns the rank order of permeation for drug-like molecules and provides an atomistic view of the permeation event. However, sampling is typically limiting, in particular for beyond rule-of-five molecules,10-12 meaning enhanced sampling methods may be required (such as metadynamics or milestoning).13,

14

Implicit

membrane models are another option to increase sampling of the small molecule in the membrane environment, at the expense of atomistic detail of the membrane.15, 16 Membrane permeation may also be approximated by the energy cost of transfer from a high dielectric (water) to a low dielectric (membrane) environment, a method that has found success for small molecules and cyclic peptides.17-19


Page 3 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this work, we report a simple descriptor (h-bonds while partitioned) that may be determined from MD simulations of membrane partitioned molecules which is found to correlate with MDCK-LE passive permeation values. Compounds are partitioned into a lipid bilayer, upon partitioning the average number of hydrogen bonds the molecule forms with either water or lipid head groups is then determined. The molecular dynamics simulations involved are a fraction of the requirements for a PMF calculation. Previous work found high correlation between h-bonds while partitioned and the free energy cost to flip-flop across the membrane core (ΔGflip-flop), which in turn correlated with membrane permeation, where “flip-flop” is defined as the free energy cost to move from the membrane partitioned position across the membrane core.5 In this work, we further examine the relationship between h-bonds while partitioned and ΔGflip-flop, finding the correlation to hold for a subset of the compounds. Finally, we find that h-bonds while partitioned correlates well with measured permeability of 49 compounds from MDCK-LE experiments, with higher performance than other molecular descriptors such as topological polar surface area (TPSA) or clogP/clogD. Furthermore, average orientations of membrane partitioned molecules provide structural information allowing understanding and medicinal chemistry optimization of passive membrane permeability.



Page 4 of 18

METHODS Dataset The dataset is composed of 49 publicly available weakly basic drug-like molecules with measured MDCK-LE passive permeation values (Papp). The compounds were selected with the criteria of adhering to Lipinski's “rule of five”20 and having a single basic center with calculated pKa in the range 8.5 – 10.5.

Table 1 - The dataset of 49 small drug-like molecules, with measured MDCK-LE permeation Papp and calculated h-bonds while partitioned (average and range over two repeats). Standard deviations in Papp values are provided in cases where multiple assay runs have been performed. Drug Name

log (Papp cm/s)

h-bonds while partitioned (AVG)

h-bonds while partitioned (RANGE)

CHEMBL471

SOTALOL

-6.02

2.72

2.72 - 2.72

CHEMBL128

SUMATRIPTAN

-6.01

2.85

2.80 - 2.89

CHEMBL649

NADOLOL

-6.00 ± 0.103

4.62

4.57 - 4.68

CHEMBL1760

TERBUTALINE

-5.98 ± 0.021

4.56

4.52 - 4.60

CHEMBL26

SULPIRIDE

-5.97 ± 0.006

3.56

3.51 - 3.61

CHEMBL142703

VILDAGLIPTIN

-5.96 ± 0.090

2.43

2.26 - 2.60

CHEMBL714

SALBUTAMOL

-5.92

3.00

2.86 - 3.15

CHEMBL6995

PRACTOLOL

-5.90 ± 0.064

2.75

2.73 - 2.76

CHEMBL24

ATENOLOL

-5.88

3.80

3.58 - 4.03

CHEMBL776

METAPROTERENOL

-5.62

3.13

3.12 - 3.14

CHEMBL483254

PANOBINOSTAT

-5.51 ± 0.037

4.03

3.94 - 4.12

CHEMBL28992

-

-5.44

1.52

1.45 - 1.59

CHEMBL178291

-

-5.38

1.36

1.13 - 1.58

CHEMBL571948

Y-39983

-5.37

2.98

2.96 - 3.00

CHEMBL27846

-

-5.33

1.42

1.34 - 1.50

CHEMBL18274

-

-5.32

1.55

1.31 - 1.79

CHEMBL263664

CCT128930

-5.29

2.53

2.40 - 2.66

CHEMBL18041

ZACOPRIDE

-5.23

1.54

1.46 - 1.61

CHEMBL563327

-

-5.11

3.43

3.42 - 3.43

CHEMBL494089

GSK-690693

-5.08

2.53

2.37 - 2.69

CHEMBL371064

BALICATIB

-5.03

2.36

2.01 - 2.71

CHEMBL21731

MAPROTILINE

-5.03

0.78

0.75 - 0.81

CHEMBL2326015

-

-4.99

2.27

2.21 - 2.34

CHEMBL500

PINDOLOL

-4.92 ± 0.068

2.04

1.95 - 2.13

CHEMBL71

CHLORPROMAZINE

-4.90

0.29

0.27 - 0.32

CHEMBL445

NORTRIPTYLINE

-4.89

0.87

0.79 - 0.96

CHEMBL72

DESIPRAMINE

-4.85

1.03

0.88 - 1.17

CHEMBL6966

VERAPAMIL

-4.84 ± 0.102

0.78

0.72 - 0.84

CHEMBL350615

N-BENZYLBENZAMIDINE

-4.83

1.40

1.36 - 1.45

CHEMBL1764

LEVOMEPROMAZINE

-4.79

0.51

0.50 - 0.52

CHEMBL669

CYCLOBENZAPRINE

-4.79

0.45

0.37 - 0.54

CHEMBL12610

BENZYDAMINE

-4.78 ± 0.008

0.92

0.91 - 0.94

CHEMBL2106741

MEZACOPRIDE

-4.78

0.59

0.55 - 0.63

CHEMBL ID


Page 5 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CHEMBL643

PROMETHAZINE

-4.78 ± 0.238

0.26

0.22 - 0.30

CHEMBL171070

-

-4.76 ± 0.111

0.40

0.36 - 0.44

CHEMBL814

FLUVOXAMINE

-4.74 ± 0.067

1.73

1.72 - 1.73

CHEMBL564

PROMAZINE

-4.74

0.50

0.45 - 0.56

CHEMBL549

CITALOPRAM

-4.73 ± 0.088

0.63

0.63 - 0.64

CHEMBL27

PROPRANOLOL

-4.66 ± 0.108

1.69

1.63 - 1.74

CHEMBL95

TACRINE

-4.64

0.83

0.77 - 0.88

CHEMBL266195

ALPRENOLOL

-4.64

1.41

1.34 - 1.47

CHEMBL134

CLONIDINE

-4.63

1.42

1.33 - 1.51

CHEMBL13

METOPROLOL

-4.61 ± 0.106

2.51

2.40 - 2.63

CHEMBL546

OXPRENOLOL

-4.60

1.13

1.04 - 1.21

CHEMBL652

FLECAINIDE

-4.58 ± 0.124

0.44

0.33 - 0.55

CHEMBL505

CHLORPHENAMINE

-4.56

0.69

0.65 - 0.73

CHEMBL1483568

-

-4.56

0.13

0.08 - 0.17

CHEMBL657

DIPHENHYDRAMINE

-4.51

0.56

0.50 - 0.61

CHEMBL11

IMIPRAMINE

-4.42 ± 0.158

0.36

0.35 - 0.38

Experimental measurement of permeation with MDCK-LE cell monolayer assay Cell culture MDCK-LE cells are a low efflux (LE) cell line developed internally at Novartis using methodology similar to that previously reported.21 MDCK-LE cells were cultured at 37°C under a 5% CO2 atmosphere, at 95% relative humidity in DMEM containing 10% FBS, penicillin-streptomycin (100 μg/mL), and 2 mM Ala-Gln. Cells were passaged weekly into an Omnitray (Nunc, Thermo Fisher Scientific, Rochester, NY) at a density of approximately 5750 cell/cm2 for continuous culture. For assay purposes cells were seeded at a density of approximately 265,000 cells/cm2 of a 96-well Transwell plate (Corning Life Sciences, Acton, MA) and cultured in the same media noted above for a period of four days. Assay The determination of the apparent permeability (Papp) was performed in the A->B direction where each compound was assayed in triplicate in a pool of three compounds, similar to that previously reported.21 The zwitterion bestatin, a poorly permeably compound, was used as marker of monolayer integrity and was included as a fourth compound in the cassette of three compounds. To initiate the assay, media was aspirated and the cells and basal chambers were washed three times with Hank’s Balanced Salt Solution (HBSS) containing 10 mM HEPES (pH 7.4). Compound test solutions were prepared in triplicate in HBSS containing 10 mM HEPES (pH 7.4) and 0.02% bovine serum albumin (BSA) to a final concentration of 10 µM and centrifuged for 2 minutes at 4000 g, then applied to the apical compartment (donor compartment) at time zero. Additionally at time zero, a 37°C solution without test articles (HBSS+10 mM HEPES (pH 7.4) plus 0.02% BSA) was added to the receiver chamber (basal) of the Transwell plate. A time zero sample of the donor solution was also sampled for further analysis. The assay was conducted for a period of 120 minutes at 37°C without shaking. At the time of assay termination samples were taken



from each donor compartment and each acceptor compartment of the Transwell plate. To each of the 0 and 120 minute sample was added an internal standard solution containing glyburide in water:acetonitrile, 50:50 (v:v)). Concentration curves were prepared using a Labcyte Echo in the same matrix noted above. Sample and concentration curve sample were analyzed by mass spectroscopy. Mass Spectroscopy Assay samples were loaded onto a RapidFire C4 cartridge by means of a RapidFire autosampler (Agilent, Santa Clara, CA). Chromatography was performed at a flow rate of 1.25 mL/min, loading with 0.1% formic acid in water and eluting in 0.1% formic acid in methanol. Mass spectroscopy was performed using an AB Sciex API5500 (Sciex, Frammingham, MA) equipped with a turbo ion spray source. The analyte concentration was calculated from the chromatographic peak area ratio of analyte to internal standard (glibenclamide, m/z 494→169), using Multiquant software V2.1(Sciex, Framingham, MA). Calculations Papp values were determined as: 𝑷𝒂𝒑𝒑 =

(

𝑽𝑨 𝑨 ∗ [ 𝑫 𝟎]

[𝑨𝟏𝟐𝟎]

)( ) ∗

𝒕

Percent recovery values were determined as: %𝒓𝒆𝒄𝒐𝒗𝒆𝒓𝒚 = 𝟏𝟎𝟎 ∗

[𝑨𝟏𝟐𝟎] + [𝑫𝟏𝟐𝟎] [ 𝑫 𝟎]

(

)

Where VA is the volume of the acceptor (mL), A is the surface area of the membrane, D0 is the donor solution concentration at t=0, D120 is the donor solution concentration at t=120, A120 is the acceptor solution concentration at t=120, and t = time (seconds). The minimum significant ratio22 averaged 1.75 from independent experiments of 45 compounds used to validate the assay and the mean relative standard deviation was 19.8%. Papp values were determined from triplicate wells using the mean acceptor and donor solution concentrations. Where independent experiments were performed standard deviations are provided to distinguish from replicate samples within a single experiment.23

Calculated physicochemical properties Calculated logP and topological polar surface area values were determined using the RDKit chemoinformatics software from the canonical isomeric SMILES molecular representation of each molecule.24 Calculated logD values were determined with the MoKa software from the canonical isomeric SMILES molecular representation of each molecule.25 A 3D polar surface area was determined with COSMOtherm26, 27 by summing the charge/area values for a low energy 3D conformer (as determined in dielectric corresponding to water) of each molecule.


Page 6 of 18

Page 7 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Membrane partitioning simulations and average hydrogen bond count For each molecule in the dataset (see Table 1) a partitioning simulation was run, whereby the drug molecule was inserted into the water layer above a 72 POPC bilayer with 60 waters/lipid and simulated for 100 ns, during which time all molecules moved into the lipid membrane, positioning just below the lipid head groups. Such a position is consistent with the free energy minima for a number of other small molecules.28 In all cases, the drug was charge neutral. Drugs were modelled with the parm@Frosst force field,29 a small molecule force field that extends AMBER ff99SB,30 and conformationally averaged AM1BCC charges,31-34 lipids were modeled using the AMBER Lipid14 force field35-37 and water using the TIP3P model.38 Molecular dynamics simulations were run with AMBER 1639 and PMEMD CUDA40-42 on GPU cards. The system was first energy minimized for 10000 steps, of which the initial 5000 steps used the steepest descent method and the remaining steps used the conjugate gradient method.43 Heating from 0 K to 100 K was then applied using Langevin dynamics44 within a 5 ps constant volume run, with restraints on the drug molecule and lipids (force constant 10 kcal/mol/Å2). The volume was then allowed to change freely and the temperature increased to 303 K (the temperature at which bulk phase properties of POPC have been validated with Lipid14) with a Langevin collision frequency of γ = 1 ps-1, and anisotropic Berendsen control of the pressure45 around 1 atm was applied by coupling the periodic box with a time constant of 2 ps for 100 ps. The same restraint of 10 kcal/mol/Å2 was maintained on the drug and lipid molecules. The pressure relaxation time was then reduced to 1 ps, all restraints were removed and the system was simulated for 100 ns in NPT. Three dimensional periodic boundary conditions with the usual minimum image convention were employed. Bonds involving hydrogen were constrained using the SHAKE algorithm,46 allowing a 2 fs time-step. PME was used to treat all electrostatic interactions47 beyond a cutoff of 10 Å. A long-range analytical dispersion correction was applied to energy and pressure. Once all simulations were complete, the position of the center-of-mass of the drug along the z-axis was analyzed. In all cases, the drug moved to a position of z ≈ 10 Å from the middle of the bilayer, demonstrating that it partitioned into the membrane, arriving at a position just below the lipid head group region. A snapshot with the drug at a membrane partitioned position was then saved. This was extended for a further 100 ns with the drug restrained at 10 Å from the bilayer center-of-mass to obtain sampling of the drug while membrane partitioned. The drug’s center-of-mass was restrained in the z-dimension only (free to diffuse in xy-plane) using a harmonic force constant of 2.5 kcal/mol/Å2. The initial 20 ns were taken as equilibration; production simulations were therefore 80 ns. This simulation was repeated with different random starting velocities to obtain an error estimate (final results are therefore the average of the two independent 80 ns drug restrained runs and range over the two repeats).



For a subset of the 49 compounds with diverse Papp values, the full free energy profile to traverse the POPC bilayer was determined using the potential of mean force method (the following eight molecules were selected: verapamil, propranolol, chlorpromazine, desipramine, atenolol, metoprolol, nadolol, terbutaline). The simulation protocol was identical to that described in reference 5. The free energy of flip-flop (∆Gflip-flop) may then be determined as the free energy cost to move from the free energy minima (in this case, the membrane partitioned position, approximately 10 Å from the bilayer center) to the membrane core. The total simulation time to compute the full PMF profile for a single compound was approximately 2000 hours using one GPU; conversely, the membrane partitioning plus two repeat membrane partitioned production runs for a single compound took 140 GPU hours. The average number of hydrogen bonds each drug molecule makes to either lipid head groups or water molecules while partitioned was determined by counting hydrogen bonds using CPPTRAJ48. The hydrogen bond definition was distance 3.5 Å and angle greater than 135°.

RESULTS h-bonds while partitioned approximates the free energy cost to membrane flip-flop

Figure 1 - Free energy cost to membrane flip-flop calculated with PMF calculations versus h-bonds while partitioned for eight small molecules. The correlation coefficient is r2=0.85. Errors in ∆Gflip-flop values from Monte Carlo bootstrap trials are given in Table S1 and are smaller than the plot markers.


Page 8 of 18

Page 9 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For a subset of eight molecules from the dataset, the full free energy of transfer profile across a POPC membrane was determined (see SI Table S1 and Figure S1). The selected compounds have a range of Papp permeation values (see Table 1). From Figure 1 it can be seen that h-bonds while partitioned has a strong correlation with the ΔGflip-flop barrier from the full PMF calculation (r2=0.85); this barrier will therefore be approximated using a membrane partitioning simulation for the molecules in the dataset.

h-bonds while partitioned correlates with permeation

Figure 2 - Plot of h-bonds while partitioned versus permeability. Values for h-bonds while partitioned are the average and range over two independent runs. The correlation coefficient is r2 = 0.63. The plot of Papp versus h-bonds while partitioned is shown in Figure 2; the correlation coefficient is found to be r2 = 0.63, indicating that this descriptor may be suitable as a virtual screening tool to optimize passive permeation. Di et al recommend Papp >2.5 x 10-6 cm/s (logPapp> -5.6) for medium to high permeability and optimal %Human Intestinal Absorption,21 from Figure 2 an approximate cut-off value of

Using Membrane Partitioning Simulations To Predict Permeability of

Recommend Documents