Revisiting the Meyer-Overton rule for drug-membrane permeabilities

ABSTRACT

Assessing the permeation rate of drug-like molecules across a lipid membrane is of paramount importance for pharmaceutical applications. While the Meyer-Overton rule relates the permeability coefficient to a few thermodynamic properties, its accuracy is limited by the homogeneity assumption. In this work, we extract an analogous relation for the permeation of a small solute through a lipid bilayer. This is obtained by systematically screening a subset of chemical space by means of high-throughput coarse-grained simulations, and relying on the accurate inhomogeneous solubility-diffusion model (ISDM). We connect the permeability coefficient of a compound to two molecular descriptors: partitioning free-energy and acid dissociation constant. In the ISDM the permeability is a functional of the potential of mean force. We discuss how – and when – combining together these profiles for a large variety of compounds can result in a smooth dependence of the permeability on the molecular descriptors. We focus on acidic molecules: for weak acids, the permeability largely depends on the main free-energy differences of the problem, which directly link to our molecular descriptors. For strong acids, the multivalued nature of the permeability requires to introduce a third variable, characterising the amphiphilicity of the small molecule.

Abbreviations: CG: coarse-grained; HTCG: high-throughput coarse-grained; PMF: potential of mean force; ISDM: inhomogeneous solubility-diffusion model; DOPC: 1; 2-Dioleoyl-sn-glycero-3-phosphocholine

GRAPHICAL ABSTRACT

KEYWORDS:

• Drug-membrane systems
• permeability coefficient
• potential of mean force
• Meyer-Overton rule

1. Introduction

The strong binding affinity of a small drug-like candidate to the specific intracellular target, alone, does not guarantee its efficacy. In addition to toxicity issues, the efficacy can be severely compromised by the presence of prohibitively large time scales for the compound to cross the lipid barrier enclosing the cell environment [1]. Understanding the transport mechanism of small molecules through lipid bilayers is therefore of fundamental importance for pharmaceutical applications.

Transmembrane transport processes are usually divided in two classes [2]: (i) carrier-mediated processes, and (ii) passive permeation, i.e. diffusion of the permeant driven only by its concentration gradient and thermal fluctuations. Although both classes can contribute to the total transport of a solute, passive permeation is believed to play a dominant role in the case of small drug-like molecules [3]. The permeability coefficient, P, quantifies the rate of passive permeation of a compound across the membrane. P is defined as [4] (1) $P=J/\mathrm{\Delta }u,$(1) where J is the steady-state flux per unit area of the solute arising from its concentration difference $\mathrm{\Delta }u$ between the two water compartments adjacent to the membrane. Equation (1(1) $P=J/\mathrm{\Delta }u,$(1) ) makes it apparent that a compound characterised by a low permeability coefficient will experience more difficulties in crossing the lipid membrane, in this way reducing the amount of substance delivered to the intended intracellular target.

Pharmacokinetic issues associated to insufficient permeation are among the major causes of high attrition rates in drug discovery [1]. Given the complexity of in vivo experiments, several in vitro methods have been developed to systematically estimate the permeability coefficient of small molecules, in order to screen them at early stages in the discovery process [5–7]. These systems have contributed to the expansion of permeability databases due to the possibility of performing high-throughput experiments. At the same time, they typically provide limited insight into the microscopic origin of permeation.

From a theoretical point of view, the first model of solute transport across a lipid bilayer is the Meyer-Overton rule [8,9]. It approximates the membrane as a homogeneous slab, for which an application of Fick's first law of diffusion provides (2) $P=KD/\sigma .$(2) In Equation (2(2) $P=KD/\sigma .$(2) ), σ is the thickness of the bilayer core, while D and K are the diffusivity and partitioning coefficient of the compound, respectively. In practice, K is usually further approximated by the partitioning coefficient between water and a bulk apolar solvent (e.g. octanol).

Despite its simplicity, the Meyer-Overton rule identifies key properties associated with the permeability: due to the dependence on the partitioning coefficient, apolar compounds are characterised by higher permeation rates [10]. However, its main limitation stems from the homogeneity assumption, at odds with the inherent asymmetry of a lipid bilayer. Accounting for the heterogeneity of the medium results in the inhomogenous solubility-diffusion model (ISDM) [11,12], which describes the diffusion process in terms of a one-dimensional Smoluchowsky equation along z, the normal distance of the compound from the bilayer midplane [13]. From the ISDM, the permeability coefficient of a compound can be calculated as an integral over the membrane extension of its potential of mean force (PMF), $G\left(z\right)$, and local diffusivity, $D\left(z\right)$, as (3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) with $\beta =1/k_{\mathrm{B}}T$. Equations (2(2) $P=KD/\sigma .$(2) ) and (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) become equivalent only in the case of a flat potential of mean force and constant diffusivity – i.e. if homogeneity holds. The ISDM paved the way for an in silico estimation of permeability coefficients, where $G\left(z\right)$ and $D\left(z\right)$ are determined by means of enhanced-sampling classical molecular dynamics simulations [13–19]. Predictions were shown to exhibit extremely high correlation with either in vitro or in vivo permeability measurements [15,20], making in silico methods a robust, complementary approach to experimental investigations. Unfortunately, the extensive computational resources required for convergence of the atomistic conformational sampling – roughly ${10}^5$ CPU hours per chemical compound – have so far limited these calculations to few small molecules with very specific chemistries, see e.g. Ref. [15,20].

