INTRODUCTION

TOP ABSTRACT INTRODUCTION SIMULATION METHODOLOGY SIMULATION RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Mechanisms of the action of disinfectants have been recently reviewed (Denyer and Stewart, 1998). These authors regard disinfectants as chemical biocides with a relative lack of selectivity. They regard the vegetative bacterial cell as presenting three broad regions for biocide interaction: the cell wall, cytoplasmic membrane, and cytoplasm. A particular problem resides in the elucidation of the mechanisms of action of chemical biocides; these studies often have to be carried out under conditions that may be remote from those used to determine anti-microbial action (Bloomfield, 1991). Although biocidal compounds derive from a variety of chemical classes, the final resulting damage may show considerable similarity. The target most frequently cited in the biocide literature is the cytoplasmic membrane (Denyer and Stewart, 1998).

Agents used in household cleaning compositions are composed of structurally diverse classes of chemicals, including surface active agents, phenols, and terpenoids. The cytoplasmic membrane is a complex structure, and modes of action could involve membrane lipids or membrane proteins. Protein denaturants will disrupt transmembrane protein structure. Low molecular weight hydrophobic materials are believed to partition into the lipophilic part of the membrane and increase its fluidity, leading to disruption and cell leakage and cell death (Denyer and Stewart, 1998).

Although nonionic surfactants have traditionally been regarded by microbiologists as mild and microbiologically inactive (Russell et al., 1992), detergent alcohol ethoxylates (that is, those surfactants capable of delivering good hard surface cleaning) have now been shown to be capable of inhibiting bacterial growth (Moore, 1997). There is evidence that nonionic surfactant can interact with in vitro lipid membranes by the formation of channels through the membrane (Schlieper and de Roberts, 1977). The occurrence of "hole" formation in bilayers of long chain surfactants ("mesh" liquid crystal phase) has been demonstrated for certain nonionic surfactants by small-angle x-ray scattering studies (Burgoyne et al., 1995). Similar structures have been found in experiments on block copolymers (Hamley et al., 1993; Zhao et al., 1996) and in simulations thereof (Groot and Madden, 1998). Finally, addition of cationic surfactants to lipid membranes leads to hole formation (Gustafsson et al., 1997, 1998). It therefore seems reasonable to enquire whether the interaction of alcohol ethoxylates and phospholipids typical of bacterial membranes would lead naturally to such mesh phase that would make the bacterial cell leaky and leading to bacterial stasis and ultimately cell death. However, no evidence has been found for structured mesh phase in deuterium NMR and x-ray diffraction studies on the interaction of a phosphatidylethanolamine extract of Escherichia coli with an alcohol ethoxylate formulation (M. V. Jones, K. L. Rabone, B. A. Timimi, and G. J. T. Tiddy, manuscript in preparation).

To gain insight in this problem by direct simulation, we need a technique that enables us to simulate a patch of lipid bilayer and its uptake of surfactant, up to a time-scale at which phase transitions occur and possible holes move around. Qualitatively different phenomena occur in lipid membranes on different time scales (Tieleman et al., 1997; Pastor and Feller, 1996). On the shortest time scale of a few picoseconds, the lipids show bond and angle fluctuations of dihedral angles within the same molecule. On a time scale of a few tens of picoseconds, trans-gauche isomerizations of the dihedrals occur. Such phenomena have already been found by Heller et al. (1993), who simulated a lipid bilayer simulation over 250 ps. In 1993, this took some 20 months' run time on a Cray 2 processor. Since then, considerable progress has been made in simulating lipid bilayers and their interaction with peptides (see e.g., Tieleman et al., 1997 for a review). These authors, however, also remark that most biological questions concern much larger length and time scales than can be investigated by straightforward molecular dynamics. For instance, the simulated structure of alamethicin at the surface of a phosphatidylcholine bilayer interacts with lipid headgroups, but does not penetrate the hydrophobic core of the bilayer within a time of 2 ns (Tieleman et al., 1999), even though this molecule is a known example of a channel-forming peptide (Bechinger, 1997; Breed et al., 1997).

On a time scale of a few nanoseconds the phospholipids rotate around their axis, and on the time scale of tens of nanoseconds, two lipids switch places within one bilayer, giving rise to lateral diffusion. Within this time scale, the individual lipids orient, and lipid membranes show protrusions (Lipowsky and Grotehans, 1993). Using parallel architecture, fully atomistic molecular dynamic simulation of liquid volumes of 10 × 10 × 10 nm³ up to 20 ns (Marrink et al., 2000) are now accessible. On this time scale, the formation of micelles can be studied. Finally, on a time scale of 100 ns, peristaltic motions and undulations occur (Lindahl and Edholm, 2000). These undulations renormalize the membrane tension, which is predicted to increase with the logarithm of the lateral size (Feller and Pastor, 1996).

