Data-driven equation for drug-membrane permeability across drugs and membranes

Drug efficacy depends on its capacity to permeate across the cell membrane. We consider the prediction of passive drug-membrane permeability coefficients. Beyond the widely recognized correlation with hydrophobicity, we additionally consider the functional relationship between passive permeation and acidity. To discover easily interpretable equations that explain the data well, we use the recently proposed sure-independence screening and sparsifying operator (SISSO), an artificial-intelligence technique that combines symbolic regression with compressed sensing. Our study is based on a large in silico dataset of 0.4 million small molecules extracted from coarse-grained simulations. We rationalize the equation suggested by SISSO via an analysis of the inhomogeneous solubility-diffusion model in several asymptotic acidity regimes. We further extend our analysis to the dependence on lipid-membrane composition. Lipid-tail unsaturation plays a key role, but surprisingly contributes stepwise rather than proportionally. Our results are in line with previously observed changes in permeability, suggesting the distinction between liquid-disordered (Ld) and liquid-ordered (Lo) permeation. Together, compressed sensing with analytically derived asymptotes establish and validate an accurate, broadly applicable, and interpretable equation for passive permeability across both drug and lipid-tail chemistry.


I. INTRODUCTION
Passive lipid-membrane permeation is and remains of great relevance for pharmaceutical applications and an improved physicochemical understanding of small-sized molecules in complex biological materials. 1,2 The technological implications of the problem has sustained the need for experiment-and simulation-free prediction of passive permeation that are rapid, inexpensive, and accurate. 3,4 Various types of surrogate models have been proposed over the years, the field having adopted machine learning early on. 5 While modern deep-learning approaches take advantage of unchallenged model expressivity to offer unprecedented prediction accuracy, they suffer from two important drawbacks: 1. Overfitting: The size of chemical space of druglike molecules is overwhelmingly large (∼ 10 60 compounds). 6 Deep surrogate models need large numbers of parameters to establish complex relationships. Unfortunately the body of experimentally available data is minuscule compared to the size of chemical space. This can lead to surrogate models that shift dangerously upon addition/removal of small numbers of compounds in the training set. The problem is aggravated by databases that are often proprietary, preventing broad availability and reproducibility. Relying on a) Electronic mail: dutta@mpip-mainz.mpg.de different measurement batches tends to also accentuate systematic errors.
2. Lack of interpretability: Surrogate models are oftentimes black-box techniques that typically cloud why a certain prediction has been made. Deep neural networks exhibit an overwhelming number of parameters and rely on highly non-linear hierarchical functions, making them nearly impossible to conceptually grasp. 7 Quantitative structureactivity relationship (QSAR) models can build multivariate models, but the combination of too many descriptors will lead to similar difficulties. Beyond predicting individual data points, we seek to gain further insight. Insight is essential, for instance as a stepping stone to solving the inverse problem, thereby establishing structure-property relationships and enabling compound design.
In this work, we address these points using a combination of approaches. Critically, we address overfitting by relying on large datasets applied to simple models. Instead of experimental approaches, we base our study on in silico measurements, taking advantage of the rapid rise of high-throughput molecular dynamics simulations. 8 All-atom simulations offer a gold standard in terms of simulation-based permeability modeling that can reach exquisite correlation against experimental reference, but their overwhelming computational costs unfortunately limit studies to tens of compounds. [9][10][11] Here instead, we use an approach based on particle-based coarse-grained (CG) simulations, making use of the transferable Martini model. [12][13][14][15] The Martini model sacrifices some chemical resolution, but retains the essential driving force, mainly the partitioning coefficient of a chemical group between different phases. This allows the CG approach to recover excellent accuracy: 1.4 kcal/mol along a potential of mean force (PMF), translating to 1 log 10 unit in the permeability coefficient, validated across an extensive set of structurally distinct compounds against both atomistic simulations and experimental measurements. 16 Critically, the accuracy of the CG model is accompanied by a three-orders-of-magnitude speedup compared to atomistic simulations. The high-throughput coarsegrained (HTCG) simulations offer unprecedentedly large databases of permeability coefficients: Menichetti et al.
reported results for 511 427 compounds. 16 The ability to screen over so many compounds results mainly from the transferable nature of the coarse-grained modelmany molecules map to the same set of CG beads, which effectively reduces the size of chemical space. 17 The benefits for an efficient computational-screening procedure outweigh the impinged degeneracy, as indicated by the above-mentioned accuracy. Overall, the database contains a nearly exhaustive subset of small organic molecules in the range 30-160 Da, thereby ensuring a dense coverage of the chemical space in this subset. A deeper analysis of this large database is the subject of this study. As for the data-driven model, we explicitly avoid building and using a black-box model, and instead turn to learning an equation. In particular, we will rely on recently proposed data-driven techniques to discover equations. 18,19 Equations relevant to physical problems often display simplifying properties, such as symmetries or separability easing both their data-driven discovery and generalization beyond the training set. Several studies have demonstrated the ability to (re)discover physics equations. 20,21 Generalization capabilities are critical, because typical training datasets are minuscule relative to the size of chemical compound space, such that overfitting can easily prevail. How do we discover simple equations from the data? To this end, we follow Occam's razor and limit the complexity of the equations we consider. The combination of descriptors through various mathematical operators will lead to an overwhelming number of trial equations, many of which may fit the data similarly well. Following ideas from compressed sensing, we select simple equations by applying the 0 regularization. Identifying simple equations both benefits generalization aspects, but also interpretability, i.e., rationalizing the derived equation.
Passive drug permeation measures the propensity of a solute to spontaneously cross a lipid membrane. In this paper, we exclude transporter-mediated uptake to focus solely on thermal diffusion. 22 Upon doing so, the solute interacts with a great variety of physicochemical environments-from an aqueous phase to a hydrophobic membrane core. Permeation is quantified by means of its coefficient, P , as the steady state flux of the solute across the soft interface. Early on, Meyer 23 and Overton 24 mod-eled passive permeation as diffusion across a homogeneous slab via P = KD/σ, where K and D are the water/membrane partitioning coefficient and diffusivity of the compound, respectively, and σ is the thickness of the bilayer core. K is typically approximated by a simpler proxy, namely the partitioning coefficient between water and octanol. Water/octanol partition, or more generally hydrophobicity, has long been identified as strongly correlating with membrane permeability. 25 Notable refinements to the homogeneous Meyer-Overton rule include the inhomogeneous solubility-diffusion model (ISDM), estimating the permeability coefficient via an integral over the membrane extension, z, of its PMF, G(z), as where R(z) and β = 1/k B T correspond to an associated resistivity and the inverse temperature, respectively. 26 The competition between different protonation states naturally follows the sum of inverse resistivities, analogous to the total resistance in a parallel electrical circuit. 10 PMFs are shifted according to the difference between the solution's pH and the compound's acid dissociation constant, pK a . In the following, we take the perspective of a neutral compound, which can deprotonate (acidic, apK a ) or protonate (basic, bpK a ). Knowledge of the PMF(s) and the diffusivity thereby fully determines the permeability coefficient. Unfortunately, these (i) are so far only available by computational techniques and (ii) typically require extensive calculations (∼ 10 5 CPU-hour per system at an atomistic resolution). 9,11 Even at a CG resolution, Eq. 1 still requires free-energy calculations to determine the PMF. The objective of this study is to gain further insight into the key physical determinants of the permeability coefficient. Beyond the widely known effect of hydrophobicity, we focus on incorporating the role of acidity via the compound's protonation state. 27 The role of acidity is crucial, partly because of the sheer number of ionizable drugs: they make ≈ 62% of the World Health Organization's list of essential drugs. 28,29 It has long been hypothesized that the neutral form of an ionizable drug is the only contributor to its permeability, known as the pH-partition hypothesis. 30,31 However, this hypothesis is limiting for two reasons. First, permanently charged compounds are known to permeate lipid membranes. 32,33 Second, differences between a compound's pK a and the surrounding pH can lead to a protonation-coupled permeation mechanism, which calls for the combined contributions of the neutral and ionized forms. 34,35 Clearly the role of acidity is expected to couple with hydrophobicity. Establishing a functional relationship connecting the two quantities is the objective of this work. To this end, we rely on modern data-driven techniques to discover an equation. Relating acidity to the permeability coefficient would not only help establishing rapid estimates for ionic compounds, but also offer insight into the coupling of these physic-ochemical properties that are valid for a wide class of compounds.
Limitations in the number of candidate descriptors and correlations between features have recently been addressed by sure-independence screening and sparsifying operator (SISSO), 36,37 which is an artificialintelligence technique that combines symbolic regression with compressed sensing. 38-40 SISSO provides a datadriven framework to discover equations-mathematical relations between input variables that best correlate with the target property. We will discuss several equations of various complexities to illustrate the balance between accuracy and interpretability. We root the simplest variant in the underlying physics by a comparison against analytical acid-base asymptotic regimes. This simple model incorporating both hydrophobicity and acidity allows us to easily extend our analysis to different lipid membranes starting from limited information. From knowledge of neutral species alone, we predict the change in passive permeability in various lipid membranes. We finally discuss the change in permeability in the context of membrane-phase behavior.

