Quantifying Small-Molecule Association with Lipid Membranes: Methods, Models, and Limitations

1. Introduction

The association of drugs with biological membranes is one of the key steps that determines drug bioavailability and pharmacokinetics. The importance of this property is well recognized and is usually estimated from drug's partition between water and octanol due to its simpler experimental evaluation. However, octanol is not able to capture the distinctive features of biological membranes. On one hand, the polarity of octanol is much higher than that of the non-polar bilayer center and much lower than that of the polar headgroups, lacking the charged groups present in phospholipids, the most abundant lipids in biomembranes. Octanol is therefore unable to mimic the membrane-interaction of ionized drugs, nor any process that is influenced by membrane electrostatic potentials [1,2]. Moreover, while octanol polarity is comparable to the average polarity of lipid bilayers, and similar to the region of intermediate polarity where the energy minimum is located for most drugs, [3,4,5,6,7,8,9,10,11] the H bonding ability is different. Together with the absence of charges, this leads to different relative stabilities of the protonated/unprotonated species of drugs with ionizable groups, which strongly influences their partition and permeability.

Other water-immiscible organic solvents have been attempted to mimic the membrane-affinity of small uncharged solutes [12]. However, no organic solvent or their combinations may capture the solvent properties of lipid bilayers due to the distinct structural features. Membrane’s lipid bilayer is a structured, very thin and strongly anisotropic media, changing from a highly polar and charged headgroup region in contact with water to the non-polar center of the bilayer in just a few nanometers [13,14]. Because the thickness of each monolayer is comparable to the dimensions of most bioactive molecules, they span membrane regions of very distinct solvation properties, adopting conformations, locations and orientations that maximize favorable interactions with the membrane [5,8,10,15,16,17,18,19,20,21,22,23]. As well recognized for many decades (e.g. [24,25]), lipid bilayers are interfacial phases and membrane-partition is of a fundamentally different nature than partition into bulk phases. In this context it is worth noting that octanol is not truly a homogeneous medium at the nanoscopic scale, being organized as inverted micelles with H-bonds between the hydroxyl groups stabilizing the micelle core [26,27,28]. The size of the micelle polar core increases in water-saturated octanol which at room temperature contains over 2 M water [29,30]. The high water content and local structure increases the solubility of polar molecules [30], and is likely the reason for its moderate success in predicting associations with proteins and lipid membranes [31], although the lack of charged groups and distinct H-bonding ability limits its mimicking capacity.

The inadequacy of non-polar solvents as models of biomembranes is one of the reasons for the poor predictivity of drugs pharmacokinetics, namely their ability to cross tight endothelia such as the blood-brain-barrier [32,33,34,35,36,37,38,39]. Membrane partition is also a major player in the pharmacodynamics of many drugs, with their membrane-affinity as well as their location and orientation in the membrane strongly influencing their interaction with membrane proteins.

Over the last century, a large number of studies have been performed regarding drug-membrane interactions. It is of the highest importance to take advantage of this data to improve the predictive ability of membrane affinity from drug's structure and molecular properties. What is delaying this important task is the lack of uniformization of the formalisms used to quantitatively describe the drug/membrane interaction, [40,41,42] and their organization in a publicly available database, containing the relevant information, manually curated, reliable and ready to use.

There are several publications where some specific methodologies to study drug/membrane association have been revised. Among them it is worth noting the work by Santos et al [42] which gives an extensive description of the optical spectroscopic methods, with important information regarding experimental details and possible artefacts. Other important reviews can also be found for the application of optical spectroscopic methodologies, [43,44,45], calorimetry [46,47], potentiometry, dialysis and NMR [48], and for the use of molecular dynamics simulations [49,50]. Each specific approach has its own particularities, namely regarding the membrane properties and the concentration of drug, but also on the parameter to quantify the drug-membrane affinity. To allow the use of the information collectively obtained by the scientific community, it is necessary to establish consensus parameters and conditions, as well as quantitative relationships to convert the information directly obtained into the consensus parameters. This will maximize the information available and allow its use for artificial intelligence applications, to find rules and quantitative relationships to predict the membrane affinity from the drug’s molecular properties.

The objective of this review is to join the information from all available methodologies to characterize the association of small molecules with lipid membranes. Each methodology is documented with some applications and their strengths and weakness are summarized. We first discuss the formalisms used to quantify the membrane affinity, establishing the quantitative relations between them and their application limits. This allows the comparison between results obtained with different approaches and guides new researchers in the field in their selection for the most appropriate approach. Key references are provided for the distinct methodologies and formalisms. The list of references, however, is not exhaustive, and the authors apologize in advance for the omission of many important works. This review will set the framework for the implementation of all available data on the interaction of small molecules with lipid membranes in its future inclusion in a publicly available database.

2. Parameters to Quantify the Affinity for the Membrane

Several parameters can be used to quantify the affinity of a small molecule for lipid membranes. To facilitate comparison between this measure of molecular lipophilicity and the commonly used water/octanol partition coefficient, the reference parameter adopted in this review is the ratio of the solute molar concentrations in the two phases (aqueous and membrane), as defined in Equation (1):

$$\begin{array}{c}{n}_{{\mathrm{S}}_{\mathrm{W}}}\begin{array}{c}{K}_{\mathrm{P}}{V}_{\mathrm{M}}\\ \leftrightarrow \end{array} n_{{\mathrm{S}}_{\mathrm{M}}}\\ \\ K_{\mathrm{P}}={\displaystyle \frac{n_{{\mathrm{S}}_{\mathrm{M}}}/V_{\mathrm{M}}}{n_{{\mathrm{S}}_{\mathrm{W}}}/V_{\mathrm{W}}}}\end{array}$$

(1)

In this equation, $n_{{\mathrm{S}}_{\mathrm{M}}}$is the moles of solute in the membrane phase which has a total volume $V_{\mathrm{M}}$, and $n_{{\mathrm{S}}_{\mathrm{W}}}$is the moles of solute in the aqueous phase with the total volume $V_{\mathrm{W}}$.

Other approaches have been considered in literature for the quantification of the affinity of the solute for the membrane. The most common include the partition in terms of solute concentration expressed in molality ($K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}$, equation (2)), molar fraction ($K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$, equation (3)), or mass ratio ($K_{\mathrm{P}}^{\mathrm{m}/\mathrm{m}}$, equation (4)) of solute in each phase. The distinct parameters are usually reported as partition coefficients without explicit indication in the parameter name of the formalism followed. Care must therefore be taken in the identification of the parameter being reported through the analysis of the equation used for its estimation. This is especially important given that all these different approaches lead to dimensionless affinity parameters, requiring careful reading of the original report for correct identification.

$$K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}={\displaystyle \frac{n_{{\mathrm{S}}_{\mathrm{M}}}/\left(m_{\mathrm{L}}+m_{{\mathrm{S}}_{\mathrm{M}}}\right)}{n_{{\mathrm{S}}_{\mathrm{W}}}/\left(m_{\mathrm{W}}+m_{{\mathrm{S}}_{\mathrm{W}}}\right)}}$$

(2)

$$K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}={\displaystyle \frac{n_{{\mathrm{S}}_{\mathrm{M}}}/\left(n_{\mathrm{L}}+n_{{\mathrm{S}}_{\mathrm{M}}}\right)}{n_{{\mathrm{S}}_{\mathrm{W}}}/\left(n_{\mathrm{W}}+n_{{\mathrm{S}}_{\mathrm{W}}}\right)}}$$

(3)

$$K_{\mathrm{P}}^{\mathrm{m}/\mathrm{m}}={\displaystyle \frac{m_{{\mathrm{S}}_{\mathrm{M}}}/\left(m_{{\mathrm{L}}_{\mathrm{M}}}+m_{{\mathrm{S}}_{\mathrm{M}}}\right)}{m_{{\mathrm{S}}_{\mathrm{W}}}/\left(m_{\mathrm{W}}+m_{{\mathrm{S}}_{\mathrm{W}}}\right)}}$$

(4)

where $m_{\mathrm{L}}$, $m_{\mathrm{W}}$, $m_{{\mathrm{S}}_{\mathrm{M}}}$ and $m_{{\mathrm{S}}_{\mathrm{W}}}$ are the mass of lipid in the membrane phase, water in the aqueous phase, solute in the membrane phase, and solute in the aqueous phase, respectively, and $n_{\mathrm{L}}$ and $n_{\mathrm{W}}$ are the number of moles of lipid in the membrane and the number of moles of water in the aqueous phase, respectively.

If both the aqueous and membrane phases are away from saturation with solute, the different parameters are interconvertible. In this case, the solute concentration is usually much lower than that of the lipid, and the total volume (or mass) of the membrane phase may be represented by that from the lipids only. The different parameters are in this case related by:

$$K_{\mathrm{P}}=K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}{\displaystyle \frac{{\rho }_{\mathrm{M}}}{{\rho }_{\mathrm{W}}}}\mathrm{ }=K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}{\displaystyle \frac{1}{\overline{V_{\mathrm{L}}}\left[W\right]}}$$

(5)

where ${\rho }_{\mathrm{M}}$ and ${\rho }_{\mathrm{W}}$ are the densities of the membrane and aqueous phase respectively. The transversal density of the lipid bilayer is very heterogeneous, being higher in the region of the polar head group, and smaller in the midplane of both membrane leaflets [51,52]. The overall density is however close to 1 g/cm³ [51,52,53,54], and therefore$K_{\mathrm{P}}\cong K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}$. For a given affinity parameter, the fraction of solute associated with the membranes (in conditions away from saturation of both the aqueous phase and the membrane) may be calculated using equations (6):

$${\displaystyle \frac{n_{{\mathrm{S}}_{\mathrm{M}}}}{n_{\mathrm{T}}}}={\displaystyle \frac{K_{\mathrm{P}}\overline{V_{\mathrm{L}}}{\left[L\right]}_{\mathrm{Molar}}}{1+K_{\mathrm{P}}\overline{V_{\mathrm{L}}}{\left[L\right]}_{\mathrm{Molar}}}}={\displaystyle \frac{K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}{\left[L\right]}_{\mathrm{molal}}}{1+K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}{\left[L\right]}_{\mathrm{molal}}}}={\displaystyle \frac{K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}{\left[L\right]}_{\mathrm{Molar}}}{\left[W\right]+K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}{\left[L\right]}_{\mathrm{Molar}}}}$$

(6)

