Medicine Chemistry
Chemistry,1970-2014-Chemistry-2011.jsonl.gz.gz
Chemistry Biology
45da14143fd68e9c0ec9ef318032af97323a4e29
Thermodynamic study of benzocaine insertion into different lipid bilayers

Thermodynamic study of benzocaine insertion into different lipid bilayers J. J. López Cascales,1,a) S. D. Oliveira Costa,1 and R. D. Porasso2 1Universidad Politécnica de Cartagena, Grupo de Bioinformática y Macromoléculas (BIOMAC), Aulario II, Campus de Alfonso XIII, 30203 Cartagena, Murcia, Spain 2Instituto de Matemática Aplicada San Luis (IMASL) Departamento de Física, Universidad Nacional de San Luis/CONICET, D5700HHW, San Luis, Argentina


I. INTRODUCTION
The existence of certain drugs with anesthetic properties has been known for over a century since the pioneering studies of Freud 1 with cocaine. Some of these anesthetic molecules are widely used in clinical anesthesiology for pain relief and during surgical operations. However, the molecular mechanism of these drugs is still a subject of debate, although there is a broad consensus on the role played by sodium channels. [2][3][4][5] There is also a tendency to concentrate studies on specific anesthetic-protein interactions, with less attention paid to the effects of the anesthetic on the protein environment, e.g., the lipids that surround the protein membranes.
In this sense, several authors 6,7 have suggested a way of action based on a ligand-receptor binding mechanism, in which the anesthetics could bind to the protein forming the sodium channels, which would block ion transduction. However, alterations in the lipid constituents of the membrane could also influence the microstructure and function of these ion channels. The requirement of proteins to bind with specific lipids in annular clefts or interactions with certain phospholipid head groups means that any displacement of these lipids due to the presence of anesthetics could also alter the protein functions. In this regard, studies based on gramicidine 8 have demonstrated that an anesthetic dissolved into the membrane induces changes in the overall dynamics of the protein, altering its function in the process. In addition, a more complex and sophisticated understanding at molecular level is required since the cellular membrane is a raft of lipids a) Author to whom correspondence should be addressed. Electronic mail: javier.lopez@upct.es.
of different compositions that incorporate certain membrane proteins. In this context, even when the argument for a lipid site action might be revived based on recent developments and research, 9 almost a century ago Meyer and Overton 10, 11 associated the potency of anesthetics with their solubility in a lipid bilayer, at a time in which the lipid bilayer structure had not even been described. It has been suggested that the solubility of anaesthetics in a lipid bilayer may cause a disturbance in membrane fluidity 12-17 that could affect the protein structure of sodium channel formation, and as a consequence of this perturbance, a conformational change in the ion channel could take place, resulting in the blockage of ion transduction. 18 Thus, rather than considering the cell membrane as a homogeneous entity, it should be considered as a structure with different domains formed by different components, in which lipid rafts of different composition are spread over the membrane, and certain protein functions might be moderated by specific lipids. 19 Thus, perturbations in a part of the membrane bilayer could affect other specific regions of membrane proteins. Therefore, studies on alterations in membrane protein signaling should look at complex processes rather than a mere alteration in the binding site of a membrane protein.
In summary, it would seem essential to characterize in the first instance the thermodynamic parameters associated with the insertion process of a local anesthetic into the bilayer according to its lipid composition, with the aim of providing insight into this process at molecular level. For this purpose we used benzocaine, a local anesthetic widely used in topical applications, which exists only in its neutral form under physiological conditions. 14, 20 The lipid bilayer was formed by DPPC (dipalmitoylphosphatidylcholine) and DPPS (dipalmitoylphosphatidylserine), where DPPC and DPPS correspond to an uncharged (neutral) and charged (negative) phospholipid in physiological conditions, respectively. Different lipid bilayers were generated by varying the lipid percentage that constitutes the membrane.
To provide insight into the process of insertion of the benzocaine into a lipid bilayer, molecular dynamics simulation was used. Based on previous studies of molecular insertions into lipid bilayers, 21,22 the study of the thermodynamic insertion of benzocaine into a lipid bilayer was carried out with atomic resolution. Information regarding the spontaneity of the insertion ( G) and the thermodynamic parameters involved in this process ( H and S) is provided.