The importance of the Meyer-Overton rule lied in its ability to relate the permeability of a compound to a key molecular descriptor directly linked to the thermodynamics of the process. However, the rule has been recently shown to overestimate the permeabilities calculated through the ISDM, which were observed to be nonlinear in the partitioning coefficient [19]. Moreover, in contrast to the ISDM the Meyer-Overton rule neglects the impact of the protonation state of a compound on its permeation rate [15].

Given the accuracy of the ISDM, being able to connect its predictions to a set of simple thermodynamic properties would be of fundamental importance to obtain rapid estimates of permeability coefficients. Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) highlights that P is a functional of the diffusivity and potential of mean force. For small molecules $D\left(z\right)$ was shown to be largely insensitive to the chemical detail [15], so that the permeation rate is mainly affected by variations in the PMF. Extracting the link between thermodynamics and permeability in the ISDM would require to extensively probe at least a representative subset of the spectrum of possible drug-like molecules, to systematically investigate the impact of the PMF on the permeability coefficient. The computational resources associated to PMF calculations at an atomistic resolution, combined with the estimated size of drug chemical space – ${10}^{60}$ compounds [21] – make this approach unfeasible in practice.

The overwhelming number of possible compounds thus calls for the introduction of new strategies that allow to more comprehensively cover large regions of chemical space. In this perspective, we recently proposed the use of coarse-grained (CG) models as a powerful instrument to preserve accurate thermodynamics while ‘clustering’ chemical space, in this way enabling a more systematic and efficient in silico exploration of chemical diversity [22,23].

In CG models, subsets of atoms are grouped together into elementary units (beads), with interactions tuned to capture a set of structural and/or thermodynamic properties of the original system [24]. In the Martini CG force field, the fundamental building blocks are a few bead types that aim at reproducing the partitioning free energy of small organic fragments between solvents of different polarity [25,26]. Due to the nature of the properties chosen as target, Martini proved to be very successful in reproducing the thermodynamics of insertion of small molecules in lipid bilayers [27,28], at the same time significantly reducing the computational investment [29].

The finite number of interaction types in the model is such that many chemically different small molecules with similar shape and polarity are mapped onto the same CG representation [22]. In contrast to traditional in silico methods at an atomistic resolution, the small number of bead types further enables a more systematic exploration of the possible variations in the thermodynamics of a small molecule: this can be obtained by combinatorially constructing all coarse-grained compounds up to a certain size (Figure 1(a)).

Figure 1. (a): we employ high-throughput coarse-grained simulations to screen a representative subset of the chemical space of small organic molecules. The systematic insertion of this large variety of compounds with very different polarities in a lipid bilayer results in an ensemble of PMFs, which dictate their permeability coefficient. (b): we highlight PMFs from the two main classes of compounds analyzed in this work: polar (top) and apolar (bottom). (c): Construction of the effective PMF $G_{\mathrm{T}}\left(z\right)$ (yellow area) contributing to the total resistivity $R_{\mathrm{T}}\left(z\right)$ starting from the neutral (red) and charged (blue) PMFs of a compound. We highlight the main free-energy differences of the problem: the water membrane free-energy $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ of the neutral compound, and the acid/base free energy difference $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$ dictating the shift between the two PMFs, directly linked to the ap$K_{\mathrm{a}}$.

We recently took advantage of these features to shed new light on the thermodynamics of insertion of a small solute in a lipid membrane. In Ref. [22] we employed high-throughput coarse-grained (HTCG) simulations to sample a representative subset of the chemical space of small organic molecules in the range 30−160 Da. Our comprehensive investigation highlighted that the overall shape of $G\left(z\right)$ can be reconstructed only from the partitioning free energy $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ – i.e. the free-energy difference of transferring the solute from bulk water to the bilayer midplane, see Figure 1(c).

In Ref. [23], by feeding HTCG results into Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) and accounting for the possible protonation states of a compound we identified a smooth relation expressing the permeability coefficient as a function of two molecular descriptors: $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and acid dissociation constant, p$K_{\mathrm{a}}$. This enabled permeability predictions for an unprecedented 500,000 small molecules including acidic, basic, and zwitterionic compounds, whose accuracy was validated against experimental as well as atomistic simulation results.

