Permeability Coefficients of Lipophilic Compounds Estimated by

Subscriber access provided by La Trobe University Library

Article

Permeability coefficients of lipophilic compounds estimated by computer simulations Zhaleh Ghaemi, Domenico Alberga, Paolo Carloni, Alessandro Laio, and Gianluca Lattanzi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01126 • Publication Date (Web): 08 Jul 2016 Downloaded from http://pubs.acs.org on July 9, 2016

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 20

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 Theory and Computation

Permeability coefficients of lipophilic compounds estimated by computer simulations Zhaleh Ghaemi,† Domenico Alberga,‡ Paolo Carloni,∗,¶ Alessandro Laio,† and Gianluca Lattanzi§ †SISSA, Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 265, 34136 Trieste, Italy ‡Dipartimento Interateneo di Fisica "M. Merlin", University of Bari "Aldo Moro", TIRES & INFN, via Orabona 4, 70126 Bari, Italy ¶Computational Biophysics, German Research School for Simulation Sciences, D-52425 Julich, Germany and Institute for Advanced Simulation, Forschungszentrum Julich, D-52425 Julich, Germany §Dipartimento di Medicina Clinica e Sperimentale and INFN - Sez. di Bari, Viale Pinto, 71122 Foggia, Italy E-mail: [email protected]

Abstract The ability of a drug to cross the intestine-blood barrier is a key quantity for drug design and employment, and is normally quantified by the permeability coefficient P , often evaluated in the so called Caco-2 assay. This assay is based on measuring the initial growth rate of the concentration of the drug beyond the cellular barrier, but not its steady state flux through the membrane. This might lead to confusion since, in the

1

ACS Paragon Plus Environment


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 20

case of lipophilic drugs, the initial slope is strongly affected by the retention of the drug in the membrane. This effect is well known, but seldom considered in the assay. Here, we exploit all-atoms molecular dynamics and bias exchange metadynamics to calculate the concentration of two lipophilic drugs across a model membrane as a function of time. This allows estimating both the steady state flux and the initial slope of the concentration growth, and comparing Caco-2 and steady state estimates of P . We show that our computational procedure is able to reproduce the experimental values, although these may differ from the permeability coefficients by orders of magnitude. Our findings are generalized by a simplified one-dimensional model of the permeation process, that may act as a roadmap to assess which measure of membrane permeability would be more appropriate and, consequently, whether retention corrections should be included in estimates based on Caco-2 assays.

Introduction A significant fraction of drug candidates and their metabolites must cross the cell membrane to reach their targets. 1 In particular, drugs that are taken orally (the most preferred way of drug administration), must pass several membranes before being absorbed in the blood. Membrane crossing may be assisted by carrier-mediated transport 2 or occur passively in a process fully dominated by drug diffusion. 3 In the latter case, determining the permeability coefficient (P ) of a drug, which describes its permeation speed through the membrane, is of paramount importance for its pharmacological potential. P is defined as the ratio between the flux through the membrane (J) and the drug concentration difference (∆C) across it:

P =

J ∆C

(1)

The most widespread method to measure P for pharmaceutical purposes is the so-called carcinoma colorectal cell based (Caco-2) assay. 3 Within this method, a drug is prepared in a “donor” compartment at an initial concentration C0 and allowed to permeate through a 2


Page 3 of 20

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


monolayer of human colon epithelial cancer cells towards a “receiver” compartment where the concentration of the drug CR (t) is measured as a function of time. At large drug concentrations, i.e. C0 >> CR (t), the drug can be assumed to flow only in one direction. P is then estimated as: 4 VR dCR (t) P = dt t=0 AC0

(2)

where VR is the volume of the receiver compartment and A is the area of the layer separating the two compartments. 5 In what follows we denote Jini =

dCR (t) . dt

If the intestine-blood barrier