where ${\left[L\right]}_x$ is the concentration of lipid relative to the total solution (the subscript $x$ denoting the concentration units), $\left[W\right]$is the concentration of water in the aqueous phase (approximately 55.5 M), and $\overline{V_{\mathrm{L}}}$ is the molar volume of the lipids in the membrane. The latter is the key property that allows the conversion of the affinity parameters into the partition coefficient$K_{\mathrm{P}}$. This property depends on the lipid composition, namely on the lipid [58–62 molar mass and the membrane phase. For a certain lipid organized as a lipid bilayer, the molar volume is larger when the membrane is in the liquid-crystalline phase, and smaller when in the gel or liquid ordered phases. The specific volume of phosphatidylcholines in the gel phase is around 0.95 mL/g, increasing to around 1.0 mL/g in the fluid phase, [55,56] corresponding to a molar volume of 0.70 and 0.74 dm³/mol for DPPC in the gel and fluid phases respectively. Increasing the length of the acyl chains adds 27 ų per methylene group, [57] which corresponds to a rise of 0.02 dm³/mol, while the introduction of cis double bonds results in only a slight increase in the molar volume (smaller than 0.01 dm³/mol per double bond [56]). The most abundant acyl chains in the phospholipids present in biomembranes contain 16 to 20 carbons, [58,59,60,61,62] leading to an average value of $\overline{V_{\mathrm{L}}}$ close to 0.8 dm³/mol. If the membranes contain high amounts of sterols, the molar volume is significantly lower due to the smaller size of the sterols and to the condensing effect induced by cholesterol in the lipid acyl chains. As observed for the effect of phase transition on the lipid molar volume, the effect of cholesterol condensation itself is relatively small (a decrease of ca. 2% for saturated acyl chains and even smaller for lipids with cis double bonds [56]). Therefore, the effects on the average lipid molar volumes are mostly due to changes in the volumes of each lipid and may be considered additive. The molar volume estimated for cholesterol in a lipid membrane is 0.38 dm³/mol [63]. At 20 mol% cholesterol in the membrane, the average lipid molar volume is therefore decreased by 10 %, and at cholesterol saturating concentrations (ca. 50 mol % [64,65]) it is decreased by 26 %, leading to $\overline{V_{\mathrm{L}}}$= 0.6 dm³/mol. For simplicity, in the conversion of partition coefficients from molar fraction ($K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$) into the consensus partition coefficient ($K_{\mathrm{P}}$), $\overline{V_{\mathrm{L}}}$= 0.8 dm³/mol is suggested for lipid membranes in the gel or in the liquid disordered phase (less than 20 mol% cholesterol) and $\overline{V_{\mathrm{L}}}$= 0.6 dm³/mol for lipid membranes in the liquid ordered phase (more than 30 mol% cholesterol).

The treatment of the lipid membrane as a phase immiscible with water leads to the above formalisms, where the association of the solute with the lipid membranes is treated as a partition coefficient. This is generally valid for the case of small molecules at low concentrations. However, for some solutes and methodologies the association is best described by a formalism considering solute binding to the membranes instead of partition between the two phases. This is particularly relevant when the solutes are large and present several charged groups, as is the case of polycationic peptides and proteins [66]. In this situation, it may be relevant to consider the presence of binding sites in the membrane, as well as the decrease in the binding capacity due to solute binding. The kinetic scheme and the equilibrium binding constant for solute binding to independent binding sites in the membrane may be defined by equation (7),

$$\begin{array}{c}{S}_{\mathrm{W}}+B \begin{array}{c}{K}_{\mathrm{b}}\\ \leftrightarrow \end{array} S_{\mathrm{M}}\\ \\ K_{\mathrm{b}}={\displaystyle \frac{\left[S_{\mathrm{M}}\right]}{\left[S_{\mathrm{W}}\right]\left[B\right]}}\hspace{0.33em};\hspace{1em}\left[S_{\mathrm{W}}\right]+\left[S_{\mathrm{M}}\right]=\left[S_{\mathrm{T}}\right]\hspace{0.33em};\hspace{1em}\left[B\right]={\displaystyle \frac{\left[L_{\mathrm{T}}\right]}{{\#}_{\mathrm{L}}}}-\left[S_{\mathrm{M}}\right]\end{array}$$

(7)

where $K_{\mathrm{b}}$ is the binding constant, $\left[S_{\mathrm{W}}\right]$and $\left[S_{\mathrm{M}}\right]$ correspond to the solute concentration free in the aqueous phase and associated with the lipid membrane, respectively, $\left[B\right]$ is the concentration of free binding sites (all with respect to the total volume of the solution), $\left[L_{\mathrm{T}}\right]$ is the total concentration of lipid, and ${\#}_{\mathrm{L}}$ is the number of lipid molecules per binding site. Equations (7) lead to a quadratic equation, where the physically meaningful solution is always $x_{-}$, equation (8).

$$\left[{\mathrm{S}}_{\mathrm{M}}\right]={\displaystyle \frac{-b-\sqrt{b^2-4ac}}{2a}};a=K_{\mathrm{b}}; b=-1-K_{\mathrm{b}}\left({\displaystyle \frac{\left[L_{\mathrm{T}}\right]}{{\#}_{\mathrm{L}}}}+\left[S_{\mathrm{T}}\right]\right);c=K_{\mathrm{b}}\left[S_{\mathrm{T}}\right]{\displaystyle \frac{\left[L_{\mathrm{T}}\right]}{{\#}_{\mathrm{L}}}}$$

(8)

In the case of low affinity binding where only a negligible fraction of solute is associated with the membrane, and when the affinity is evaluated from the dependence of $S_{\mathrm{M}}$ with $S_{\mathrm{T}}$, $K_{\mathrm{b}}$ may be estimated from the simplified equation (9):

$$K_{\mathrm{b}}\cong {\displaystyle \frac{\left[S_{\mathrm{M}}\right]}{\left[S_{\mathrm{T}}\right]\left(\frac{\left[L_{\mathrm{T}}\right]}{{\#}_{\mathrm{L}}} - \left[S_{\mathrm{M}}\right]\right)}}; \left[{\mathrm{S}}_{\mathrm{M}}\right]\cong {\displaystyle \frac{K_{\mathrm{b}}\left[S_{\mathrm{T}}\right]\frac{\left[L_{\mathrm{T}}\right]}{{\#}_{\mathrm{L}}}}{1+K_{\mathrm{b}}\left[S_{\mathrm{T}}\right]}}$$

(9)

The binding constant is related to the partition coefficient by equation (10) [67].

$$K_{\mathrm{P}}={\displaystyle \frac{K_{\mathrm{b}}}{{\#}_{\mathrm{L}}\overline{V_{\mathrm{L}}}}}\hspace{0.33em}$$

(10)

In most publications following this approach, the number of lipids per binding site is not explicitly assumed. This is equivalent to the assumption of ${\#}_{\mathrm{L}}$=1, and the resulting binding constant is usually referred to as $K_{\mathrm{L}}$.

Another association parameter that has been used to quantify the affinity of solutes for lipid membranes results from considering the whole lipid vesicle as the binding agent. This formalism is used when the equilibrium affinity is obtained from the ratio of the association and dissociation rate constants, calculated from the dependence of the rate of association with the concentration of lipid vesicles, and is relevant to evaluate if solute-membrane association is a diffusion-limited process [68,69,70,71,72]. In this case, the equilibrium constant is given by equations (11) and related with the corresponding partition coefficient by equation (12).

$$\begin{array}{c}{S}_{{\mathrm{W}}^{\mathrm{o}}}+L_{\mathrm{V}} \begin{array}{c}{k}_{+}^{L_V} \left[L_V\right]\\ \rightleftharpoons \\ k_{-}\end{array}{ S}_{{\mathrm{L}}^{\mathrm{o}}}\\ K_{{\mathrm{L}}_{\mathrm{V}}}={\displaystyle \frac{k_{+}^{{\mathrm{L}}_{\mathrm{V}}}}{k_{-}}}={\displaystyle \frac{\left[S_{{\mathrm{L}}^{\mathrm{o}}}\right]}{\left[S_{{\mathrm{W}}^{\mathrm{o}}}\right]\left[L_{\mathrm{V}}\right]}} ;\hspace{1em}\left[L_{\mathrm{V}}\right]={\displaystyle \frac{\left[L_{\mathrm{T}}\right]}{n_{\mathrm{L}/{\mathrm{L}}_{\mathrm{V}}}}}\end{array}$$

(11)

$$K_{\mathrm{P}}={\displaystyle \frac{K_{{\mathrm{L}}_{\mathrm{V}}}}{{\overline{V}}_{{\mathrm{L}}_{\mathrm{V}}^{\mathrm{o}}}}}={\displaystyle \frac{K_{{\mathrm{L}}_{\mathrm{V}}}}{n_{\mathrm{L}/{\mathrm{L}}_{\mathrm{V}}^{\mathrm{o}}}{\overline{V}}_{\mathrm{L}}}}$$

(12)

where $S_{{\mathrm{W}}^{\mathrm{o}}}$ is the solute in the aqueous media outside the vesicles, $S_{{\mathrm{L}}^{\mathrm{o}}}$ is the solute associated with the outer leaflet of the vesicle membrane, $\left[L_{\mathrm{V}}\right]$ is the concentration of the lipid vesicles, ${\overline{V}}_{{\mathrm{L}}_{\mathrm{V}}^{\mathrm{o}}}$ is the volume of lipid in the outer leaflet of one mole of vesicles, which is equal to the product of the number of lipid molecules in the leaflet ($n_{\mathrm{L}/{\mathrm{L}}_{\mathrm{V}}^{\mathrm{o}}}$) and the lipid molar volume. The latter is calculated from the size of the vesicles and the area per lipid, which in the fluid phase is 60 to 70 Ų for phosphatidylcholine, sphingomyelin and phosphatidylserine [73,74,75], and 55 to 65 Ų for phosphatidylethanolamine [76,77], the area increasing with the number of cis unsaturations in the acyl chains [74,78], and being smaller when the membrane is in the gel or liquid ordered phases [74,76,78,79]. A cross-sectional area of 39 Ų is usually considered for cholesterol, independently of the membrane phase [80].

In early work there was additionally another measure used to quantify the affinity of the solute to the membrane, which was defined as the thickness of the aqueous media in contact with the surface of the membrane which contained the same moles of solute as found associated with the membrane [81]. This affinity parameter has units of length as is related with the partition coefficient by equation (13),

$$K=K_{\mathrm{P}}\hspace{0.33em}\mathcal{l}$$

(13)

where l is the thickness of the membrane to which the solute equilibrates (the whole bilayer if the solute equilibration between the two membrane leaflets is fast). Typical values for the thickness of fully hydrated bilayers composed of phospholipids with 16 carbons in the acyl chains are close to 5 nm, being somewhat larger in the gel than in the fluid phase [76]. Shorter/longer acyl chains lead to bilayers with lower/higher thickness [82]. The exact value reported for the bilayer thickness depends on the lipid groups considered, methodology used, and lipid composition [51,57,73,83,84,85,86,87,88]. Nevertheless, 5 nm can usually be considered as a reasonable estimate for bilayers prepared from the phospholipids usually found in biomembranes. The use of K as a quantitative parameter for the drug-membrane affinity may be of particular interest when evaluating the effect of lipophilicity on permeability and using membrane systems with large unstirred water layers (e.g. permeability assays using cell monolayers [89,90,91,92,93,94]).

2.1. Effect of Electrostatic Interactions and Ionic Strength on the Observed Solute-Membrane Affinity

In the partition formalisms described above, it is assumed that the solute does not affect the membrane properties and therefore the solute affinity for the membrane is not dependent on its concentration. In this case, it is not necessary to include changes in the membrane volume (or mass), nor in membrane’s charge or solvation properties. This is not the case at high solute concentrations, especially when the solute is a lipid like molecule that becomes a constituent of the membrane and may achieve very high local concentrations [47]. Due to the changes in membrane’s properties, the partition coefficient obtained will depend on solute concentration. For the comparison of different solutes and rationalization of the relation between the affinity for the membrane and solute structural properties, it is necessary to consider the affinity for an unperturbed membrane. In the case of charged solutes, it may also be relevant to distinguish between the contribution of electrostatic interactions (repulsive or attractive interactions with the charged lipid headgroups or solute-solute repulsive interactions) and hydrophobic interactions established with the non-polar portion of the membrane. Although both types of interaction contribute to the affinity observed, they depend differently on the solute concentration and on system properties such as the lipid type, pH or ionic strength. By knowing both contributions separately, well-established formalisms such as the Gouy-Chapman theory, may be used to calculate the partition coefficient at conditions other than those evaluated experimentally.

The observed partition coefficient ($K_{\mathrm{P}}$) is related to the partition coefficient in the absence of electrostatic effects ($K_{\mathrm{P}}^0$), by the membrane surface potential (${\psi }_0$) and the solute charge ($z_{\mathrm{S}}$), equation (14)

$$K_{\mathrm{P}}=K_{\mathrm{P}}^0\hspace{0.33em}{e}^{-{\displaystyle \frac{z_{S}F{\psi }_0}{RT}}}\equiv K_{\mathrm{P}}^0=K_{\mathrm{P}}\hspace{0.33em}{e}^{{\displaystyle \frac{z_{S}F{\psi }_0}{RT}}}$$

(14)

where $F$ is the Faraday constant (96485 C/mol), $R$ is the molar gas constant (8.314 J K^-1 mol^-1), and $T$ is temperature in Kelvin. This dependence of the observed affinity on the electrostatic interaction of the solute with the membrane is due to the accumulation (or depletion) of solute at the water-membrane interface, the aqueous media in the immediate vicinity of the membrane from which the solute partitions into the membrane, ${\left[S_{\mathrm{W}}\right]}_{\mathrm{int}\mathrm{er}}$, [95] equation (15):

$${\left[S_{\mathrm{W}}\right]}_{\mathrm{inter}}={\left[S_{\mathrm{W}}\right]}_{\mathrm{bulk}}{e}^{-{\displaystyle \frac{z_{S}F{\psi }_0}{RT}}}$$

(15)

After calculation of $K_{\mathrm{P}}^0$ for the specific experimental conditions used to characterize the solute affinity, the observed partition coefficient ($K_{\mathrm{P}}$) may be calculated at any other condition provided that the electrostatic interaction is the only different variable (namely a different solute concentration while at non-perturbing conditions, a different density of charged lipids in the membrane, or at a distinct ionic strength in the aqueous media).

To calculate $K_{\mathrm{P}}^0$ from $K_{\mathrm{P}}$ using equation (14) it is necessary to know the membrane surface potential. This property is related to the density of charged species in the membrane (${\sigma }_0$) by equation (16) [95,96,97]:

$${\sigma }_0={\displaystyle \frac{{\psi }_0}{\left|{\psi }_0\right|}}\sqrt{2000{ϵ}_0{ϵ}_{r}RT\sum _{\mathrm{j}}{C}_{\mathrm{j}}\left(e^{-{\displaystyle \frac{z_{\mathrm{j}}F{\psi }_0}{RT}}}-1\right)}$$

(16)

where ${ϵ}_0$ is the vacuum permittivity (8.85×10^-12 F/m) and ${ϵ}_{\mathrm{r}}$ is the relative permittivity of the medium where the surface potential is calculated (usually considered to be equal to that of the aqueous media, 78 at 25 °C). The sum in equation (16) is extended to all ions present in the aqueous media, with a concentration $C_j$ and charge $z_j$.

The membrane charge density may also be calculated directly from the charges present at the membrane surface,

$${\sigma }_0=e_0{\displaystyle \frac{\sum _i{z}_{\mathrm{i}}{n}_{\mathrm{i}}}{\sum _i{A}_{\mathrm{i}}{n}_{\mathrm{i}}}}$$

(17)

where $e_0$ is elemental electrostatic charge (1.60×10^-19 C), the sum $\sum _i{z}_{\mathrm{i}}{n}_{\mathrm{i}}$ extends to all charges in the membrane (lipids, solute and eventual ions adsorbed), and $\sum _i{A}_{\mathrm{i}}{n}_{\mathrm{i}}$ is the membrane surface area (including only the lipids if the local solute concentration is low, or both lipids and solutes if the solute molar fraction is significant [3,98]). The amount of solute in the membrane is in turn calculated from equation (6) and the value of $K_{\mathrm{P}}$ calculated from equation (14), providing an independent link between ${\sigma }_0$ and ${\psi }_0$. The charge density at the membrane surface depends on the affinity of the solute towards the membrane and in turn, the solute affinity depends on the membrane surface charge density. The membrane affinity of the solute must therefore be obtained from the convergence of the values of ${\sigma }_0$ calculated from equations (16) and (17) at different values of ${\psi }_0$ and $K_{\mathrm{P}}^0$ (and corresponding $K_{\mathrm{P}}$).

For small surface potentials ($\left|{\psi }_0\right|$ < 25 mV) and electrolyte solutions mostly composed by monovalent ions, equation (16) may be simplified to equation (18) [99,100], where ${\lambda }_{\mathrm{D}}$ is the Debye length, and $I$ is the ionic strength.

$${\psi }_0={\sigma }_0{\displaystyle \frac{{\lambda }_{\mathrm{D}}}{{ϵ}_0{ϵ}_r}}\mathrm{ };\mathrm{ }{\lambda }_{\mathrm{D}}=\sqrt{{\displaystyle \frac{RT{ϵ}_0{ϵ}_r}{2F^{2}I}}}$$

(18)

This allows the direct calculation of the surface potential from the surface charge density obtained by equation (17), without the need for iterative procedures until convergence [3].

When collecting partition coefficients from literature with the goal of establishing rules between the solute-membrane affinity and solute structural properties, both $K_{\mathrm{P}}^0$ and $K_{\mathrm{P}}$ are relevant. The former is independent on the solute and membrane charges (allowing the evaluation of the effects of solute properties other than its charge), while the latter reflects the effective affinity for some specific conditions (membrane charge and solute ionization and concentration).

Correction for electrostatic effects is usually necessary when the membrane affinity of charged molecules is characterized in vitro due to the relatively high local concentrations of solute in the membrane achieved. To evaluate the necessity to correct for the effects of electrostatics and solute concentration when considering drug distribution between membrane compartments in vivo, one should consider the expected in vivo conditions. The amount of membranes in the blood or in the tissues is quite high, with the concentration of phospholipids being in the mM range and up to 40 mM in the brain [39,101,102,103]. On the other hand, the total concentration of a given drug is low, being typically a few tens of mg per dose which leads to an overall concentration in the µM range (considering an average molar mass of 200 g/mol for the solute, 70 kg for the whole system weight with a density of 1 g/mL). Therefore, at equilibrium, the local concentration of a specific drug in biomembranes is always very low, although moderate local concentrations may be attained transiently due to slow equilibration kinetics [104,105]. Another important parameter that should be considered when evaluating the effects of charged solutes on their membrane affinity is the physiological ionic strength. This parameter depends on the specific concentration of ions in the aqueous compartments but is always high, 0.15 M being a reference value. In Figure 1 it is presented the expected effect of solute concentration on the apparent affinity of a charged solute to a lipid membrane. The lipid concentration considered was 1 mM, of which 10% are negatively charged lipids. The solute was considered to have a global charge of +1, and a moderate affinity for an uncharged membrane ($K_{\mathrm{P}}^0$=10³).

For very low ionic strength, the membrane surface potential in the absence of solute is large (79 mV at $I$= 0.01 M), and it approaches zero as the ionic strength increases (being -27 mV at $I$= 0.15 M). The electrostatic attraction between the positively charged solute and the negatively charged membrane leads to an increase in the observed affinity ($K_{\mathrm{P}}>K_{\mathrm{P}}^0$). However, as the solute concentration increases, the charges in the membrane are neutralized and the surface potential tends to zero. This leads to a decrease in $K_{\mathrm{P}}$ that approaches $K_{\mathrm{P}}^0$. At very high solute concentrations, the fraction of charged solute in the membrane becomes higher than that of the charged lipid, leading to a positive surface potential. At these conditions, $K_{\mathrm{P}}$ becomes lower than $K_{\mathrm{P}}^0$ due to the electrostatic repulsion between the solute and the charged solute-containing membrane.

For the high ionic strength characteristic in biological fluids and solutes with a single charge, the observed solute affinity is not significantly affected by the solute concentration and is similar to that observed for an uncharged membrane (within a factor of 2, red curves). In these simulations it is assumed that the solute does not influence the membrane properties other than its overall surface charge density and corresponding surface potential. This assumption is reasonable given the moderate local solute concentrations achieved ($\lesssim 10\%$ relative to the lipid). However, due to the low sensitivity of some methodologies, high solute concentrations may be required to characterize its membrane affinity leading to changes in the membrane properties. It is therefore important to estimate $K_{\mathrm{P}}^0$ from the observed $K_{\mathrm{P}}$ obtained experimentally, and the former should be the reference parameter from which the affinity at the relevant specific conditions may be calculated.

2.2. Effect of the Interaction Kinetics on the Apparent Solute-Membrane Affinity

In the formalisms presented in the previous sections, the membrane volume (and the corresponding lipid concentration) refers to the membrane portion accessible to the solute during the solute-membrane equilibration time. The equilibration of the solute between the aqueous phase and the membrane involves several steps that will limit the extent of the association during a certain time window. The effect of the kinetics of the distinct steps on the observed amount of solute associated with the membrane is addressed in this section.

The distinct steps in the equilibration of the solute between the aqueous phase and the membrane are schematically represented in Figure 2 for the case of unilamellar liposomes. The solute may be dissolved in the aqueous media outside and inside the liposomes ($S_{\mathrm{W}\mathrm{o}}$ or $S_{\mathrm{W}\mathrm{i}}$, respectively) or may be associated with the outer or inner leaflet of the lipid bilayer ($S_{\mathrm{L}\mathrm{o}}$ or $S_{\mathrm{L}\mathrm{i}}$, respectively). Depending on the solute size and polarity, the rate limiting step may be insertion in the membrane ($k_{\mathrm{i}\mathrm{n}}$), desorption from the membrane ($k_{\mathrm{o}\mathrm{u}\mathrm{t}}$), or translocation (flip-flop) between the membrane leaflets ($k_{\mathrm{f}}$). Hypothetical free energy profiles are represented by the green line, with ${\mathsf{\Delta }G}_{\mathrm{P}}^{\mathrm{o}}$ indicating the free energy variation of the solute between the aqueous phase and its equilibrium position in the membrane, and ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{X}}^{\mathrm{o}}$ the free energy variation between the distinct energy local minimum and the high energy transition states associated with the solute’s insertion ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{i}\mathrm{n}}^{\mathrm{o}}$, translocation ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{f}}^{\mathrm{o}}$, and desorption ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{o}\mathrm{u}\mathrm{t}}^{\mathrm{o}}$ steps.

In the case of relatively polar solutes, the dissolution in the non-polar center of the bilayer is the less favorable state, corresponding to an absolute energy maximum at this location (left plot), and translocation between the two membrane leaflets is the rate limiting step with an energy barrier ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{f}}^{\mathrm{o}}$. On the other hand, for very lipophilic solutes, the less favorable state is their dissolution in the very polar region of the membrane/water interface, implying that, in this case, desorption from the membrane into the aqueous phase is the rate limiting step with an energy barrier ${\mathsf{\Delta }}^{‡}{G}_{\mathrm{o}\mathrm{u}\mathrm{t}}^{\mathrm{o}}$.

In most studies of solute partition to membranes, the solute is added to the aqueous phase outside unilamellar liposomes. For liposomes with a diameter $\lesssim$ 100 nm, the fraction of aqueous media inside the liposomes is negligible, and the final step of equilibration between the inner membrane leaflet and the aqueous media does not influence significantly the overall solute distribution. In this case, the time needed to reach equilibrium (or the deviation from equilibrium at a given incubation time, $t_{\mathrm{i}\mathrm{n}\mathrm{c}}$) depends only on the rate of insertion (with a characteristic time ${\tau }_{\mathrm{i}\mathrm{n}}=1/k_{\mathrm{i}\mathrm{n}}$) and on the rate of translocation or flip-flop (with a characteristic time ${\tau }_{\mathrm{f}}=1/k_{\mathrm{f}}$). For small solutes, the energy barrier for insertion is usually small, and this step is typically completed within seconds at both the lipid concentrations used in vitro and found in vivo [68,106,107,108,109]. Depending on the solute charge and on the strength of the interactions established with the hydrated membrane, the characteristic time for solute translocation between the two membrane leaflets may vary between sub-seconds to several hours [110,111,112,113,114,115,116,117,118,119,120,121], controlling the extent of solute-membrane equilibration. For short incubation times ($t_{\mathrm{i}\mathrm{n}\mathrm{c}}$>> ${\tau }_{\mathrm{i}\mathrm{n}}$, but $\ll {\tau }_{\mathrm{f}}$), equilibration will occur only with the membrane’s outer leaflet. In this case, only the lipid in the outer leaflet should be considered in the equations (1) to (13). In the case of large unilamellar vesicles this corresponds to half the total lipid, and therefore the fraction of lipid accessible to the solute (${\gamma }_{\mathrm{L}}$) would be equal to 0.5. On the other extreme, if the system is allowed to equilibrate during at least twice the characteristic time for translocation ($t_{\mathrm{i}\mathrm{n}\mathrm{c}}\gtrsim 2 {\tau }_{\mathrm{f}}$), almost full equilibration will be achieved and ${\gamma }_{\mathrm{L}}$ ≅ 1.

Partial equilibration leads to an apparent partition coefficient ($K_{\mathrm{P}}^{\mathrm{a}\mathrm{p}\mathrm{p}}$) that underestimates the true $K_{\mathrm{P}}$. This is not a major concern when using large unilamellar liposomes where only half of the lipid is in the inaccessible inner leaflet. In this case, an incubation time equal to ${\tau }_{\mathrm{f}}$ leads to $K_{\mathrm{P}}^{\mathrm{a}\mathrm{p}\mathrm{p}}\cong 0.8 K_{\mathrm{P}}$, and $K_{\mathrm{P}}^{\mathrm{a}\mathrm{p}\mathrm{p}}\cong 0.55 K_{\mathrm{P}}$ for $t_{\mathrm{i}\mathrm{n}\mathrm{c}}={0.1\times \tau }_{\mathrm{f}}$. The solute-membrane affinity is thus underestimated by less than a factor of 2.

The underestimation of $K_{\mathrm{P}}$ may however be very significant when using multilamellar liposomes because the fraction of lipid accessible to the solute added to the outer aqueous media may be very small. As an example, if only 10 % of the lipid is in the outer leaflet of the outer lamella, $t_{\mathrm{i}\mathrm{n}\mathrm{c}}={0.1\times \tau }_{\mathrm{f}}$ leads to $K_{\mathrm{P}}^{\mathrm{a}\mathrm{p}\mathrm{p}}\cong 0.1 K_{\mathrm{P}}$. Moreover, in this case, equilibration is also dependent on the rate of desorption from the membrane. Incubation for $t_{\mathrm{i}\mathrm{n}\mathrm{c}}={2\times \tau }_{\mathrm{f}}$ leads to full equilibration with the outer lamella, but depending on the rate of desorption it may be insufficient for equilibration with the inner lamellae. For very lipophilic solutes, the characteristic time for desorption may reach several hours or even days [5,68,69,71], making incubation until full equilibration an impossible task. It is therefore of high importance to use unilamellar liposomes when characterizing the partition coefficient.

3. Methods to Obtain Solute-Membrane Affinity

3.1. Experimental Methodologies

There is a large variety of methods that may be used to characterize the association between small molecules and lipid membranes. When the membrane phase is organized as liposomes, if the aqueous media outside the liposomes may be sampled directly, this may be used to estimate the partition coefficient. Other approaches are based on changes in the properties of the solute upon association with the membrane, on the changes in the properties of the membranes due to solute association, or on changes in the properties of the whole system. Each approach has advantages and disadvantages, and there are specific formalisms for the analysis of the results to obtain the association affinity. In this section an extensive revision of the different approaches will be done, with the presentation of the formalisms required, limits of applicability and common pitfalls. The main objective is to provide the framework for a quantitative and critical interpretation of the results reported in literature for solute-membrane affinity and to facilitate its inclusion in comprehensive databases. It is also intended to guide new researchers in the field in the important task of obtaining accurate estimates for solute-membrane affinity in the case of new or poorly characterized solutes and/or membranes. Some references are provided for each methodology presented, focusing on the early works where the methodologies and analysis formalisms were developed and including also some more recent representative references.

3.1.1. Direct Quantification of the Solute in at Least One of the Phases

For an accurate quantification of the distribution of solute between the aqueous phase and the membrane it is better to quantify the amount of solute in both phases, allowing the evaluation of eventual solute degradation and/or adsorption to the apparatus. This procedure cannot however be followed to characterize the distribution towards lipid membranes because this phase is organized as liposomes which always contains large amounts of the aqueous medium. Nevertheless, the aqueous phase may often be physically separated from a membrane-enriched solution, allowing the direct quantification of the solute. The most common methods are schematically represented in Figure 3.

Separation by centrifugation requires a higher density of membranes when compared to aqueous phase. This approach has the advantage that no separation material is required, avoiding solute loss by adsorption. However, the density of phospholipid membranes is usually only slightly higher than that of pure water, requiring the use of ultracentrifugation for efficient sedimentation of unilamellar vesicles. Also, the presence of salt in the aqueous medium increases the density of this phase, decreasing the density difference or even with the aqueous phase becoming denser than the lipid vesicles. For these reasons, separation by centrifugation is usually followed when the membrane is assembled as multilamellar vesicles [122,123,124,125,126,127,128]. This may however lead to an underestimation of the membrane affinity due to incomplete equilibration (see section 1.2). More recently, unilamellar liposomes encapsulating aqueous media with high density [129,130,131], or with dense particles covered by a phospholipid membrane [132,133,134], have been used instead of MLVs, thus, allowing the easy use of centrifugation as the separation method, while minimizing the uncertainty in the fraction of lipid accessible.

Equilibrium dialysis permits the physical separation between a membrane-enriched compartment composed of regular unilamellar liposomes and is among the first methods used, still being considered as a reference [48,135,136]. The porous membrane allows equilibration of the solute with the aqueous media outside the dialysis bag while the liposomes are retained. The simplicity of this method justifies its use in most seminal works on solute-membrane interactions [48,135,137,138,139,140,141,142,143,144,145,146]. There are however important limitations, including very slow equilibration and a significant adsorption of the solute to the dialysis membrane.

Alternative approaches used to separate the membrane-enriched solution include filtration and ultrafiltration, usually performed in conjunction with centrifugation for a faster separation [123,125,132,133,134,147]. When compared with centrifugation alone, this procedure has the advantage that it does not require vesicles with a higher density than the surrounding aqueous medium. It has however the disadvantage of possible adsorption of the solute to the filter, a drawback that is shared with the most established method of equilibrium dialysis.

Size exclusion chromatography is another methodology that could be used to separate the membrane vesicles from the solute in the aqueous medium. This separation methodology has been thoroughly used to evaluate the aggregation behavior of surfactants [148,149], the interaction of solutes with small micelles and proteins [150,151], or the interaction of proteins with liposomes [152]. It has however been used only occasionally to evaluate the binding of small molecules to lipid membranes [153], due to the possible removal of solute from the liposomes during elution through the column.

Whatever the method followed for the physical separation of the aqueous phase and the membrane-enriched solution, the calculation of the partition coefficient should include explicitly the aqueous media contained in the membrane-enriched solution, [37] equation (19).

$$K_{\mathrm{P}}={\displaystyle \frac{n_{\mathrm{S}1}-n_{\mathrm{S}2}\frac{V_{\mathrm{W}1}}{V_2}}{V_1-V_{\mathrm{W}1}}}/{\displaystyle \frac{n_{\mathrm{S}2}}{V_2}}$$

(19)

The subscript 1 refers to the membrane-enriched solution and subscript 2 refers to the aqueous phase; $V_{\mathrm{i}}$ and $n_{\mathrm{S}\mathrm{i}}$ are respectively the volume and the mole number of solute in each fraction; and $V_{\mathrm{Wi}}$ is the volume of aqueous medium in the membrane-enriched solution. In Equation (19) it is assumed that the aqueous compartment has been efficiently separated from the vesicles and therefore that compartment 2 does not contain membranes.

The volume of aqueous media in the membrane-enriched compartment is relatively small when separation is performed by centrifugation. However, it is very significant when separation is performed by dialysis or filtration, usually several orders of magnitude larger than the volume occupied by the membranes themselves, and the correction introduced by equation (19) must always be taken included. An additional common source of artefacts when using dialysis is the increase in volume inside the dialysis container during equilibration, especially if it is a dialysis bag.

The most accurate procedure to calculate the partition coefficient after separation of the two solutions is to quantify the amount of solute on both. However, in some situations only the solute in the aqueous compartment is measured and the amount of solute in the membrane-enriched phase is calculated from the mass balance. This procedure overlooks eventual adsorption of the solute to the container and filtration material which may be significant for the case of non-polar solutes and leads to an overestimation of $K_{\mathrm{P}}$. When following this simplified procedure, control experiments should be performed with distinct amounts of lipid vesicles, to evaluate for consistency in the calculated partition coefficient.

The assumption that the aqueous phase does not contain lipids may also not hold, especially when the lipids tend to associate as small micelles, as is the case for lysophospholipids and fatty acids. The presence of residual amounts of lipid in the aqueous phase may significantly increase the amount of solute in the “membrane-free” compartment, leading to an underestimation of $K_{\mathrm{P}}$. This is particularly relevant for solutes with very large affinity for the membranes.

In addition to the methods indicated above, where the quantification of the solute in the aqueous phase requires physical separation from a membrane-enriched solution, there are a few methods where it may be performed with the membranes dispersed in the aqueous phase. One example is the case of volatile solutes, where their vapor pressure above the sample provides an accurate measure of the concentration in the aqueous phase. This approach has been followed to obtain the partition coefficient of ethanol to membranes of distinct lipid compositions, with the concentration of ethanol being quantified by headspace gas chromatography [154]. Other methods take advantage of a second equilibrium involving the solute in the aqueous phase. This is the case for the association of fatty acids or lysolipids with the fluorescent fatty acid binding protein ADIFAB [155,156,157,158,159,160,161,162,163]. In this approach, it must be verified if the association of the solute to the reporter binding agent does not displace the solute from association with the membranes. Otherwise, the two equilibria must be considered when calculating the solute-membrane affinity.

3.1.2. Indirect Quantification of the Solute from Changes in Its Properties

When an easily measurable property of the solute changes upon association with the membranes, this may be used to indirectly quantify the amount of solute in the membranes. This approach has been followed by a large diversity of methods based on changes of very distinct properties. The major advantage in comparison with the methodologies presented in section 2.1.1 is that the equilibrium of the solute between the aqueous phase and the membranes is not disturbed by the physical separation of the two phases. The major disadvantage is that it is an indirect measurement and therefore some general methodological aspects must be taken into consideration.

The master equation when following the variation in a solute’s property is given in equation (20):

$$P_{\mathrm{S}}=f_{S_{\mathrm{W}}}{P}_{{\mathrm{S}}_{\mathrm{W}}}+f_{S_{\mathrm{M}}}{P}_{{\mathrm{S}}_{\mathrm{M}}}$$

(20)

where $P_{\mathrm{S}}$ is the property of the solute measured at a given lipid concentration, $P_{{\mathrm{S}}_{\mathrm{W}}}$ is the same property when all solute is in the aqueous phase, and $P_{{\mathrm{S}}_{\mathrm{M}}}$ when the solute is associated with the membrane phase;$f_{S_{\mathrm{W}}}$, and $f_{S_{\mathrm{M}}}$ correspond to the fraction of solute in the aqueous phase and associated with the membrane respectively. The fraction of solute in each phase is dependent on the amount of membrane in solution and on the equilibrium association parameter, which may be defined in different ways as described in section 1.

The validity of equation (20) relies on several aspects that should be verified:

• The fraction of membrane accessible to the solute during the equilibration time must always be considered. As indicated in section 1.2, if the membrane is organized as unilamellar vesicles this introduces a relatively small uncertainty (up to a two-fold underestimation of the affinity), but corresponds to a major uncertainty if multilamellar vesicles are used.
• It is assumed that the property of interest for the solute in the aqueous media and associated with the membrane is not dependent on the lipid concentration nor on the local concentration of solute; that is $P_{{\mathrm{S}}_{\mathrm{M}}}$ and $P_{{\mathrm{S}}_{\mathrm{W}}}$ are constants in the concentration range characterized. This is a major requirement and is not always guaranteed. Frequent problems are due to attaining very high local concentrations of solute in the membrane. This may lead to changes in the properties of the membranes, to alterations in the position of the solute in the membrane, and in the interactions established. An indication regarding the validity of this assumption may be obtained by calculating the lipid/solute ratio considering the partition coefficient obtained from the best fit to the changed property. For small and uncharged solutes, a ratio larger than 20 at all lipid concentrations tested is a good indication that the membrane remains unperturbed [97,164]. An exception is when fluorescence is the property being followed. In this case the ratio should always be higher than 100 due to possible fluorescence self-quenching even at solute concentrations as low as 1 mol% [165,166,167]. A higher lipid/solute ratio should also be used for the case of large solutes or when their global charge is very high, because each solute molecule influences the properties of many lipids or leads to a strong variation in the membrane electrostatic properties, respectively. Another common problem when using very hydrophobic solutes is the use of solute concentrations above their critical aggregation concentration in the aqueous media. In this case there are two equilibria in the system, and if the aggregation equilibrium in the aqueous phase is not included in the analysis it will lead to an underestimation of the partition coefficient. The contribution from this problem may be evaluated through performing experiments at different solute concentrations, and/or by following the time evolution of the property at a given lipid concentration (see e.g. [5,168]).
• The properties of the solute when associated with the membranes ($P_{{\mathrm{S}}_{\mathrm{M}}}$) must be well characterized. This requires evaluating the property for a fractional volume of membrane above $1/K_{\mathrm{P}}$, where most of the solute is associated with the membrane. For solutes with small affinities for the membranes ($K_{\mathrm{P}}\le {10}^2$) this requires going to very high lipid concentrations ($\left[L\right]\gg 10\hspace{0.33em}\mathrm{mM}$). The upper limit of lipid concentration with the lipid phase organized as unilamellar vesicles of 100 nm diameter prepared by extrusion is around 100 mM (corresponding to 1/3 of the solution volume being occupied by the LUVs), and those high concentrations can only be achieved with charged membranes to decrease the tendency to form multilamellar vesicles. The use of liposomes to characterize solute-membrane affinity is therefore limited to solutes with $K_{\mathrm{P}}>10$. At these high lipid concentrations, the solution presents high turbidity and may lead to severe artefacts when an optical property is followed. Corrections for background signals do not always solve the problem because the contributions may be non-additive. One possible approach to evaluate the magnitude of those effects is to use distinct excitation and/or emission pathlengths and check for consistency. Alternatively, the effects may be corrected through control experiments using solutes that have been previously well characterized (e.g. [169]).

In the following sections, the most common methodologies to obtain the partition coefficient from changes in the properties of the solute are presented, with focus on method-specific aspects that are not thoroughly addressed in available reviews.

3.1.2.1. Quantification of the Solute from Changes in Its UV-Vis Spectroscopic Properties

UV-vis absorption is a convenient methodology to follow the transfer of the solute between the aqueous and the membrane phases, given the easy availability of the required equipment and the high absorptivity of most bioactive molecules at least in the UV region. However, this property depends weakly on the environment properties, and most molecules show little variations upon association with the membrane. Another important drawback when following absorption in the UV region is the strong light scattering by the membrane vesicles. These limitations are somewhat overcome when using derivative spectrophotometry [170,171], and this approach has been successfully used to characterize the partition coefficient of several molecules [3,170,171,172,173,174,175,176,177,178,179,180,181]. To avoid artifacts introduced by non-additivity of the signals due to saturation, it is highly recommended that the equipment baseline is performed with a non-absorbing solution, with the changes in the absorption of the solute due to interaction with the membrane being obtained by subtracting the absorption from the blank solutions (containing all components except the solute). When the experiments are performed in the absence of signal saturation (e.g. using low to moderate concentration range of lipid membrane), the second derivative numerical method, straightforwardly eliminates the contribution from light scattering by the liposomes suspensions, since the background scatter is nearly described by a linear function of the wavelength whose second derivative is zero [177].

In contrast with UV-vis absorption, fluorescence emission is a property that is strongly dependent on the environment. The polarity and viscosity of the lipid membrane differ significantly from that of the aqueous medium and therefore, the fluorescence properties of the solute usually change upon association with the membrane. Those properties are very convenient to evaluate the solute-membrane association due to the high sensitivity and low interference from other compounds present in the system. This approach is particularly adequate to characterize fluorescent probes due to their high fluorescence quantum yield. Nevertheless, a large fraction of bioactive molecules has aromatic groups in their structure, providing them with some fluorescence at least in the UV region.

The most commonly followed property is the variation in the fluorescence quantum yield and/or spectra [16,69,71,107,123,169,182,183,184,185,186,187,188,189,190,191,192,193], although changes in fluorescence lifetime have also been used [184,194,195]. The fluorescence quantum yield and lifetime are also influenced by the presence of other solutes in the system, and this has been used to follow the partition of the solute to the membrane. An example is the variation of the fluorescence intensity due to quenching and/or fluorescence resonance energy transfer occurring only in one of the phases [16,183,188,196,197,198,199,200]. In addition to providing information regarding the fraction of solute in each phase, this approach may also give valuable information regarding the transversal location of the solute in the membrane [183,196,201,202,203,204].

Due to the higher viscosity of the membrane when compared to the aqueous medium, the fluorescence anisotropy of the solute usually increases when associated with the membrane. This property is therefore also commonly used to follow solute-membrane association [183,184,185,205]. An important advantage of this approach is that it does not depend on solute concentration, at least if maintained low enough to prevent significant depolarization due to homo-FRET in the case of overlapping absorption and emission spectra [206,207]. It is, however, strongly influenced by scattered light, which represents an important drawback if using large vesicles and fluorescence in the UV region.

Fluorescence Correlation Spectroscopy (FCS) is another methodology that may be used to follow solute association with the membrane, relying on the much smaller diffusion coefficient of the large liposomes [131,208].

There are several reviews in literature where the specificities of each fluorescence-based approach are described in detail [42,43,206,207,209], and they will not be addressed in this review. The focus of this review is to complement the information available with a discussion on the validity limits of the Master equation (20), namely regarding the additivity of the contribution from the solute in each phase. With this goal, the results expected at some typical situations were simulated. The dependence of several fluorescence properties with the lipid concentration and the partition coefficient obtained from the best fit of the Master equation are shown in Figure 4. The variable properties considered were the fluorescence quantum yield (${\varphi }_{\mathrm{F}}$) and the wavelength at maximum fluorescence emission (${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$). The conditions simulated are: i) large (10-fold) increase in ${\varphi }_{\mathrm{F}}$ upon partition to the membrane without changes in ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$; ii) moderate (2.5 fold) increase in ${\varphi }_{\mathrm{F}}$ accompanied by a moderate variation in ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (50 nm blue-shift); iii) moderate variation in ${\varphi }_{\mathrm{F}}$ (2.5-fold) and large variation ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (100 nm blue-shift); and iv) no variation in ${\varphi }_{\mathrm{F}}$ and large variation ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ (100 nm blue-shift). For simplicity it was assumed that the solute absorption spectrum (shape and molar absorptivity) is the same when in the aqueous or membrane phases. It is also assumed that the solute does not change its ionization state upon association with the membranes, which means that no additional equilibria must be considered other than solute equilibration between the aqueous and membrane phases. The ratio between the partition coefficient estimated from the best fit of the simple master equation ($K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$), with $f_{S_{\mathrm{M}}}$ calculated from equation (6), and the true $K_{\mathrm{P}}$ (fixed at 2×10³), is included in Figure 4.

The first observation is that the dependence of the fluorescence property with the lipid concentration is usually well described by the master equation (20) for all properties evaluated. However, the value of $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ obtained from the best fit may deviate from the true $K_{\mathrm{P}}$ by more than an order of magnitude. Therefore, the quality of the best fit provides no evidence regarding the validity of the approximations assumed in the derivation of the master equation.

The first property analyzed was the fluorescence intensity at a fixed wavelength ($I_{\mathrm{F}}$), first set of plots below the spectra. This is a very robust property, with $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ being always equal to the true $K_{\mathrm{P}}$ irrespective of the wavelength selected.

The adequacy of the maximum fluorescence intensity ($I_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$) was also evaluated (black curves in the middle plots). In the absence of a spectra shift this is equivalent to the intensity at a fixed wavelength, leading to $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}=K_{\mathrm{P}}$. However, as the spectral shift becomes significant, the overlap between the spectra of the solute in the aqueous media and the membrane decreases, and this property leads to incorrect estimations of the partition coefficient. If the spectra overlap is very small (as in case iv), the shape of the variation is not well described by the Master equation, with the results at low and high lipid concentrations leading to different $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$.

The selection of intensive fluorescence properties may be of high value to characterize the solute-membrane affinity using high throughput methodologies where solute concentration may show some variations along the distinct datapoints. In the case of significant spectral shift, a convenient property could be the wavelength of maximum fluorescence intensity (${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$). The dependence of ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is shown in the middle plots (magenta curves). For small spectral shifts maintaining a large overlap between the emission spectra of the solute in water and in the membrane (case ii), a reasonably good fit is obtained but $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ leads to an over-estimation of $K_{\mathrm{P}}$. This may however be corrected if weighed by the relative fluorescence intensity of the solute in both media. In this case $I_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\mathrm{S}}_{\mathrm{M}}\right)/I_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\mathrm{S}}_{\mathrm{W}}\right)$ = 2.5, and an accurate estimate of $K_{\mathrm{P}}$ is obtained when $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ is divided by this ratio. However, as the spectral shift increases (cases iii and iv), poor fits are obtained. In the simulations shown, considering lipid concentrations that lead up to 90% of the solute in the membrane, although the quality of the best fit is not good, $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ leads to a reasonable estimate of $K_{\mathrm{P}}$ when divided by the relative intensities in both media ( $I_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\mathrm{S}}_{\mathrm{M}}\right)/I_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\mathrm{S}}_{\mathrm{W}}\right)$= 2.5 for case iii and 1.0 for case iv). However, in a real situation, the poor quality of the best fit will lead to $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ values that vary with the lipid concentration used. Showing that this property is not robust enough for the estimation of partition coefficients.

In the second plots from the bottom, the variation of the signal at a given wavelength is represented after normalization of the different spectra at their own maxima, $I_{\mathrm{F}} \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{a}$, with the color codes corresponding to the emission wavelengths indicated by the arrows in the top plots. As expected, no variation is observed in the absence of a spectral shift (case i). For a small spectral shift and moderate increase in the fluorescence quantum yield (case ii), the shape of the curve is well described by the Master equation at all wavelengths. However, the best fit leads to incorrect estimations of $K_{\mathrm{P}}$, usually leading to an overestimation. A good agreement was obtained when the wavelength of maximum emission in the aqueous phase was selected (blue curves), but this adequacy is not robust enough to be used for the characterization of unknown partition coefficients. As the spectral shift increases (cases iii and iv), the curve shape is no longer well described by the Master equation. In the case of a large spectral shift without changes in quantum yield (case iv), selecting the isoemissive wavelength reveals two distinct regimes: an increase in normalized fluorescence intensity at low lipid concentrations and a decrease at high lipid concentrations. Moreover, each regime leads to a different fitted partition coefficient (green dark solid line, and light green dashed line).

Fluorescence anisotropy is another intensive parameter that was evaluated in Figure 4, lower plots. The shape of the curve is always very well described by the Master equation. However, the value obtained for $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ deviates significantly from $K_{\mathrm{P}}$ at some of the conditions considered. The inadequacy of variations in the fluorescence anisotropy as a direct reporter of solute partition to the membrane has been clearly stated before (e.g. [42]) and is due to the fact that the contribution of the distinct solute species (in the aqueous or membrane compartment) to the overall parameter is not proportional to their fraction. This may be corrected by weighting the fraction of the species by their contribution (${\alpha }_{\mathrm{W}}$ or ${\alpha }_{\mathrm{M}}$), equation (21), [42,184,210].

$$P_{\mathrm{S}}={\displaystyle \frac{{\alpha }_{\mathrm{W}}{ f}_{S_{\mathrm{W}}}{P}_{{\mathrm{S}}_{\mathrm{W}}}+{\alpha }_{\mathrm{M}}{ f}_{S_{\mathrm{M}}}{P}_{{\mathrm{S}}_{\mathrm{M}}}}{{\alpha }_{\mathrm{W}}{ f}_{S_{\mathrm{W}}}+{\alpha }_{\mathrm{M}}{ f}_{S_{\mathrm{M}}}}}$$

(21)

The contribution of each species depends on the intensity of light absorbed, which is proportional to their molar absorptivity at the excitation wavelength (${\epsilon }_{\lambda \mathrm{exc}}^{{\mathrm{S}}_{\mathrm{x}}}$) if absorption is lower than 0.1, and on their fluorescence emission intensity at the emission wavelength considered, which depends on the fluorescence quantum yield (${\varphi }_{\mathrm{F}}^{{\mathrm{S}}_{\mathrm{x}}}$) and on the fraction of emitted light at the collected wavelength ($g_{\lambda \mathrm{em}}^{{\mathrm{S}}_{\mathrm{x}}}$) corresponding to the fluorescence intensity at that wavelength on a fluorescence spectra normalized by its area, equation (22).

$${\alpha }_{\mathrm{W}}={\epsilon }_{\lambda \mathrm{exc}}^{{\mathrm{S}}_{\mathrm{W}}}\hspace{0.17em}{\varphi }_{\mathrm{F}}^{{\mathrm{S}}_{\mathrm{W}}}\hspace{0.17em}{g}_{\lambda \mathrm{em}}^{{\mathrm{S}}_{\mathrm{W}}}\hspace{0.33em};\hspace{1em}{\alpha }_{\mathrm{M}}={\epsilon }_{\lambda \mathrm{exc}}^{{\mathrm{S}}_{\mathrm{M}}}\hspace{0.17em}{\varphi }_{\mathrm{F}}^{{\mathrm{S}}_{\mathrm{M}}}\hspace{0.17em}{g}_{\lambda \mathrm{em}}^{{\mathrm{S}}_{\mathrm{M}}}\hspace{0.33em}$$

(22)

When the anisotropy obtained from the simulated results is analyzed with the equation resulting from (21) and (22), the value of $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}$ is essentially equal to $K_{\mathrm{P}}$, reflecting the fact that anisotropy is an additive property if weighted by the fractional fluorescence intensity from each species [211]. In contrast, the related property fluorescence polarization is not additive even when weighted by the fractional fluorescence intensity [211], and it cannot therefore be used to obtain the affinity of the solute to the membrane.

The very large number of parameters in equation (21) may discourage the best intended researcher and may justify the inadequate use of the simple master equation in many publications. It should however be noted that the individual parameters do not need to be known, only the ratio between their product when the solute is in the membrane and in water, ${\alpha }_{\mathrm{M}}/{\alpha }_{\mathrm{W}}$. These parameters may be easily obtained through the relative fluorescence intensity in the aqueous phase and in the presence of sufficiently high amounts of membrane phase, permitting the use of the correct formalism without too much effort (as done e.g. in [212]). When the fluorescence intensity from both species is similar at the wavelength selected (namely at the isoemissive wavelength, shown in green), $K_{\mathrm{P}}$ may in fact be obtained directly from the best fit of equation (20).

In Figure 4 it was assumed that the absorption spectrum is the same in both media. The effect of changing the excitation wavelength when there are variations in the absorption spectra is explored in Figure 5. A common situation was considered where the absorption and the emission spectra are shifted to higher energies (lower wavelengths) when the polarity of the medium decreases (${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{A}}^{\mathrm{W}}$=320 nm and ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{A}}^{\mathrm{M}}$=300 nm, and ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{F}}^{\mathrm{W}}$=450 nm and ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{F}}^{\mathrm{M}}$=400 nm). It was also assumed that the molar absorptivity at the maximum absorption was unchanged, while a higher fluorescence quantum yield was considered when associated with the membrane (${\mathsf{\Phi }}_{\mathrm{F}}^{\mathrm{M}}/{\mathsf{\Phi }}_{\mathrm{F}}^{\mathrm{W}}$= 2), the anisotropy was the same as in the simulations shown in Figure 4, $r_{\mathrm{W}}$= 0.01 and $r_{\mathrm{M}}$= 0.2.