The key question pertinent to the problem that we want to address is: can we predict the dynamic structure of a lipid membrane and its phase changes in the presence of surfactant? More generally, we may wonder: how far we can get by pushing the hardware and by combining molecular dynamics with Monte Carlo methods (Forrest and Sansom, 2000)? The explicit simulation of the formation of a lamellar bilayer, and the analysis of the lateral correlation function to extract the bending rigidity took 18 months of CPU on an R4400 processor (Goetz et al., 1999). Extra simulation speed was gained here by using a united atom approach, where CH₂ groups are represented by a single Lennard-Jones particle (Goetz and Lipowsky, 1998). By virtue of parallelization over several processors or PC clusters, hardware developments have now pushed the limit of molecular simulation to 100 ns (Lindahl and Edholm, 2000). Nevertheless, there is a limit beyond which hardware developments cannot help us. For instance, phenomena like the cooperative motion in phase transitions, the insertion of large molecules like proteins into membranes, or membrane fusion occur on much larger time scales and are well outside the range of current simulation power. This requires simulation of the microsecond range, and a new set of phenomena could be studied if we could address the millisecond time scale.

To simulate on these time scales, there is no alternative but to simplify the model. Indeed, it is conceivable that, for large-scale collective motions, or rare events that happen on long time scales, not all atomistic details of the model are essential. For instance, biased sampling Monte Carlo simulation has been used to determine the rate-determining step in protein crystallization (TenWolde and Frenkel, 1997) on a very simplified model. The key question thus becomes: how to simplify the model such that we retain the essential physics? One possible simplification that has been applied to membrane/protein interactions is to model the protein in full atomic detail, but to represent the membrane by a mean-field potential (Biggin et al., 1997). This approach goes back to self-consistent field theory of lipid membranes (Leermakers and Scheutjens, 1988; Barneveld et al., 1992). In this technique, a flat bilayer is assumed explicitly, hence all correlated motion of molecules in the membrane and all hydrodynamic interaction are neglected at the start. For the present purpose, this method is too coarse an approximation, because modeling uncorrelated holes in a layer is impossible with this technique. For this reason, a new approach is necessary, which is more detailed than a plain self-consistent field approach, but yet is faster than simulating the united atom model based on the Lennard-Jones potential (Goetz and Lipowsky, 1998).

A promising simulation technique that might bridge the gap between atomistic and mesoscopic simulation is dissipative particle dynamics (DPD) (Groot and Warren, 1997). This method has been applied previously to single-chain surfactant membranes (Venturoli and Smit, 1999), and to the related problem of block-copolymer mesophase formation (Groot and Madden, 1998; Groot et al., 1999), where the lamellar phase is formed spontaneously in the simulation. Finally, the spontaneous formation of surfactant micelles and the formation of polymer-surfactant complexes has been simulated recently up to a time of 16 µs using the DPD technique (Groot, 2000). However, in these earlier simulations, the parameters were not related to molecules of specific chemistry. This problem is addressed here. The nature of the DPD simulation technique, and the parameterization required to simulate specific lipids and surfactants in sufficient chemical detail, is outlined in the next section. A cross-validation of the method and the results for the structural phases of lipid/surfactant mixtures are presented thereafter.

	SIMULATION METHODOLOGY

TOP ABSTRACT INTRODUCTION SIMULATION METHODOLOGY SIMULATION RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Outline of simulation method

The strategy to simulate molecular motions on length and time scales that are much larger than what can be achieved with ordinary molecular dynamics (MD) simulations is based on two main ingredients. First, atoms are lumped together into "united atoms" that describe more than one atom. In the case of hydrocarbons, we lump together three CH₂ groups into one new entity. This leads to a coarse-grained description of the molecular structure. Such a coarse graining is not a problem, because we are not really interested in where each hydrogen atom goes, but we are interested in the mobility of whole molecules, and the motion of the bilayer as a whole. For this purpose, three CH₂ groups lumped together into one bead is enough detail. The advantage of this coarse graining is that fewer particles have to be considered than there are atoms in the system, which translates into higher simulation speed.

The second ingredient used is that these new particles interact with each other by rather soft forces as the positions of the underlying atoms are smeared out. Because we want to describe the correct thermodynamics (and dynamics) on a larger length scale than an atom, we only need to reproduce the correct compressibility of the liquid and the correct solubilities of the various components into each other (Groot and Warren, 1997). To arrive at this goal, we have the freedom to choose the effective interaction as a rather soft repulsion, provided that we satisfy the above criteria. This means that we can leave out the hard-core repulsive interaction between the atoms. Because it is the hard-core interaction that forces the use of small time-steps (10¹⁵ s), the removal of this core allows a considerable increase of the time step (typically 5 × 10¹² s). The particle positions are evolved in time using the DPD method.

In the DPD method, a set of interacting particles is considered, whose time evolution is governed by Newton's equation of motion. Hence, at every time step, the set of positions and velocities, {r_i, v_i}, follows from the positions and velocities at earlier time. The force acting on a particle is given by the sum of a conservative force, a drag force and a pair-wise additive random force, i.e., f_i = _j (F + F + F), where the sum runs over all neighboring particles within a certain distance R_c. All forces depend on coordinate differences. The conservative force is given by (Hoogerbrugge and Koelman, 1992; Groot and Warren, 1997)