A. Setting up the simulation conditions
The lipid bilayer was assembled from 72 lipids (36 lipids per leaflet). Further information on how the three-dimensional computational box were generated, may be consulted in previous works. [23][24][25] Assuming that asymmetric models of lipid bilayers 25 mimics much better a biological plasma membrane than the models of symmetric lipid bilayer; however, with the aim of studying the thermodynamic properties involved in the benzocaine insertion into binary lipid bilayers of different compositions, symmetric bilayers of different DPPC/DPPS ratios were simulated. This substitution was made on both sides of the bilayer so that the bilayer symmetry was maintained in all simulations. Thus, the parameter χ was defined as follows: where χ represents the fraction of DPPS with respect the total number of lipids that form the membrane. In this respect, χ = 0 corresponds to a pure membrane of DPPC, and χ = 1 represents a membrane constituted only by DPPS. Thus, the following membranes were assembled as a function of the DPPC:DPPS composition: 72:0, 60:12, 48:24, 24:48, and 0:72, corresponding to χ = 0, 0.2, 0.3, 0.7, and 1, respectively. The system corresponding to χ = 0 contained 2517 molecules of the single point charge (SPC) water model 26 in the simulation box. In the case of χ = 0, to balance the negative charge associated with each DPPS lipid, the corresponding sodium ions Na + were introduced by substitution of water molecules, to balance the charge associated to the DPPS. Thus, 12, 24, 48, and 72 sodium ions were introduced into the system to balance the number of DPPS in the bilayer. As regards the anesthetic, two benzocaine molecules were introduced into the system at different positions along the z-axis perpendicular to the lipid surface. A scheme of the benzocaine molecule used in this work is shown in Figure 1. The force field parameters and electrostatic charge distribution of the DPPC and DPPS were the same as those used by Egberts et al. 27 and Cascales et al. 28 in previous simulations of pure bilayers of DPPC and DPPS in their liquid crystalline state. The charge distribution and force field parameters for benzocaine were taken from a previous article. 22 The MD simulations were carried out using the GROMACS-3.3.3 package. The system was coupled to a weak temperature and pressure coupling bath algorithm proposed by Berendsen,29 with coupling constants of 0.1 ps and 1 ps for the temperature and pressure, respectively. A time step of 2 fs was used as integrator time to generate the trajectories. The electrostatic interactions were evaluated using the fast particle-mesh Ewald method 30,31 and the Lennard-Jones potential was used for the long range interactions, considering a cut-off of 1 nm. Finally, all the bond lengths were constrained using LINCS. 32 All the simulations were performed at 350 K and 1 atm. This temperature was chosen above the transition temperature of 314 K and 326 K for the DPPC and DPPS, respectively, 33,34 to ensure that the bilayers were always in their liquid crystalline phase. In addition, this temperature provided a window of temperatures to avoid any change in the lipid phase, within the range of temperatures (between 345 K and 355 K) used to calculate the thermodynamic properties of the system.