The first observation from Figure 5 is that if the combination of the excitation and emission wavelengths favors the species in the aqueous media, the apparent partition coefficient obtained when the anisotropy variation is fitted with the unweighted equation (20) is lower than $K_{\mathrm{P}}$ (wavelengths and anisotropies marked in blue). Conversely, if the species in the membrane are favored then $K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}>K_{\mathrm{P}}$ (wavelengths and anisotropies marked in orange). On the other hand, if the fluorescence emission is collected at a wavelength where the fluorescence intensity is similar when the solute is in the water or associated with the membrane, ${\alpha }_{\mathrm{M}}/{\alpha }_{\mathrm{W}}\cong$1, a good estimate of $K_{\mathrm{P}}$ is obtained from the fit of the unweighted equation (wavelengths and anisotropy variations marked in green). The wavelength where this is observed depends on the excitation wavelength − closer to ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{F}}^{\mathrm{W}}$ at excitation wavelengths that preferentially excite the solute in the membrane, bottom to top plots. Therefore, even if the absorption spectra, molar absorptivities and/or fluorescence quantum yields are unknown, an accurate estimation of $K_{\mathrm{P}}$ may be obtained provided that the fluorescence emission is collected at a wavelength where both species show similar fluorescence intensity.

In any case, $K_{\mathrm{P}}$ may always be calculated from $K_{\mathrm{P}}^{\mathrm{a}\mathrm{p}\mathrm{p}}$ provided that ${\alpha }_{\mathrm{M}}/{\alpha }_{\mathrm{W}}$ is known, equation (23), which is obtained from the equality of equation (24) and (25).