Importantly, the relation in Ref. [23] connects the permeation rate to only two thermodynamic properties, providing an alternative to the Meyer-Overton rule for predicting permeability coefficients in the framework of the accurate ISDM. Results were obtained by combining together permeabilities for compounds with highly different polarities, characterised by a large variety of neutral and charged potentials of mean force (Figure 1(a)). In this work, we analyze in more detail the impact of the molecular descriptors on the permeation rate, resulting – when possible – in a smooth dependence of the permeability coefficient on these thermodynamic variables. We here focus on acidic compounds, described by the acidic p$K_{\mathrm{a}}$ or ap$K_{\mathrm{a}}$, separately discussing neutral compounds ($\mathrm{ap}{K}_{\mathrm{a}}⪆8$), weak acids ($4⪆\mathrm{ap}{K}_{\mathrm{a}}\lesssim 8$) and strong acids ($\mathrm{ap}{K}_{\mathrm{a}}\lesssim 4$). For neutral compounds, the permeability is largely insensitive to the detailed shape of the underlying PMF, but mainly dictated by $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ (Figure 1(c)) [16]. This free-energy difference is not correctly accounted for in the Meyer-Overton rule, thus generating the observed discrepancies with ISDM predictions [19]. Moving towards weak acids, the ap$K_{\mathrm{a}}$ comes into play. However, the insensitivity of the permeability on the shape of the PMF suggests that the permeation rate mainly depends on a combination of two free-energy differences: $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$, where the latter is the free energy of neutralising the compound in bulk water, directly linked to the ap$K_{\mathrm{a}}$ (Figure 1(c)). This results in the smoothness of the permeability coefficient observed as a function of our thermodynamic variables. Finally, for strong acids, at fixed ap$K_{\mathrm{a}}$ the permeability coefficient does no longer smoothly depend on $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. This requires to introduce a third variable characterising the degree of amphiphilicity of the small molecule. In this case, our results provide a lower bound for the permeation rate of a compound: at a fixed ap$K_{\mathrm{a}}$, the permeability coefficient cannot be lower than a value dictated by the combination of the two free-energy differences, $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$.

2. Methods

2.1. Molecular dynamics simulations

Molecular dynamics simulations of the Martini force field [25,26,30] were performed in Gromacs 4.6.6 [31]. We employed an integration timestep $\delta t=0.02\tau$, where τ is the model's natural unit of time. Control over the system temperature ($T=300\mathrm{K}$) and pressure ($P=1\mathrm{bar}$) was obtained by means of a Parrinello-Rahman barostat [32] and a stochastic velocity rescaling thermostat [33], with coupling constant ${\tau }_P=12\tau$ and ${\tau }_T=\tau$ respectively. A membrane of $36{\mathrm{nm}}^2$ containing N=128 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipids (64 per layer), $N^{\prime }=1890$ water molecules, $N^{″}=190$ antifreeze particles [26] and enough counterions to ensure charge neutrality was generated by means of the Insane building tool [34], and subsequently minimised, heated up, and equilibrated.

We determined potentials of mean force $G\left(z\right)$ by means of umbrella sampling simulations [35]. For each analyzed compound, we employed 24 windows with harmonic biasing potentials ($k=240\mathrm{kcal}/\mathrm{mol}/{\mathrm{nm}}^2$) centred every $0.1\mathrm{nm}$ along the normal to the bilayer midplane. In each window, two solute molecules were placed in the membrane in order to increase sampling and alleviate leaflet-area asymmetry [36,37]. Each production simulation was run for $1.2\cdot {10}^5\tau$, and potentials of mean force were extracted by means of the weighted histogram analysis method [38–40].

2.2. Permeability coefficients

The permeability coefficient of a compound is obtained as an integral of its resistivity $R\left(z\right)$ over the membrane normal, see Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ). While CG simulations provide access to the potential of mean force $G\left(z\right)$, the acceleration of the dynamics introduced by CG models [41] makes the estimation of the inhomogenous diffusivity $D\left(z\right)$ a priori difficult. However, atomistic simulations of several small molecules in a DOPC membrane showed that $D\left(z\right)$ is largely insensitive to the chemical detail [15]. Moreover, $D\left(z\right)$ only provides a logarithmic correction to the ${\log }_{10}P$ of a compound, see Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ). In our calculations we thus relied on a unique diffusivity profile, and parametrised it as (4) $D\left(z\right)=\alpha +\frac{\beta }{{\mathrm{e}}^{-\gamma (x-\delta )}+1},$(4) with $\alpha =0.85\cdot {10}^{-6}{\mathrm{cm}}^2/\mathrm{s}$, $\beta =9.15\cdot {10}^{-6}{\mathrm{cm}}^2/\mathrm{s}$, $\gamma =7.5{\mathrm{nm}}^{-1}$, $\delta =3\mathrm{nm}$, which correctly captures the main features of its atomistic counterpart [15].

In the case of acidic compounds, accounting for deprotonation requires to slightly modify Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) to include the flux of both neutral and charged species. The total resistivity reads [15] $R_{\mathrm{T}}(z{)}^{-1}=R_{\mathrm{N}}(z{)}^{-1}+R_{\mathrm{C}}(z{)}^{-1}$, where $R_{\mathrm{N}}$ and $R_{\mathrm{C}}$ are the resistivities of the neutral and charged forms of the compound, respectively. The corresponding potentials of mean force $G_{\mathrm{N}}\left(z\right)$ and $G_{\mathrm{C}}\left(z\right)$ must be offset by the acid/base free-energy difference in bulk water [36] (5) $\mathrm{\Delta }{G}_{\mathrm{b}\to a}=G_{\mathrm{acid}}-G_{\mathrm{base}}=k_{\mathrm{B}}T(\mathrm{pH}-\mathrm{ap}{K}_{\mathrm{a}})\ln 10,$(5) where we consider the case of neutral $\mathrm{pH}=7.4$. Assuming that the diffusivity does not depend on the protonation state of the compound, the combination of the two PMFs in the total resistivity $R_{\mathrm{T}}$ can be represented by an effective potential $G_{\mathrm{T}}\left(z\right)=min\left[{\overline{G}}_{\mathrm{N}}\right(z),{\overline{G}}_{\mathrm{C}}(z\left)\right]$ for all z, see Figure 1(c), from which $R_{\mathrm{T}}\left(z\right)\propto \exp \left[\beta G_{\mathrm{T}}\right(z\left)\right]$. ${\overline{G}}_i\left(z\right),i=\mathrm{N},\mathrm{C}$ indicate the shifted counterparts of the neutral and charged PMFs, respectively.