behaves as a barrier, the flux at the onset of the permeation Jini and the steady-state flux J entering Eq. 1 are very similar, and Eq. 2 provides a correct estimate of the permeability coefficient. This is actually always the case when the drug is hydrophilic, or the transport predominantly takes place via active carriers. However, for a lipophilic drug that permeates passively, the flux Jini can differ from the steady-state flux J by several orders of magnitude. This is because the blood-intestine “barrier” for a lipophilic compound is actually a sink. The flux at the onset of the permeation process, Jini , is limited by the time required to escape from this sink. At the steady state, the membrane is saturated by the drug, and the flux becomes much larger, since it is limited only by diffusion and not by the escape process anymore. This effect is actually well known, and it has been previously reported in literature. 6 Possible solutions included the development of less biased assays 7 or introducing empirically the so-called retention correction. 8 In the latter case, the estimate of the permeability coefficient is modified as P =

dCR (t) VR 1 , dt AC0 1−R

where R is the fraction of compound retained in

the membrane. 8 This fraction can be very close to one for strongly hydrophobic compounds, significantly boosting the value of P with respect to the value predicted by Eq. 2. Also molecular dynamics calculations 9 have shown that the permeation properties differ considerably for different penetrants in particular when the membrane acts as a trap rather than a barrier. 10,11 In a previous publication, we derived an expression for estimating P 3



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 molecules with several degrees of freedom using simulations that could reproduce with a reliable accuracy the permeability coefficient of ethanol. 12 However, a computational protocol that allows estimating P for compounds with different chemical nature is still missing. Here we demonstrate that Caco-2 essays of lipophilic compounds do not actually measure the permeability coefficient, but rather a quantity that may significantly differ from it. To this aim, we employ a general computational protocol that simulates the permeation process of drugs. The protocol allows to estimate the membrane permeability and to calculate the time-dependence of the concentration of a drug across a model membrane barrier, thus simulating an idealized Caco-2 experiment. We consider the permeation process of a few compounds through a model palmitoyloleoylglycerophosphocholine (POPC) membrane. We perform molecular dynamics simulations, exploiting bias exchange metadynamics, 13 an approach that is able to reconstruct the multidimensional free energy landscape for complex biochemical processes. We describe a method to derive a master equation from the analysis of the detailed molecular simulations. The obtained master equation allows modeling the permeation kinetics on the time scale of seconds and for a large number of permeating molecules. In this manner, we are able to simulate a Caco-2 experiment, and compute the value of P by both Eqs. 1 and 2. Within this approach, we verified that the P values estimated from the two equations are consistent for a hydrophilic molecule, ethanol, but they differ by 4 and 8 orders of magnitudes for two FDA approved anti-HIV-1 drugs, namely efavirenz (EFV) and etravirine (ETR), respectively. Both molecules were referred to as lipophilic in the pharmaceutical literature. 14,15 It may be argued that, even in this case, the initial uptake is still the relevant biochemical quantity, rather than the permeation coefficient in the steady state flux. However, the discrepancy between these two quantities clearly points out that the two methods are not equivalent and suggests that lipophilic drugs might have been discarded just because of their initial retention in the membrane. The Caco-2 assay employs living cells, with membranes containing a large variety of lipids and active transporters. Clearly, any model (including the one presented here) that

4


Page 4 of 20

Page 5 of 20

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