$$K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}=K_{\mathrm{P}}{\displaystyle \frac{{\alpha }_{\mathrm{M}}}{{\alpha }_{\mathrm{W}}}}$$

(23)

$$r={\displaystyle \frac{{\alpha }_{\mathrm{W}}{f}_W r_{\mathrm{W}}+{\alpha }_{\mathrm{M}} f_{\mathrm{M}} r_{\mathrm{M}}}{{\alpha }_{\mathrm{W}}{f}_W+{\alpha }_{\mathrm{M}} f_{\mathrm{M}}}},\hspace{0.33em}\mathrm{with} f_{\mathrm{W}}={\displaystyle \frac{1}{V_{\mathrm{W}}+K_{\mathrm{P}}{V}_{\mathrm{M}}}} \mathrm{and} f_{\mathrm{M}}={\displaystyle \frac{K_{\mathrm{P}}{V}_{\mathrm{M}}}{V_{\mathrm{W}}+K_{\mathrm{P}}{V}_{\mathrm{M}}}}\hspace{1em}\hspace{0.33em}$$

(24)

$$r=f_{\mathrm{W}}^{\mathrm{F}\mathrm{i}\mathrm{t}} r_{\mathrm{W}}+f_{\mathrm{M}}^{\mathrm{F}\mathrm{i}\mathrm{t}} r_{\mathrm{M}},\hspace{0.33em}\mathrm{with} f_{\mathrm{W}}^{\mathrm{F}\mathrm{i}\mathrm{t}}={\displaystyle \frac{1}{V_{\mathrm{W}}+K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}{V}_{\mathrm{M}}}} \mathrm{and} f_{\mathrm{M}}^{\mathrm{a}\mathrm{p}\mathrm{p}}={\displaystyle \frac{K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}{V}_{\mathrm{M}}}{V_{\mathrm{W}}+K_{\mathrm{P}}^{\mathrm{F}\mathrm{i}\mathrm{t}}{V}_{\mathrm{M}}}}\hspace{1em}\hspace{0.33em}$$

(25)

In conclusion, the best property to report on the association of fluorescent solutes to membranes is the fluorescence intensity at a fixed wavelength. This property is robust and provides good direct estimates in all cases. If the conditions of the experiment introduce variability in the concentration of solute, fluorescence anisotropy may be used, but it must be weighted by the relative fluorescence intensity of the solute in the distinct phases. Unfortunately, simple corrections with general validity cannot be applied to the convenient properties of spectral shift or $I_{\mathrm{F}} \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{a}$ to estimate the true $K_{\mathrm{P}}$. Spectral shifts may be used but only if the spectra of the solute in water and in the membrane strongly overlaps (case ii in Figure 4). In this case a good fit of the Master equation is expected and the true $K_{\mathrm{P}}$ may be obtained from equation (23). These properties should therefore be used with caution, preferably supported by alternative independent approaches.

3.1.2.2. Quantification of the Solute from Changes in Its Ionization Equilibria

The variation of solute ionization equilibria is a property that has also been used often to characterize the association of weak acids and bases with membranes [48,144,213,214,215,216,217,218]. This may be obtained directly by potentiometric titration or through variations in the solute properties at different pH values in the presence of varying concentrations of lipid.

The equilibria that must be considered includes the ionization of the solute in the aqueous phase and when associated with the membrane, as well as the partition coefficient of all solute species for the membrane. In the case of solutes with a single ionization equilibrium, the relevant equilibria are shown in Figure 6.

The four equilibria establish a thermodynamic cycle, so the equilibrium constants are constrained by the relation:

$$K_{\mathrm{P}}^{\mathrm{SH}}{ K}_{\mathrm{a}}^{\mathrm{M}}=K_{\mathrm{P}}^{\mathrm{S}}{ K}_{\mathrm{a}}^{\mathrm{W}} ; \mathrm{p}{K}_{\mathrm{a}}^{\mathrm{M}}=\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{W}}-\mathrm{l}\mathrm{o}\mathrm{g}\left({\displaystyle \frac{K_{\mathrm{P}}^{\mathrm{S}}}{K_{\mathrm{P}}^{\mathrm{SH}}}}\right)$$

(26)

At a given pH and volume of membrane, the solute is distributed between the distinct species and compartments, showing a global, observed, ionization constant ($K_{\mathrm{a}}^{\mathrm{o}\mathrm{b}\mathrm{s}}$) given by:

$$K_{\mathrm{a}}^{\mathrm{o}\mathrm{b}\mathrm{s}}={\displaystyle \frac{{nS}_{\mathrm{W}}+{nS}_{\mathrm{M}}}{{nSH}_{\mathrm{W}}+{nSH}_{\mathrm{M}}}}{ \left[H^{+}\right]}_{\mathrm{T}};$$

(27)

$$\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{obs}}=\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{W}}+\mathrm{l}\mathrm{o}\mathrm{g}\left({\displaystyle \frac{1+K_{\mathrm{P}}^{\mathrm{SH}}\nu }{1+K_{\mathrm{P}}^{\mathrm{S}}\nu }}\right)$$

(27)

where $\nu =V_{\mathrm{M}}/V_{\mathrm{W}}\cong V_{\mathrm{M}}/V_{\mathrm{T}}$, and the observed, global, partition coefficient ($K_{\mathrm{P}}^{\mathrm{o}\mathrm{bs}}$) is given by:

$$K_{\mathrm{P}}^{\mathrm{obs}}={\displaystyle \frac{\left({ \left[H^{+}\right]}_{\mathrm{T}}+K_{\mathrm{a}}^{\mathrm{M}}\right)K_{\mathrm{P}}^{\mathrm{SH}}}{\left({ \left[H^{+}\right]}_{\mathrm{T}}+K_{\mathrm{a}}^{\mathrm{W}}\right)}}$$

(28)

Simulations for the dependence of $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{obs}}$ on $\nu$ for some typical situations are shown in Figure 7.

The different lines in each plot correspond to different magnitudes of stabilization of the partition coefficient of the neutral species ($AH$ or $B$, for acids and bases, respectively). In the case of equal membrane affinities ($K_{\mathrm{P}}^{\mathrm{A}\mathrm{H}}=K_{\mathrm{P}}^{{\mathrm{A}}^{-}}$, or $K_{\mathrm{P}}^{\mathrm{B}}=K_{\mathrm{P}}^{{\mathrm{B}\mathrm{H}}^{+}}$), the ionization equilibrium is unaffected by partition to the membrane, and no variation is observed in the $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{obs}}$ (grey lines). Stabilization of the neutral species by the membrane leads to an increase in its fraction at higher lipid concentrations, leading to an increase in $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{obs}}$ in the case of acids (upper plots) and a decrease in the case of bases (lower plots). The effect of the lipid concentration on $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{obs}}$ depends on the lipophilicity of the solute, with less than one $\mathrm{p}{K}_{\mathrm{a}}$ unit at 4 mM lipid for solutes with moderate lipophilicity ($\mathrm{l}\mathrm{o}\mathrm{g}{K}_{\mathrm{P}}^{\mathrm{A}\mathrm{H}}$= 3 or $\mathrm{l}\mathrm{o}\mathrm{g}{K}_{\mathrm{P}}^{\mathrm{B}}$= 3, rigth plots), and up to several $\mathrm{p}{K}_{\mathrm{a}}$ units for very lipophilic solutes ($\mathrm{l}\mathrm{o}\mathrm{g}{K}_{\mathrm{P}}^{\mathrm{A}\mathrm{H}}$= 5 or $\mathrm{l}\mathrm{o}\mathrm{g}{K}_{\mathrm{P}}^{\mathrm{B}}$= 5, left plots).

One difficulty of using potentiometric titration to characterize the interaction of ionizable solutes with membranes is the high concentrations of solute that must be used to allow a direct pH titration. This leads to very high local concentrations in the membrane that may alter the surface electrostatic potential of the membrane and influence the affinity observed. Corrections for electrostatic effects are therefore usually required to characterize the association of the charged solute species with the membranes. This is achieved by substituting the partition coefficients in equation (27) by the intrinsic partition coefficients that depend on the electrostatic surface potential of the membrane, see equations (14) to (18). This was performed in reference [216] to obtain the intrinsic partition coefficients of several weak acids and bases. Some additional challenges with the application of this methodology include the presence of ionizable groups in the lipids themselves, and the possible displacement of the ionization equilibria in the aqueous media at the water/membrane interface due to accumulation or depletion of ${\mathrm{H}}^{+}$ and/or ionized solute species (see e.g. reference [48] for details).

Several additional methodologies may be used to indirectly follow solute ionization in the presence of the membrane and thus obtain the partition coefficient of the distinct species and $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{M}}$. This may be achieved by obtaining $K_{\mathrm{P}}^{\mathrm{o}\mathrm{b}\mathrm{s}}$ at distinct pH values, or through the characterization of $\mathrm{p}{K}_{\mathrm{a}}^{\mathrm{o}\mathrm{b}\mathrm{s}}$ at distinct lipid concentrations. The only requirement is that the methodology is able to distinguish between the different solute ionization species when in the aqueous medium and associated with the membrane. All concerns discussed in the preceding sections will apply, namely the effects of the solute on membrane properties, effects of the interaction kinetics, and liposome multilamellarity (discussed in section 1), as well as issues with additivity of the measured property (as discussed in section 2.1.2.1).

Isothermal Titration Calorimetry (ITC) is of particular relevance to characterize partition and ionization and will be discussed in detail in section 2.1.3.1. Fluorescence based methodologies are also useful since ionization usually influences the fluorescence spectra and/or quantum yield and because its high sensitivity allows the use of very low solute concentrations providing the parameters in unperturbed membranes. However, this methodology is difficult to implement due to changes in the fluorescence properties of the solute when associated with the non-polar membrane, the full characterization requiring knowledge of the absorption and fluorescence spectra of all species in the aqueous medium and associated with the membrane [212]. Other relevant methodologies are briefly discussed in the next section.

The characterization of the partition and ionization equilibria for ionizable solutes is of major practical importance when the goal is to predict the behavior of the solutes in vivo, namely their Absorption, Distribution, Metabolism, Elimination and Toxicity (ADME/Tox) as well as their pharmacodynamics. This is because all those processes are influenced by the localization of the solute in the complex in vivo system (namely the aqueous compartments and biomembranes), and because the pH varies in the distinct aqueous compartments of the cells, tissues and organs. The charge of the solute when associated with the biomembranes has a strong influence on its rate of membrane permeation (with direct effects on its pharmacokinetics), and on its transversal location in the membrane (influencing their interaction with membrane proteins and therefore on their pharmacodynamics). Due to its complexity, there is however a very limited number of studies where the interplay between solute ionization and membrane association has been characterized in detail. Hopefully, the computational tools more recently available will facilitate the analysis and the interpretation of complex sets of data and will contribute to overcoming these limitations.

3.1.2.3. Other Methodologies Based on Changes in Solute Properties

Nuclear Magnetic Resonance (NMR) is a powerful technique that has been used to study solute-membrane interactions. This is a powerful and versatile methodology where the resonance may be followed from distinct nuclei (most commonly ¹H, ²H, ¹³C, ³¹P, ¹⁹F), in solution and solid-state NMR, with many variations in the pulse sequences used for the spectra acquisition including 1-dimensional and multi-dimensional NMR. A major advantage of this methodology is the possibility of obtaining the localization of the solute in the membrane [219,220,221]. Effects of the solute on membrane properties can also be obtained, including changes in the lipid order(e.g. [24]) or membrane curvature (e.g. [222]) due to intercalation of the solute, and surface electrostatics due to changes in the head group density or association of charged solutes (e.g. [223]). The resonances from the solute nuclei are not only sensitive to the polarity of the environment but also to changes in ionization, thus allowing the characterization of changes in pKa upon association with the membrane [224]. In some conditions it is possible to distinguish between the solute associated with the outer and inner leaflet of the membrane, thus allowing the characterization of the rate of permeation [225,226]. Information about the rate of solute equilibration with the membrane can also be obtained through the dependence of the solute resonance frequency and width with the lipid concentration, with slow equilibration originating resonances at two distinct frequencies for the solute in water and in the membrane [227,228,229,230] while in fast equilibration resonance at a single intermediate frequency is obtained [226,230,231,232,233,234,235].

In spite of the immense potential, there are very few reports on the use of NMR to obtain partition coefficients [228,232,233,234,235,236,237], most using ¹⁹F [228,233,234] but also ¹H [232,235,237] and ²H [229,236]. In most of the publications, only qualitative or semi-quantitative information is obtained [238,239,240,241,242,243,244,245]. The major difficulty is signal broadening due to the slow tumbling of the large lipid aggregates. This is overcome through the use of small lipid aggregates such as small unilamellar vesicles (SUVs) [232,234] or lipid-surfactant bicelles [237], by following only the unbound solute [233,235], or through the use of large lipid vesicles in solid-state NMR with fast spinning [228]. All strategies have their limitations, with the lipid small aggregates being poor model systems for biomembranes, existence of non-linear effects when following only the unbound solute, and the use of lipid membranes not fully hydrated when using solid-state NMR. Another significant challenge is the inherently low sensitivity of NMR, which often requires solute concentrations in the millimolar range. This can cause solubility issues in the aqueous phase and lead to high local concentration in the membrane, potentially perturbing its properties. Sensitivity improves with increasing NMR magnetic field strength and with the use of cryogenic probes, and significant advances are expected in the coming years for applying this technique to the characterization of solute–membrane interactions.