(1)

where a_ij is a maximum repulsion between particle i and particle j, r_ij = r_j  r_i and _ij = r_ij/|r_ij|. Between neighboring particles on a chain, an extra spring force is defined that binds the particles together, given by

(2)

The drag force F and the random force F act as heat sink and source, respectively, so their combined effect is a thermostat. They are given by the random force F = (r_ij)r_ij/ and the drag force F = 1/2²(r_ij)²/kT(v_ij · r_ij)r_ij, where  is a random variable with zero mean and variance 1, and (r) = (1  r) for r < 1 and  = 0 for r > 1. Following Groot and Warren (1997), we use the fixed noise amplitude  = 3. This particular thermostat is special in that is conserves (angular) momentum, which leads to a correct description of hydrodynamics (Español, 1995). We choose the particle mass and the temperature and interaction range as units of mass, energy, and length, hence m = kT = R_c = 1, and the simulated time is expressed in the natural unit of time

(3)

The DPD method, in general, has been shown to produce a correct (N, V, T) ensemble if the fluctuation-dissipation relation is satisfied (Español and Warren, 1995; Groot and Warren, 1997). At every time step, the set of positions and velocities, {r_i, v_i} is updated from the positions and velocities at earlier time using a modified version of the velocity-Verlet algorithm:

(4)

The masses of the particles are put at 1, so that the force acting on a particle equals its acceleration. The force is updated once per iteration. Because the force depends on the velocities, the velocity in the next time step has to be estimated by a predictor method. This is done in the second step of our algorithm. The velocity is corrected in the last step. If the parameter  is put at  = 0.5, this scheme equals the velocity-Verlet algorithm (Allen and Tildesley, 1987). However, here we use  = 0.65, where we find a very accurate temperature control, even at the time step t = 0.06 that is used. A more systematic study into the influence of parameter  was presented by Den Otter and Clarke (2000).

In the simulations, a bilayer of phospholipids is simulated in a periodic cell, where the bilayer is oriented perpendicular to the x-axis. Therefore, the local density of each component is measured in thin slabs perpendicular to the x-axis, and the stress tensor is averaged locally and over the whole system. The stress tensor leads to the surface tension via

(5)

where A is the area of the yz-plane, and F_ij is the total conservative force between particles i and j. To simulate a bilayer at a predefined surface tension (usually zero tension), sub-averages of the surface tension are taken during the run. To control the surface tension, frequent transformations are made where all y and z-coordinates are stretched while the x-coordinates are shrunk, or vice versa. Hence,

(6)

where the amount of stretch  is determined from the deviation of the actual stress from the externally applied stress. Note that this transformation conserves the total volume of the sample. Generally, the procedure must be repeated a number of times until the desired surface tension is reached. From then on, the area should be fixed, because a simulation with repeated transformations does not correspond to a correct Hamiltonian system. To have a Hamiltonian system, one would have to introduce an external force that couples to the instantaneous surface tension, but this has not been implemented here. An alternative approach is to apply the transformation in Eq. 6 at random, but to accept or reject a transformation according to the Monte Carlo algorithm (Venturoli and Smit, 1999).

Simulated system and physical length and time scales

The system that is simulated is composed of water, phosphatidylethanolamine (PE), and a nonionic surfactant, either C₁₂E₆ or C₉E₈. The structure of PE contains two C₁₅ chains connected to a glycerol group. The molecule is shown in Fig. 1, together with its mapping on the coarse-grained DPD model. To construct a mesoscopic model, we first need to determine the volume of the simulation beads, and hence determine the length scale. We choose a coarse-graining where three carbon atoms are taken together and grouped into one bead.

View larger version (20K): [in this window] [in a new window]   FIGURE 1   The simulated phospholipid and its mapping on the DPD model. Each bead represents roughly the same liquid volume.

Within the same model the surfactant C₁₂EO₆ is represented by the DPD bead structure c₄e₄. Hence, whereas, a c-bead represents three carbon atoms, an e-bead represents 1.5 EO group, or the glycerol-linking unit in the phospholipid. This mapping is justified by studying the partial volumes of (CH₂)₃ and EO groups. From Lu et al., (1993), we find that the volume of C₁₂H₂₅ is 350 Å³, the volume of (OC₂H₄)₆OH is 380 Å³, and the volume of a water molecule is 30 Å³. If we subtract the H₂O volume from the E₆OH volume, we find exactly the same volume for 6 EO groups as for C₁₂H₂₅. Hence each (e- or c-) bead represents a liquid volume of 90 Å³. Because a water molecule has a volume of 30 Å³, and the water beads are to represent the same volume of 90 Å³, the water beads (w) must represent three water molecules, whose volume also adds up to 90 Å³. Because the simulated bead density is R = 3, a cube of R contains three beads and therefore corresponds to a volume of 270 Å³. Thus, we find the physical size of the interaction radius,

(7)

