2.1. Partition to POPC LUVs
As explained in detail in the
Section 3.5, the partition of H33342 to POPC LUVs was characterized from the variation in the fluorescence intensity and from the variation in the fluorescence anisotropy. At pH ≥ 6, association with the membrane leads to a large increase in the fluorescence intensity, allowing to directly obtain the amount of H33342 in both media. However, for pH ≤ 6 the major species in the aqueous media has a large fluorescence quantum yield, and only a small variation is observed upon association with the membrane. Under these conditions, the quantification of the partition coefficient from the variation of the fluorescence intensity has a large uncertainty. On the other hand, interaction with the membrane imposes constraints in H33342 mobility, leading to a significant increase in the fluorescence anisotropy. The variation observed in the fluorescence anisotropy at a given volume of membrane is however not a direct measure of the amount of H33342 in each media. Instead, it reflects the relative fluorescence intensity from the probe molecules in both media. To validate the methodology based on the fluorescence anisotropy, both approaches were followed at pH=6. Typical results obtained at this pH are shown in
Figure 2. The fluorescence anisotropy depends much strongly on the POPC concentration because the anisotropy is dominated by the H33342 associated with the membrane. That is, the ratio between the fluorescence intensity from the probe in the membrane and that from the probe in the aqueous phase (
) is much larger than 1. This parameter is larger at shorter wavelengths because the spectrum of probe associated with the membrane is blue-shifted when compared to that for the probe in the aqueous media. Thus, the fluorescence anisotropy increases more abruptly with the lipid concentration for fluorescence detection at lower wavelengths (
Figure 2b). When the variation in the fluorescence anisotropy is weighted by
(
=13 in the range 450-470 nm, 8 in the range 470-490 nm, and 5 in the range 490-510 nm), the obtained partition coefficient is essentially independent on the detection wavelength. Statistically equivalent values are also obtained from the analysis of the fluorescence intensity (equation (2)) and the fluorescence anisotropy (equation (3)).
Typical results obtained at all pH values studied are shown in
Figure A1 and
Figure A2 of
Appendix A, and the average values of
KP obtained are provided in
Table 1. The partition coefficient obtained at neutral pH is in excellent agreement with the values reported in literature for partition to DMPC membranes in the liquid disordered phase [
5], and to the lipid bilayer of inside-out vesicles of
Lactococcus lactis [
33]. A maximal value of
KP is observed at pH=10, in agreement with the higher abundance of the neutral form, predicted in literature studies [
7,
34,
35] and also by the MarvinSketch software [
36]. The affinity for the POPC bilayers is essentially unaffected by the pH when the latter varies between 7 and 11, although it is predicted that at pH=7 the most abundant species is the monocation. A significant decrease is observed for acidic pH values. However, this variation is still relatively modest, with only one order of magnitude decrease between pH=10 and 3. This result is unexpected, given that the most abundant species is +3 at pH=3.
The partition of H33342 to the POPC membranes was also characterized by ITC to obtain the enthalpy of interaction. The enthalpy variation obtained directly from the thermogram includes contributions from the interactions established between the solute and the medium where it is dissolved (aqueous medium or membrane), as well as from changes in solute ionization upon partition, the latter also including contributions from changes in the ionization of the buffer [
21,
38,
39]. The phosphate buffer used in the ITC experiments has a very low ionization enthalpy [
40], thus with a negligible contribution to the observed enthalpy variation. Therefore, the overall value obtained for the enthalpy variation reflects only changes in H33342, interactions established with the aqueous medium and the lipid membrane, and possible changes in ionization upon partition.
Figure 3 shows typical titration curves obtained at different pH values. The partition coefficients obtained are in good agreement with those obtained from changes in H33342 fluorescence properties, although somewhat lower. This may be due to the higher concentration used in the ITC experiments, to incomplete equilibration of H33342 between the two bilayer leaflets during the titration [
20], or to a small preference for membranes with a higher curvature [
41]. This was not further explored because the enthalpy variation is the most relevant parameter from those experiments.
It is observed that the interaction of H33342 with the POPC membranes shifts from being stabilized by enthalpy at acidic pH values (ΔHo = -54 ± 7 kJ/mol and TΔSo = -38 kJ/mol at pH=5.3), to equal contribution from enthalpy and entropy at alkaline pH values (ΔHo = -9 ± 2 kJ/mol and TΔSo = 9 kJ/mol at pH=8.2). This suggests that the association of the neutral form of H33342 with POPC membranes is mostly stabilized by the hydrophobic effect, while stabilizing electrostatic interactions are established between the cationic forms of H33342 and the membrane, even in the case of membranes containing only zwitterionic lipids such as POPC. To further elucidate on the nature of the interactions established, computational studies were performed with H33342 in distinct protonation states. The results obtained are presented in the next sections.
2.3. Unrestrained MD Simulations
As described in the methods section, the single most stable structures obtained from the quantum chemical calculations for H(0), H(+1) and H(+2) were parameterized for MD simulations and simulated in triplicate for 200 ns, in the presence of a POPC bilayer.
Figure 5 depicts final snapshots in each simulation (run 1 and 2 refers to the system with probe molecules initially placed in the water phase, while in run 3 the probe molecules were initially placed in the bilayer midplane). Even though they only depict a single configuration state, some general features can already be observed. All molecules, irrespective of the ionization state, interact with the lipid bilayer. H(+2) molecules often adopt an orientation perpendicular to the membrane plane, where they insert the region of the molecule bearing no protonated sites (ring 3, ethoxy substituent) into the bilayer, while keeping the protonated region (ring 1, piperazine) in the upper headgroup region or even protruding into the water. This tendency is progressively lost for lower ionization states.
Figure 5 reveals that two molecules in run 1 of H(+1) aggregate (average distance between centers of mass is 0.90 nm in the last 150 ns of simulation, with minimal distance of 0.36 nm; the dimer is clearly visible in the lower leaflet). Aggregation markedly affects probe interactions with the bilayer, and because we believe that this is a result of the relatively high probe:lipid ratio (1:32) necessary for simulating an adequate number of solute molecules, we excluded these two molecules from subsequent analyses. Therefore, whereas in the H(0) and H(+2) analyses twelve molecules were used for averaging, only ten were considered for H(+1).
The time evolution of the center of mass vertical position (
z) for the four individual molecules in each run (
Figure S1) shows that
z converges rapidly in the simulations that start with the molecules in the bilayer midplane (
z = 0, replicate 3), mostly to values in the 1 nm < |
z| < 1.5 nm range. However, in replicates 1 and 2, where H33342 molecules were initially placed in the water/headgroup interface, insertion of some molecules is still occurring after 50 ns. Some molecules in these systems appear to not actually having inserted by the end of the simulations, residing at |
z| ≈ 2 nm. Despite these different behaviors, all molecules (except the aggregated pair in run 1 of H(+1), as previously explained) were included for subsequent analyses, after discarding the first 100 ns of each simulation.
The mass density of H33342 along the
z direction is shown in
Figure 6 for the three ionization states. Overall, no pronounced differences are observed, with the preferred location being around 1.5 nm from the bilayer center for all forms. This was already suggested by the similar position of the center of mass (
Figure S1). Curiously, H(+1) shows a more internal peak than both H(0) and H(+2), but it has a wider density profile.
Despite this similarity, the molecules clearly change their average orientation in the bilayer in response to varying their degree of protonation. We looked at orientation differences in two ways: by calculating average positions of different structural groups (
Figure 7) and angular distributions for the H33342 long axis (defined in
Figure 1a) tilt relative to the bilayer normal (
Figure 8).
Figure 7 shows average transverse positions of groups and atoms of H33342, in comparison with those determined for selected atoms of POPC in a probe-free system (as shown in
Figure S2, the presence of H33342 leads to very slight changes in POPC atom locations, <1.5% for those included in
Figure 7). Globally, it illustrates that the probe has similar vertical localization in the membrane for all ionization states, corresponding to the upper acyl chain/carbonyl/glycerol region of the bilayer. The center of mass position appears to show a very slight dependence on the ionization state, and slightly more external locations are observed for progressively increased overall charge. However, the difference is ultimately very small and well within the estimated uncertainty, in accordance with the density profile analysis above. Still, there is a larger difference in positions of particular H33342 atoms/groups, which, while not entirely significant at the 95% confidence level considered in the figure error bars, appears to be sufficiently consistent overall to allow the following considerations. The uncharged species prefers to adopt an almost horizontal orientation, in the plane of the bilayer. The average positions of the different atoms of H(0) lie in a narrow bilayer region (between ~1.2 and 1.5 nm) for this ionization state, with the piperazine ring (N2) slightly more internal than other parts of the molecule. This contrasts with the behavior observed for H(+2). For this species, the most internally located part is the ethoxy group, with the O41 atom located on average at
z = (0.89 ± 0.21) nm. Conversely, the opposite end of the molecule (the piperazine ring) has a very external transverse location in the headgroup region (
z = (1.94 ± 0.19) nm), implying that in this ionization state, H33342 spans ~1 nm in the direction normal to the bilayer plane. An intermediate scenario is observed for H(+1), with specific atoms/groups located on average between
z = (1.12 ± 0.17) nm (for O41) and (1.71 ± 0.54) nm (for N2).
The long axis tilt distributions (
Figure 8) confirm these patterns. All distributions are very wide, illustrating the flexibility in the orientation of H33342 in the bilayer as well as distinct behaviors of different molecules (some of which actually do not insert properly in the course of the simulation, as pointed out above). Still, the widest distribution is clearly that corresponding to H(0), which displays significant probability density both for θ < 90º (piperazine pointing towards water) and θ > 90° (piperazine pointing towards center of bilayer), though the latter dominates slightly, and the average tilt angle for this state is <θ> = 100°. The opposite behavior is verified for H(+2). This form favors orientation of piperazine towards water and insertion of the opposing end, with <θ> = 52°. Again, H(+1) displays intermediate behavior (<θ> = 74°).
To obtain further details on the interactions established between H33342 and the lipid molecules, and how they may depend on the probe ionization state, we studied the fractional frequency of hydrogen bonding between NH donor groups of H34442 and oxygen acceptor atoms of POPC (
Figure 9). Globally, it is observed that H bonding to glycerol O atoms (O14 and O35 for
sn-2 and
sn-1 lipid acyl chains, respectively) is low in all cases (<2% of total for the three forms), and only bonding to POPC O atoms from the phosphate (O7, O9, O10, and O11, mainly the three latter atoms) and carbonyl (O16 and O37 for
sn-2 and
sn-1 lipid acyl chains, respectively) groups need to be considered.
It must be noted that the higher the probe charge, the larger the number of NH donors available for H bonding. H(0) has two NH groups, one in each of the benzimidazole rings. H(+1) has one more NH group, in the piperazine ring, while H(+2) has that group and an additional one, in the benzimidazole ring closest to the piperazine. The fractional frequency of H-bonding summed across all POPC acceptor atoms reflects this obvious fact, and increases in the order 0.98 (H(0)) < 1.31 (H(+1)) < 2.53 (H(+2)). When these values are divided by the actual number of NH donor groups, similar values are obtained for H(0) and H(+1) (0.49 and 0.44, respectively), whereas the one for H(+2) is slightly higher (0.63). Overall, these results indicate that H33342 NH groups are roughly involved in H bonding to lipid O acceptors 40-60% of the time. The fractions of bonding to the phosphate (39%, 25% and 41% for H(0), H(+1) and H(+2), respectively) or carbonyl POPC O atoms (59%, 74% and 57% for H(0), H(+1) and H(+2), respectively) are similar for the three forms. The fact that maximal bonding to the phosphate occurs for H(+2) (41%) reflects perhaps the more external position of piperazine and benzimidazole ring 1, all with NH donors. However, in general, the transverse location of the NH donor groups is around
z ~1.5 nm (
Figure 7), enabling efficient H bonding to carbonyl POPC atoms, which are the dominant interactions.
The higher number of hydrogen bonds from H33342 NH to POPC acceptor atoms, as H33342 protonation increases, may be correlated with the increase observed in the interaction enthalpy as the pH is decreased (
Figure 3). Those favourable interactions can also explain the relatively small variation observed in the affinity for the membrane, despite the large decrease in the hydrophobicity of the protonated charged forms of H33342.
2.4. Free Energy Profiles from Umbrella Sampling MD
From umbrella sampling simulations, free energy profiles were calculated for the three considered H33442 forms, as described in detail in the Methods section. The potential of mean force (PMF) curves are shown in
Figure 10, whereas
Figures S3-S5 depict uncertainty estimates and convergence analysis (PMF profiles obtained using different regions of the 0 <
t < 120 ns simulation time range, following the procedure of reference [
42]) for the different curves.
As expected for amphiphilic solutes [
42,
43,
44,
45,
46], the free energy is reduced when they interact with the lipid bilayer, coming from the aqueous medium. Energy minima are observed for the three ionization states at depths close to those obtained from the unrestrained simulations and shown in
Figure 7, again similar for all forms (slightly more interior for the neutral form). The minimal ordinate values are also similar for the three states. Whereas the value calculated for H(+2) form is lower than those obtained for the other ones by ~5-6 kJmol
-1, this difference is smaller than the uncertainty of approximately ±10 kJmol
−1 indicated from bootstrapping analysis (Figs. S3E-S5E). In the absence of clearcut distinction in the regions around the equilibrium locations, the most striking difference in the three free energy profiles occurs for the regions nearest the center of the bilayer (
z ≈ 0). Neutral H33442 displays a 24 kJmol
-1 energy barrier for crossing from the free energy minimum to a maximum located at
z = 0.26 nm, before decreasing by 3.5 kJmol
-1 between the latter location and
z = 0. In recent studies by our group, free energy plateaus or slight decreases near the bilayer center have been obtained for neutral solutes such as cholesterol [
43], neutral rhodamine dyes [
44] or undissociated fatty acids [
45]. This feature is not observed in the profiles of the charged H(+1) and H(+2), for which maxima are observed at
z = 0, with free energy barriers (of 40 and 78 kJmol
-1, respectively) that increase with the charge of the solute.
Following the procedure used in recent studies [
44,
47,
48,
49,
50,
51], these calculated free energy profiles may be used to obtain computational estimates of the partition coefficient:
Because of different underlying reference states in the experimental and calculated
KP, as well as arbitrariness in the location of the bulk water/lipid interphase in the simulations, here considered to be
a = 4.0 nm, no attempt is made to calculate absolute
KP values. Still, we can compare relative values of computational estimates of
KP for the three simulated forms. From the PMF profiles, the lowest
KP value was obtained for the H(+1) state, below the corresponding ones for H(0) and H(+2) by 0.3 and 1.0 log units, respectively. This order does not fully coincide with the one resulting from the experimental determinations (
Table 1), which point to a steady ~2-fold (or 0.4 log units) increase in
KP as pH is increased from 5 to 9, with dominant protonation state varying between +2 and 0 [
7,
34,
35] Still, both experimental and computational values agree in the fact that
KP shows low sensitivity across a wide range of pH values.