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


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.
as [4] P = J/ u, ( 1 ) where J is the steady-state flux per unit area of the solute arising from its concentration difference u between the two water compartments adjacent to the membrane. Equation (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][6][7]. These systems have contributed to the expansion of permeability databases due to the possibility of performing highthroughput 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 In Equation (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(z), and local diffusivity, D(z), as with β = 1/k B T. Equations (2) and (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(z) and D(z) are determined by means of enhanced-sampling classical molecular dynamics simulations [13][14][15][16][17][18][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) highlights that P is a functional of the diffusivity and potential of mean force. For small molecules D(z) 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]. contributing to the total resistivity R T (z) 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 G W→M of the neutral compound, and the acid/base free energy difference G b→a dictating the shift between the two PMFs, directly linked to the apK a .
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)).
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(z) can be reconstructed only from the partitioning free energy G W→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) 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: G W→M and acid dissociation constant, pK 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 pK a or apK a , separately discussing neutral compounds (apK a 8), weak acids (4 apK a 8) and strong acids (apK a 4). For neutral compounds, the permeability is largely insensitive to the detailed shape of the underlying PMF, but mainly dictated by G W→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 apK 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: G W→M and G b→a , where the latter is the free energy of neutralising the compound in bulk water, directly linked to the apK 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 apK a the permeability coefficient does no longer smoothly depend on G W→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 apK a , the permeability coefficient cannot be lower than a value dictated by the combination of the two free-energy differences, G W→M and G b→a .

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 δt = 0.02 τ , where τ is the model's natural unit of time. Control over the system temperature (T = 300 K) and pressure (P = 1 bar) was obtained by means of a Parrinello-Rahman barostat [32] and a stochastic velocity rescaling thermostat [33], with coupling constant τ P = 12 τ and τ T = τ respectively. A membrane of 36 nm 2 containing N = 128 1,2-dioleoylsn-glycero-3-phosphocholine (DOPC) lipids (64 per layer), N = 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(z) by means of umbrella sampling simulations [35]. For each analyzed compound, we employed 24 windows with harmonic biasing potentials (k = 240 kcal/mol/nm 2 ) centred every 0.1 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 · 10 5 τ , and potentials of mean force were extracted by means of the weighted histogram analysis method [38][39][40].

Permeability coefficients
The permeability coefficient of a compound is obtained as an integral of its resistivity R(z) over the membrane normal, see Equation (3). While CG simulations provide access to the potential of mean force G(z), the acceleration of the dynamics introduced by CG models [41] makes the estimation of the inhomogenous diffusivity D(z) a priori difficult. However, atomistic simulations of several small molecules in a DOPC membrane showed that D(z) is largely insensitive to the chemical detail [15]. Moreover, D(z) only provides a logarithmic correction to the log 10 P of a compound, see Equation (3). In our calculations we thus relied on a unique diffusivity profile, and parametrised it as with α = 0.85 · 10 −6 cm 2 /s, β = 9.15 · 10 −6 cm 2 /s, γ = 7.5 nm −1 , δ = 3 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) to include the flux of both neutral and charged species. The total resistivity reads [15] where R N and R C are the resistivities of the neutral and charged forms of the compound, respectively. The corresponding potentials of mean force G N (z) and G C (z) must be offset by the acid/base free-energy difference in bulk water [36] (5) where we consider the case of neutral 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 T can be represented by an effective potential 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 pK 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.

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 N (z) for all neutral unimers and dimers, 119 in total. For each neutral compound, we then determined G C (z) 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 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 N (z) and G C (z) in the total resistivity R T (z), we calculated the permeability coefficient of every compound as a function of the apK a every 0.2 apK a unit. The set of permeability data was then projected on the (apK a , G W→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 (apK a 4) and strong acids (apK a 4). Figure 2 displays a two-dimensional projection of the permeability data on the (log 10 P, G W→M ) plane for neutral compounds and weak acids, apK a 4. We present results for a subset of apK 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].