Because noise and friction are included in the simulation method, the hydrodynamic regime is simulated already with few particles and time steps (Español, 1995). The consequence of this strategy, however, is that we have lost track of our physical unit of time. To gauge our time unit, we match the long-time diffusion constant of water.

Some care must be taken here. The self-diffusion constant of a water-bead is not the same as the self-diffusion constant of water, because the bead represents three water molecules. When these move over the vectors R₁, R₂, and R₃, their center of mass moves over the vector R_w = (R₁ + R₂ + R₃)/3. Hence, the ensemble average of the mean square displacement of the water beads is

(8)

where R² is the mean square displacement of a water molecule. Because the mean square displacement of the water beads is one-third of that of the water molecules, the diffusion constant of the beads is one-third of that of water.

The diffusion constant of the water beads was obtained by averaging the mean square displacement over three runs of 50,000 time steps each, and determining the slope of R(t)/6 against time. At the noise and repulsion parameters used here, this led to the bead-diffusion constant,

(9)

Equating this to the experimental diffusion constant of water (Partington et al., 1952), D_water = (2.43 ± 0.01) × 10⁵ cm²/s, leads, together with Eq. 8, to the time scale

(10)

It is particularly instructive to see how length and time scales change if we would incorporate more atoms in one bead. Let, in general, a bead correspond to N_m water molecules. The number N_m can be viewed as a real-space renormalization factor. The case discussed above thus corresponds to the choice N_m = 3. Hence, a cube of volume R represents N_m water molecules, where  is the number of DPD beads per cubic R_c. Because the physical volume of this cube equals 30 N_m ų, the length scale R_c follows as

(11)

Following the reasoning above to match the diffusion constant of pure water, we find that Eq. 10 in the general case becomes

(12)

In this equation, it is implicitly assumed that the repulsion parameter between equal beads is fixed at the value of a = 78, and that the bead density is fixed at  = 3. Following Groot and Warren (1997), a modified velocity-Verlet algorithm is used that allows for time steps of t = 0.06 at this repulsion parameter, hence, time-steps of 5.3 ps are taken.

At this point, we can understand why the DPD method is so much faster than straight-forward molecular dynamics. There are two combined effects that lead to speed-up. The first contribution comes from the low Schmidt number in the simulation (Groot and Warren, 1997). The Schmidt number is the ratio between viscosity and the self-diffusion constant, Sc = /D. In an ordinary liquid like water, this ratio is roughly Sc  1000, whereas, in the DPD method, we have Sc  1. The origin of this difference can be traced back to the removal of the hard core from the interaction potential. This hard core leads to a caging effect, so that an atom undergoes many collisions before it is actually transported. The soft potential used here removes this caging affect, so that the diffusion of particles is increased by a factor of 1000. The second factor leading to fast simulation is the scaling of the physical time with the renormalization factor N_m as in Eq. 12. On top of the power , by which the physical time scale increases, the amount of CPU time will decrease inversely proportional to N_m if we want to simulate a given volume simply because we have to update the position of fewer objects. Thus, for a given system volume, DPD can be expected to be faster than MD by a factor of roughly 1000 N  2 × 10⁴ for N_m = 3 and about 10⁵ for N_m = 6, independent of hardware and disregarding the CPU time spent on evaluating the (relatively long-ranged) Lennard-Jones potential. A direct comparison with the early work by Heller et al. (1993) showed that the present simulation is 10⁷ times faster. Since 1993, software and hardware improvements have accelerated computational speed by roughly a factor of 3000 (see e.g., Marrink et al., 2000), so that the actual method is indeed some 10⁴ times faster than MD when N_m = 3 is used.

Interaction parameters

To find in practice the interaction parameters for this model, we need to match the liquid structure function in the limit k  0, because this determines the free energy change associated to density fluctuations. This, in turn, is related to the compressibility and solubilities. Note that the pressure itself drops out in an (N, V, T) ensemble, because this is a linear variation of the free energy. It was previously proposed that the following relation should hold (Groot and Warren, 1997):

(13)

where  is the bead density in the simulation, and n is the density of, e.g., water molecules in liquid water. However, this relation only holds if one DPD bead corresponds to one water molecule. In general, the system should satisfy

(14)

where N_m is the number of water molecules per DPD bead. In the previous section, N_m has been chosen at N_m = 3. For this value, the compressibility of water at room temperature is matched if the repulsion parameter in Eq. 1 is determined at (Groot and Warren, 1997)

(15)

where a_ii is the repulsion parameter between particles of the same type. Note that it is taken the same for all liquid components, because we actually simulate equal liquid volumes for all components.

The next observables to match are the mutual solubilities. In polymer chemistry, this is usually expressed by specifying the Flory-Huggins -parameters. This parameter represents the excess free energy of mixing in the Flory-Huggins model. This is a cell model, where every cell is filled by a fraction  of A molecules and by a fraction 1   of B molecules. Hence, the lattice is completely filled. If A is a polymer that occupies N_A cells, and B is a solvent that occupies N_B cells, then the free energy per cell (disregarding constants and terms linear in ) can be written as