Permeability coefficients predicted by the CG model were validated in Ref. [23] against atomistic simulations results as well as experimental measurements for several chemically different small molecules covering a wide range of p$K_{\mathrm{a}}$. In both cases we observed a mean absolute error well within one ${\log }_{10}$ unit, demonstrating the accuracy of the CG model in reproducing the permeation rate of a compound.

3. Results

We considered all possible CG compounds consisting of a single Martini bead (‘unimers’) as well as those generated by binding together all possible combinations of two beads (‘dimers’), representative of a subset of the chemical space of small organic molecules in the range 30–160 Da [22]. We calculated the PMF $G_{\mathrm{N}}\left(z\right)$ for all neutral unimers and dimers, 119 in total. For each neutral compound, we then determined $G_{\mathrm{C}}\left(z\right)$ for its deprotonated counterpart. At the CG level, deprotonating a chemical group amounts to replacing the original bead type with a negatively charged one. We always represent deprotonated fragments by means of a $Q_{\mathrm{a}}$ type, following the standard Martini notation. In the case of unimers, this choice automatically determines the CG representation of the deprotonated form of a compound. In the case of dimers, we initially assumed that deprotonation always occurs in the chemical fragment represented by the bead of higher polarity, but relax this assumption along the discussion, vide infra. By combining $G_{\mathrm{N}}\left(z\right)$ and $G_{\mathrm{C}}\left(z\right)$ in the total resistivity $R_{\mathrm{T}}\left(z\right)$, we calculated the permeability coefficient of every compound as a function of the ap$K_{\mathrm{a}}$ every 0.2 ap$K_{\mathrm{a}}$ unit. The set of permeability data was then projected on the $(\mathrm{ap}{K}_{\mathrm{a}},\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}})$ plane.

To analyze how the resulting permeabilities are affected by a change in the two molecular descriptors, we separately discuss the case of neutral compounds/weak acids ($\mathrm{ap}{K}_{\mathrm{a}}⪆4$) and strong acids ($\mathrm{ap}{K}_{\mathrm{a}}\lesssim 4$).

3.1. Neutral compounds/weak acids: revisiting the Meyer-Overton rule

Figure 2 displays a two-dimensional projection of the permeability data on the $({\log }_{10}P,\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}})$ plane for neutral compounds and weak acids, $\mathrm{ap}{K}_{\mathrm{a}}⪆4$. We present results for a subset of ap$K_{\mathrm{a}}$ values, and colour the data accordingly. Panels (a) and (b) distinguish between the permeabilities obtained for Martini unimers and dimers, effectively amounting to a segregation in molecular weights [22].

Figure 2. Two-dimensional projection of the permeability data on the $({\log }_{10}P,\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}})$ plane in the case of neutral compounds and weak acids, $\mathrm{ap}{K}_{\mathrm{a}}⪆4$. Data are coloured according to the value of the ap$K_{\mathrm{a}}$. Panels (a) and (b) distinguish between results obtained for Martini unimers and dimers. For neutral compounds with $\mathrm{ap}{K}_{\mathrm{a}}⪆8$, the permeability becomes independent on the value of the acid dissociation constant. Accordingly, we only show data up to $\mathrm{ap}{K}_{\mathrm{a}}=10$. We present permeability coefficients obtained from HTCG simulations (points) and the corresponding interpolated permeability profiles (lines). For neutral compounds, we include predictions from the Meyer-Overton rule in the presence of flat potentials of mean force with heights equal to $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ (‘Overton-mem’, dashed lines), see text.

Due to the finite number of interaction types in the Martini model, simulation results are located at a discrete set of partitioning free energies. In both cases, the regularity of the points at fixed ap$K_{\mathrm{a}}$ as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ allows to reconstruct continuous permeability profiles. We include these interpolations in the corresponding panels of Figure 2. The independence of the interpolated curves on the CG representation, together with the dense covering in $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ generated by CG dimers, highlight the smoothness of the permeability coefficient as a function of ap$K_{\mathrm{a}}$ and $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. These results enable the calculation of the permeation rate only starting from two thermodynamic properties, providing a phenomenological alternative to the Meyer-Overton rule in the framework of the accurate ISDM.