Another technique that has been occasionally used to evaluate solute-membrane interactions is Electron Paramagnetic Resonance (EPR). Because most biologically relevant small molecules do not present paramagnetic properties, these studies are mostly limited to lipophilic electron spin probes [18,245,246,247,248,249,250], or solutes labelled with electron spin probes [251,252,253]. Indirect information on solute-membrane interactions has also been obtained from the effect of paramagnetic solutes on the NMR signal of membrane lipids [226,254], or solute effects on membrane properties reported by EPR probes [244,245,255,256,257,258]. However, these indirect studies are only qualitative or semi-quantitative and reflect association at membrane perturbing conditions.

3.1.3. Indirect Quantification of the Solute from Changes in the Properties of the Whole System and/or in the Membrane Properties

Due to the low concentration of solute expected in vivo (see section 1.1) and in order to enable comparisons among different solutes, membrane compositions, and experimental conditions, it is generally preferable to characterize membrane affinity using very low solute concentrations, thereby avoiding membrane perturbation. However, most indirect methods rely on detecting changes in membrane properties and therefore typically require high local solute concentrations, often approaching membrane saturation.

This section presents the most common methods used to characterize solute–membrane partitioning that rely on changes in the properties of the membranes and/or the aqueous medium. The methods discussed range from isothermal titration calorimetry (ITC), which allows determination of affinity at both low and high solute concentrations, to techniques based on changes in membrane properties, such as electrostatic characteristics, liposome size, and solid–liquid phase transitions. Finally, methods based on changes in the properties of supported lipid membranes are briefly discussed.

3.1.3.1. Isothermal Titration Calorimetry (ITC)

In ITC, the property measured is the heat variation upon transfer of the solute from the aqueous media to the membrane, which depends on solute-solute, solvent-solvent, and solute-solvent interactions. In the more conventional approach, aliquots of a concentrated membrane solution are added to a solute solution and the heat evolved after each addition is recorded. If the solvent in the membrane and solute solutions is perfectly matched and do not alter their properties in the concentration range evaluated (no heat of dilution), the heat evolved due to the addition of the membrane solution aliquot $i$, is proportional to the enthalpy variation of the association (${\Delta H}^{\mathrm{o}}$) and to the variation in the number of moles of solute associated with the membrane due to the increase in the lipid concentration, ${∆}^{\mathrm{i}}{n}_{\mathrm{SL}}$, equation (29).

$$q_{\mathrm{i}}={\mathsf{\Delta }H}^{\mathrm{o}}\hspace{0.33em}{\mathsf{\Delta }}^{\mathrm{i}}{n}_{\mathrm{SL}}$$

(29)

The variation in the moles of solute associated with the membrane may be calculated from the dependence of $q_{\mathrm{i}}$ with the lipid concentration, allowing calculation of the association affinity. Some methodological details may be consulted for example in references [97,164,259,260].

The sensitivity of the initial ITC instruments was in the range of μcal s^-1 and therefore high concentrations of solute were needed, leading to very high local concentrations of solute in the membrane. In fact, most studies were dedicated to the study of membrane solubilization by high concentrations of amphiphilic solutes (e.g. [261]). However, modern ITC equipment have sensitivity in the range of ncals^-1, allowing the use of much lower solute concentrations and thus the characterization of solute association to essentially unperturbed membranes [97,164]. It should be noted that membrane perturbation may be particularly problematic in experiments with the successive addition of membrane-enriched aliquots because the different titration steps are not independent and because higher local concentrations of solute are attained in the first titration points. Membrane perturbation is therefore higher in the beginning of the titration, with the perturbation propagating to all the titration curve.

The simpler model to describe the dependence of ${∆}^{\mathrm{i}}{n}_{\mathrm{SL}}$ with the lipid concentration for small and uncharged solutes is considering a simple partition, equation (1) (e.g. [164]). If the solute is charged, electrostatic effects may need to be included in the model, allowing the characterization of the partition coefficient at the experimental conditions followed ($K_{\mathrm{P}}$), and the intrinsic partition coefficient in the absence of electrostatic interactions ($K_{\mathrm{P}}^0$), equation (14) (e.g. [96,242]). If high local solute concentrations in the membrane are attained during the titration and/or the ionic strength of the aqueous media is low, the electrostatic potential at the membrane surface changes throughout the titration and the data must be analyzed with the complete model including electrostatic effects (details in section 1.1 and analysis applications provided in references [97,98,164]). However, if the aqueous media contains high concentrations of ions (e.g. 0.15 M NaCl), if the solutes have a single charged group, and solute local concentrations in the membrane are lower than 5 mol%, the membrane surface potential is small and maintained throughout the titration, allowing the use of the simple partition model to analyze the data and the calculation of the intrinsic partition coefficient using equation (14) [97].

The partition models discussed in section 1 are not available in the analysis software provided by the ITC manufacturer’s, which is focused on ligand-protein binding. For that reason, some authors analyze the results assuming binding of the solute to well defined binding sites in the membrane (e.g. [262]). Although this model does not describe adequately the association of small solutes with lipid membranes, the quality of the best fit obtained is usually good, providing estimates of the binding constant ($K_{\mathrm{b}}$) and the size of the binding sites in the membrane (number of lipid molecules per binding site, ${\#}_{\mathrm{L}}$), from which the partition coefficient may be calculated using equation (10). A major concern when following this approach is that the estimates obtained for $K_{\mathrm{b}}$ and ${\#}_{\mathrm{L}}$ correspond to a set of values in a wide range of possible parameter combinations. This is illustrated in Figure 8 for the case of association of the antibiotic rifampicin (data from reference [263]) and the surfactant deoxycholic acid (data from reference [119]) with POPC LUVs. If the number of binding sites available in the membrane is much larger than the number of solute molecules bound, the quality of the best fit is independent on ${\#}_{\mathrm{L}}$, and only the ratio $K_{\mathrm{b}}/{\#}_{\mathrm{L}}$ may be obtained. The value of $K_{\mathrm{b}}$ cannot therefore be used as a direct estimate of the solute-membrane affinity. However, this procedure allows the calculation of accurate estimates for $K_{\mathrm{P}}$, using equation (10), and $\mathsf{\Delta }{H}^{\mathrm{o}}$. The yellow box in the right plots of Figure 8 highlights the range of ${\#}_{\mathrm{L}}$ values that lead to a best fit quality similar to that obtained considering the simple partition model (${\chi }^2\le 1.2 {\chi }_{\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}^2$, corresponding to the upper limits of a 80% confidence interval). For the case of titration of 10 µM rifampicin this interval corresponds to ${\#}_{\mathrm{L}}$ up to 30, while the maximum number of lipids per binding site is decreased to 15 in the case of titration of 25 µM deoxycholic acid. The different intervals reflect the distinct maximum local concentrations achieved, with 90 lipid molecules per solute bound in the case of rifampicin, and 50 lipid molecules per solute bound in the case of deoxycholic acid. At the limit of the 80% confidence interval, the value of $K_{\mathrm{P}}$ is overestimated by approximately two-fold and $\mathsf{\Delta }{H}^{\mathrm{o}}$ is underestimated by about 50%. The deviations from the true parameters are comparable to other experimental and model uncertainties such as the fraction of lipid accessible due to liposome multilamellarity and/or slow solute translocation (see section 1 for details). Therefore, provided that the solute-membrane affinity is reported as ${n K}_{\mathrm{b}}$, $n$ being the number of binding sites per lipid, the best fit using the one set of sites model available in the analysis software provided by the ITC manufacturer is a reasonable approach to describe the association of small molecules with lipid membranes. Given the comparable size of the solutes and lipids, it is more meaningful to define the binding sites in terms of the number of lipid molecules in each binding site (${\#}_{\mathrm{L}}$), and not as the number of binding sites per lipid ($n$), with ${\#}_{\mathrm{L}}=1/n$. In any case, this parameter is usually not very informative given that best fits with similar quality are obtained from distinct sets of parameter values ($n$ and $K_{\mathrm{b}}$). The partition coefficient ($K_{\mathrm{P}}$) may however be calculated from $n K_{\mathrm{b}}/{\overline{V}}_{\mathrm{L}}$ or the equivalent $K_{\mathrm{b}}/\left({{\#}_{\mathrm{L}}\overline{V}}_{\mathrm{L}}\right)$, equation (10).

The Langmuir adsorption model (Langmuir isotherm) has also been used to describe solute association with lipid membranes (e.g. [264]). This model is in fact equivalent to the one set of sites model with the size of the binding sites described as a surface area.

In addition to providing acceptable estimates of the partition coefficients and association enthalpy variations for small solutes, the explicit consideration of binding sites in the membrane may in fact be the most adequate models for the interaction of some specific solutes and membranes, namely in the case of large solutes with multiple charged groups and membranes with high surface charge density (e.g. [66]). It has nevertheless been shown that these systems may usually still be well described by partition models provided that the contributions from electrostatic interactions are included (e.g. [265]).

Complementing the most usual approach of titrating the solute with increasing concentrations of lipid, some applications are based on the titration of the membrane with increasing concentrations of solute. This approach has mostly been used with only a few titration steps to avoid large variations in the solute concentration with the goal of obtaining the interaction enthalpy and the assumption of full solute binding (e.g. [266]). When a full titration is performed, the heat evolved usually decreases with the increase in the solute concentration, and this has been interpreted as saturation of independent binding sites in the membrane, and analyzed with the one set of sites model, scheme (7) and equation (8) (e.g. [66,262,267]. However, the assumption of independent binding sites in the membrane is not usually valid, and while for solute titration with lipid this analysis may provide reasonable estimates of the association affinity, the decrease in the amount of solute that associates with the membrane at high solute concentrations is in fact due to changes in the membrane properties.

A combination of these two approaches (lipid addition to solute in the cell or solute addition to lipid in the cell) has been proposed to characterize the membrane-association of very low affinity solutes [263], and the data obtained from both sets of titrations may provide a critical evaluation of the use of binding sites to characterize solute-membrane affinity. In that approach, several independent titrations are performed by adding the solute to distinct lipid concentrations in the cell. In each titration, increasing solute concentrations are achieved, providing information regarding membrane saturation from the dependence of the heat evolved ($q_{\mathrm{i}}$) with solute concentration. The association to unperturbed membranes, may however be obtained from the dependence of $q_{\mathrm{i}}$ of the first solute addition for independent titrations with distinct lipid concentrations. This approach has two advantages. First, the lipid concentrations that may be attained in cell are much higher than when following the usual approach of adding the lipid membranes to the solute in the cell. Second, all titration data points are independent and report the interaction with unperturbed membranes. A drawback is however the need to perform several independent titrations to obtain a single estimate of the partition coefficient. To assess the adequacy of obtaining solute affinity from the variation of the heat with the solute concentration (that is, from the saturation curve in each titration curve), data obtained in reference [263] for the interaction of rifampicin with some representative concentrations of POPC were reanalyzed and are shown in Figure 9.

The left plot shows the best fit of a model considering binding sites in the POPC membrane, equation (8), but fixing the partition coefficient and enthalpy variation at the intrinsic parameter values obtained from the dependence of $q_{\mathrm{i}}$ with the lipid concentration ($K_{\mathrm{P}}$=3.1×10³ and ${\mathsf{\Delta }H}^{\mathrm{o}}=-$4.3 kJ/mol). A reasonable fit is obtained, especially for low lipid concentrations, providing an estimate of the number of lipids per binding site (${\#}_{\mathrm{L}}$) and binding constant ($K_{\mathrm{b}}$, obtained from equation 10). The right plot shows the best fit obtained when all parameters are allowed to adjust freely. As expected, the quality of the best fit improves, leading to excellent fits. However, the binding constant ($K_{\mathrm{b}}$) is the only parameter that is reasonably consistent for the different lipid concentrations, whereas the estimates of $K_{\mathrm{P}}$ and ${\mathsf{\Delta }H}^{\mathrm{o}}$ are significantly different from those obtained for the interaction with unperturbed membranes and vary widely for the different lipid concentrations. Also, in most cases, an accurate estimate of the parameters could not be obtained, only an upper or lower limit.

The analysis shown in Figure 9 shows that the solute saturation curves obtained by titrating the membrane with increasing solute concentrations do not lead to adequate parameters when analyzed with the assumption of well-defined binding sites in the membrane. Therefore, although this model may be justified in the case of association of large multivalent ions with membranes of opposite charge, where binding is mostly due to electrostatic interactions [66,268], the results may still be described by a simple partition regarded that the contribution from electrostatic interactions is included [268,269]. Moreover, distinct association parameters are obtained when following both approaches (addition of lipid to ligand or vice versa) [268], and more complex models may be required to adequately describe the association of solutes with membranes, including solute-membrane and solute-solute interactions (e.g. [194,270]). The parameters obtained from solute saturation curves should therefore be interpreted with caution, at least for the interaction of small ligands when followed by ITC.

The application of ITC to characterize solute-membrane interactions goes far beyond obtaining the interaction affinity. A major advantage is the full characterization of the variation in the thermodynamic parameters of the system associated with the transfer of the solute from the aqueous media to the membrane ($\mathsf{\Delta }{G}^{\mathrm{o}}$, $\mathsf{\Delta }{H}^{\mathrm{o}}$, and $\mathsf{\Delta }{S}^{\mathrm{o}}$), providing important insight into the type of interactions established, and/or changes in solute/membrane properties [106,271,272]. For this goal, the definition of the reporter affinity parameter is of major importance because this determines the magnitude of $\mathsf{\Delta }{G}^{\mathrm{o}}$, from which $\mathsf{\Delta }{S}^{\mathrm{o}}$ is calculated, equations (30).

$$\mathsf{\Delta }{G}^{\mathrm{o}}=-RT \mathrm{l}\mathrm{n} K=\mathsf{\Delta }{H}^{\mathrm{o}}-T\mathsf{\Delta }{S}^{\mathrm{o}}$$

(30)

From studies of hydrocarbon partition between water and organic solvents, evidence has been obtained for a higher adequacy of partition coefficients defined as volume fraction corrected for the distinct volumes of the species involved to account for variations in the mixing entropy (Flory-Huggins theory) [273,274]. However, controversy remains regarding the most adequate formalism for partition to lipid membranes [271,275]. The impact of the definition of the partition coefficient on $\mathsf{\Delta }{G}^{\mathrm{o}}$ − and correspondingly on $\mathsf{\Delta }{S}^{\mathrm{o}}$ − is given in equation (31) for $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}$ and $K_{\mathrm{P}}$.

$$\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}^{\mathrm{n}/\mathrm{m}}\right)=\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}\right)-RT\ln {\displaystyle \frac{{\rho }_{\mathrm{W}}}{{\rho }_{\mathrm{M}}}}\mathrm{ }\cong \mathrm{ }\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}\right)$$

(31)

$$\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}\right)=\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}\right)-RT\ln \left({\overline{V}}_{\mathrm{L}}\left[W\right]\right)\mathrm{ }\cong \mathrm{ }\mathsf{\Delta }{G}^{\mathrm{o}}\left(K_{\mathrm{P}}\right)-10\mathrm{ }\mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{ }$$

(32)

In the case of a solute with $K_{\mathrm{P}}=1\times {10}^3$ (leading to $\mathsf{\Delta }{G}^{\mathrm{o}}=-18 \mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$), a value of $\mathsf{\Delta }{H}^{\mathrm{o}}=-20 \mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$ would predict $T\mathsf{\Delta }{S}^{\mathrm{o}}=-2 \mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$. If the same membrane affinity is expressed as a molar fraction, $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}=4\times {10}^4$ (leading to $\mathsf{\Delta }{G}^{\mathrm{o}}=-28 \mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$), the same $\mathsf{\Delta }{H}^{\mathrm{o}}$ would predict $T\mathsf{\Delta }{S}^{\mathrm{o}}=8 \mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$. Thus, although in both cases it is concluded that the solute-membrane association is driven by enthalpy, the use of $K_{\mathrm{P}}$ predicts a small decrease in the system entropy, while a significant increase would be predicted if the partition coefficient is reported as $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$. The estimated contribution of entropy variation to the solute-membrane affinity will depend on $\mathsf{\Delta }{G}^{\mathrm{o}}$ and $\mathsf{\Delta }{H}^{\mathrm{o}}$, with $T\mathsf{\Delta }{S}^{\mathrm{o}}$ being always 10 kJ/mol higher (more positive or less negative) if $\mathsf{\Delta }{G}^{\mathrm{o}}$ is calculated from $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$ instead of $K_{\mathrm{P}}$.

Another important application of ITC is in the characterization of changes in ionization due to solute-membrane association. The enthalpy variation directly measured in an ITC titration reflects changes in the system upon transfer of the solute from the aqueous phase to the membrane. It includes changes in the interactions established by the solute with both solvents, as well as other concomitant processes involving heat variations. These include solvent-solvent interactions, which contribute directly to the solutes’ partition affinity, but also other indirect changes in the system components. An important example is the case of changes in the ionization of solute or membrane lipids, which is accompanied by H⁺ release or capture by the pH buffer in the aqueous media. Depending on the enthalpy variation associated with H⁺ dissociation from the buffer ($\mathsf{\Delta }{H}_{\mathrm{b}\mathrm{u}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}}^{\mathrm{o},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}$), the enthalpy variation directly obtained may differ significantly from that associated with solute transfer between the aqueous media and the lipid membrane ($\mathsf{\Delta }{H}_{\mathrm{P}}^{\mathrm{o}}$), the difference being related with the number of H⁺ captured by the solute/lipid (released by the buffer) due to solute transfer from the aqueous media to the membrane, $\mathsf{\Delta }{n}_{{\mathrm{H}}^{+}}$, equation (33).

$$\mathsf{\Delta }{H}^{\mathrm{o}}=\mathsf{\Delta }{H}_{\mathrm{P}}^{\mathrm{o}}+\mathsf{\Delta }{n}_{{\mathrm{H}}^{+}}\mathrm{ }\mathsf{\Delta }{H}_{\mathrm{b}\mathrm{u}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}}^{\mathrm{o},\mathrm{ }\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}$$

(33)

The dependence of $\mathsf{\Delta }{H}^{\mathrm{o}}$ on the pH buffer properties is both an inconvenience and an opportunity. On the one hand, to obtain a full characterization of the thermodynamic parameters associated with solute partition from a single titration it is mandatory to use buffers with small $\mathsf{\Delta }{H}_{\mathrm{b}\mathrm{u}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}}^{\mathrm{o},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}$ such as inorganic phosphate [276,277]. On the other hand, by performing distinct titrations using buffers with distinct ionization enthalpies it is possible to obtain $\mathsf{\Delta }{n}_{{\mathrm{H}}^{+}}$. When the membrane is formed from non-ionizable lipids and/or phosphatidylcholines only, the H⁺ transferred may be directly attributed to the solute and allow characterizing changes in solute ionization upon association with the membrane. This approach has been followed to obtain the charge of several drugs and natural biomolecules when associated with lipid membranes, providing important insights regarding the mechanisms of their bioavailability and bioactivity [3,119,192,265,266,278,279].