(16)

Different polymers usually tend to segregate, and, to represent this behavior, we impose a larger repulsion between unlike beads than between beads of the same type. It has been established that the -parameter is linearly related with the excess of the AB repulsion over the AA repulsion (Groot and Warren, 1997). Following their procedure for the present repulsion between equal beads, simulations have been run for mixtures of A₃ and B₃ chains (6000 beads in total) in a box of size 10 × 10 × 20R. The volume fraction of A in the majority B phase has been measured. This is inserted in the mean-field expression for the binodal for N_A = N_B = N = 3:

(17)

which should be valid far away from the critical point. Thus the Flory-Huggins -parameters have been obtained for AB repulsion a_AB = 82, 83, 84, and 85, which led to the correspondence

(18)

where a = a_AB  a_AA is the excess repulsion.

The pertinent -parameters are determined by matching relevant thermodynamic data to the same Flory-Huggins model. To obtain the parameter between hydrocarbon and water, the binodal is determined numerically from the set of equations for molecules of unequal volume:

(19)

where µ^I, µ^II, p^I, and p^II are the chemical potentials and osmotic pressures in the dilute and in the dense phase, respectively. This binodal is subsequently matched to the solubility data of hexane, heptane, and octane in water. The volume fraction of oil in water for these hydrocarbons is 1.8 × 10⁵, 3.5 × 10⁶, and 9 × 10⁷, respectively, whereas the volume fraction of water in these oils is 7.3 × 10⁵, 6.2 × 10⁵, and 5.6 × 10⁵, respectively (Shaw, 1989). For the case where three water molecules are represented by one DPD bead, the interaction parameter is found as _{hydrocarbon-water}  6.0, and appears to be relatively independent of temperature. Because this parameter scales linearly with the bead volume, the value 6.0/N_m = 2.0 should be compared to values cited in the literature for the -parameter per carbon atom: Barneveld et al. (1992) also derived  = 2.0 by matching the critical micelle concentration of surfactant solutions, whereas Leermaker and Scheutjens (1988) used  = 1.6. It should be noted that there is an inconsistency in the Flory-Huggins model: when the -parameter between hydrocarbon and water is determined by matching the solubility of water in oil, a higher value _{hydrocarbon-water} = 9.3 is obtained, which does depend on temperature. A more rigorous parameterization is possible using a closed expression for the binodal that follows from the DPD simulations directly. This avoids using mean-field theory (Wijmans et al., 2001), but the present study is based on _cw  6.0.

The second -parameter to match is the one between polyethyleneoxide (PEO) and water. The problem here is that PEO and water at room temperature mix in all ratios, hence the solubility does not lead to a parameter value. An alternative way is to fit experimental adsorption data. To model PEO adsorption on polystyrene latex, Cohen Stuart et al. (1984) used _{ethyleneoxide-water} = 0.45, but they fail to give either a source for this number, or the precise data that they fitted to obtain it. At elevated temperatures (T > 100°C), water and PEO no longer mix ideally. Hence, an alternative route is to describe this demixing at higher temperatures, and to extrapolate the temperature dependence of the -parameter back to room temperature. This way, Barneveld (1991) estimated this -parameter as a function of temperature and found the value _ew  0.30-0.38 at room temperature. However, he only took the cloud point into account where mean-field theory is least reliable. Taking the whole shape of the binodal into account and extrapolating back to room temperature, we find, from the experimental data by Seaki et al. (1976), _ew  0.30 ± 0.04, which is close to the value obtained by Barneveld.

A third important -parameter is the interaction between PEO and hydrocarbons. To estimate this -parameter, only sparse experimental data is available. Ideally, we would need the solubility of EO₆ in C₁₂H₂₅, and vice versa. What is available is neutron-scattering data of C₁₂E₆ at the air-water interface (Lu et al., 1993). Assuming that this is not too dissimilar from the oil-water interface, these data can be compared to DPD simulation data of an oil-surfactant-water system. The experiment shows a significant overlap between the surfactant head group and the surfactant tails, and both distributions can be described by Gaussian peaks with full-width-at-half-height  = 13.5 ± 1.0 Å. Furthermore, the sum of head and tail peaks was found to have a full-width-at-half-height of 18.3 ± 1.0 Å (Lu et al., 1993).