can be studied by atomistic simulations with the computational resources available nowadays cannot take this complexity into account. However, Eq. 2 is based on the slope of the flux at the onset of the permeation. This, in turn, is determined, both in the Caco-2 assay and in our simulations, by the escape rate of the molecules trapped in the membrane that is in direct contact with the water solution in the receiver compartment. In this respect, the number of membranes separating the two compartments is not relevant, since only the last membrane will release the drug whose concentration is measured by the detector. Also the exact microscopic composition of the membrane is, at least in a first approximation, irrelevant, since the escape rate is primarily determined by the depth of the sink associated with the membrane, that is approximately similar for different fatty acids, since it is primarily determined by the lipophilicity. Hence, we expect agreement between our computational estimates of permeability based on Eq. 2 and the experimental values obtained in Caco-2 assays. Building on these results, we introduce a simple mesoscopic model of the permeation process aimed at predicting when and why discrepancies between the estimates of P provided by Eqs. 1 and 2 can be expected. In short, this model shows that while the estimate of Eq. 1 is determined by the free energy of the highest barrier, the estimate of Eq. 2 is determined by the free energy difference between the highest barrier and the deepest free energy well. This explains why the differences between the two predictions is particularly significant for lipophilic drugs.

Methodology We performed molecular dynamics simulations on the permeation process of one EFV and one ETR molecule through a 2×16 POPC lipid bilayer solvated with 1448 (1441) TIP3P water molecules, 16 respectively. POPC lipids were parametrized as in Ref. 17 and a standard procedure was used for the parametrization of the drugs: molecule charges were computed

5



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

via the RESP method, 18 by fitting the electrostatic potential calculated from Density Functional Theory with 6-31G∗ basis set and a hybrid exchange-correlation functional (B3LYP). Geometry-related parameters were the same as in GAFF. 19 The charges for the two drug molecules are reported in Supplementary Information. The permeation processes of these lipophilic compounds are rare events on the molecular time scale: a similar computational procedure was already employed to describe accurately the permeation of ethanol and is fully reported in ref. 12 The permeation process is satisfactorily described at the mesoscopic level in terms of the collective variable (CV) z, defined as the z-component of the vector pointing from the center of mass of the membrane to the center of mass of the permeating drug (the z-axis of the simulation cell being orthogonal to the membrane surface). However, in ref., 12 we showed that neglecting the multidimensionality inherently associated with the molecular nature of the permeation process leads to a slowly decaying memory kernel in the dynamics with respect to the z variable. We therefore build a kinetic model in a larger CV, in which the dynamics is Markovian. Briefly, we first reconstructed the multi-dimensional free energy surface associated with the permeation process by using bias-exchange metadynamics (Be-Meta) 13 with 13 and 14 collective variables (CVs), for EFV and ETR, respectively. The first CV is z while 9 (12) other CVs for ETR (EFV) measure the number of contacts between the permeating molecule and the water molecules, the phosphate groups, and the aliphatic carbons of the lipid tails; finally, other 4 CVs describe the conformation of the internal dihedrals of the ETR (see Supplementary Information for details). The values of the Gaussian widths were chosen upon the procedure described in ref., 20 namely they were taken as the standard deviations of the corresponding CVs obtained in a preliminary unbiased simulation. We performed a total of 3.2 and 2.6 µs of Be-Meta simulations for ETR and EFV, respectively. All simulations were performed using the GROMACS 4.0 package. 21 , 22 using a NoséHoover thermostat 23,24 with a coupling constant of 1 ps. All bonds were constrained to their equilibrium values by the LINCS algorithm. 25 Lennard-Jones interactions were cut off at a

6


Page 6 of 20

Page 7 of 20

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