ITC can also provide important information regarding the kinetics of solute-membrane interactions. The uptake and release protocol proposed by Heerklotz and co-workers provides semi-quantitative information on the rate of solute translocation between the membrane outer and inner leaflets through the combination of two types of titrations [280]. The uptake titration follows the usual ITC procedure, with liposomes added to the solute in the ITC cell, while in the release titration liposomes previously equilibrated with the solute are added to the aqueous buffer, thus reporting on the heats evolved due to the release of solute from the liposomes. The global analysis of both titrations provides the association affinity, the enthalpy variation and the fraction of lipid accessible to the solute during each titration step. This information is related to the rate of solute equilibration between the two membrane leaflets and also to possible multilamellarity on the liposomes. Additional limitations that prevent the use of this approach to quantitatively obtain the rate of solute translocation include the low rate of titrant addition in the ITC experiment (typically 0.5 µL/s, preventing the characterization of fast translocation processes) and the propagation of the heat evolved in the first titration steps throughout the whole titration in the case of slow kinetics. Nevertheless, this approach allows classification in terms of slow, intermediate or fast permeating solutes [97,119,192,279,281,282,283]. Transformations that occur on the time scale of a few minutes may also be characterized directly through the time dependence of the heat evolved after each titration step. This approach has been used in the case of the interaction of chlorpromazine with liposomes of distinct lipid compositions [97]. After the addition of the liposome aliquot, a large heat variation was observed that reflected the fast association of the solute with the membrane. The full equilibration of the solute with the membrane (that is, the return to baseline) involved however an additional slow step that occurred on the time scale of several minutes and could therefore be quantitatively characterized. After exclusion of contributions from membrane multilamellarity, this slow step was attributed to solute translocation between the outer and inner membrane leaflets. Contributions from other slow processes involved in the full equilibration of charged solutes across lipid membranes cannot however be discarded [284,285].

As a final remark, despite its potentially high utility, the characterization of solute–membrane interactions by ITC is not always straightforward. This difficulty arises primarily from the small enthalpy changes expected for associations that are largely driven by the hydrophobic effect. Although the establishment of electrostatic interactions between the solute and lipid polar groups can produce measurable heat changes, these variations are typically limited to at most a few tens of kJ/mol. Moreover, such interactions are often of low to moderate affinity, resulting in the distribution of the total heat change over many titration steps and, consequently, a further reduction in the heat per injection. To accurately characterize the interactions of small molecules with lipid membranes, high-purity lipids, perfectly matched solvents, and ITC instrumentation in optimal condition are required. A standard operating procedure has recently been published that provides important recommendations, including minimization of solvent mismatch and detailed equipment cleaning and calibration procedures [286].

3.1.3.2. Methods Based on Changes in the Membrane Electrostatic Properties

The membrane-association of charged solutes may be characterized through variations of the membrane surface charge density (${\sigma }_0$), and corresponding surface potential (${\psi }_0$), at increasing solute concentrations. This property may be obtained from the liposome electrophoretic mobility ($\overline{u}$) or zeta potential ($\zeta$) [3,45,140,179,192,213,287,288,289,290], knowing the ionic strength and viscosity of the aqueous media and the distance between the membrane surface and the liposome slipping plane ($d_{\zeta }$), equations (34), [213,288]

$$\zeta ={\displaystyle \frac{2kT}{e}}\mathrm{l}\mathrm{n}{\displaystyle \frac{1+\alpha \mathrm{ }\mathrm{e}\mathrm{x}\mathrm{p}\left(-\kappa {\mathrm{ }d}_{\zeta }\right)}{1-\alpha \mathrm{ }\mathrm{e}\mathrm{x}\mathrm{p}\left(-\kappa {\mathrm{ }d}_{\zeta }\right)}}\mathrm{ };\mathrm{ }\alpha ={\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(e{\mathrm{ }\psi }_0/2\mathrm{ }kT\right)-1}{\mathrm{e}\mathrm{x}\mathrm{p}\left(e{\mathrm{ }\psi }_0/2\mathrm{ }kT\right)+1}}$$

(34)

where $\kappa$ is the inverse of the Debye length (${\lambda }_{\mathrm{D}}$, given by equation (18)).

Evidence for $d_{\zeta }$=2 Å has been obtained for the case of membranes from pure phosphatidylserine and mixtures with phosphatidylcholine, in aqueous media containing monovalent salts [100,213], and this value is commonly used to calculate the surface potential from the measured zeta potential. The membrane electrostatic properties can also be obtained using membrane bound fluorescent sensors (e.g. fluorescein phosphatidylethanolamine or di-8-ANEPPS to report on the membrane surface and dipole potential, respectively), and this has been used to characterize the membrane-affinity of charged peptides [291].

Because the affinity is obtained from changes in the membrane electrostatic properties, to obtain the intrinsic affinity it is always necessary to correct for the effects of the surface potential on the association affinity. This is usually done using the Gouy-Chapman theory (section 1.1).

An important limitation of these methods is the assumption of additivity for the contribution from the solute and lipids, and that their electrostatic properties are independent of the local concentrations in the membrane. Another concern is that they require knowing the average charge of the solute when associated with membrane, a parameter that is not always well known due to changes in the solute ionization equilibria when associated with the membrane. On the other hand, when the membrane affinity is known from independent experiments performed at the same conditions, these methods allow the characterization of the solute average charge when associated with the membrane (e.g. [3]).

3.1.3.3. Methods Based on Changes in the Liposome Size

Thermophoresis is a property based on the directed motion of molecules or supramolecular entities induced by temperature gradients, which depending on the molecular properties and temperature lead to a depletion (thermophobicity) or accumulation (thermophilicity) in the region with a higher temperature [292]. The concentration gradient depends on the diffusion coefficient which allows the characterization of the size of the diffusing entity [293], and in turn may be used to characterize molecular interactions. This methodology is increasingly being used to evaluate ligand-protein binding affinities [294,295,296], and has also been applied to characterize the membrane-binding affinity of proteins and peptides, [297,298] and other bioactive molecules [299,300,301]. To eliminate contributions from solvent flow, the sample is placed in thin capillaries of 100 µm diameter. This methodology requirement is in fact one of the major advantages of this approach, requiring the use of very small amounts of sample. The concentration variations in the heated region of the sample are usually visualized by fluorescence, the method thus requiring the use of fluorescent binding agents, either resulting from the presence of intrinsic fluorescent moieties or after their modification with fluorescent probes. In a typical experiment, the concentration of binding agent is maintained constant, and the dependence of the fluorescence intensity on the ligand concentration is followed overtime during the establishment of the temperature gradient. The saturation fraction of the binding agent with the ligand ($f$) is then calculated from the variation in the fluorescence intensity at a given ligand concentration (${\mathsf{\Delta }c}_{\mathrm{i}}$), which depends on the variation observed in the absence of ligand (${\mathsf{\Delta }c}_{\mathrm{u}\mathrm{n}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$) and at full saturation of the binding agent (${\mathsf{\Delta }c}_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$), equation (35). The binding affinity may then be calculated using the most adequate binding formalism which may include a set of several distinct binding sites with or without binding cooperativity [67,302].

$${\mathsf{\Delta }c}_{\mathrm{i}}=f{\mathsf{\Delta }c}_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}+\left(1-f\right){\mathsf{\Delta }c}_{\mathrm{u}\mathrm{n}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$$

(35)

When applied to the characterization of ligand-membrane binding, this methodology suffers from the same limitations as previously indicated for all methods that rely on membrane saturation with ligand, and thus do not provide the parameters for association to unperturbed membranes.

In the case of large solutes and/or at high local solute concentrations, the association of solute with the membrane increases the liposome size and this may be followed directly by Dynamic Light Scattering (DLS). As indicated above, if the liposome size is followed at increasing solute concentrations the results provide only information regarding the membrane saturation with the solute. However, data can also be obtained for a fixed solute concentration at increasing liposome amounts. In this case, the liposome size is larger at low lipid concentrations (due to the high local concentration of solute) and tends towards the size of pure liposomes as the solute-lipid ratio decreases. This approach has been followed to characterize the association of a lipidic α-amino acid with POPC liposomes, leading to a partition coefficient comparable to that obtained by ITC [192]. In this work, the variation in the liposome size was however very small (less than 2 nm radius for liposomes with a radius of 55 nm) and required the use of high local concentrations of solute (above 10 mol%), highlighting the low sensitivity of the method for small solutes. Although with limited application to characterize quantitatively the solute-membrane affinity, DLS is an important complementary technique, providing information regarding the effects of the solutes on membrane properties (e.g. [47,303,304,305]).

3.1.3.4. Methods Based on Changes in the Lipid Phase Transitions

The temperature of the phase transition of lipid membranes depends on the presence of foreign molecules [306], similarly to solute effects on the freezing temperature depression in homogeneous media, and this may be used to follow the association of solutes with membranes. A major advantage of this method is that it may be used with any solute, without requiring any specific solute property such as charge or fluorescence. The changes in the lipids phase transition may be followed by differential scanning calorimetry (e.g. [240,307,308,309,310]) but also by turbidimetry due to the increase in the refractive index of membranes in the gel phase (e.g. [311,312,313,314]), or by changes in the fluorescence properties of membrane embedded probes such as in the anisotropy of diphenylhexatriene (e.g. [178,179,315,316]).

The formalism to analyze the decrease in the melting temperature is based on the van’t Hoff equation for freezing-point depression [317]. If the solute associates only with the fluid membrane, the decrease in the membrane phase transition temperature is given by equation (36),

$$\mathsf{\Delta }{T}_{\mathrm{m}}=-{\displaystyle \frac{R{\left(T_{\mathrm{m},0}\right)}^2}{\mathsf{\Delta }{H}_{\mathrm{m},0}^{\mathrm{o}}}}{\displaystyle \frac{K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}}{\left[W\right]+K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}\left[L\right]}}\left[S_{\mathrm{W}}\right]$$

(36)

where $T_{\mathrm{m},0}$ and $\Delta H_{\mathrm{m},0}^{\mathrm{o}}$ are the melting temperature and the melting enthalpy of the lipid in the absence of solute, $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$ is the partition coefficient of the solute defined in terms of mole fraction, equation (3), and $\left[W\right]$, $\left[L\right]$ and $\left[S_{\mathrm{W}}\right]$ are respectively the concentration of water, of lipid in the fluid phase (half the total lipid at $T=T_{\mathrm{m}}$) and the concentration of solute in the bulk aqueous phase, all defined as molar fraction.

If the fraction of solute associated with the membrane is negligible − obtained at low lipid concentrations or for solutes with weak membrane-affinity − equation (36) may be simplified to equation (37) and a linear dependence is observed between $\mathsf{\Delta }{T}_{\mathrm{m}}$ and the total solute concentration $\left[S\right]$.

$$\mathsf{\Delta }{T}_{\mathrm{m}}=-{\displaystyle \frac{R{\left(T_{\mathrm{m},0}\right)}^2}{\mathsf{\Delta }{H}_{\mathrm{m},0}^{\mathrm{o}}}}{K}_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}\left[S\right]$$

(37)

The assumption that solutes dissolve only in the melted (fluid) phase does not always hold for lipid membranes. Some lipids, such as phosphatidylcholines—the main lipids in biomembranes—can form an intermediate phase between the gel and fluid states, which can incorporate substantial amounts of solute [318,319]. When the partition coefficient is estimated from the melting temperature depression using equation (36), what is being reported is the difference between solute partitioning into the membrane above and below the melting temperature. If solute binding to the phase below the melting point is significant, the method underestimates the true partition coefficient between the aqueous phase and the fluid membrane [320].

The general equation for the decrease in the melting temperature is given in equation (38) [317].

$$\Delta T=-{\displaystyle \frac{R{\left(T_{\mathrm{m},0}\right)}^2}{\Delta H_{\mathrm{m},0}^{\mathrm{o}}}}\hspace{0.33em}{\displaystyle \frac{K_1^{\mathrm{n}/\mathrm{n}}-K_2^{\mathrm{n}/\mathrm{n}}}{\left[W\right]+\left\{\alpha K_1^{\mathrm{n}/\mathrm{n}}+\left(1-\alpha \right)K_2^{\mathrm{n}/\mathrm{n}}\right\}\left[L\right]}}\hspace{0.33em}\left[S\right]$$

(38)

where $K_{\mathrm{i}}^{\mathrm{n}/\mathrm{n}}$ is the partition coefficient between the aqueous media and membrane phase just above$\left(i=1\right)$, and below $\left(i=2\right)$ the melting temperature, $\alpha$ is the mole fraction of lipid in the melted phase at the temperature evaluated ($\alpha$= 0.5 at the midpoint of the transition).

Using equation (38) to analyze the dependence of the melting temperature on the solute concentration, at different concentrations of lipid, the partition coefficients of several drugs between the aqueous media and DPPC membranes in the fluid and in the intermediate phase (the rippled phase) were obtained [317]. The ratio between the two partition coefficients were equal to two or less, indicating that the partition coefficients for the fluid phase are underestimated by 50% or more if the data is analyzed with the simplified equation (36). The partition into DPPC above and below the main phase transition (from rippled to gel) was characterized by Kawamura et al for a homologous series of benzyl alkanols from phenol to 8-phenyl-1-octanol [321]. It was found that the relative partition into the less fluid rippled phase increased with the length of the alcohol (being negligible for $n\le 4$ and increasing up to 30 % for$n=8$). The authors also analyzed the effect of the solutes on the temperature of the pre-transition (from gel to rippled) and found negligible association of the solutes with the more ordered gel phase. These studies show that variations in the membrane main phase transition may be used to report quantitatively the association of solutes with membranes where the main transition is from gel to fluid phases and when the solute shape is much different from that of the membrane lipids. However, it must be used with caution when the gel-to-fluid transition involves phases with intermediate properties and/or when the solutes have lipid-like structural properties. For this reason, the measurement of changes in the transition temperature is a methodology commonly used, but only for a qualitative assessment of solute-membrane association (e.g. [90,243,322,323,324,325,326,327,328]).

3.1.3.5. Methods Based on Interactions with Solid-Adsorbed Membranes

Nanoplasmonic sensors with application in lipid membrane studies include Surface Plasmon Resonance (SPR) and Quartz Crystal Microbalance (QCM). The details of these methodologies have been reviewed recently including the principles involved and most common applications [329]. A major advantage is that the solute dissociation rate constant from the membrane may be easily characterized, providing information on the interaction kinetics.

Common limitations of these methodologies include the use of membranes attached to solid supports, with possible effects on membrane structural and dynamic properties, and that titrations are performed at increasing concentrations of solutes, where the solute-membrane affinity is obtained from membrane saturation curves. The use of immobilized membranes greatly simplifies the experimental approach and allows the characterization of a large number of solutes with the same membrane preparation. Therefore, although most commonly used in the characterization of solute-macromolecule interactions, these methodologies have also been used to characterize solute-membrane interactions [194,330,331,332,333,334,335,336,337,338,339,340].

The formalism to obtain partition coefficients from the variations in the SPR signal at increasing solute concentrations has been recently derived, [331] and given by equation (39), where ${RU}_{\mathrm{L}}$ and ${RU}_{\mathrm{S}}$ corresponds to the SPR signal in the absence and presence of solute, respectively, $M_{\mathrm{L}}$ and $M_{\mathrm{S}}$ is the molar mass of lipid and solute, and $\sigma$ is the lipid-to-solute ratio at saturation (equivalent to the number of lipids per binding site, ${\#}_{\mathrm{L}}$).

$${\displaystyle \frac{{RU}_{\mathrm{S}}}{{RU}_{\mathrm{L}}}}=\hspace{0.33em}{\displaystyle \frac{\overline{V_{\mathrm{L}}}{K}_{\mathrm{P}}\frac{M_{\mathrm{S}}}{M_{\mathrm{L}}}\left[S_{\mathrm{T}}\right]}{1+\sigma \overline{V_{\mathrm{L}}}{K}_{\mathrm{P}}\left[S_{\mathrm{T}}\right]}}\hspace{0.33em}$$

(39)

The assumptions in the derivation of this equation are: i) a linear dependence of the SPR signal with the mass of the membrane in the sensor (including the membrane and the membrane-associated solute); ii) that the fraction of solute associated with the membrane is negligible ($\left[S_{\mathrm{W}}\right]\cong \left[S_{\mathrm{T}}\right]$); and iii) that the solute affinity for the membrane is maintained ($K_{\mathrm{P}}$ independent on $\left[S_{\mathrm{M}}\right]$) with the volume of membrane available for solute partition decreasing with the amount of solute associated with the membrane, and being given by $\sigma \overline{V_{\mathrm{L}}}\left[L_{\mathrm{T}}\right]$. While assumptions i) and ii) are generally considered in SPR, assumption iii) assumes that the solute interacts with well-defined binding sites in the membrane, and its validity may depend on the solute and membrane. This approach has been followed to characterize the interaction of a protein (${\mathrm{M}}_{\mathrm{W}}$=16.44 kDa), and two peptides (${\mathrm{M}}_{\mathrm{W}}$=4.492 and 0.337 kDa) with POPC membranes, leading to $\sigma$ values of 9.3, 50, and 69, respectively [331]. The observation that the ratio of lipid-to-solute at saturation (that is, the size of the binding site) is inversely related with the solute size, points towards some possible pitfalls in the application of this approach to obtain the solute-membrane affinity parameters. Nevertheless, the partition coefficients obtained are in reasonable agreement with those obtained independently by changes in the fluorescence quantum yield of the solute, highlighting the potential of this methodology at least to obtain relative membrane affinities, with the advantage of providing information on the interaction kinetics as well [341,342,343,344,345,346,347].

3.2. Computational Approaches

3.2.1. Molecular Dynamics Simulations

Computational tools are becoming a standard methodology to evaluate solute-membrane interactions, with molecular dynamics (MD) simulations under NPT conditions being of particular relevance (e.g. [348,349,350,351]). Local partition coefficients, $K\left(z\right)$, are usually calculated from the excess free energy at the $z$ coordinate relative to a reference position, $\mathsf{\Delta }\mathsf{\Delta }{G}^o\left(z\right)$, equation (40).

$$K\left(z\right)=e^{{\displaystyle \frac{-∆∆G^{\mathrm{o}}\left(z\right)}{RT}}}$$

(40)

The solute-membrane affinity parameter in equation (40) was intentionally not defined. It is however necessary to define whether it is $K_{\mathrm{b}}$, or any specific partition coefficient, namely $K_{\mathrm{P}}$ or $K_{\mathrm{P}}^{\mathrm{n}/\mathrm{n}}$, because their quantitative relation with $\mathsf{\Delta }{G}^{\mathrm{o}}$ is different and would lead to the prediction of different solute-membrane affinities, equation (32). To select between the two reporters for solute-membrane affinity it is necessary to understand how is $\mathsf{\Delta }{G}^{\mathrm{o}}$ obtained. In unrestrained simulations, what is characterized directly is the time that the solute stays associated with the membrane or in the water. The fractional time spent in each position along the $z$ coordinate gives directly the probability density function, by calculating a histogram with a discrete number of bins of equal width along the z-coordinate (that contains the membrane and water media). Since the bins have equal widths and $xy$ areas, the local solute densities obtained for each bin can be treated as relative local concentrations in units of moles per unit volume. This suggests that the partition coefficient obtained directly from the relative time spent at each $z$ coordinate is defined as $K_{\mathrm{P}}$, equation 1.

The relation is not so obvious when the free energy profile is obtained from enhanced sampling simulations such as Umbrella Sampling, where the probability density is obtained under the application of a force [352]. Because $\mathsf{\Delta }{G}^{\mathrm{o}}$ is related with the partition coefficient, equation (40), the definition of $K\left(z\right)$ has an impact on $\mathsf{\Delta }{G}^{\mathrm{o}}$, and vice versa. So, before this issue is clarified, the corresponding uncertainty for the predicted solute-membrane affinities obtained from biased MD simulations must be included in the comparison with experimentally measured affinities. For simplicity, in the following analysis it will be assumed that the solute-membrane affinity obtained from MD simulations (unrestrained or biased) is defined as $K_{\mathrm{P}}$.