The permeability coefficient weights the contributions of both the charged $G_{\mathrm{C}}\left(z\right)$ and neutral $G_{\mathrm{N}}\left(z\right)$ of a compound in an effective potential $G_{\mathrm{T}}\left(z\right)$ which approximately dictates the shape of the total resistivity, $\ln R_{\mathrm{T}}\left(z\right)\sim \beta G_{\mathrm{T}}\left(z\right)$ (Figure 1(c)). Understanding the smooth dependence of the permeability on both descriptors thus requires to more carefully analyze how the underlying potentials of mean force combine in $G_{\mathrm{T}}\left(z\right)$. For this purpose, within the class of weak acids we separately discuss the case of neutral compounds – or ‘extremely weak’ acids – with $\mathrm{ap}{K}_{\mathrm{a}}⪆8$, and weak acids, $4\lesssim \mathrm{ap}{K}_{\mathrm{a}}\lesssim 8$.

3.1.1 Neutral compounds: $ap{K}_{\mathrm{a}}⪆8$

For neutral compounds, the permeability is independent of the acid dissociation constant. This stems from the fact that by increasing the ap$K_{\mathrm{a}}$ beyond $\mathrm{ap}{K}_{\mathrm{a}}\approx 8$ the charged repulsive PMF $G_{\mathrm{C}}\left(z\right)$ is shifted towards increasingly higher values, and does not contribute to the effective PMF (Figure 1(c)). One obtains $G_{\mathrm{T}}\left(z\right)\simeq G_{\mathrm{N}}\left(z\right)$, such that $R_{\mathrm{T}}\left(z\right)\sim \exp \left[\beta G_{\mathrm{N}}\right(z\left)\right]$.

For $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}\lesssim 0$, Figure 2 shows that ${\log }_{10}P$ roughly plateaus. The reason is straightforward: apolar compounds are characterised by a negative, attractive well of $G_{\mathrm{N}}\left(z\right)$ within the membrane core (Figure 1(b)). This well gives a negligible contribution to the integral in Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) compared to the regions in which the PMF reaches saturation – i.e. as the compound approaches the bulk water environment.

Moving towards $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}⪆0$ the repulsive nature of $G_{\mathrm{N}}\left(z\right)$ (Figure 1(b)) increases the resistivity $R_{\mathrm{T}}\left(z\right)$, resulting in a reduction of the permeation rate. In this case, Figure 2 highlights a linear dependence of ${\log }_{10}P$ as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. It was recently observed that for polar compounds the permeability is mainly dictated by the height of the repulsive PMF (i.e. $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$) irrespective of its detailed shape – e.g. width, presence of attractive tails [16]. This suggests the linearity in Figure 2 to be a consequence of this insensitivity.

We thus approximated the membrane as a bulk homogeneous solution, and employed the Meyer-Overton rule to calculate the corresponding permeabilities. In Equation (2(2) $P=KD/\sigma .$(2) ), we assumed $K=K_{\mathrm{WM}}=\exp [-\beta \mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}]$ and set $\sigma ,D$ equal to the average thickness and diffusivity values. The use of $K_{\mathrm{WM}}$ amounts to replacing an ensemble of repulsive PMFs smoothly approaching zero by square-well potentials with heights equal to $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. We stress that $K_{\mathrm{WM}}$ differs from the partitioning coefficient K in Equation (2(2) $P=KD/\sigma .$(2) ), which is instead calculated accounting for the inhomogeneous nature of the membrane, vide infra.

Permeabilities obtained from this homogenous model are included in both panels of Figure 2 (dashed lines). In the case of polar compounds results only slightly underestimate ISDM ones, the maximum deviation being roughly one ${\log }_{10}$ unit. This confirms that for $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}⪆0$ the permeability mainly depends on the height of the repulsive PMF [16], and explains the linearity of ${\log }_{10}P$ as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. On the contrary, for $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}\lesssim 0$ the homogeneous model strongly overestimates the permeability coefficient: this stems from the fact that the square-well potentials neglect the saturation of the PMFs that occur by approaching the bulk water environment – the only region that for apolar compounds contributes to the resistivity.

This homogeneous model differs from the Meyer-Overton rule, in which the partitioning coefficient K of a compound reads (6) $K=\frac{1}{\sigma }\int \mathrm{d}z\exp [-\beta G_{\mathrm{N}}(z\left)\right].$(6) Predictions from the rule were recently shown to overestimate ISDM ones [19], and the permeability coefficient was observed to be multivalued in K. By approximating $D\left(z\right)$ in Equation (3(3) $P^{-1}=\int \mathrm{d}zR\left(z\right)=\int \mathrm{d}z\frac{\exp \left[\beta G\right(z\left)\right]}{D\left(z\right)},$(3) ) with its average value, a comparison with Equation (2(2) $P=KD/\sigma .$(2) ) highlights that in the latter we are replacing the inverse of the integral of $\exp \left[\beta G_{\mathrm{N}}\right(z\left)\right]$ with the integral of its inverse. The regions contributing to the two integrals are essentially swapped: for apolar compounds the negative well becomes the major contributor to Equation (6(6) $K=\frac{1}{\sigma }\int \mathrm{d}z\exp [-\beta G_{\mathrm{N}}(z\left)\right].$(6) ), while in the ISDM its contribution is negligible compared to the region of saturation of the PMF. The opposite occurs in the case of polar compounds, where the free-energy barrier $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ that dictates permeation rate in the ISDM does now not contribute to Equation (6(6) $K=\frac{1}{\sigma }\int \mathrm{d}z\exp [-\beta G_{\mathrm{N}}(z\left)\right].$(6) ), which is instead governed by the regions where the PMF reaches saturation or exhibits attractive tails [19]. This lies at the origin of the failure of the Meyer-Overton rule in reproducing ISDM permeabilities for neutral compounds.