distance of 1.4 nm and the time step was set to 2 fs. Long-range electrostatic interactions were computed with the particle-mesh Ewald algorithm: 26 the numbers of grid points along x, y and z were 30, 26 and 66, respectively, the interpolation order was set to 4, the relative accuracy of direct/reciprocal space was 10−5 and the distance for the Coulomb cut-off was set at 1.4 nm. The systems were coupled to a heat bath at a temperature of 323 K, well above the POPC phase transition temperature, 22 using a Nosé-Hoover thermostat 23,24 with a coupling constant of 14 ps. First, the systems underwent molecular dynamics simulations for 5 ns in the NPT ensemble. The pressure was controlled in the direction normal to the bilayer surface by a Parrinello-Rahman barostat 27 with a coupling time constant of 1 ps. Then, the systems underwent Be-Meta in the NVE ensemble. The obtained trajectories were processed to reconstruct a thermodynamic and kinetic model within a suitably chosen subset of CVs (see Supplementary Information), as described in Ref. 28 To this aim, we exploited the low-dimensional free energies obtained from Be-Meta to estimate by a weighted histogram procedure the free energies Fα of a finite number of microstates, labeled by α, that are representative of all the configurations explored by the system. Then, we determined the position-dependent diffusion matrices Dα associated to the transition processes between these microstates. We here assume that the diffusion matrix takes three different values within the water region (WAT), the region under the surface of each membrane leaflet (HB), and the hydrophobic core of the membrane (MID), and is constant within these regions. Moreover, we assume that D is diagonal. Both these assumptions can induce small artifacts in the dynamics, and could be relaxed, but only at the price of optimizing the likelihood by much longer unbiased molecular dynamics simulations. The values of the diffusion matrices are reported in Supplementary Information. The dynamics of the system is then described by a master equation for the probabilities of the microstates (pα ): dpα X = (kβα pβ − kαβ pα ) , dt β 7


(3)


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 8 of 20

in which the transition rates are:

kαβ where Dαβ =

1 2

1 Dαβ exp − = (Fβ − Fα ) 2kB T |sα − sβ |2

(4)

(Dα + Dβ ), sα is the position of the center of microstate α in the space

of the CVs, T is the temperature and kB is Boltzmann’s constant. According to Eq. 4 the prefactors in the kαβ transition rates depend on the inverse of the square of the size of the microstates. However, this is essential to ensure that the observable (the permeability coefficient in our case) is, at least approximately, independent on the size itself. We have checked this hypothesis by increasing up to 20% the dimensions of the microstates in each direction of the manifold defined by the chosen CVs: in all the examined cases, the permeability coefficients calculated with both equations did not change significantly. The estimate of P corresponding to Eq. 1 was obtained by the numerical solution of Eq. 3 subject to the conditions that pα = C, where C is a fixed concentration, for all microstates α for which zα < −zm and pα = 0 for zα > zm , and zm = 2.1 nm is the value of z at which the free energy becomes constant in z within the error bars1 . Imposing these boundary conditions induces a stationary current J from the source (microstates with zα < −zm ), P towards the sink (microstates with zα > zm ). J is estimated as Jα = β kβα pβ for α such that zα = zm and the concentration difference in Eq. 1 is given by ∆C = C.

The estimate of P corresponding to Eq. 2 was obtained by numerical solution of Eq. 3 with initial condition of pα = CS if zα < −zm and zero otherwise. We then computed the P concentration in the receiver compartment CR (t) = α:zα >zm pα (t). The values of CR as a function of time for EFV, ETR and ethanol are reported in Supplementary Information. dCR (t) dt

in Eq. 2 corresponds to the early slope of the concentration relaxation in the receiver

compartment, calculated with a linear fit by taking the maximum number of points for which r 2 = 0.99. 1 Different choices for the membrane boundaries lead to a slight decrease of P when the membrane thickness is increased to 5 nm. P becomes practically independent from the membrane thickness above 5 nm (see Supplementary Information).

8


Page 9 of 20

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


The simulations and the free energy calculations are necessarily performed on a small simulation box, in which the donor and the receiver compartments are 1-2 nm thick. However, in the rate model one can easily increase the volume of the donor (or the receiver) compartment by adding as many microstates as necessary along the z direction, with a constant free energy value. This allows analyzing the influence of the volume of the donor compartments on the estimated value of P . The result of this analysis is reported in the Supplementary Information.