A. Database of drug-membrane thermodynamics
Our analysis is based on the passive-permeability database provided by Menichetti et al. 16 Reference information includes water/octanol partitioning free energy (∆G W→Ol ), acid dissociation constant for acids and bases (apK a and bpK a ), and simulation-based permeability coefficient (expressed by its order of magnitude, log 10 P ) for solutes through a single-component bilayer made of 1,2dioleoyl-sn-glycero-3-phosphocholine (DOPC). The water/octanol partitioning free energies were predicted using the neural network ALOGPS, 41 and acid dissociation constants pK a were predicted from ChemAxon Marvin. 42 Filtering of the Generated DataBase (GDB) 43 for compounds that mapped to a one-or two-bead coarse-grained Martini representation, i.e., a monomer or a dimer, led to 511 427 small organic molecules (30-160 Da). Enhancedsampling molecular dynamics simulations yielded the PMFs for both neutral and (de)protonated species, and Eq. 1 was used to compute the permeability coefficient, P . The diffusivity profile was not extracted from coarsegrained simulations, but instead from previous atomistic studies, taking advantage of its relatively small chemical dependence and its logarithmic impact on the permeability. 10 The pK a of a chemical group can be either acidic (apK a ) or basic (bpK a ), in that a neutral compound can either deprotonate or protonate (see SI for definitions). While the ionization constant of conjugated acid/base pairs typically lie between 10 −2 and 10 16 , we considered compounds with pK a values between −10 and 20. 44 This led to a dataset of 418 828 compounds used as part of this work. A follow up work to Menichetti et al.