B. Free energy profile
The umbrella sampling method 35 was used to evaluate the free energy profile associated with benzocaine insertion into a lipid bilayer. The relationship between the free energy profile and the potential mean force is depicted in the following equation: where G b (z) represents the free energy profile along the zaxis perpendicular to the bilayer, R the constant of the gases, T the temperature, C eq b (z) the concentration profile of the benzocaine along the z-axis, C * b the concentration of benzocaine in the solution (aqueous solution in our case), and P MF (z) the potential of mean force computed with respect to a reaction coordinate equal to the z coordinate of the center of mass of the benzocaine molecule, using the umbrella algorithm and W H AM method (weight histogram analysis method). 36 Since an accurate profile of C eq b (z) could not be obtained during the simulation time, a biasing potential must be applied to the benzocaine. In this respect, 36 independent simulations of 50 ns each were performed, for each of the five systems studied in this work, which amount a total of 9000 ns of simulation time. In each of the 36 simulations for any of the systems studied, the benzocaine was restrained at a given depth in the bilayer by a harmonic potential on the z-coordinate perpendicular to the membrane plane, leaving it to move freely on the membrane plane. The force constant of this harmonic spring was 3000 kJ mol −1 nm −2 for all the restrained positions of the benzocaine across the membrane. Thus, benzocaine was shifted by 0.1 nm between consecutive simulation windows. As in previous works, 21, 22, 37 two benzocaine molecules were introduced into the system to save computing time, with a separation of 3.5 nm along the z-axis to ensure the absence of relations between both molecules. Figure 2 depicts four snapshots of the system corresponding to χ = 0, where one can appreciate the benzocaine molecules placed at different positions along the z-axis.

A. Surface area and deuterium order parameters
In order to validate the simulation parameters and force field used in this work, two different properties were studied, the surface area per lipid and the deuterium order parameter for the lipids forming the bilayer in the absence of benzocaine. Figure 3 shows the average running surface area per lipid for three of the five bilayers simulated in this work, in the absence of benzocaine, as a function of their lipid compositions, χ . From Figure 3, an equilibration time of 10 ns was estimated, for all the bilayers studied in this work. Table I depicts the average surface area per lipid calculated from the last 40 ns of the trajectory length.
Unfortunately, only the surface area corresponding to pure lipid bilayers (formed exclusively by DPPC or DPPS) can be compared with experimental data. In the case of a pure DPPC bilayer, corresponding to χ = 0, a surface area per lipid of 0.683 ± 0.006nm 2 was measured, which is in good agreement with the experimental data of 0.70 nm 2 proposed by Phillips et al. 38 and Janiak et al. 39 for a DPPC bilayer at 350 K. In the case of the DPPS bilayer, corresponding to χ = 1, a surface area of 0.531 ± 0.005 nm 2 was calculated, which agrees with the experimental data of DPPS in its liquid crystalline phase, which range from 0.45 to 0.55 nm 2 . 34,40 In general, the surface area per lipid was seen to decrease as χ increased (in other words, with the increasing content of DPPS in the lipid bilayer) in our simulations.
Another property widely used to validate simulations involves the deuterium order parameter of the methylene groups along the lipid hydrocarbon tails, −S CD . Experimentally, −S CD can be calculated from H-NMR experiments of specific deuterated methylene groups along the hydrocarbon lipid tails with the axis normal to the lipid bilayer. 17,41,42 By simulation, the order parameters may be calculated from the methylene hydrogen orientation with respect the z-axis normal to where x, y, z are the local coordinates of the system, θ a is the angle made by the molecular axis with the bilayer normal, and δ ab is the Kronecker delta. From simulation, the order parameter −S CD can be determined using the relation proposed by Egberts and Berendsen, 43 where S xx and S yy are the terms of the order parameter tensor of Eq. (3).
As it was explained in previous works, 23, 25 a procedure was implemented in the S CD calculations to circumvent the fact that hydrogens are not explicitly considered in the unit atoms that form the lipid molecules in our simulations.
Thus, Fig. 4 shows the deuterium order parameters of three of the systems studied in this work. From this figure, one can appreciate the good agreement between simulation and experimental data for pure membranes of DPPC and DPPS (the only two cases for which experimental data are available), and the enhancement of the order parameters of the lipid tails with increasing χ , accompanied by shrinking of the lipid surface with increasing DPPS content.
In this regard, to explore the effect of benzocaine on the order parameter of the lipids that form the bilayer, Fig. 5 depicts the variation of the order parameters of a pure bilayer  of DPPC with 2 and 10 benzocaine molecules that can move freely inside the bilayer, corresponding to benzocaine concentrations of 4.7 and 23.3 g/cm 3 , respectively. Thus, Fig. 5 shows a slight increasing in the disorder at the end of the lipid tails with the increasing of the benzocaine (BZCA) concentration (although this effect is small to be significant considering the error bars of our simulations), in good agreement with the experimental data of Kuroda et al., 42 where they reported how the disorder increased beyond the carbon 8 in presence of benzocaine. In line with these results, almost no variation in the surface area per lipid and thickness of the lipid bilayer was measured from our simulations. Figure 6 depicts the free energy profile associated with the insertion of benzocaine into a lipid bilayer as a function of the lipid composition, χ . Note the existence of a minimum in the free energy profile located in the hydrocarbon region within the lipid bilayer. These values correspond to   Figure 6 shows how increases in the DPPS content (higher values of χ ) lead to a widening of the free energy profile along the z-axis. This phenomenon is associated with the increasing thickness of the lipid bilayer as a result of the increase in the order parameters of the lipids, as was discussed above. Figure 7 shows how the maximum in free energy associated with benzocaine insertion into a lipid bilayer fits an exponential function, which can be expressed in its linear form as follows: ln( G max )/(kJ mol −1 ) = −2.93χ + 1.48 (5) or expressed in its exponential form as follows:

B. Free energy profile, G(z)
expressed G in both equations in kJ/mol. This result confirm that an increase in the DPPS content of the lipid bilayer facilitates the insertion of the benzocaine into the lipid bilayer by weakening the barrier effect of the free energy located at the lipid/solution interface. This result may be related with the experimental evidence that the activity of certain local anesthetics is modulated by the DPPS composition of the cellular membrane (as corresponds to most of the cells of the central nervous system compared with other eukariotic cells). Thus, Baenziger et al. 44 reported how the presence of charged lipids in the membrane enhances the drug activity by increasing the partitioning of tetracaine into the lipid bilayer, which is in a good agreement with our findings.
With the aim of validating the free energy profiles obtained by the umbrella sampling method, Fig. 8 shows the free energy profile G(z) across the lipid bilayer, and the normalized distribution of probability of two free benzocaine molecules in the interior of a lipid bilayer of χ = 0.3. This figure underlines the consistency of the results obtained using two different computational techniques, where the minimum of the free energy coincides with the maximum probability of finding benzocaine in the interior of the lipid bilayer. Analogous results were obtained for the rest of systems studied in this work.

C. Enthalpy and entropy associated with the insertion process
One aspect of undoubted interest in any thermodynamic study is related with the identification of the enthalpy and entropy contribution to the free energy of the process. In this respect, from classical thermodynamics it is known that entropy S, enthalpy H, and the free energy G are related by the following expression: where H and S correspond to the variation in enthalpy (energy) and entropy (disorder) associated with the thermodynamic process. Thus, experimentally, the variation in entropy associated with a certain thermodynamic process can be obtained from the difference in free energy measured at two different temperatures, as follows: This differential equation can be solved numerically as follows: Thus, from the free energy profile of G at two different temperatures, the entropy associated with the insertion of benzocaine into the bilayer can be estimated. Once S has been calculated (from Eq. (7)), the enthalpic contribution to the free energy can be estimated as well. In our case, additional free energy profiles corresponding to χ = 0,0.3 and 1 were obtained at temperatures of 345 K, 350 K, and 355 K. To reduce the error bars, the trajectory lengths were extended to 80 ns each, which amounted to a total of 17 280 ns of extra simulation time. Figure 9 shows the enthalpy and entropy profiles associated with benzocaine insertion into the lipid bilayer, corresponding to the membranes of composition χ = 0 (a pure DPPC bilayer), χ = 0.3 (a mixture of DPPS + DPPC), and χ = 1 (a pure DPPS bilayer).
Unfortunately, even when the simulation trajectories were lengthened to 80 ns each (compared with the previous simulations of 50 ns), the error bars of S and H at the lipid/solution interfaces show long error bars, which limits the accuracy of our results. However, the error bars were smaller within the bilayer for both enthalpy and entropy. From the results depicted in Fig. 9, and based on the fact that benzocaine insertion into the membrane is a spontaneous process, it is observed how entropy is the driving force of the process, considering that benzocaine insertion into the bilayer is an endothermic process.
Thus, the shape associated with the entropy from bulk solution to the interior of the lipid bilayer may be explained as follows: 1. In bulk solution, the benzocaine dipole presents a homogeneous orientational distribution function, as evidenced in a previous work. 22 However, in the vicinity of the lipid/water interface, this distribution is strongly perturbed due to the fact that the benzocaine dipole preferably points toward the bilayer surface. This perturbation in the orientational distribution function of benzocaine is enhanced by the presence of DPPS in the bilayer. As a consequence, increasing the anisotropy of the benzocaine dipole orientation by increasing the DPPS content in the lipid bilayer should lead to a diminution in the entropy (or disorder), which is in good accordance with our simulation data. 2. Within the bilayer, an increase in entropy is expected with increasing disorder of the membrane (chaos increases, i.e., the entropy). Thus, from the results of the deuterium order parameters (see Sec. III A), an increase in the DPPS fraction in the lipid bilayer increases the order of the lipid tails (i.e., decreased disorder). Thus, an increase in order in the interior of the lipid bilayer correlates with a diminution in the entropy of the same as reflected in Fig. 9.
Unfortunately, as far as we know, no experimental data are available to make a straight comparison with our simulation results. However, Lozano, and Martínez 45 studied the transference process associated with cetoprofene from aqueous solution to the interior of a DPPC liposome. These authors proposed the following thermodynamic data for this process: G= −32.72±0.09 kJ/mol, H = 56.2 ± 1.5 kJ/mol, and S = 0.298 ± 0.008 kJ/mol.K, results that are in good agreement with our simulation data for a pure bilayer of DPPC with G = −24 ± 0.14 kJ/mol, H = 73.82 ± 0.02 kJ/mol, and S = 0.29 ± 0.01 kJ/mol.K.
Thus, considering the H and S profiles depicted above, the free energy profile G(z) can be obtained for a range of T in the absence of phase transition. Figure 10 shows the free energy profile in the range 340-360 K, for a DPPC bilayer (χ = 0). From this figure, it can be observed how an increase the temperature produces an increase in the free energy in the middle of the lipid bilayer, as might be expected from an entropic process. However, almost no variation in maximum free energy was observed at the lipid-solution interface in the range of temperatures studied. An analogous trend was shown by the other bilayers of different χ (not included in this work). This enhancement in the free energy with temperature is in a good agreement with the experimental results provided by Matsuki et al.,20 where an increase in the partition constant for benzocaine between a water/DPPC bilayer was measured.

IV. CONCLUSIONS
MD simulations of benzocaine insertion into lipid bilayers of different composition pointed to the spontaneity of this process, regardless of the lipid composition of the bilayer, with values ranging from −24 to −28 kJ/mol, for χ = 0 (pure DPPC) and χ = 1 (pure DPPS). However, even though this is a spontaneous process (from a thermodynamic point of view), as it penetrates into the bilayer, benzocaine must pass through a barrier of free energy that depends exponentially on the DPPS content of the bilayer. Thus, an increase in the DPPS fraction of the lipid bilayer facilitates the insertion of the benzocaine into the bilayer, an observation that could be related with the activity of certain drugs that depend on the lipid composition of the cell membrane.
As regards the nature and spontaneity of benzocaine insertion, it was seen how this insertion process is mainly dominated by the entropic nature of the process. Thus, as the order of the lipid bilayer associated with the presence of DPPS increases (higher χ values), so the entropic contribution to the insertion process diminishes.