Results and discussion In order to estimate the permeation coefficients obtained by Eqs. 1 and 2, we first considered a minimalist model of the permeation process, in which the permeating molecule has to cross a bilayer of POPC molecules. We here considered two lipophilic drugs, Efavirenz (EFV) and Etravirine (ETR), while a hydrophilic molecule, ethanol, was already studied with the same setup in Ref. 12 We model the permeation process by accurate atomistic simulations, using the approach described in detail in Ref. 12 In short, we first compute the multi-dimensional free energy surface as a function of a large set of collective variables associated with the permeation process by using bias-exchange metadynamics (Be-Meta). 13 Afterwards, we identify a set of microstates on this free energy surface. A microstate includes configurations with similar values of all collective variables. Finally, we estimate the transition probabilities between all microstates. This approach, described in details in the Methodology section and SI, allows predicting the long-time evolution of any explicit function of the collective variables by solving a master equation in the probability distribution of the microstates. Hence, following this methodology, we are able to estimate the current through the membrane at the beginning of the permeation process (Jini ) entering in Eq. 2 and the steady-state flux J entering in Eq. 1, and compute the corresponding values of the permeability coefficient.

9



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

Permeation free energy surfaces Fig. 1 shows the permeation free energies for ETR, EFV and ethanol, projected along the z component of the vector pointing from the center of mass of the membrane to the center of mass of the permeating drug (z). The permeation of a hydrophilic molecule such as ethanol is associated with a barrier while for the two lipophilic molecules, ETR and EFV, the free energies display minima inside the membrane.

Figure 1: Free energy profiles as a function of the z component of the vector pointing from the center of mass of the membrane to the center of mass of the permeating drug: for lipophilic molecules, ETR and EFV, the profiles include significantly deep wells while for the hydrophilic molecule ethanol, the profile is characterized by the presence of a high barrier. Free energy profiles are obtained by averaging the one-dimensional history dependent potentials of Be-Meta simulations. The calculations for ethanol are reported in Ref. 12

10


Page 10 of 20

Page 11 of 20

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


Representative states: ETR The red profile in Fig. 1 shows the free energy profile of ETR permeation as an average of history dependent potentials of Be-Meta simulations along the collective variable z. The free energy minimum with respect to this coordinate (HBET R ) corresponds to the representative structure shown in Fig. 2A. In this state, the nitrogen atoms are involved in hydrogen bonds as well as in electrostatic interactions with the carbonyl atoms of POPC molecules. Although this state corresponds to a global free energy minimum, it is clearly rate-limiting at the onset of the permeation process, since these interactions would keep ETR in the close proximity of the membrane surface. The middle of the membrane represents a local maximum for the ETR free energy profile. In this state (MIDET R ), the aliphatic chains of POPC molecules produce steric hindrances that force ETR to assume two possible orientations, well described by the values of the dihedrals 2 and 3 as shown in SI Fig. 2. At the center of the membrane, the extended (Fig. 2B) and compact (Fig. 2C) conformations display a probability of 66% and 34%, respectively. This is consistent with the low mass density of the membrane in the center (SI Fig. 4). However, going away from the center, as the density of the POPC membrane increases (SI Fig. 4) and because of the presence of the elongated aliphatic chains, the trend reverses and the ETR compact conformation becomes more probable (70% of analysis time). ETR conformational changes are therefore the rate-limiting steps for the molecule translocation across the membrane. Representative states: EFV The same analysis performed on the EFV free energy profile (blue line in Fig. 1) points out similar results. Representative structures for the HBEF V state are shown in Fig. 2D-E, corresponding to the presence of hydrogen bonds with the phosphate group (Fig. 2D) or the carbonyl (Fig. 2E) of the phospholipid molecule. Again, these hydrogen bonds would tend to favor the interaction of EFV with the hydrophilic heads of the POPC molecules, with the 11



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