The local partition coefficients obtained from equation (40) are required to calculate the permeability coefficients from the aqueous media on both sides through the membrane, according to the inhomogeneous solubility-diffusion mechanism [353]. The local partition coefficients cannot however be directly compared with those obtained experimentally which report on the solute association with the whole membrane.

Different approaches have been proposed to calculate the overall partition coefficient ($K$, which following the discussion above is assumed to correspond to $K_{\mathrm{P}}$) from the free energy profile obtained in MD simulations. The simpler approach is to consider the partition between the aqueous media and the equilibrium position in the membrane ($z_{\mathrm{E}\mathrm{q}}$), corresponding to the minimum in the free energy profile [354], equation (41).

$$K_{\mathrm{P}}=e^{-{\displaystyle \frac{\Delta G^{\mathrm{o}}\left(\mathrm{W}\to z_{\mathrm{E}\mathrm{q}}\right)}{RT}}}=e^{-{\displaystyle \frac{\left(\Delta G^{\mathrm{o}}\left(z_{\mathrm{E}\mathrm{q}}\right)-\Delta G^{\mathrm{o}}\left(z_{\mathrm{W}}\right)\right)}{RT}}}$$

(41)

This procedure has been successfully applied to obtain partition coefficients between two homogeneous solvents, [6] where the free energy profile is mostly invariant within each phase. This is not however the case for lipid membranes, with $\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)$ varying significantly with $z$ (the transversal position across the lipid bilayer). The partition coefficient calculated from equation (41) is related to the maximal membrane-affinity, to a specific location in the membrane, and as expected, usually leads to an overestimation of the partition coefficients. Moreover, deviations will depend on the width of the equilibrium position [6,355,356], introducing uncertainty even in semi-quantitative analysis such as in the ordering of distinct solutes. Those limitations have been pointed-out since the first use of this approach [354], but are often overlooked [357].

A more realistic comparison between the partition coefficient obtained from the free energy profile and that obtained experimentally is obtained when the whole membrane is considered. This has first been done by MacCallum and Tieleman [358] for the distribution of hexane between water and a DOPC bilayer, Equation (42).

$$\begin{array}{c}{K}_{\mathrm{P}}={\displaystyle \frac{1}{5}}{\displaystyle {\int }_{-2.5}^{2.5}}K\left(z\right) dz={\displaystyle \frac{1}{5}}{\displaystyle {\int }_{-2.5}^{2.5}}{e}^{-{\displaystyle \frac{\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)}{RT}}}dz\\ \Delta \Delta G^{\mathrm{o}}\left(z\right)=\Delta G^{\mathrm{o}}\left(z\right)-\Delta G^{\mathrm{o}}\left(W\right)\end{array}$$

(42)

where z = 0 was defined for the non-polar center of the membrane, and the limits of the integration correspond to the extreme coordinates of hexane distribution found experimentally (-2.5 nm < z < 2.5 nm) [359]. The integral is divided by the width considered for the bilayer (5 nm) to make the partition coefficient independent of the size of the system considered. This is the method that is being followed in most publications. However, the limits of the integration and the normalization factor have been modified by several authors. In most publications, the integration is only performed for one membrane leaflet, justified by the symmetry of the membrane with respect to both leaflets [193,356,360]. Regardless of the integration being for the whole bilayer or for one leaflet only − leading to equivalent results in the case of symmetric membranes − some authors use different integration limits. In some case this is done to include thicker membranes, with the integration being extended until a flat free energy profile in the aqueous media is obtained [193]. However, in some publications the choice of the integration limits seems rather arbitrary [356,360], and surprising at first, the normalization factor does not always correspond to the width of the integration [356].

Piasentin et al. [361] proposed an alternative equation to calculate the partition coefficient from the free energy profile, Equation (43), where the simulation box is divided into several layers parallel to the membrane surface, $n_{\mathrm{W}}$ corresponding to layers containing water, and $n_{\mathrm{L}}$ to layers containing the lipid membrane. The distinction between both media is done on the basis of the solute free energy profile, with the cut-off set at $\lceil ∆G^{\mathrm{o}}\rceil$ > 0.2 kJ/mol.

$$K_{\mathrm{P}}={\displaystyle \frac{n_{\mathrm{W}}{\displaystyle {\sum }_{i=0}^{n_{\mathrm{L}}}}{e}^{-\frac{∆G\left(z_i\right)}{RT}}}{n_{\mathrm{L}}{\displaystyle {\sum }_{j=0}^{n_{\mathrm{W}}}}{e}^{-\frac{∆G\left(z_j\right)}{RT}}}}$$

(43)

The box layers in equation (43) containing bulk water do not in practice contribute to the calculated $K_{\mathrm{P}}$ because $∆G^{\mathrm{o}}\cong 0$ and $n_{\mathrm{W}}/{\displaystyle {\sum }_{j=0}^{n_{\mathrm{W}}}}{e}^{-\cong 0}\cong$1. Thus, this equation is essentially equivalent to Equation (42) with the integration performed up to $\lceil ∆G\rceil$ > 0.2 kJ/mol and a pre-exponential factor corresponding to the width of the integration.

The free energy profile and corresponding local partition coefficient of a set of hypothetical solutes is shown in Figure 10 to illustrate the effect of the different methods discussed above for the determination of $K_{\mathrm{P}}$. To facilitate the comparison, the minimum of $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)$ is the same for all solutes, only the width of the free energy well and the coordinate of the minimum are changed. The curves shown in black represent a small and very hydrophobic solute whose center of mass (COM) may be located at any of the $z$ coordinates of the lipid non-polar acyl chains. On the other extreme, the curves in blue represent a polar solute that interacts only with the membrane polar region with its COM located at the membrane/water interface. The additional two solutes (magenta and green curves) represent solutes with intermediate polarity where the width and location of the $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)$ minimum is varied.

Because the minimum of $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)$ is the same for all solutes considered, the value obtained for $K_{\mathrm{P}}$ when using equation (41) is also the same, being equal to 4.4×10⁴ (Log $K_{\mathrm{P}}$=4.6). However, a higher membrane affinity would be expected for the solutes with a broader distribution in the membrane, and this is captured by the overall partition coefficients calculated from the integration of $K_{\mathrm{P}}\left(z\right)$. This calls attention to the limitations of using equation (41) for the calculation of the solute-membrane affinity. The exact value obtained for $K_{\mathrm{P}}$ depends on the integration limits and pre-exponential factor. While the integration limits should include all the $z$ coordinates with $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)\ne 0$, the region where the solute free energy is affected by the presence of the lipid membrane, the choice of the pre-exponential factor deserves some discussion.

If the goal is to obtain the average concentration of solute when associated with the membrane, then a pre-exponential factor equal to the width of the integration interval is the most adequate (as performed in equation (42) and (43). This may however be inadequate if the goal is to compare with the partition coefficient obtained experimentally. For a quantitative comparison with experimental data, the pre-exponential factor should be the thickness of the membrane considered in the experimental characterization of the partition coefficient. In MD simulations the membrane is reduced to one dimension, justified by the assumption that the membrane in MD simulations is flat and laterally homogeneous. Therefore, for the quantitative comparison between the partition coefficients obtained by integration of $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)$ and the experimental $K_{\mathrm{P}}$, the pre-exponential factor should be the membrane thickness ($h$) that corresponds to the lipid molar volume considered experimentally (${\overline{V}}_{\mathrm{L}}$). The thickness of each membrane leaflet may be calculated with equation (44),

$$h\left(\mathrm{n}\mathrm{m}\right)={\displaystyle \frac{{\overline{V}}_{\mathrm{L}}\left({\mathrm{d}\mathrm{m}}^3{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}\right)}{a_{\mathrm{L}}\left({\mathrm{n}\mathrm{m}}^2\right)\mathrm{ }{N}_{\mathrm{A}}\left({\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}\right)\mathrm{ }{10}^{-24}\left({{\mathrm{d}\mathrm{m}}^3\mathrm{ }\mathrm{n}\mathrm{m}}^{-3}\right)}}$$

(44)

where $a_{\mathrm{L}}$ is the area per lipid, and $N_{\mathrm{A}}$ is the Avogadro constant. For membranes composed of phosphatidylcholine lipids in the fluid phase and acyl chains with 16 or 18 carbons (e.g. POPC), ${\overline{V}}_{\mathrm{L}}$≅ 0.8 dm³ mol^-1 [57] and $a_{\mathrm{L}}$=0.64 nm² at 30 °C [73]. This leads to $h$ ≅ 2 nm for each membrane leaflet, in very good agreement with the bilayer thickness measured experimentally [73].

The integration limits and pre-exponential factors used in equation (42) that provide consistency with the partition coefficients measured experimentally are therefore integration through the membrane plus interface until the water regions where $\mathsf{\Delta }{G}^{\mathrm{o}}$ is constant (where $\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}=0$) on both sides of the bilayer and the pre-exponential should correspond to the sum of both monolayer thickness calculated from equation (44) considering the molar volumes used to obtain the experimental partition coefficient, equation (45).

$$K_{\mathrm{P}}={\displaystyle \frac{1}{h_1+h_2}}{\displaystyle {\int }_{-{\mathrm{z}}_{\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}=0}}^{{\mathrm{z}}_{\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}=0}}}K\left(z\right) dz$$

(45)

In the case of symmetric membranes, this is equivalent to integrate only on one of the membrane leaflets, from the bilayer midplane to water, equation (46).

$$K_{\mathrm{P}}={\displaystyle \frac{1}{h}}{\displaystyle {\int }_{\mathrm{z}=0}^{{\mathrm{z}}_{\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}=0}}}K\left(z\right) dz$$

(46)

The partition coefficients obtained from equation (46) for the hypothetic solutes are represented in Figure 10. As expected, the partition coefficient obtained for the small and very hydrophobic solute that may be located anywhere in the membrane non-polar region is close to the partition coefficient obtained from the minimum energy, equation (41) − $\mathrm{L}\mathrm{o}\mathrm{g}{K}_{\mathrm{P}}=$ 4.5 and 4.6, respectively. However, for the solutes that are located at a well-defined depth in the membrane, the two estimates of $K_{\mathrm{P}}$ may differ by orders of magnitude, with equation (41) overestimating the solute-membrane affinity.

A note should be given regarding the fraction of lipid accessible to the solute, section 1.2. While all lipids are accessible in MD simulations, limited accessibility may be observed experimentally depending on the rate of solute equilibration and if the liposomes present multilamellarity, leading to an underestimation of the solute-membrane affinity obtained experimentally.

An additional concern when comparing the solute-membrane affinities obtained experimentally and by MD simulations is related to the assumption that the simulated membrane is flat and laterally homogeneous, allowing the quantification of the relative amount of solute in the membrane simply from its probability density at a given $z$ coordinate. In fact, significant lateral heterogeneity and local membrane curvature may be observed when performing biased simulations for very polar or strongly amphiphilic solutes [362,363,364,365,366,367]. When enhanced sampling methods are used that force solutes with polar groups to sample the non-polar region of the membrane, the flexible membrane deforms locally maintaining the interactions between the solute and the lipid polar groups in spite of the low average value of the z coordinate. In turn, when the COM of amphiphilic solutes sample z coordinates which on average correspond to the water medium, the membrane deforms to shield the solute non-polar regions from the water. These effects lead to artificial free energy profiles which may severely influence the solute-membrane affinity calculated.

An example is illustrated in Figure 11 for the case of a homologous series of NBD-Cn amphiphiles with a polar region (the NBD group) and alkyl chains of different length (Cn). The results from MD simulations are taken from reference [362], and compared with experimental results from reference [368]. The results shown in the upper plots were obtained with the coordinate $z$ calculated with respect to the membrane mid-plane in a small cylinder centered in the solute, with this approach being followed to prevent membrane deformation when the solute is at the membrane/water interface. For the less hydrophobic solute (with a 4 carbon alkyl chain, NBD-C4), $K_{\mathrm{P}}$ calculated with equation (46) is in reasonable agreement with that obtained experimentally, while equation (41) leads to an overestimation as expected − ${\mathrm{L}\mathrm{o}\mathrm{g}K}_{\mathrm{P}}=$ 3.0, 3.6, and 4.2 for experimental results and calculated from equation (46) and (41), respectively. The agreement between the membrane affinity calculated from the free energy profile obtained in the MD simulations and that experimentally obtained progressively worsens as the amphiphiles hydrophobicity increases (longer alkyl chains), and for NBD-C16 even equation (46) overestimates $K_{\mathrm{P}}$ by more than 4 orders of magnitude. When the $z$ coordinate is calculated with respect to the mid-plane of the whole membrane, severe membrane deformations are observed and the free energy profiles are strongly dependent on the pulling method (CW – from the membrane center to water, or WC – from water to the membrane center), lower plots. The positive curvature generated in the membrane when following the CW pulling method leads to an increase in the calculated $K_{\mathrm{P}}$, while a membrane thinning and a negative curvature is observed in the WC pulling method that leads to a decrease in the calculated $K_{\mathrm{P}}$.

The example shown in Figure 11 illustrates the importance of preventing membrane deformation when performing biased MD simulations to obtain free energy profiles and partition coefficients. The best approach in this respect would be to perform unrestrained MD simulations. However, adequate sampling at $z$ coordinates with $\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)-\mathsf{\Delta }{G}^{\mathrm{o}}\left(z_{\mathrm{E}\mathrm{q}}\right)\gg 10 RT$ would require unfeasible long simulation times. This approach is therefore limited to solutes with ${\mathrm{L}\mathrm{o}\mathrm{g}K}_{\mathrm{P}}\cong 0 \pm 1$, as performed recently for the interaction of the MRI contrast agent [Gd(DOTA)]⁻ with POPC membranes [369].

In unrestrained simulations, the probability density of the solute at each $z$ coordinate is obtained directly from the fractional time spent by the solute at that coordinate. The local partition coefficient for a specific depth in the membrane (z) may then be obtained directly from the relative solute probability density ($p_{\mathrm{S}}$) or the relative time spent at that $z$ coordinate, equation (47) [353,369,370].

$$K\left(\mathrm{r}\mathrm{e}\mathrm{f}z\right)={\displaystyle \frac{p_{\mathrm{S}}\left(z\right)}{p_{\mathrm{S}}\left(z_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}}={\displaystyle \frac{\mathsf{\Delta }t\left(z\right)}{\mathsf{\Delta }t\left(z_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}}$$

(47)

These local partition coefficients are equivalent to those obtained with equation (40). The excess Gibbs free energy at the $z$ coordinate may also be calculated from the relative solute probability density, equation (48), and all previously discussed approaches may be used to calculate the partition coefficient from the Gibbs free energy profile.

$$\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\left(z\right)=-RTln\left({\displaystyle \frac{p_{\mathrm{S}}\left(z\right)}{p_{\mathrm{S}}\left({\mathrm{z}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}}\right)$$

(48)

The overall partition coefficient may also be calculated directly from the probability density, in analogy with what is done experimentally, equation (49) [369], where $p_{\mathrm{S}}\left(\mathrm{M}\right)$and $p_{\mathrm{S}}\left(\mathrm{W}\right)$ are the probability density of finding the solute in any coordinate within the membrane or in the water, respectively, and $\langle z_{\mathrm{M}}\rangle$ and $\langle z_{\mathrm{W}}\rangle$ correspond to the thickness of the membrane and water layers in the simulation box.

$$K_{\mathrm{P}}={\displaystyle \frac{\frac{p_{\mathrm{S}}\left(\mathrm{M}\right)}{\langle z_{\mathrm{W}}\rangle }}{\frac{p_{\mathrm{S}}\left(\mathrm{W}\right)}{\langle z_{\mathrm{W}}\rangle }}}={\displaystyle \frac{p_{\mathrm{S}}\left(\mathrm{M}\right)\langle z_{\mathrm{W}}\rangle }{p_{\mathrm{S}}\left(\mathrm{W}\right)\langle z_{\mathrm{M}}\rangle }}$$

(49)

The definition of the membrane/water interface may be defined in terms of the dependence of $p_{\mathrm{S}}\left(z\right)$ with the $z$ coordinate, with water corresponding to the region where $p_{\mathrm{S}}\left(z\right)$ is independent on $z$ (that is, where the solute is not affected by the presence of the membrane). In the case of the very polar contrast agent [Gd(DOTA)]⁻ the transition occurs only at $z=$ 4 nm, with a local maximum in $p_{\mathrm{S}}\left(z\right)$ being observed at $z\cong$ 2.8 nm corresponding to the solute adsorbed at the membrane/water interface, and a global maximum in $p_{\mathrm{S}}\left(z\right)$ observed at $z\cong$ 1.6 nm corresponding to the solute inserted in the membrane [369]. Both solute populations should be included in the calculation of the overall $K_{\mathrm{P}}$ because both contribute to the amount of solute associated with the membrane obtained experimentally. However, as discussed above, the division of $p_{\mathrm{S}}\left(\mathrm{M}\right)$ by the width of the total region where the solute is affected by the membrane ($\langle z_{\mathrm{M}}\rangle$) is not in agreement with that assumed experimentally. Following the same reasoning as before, the width considered for the membrane should be that calculated from equation (44). In any case, the different values considered for $\langle z_{\mathrm{M}}\rangle$ introduce only small variations in the calculated $K_{\mathrm{P}}$, up to a factor of 2 in this case.

It should be mentioned that equations (43), (46) and (49) are equivalent, they are simply expressed with regard to distinct solute and system properties [369]. In all cases, to obtain an estimate of the fraction of solute associated with the membrane it is necessary to sum up from $z$=0 until the $z$ coordinate where the solute is no longer influenced by the presence of the membrane ($\left|\mathsf{\Delta }\mathsf{\Delta }{G}^{\mathrm{o}}\right|\cong 0$). However, if the objective is to compare with experimentally obtained partition coefficients, the solute concentration in the membrane should be calculated assuming the same membrane volume as considered experimentally. Therefore, when using equations (43) and (49), the number of membrane windows or its corresponding width should be equal to that assumed experimentally. This is taken into consideration in equation (46) by using the appropriate monolayer thickness $h$ in the pre-exponential.

As a final remark, to obtain a quantitative agreement between the partition coefficient estimated from in silico approaches and those obtained in vitro or in vivo, several fundamental aspects may also need to be addressed. Namely, the need to include contributions from the distinct translational entropy in the two media [371,372,373,374,375], the need to equilibrate the system regarding solute orientation and conformation when associated with the membrane [6,15,360,376,377,378,379]. It is also necessary to consider the changes in the charge distributions and ionization of the solute in the membrane highly heterogeneous medium, [380,381,382,383,384,385] this last aspect being addressed in section 3.2.3.

3.2.2. Implicit Solvation Models

Despite their obvious utility, molecular dynamics simulations with explicitly defined lipid and solvent molecules involve considerable computational cost, namely when calculation of the free energy for the interaction of solutes with the membrane is required, as is the case for partition. For this reason, methods that employ a simplified implicit representation of water and/or lipid molecules have been proposed. Two such approaches are described here, based on the Generalized Born (GB) and Conductor-like Screening (COSMO) solvation models.

Using a mean-field approach to represent the environment, the electrostatic component of solvent-solute interactions can be calculated by numerically solving the Poisson-Boltzmann equation, or slightly less accurately but more efficiently, using analytical approximations such as the generalized Born (GB) method [386]. According to this formalism, the electrostatic contribution $\mathsf{\Delta }{G}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}^{\mathrm{o}}$ to the solvation free energy $\mathsf{\Delta }{G}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}^{\mathrm{o}}$ of a given solute is described by a sum of pairwise terms of interaction between atoms $i$ and $j$, separated by a distance $r_{\mathrm{i}\mathrm{j}}$, with charges $q_{\mathrm{i}}$ and $q_{\mathrm{j}}$, respectively [387].

$$\mathsf{\Delta }{G}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}^{\mathrm{o}}=-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{{ϵ}_{\mathrm{r}}^1}}-{\displaystyle \frac{1}{{ϵ}_{\mathrm{r}}^2}}\right){\sum }_{i,j}{\displaystyle \frac{q_{\mathrm{i}}{q}_{\mathrm{j}}}{f_{\mathrm{i}\mathrm{j}}^{\mathrm{G}\mathrm{B}}}}$$