To reproduce these data, a DPD simulation was performed containing 440 C₁₂H₂₅ molecules, 140 C₁₂E₆ molecules, and 1620 water beads (=4860 molecules) in a box of size 18.72 × 8.95 × 8.95R (=121.0 × 57.9 × 57.9 Å³), where _ew = 0.3 and _cw = 6 were used. To arrive at the same amount of overlap in the simulation as in the experiment, a _ce parameter much smaller than the hydrocarbon-water parameter needs to be used. Because an EO group consists of of hydrocarbon and of oxygen, it is not unreasonable to assume that the -parameter between C and EO is of that between C and water. Therefore, simulations were performed at _ce = 2. The box size was selected such that the surface tension vanishes within the simulation error. The area per molecule thus found in the simulation (48 Å²) is slightly smaller than the experimental value (55 Å²). In this simulation, we find two Gaussian peaks of the same half-height width and the same overlap as seen in the experiment (see Fig. 2). In the simulation, the tail peak is 13.3 ± 0.7 Å wide and the head-peak width is 13.3 ± 0.8 Å, which compare very well with the experimental 13.5 ± 1.0 Å. The width of the sum of head and tail peaks is 18.1 ± 0.7 Å in the simulation, and 18.3 ± 1.0 Å in the experiment. Hence, the assumed value for -parameter between hydrocarbons and EO beads, _ce = 2.0, is consistent with the data.

View larger version (21K): [in this window] [in a new window]   FIGURE 2   DPD simulation of a C12H25/C12E6/ water interface. The density profiles of the surfactant head and tail groups are indicated by the dashed curves.

Finally, the -parameters describing the head group of the lipid molecule have to be defined. Because these groups contain more oxygen than EO does, and also have partial charges, they have been treated as if they were water, with respect to C and EO. To represent the partial charges, one might increase the repulsion between the head groups mutually, and reduce the repulsion of the head groups with water. To study the influence of these variations, two sets of parameters were studied, listed below both in terms of -parameters and as the actual repulsion parameters.

(20)

The resulting bilayers have been compared with molecular dynamic simulations on 1-palmitoyl-2-oleoyl-sn-glycero-3-phophatidylcholine (POPC) (Heller et al., 1993). This comparison is described in more detail in the next section. It is found that the two -parameter sets lead to very similar density profiles and areas per lipid head group. Hence, although the head-group parameters that have been varied here are not very well known, this is, in practice, not a problem, because these parameters are not critical for the result.

	SIMULATION RESULTS

TOP ABSTRACT INTRODUCTION SIMULATION METHODOLOGY SIMULATION RESULTS DISCUSSION AND CONCLUSIONS REFERENCES

Validation of the DPD simulation

The density profiles obtained from DPD simulations of lipid bilayers containing 200 lipid molecules are compared to previous MD simulation results (Heller et al., 1993) in Fig. 3. Each DPD run of 50,000 time steps took 1 h of CPU on a single processor R8000 workstation. The difference between the results of the two parameter sets is so small that only the second set is shown. The two peaks in the CH₂ density profile arise because the CH₃ chain endpoint is localized in the center of the bilayer, and is not included in the average. In the DPD simulation, the average is over all but the last c-bead, i.e., the last three carbon atoms are excluded from the average. This explains the difference between the DPD and MD gap in the middle of the CH₂ density profile. A further difference arises because the e-beads in the DPD simulation contain ester bonds that are not included in the glycerol backbone. The difference in parameter sets 1 and 2 lies in the (not very well known) parameters for the phosphatidyl groups. The fact that this variation in the parameters leads to only slight differences in the density profiles demonstrates that the system is not very sensitive to these parameters. The good correspondence with the MD simulation for either set thus gives confidence in the reliability of the model and parameters.

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Molecular dynamics result for CH2 density profile () and glycerol backbone () taken from Fig. 18 of Heller et al. (1993) compared to the DPD result for parameter set 2.

The difference in the two parameter sets is reflected in the area per lipid molecule in the bilayer. Experimental values for the area per head group of DPPC in the L phase are reported to be 70 Å² (Rand and Parsegian, 1989), or to vary as 57.6-70.9 Ų (Nagle and Wiener, 1988). For DOPC, 60 Å² is reported (Wiener and White, 1992a,b), and finally 65.5 Å² is reported for simulations of POPC (Heller et al., 1993). See Nagle and Tristram-Nagle (2000) for a recent overview of the experimental data. For comparison, our parameter set 1 leads to an area of 62.1 ± 0.1 Å² per lipid, and parameter set 2 leads to 66.8 ± 0.1 Å². These values cover the range of experimental data. For the purpose of the present simulations, we did not take into account the renormalization of the surface tension with the size of the membrane (Feller and Pastor, 1996). Recent simulations show that the area per lipid increases from 63 Å² at 256 lipids to 63.5 Å² at an infinite system (Lindahl and Edholm, 2000). Because this difference falls within the uncertainty of the present parameterization, it is not relevant for our simulations. Related to the area per lipid is the stress profile through the membrane. An example of the stress profile is provided in Fig. 4. It compares well with the stress profile simulated by Venturoli and Smit (1999), but it is at variance with the stress profile predicted by Goetz and Lipowsky (1998), which shows several oscillations in the stress profile. These may well be caused by the use of short hydrophobic tails with explicit hard (Lennard-Jones) interaction. Figure 4 also shows the volume fraction profiles for the lipid-chain end points, water, and head groups. It should be noted that, although the chain ends tend to peak around the center of the bilayer, the distribution is not as narrow as that of the CH₃ groups seen in MD. A straight comparison is not possible, however, because the last bead not only models the CH₃ group, but also the last two CH₂ groups. Note further that the main positive contribution to the stress comes from the water-hydrocarbon interface. These peaks (the shaded areas in Fig. 4) integrate to surface tensions of 20 and 21 mN/m, respectively, and are compensated to zero by the negative tension in the head groups and the tails. These simulated values compare well to the surface tension between medium chain triglyceride oil and water, which is 23.5 mN/m (Hugelshofer et al., 2000).