Figure 2: Representative snapshots of ETR and EFV permeation. (A) HBET R state: the interaction between the nitrogen atoms of ETR and the polar heads of phospholipid molecules results in a free energy minimum. (B) In the compact conformation of ETR, both dihedrals 2 and 3 are in the trans conformation. (C) In the extended conformation of ETR, dihedral 2 is in the cis conformation, while dihedral 3 is in trans. In the middle of the membrane, there is a net preference for the extended conformation, while the compact conformation is favored away from the centre, where the mass density of the membrane increases. (D-E) HBEF V state: the nitrogen atoms of EFV form hydrogen bonds with (D) the phosphate group and (E) the carbonyl of the phospholipid molecules. These interactions result in a free energy minimum. (F) Rate limiting step in the EFV translocation across the model POPC membrane corresponds to an effective chain of water molecules that protrudes in the membrane. In all panels, POPC molecules coordinating the drugs are represented in licorice, the drugs EFV in CPK, other POPC molecules as gray lines.

12


Page 12 of 20

Page 13 of 20

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


whole molecule inside the membrane but in the close proximity of its surface. The MIDEF V state corresponds to the loss of the hydrogen bonds formed with the POPC molecules. However EFV is able to bind water molecules that can be dragged to the middle of the membrane, as clearly shown in the representative snapshot in Fig. 2F. This behavior is quite similar to other previously reported computational studies on small molecule permeation through model POPC membranes. 29 The insertion of water molecules inside the membrane is usually accompanied by a disruption of the membrane surface that results in a higher free-energy barrier. This is also visible in the blue profile of Fig. 1.

Estimate of the permeability coefficients We estimated P using two approaches: one from Eq. 1, which relies on the definition of P , and the other one from Eq. 2, corresponding to typical Caco-2 experiments. Fig. 3 reports P values for the molecules versus their free energy differences between the solution and the center of the membrane. Hence, the free energy difference will be positive for hydrophilic molecules like ethanol. On the contrary, lipophilic molecules will be characterized by negative free energy differences, corresponding to their effective preference for the membrane. As shown in Fig. 3, the P values estimated by atomistic simulations applying Eqs. 1 and 2 are consistent only for ethanol, while they substantially differ for both lipophilic drugs. Importantly, the P values obtained from our simulations are close to those reported in Caco2 assays for EFV 15 and ETR 30 only when the estimates are based on Eq. 2, indicating that for these two drugs our computational approach is able to predict the experimental outcome. This implies that the experimental estimates of P reported in Refs. 15 and 30 do not correspond to the permeability coefficient as defined by Eq. 1.

One-dimensional model of membrane translocation Our analysis shows that it is possible to simulate the time dependence of the concentration of a permeant beyond a model membrane, and compute the permeability coefficient by both 13



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

Figure 3: Estimates of P based on computed free energy differences for the permeation of EFV, ETR and ethanol through a model POPC membrane. The computed ∆F values are evaluated by taking into account the full reconstructed multidimensional free-energy profiles and hence may slightly differ from the heights of the barriers reported in Fig. 1. Experimental data from Caco-2 assays are reported for comparison, when available. Eqs. 1 and 2. Which equation is suitable for the estimate does not matter in the case of hydrophilic drugs, but is extremely relevant in the case of lipophilic drugs. It is therefore important to establish which quantity is relevant in each specific case, the steady state flux or the flux at the onset of the permeation process. To further quantify the difference between these two quantities, we introduce here a simplified model system in which a molecule diffuses on a one-dimensional free energy landscape at a temperature T = 1. A membrane, or a sequence of membranes, in this model is simply represented by a series of barriers and wells on which the drug will diffuse. Some examples of the profiles we considered are shown in panels A-D in Fig. 4. For each of the 4 shapes, we considered several choices of barrier height ∆Gb and maximum well depth ∆Gw . The chosen values for ∆Gb and ∆Gw are listed in the panels. The value of P for these systems is calculated by solving numerically the associated onedimensional Fokker-Plank equations with the appropriate boundary conditions, as reported in the Methodology section. The results are presented in the top and bottom panels on the left. Consistently with the results obtained by atomistic simulations, the values of P 14


Page 14 of 20