(50)

In equation (50), ${ϵ}_{\mathrm{r}}^1$ and ${ϵ}_{\mathrm{r}}^2$ are the relative solute and solvent dielectric constants, respectively, the absolute permittivity and additional pre-factors [388] being absorbed into the Coulomb constant used by the force field, and

$$f_{\mathrm{i}\mathrm{j}}^{\mathrm{G}\mathrm{B}}\left(r_{\mathrm{i}\mathrm{j}}\right)=\sqrt{r_{\mathrm{i}\mathrm{j}}^2+R_{\mathrm{i}}{R}_{\mathrm{j}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-r_{\mathrm{i}\mathrm{j}}^2/\left(4R_{\mathrm{i}}{R}_{\mathrm{j}}\right)\right\}}$$

(51)

In the latter equation, $R_{\mathrm{i}}$ and $R_{\mathrm{j}}$ are distance parameters denoted effective Born radii for atoms $i$ and $j$, respectively. In equation,

Tanizaki and Feig [389] adapted this formalism to a “heterogeneous dielectric generalized Born” (HDGB) model of a lipid bilayer. In brief, this extension was made by considering that $ϵ$ in the bilayer phase depends on the transverse location $z$, and considering an apparent local dielectric constant profile, ${ϵ}_{\mathrm{r}}^{\text{'}}\left(z\right)$). Thus, the electrostatic interaction for solute atom pair $ij$ is characterized by a dielectric constant given by

$${ϵ}_{\mathrm{i}\mathrm{j}}={\displaystyle \frac{{ϵ}_{\mathrm{r}}^{\mathrm{\text{'}}}\left(z_{\mathrm{i}}\right)+{ϵ}_{\mathrm{r}}^{\mathrm{\text{'}}}\left(z_{\mathrm{j}}\right)}{2}}$$

(52)

The ${ϵ}_{\mathrm{r}}^{\text{'}}\left(z\right)$ function is computed numerically by solving the Poisson equation as a probe ion moves across an implicit membrane/water model, most efficiently made up of three regions with different dielectric constant values (three-dielectric model; e.g. ${ϵ}_{\mathrm{r}}=2$ in the bilayer core, 7 near the lipid/water interface, and 80 in bulk water).

The electrostatic component is then added to the nonpolar contribution to the solvation free energy ${\mathsf{\Delta }G}_{\mathrm{n}\mathrm{p}}^{\mathrm{o}}$, the latter made up by summing the van der Waals solute-solvent interactions (${\mathsf{\Delta }G}_{\mathrm{v}\mathrm{d}\mathrm{w}}^{\mathrm{o}}$) and the cost of cavity formation (${\mathsf{\Delta }G}_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}^{\mathrm{o}}$), both assumed to depend linearly on the solvent-accessible area (${SA}_{\mathrm{i}}$ for atom $i$):

$${\mathsf{\Delta }G}_{\mathrm{n}\mathrm{p}}^{\mathrm{o}}={\mathsf{\Delta }G}_{\mathrm{v}\mathrm{d}\mathrm{w}}^{\mathrm{o}}+{\mathsf{\Delta }G}_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}^{\mathrm{o}}=\gamma {\sum }_{i}S\left(z_{\mathrm{i}}\right)S{A}_{\mathrm{i}}$$

(53)

In the latter equation, $\gamma$ is an empirical surface tension parameter and $S\left(z\right)$ is a profile function that accounts for the variation of the surface tension along $z$. In the original paper [389], this was obtained by fitting an analytical differentiable function to an MD-derived free energy profile of O₂ insertion, assuming that $SA$ is independent of $z$ for this solute. The HDGB model was then applied to the calculation of the free energy profile of water and amino acid side chain analogs (averaging ${\mathsf{\Delta }G}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}^{\mathrm{o}}$ for different orientations at a given fixed center-of-mass $z$) and also incorporated in Langevin MD simulations of melittin and bacteriorhodopsin.

More recently, an extension of the HDGB was proposed, to account for the possibility of membrane deformation upon insertion of charged solutes [390]. This introduced an additional term ${\mathsf{\Delta }G}_{\mathrm{m}\mathrm{e}\mathrm{m}}^{\mathrm{o}}$, to be added to ${\mathsf{\Delta }G}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}^{\mathrm{o}}$ and ${\mathsf{\Delta }G}_{\mathrm{n}\mathrm{p}}^{\mathrm{o}}$. The term ${\mathsf{\Delta }G}_{\mathrm{m}\mathrm{e}\mathrm{m}}^{\mathrm{o}}$ can be computed by integrating over the membrane plane the contributions from compression, bending and surface tension, characterized by moduli $K_{\mathrm{a}}$, $K_{\mathrm{c}}$ and $\alpha$, respectively [391,392]:

$$\mathsf{\Delta }{G}_{\mathrm{m}\mathrm{e}\mathrm{m}}^{\mathrm{o}}=\iint {\displaystyle \frac{1}{2}}\left\{{\displaystyle \frac{K_{\mathrm{a}}}{d_0}}{u}^2+{\displaystyle \frac{1}{2}}{K}_{\mathrm{c}}{\left({\nabla }^{2}u\right)}^2+{\displaystyle \frac{1}{2}}\alpha K_{\mathrm{c}}{(\nabla u)}^2\right\}d\mathsf{\Omega }$$

(54)

Here, $u$ is the deviation of a membrane leaflet from its unperturbed width ($d_0/2$). Solutes are considered as cylindrical inclusions, with the assumption of independent (and therefore additive) deformation energies in the two membrane leaflets. $u$ is obtained by minimization of the deformation energy, leading to a differential equation which is solved with the leaflet deformation and the contact slope around the contact curve as boundary conditions, together with the assumption that the membrane is flat far from the solute. The extended formalism, termed dynamic HDGB (DHDGB) performs similarly to HDGB for noncharged amino acids, but improves agreement with MD free energy profiles of insertion of charged ones, avoiding the overestimation present in the previous model. Both HDGB and DHDGB were later applied to the calculation of free energy profiles for insertion of a variety of drug solutes, which were conjugated with diffusivity estimates for calculation of permeability coefficients [393].

COSMOmic is an alternative emerging approach for the computation of free energy profiles for solute interaction with anisotropic systems such as micelles and lipid bilayers using a continuum dielectric solvent approach. It was originally introduced by Klamt et al. [394] as an extension of the conductor-like screening model for realistic solvation (COSMO-RS) quantum chemical method. COSMO-RS uses the COSMO solvation formalism to model electrostatic solvent-solute interactions [395], combined with equilibrium statistical thermodynamics, to calculate chemical potentials in solution. While COSMO-RS is useful to estimate partition between immiscible liquids (e.g., [396]), it does not account for the transverse anisotropy of membranes. COSMOmic achieves this by considering the dependence of the chemical potential ${\mu }^{\mathrm{X}}$ of solute $X$ on its transverse position $r^X$ and orientation ${\mathbf{d}}^{\mathrm{X}}$ [394]:

$${\mu }^{\mathrm{X}}\left(r^X,{\mathbf{d}}^{\mathrm{X}}\mathrm{ }\right)={\mu }_{\mathrm{C}\mathrm{R}\mathrm{S}}^{\mathrm{X}}\left(r^X,{\mathbf{d}}^{\mathrm{X}}\mathrm{ }\right)+{\mu }_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{b}}^{\mathrm{X}}\left(r^X,{\mathbf{d}}^{\mathrm{X}}\mathrm{ }\right)+{\mu }_{\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}^{\mathrm{X}}\left(r^X,{\mathbf{d}}^{\mathrm{X}}\mathrm{ }\right)+{\mu }_{\mathsf{\zeta }}^{\mathrm{X}}\left(r^X,{\mathbf{d}}^{\mathrm{X}}\mathrm{ }\right)$$

(55)

In equation 55, the CRS term refers to the location- and orientation-dependent chemical potential calculated from the screening of charge density, which is dominant for polar solutes. The use of COSMO-RS in this calculation requires information on the distribution of the atom types over the bilayer normal, as well as the structure of a representative solute and lipid molecule. In the original COSMOmic article, lipid atom distributions and the structure of a single lipid conformer were retrieved from a sole MD simulation frame of the hydrated bilayer (note that no solute is required, since the structure of the latter is usually obtained from quantum chemical methods such as DFT). More recently, an improved implementation, using lipid atom distributions resulting from the whole MD trajectory and selection of an individual lipid molecule with a solvent accessible surface (SAS) equal to the ensemble average value of this property have been proposed [49].

The second term in the right-hand side, also present in ordinary COSMO-RS, is the so-called combinatorial contribution, which takes into account size and shape ratios of solute and solvent molecules. The third and fourth terms were purposedly introduced to account for the elastic deformation and membrane electrostatic zeta potential contributions to the chemical potential, respectively. In the original COSMOmic paper [394], the elastic term was shown to not lead to significant improvement, and therefore has been neglected in subsequent reports. On the other hand, the zeta potential term was not included in the calculations described in the original COSMOmic paper (which thus ended up being used only for neutral solutes), because of unrealistically strong zeta potentials estimated from MD simulations [394]. This issue was solved in subsequent reports with ionic solutes through the implementation of an empirical model potential with a small number of adjustable parameters [397,398].

Inclusion of bilayer transverse anisotropy is achieved by modeling the bilayer (or micelle) as a stratified liquid, with $n$ layers of typical width 0.1 to 0.2 nm. For each layer, a large number $m$ of possible solute orientations is considered. The probability $p^{\mathrm{X}}$ of finding the solute in layer $i$ (of central position $z_{\mathrm{i}}$) is given by

$$p^{\mathrm{X}}\left(z_{\mathrm{i}}\mathrm{ }\right)={\displaystyle \frac{{\displaystyle {\sum }_{j=1}^m}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta {\mu }^{\mathrm{X}}\left(z_{\mathrm{i}},{\mathbf{d}}_{\mathbf{j}}\right)\right)}{{\displaystyle {\sum }_{k=1}^n}{\displaystyle {\sum }_{j=1}^m}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta {\mu }^{\mathrm{X}}\left(z_{\mathrm{k}},{\mathbf{d}}_{\mathbf{j}}\right)\right)}}$$

(56)

where $\beta ={\left(RT\right)}^{-1}$. It follows that the free energy profile, $\mathsf{\Delta }{G}^{\mathrm{X}}\left(z_{\mathrm{i}}\right)$, is calculated as

$$\mathsf{\Delta }{G}^{\mathrm{X}}\left(z_{\mathrm{i}}\right)=-RT \mathrm{l}\mathrm{n}{\displaystyle \frac{p^{\mathrm{X}}\left(z_{\mathrm{i}}\right)}{p^{\mathrm{X}}\left(z_{\mathrm{n}}\right)}}$$

(57)

In equation (57), the position $z_{\mathrm{n}}$ is that of the outermost layer, corresponding to bulk water. From the free energy profile, the partition coefficient is then calculated using equation (43).

Generally good correlations have been obtained between experimental and COSMOmic estimates for membrane/water partition coefficients [49,394,397,398,399,400,401]. It should be mentioned that quantitative absolute agreement requires subtraction of a constant fitted offset value (ca. 0.32 log units in references [397,398]). Additionally, the partition coefficients of more hydrophilic solutes are commonly underpredicted [401,402].

Still, from the computational point of view, COSMOmic provides a highly efficient way in predicting partition coefficients between water and anisotropic media such as membranes. The partition behavior of each solute (notably including the free energy profile) takes a few seconds, provided that all necessary inputs are in place. The latter typically require MD simulations for the atom distributions and lipid conformer structure (experimental data could be used in principle, though so far this has not been the case) and quantum chemical approaches such as DFT calculations for both lipid and solute. In any case, taking into account that a single MD simulation of a specific bilayer is necessary to predict the partition of any solute in this system, this makes for a considerable gain in speed compared to enhanced-sampling MD. Thus, COSMOmic calculations constitute a very attractive approach for rapid computational estimation of free energy profiles and partition coefficients.

3.2.3. Inclusion of Changes in Solute Ionization upon Membrane-Association

The computational methods described in sections 2.2.1 and 2.2.2 do not take into consideration possible changes in the solute ionization when partitioning from the aqueous media to the lipid membrane. This is a major limitation given the common deviation observed in the solute pKa upon membrane-association (see section 2.1.2.2). To overcome this limitation, the free energy profile of the most relevant ionization forms of the solute may be obtained using independent simulations. The partition coefficient of each species may then be obtained from the respective free energy profile, following the formalisms discussed in section 2.2.1. If the ionization constants in the aqueous media ($K_{\mathrm{a}}^{\mathrm{W}})$are known, these simulations provide all the necessary parameters to describe the interaction of the solute with the lipid bilayer at a given pH value. In the case of a single ionizable group, the full characterization of solute interaction with the membrane involves 4 species (the protonated and unprotonated species in the water and associated with the membrane, ${SH}_{\mathrm{W}}$, $S_{\mathrm{W}}$, ${SH}_{\mathrm{M}}$, and $S_{\mathrm{M}}$, respectively) and 4 equilibrium constants ($K_{\mathrm{P}}^{\mathrm{S}\mathrm{H}}$, $K_{\mathrm{P}}^{\mathrm{S}}$, $K_{\mathrm{a}}^{\mathrm{W}}$, and $K_{\mathrm{a}}^{\mathrm{M}}$), Figure 6. The simulations provide estimates for $K_{\mathrm{P}}^{\mathrm{S}\mathrm{H}}$and $K_{\mathrm{P}}^{\mathrm{S}}$, and if $K_{\mathrm{a}}^{\mathrm{W}}$ is known, the thermodynamic cycle allows calculation of $K_{\mathrm{a}}^{\mathrm{M}}$, equation (26).

The free energy profiles of the distinct species may also allow obtaining the ionization equilibrium constant at distinct depths in the membrane ($K_{\mathrm{a}}^{\mathrm{z}}$). For this goal it is first necessary to shift the energy profiles of the distinct ionization forms to match their free energy difference in the aqueous media predicted from $K_{\mathrm{a}}^{\mathrm{W}}$. This procedure also allows the calculation of the ionization constant at the equilibrium position in the membrane ($K_{\mathrm{a}}^{{\mathrm{z}}_{\mathrm{E}\mathrm{q}}}=K_{\mathrm{a}}^{\mathrm{M}})$. The overall free energy profile at a given pH in the aqueous media may also be calculated, corresponding to the minimum free energy path, from which the overall partition coefficient may then be calculated and compared with that obtained experimentally at a given pH value. The overall free energy profile can also be used to calculate the rate of solute permeation through the lipid membrane, but this is beyond the scope of this revision (interested readers may refer to the recent reviews [383,403]).

The approach described above has been followed to obtain the relative stability of the neutral and ionized species of several drugs at different depths in the membrane and used to calculate the permeability coefficient at different pH values [17,384,404]. The free energy profiles of the distinct ionization forms have also been obtained for some additional solutes, including drugs [3,6,8,405,406], other bioactive molecules [17,407], and fluorescent probes [193,212,408], but the distinct free energy profiles were mostly used to obtain the relative affinity of the distinct forms for the membrane, and their relative permeability coefficients.

An alternative approach to incorporate solute ionization when characterizing its interaction with membranes is following the constant pH methodology, where changes in the solute ionization are allowed during the MD simulations [384,409,410,411]. Constant-pH molecular dynamics (CpHMD), refers to a group of methods that simultaneously determines the charge states of ionizable groups during conformational dynamics. While conventional molecular dynamics simulations generally assume fixed protonation states, a constant-pH methodology addresses the dynamic shifts in protonation states during the simulation. The CpHMD methods are categorized as discrete or continuous, based on how charge states are represented and sampled during the simulations. In discrete CpHMD [409], a standard MD simulation is periodically interrupted to add a step of charge state evaluation of each ionizable group using Monte Carlo methods, which determines if the randomly generated new charge states should be accepted or rejected. Then, MD simulation will continue with the new charge states if accepted and/or old charge states if rejected. In continuous CpHMD [412,413,414,415], an extended Hamiltonian of the system is used where the added degree of freedom is based on the λ-dynamics approach for free energy calculations, where each ionizable site is assigned to a λ value between 0 (protonated form) and 1 (deprotonated form), that propagates simultaneously with spatial coordinates. This formalism avoids the large energy jumps of discrete methods by allowing a progressive rearrangement of the local environment [416].

The CpHMD methods have been applied to the mechanistic understanding of variety of biomolecular systems. Examples include the characterization of the protonation equilibria of titratable membrane lipids, solutes at the membrane/water interface [382,417,418,419,420,421,422] and during membrane insertion and permeation [383,384,423].

Yue et al. characterized propranolol ionization equilibria and membrane permeation using the CpHMD methods and those involving the combination of distinct fixed-ionization MD simulations, with a good agreement being obtained for the lowest free energy path from both methods [384]. In that work protonation and conformation equilibria were dynamically coupled, to study the permeation of the ionizable drug through a POPC lipid bilayer. It was found that propranolol migrates into the bilayer in the charged form and deprotonates at the hydrophobic boundary, and that by involving dynamic protonation the permeation coefficient was predicted with high accuracy. This work shows that the correct determination of the location-dependent ionization is required for the correct determination of membrane partition and membrane permeation.

The use of constant-pH MD simulations is still limited despite the developments for the application in different simulation software [424,425] and force fields [385,425,426]. This is mainly due to their higher level of system preparation and computational complexity [427,428]. Fortunately, nowadays with tools such as the pHbuilder for the automatic preparation of constant-pH MD simulations in GROMACS [429], it is expected that constant pH simulations become easier to set up and run, thereby making them more accessible to the scientific community.

4. Conclusions

There is a large and varied ensemble of methodologies to quantify the association affinity of small molecules, such as drugs, for model lipid membranes. In this work, the advantages and limitations of experimental and computational approaches have been systematically analyzed and the diversity of affinity parameters reported in the literature has been rationalized. Particular attention has been given to the assumptions, formalisms, and experimental conditions underlying each methodology, as these factors ultimately determine the validity and comparability of the obtained results. By explicitly addressing these aspects, this review aims to clarify the origins of apparent discrepancies between methods and to define the conditions under which meaningful consensus parameters may be established.

Beyond a critical comparison of methodologies, this integrative perspective contributes to the establishment of robust biochemical and biophysical foundations for solute–membrane interactions. Such foundations are essential for relating experimentally accessible affinity parameters to the molecular properties of small molecules, including size, polarity, charge, and chemical functionality. In this context, the combined use of carefully designed experiments and appropriately parameterized computational models emerges as a promising strategy to bridge molecular-scale descriptions and macroscopic observables. Ultimately, the framework presented here provides a basis for the development of quantitative and predictive relationships aimed at forecasting membrane affinity from molecular structure, thereby supporting both fundamental studies and applied efforts in drug design and membrane biophysics.