The membrane models used for this study are equilibrated fully hydrated dimyristoyl-phosphatidylcholine (DMPC) bilayers. At the temperatures set for the study, i.e., 300 K, the bilayer is in the biologically relevant liquid crystal L phase. We considered several systems. First, for the bare bilayer, we studied a small system consisting of 64 DMPC units and 1,645 water molecules. To investigate size effect, a larger system was constructed by replicating the system in the x and y directions of the water-bilayer interface, resulting in an assembly of 256 DMPC units.
The system containing the peptide channel was taken from Tarek et al. (2003). It consists of a nanotube formed of eight cyclo[(L-Trp-D-Leu)3-L-Gln-D-Leu] subunits, organized in an antiparallel, ß-sheetlike channel embedded in a hydrated DMPC bilayer (218 lipid molecules). Such cyclic peptides self-assemble into hollow tubular structures by means of a network of hydrogen bonds that connect the stacked subunits (Bong et al., 2001; Ghadiri et al., 1993). The complete system (synthetic channel, 215 lipid units, and 6,469 water molecules) represented a total of 46,097 atoms.
To investigate the effect of the external field on a short DNA strand located at the lipid-water interface, we considered a well-equilibrated palitoyl-oleyl-phosphatidylcholine (POPC) lipid bilayer (288 lipids), with excess water in which a 12 basepair 5'-cgcgaattcgcg-3' ecor1 DNA duplex was placed in the aqueous phase near the lipid headgroups. The complete systems comprised the DNA duplex, 288 POPC lipids in united atoms representation, 14,500 water molecules, and 22 counterions to balance the DNA charges (total of 65,609 atoms). The number of lipid units and of water molecules considered is chosen such that the size of the system precludes interactions between the DNA strand and its images due to the use of periodic boundary conditions in the simulation. POPC was chosen for convenience as preequilibrated configurations of a united-atom model of a phosphatidylcholine lipid were at hand.
The MD simulations presented here were carried out using the program NAMD targeted for massively parallel architectures (Bhandarkar et al., 2002; Kale et al., 1999). All systems were examined in the (NPT) ensemble using three-dimensional periodic boundary conditions. The equations of motion were integrated using the multiple time-step algorithm (Izaguirre et al., 1999, 2001; Martyna et al., 1996). A Langevin piston was used to maintain the pressure at 1 atm. The temperature of the system was fixed at 300 K. Long-range, electrostatics forces were taken into account using the particle mesh Ewald approach (Darden et al., 1993; Essmann et al., 1995). The water molecules were described using the TIP3P model. Bond stretching, valence angle deformation, and torsional and nonbonded parameters of the cyclic peptides forming the nanotube, of the DNA nucleic acids, and of the DMPC and POPC lipid units were extracted from the all-atom CHARMM force field (MacKerell et al., 1998). The united atoms approximation was used for the POPC lipid acyl chains. The similarities of results presented in the following for simple lipid bilayers and those obtained by Tieleman (2004) indicate that the results are not highly dependant on the details of the simulation algorithms or of the force field used.
In experiments, a potential difference V may be imposed by means of a voltage clamp or a pulse. In simulations, due to the small sizes of the systems and the use of periodic boundary conditions, it is impossible to impose a TM voltage by simple addition of explicit ions to the solution at both sides of the bilayer (Roux, 1997; Tieleman et al., 2001) Here, we adopted the same strategy as in Tieleman et al. (2001), i.e., applied an external electric field perpendicular to the membrane plane to maintain a fixed voltage difference across the bilayer (Tieleman et al., 2003, 2004). In practice, this is done by adding a force on all the atoms bearing a charge where E is the constant external electric field.
For the analysis of the electrostatic properties of membrane, we recorded the electrostatic potential across the bilayer as the electroporation of the bilayer is taking place. This potential may be estimated using Poisson's equation and derived from the MD simulation as a double integral of the molecular charge density distributions at Z, following
neglecting therefore the explicit electronic polarization (Tieleman et al., 1997
; Tobias et al., 1997
). Phosphadilylcholine headgroups have a large dipole moment. In hydrated bilayers, the orientation of the lipid headgroups with respect to the membrane surface is often found to average around 30°. This causes a nonnegligible local electric field that is partially compensated by a specific orientation of interfacial water molecules and results in a net dipole potential across each interface (Cheng et al., 2003
; Gawrisch et al., 1992
; Shinoda et al., 1998
), i.e., between the interior of the hydrocarbon layer and the aqueous phase. For this study, and in contrast to previous simulations (Berger et al., 1997
; Tieleman et al., 1997
; Tobias, 2001
), we consider the dipole across the whole membrane. Owing to the symmetry of the bilayer and in the absence of salt, the total dipole across the bilayer is null. When an external electric field is applied, one expects that water molecules and lipid headgroups reorient, changing therefore the electrostatic properties of the membrane and hence the measured total dipole potential.
In agreement with Tieleman (2004), we find that the applied electric field E induces a voltage difference over the whole system where Lz is the size of the simulation box in the direction perpendicular to the applied field. For instance, as we will see later, in the case of the bare lipid bilayers ( Å at rest) the total potential drop across the systems is 3 and 6 V for the applied fields of intensity E = 0.5 V.nm–1 and 1.0 V.nm–1, respectively.
Except for the system containing the peptide nanotube where the initial configuration was taken from our previous work, the two other systems were first equilibrated without application of the transverse electric field to afford initial configurations. The lengths of the various simulations ranged from 5 to 10 ns, depending on the system and the trajectories as will be indicated below for each system. As we will see in the following, these timescales are long enough for the electroporation to occur.