View larger version (26K): [in this window] [in a new window]   FIGURE 4   Volume fraction profiles of DPD simulated lipid bilayer, including the profiles for the head groups and for the end points of the hydrocarbon tails and water (top). The lower panel shows the stress profile across the interface.

We now focus attention to the validation of the time-scale of the simulation. To this end, the bilayer simulation is used to obtain the lateral diffusion constant of the lipid molecules, which can be checked against experimental and MD simulated values. For POPC, the simulation result is D = 0.073 × 10⁵ cm²/s (Heller et al., 1993). Some care must be taken here, because this simulation was too short for two molecules to swap places in the layer. Experiments they cite (Cohen and Turnbull, 1959; Pfeiffer et al., 1989) indicate lateral diffusion constants of DOPC in the range D = 0.036 × 10⁵ cm²/s and D = 0.02 × 10⁵ cm²/s.

In the DPD simulation, the mean square displacements parallel (and also perpendicular) to the bilayer were averaged. These mean square displacements, divided by 2 (perpendicular) or by 4 (parallel motion) and scaled according to the proper length and time units, are shown in Fig. 5 for parameter set 2. The slopes of the lines give the respective diffusion constants. The mean displacement of the lipid molecules during the run is 8.0 nm. Because the area per molecule is 0.668 nm², a typical distance between molecular centers is 0.82 nm; hence, during this run, each lipid molecule traveled 9.8 times that distance on average. The diffusion constants derived from the motion of the CH₂ groups of the lipid and of the head-groups in the DPD simulation are both D = 0.06 × 10⁵ cm²/s, and the perpendicular diffusion vanishes within the accuracy of the simulation because the total simulation time is way beyond the time scale on which lipid protrusions from the bilayer occur. The lateral diffusion constant found here (0.06 × 10⁵ cm²/s) compares well with the experimental (0.02-0.04 × 10⁵ cm²/s) and MD simulation (0.07 × 10⁵ cm²/s) values in the literature, even though no attempt has been made to match the relative viscosity of the lipid phase. The diffusion results given here are obtained for parameter set 2; for parameter set 1, the results are the same within the simulation error. Notably, when the glycerol solubility is reduced from _ew = 0.3 to _ew = 0.4, the lateral diffusion of the lipids drops to 0.05 × 10⁵ cm²/s. Hence, lateral diffusion is sensitive to the solubility of the glycerol group.

View larger version (20K): [in this window] [in a new window]   FIGURE 5   Mean square displacement of water and lipids in the DPD bilayer simulation with set 2 parameters, divided by 2 for diffusion perpendicular to the membrane, and by 4 for diffusion parallel to it. The slope of this line gives the diffusion constant in 105 cm2/s.

The diffusion constants found for water are D = 2.21 × 10⁵ cm²/s and D = 0.004 × 10⁵ cm²/s, where the factor of 3 has already been taken into account to convert from DPD beads to water molecules. The diffusion constant of water parallel to the bilayer is reduced by 9% relative to the value in bulk water, which is not unreasonable for diffusion in a narrow gap. The error of the lateral water diffusion constant is estimated at about 1%. The water diffusion perpendicular to the layers is determined by the penetration of the bilayer by water. Because the solubility of water in oil at the present -parameter is overestimated, the real transverse water-diffusion constant is estimated at D  2 × 10⁹ cm²/s.

Diffusion of water through lipid/surfactant bilayers

To investigate the structure of mixtures of lipid/surfactant bilayers, a series of simulations was done with 540 surfactant or lipid molecules at a mole fraction of surfactant varying from 10 to 100% in steps of 10%. Simulations have been done both with C₁₂E₆ (bead structure c₄e₄) and with C₉E₈ (bead structure c₃e₅). First, we concentrate on the results regarding C₁₂E₆. The bilayer with 100% C₁₂E₆ should fall apart because the lamellar phase is unstable at the given concentration (32% volume fraction). In the course of the first 100 ns, however, the system forms a perforated lamellar sheet with holes that move around and remain meta-stable for at least 700 ns. This could be a finite size effect, but also a finite time effect because the cession of liquid bridges is a slow process (Groot et al., 1999).

The system with 90% mole-fraction C₁₂E₆ displays holes that move around. They are relatively stable; the typical lifetime is 20-40 ns, and, over the course of the run, there is hardly ever a conformation without a hole in the patch of 124 nm². Systems with 80% mole-fraction surfactant or less do not show stable perforated conformations. Occasionally, small holes do appear, but these disappear very quickly. At 80% mole-fraction surfactant the typical life time of a hole is some 0.4-1 ns (see Fig. 6).