B. Learning algorithm
To learn an interpretable model of passive permeability, we used SISSO, as implemented in Ref. 46. SISSO aims at establishing a functional relationship, y = f (Φ), between n primary features, Φ 0 = {φ 1 , φ 2 , · · · , φ n }, and a target property, y, based on N training compounds. SISSO assumes that y can be reliably expressed as a linear combination of non-linear, but closedform, functions of primary features. To construct these non-linear functions, SISSO recursively applies a set of user-defined unary and binary operators (we used on the primary features and builds up sets of candidate features. Φ q denotes the set of candidate features at each level of recursion q. The number of candidate features in Φ q increases sharply with increase in the recursion level q, the number of operators used, and the number of primary features n. For each q, SISSO selects iteratively subsets of candidate features that have the largest linear correlations with the target y and then with the subsequent residuals, i.e., each portion of y that is captured by the previous iterations (see Ouyang et al. 37 ). The number of iterations in this procedure, which equals the number of terms in the linear expansion of f (Φ), and hereby denoted as dimension of the model, is controlled by a sparsifying 0 regularization. For each q and number of dimensions, SISSO selects the model with the smallest root mean-squared error (RMSE). We also quantify model performance using the maximum absolute error (MaxAE) and the square of the Pearson correlation coefficient, r 2 .

C. Feature construction and training
We apply SISSO to three easily accessible and interpretable primary features: the water/octanol partitioning free energy, ∆G W→Ol , and the acid dissociation constants apK a and bpK a as provided by Menichetti et al. 16 We thereby seek a refinement or correction to the commonly used model based on hydrophobicity alone. 25 The mean absolute errors associated with the ∆G W→Ol and pK a predictions (0.36 kcal/mol 41 and 0.86 units, 47 respectively) make them reliable primary descriptors. To avoid constructing features with different units, we multiply the partitioning free energy by the inverse temperature: β∆G W→Ol , using T = 300 K following Menichetti et al. 16 Starting with the set of primary features Φ 0 = {β∆G W→Ol , apK a , bpK a }, we consider the construction of secondary features for up to two iterations (i.e., q = 2), where Φ 1 and Φ 2 consist of roughly 30 and 2 000 features, respectively. We limit the SISSO screening size to 500 and consider up to three-dimensional descriptors. We train on 10% of the available data (see SI for the input script), and use the remaining 90% for holdout evaluation. We draw these train/test sets uniformly at random, without replacement. To reduce variance, we report the average performance over ten independently drawn train/test sets. Ouyang et al. 36 emphasized that SISSO works reliably when the number of materials in a training set, N , sufficiently exceeds the number of candidate features. In particular, SISSO requires N ≥ kd ln(#Φ), where k ∼ 1-10, d is the dimension, and #Φ is the size of the feature space. For Φ 2 , the relevant feature space used to train our model, the relation requires N ≥ 10·2·ln(2·10 3 ) 150. By training our models with only 10% of the dataset (∼ 42 000 compounds), SISSO is well within a regime to provide meaningful and consistent results.

III. RESULTS AND DISCUSSIONS
A. Learning permeability models Table I contains the four models considered in this work: (i) f Hyd is a baseline hydrophobicity model, which linearly correlates with water/octanol partitioning free energy; and (ii-iv) f 1D to f 3D linearly correlate with one to three secondary feature(s) identified by SISSO. For each model, c m i correspond to non-zero coefficients for model m and index i, reported in Tab. I. For all models, ∆G W→Ol takes on a central role, as expected by the performance of the baseline. The simplest SISSO , is remarkably robust: it is systematically identified as the best performing 1D model across all 10 training sets. Given that we trained on only 10% (∼ 42 000 compounds) of the dataset, this highlights this model's performance compared to all other candidates (see SI for a list of candidate one-dimensional models). The stability of the model-given such a small training fraction-speaks for the robustness of the equation. Its improved accuracy compared to the baseline will be evaluated further down.
We also report more complex 2D and 3D models in Tab. I. While we will show below that these yield even better accuracy compared to f 1D , they are specifically tailored to the training set used: f 2D and f 3D are ranked as the best model in eight and five out of the 10 training sets, meaning that other models of similar complexity closely compete.

B. Model performance
We now turn to the performance of these four models. Tab. I reports their RMSE, MaxAE, and squared Pearson correlation coefficient, r 2 , averaged over the test sets. Going from the baseline to more complex SISSO models, the systematic decrease in the RMSE is accompanied by an increase in the correlation coefficient. On the other hand, the maximum absolute error does show a clear minimum for f 1D . This offers a first hint at the appealing balance between generalization and interpretability of the 1D SISSO model. The performance of these four models is depicted in Fig. 1 for the entire dataset, where we report each model against reference values. Going from baseline to SISSO models of increasing complexity, the distribution does lean increasingly toward the y = x correlation line. The presence of horizontal stripes in Fig. 1 results from the degenerate use of CG mappings for many molecules. 17 This artifact is most notable for f Hyd , which solely relies on hydrophobicity, whereas the others have chemically specific acidity information. For the 2D and 3D models, we also point out outliers at the lowest permeability values. Fig. S1 in the SI shows the distribution of compounds: these low-permeability values are scarcely populated, both in algorithmically and synthesized compound databases. 16 Here they represent only 0.07% of the dataset, and our uniform sampling of training points likely brought in only few of them. They likely result from poor extrapolation behavior of the 2D and 3D models, which notably include powers of two of several variables. To better understand the performance of each model,  we analyze their errors in greater detail. Fig. 2a shows the distribution of absolute error. The large dataset at our disposal allows us to evaluate more than 4 orders of magnitude of this distribution. The comparison between f Hyd and f 1D proves insightful: while they are remarkably close up to errors of 5 log 10 units, the baseline then displays a significant hump, while the SISSO 1D model keeps decaying monotonously. Both models rely on β∆G W→Ol , which explains the remarkable agreement early on, while the stark difference between the two curves is entirely due to the effect of acidity. This is confirmed by a further decomposition of the error as a function of acidity, showing that f Hyd leads to larger errors for stronger acids and bases (Panels b and f), while SISSO 1D significantly reduces the error in this regime (Panels c and g).
In comparison, the more complex SISSO 2D and 3D display a shift in the error distribution toward lower errors (Fig. 2a) compared to the 1D model. At low probability however, we observe a significant change in the slope of the decay, indicating worse performance for a small number of outliers. This is also illustrated when decomposing the error in terms of acidity in Panels (d-i): while the overall performance improves, we identify more outliers. These outliers mostly lie at low permeability values (reminiscent of Fig. 1), and for strong apK a or bpK a . Poor performance at large acidity values could take place if these were absent of the small training fraction.

C. Validation against atomistic simulations
The SISSO models should naturally be prone to systematic errors inherent to the training dataset. While we expect our systematic integration of the ISDM permeability coefficient (Eq. 1) to ensure robust functional relationships, systematic errors in the parameters are likely to affect the fitting coefficients. Reference permeability values were extracted from computationally efficient coarse-grained computer simulations, at the cost of forcefield accuracy. Still, a comparison of the coarse-grained simulations against atomistic computer simulations had indicated an excellent agreement of 1 log 10 unit across a limited set of small molecules 16 . Here we compare the performance of the four permeability models against the atomistic simulations of Carpenter et al. 10 This set of 12 organic compounds covers a range of molecular weights that goes significantly beyond our training set: an average of 243 Da and up to 319 Da, comparable to that of real drugs, given that more than 60% of drugs have molecular weight below 300 Da. 48 On the other hand our training HTCG database only contained compounds up to 160 Da. This thus presents a challenging test for the generalizability of the SISSO models. Fig. 3 shows the absolute error against atomistic simulations across all 12 small molecules and for our four models. For each model, we display the error as a function of both apK a and bpK a . The baseline model yields absolute errors between 1 and almost 4 log 10 units. While larger errors correlate with strong acids, they do not seem to correlate with larger bases. Unlike Fig. 2, the minuscule set of atomistic compounds prevents us from drawing conclusions that would hold across any significant subset of chemical space. Turning to SISSO 1D, we observe a small but noticeable decrease in performance, where the mean absolute error (MAE) increases from 1.55 to 2.03 log 10 units. The MAE, however, decreases against the baseline when considering the more complex SISSO 2D and 3D: 0.87 and 1.33 log 10 units, respectively. The nonmonotonic decrease of the MAE with increasing complexity in Fig. 3b-d suggests the role of the small validation dataset considered. Overall though, the incorporation of acidity does lead to an improved reproduction of the permeability coefficient. It validates the SISSO-derived equations on permeability coefficients derived using independent methods and outside the chemical space of the training data.

D. Acid-base asymptotes
The analysis so far highlights how model complexity impacts accuracy. Missing from the analysis so far is the consideration of interpretability. The two one-dimensional models-the baseline and SISSO 1Dhighlight a simple mechanism as to the functional dependence of the permeability coefficient on both hydrophobicity and acidity. Focusing on SISSO 1D specifically, we rewrite the model in terms of two contributions Fig. 4 displays the permeability coefficient as a function of these two contributions. The symmetric contribution of β∆G W→Ol in the two terms of Eq. 2 indicates that the baseline hydrophobicity model manifests itself along the diagonal. Notably missing from the diagonal behavior are the dark vertical and horizontal basins. They localize at apK a − β∆G W→Ol ∼ 0 and −bpK a − β∆G W→Ol ∼ −15, and represent strong-acid and strong-base regimes. In what follows, we provide an asymptotic rationalization of the functional form of Eq. 2.
To rationalize Eq. 2, we first outline the role of our three primary descriptors in the ISDM model (Eq. 1). Fig. 5 illustrates the well known interplay between PMF and acidity, in particular how the latter shifts the PMFs of the neutral and (de)protonated species. In the following, we will denote the PMFs of the neutral, protonated, and deprotonated species as G N (z), G B (z), and G A (z), respectively.
The difference between apK a or bpK a and pH dictate the propensity for the PMFs to cross each other. The ISDM relies on a total resistivity (defined in Eq. 1), R T , such that R −1 , analogous to the total resistance in a parallel electric circuit.
The PMFs of the neutral, protonated, and deprotonated species can be linked in water thanks to their apK a and bpK a values, as well as the pH of the environment, through the following equations where G(∞) indicates that the compound is located in bulk water. Eqs. 3 and 4 effectively link the difference between the pH of the environment with the pK a of the compound to a shift in the PMFs. Without loss of generality, we will shift all free energies such that zero corresponds to the most favorable compound in bulk water. Eqs. 3 and 4 allow us to explicitly link pK a information with the total resistivity where we assume that all protonation states yield identical diffusivity. 16 Because Eq. 5 takes on a relatively complex form, we will consider only asymptotic regimes: • Neutral compounds entail no strong acid or base characteristic, i.e., apK a pH and bpK a pH, such that the compound is effectively unable to (de)protonate, and G N (∞) = 0. Eq. 5 can be sim- • Strong acids consist of solutes that display at least one functional group for which apK a pH. For neutral pH (≈ 7), this would correspond to apK a < 4. We set G A (z → ∞) = 0. In this regime the third exponential in Eq. 5 would dominate the other two, leading to R T (z) ≈ D −1 (z) exp [βG N (z) + (pH − apK a ) ln 10].
• Strong bases would display at least one functional group where bpK a pH. For neutral pH, this would correspond to bpK a > 10. We set G B (z → ∞) = 0.
Using a similar argument, Eq. 5 would be dominated by  The total resistivities still require integration over the reaction coordinate z, which we simplify to the largest contribution of the PMF. 11 The effective resistivity model is equivalent to choosing the lowest PMF at any value of z: G eff (z) = min i G i (z), where i runs over the neutral, protonated, and deprotonated species. In addition, the dominating contribution of the effective permeability will come from its maximum value, corresponding to a posi- tion z * = arg max z G eff (z). Assuming that the largest contribution of the permeability arises from the total re-sistivity at z * , we obtain P ≈ R −1 T (z * ). This yields the following acid-base asymptotic regimes if apK a pH and bpK a pH.
To numerically test Eq. 6, we identify datapoints corresponding to the three asymptotic regimes: the neutral compounds (apK a > 10 and bpK a < 4), strong acids (0 < apK a < 4 and bpK a < 4), and strong bases (10 < bpK a < 14 and apK a > 10). For simplicity, we only considered non-zwitterionic compounds. We assume that log 10 D(z) yields no significant, chemically specific effect, and uniformly shift the permeability coefficient across the chemical space of compounds considered. Fig. 6 shows the agreement between Eq. 6 and the reference permeability coefficients. All show strong linear correlation: for neutral compounds (r 2 = 0.998), strong acids (r 2 = 0.959), and strong bases (r 2 = 0.986). These results numerically validate the asymptotes of Eq. 6.
More importantly, the asymptotes provide a physically motivated rationale for the two contributions of Eq. 2: apK a − β∆G W→Ol and −bpK a − β∆G W→Ol . We first note that ∆G W→Ol is related to G N (z * ). The depth at which the effective PMF is the highest, z * , will almost always be close to the membrane midplane: z * ≈ 0. The main exception to this are hydrophobic compounds, for which the highest point in the PMF is in water (Fig. 5c). Furthermore, G(z * = 0), which corresponds to the transfer free energy from water to membrane midplane, has been shown to linearly relate to the water/octanol partitioning free energy, ∆G W→M ∝ ∆G W→Ol . 50 This connects G N (z * ) present in Eq. 6 to ∆G W→Ol in Eq. 2. We then observe that the asymptotes and Eq. 2 have the same signs and exponents of apK a and bpK a . For an acidic or a basic compound, if we consider the compoundspecific pK a and substitute ∆G W→Ol by the compound's G N (z * ) in the relevant among the two contributions of Eq. 2, we indeed recover one the asymptotes. As for neutral compounds, G N (z * ) is the only relevant quantity while estimating permeability.

E. Drug-membrane permeability across membranes
The following analyzes how drug-membrane permeability changes according to membrane composition. We hypothesize that the functional form of SISSO 1D is applicable to other lipid membranes, and use it as a starting point. We take advantage of the above-mentioned asymptotic regimes to limit the amount of information needed from new membranes. The regime of neutral com-pounds described in Eq. 6 can be used advantageously because it only requires information on neutral PMFs. We rely on the dataset of Hoffmann et al., which precisely contains PMF information-but no permeabilityin various membranes, and only for neutral compounds. 45 We specifically analyze the change in permeability when turning to phosphocholine (PC) membranes made of different lipids, varying in both tail length and level of unsaturation. Fig. 7 shows the relation of −βG N (z * )-the dominant term for the permeability of neutral compounds (Eq. 6)between the original membrane used in this work, DOPC, and others. All curves follow a line, indicating that the asymptotic regime for neutral compounds of Eq. 6 holds for all membranes. We find two families of lines with different intercepts: DLPC, DPPC, and POPC show an intercept with DOPC that is roughly 0, while DIPC and DAPC have an intercept that is approximately 1.4.
To better understand these results, we first recall that G N (z * ) corresponds to the highest value of the neutral PMF. We can safely ignore contributions of the charged PMFs, such that G N (z * ) denotes the highest point of the effective PMF. The excellent agreement between DLPC and DPPC indicates that tail length (3 and 4 beads long, respectively) does not impact the permeability. This is expected, given that tail length is only expected to change the length, but not the height, of the hydrophobic plateau. On the other hand, the further agreement between them and POPC and DOPC indicate a lack of dependence on tail saturation for these lipids (1 and 2 unsaturated beads, respectively). Remarkably, the shift in the intercept appears only for DIPC and DAPC-lipids that exhibit more unsaturations: 4 and 8, respectively. Notably, they display the same shift in −βG N (z * ). As such, the level of unsaturation is a determining factor, but does not gradually impact −βG N (z * ).
Interestingly, previous Martini studies of ternary membranes with DPPC and cholesterol have shown that DIPC and DAPC strongly phase separate into a liquiddisordered (Ld) phase. 51-53 Inconsistent conclusions were drawn from different studies with DOPC, pointing to a thermodynamic drive that is weak, at best. 52,54,55 As for POPC, there is no sign of phase transition. 52 The results on hinges on the presence of DPPC and cholesterol, which are notably absent from our reference simulations. 45 The trends are surprising in that they show a dependence on lipid-tail unsaturation that is stepwise rather than proportional. While we defer a more detailed study to future work, we suggest the role of the entropic character of the lipid-tail fluctuations.
Other studies before us have reported a clear change in permeability between Ld and Lo domains: Ghysels et al. used both atomistic simulations and electron paramagnetic resonance spectroscopy experiments on the permeation of water and oxygen, and found a permeability ratio P (Ld)/P (Lo) ≈ 3. 56 Here the CG simulations yield a shift in −βG N (z * ), which translates to a ratio of the permeability coefficients of ≈ 25. Our estimates are thus within one log 10 unit of the results of Ghysels et al. for their specific compounds. The mechanism remains to be clearly identified, although this could be consistent with the proposed role of local membrane surface density (i.e., its propensity to form transient holes). 57 To summarize, our use of single-component lipid membranes only allows us to speculate as to the shift in Fig. 7. The link to Lo/Ld domain formation is in line with prior atomistic simulations and experiments for specific compounds. 56 In the broader context, our results could help generalize their results across the chemical space of drugs. The results also suggest simple additive corrections to our effective equation when considering different membranes.

IV. CONCLUSIONS
We propose to learn a functional relationship between hydrophobicity and acidity as a simple surrogate for passive membrane permeability. Our approach, combining symbolic regression and compressed sensing, is data-driven and interpretable, and based on large databases of high-throughput coarse-grained simulations. Sure-Independence Screening and Sparsifying Operator (SISSO) builds a hierarchy of models of increasing complexity. Models prove increasingly accurate, yet more complex models are more prone to discrepancies for a few outliers. Our SISSO 1D model offers improved accuracy compared to the hydrophobicity baseline, and yet excellent interpretability. We identify the simple and interpretable equation f 1D = c 1D 0 + c 1D 1 (apK a − bpK a − 2β∆G W→Ol ), where apK a and bpK a characterize acidity, ∆G W→Ol is the water/octanol partitioning coefficient, and c 1D 0 and c 1D 1 are the only two fitting parameters. We rationalize the model by an analysis of the asymptotic regimes of the inhomogeneous solubility-diffusion model (ISDM). The asymptotes are validated numerically and confirm the SISSO 1D equation, implicitly testifying to the accuracy of the underlying HTCG resolution. Broad agreement numerically validates the use of a single bulk hydrophobicity measure to effectively replace the potential of mean force, which has been exploited by others before. 11 Importantly, the interplay of hydrophobicity together with acidity leads to a significant improvement in model accuracy of the SISSO 1D equation for ionizable groups. The SISSO equations show improvements over a challenging set of compounds that are much larger than the training set. Critically, our work refines the common role of hydrophobicity in passive permeation to relate it functionally with acidity.
The separation of the ISDM in asymptotic regimes allows us to build drug-permeability models across membranes with limited information only. Using only the potential of mean force of neutral solutes, we infer the change in permeability for membranes with lipids of varying tail length and level of unsaturation. We observe a surprising change in permeability coefficient: lipid-tail unsaturation contributes stepwise rather than proportionally. Our findings are in line with recent atomistic simulations and electron paramagnetic resonance spectroscopy experiments, highlighting the distinction between lipids primarily involved in liquid-disordered (Ld) and liquid-ordered (Lo) domains. The approach offers a data-driven, interpretable analysis of drug-membrane passive permeability across both drugs and membranes.

SUPPLEMENTARY MATERIAL
See supplementary material for definitions of dissociation constants (apK a and bpK a ), distribution of compounds across permeability, Tab. I along with the error values, list of best one-dimensional descriptors, and the SISSO input script.

DATA AVAILABILITY
In this study, we used the database provided by Menichetti et al. 16 It is openly available at https://doi. org/10.1021/acscentsci.8b00718.
To comply with the unified definition of pK a of ChemAxon-it is the ratio of conjugate acid to conjugate base-the corresponding basic pK a , denoted as bpK a , is defined as bpK a = − log 10 K a = pH + log 10 The usefulness of this definition is that now both apK a and bpK a are written as pK a = pH + log 10 [conjugate acid] [conjugate base] . Strong acids, as defined in the Acid-base asymptotes section of the paper, have low (apK a ≤ 4). At pH = 7, from Eq. 4 we find that [AH] [A − ] = 10 apKa−7 ≤ 10 −3 .     S2: Best one-dimensional descriptors for ten training sets as predicted by SISSO. Each column corresponds to a particular training set. The data demonstrates the robustness of the predictions across training sets-only twelve unique descriptors are present the best ten descriptors from all training sets. The top three descriptors do not change. The best one-dimensional descriptor ((apK a − β∆G W→Ol ) − (bpK a + β∆G W→Ol )) and the baseline hydrophobicity descriptor have been highlighted for reference.