Page 15 of 20

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


calculated with Eqs. 1 and 2 coincide only for “hydrophilic” molecules, for which the profile is characterized by a significant barrier and a small well. Instead, the two estimates differ substantially in the presence of a significant free-energy minimum. In the upper panel, the values of P predicted by the two equations are presented as a function of ∆Gb . The full points, that correspond to the estimate provided by Eq. 1, are approximately fit by a line of equation log P = A −

∆Gb , T

where A is a constant. This indicates that, regardless from the

depth of the well, the value of P provided by Eq. 1, is determined solely by the exponential of the barrier height. The empty points, corresponding to the estimate provided by Eq. 2 are always below the full points, indicating that the estimate of P provided by this equation is systematically lower. In the upper panel, the values of P are presented as a function of ∆Gb + ∆Gw . In that case, the empty points are approximately fit by a line of equation w , while in this representation the full points cannot be fitted. This result log P = A− ∆Gb +∆G T

is in agreement with Ref., 10 where it is also observed that the maximum well depth must be taken into account in order to predict the value of the permeability coefficient reported in Caco-2 assays. This indicates that the value of P provided by Eq. 2 is approximately determined by the exponential of the sum of the barrier height and of the well depth. This analysis explains why the two estimates are different, and allows relating the estimates to microscopic quantities characterizing the membrane.

Conclusions In this work we introduced a computational procedure that allows simulating the time dependence of the concentration of a permeant (for example a drug) beyond a model membrane barrier. The procedure allows, within the same computational framework, estimating the permeability coefficient and simulating a Caco-2 permeability assay. We have shown that our approach is able to provide quantitative estimates of the initial slope of the concentration that can be compared with the experimental values, and reported a remarkable consistency

15



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

Figure 4: Upper panel: logP values estimated from Eq. 1 depend on the height of the freeenergy barrier, ∆G for the free energy profiles reported in the panels A-D. Lower panel: logP values estimated from Eq. 2 depend on the maximum free-energy difference ∆G + ∆H for the free energy profiles reported in the panels A-D. (A) A Gaussian barrier of positive/negative height G mimics the effect of the membrane on a hydrophilic/lipophilic molecule, respectively. (B) A more complicated model of the free-energy profile for the permeation of a lipophilic compound, with free-energy barriers and wells. (C-D) Peak-trough profiles.

16


Page 16 of 20

Page 17 of 20

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


that is apparently not affected by the complexity of the experimental setup. In particular, in the case of lipophilic drugs (ETR and EFV), the permeation process is not associated with a significant free energy barrier. In these cases, our calculations clearly show that the value of P , defined in terms of the steady state flux J (Eq. 1), differs significantly from that based on the initial rate of uptake, Eq. 2. Not surprisingly, the experiments are consistent with the estimate provided by Eq. 2. Therefore, the Caco-2 assay does not provide an estimate of the permeability coefficient defined in Eq. 1, unless the retention factor is correctly included in the estimate. In the case of hydrophilic ligands, both Eqs. 1 and 2 yield consistent values that are fully accounted for by our computational procedure. These conclusions have been generalized by our simplified one-dimensional Fokker-Planck formulation of the permeation process. The different free energy profiles, though not exhaustive, still represent a variety of cases that may be encountered in the permeation of several drug molecules. This simplified model is therefore useful to classify the permeation process of any compound and thus assess whether it would be desirable to use Eq. 1 or Eq. 2 to obtain its membrane permeability, by measuring it in an experiment, or by computation in all atoms molecular dynamics simulations.

Supporting Information Available Parametrization of EFV and ETR. Details of the computational setup, of the chosen collective variables, microstate determination, mass density of the membrane. PDB files of the representative structures in Fig. 2. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) K., M. Brody’s Human Pharmacology: Molecular to Clinical., 4th ed.; Elsevier, 2004. 17



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 18 of 20