View larger version (75K): [in this window] [in a new window]   FIGURE 6   Hydrophobic part of mixed bilayers containing 80% mole-fraction C12E6 (left) and 90% (right). The surfactant C12 chains are represented in gray, the lipid C15 chains are black. The hole in the left conformation is transient, at the right they are stable.

The diffusion of water perpendicular to the layers provides a useful way of testing the porosity of the bilayers. To interpret these result, we first solve the diffusion equation for the mean square displacement of water in a narrow slit of impenetrable walls at distance L. For diffusion constant D, the mean square displacement perpendicular to the walls follows analytically as

(21)

The first expression is exact, the second is a good approximation that can be used for curve fitting. Note that, in discarding all terms higher than n = 3, for the last term, we have replaced 96/81⁴  0.0122 by (1  96/⁴)  0.0145, otherwise, the right-hand side of Eq. 21 would not vanish in t = 0. To find a reliable value for the transverse diffusion coefficient, the mean square displacement of water normal to the bilayers is recorded, and the transverse diffusion is obtained by fitting the data to

(22)

where R is the diffusion in a slit given by Eq. 21, and a and  are free-fit parameters. The fits are shown in Fig. 5, together with the raw data for the surfactant-free bilayer. The resulting transverse diffusion of water is shown in Fig. 7 as a function of the surfactant mole fraction in the layer.

View larger version (20K): [in this window] [in a new window]   FIGURE 7   Transverse diffusion of water through the bilayer as a function of surfactant volume fraction for C9E8 () and C12E6 (). The bars give the regions where no holes, fluctuating holes, and stable holes occur, respectively.

It is found that up to and including a mole fraction of 50% C₁₂E₆, the diffusion of water through the bilayer is independent of the amount of added surfactant. From that point onward, the water diffusion through the layers steadily increases. This increase is attributed to the formation of small holes that open and close on a time scale of 0.5 ns or less, depending on the surfactant content. Permanent holes occur from 90% surfactant. Consequently, the transverse diffusion constant increases sharply above 80% surfactant.

The transverse diffusion of water through the PE/C₉E₈ bilayer is also shown in Fig. 7, for comparison with the C₁₂E₆ results. Again, here, we see the same qualitative picture involving the creation and annihilation of small transient holes, followed by a phase with stable holes at higher surfactant content. However, the range of surfactant concentrations over which transient holes occur is much shorter than for C₁₂E₆. Also, the point where stable holes are formed is reached already at 70% mole-fraction surfactant, rather than 90%. The real transition points must be between 60 and 70% for C₉E₈ and between 80 and 90% for C₁₂E₆, i.e., some 56% and 78% by weight, respectively. The approximate ranges where no holes, fluctuating holes, and stable holes are found are indicated in Fig. 7 by the black and grey bars.

Rupture of bilayers under strain

So far, we have been concerned only with membranes at vanishing surface tension. However, in many cases, it was found that the dividing cells are especially vulnerable. Dividing cells are not necessarily in a state of vanishing membrane tension. For instance, when yeast cells divide, their cell membrane buds out of the cell wall and is no longer protected by it. Instead, the membrane is exposed to the solution. The osmotic pressure difference  between the inside and outside of the cell then leads to a finite surface tension on the membrane, which is given by

(23)

where R is the mean radius of curvature of the bud. The membrane will react to this osmotic pressure by expanding, which is an obvious prerequisite for cell division when the budding mechanism is pertinent. If the membrane cannot withstand this expansion, the cell will die. For these reasons, it is prudent to simulate cell membranes under strain, rather than to study them at zero surface tension, as far as the mechanism for cell death is concerned. Simulations of mixed membranes of lipid and C₁₂E₆ were undertaken in which the membrane is stretched over time, leading to increasing tension and, ultimately, to rupture. An example of this process is shown in Fig. 8, where the actual creation and expansion of holes is monitored. The successive frames are taken at time intervals of 1.2 ns, and the patches are 17 × 17 nm² across. This membrane consists of 70% PE and 30% C₁₂E₆. It ruptures when its area is increased by 74%.

View larger version (53K): [in this window] [in a new window]   FIGURE 8   Sequence of frames showing the rupture of a 30% C12E6/lipid bilayer. The time between two frames is roughly 1.2 ns. The size of the slab is 17 × 17 nm2.

A full stress history is shown in Fig. 9 for the membrane containing a mole fraction of 50% surfactant. The system is left to equilibrate over 1000 time steps (5.3 ns), after which time the y- and z-coordinates are expanded by a factor 1.03, and the x-coordinates are contracted by a factor 0.94. This cycle is repeated 12 times. The curve in Fig. 9 gives the running average of the surface tension, averaged from each expansion point. Hence, the last point in each block is the most accurate estimate of the surface tension at that particular area. These data points are collected in Fig. 10, together with the data for 10 and 80% surfactant. The stress is plotted against the area of the bilayer, divided by the area at vanishing surface tension. Each simulation shows a clear rise in surface tension, up to a critical point where<