3.1.2 Weak acids: $4\lesssim ap{K}_{\mathrm{a}}\lesssim 8$

At fixed $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, Figure 2 shows that weak acids with $4\lesssim \mathrm{ap}{K}_{\mathrm{a}}\lesssim 8$ are characterised by a reduction in the permeability coefficient with respect to neutral compounds. In addition to polarity, the deprotonation of a small molecule affects its permeation rate.

By decreasing the ap$K_{\mathrm{a}}$, both the neutral and charged PMFs of a compound contribute to the total resistivity. The neutral PMF $G_{\mathrm{N}}\left(z\right)$ is shifted towards increasingly higher values by a factor of $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$, and the minimum between the two profiles dictates the shape of $G_{\mathrm{T}}\left(z\right)$ (Figure 1(c) and Section 2.2), with $R_{\mathrm{T}}\left(z\right)\sim \exp \left[\beta G_{\mathrm{T}}\right(z\left)\right]$. We performed this construction for all compounds investigated in this work, resulting in the combination of a large variety of PMFs. However, Figure 2 highlights that at fixed ap$K_{\mathrm{a}}$ the only differences between the permeabilities of acidic and neutral compounds as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ consist in the length of the initial plateau and in the presence of a negative offset for the linear regime.

To understand the origin of this behaviour, Figure 3 displays two examples of the construction of $G_{\mathrm{T}}\left(z\right)$ for a polar and an apolar compound with $\mathrm{ap}{K}_{\mathrm{a}}=5$. For $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}⪆0$ (main plot in Figure 3), combining the neutral and charged PMFs results in a repulsive $G_{\mathrm{T}}\left(z\right)$ with an effective barrier given by the sum of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$. As for polar compounds the permeability coefficient is mainly dictated by the height of the PMF, increasing this by a factor $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$ – which is constant at fixed ap$K_{\mathrm{a}}$ and independent of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, see Figure 1(c) – preserves the linearity of ${\log }_{10}P$ in $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ but generates the offsets observed in Figure 2.

Figure 3. Construction of $G_{\mathrm{T}}\left(z\right)$ (blue) for an acidic compound from its charged $G_{\mathrm{C}}\left(z\right)$ (red) and neutral $G_{\mathrm{N}}\left(z\right)$ (orange) PMFs. In calculating $G_{\mathrm{T}}\left(z\right)$, $G_{\mathrm{N}}\left(z\right)$ is vertically shifted by $\mathrm{\Delta }{G}_{\mathrm{b}\to a}=k_{\mathrm{B}}T(\mathrm{pH}-\mathrm{p}{K}_{\mathrm{a}})\ln 10$, and we consider $\mathrm{ap}{K}_{\mathrm{a}}=5$. We present results for a polar (main plot) and an apolar (inset) dimer.

For $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}\lesssim 0$ (inset in Figure 3), the effect of a small vertical shift of $G_{\mathrm{N}}\left(z\right)$ for a fixed ap$K_{\mathrm{a}}$ is to reduce the attractive well that in the neutral case was not contributing to the resistivity. By increasing $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, the well crosses zero and $G_{\mathrm{T}}\left(z\right)$ becomes repulsive. At this stage, the permeability starts following the linear regime. This explains the reduction in the length of the plateau observed in the permeability profiles for $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}\lesssim 0$ by increasing the strength of the acid.

Overall, our results suggest that in the case of weak acids the permeability is mainly dictated by the combination of two free-energy differences: the water/membrane partitioning free energy $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, and the acid/base free-energy difference $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$, which directly links to the ap$K_{\mathrm{a}}$ (Figure 1(c)). This is observed irrespective of the detailed shape of the underlying PMFs, and explains the smooth dependence of the permeability observed as a function of the two molecular descriptors.

The construction of $G_{\mathrm{T}}\left(z\right)$ allows for a further discussion about the role of the charged PMF of a compound in the permeation rate. Figure 3 shows that if $G_{\mathrm{C}}\left(z\right)$ has reached the shifted neutral PMF within the membrane core – a saturation condition always satisfied for $\mathrm{ap}{K}_{\mathrm{a}}⪆4$ – an increase in its amplitude does not affect $G_{\mathrm{T}}\left(z\right)$. It follows that the permeability becomes largely independent on the polarity of the deprotonated form of the compound.

For a small molecule mapping to a CG dimer we originally assumed that deprotonation always occurred in the chemical fragment represented by the bead of higher polarity. We thus recalculated permeability coefficients relying on the opposite assumption – i.e. deprotonating the less polar bead. The resulting charged PMFs are characterised by a stronger repulsion within the membrane core, an example being shown in the inset of Figure 4.