(2) Smith, D. A. Mol. Pharmaceutics 2013, 10, 1162–1170. (3) Sugano, K.; Kansy, M.; Artursson, P.; Avdeef, A.; Bendels, S.; Di, L.; Ecker, G. F.; Faller, B.; Fischer, H.; Gerebtzoff, G.; Lennernaes, H.; Senner, F. Nat. Rev. Drug Discovery 2010, 9, 597–614. (4) Artursson, P.; Karlsson, J. Biochem. Biophys. Res. Commun. 1991, 175, 880–885. (5) Avdeef, A. Absorption and Drug Development: Solubility, Permeability, and Charge State; John Wiley & Sons, 2003. (6) Krishna, G.; Chen, K.; Lin, C.; Nomeir, A. A. Int. J. Pharm. 2001, 222, 77–89. (7) Fossati, L.; Dechaume, R.; Hardillier, E.; Chevillon, D.; Prevost, C.; Bolze, S.; Maubon, N. Int. J. Pharm. 2008, 360, 148–155. (8) Heikkinen, A. T.; Monkkonen, J.; Korjamo, T. J. Pharmacol. Exp. Ther. 2009, 328, 882–892. (9) Marrink, S. J.; ; Berendsen, H. J. C. J. Phys. Chem. 1996, 100, 16729–16738. (10) Orsi, M.; ; Essex, J. W. Soft Matter 2010, 6, 3797–3808. (11) Paloncýová, M.; Berka, K.; Otyepka, M. J. Phys. Chem. B 2013, 117, 2403–2410. (12) Ghaemi, Z.; Minozzi, M.; Carloni, P.; Laio, A. J. Phys. Chem. B 2012, 116, 8714–8721. (13) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562–12566. (14) Cada, D. J.; Levien, T.; Baker, D. E. Hosp. Pharm. 2008, 43, 498–510. (15) Takano, R.; Sugano, K.; Higashida, A.; Hayashi, Y.; Machida, M.; Aso, Y.; Yamashita, S. Pharm. Res. 2006, 23, 1144–1156. (16) Jorgensen, W.; J.D., M. J. Am. Chem. Soc. 1983, 105, 1407–1413.

18


Page 19 of 20

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


(17) Jojart, B.; Martinek, T. J. Comput. Chem. 2007, 28, 2051–2058. (18) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269– 10280. (19) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; ; Case., D. A. J. Comput. Chem. 2004, 25, 1157–1164. (20) Laio, A.; Gervasio, F. Rep. Prog. Phys. 2008, 71, 126601. (21) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (22) Epand, R. M.; Bottega, R. Biochim. Biophys. Acta, Biomembr. 1988, 944, 144 – 154. (23) Nosé, S. Mol. Phys. 1984, 52, 255–268. (24) Hoover, W. Phys. Rev. A 1985, 31, 1695–1697. (25) Hess, B.; Bekker, H.; Berendsen, H.; Fraaije, J. J. Comput. Chem. 1997, 18, 1463–1472. (26) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10093. (27) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182–7190. (28) Marinelli, F.; Pietrucci, F.; Laio, A.; Piana, S. PLoS Comput. Biol. 2009, 5, e1000452. (29) Minozzi, M.; Lattanzi, G.; Benz, R.; Costi, M. P.; Venturelli, A.; Carloni, P. PLoS One 2011, 6, e23187. (30) Scholler, M.; Hoetelmans, R.; Beets, G.; Vandermeulen, K.; Peeters, M.; Bastiaanse, L.; Leemans, R.; Debroye, C.; Woodfall, B. Poster Exhibition: The 3rd IAS Conference on HIV Pathogenesis and Treatment. Rio De Janeiro, Brazil. Abstract no. TuPe3.1B11 July 24- 27, 2005,

19



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

Graphical TOC Entry

20


Page 20 of 20

Permeability Coefficients of Lipophilic Compounds Estimated by

Recommend Documents