Neutral compounds/weak acids: revisiting the Meyer-Overton rule
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 apK a as a function of G W→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 G W→M generated by CG dimers, highlight the smoothness of the permeability coefficient as a function of apK a and G W→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 C (z) and neutral G N (z) of a compound in an effective potential G T (z) which approximately dictates the shape of the total resistivity, ln R T (z) ∼ βG T (z) (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 T (z). For this purpose, within the class of weak acids we separately discuss the case of neutral compounds -or 'extremely weak' acids -with apK a 8, and weak acids, 4 apK a 8. Figure 2. Two-dimensional projection of the permeability data on the (log 10 P, G W→M ) plane in the case of neutral compounds and weak acids, apK a 4. Data are coloured according to the value of the apK a . Panels (a) and (b) distinguish between results obtained for Martini unimers and dimers. For neutral compounds with apK a 8, the permeability becomes independent on the value of the acid dissociation constant. Accordingly, we only show data up to apK 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 G W→M ('Overton-mem', dashed lines), see text.

Neutral compounds: apK a 8
For neutral compounds, the permeability is independent of the acid dissociation constant. This stems from the fact that by increasing the apK a beyond apK a ≈ 8 the charged repulsive PMF G C (z) is shifted towards increasingly higher values, and does not contribute to the effective PMF (Figure 1(c)). One obtains For G W→M 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 N (z) within the membrane core (Figure 1(b)). This well gives a negligible contribution to the integral in Equation (3) compared to the regions in which the PMF reaches saturation -i.e. as the compound approaches the bulk water environment.
Moving towards G W→M 0 the repulsive nature of G N (z) (Figure 1(b)) increases the resistivity R T (z), 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 G W→M . It was recently observed that for polar compounds the permeability is mainly dictated by the height of the repulsive PMF (i.e. G W→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), we assumed K = K WM = exp[−β G W→M ] and set σ , D equal to the average thickness and diffusivity values. The use of K WM amounts to replacing an ensemble of repulsive PMFs smoothly approaching zero by square-well potentials with heights equal to G W→M . We stress that K WM differs from the partitioning coefficient K in Equation (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 G W→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 G W→M . On the contrary, for G W→M 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 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(z) in Equation (3) with its average value, a comparison with Equation (2) highlights that in the latter we are replacing the inverse of the integral of exp[βG N (z)] 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), 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 G W→M that dictates permeation rate in the ISDM does now not contribute to Equation (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.

Weak acids: 4 apK a 8
At fixed G W→M , Figure 2 shows that weak acids with 4 apK a 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 apK a , both the neutral and charged PMFs of a compound contribute to the total resistivity. The neutral PMF G N (z) is shifted towards increasingly higher values by a factor of G b→a , and the minimum between the two profiles dictates the shape of G T (z) (Figure 1

(c) and Section 2.2), with R T (z) ∼ exp[βG T (z)].
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 apK a the only differences between the permeabilities of acidic and neutral compounds as a function of G W→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 T (z) for a polar and an apolar compound with apK a = 5. For G W→M 0 (main plot in Figure 3), combining the neutral and charged PMFs results in a repulsive G T (z) with an effective barrier given by the sum of G W→M and G b→a . As for polar compounds the permeability coefficient is mainly dictated by the height of the PMF, increasing this by a factor G b→a -which is constant at fixed apK a and independent of G W→M , see Figure 1(c) -preserves the linearity of log 10 P in G W→M but generates the offsets observed in Figure 2.
For G W→M 0 (inset in Figure 3), the effect of a small vertical shift of G N (z) for a fixed apK a is to reduce the attractive well that in the neutral case was not contributing to the resistivity. By increasing G W→M , the well crosses zero and G T (z) 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 G W→M 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 Figure 3. Construction of G T (z) (blue) for an acidic compound from its charged G C (z) (red) and neutral G N (z) (orange) PMFs. In calculating G T (z), G N (z) is vertically shifted by G b→a = k B T(pH − pK a ) ln 10, and we consider apK a = 5. We present results for a polar (main plot) and an apolar (inset) dimer. partitioning free energy G W→M , and the acid/base free-energy difference G b→a , which directly links to the apK 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 T (z) 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 C (z) has reached the shifted neutral PMF within the membrane core -a saturation condition always satisfied for apK a 4 -an increase in its amplitude does not affect G T (z). 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.
Permeability coefficients obtained from the two deprotonation prescriptions are presented in Figure 4 (main plot). Results are essentially indistinguishable. Slight deviations are observed for apK a = 4, in which the condition of saturation of G C (z) is not always satisfied as we are moving towards strong acidsvide 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.

Strong acids: the role of amphiphilicity
Panel (a) of Figure 5 displays permeability data of CG dimers as a function of G W→M for apK a = 2, together with some previously analyzed results for apK a 4. At apK a = 2, data initially follow a linear trend analogous to the one of weak acids due to the saturation of G C (z) in the effective PMF. However, for G W→M 0 we observe a spread of the points over the (log 10 P, G W→M ) plane: the permeability coefficient at fixed apK a becomes multivalued in G W→M . This does not occur in the case of CG unimers, whose log 10 P continue to behave linearly as G W→M increases (data not shown).
In strong acids, G b→a is such that neutral PMFs are shifted towards high values with respect to their charged counterparts (Figure 1(c)): only G C (z) contributes to the effective PMF, and R T (z) ∼ exp[βG C (z)]. Different permeabilities for a given G W→M thus arise from compounds presenting very different polarities in their depronotonated forms, despite having similar ones in the neutral case -i.e. similar G W→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 G W→M : while these will lead to roughly the same free-energy difference for the neutral PMFs G N (z), 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 Two-dimensional projection of the permeability data on the (log 10 P, G W→M ) plane in the case of strong acids with apK a = 2. We include as reference a subset of the data obtained for weak acids, apK 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 N (z) and G C (z) 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. asymmetry in polarity between the two beads. Due to additivity, the compounds are characterised by similar neutral PMFs G N (z) (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 T (z) ∼ exp[βG C (z)], 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, G asy = G 2 − G 1 , where beads are arranged in order of decreasing polarity. For a given G W→M , increasing G asy amounts to decreasing/increasing the Figure 6. Permeability coefficients (log 10 scale) for CG dimers in the case of strong acids with apK a = 2. We display results as a function of G W→M and the asymmetry G asy . G 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.
polarity of the more apolar/polar chemical group. As such, G asy is essentially quantifying the amphiphilicity of the small molecule.
Permeability data for apK a = 2 as a function of G W→M and G asy are presented in Figure 6. The introduction of G asy completely solves the degeneracy in G W→M observed in Figure 5(a): for a given G W→M of the neutral form of the acid, increasing G 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 C (z). 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 T (z). 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 apK a . For strong acids, however, our results provide a lower bound for the permeation rate: at a fixed apK a , the permeability coefficient cannot be lower than the linear relation dictated by combining two free-energy differences, G W→M and G b→a .

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, G W→M , and acid dissociation constant, apK 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 apK a and smoothly depends on G W→M . This is due to the fact that the permeability coefficient is mainly dictated by G W→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 G W→M and the acid/base free-energy difference G b→a , the latter being a linear function of the apK 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 G W→M and G b→a .