Figure 4. Main plot: Permeability data of CG dimers (${\log }_{10}$ scale) calculated by relying on different CG representations of the deprotonated form of a compound. Starting from a neutral dimer, coloured points are obtained by assuming that deprotonation always occur in the chemical fragment encapsulated in the bead of higher polarity. In dashed lines (‘Apolar-dep’) we always deprotonate the bead of lower polarity. Inset: Comparison of the charged PMFs $G_{\mathrm{C}}\left(z\right)$ arising from the two different CG representations of the deprotonated form of the neutral polar dimer in Figure 3. The PMF of the compound obtained by deprotonating the bead of lower polarity is characterised by a stronger repulsion.

Permeability coefficients obtained from the two deprotonation prescriptions are presented in Figure 4 (main plot). Results are essentially indistinguishable. Slight deviations are observed for $\mathrm{ap}{K}_{\mathrm{a}}=4$, in which the condition of saturation of $G_{\mathrm{C}}\left(z\right)$ is not always satisfied as we are moving towards strong acids – vide infra. Importantly, this analysis shows that while the permeation rate could in principle depend on which chemical fragment of the small molecule is subject to deprotonation, this does not matter in practice.

3.2. Strong acids: the role of amphiphilicity

Panel (a) of Figure 5 displays permeability data of CG dimers as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ for $\mathrm{ap}{K}_{\mathrm{a}}=2$, together with some previously analyzed results for $\mathrm{ap}{K}_{\mathrm{a}}⪆4$. At $\mathrm{ap}{K}_{\mathrm{a}}=2$, data initially follow a linear trend analogous to the one of weak acids due to the saturation of $G_{\mathrm{C}}\left(z\right)$ in the effective PMF. However, for $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}⪆0$ we observe a spread of the points over the $({\log }_{10}P,\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}})$ plane: the permeability coefficient at fixed ap$K_{\mathrm{a}}$ becomes multivalued in $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. This does not occur in the case of CG unimers, whose ${\log }_{10}P$ continue to behave linearly as $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ increases (data not shown).

Figure 5. Panel (a): Two-dimensional projection of the permeability data on the $({\log }_{10}P,\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}})$ plane in the case of strong acids with $\mathrm{ap}{K}_{\mathrm{a}}=2$. We include as reference a subset of the data obtained for weak acids, $\mathrm{ap}{K}_{\mathrm{a}}⪆4$. Points were obtained by assuming that starting from a neutral dimer, deprotonation always occur in the chemical fragment encapsulated in the bead of higher polarity. On the contrary, in dashed lines (‘Apolar-dep’) we always deprotonate the bead of lower polarity. Panel (b): Comparison of the neutral (main plot) and charged (inset) PMFs $G_{\mathrm{N}}\left(z\right)$ and $G_{\mathrm{C}}\left(z\right)$ of a $P_5{C}_1$ and a $P_2{C}_5$ dimers. In both cases, charged PMFs are obtained by assuming that deprotonation always occurs in the bead of higher polarity.

In strong acids, $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$ is such that neutral PMFs are shifted towards high values with respect to their charged counterparts (Figure 1(c)): only $G_{\mathrm{C}}\left(z\right)$ contributes to the effective PMF, and $R_{\mathrm{T}}\left(z\right)\sim \exp \left[\beta G_{\mathrm{C}}\right(z\left)\right]$. Different permeabilities for a given $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ thus arise from compounds presenting very different polarities in their depronotonated forms, despite having similar ones in the neutral case – i.e. similar $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$.

Martini assumes additivity in the partioning free energy of its constituting beads. Many possible combinations of two neutral beads can result in the same $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$: while these will lead to roughly the same free-energy difference for the neutral PMFs $G_{\mathrm{N}}\left(z\right)$, this does not necessarily hold for their deprotonated counterparts. As an example we consider a $P_5{C}_1$ dimer, composed by the most polar (P) and most apolar (C) Martini bead types, and a $P_2{C}_5$ dimer, in which we reduce the asymmetry in polarity between the two beads. Due to additivity, the compounds are characterised by similar neutral PMFs $G_{\mathrm{N}}\left(z\right)$ (main plot in Figure 5(b)). Assuming that deprotonation always occurs in the bead of higher polarity, however, the charged form of the first compound shows a weaker repulsion from the membrane core (inset of Figure 5(b)). As $R_{\mathrm{T}}\left(z\right)\sim \exp \left[\beta G_{\mathrm{C}}\right(z\left)\right]$, $P_5{C}_1$ is characterised by a higher ${\log }_{10}P$: $-4.20$ compared to $-6.0$ for $P_2{C}_5$.

The sensitivity of the permeability coefficient on the asymmetry in polarity between the chemical fragments composing the molecule suggests to introduce a third molecular descriptor. The most natural choice is the difference between the partitioning free energies of the two beads, $\mathrm{\Delta }{G}_{\mathrm{asy}}=\mathrm{\Delta }{G}_2-\mathrm{\Delta }{G}_1$, where beads are arranged in order of decreasing polarity. For a given $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, increasing $\mathrm{\Delta }{G}_{\mathrm{asy}}$ amounts to decreasing/increasing the polarity of the more apolar/polar chemical group. As such, $\mathrm{\Delta }{G}_{\mathrm{asy}}$ is essentially quantifying the amphiphilicity of the small molecule.

Permeability data for $\mathrm{ap}{K}_{\mathrm{a}}=2$ as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{asy}}$ are presented in Figure 6. The introduction of $\mathrm{\Delta }{G}_{\mathrm{asy}}$ completely solves the degeneracy in $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ observed in Figure 5(a): for a given $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ of the neutral form of the acid, increasing $\mathrm{\Delta }{G}_{\mathrm{asy}}$ generates an increase in the permeation rate. On the one hand, this is not surprising, meaning that the reduction in permeability generated by the deprotonation of a functional group can be mitigated by coupling it with a highly apolar one. However, in analyzing these results one has to account for the tendency of Martini to underestimate the PMF of charged compounds [27,42]. We thus repeated our calculations by assuming that deprotonation always occur in the bead of lower polarity, which as discussed in Section 3.1 generates an increase in the repulsive barrier of $G_{\mathrm{C}}\left(z\right)$. Results are included in Figure 5(a) (dashed line). In this case, the permeability coefficient follows a linear regime analogous to the case of weak acids due to saturation of the charged PMFs in $G_{\mathrm{T}}\left(z\right)$. We conclude that if the permeability coefficient depends on the amphiphilicity of the small molecule, due to the underestimation of the PMF of charged compounds in Martini this will be observed at extremely low values of the ap$K_{\mathrm{a}}$. For strong acids, however, our results provide a lower bound for the permeation rate: at a fixed ap$K_{\mathrm{a}}$, the permeability coefficient cannot be lower than the linear relation dictated by combining two free-energy differences, $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$.

Figure 6. Permeability coefficients (${\log }_{10}$ scale) for CG dimers in the case of strong acids with $\mathrm{ap}{K}_{\mathrm{a}}=2$. We display results as a function of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and the asymmetry $\mathrm{\Delta }{G}_{\mathrm{asy}}$. $\mathrm{\Delta }{G}_{\mathrm{asy}}$ is calculated as the difference between the water/membrane partitioning free energy of the beads composing the dimer, sorting them in order of decreasing polarity. In all cases, the charged form of the compound was obtained by assuming that deprotonation always occur in the bead of higher polarity.

4. Conclusions

While the Meyer-Overton rule links the permeability of a compound across a lipid bilayer to the partitioning coefficient, its validity is limited by the homogeneity assumption. An analogous relation in the framework of the accurate inhomogeneous solubility-diffusion model (ISDM) was missing up to now. In the ISDM, the permeability coefficient is a functional of the potential of mean force (PMF). Identifying this relation would require to probe large regions of chemical space to systematically investigate the impact of the PMF on the permeation rate. This is hampered by the overwhelming number of possible compounds.

In this work, we rely on CG simulations to systematically explore the permeability variations in a subset of the chemical space of small organic molecules, focussing on acidic compounds. We connect permeability coefficients calculated from the ISDM to two thermodynamic properties: the water/membrane free energy, $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$, and acid dissociation constant, ap$K_{\mathrm{a}}$. Results are obtained by probing a large variety of small molecules with very different PMFs, both neutral and charged. We carefully analyze how – and when – combining together these profiles can result in a smooth dependence of the permeability on the two thermodynamic variables.

For neutral compounds, the permeation rate is independent on the ap$K_{\mathrm{a}}$ and smoothly depends on $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$. This is due to the fact that the permeability coefficient is mainly dictated by $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ irrespective of the detailed shape of the neutral PMF. The deviation of the Meyer-Overton rule from ISDM predictions stems from incorrectly accounting for this free-energy difference. For weak acids, the deprotonation of a chemical group affects the permeation rate of a compound. In this case, charged and neutral PMFs combine in an effective profile whose water/membrane free-energy difference is the sum of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and the acid/base free-energy difference $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$, the latter being a linear function of the ap$K_{\mathrm{a}}$. The insensitivity of the permeability coefficient on the detailed shape of the effective PMF results in a smooth dependence of the permeation rate on our thermodynamic variables. Interestingly, we show that in weak acids the permeability is independent of the chemical fragment subject to deprotonation. For strong acids, the permeability is no longer a smooth function of the two original molecular descriptors. This requires to introduce a third orthogonal variable characterising the amphiphilicity of a small molecule. In this case, our results provide a lower bound to the permeation rate of a compound: the permeability cannot be smaller than a linear relation generated by the combination of $\mathrm{\Delta }{G}_{\mathrm{W}\to \mathrm{M}}$ and $\mathrm{\Delta }{G}_{\mathrm{b}\to a}$.

Acknowledgments

The authors thank Alessia Centi and Omar Valsson for critical reading of the manuscript and useful comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Roberto Menichetti  http://orcid.org/0000-0002-5961-6438

Tristan Bereau  http://orcid.org/0000-0001-9945-1271

Additional information

Funding

This work was supported by the Emmy Noether program of the Deutsche Forschungsgemeinschaft (DFG). The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA at Jülich Supercomputing Centre (JSC).