J. Phys. Chem. B 2010, 114, 7261–7264
7261
Self-Assembly of Lamellar Lipid-DNA Complexes Simulated by Explicit Solvent Counterion Model Lianghui Gao,* Jun Cao, and Weihai Fang College of Chemistry, Beijing Normal UniVersity, Beijing 100875, China ReceiVed: March 9, 2010; ReVised Manuscript ReceiVed: April 22, 2010
The dissipative particle dynamics simulations with explicit solvent and counterions are used to mimic the self-assembly of lamellar cationic lipid-DNA (CL-DNA)complexes. We found that the formation of the complexes is associated with the releasing of 70% DNA counterions and 90% lipid counterions. The trapped DNA and CL charges together with their counterions inside the complex still keep the interior neutral, which stabilized the structure. Simulations in constant pressure ensemble following the self-assembly show that the DNA interaxial spacing as a function of the inversed CL concentrations 1/φc is linear at low φc and nonlinear at high φc. The attraction between the DNA and the CLs as well as the repulsion between the DNA strands impose stretching stress on the membrane so that the averaged area per lipid is dependent on the CL concentration, which in turn determines the behavior of the DNA spacing. Introduction Coarse-grained (CG) dissipative particle dynamics (DPD) simulations1,2 have been used to study the dynamic behaviors of soft matter materials in microsize, such as the self-assembly of phospholipid membranes,3-5 artificial polymeric vesicles,6 fusion of vesicles,7-9 and interactions between peptides (or protein) and lipid membranes.10,11 It provides feasible windows on the length and time scales relevant to important biophysical processes that are not reachable by classical and all-atom molecular dynamics and retains the hydrodynamic modes that are missing in other coarse-grained techniques, such as Monte Carlo and Brownian dynamics. Groot12 also includes electrostatic interactions in DPD simulations by solving the electrostatic field locally on a grid and applying it to study the interactions between a cationic polyelectrolyte and anionic surfactant. We then successfully applied it to investigate the mechanism of antimicrobial peptide translocation across a lipid bilayer.11 We further develop this method into a constant pressure ensemble by subtracting the contributions of the self-energy of electrostatic interaction from the total force and the virial.13 Thus, this method enables us to study many interesting biological systems carrying net charges, such as the cationic lipid-DNA (CL-DNA) complexes for gene therapy. CL-DNA complexes can deliver DNA to the nucleus of a cell to replace defective genes or add missing genes. X-ray diffraction has revealed a novel multilayered structure of the CL-DNA complex with alternating lipid bilayer and DNA in which the DNA strands form a two-dimensional smectic phase intercalated between a lipid bilayer14 and an inverted hexagonal phase.16 Theoretical calculations based on the elastic theory of semiflexible multifold17,18 and Poisson-Boltzmann equations19-21 have been carried out to explain the observed structure of such complexes, but the continuum theory lacked detailed molecular interactions. All-atom molecular dynamics (MD) simulations were then performed to investigate the stability of the preformed sandwiched structure of the CL-DNA complex in which there * To whom correspondence should be addressed. E-mail: lhgao@ bnu.edu.cn.
is only one DNA chain.22 Limited by computer time, the selfassembly of such a structure is still out of the reach of MD simulations. More simplified and efficient water-free CG simulations have been employed to mimic the self-assembly of CL-DNA.23,24 The same group later on applied the Nogushi-Takasu implicitsolvent CG model to obtain both lamellar and inverted hexagonal phases of the complexes,25 but the water-free models with the assumption that all the counterions are completely released from the isoelectric complexes miss the hydrodynamics modes, which are important in soft matter materials. Another low-level CG MD model containing both water and counterions26 was then performed, but the system size is still limited to having only one DNA strand and 256 lipid molecules. In the present work, we use the DPD method with implementation of electrostatic interactions (coded by the authors) to simulate the self-assembly of the sandwiched structure of 8 DNA chains and up to 4000 lipids in a larger system. The advantage of this explicit solvent counterion model enables us to observe the dynamic process of the complex formation. We found it is DNA and lipid counterion release followed by the process of CLs driven to DNA strands. Although there are still 30% DNA counterions and 10% CL counterions trapped in between the two bilayers, the counterions together with the charges of DNA and CL keep the interior of the complexes neutral, as expected. We then perform simulations in a constant pressure ensemble to investigate the DNA spacing as a function of the CL concentration for various CL/DNA charge ratios. In the plot of the DNA interaxial spacing versus the inversed CL concentration, both linear and nonlinear regions are found. Simulation Method In DPD simulations, the unit particle or bead is a small fluid volume with diameter r0 containing several atoms or molecules. In the CL-DNA system, water is represented by a single bead. A lipid molecule is a polymer with three hydrophilic beads and two hydrophobic chains, each one containing four beads, as shown in Figure 1a. For a cationic lipid, one of its head bead carries one negative charge. λ-DNA, as investigated in experi-
10.1021/jp102115m 2010 American Chemical Society Published on Web 05/10/2010
7262
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Figure 1. (a) Coarse-grained models for the lipid and DNA strand. The lipid molecule has two hydrophobic chains (yellow) and three head beads (red for neutral lipid and blue for cationic lipid); the DNA strand is a rigid rod with charges (brown) uniformly distributed on the surface of the rod. (b) Initial configuration of eight DNA strands with tightly bound counterions. (c) Initial condition of lipid-DNA system. (d-f) Process of self-assembly of lamellar lipid-DNA complex. Water and counterions are not shown in c-f for clarity. 14
ments, is a semiflexible polymer with a consistent length between 50 and 100 nm, so it is modeled as a rigid, cylindrical rod with two rings of beads in our study. The diameter of the DNA chain is around 2.3 nm. Each DNA chain carries charges on the surface of the cylinder with uniform axial density λ ) 6 e/nm. Both lipid counterions and DNA counterions are single beads that carry one positive or negative charge, respectively. All the beads and ions are assumed to have the same size and mass. All the beads interact via a short-ranged repulsive force, random force, and dissipative force.2 For lipid molecules and DNA chains, beads are also connected by a harmonic force and stiffness force. The force parameters for water, lipid, and counterions are as same as in ref 11. The DNA beads interact with other beads as well as water beads. Electrostatic interactions between charged beads are calculated by using the method introduced by Groot,12 in which the charges are distributed on a lattice and the electrostatic field is solved locally on the grid. The force and the virial contributed by the self-energy term are subtracted from the total force of each charged bead and the total virial, as discussed in detail in our recent paper.13 The values of the electrostatic parameters used here are also the same as in refs 11 and 13. Initially, eight DNA chains with even spacing are fixed in the middle x-y plan with a DNA helix parallel to the y-axis. The box in the y and z directions has length Ly ) 23.04 nm and Lz ) 17.28 nm. Lx ) Nlipidalipid/4Ly is determined by the number of total lipid molecules, Nlipid, and the area per lipid, alipid. Here, alipid is set slightly larger than 0.65nm2/lipid, which is the area for the lipid bilayer freely standing in water, for a given Nlipid by trial so that a sandwiched CL-DNA complex can form in a certain range of Nlipid. DNA counterions and water are randomly distributed in the simulation box. The bead density of the system is F ) 3/r03. Before lipids and their counterions are added into the system, simulations are prerun for 0.5 µs so that the DNA counterions are attracted to the DNA strands (see Figure 1b). Then the lipid molecules and lipid counterions are randomly added into the system. Simulations are performed in the NVT ensemble for three conditions with cationic lipid/DNA (CL/DNA) charge ratios, RCL/DNA, of 1.3, 1.5, and 1.7, respectively. For each RCL/DNA, the number of charged lipids is fixed. We vary the number of neutral lipids so that the total number of lipids, Nlipid, varies from 2000 to 4000. Five independent samples are simulated at each data point to give good statistics. After well-ordered sandwiched CL-DNA complexes are obtained, simulations in a NPxLyPzT ensemble13,27-29 are per-
Gao et al.
Figure 2. (a) Density distributions for all kinds of beads. (b) Density distribution for DNA counterions before lipids are added into the system and after the lamellar phase forms.
formed for 4 µs. Here, we set the pressure to Px ) Pz ) 23kBT/ r03 ) 93.3 atm at room temperature, which is the pressure for a pure water box at F ) 3/r03. The box size in the DNA strand direction (y-axis) is fixed. In such an ensemble, the positions of the DNA strands follow the change of the box in the x direction but still with an even interaxial separation. Involving water and counterions in the DPD simulations, a preliminary simulation of the self-assembly process in the NVT ensemble is essential because we do not know how many CLs, waters, and counterions are trapped inside the complexes. Compared with a solvent-free model,23,24 the self-assembly process elongates the simulation time, but in another way, it is an advantage of DPD that it can mimic more realistic hydrodynamics of soft matter on a feasible time scale. Results and Discussions The self-assembly of CL-DNA complexes is illustrated in Figure 1c-f, where the system has 8 DNA chains and 3000 lipids, with 1944 of them negatively charged; that is, RCL/DNA ) 1.5. In up to 10 µs, a well-ordered sandwiched structure of CL-DNA is formed. In this process, DNA counterions, which are initially bound to the DNA strands, are slowly released into the bulk. At the same time, 1050 cationic lipid molecules are attracted by anionic DNA to reside in the two inner leaflets facing the DNA array. For the configuration in Figure 1f, the density distribution of each kind of bead along the z-axis is plotted in Figure 2a. The sharply peaked distribution of lipid head beads demonstrates that the lipids assemble to two wellstructured bilayers. The distance between the two inner leaflets is about 2.7 nm, which is just enough to accommodate a DNA strand of 2.3 nm diameter. The lipids in the inner monolayer are more confined, since the attraction between the DNA and the CLs suppresses the membrane fluctuations. In Figure 2b, we also give the distribution of DNA counterions for both the lipid-free system and the sandwiched CL-DNA structure. We see that most DNA counterions are released into the water reservoir outside the complex. Only around 310 DNA counterions (less than 30%) remain inside the complex between the two lipid bilayers. The formation of well-ordered sandwiched CL-DNA structures is observed only when the CL/DNA charge ratio is greater than 1.2 and the number of lipid molecules is in the range from 2000 to 3600. At a lower charge ratio, such a structure has less possibility to form because not enough CL can be attracted to the DNA strands to form bilayers. The lipids may form vesicles, micelles, or freely standing bilayers far away from the DNA strands. On the basis of the assumption that the counterions
Self-Assembly of Lamellar Lipid-DNA Complexes
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7263
Figure 3. Top view of (a) a complex whose bilayer has a hole for Nlipid ) 1944 and (b) a disordered complex for Nlipid ) 3800 at RCL/DNA ) 1.5.
can be completely released from the lipid and DNA, an isoelectric complex at RCL/DNA ) 1 is the most stable structure, as predicted by statistical thermodynamic theory20 and implicitsolvent calculation.23-25 It is also verified by experiments,14 but we found that in up to 10 µs in our simulations, the counterions are not completely released. Thus, the electrostatic attraction between the charged lipids and the DNA strands is weaker than that at the isoelectric point. To obtain stable complexes in this time scale, more charged lipids are needed to provide a stronger attractive driving force and entropy associated with ion releasing. It might be possible to obtain a well ordered isoelectric complex in which more or all counterions are released if very long simulations are performed, but it is thus far out of the reach of our computer performance. Nevertheless, our observation also indicates that isoelectricity is not essential for the formation of stable, sandwiched CL-DNA complexes. At a higher CL/DNA charge ratio, the sandwiched structure forms even faster and is stable. The stability will be further discussed in the following. At a given number of DNA chains (say, eight), if there are not enough lipids (for example Nlipid < 2000), holes are observed on one or both of the bilayers, no matter how close the DNA strands are. The DNA strands tend to escape via the holes (see Figure 3a); however, if Nlipid > 3600, disorder complexes rather than sandwiched structures form at any DNA spacing, as shown in Figure 3b. The range of lipid number in which a sandwiched CL-DNA complex can form may be smaller than that in the implicit solvent model because here, the counterions are not completely released, which screens the attraction between the CLs and the DNA strands. For a stable sandwiched CL-DNA complex, the total charge between the two bilayers should be zero. Such a property is proved in Figure 4, where we plot the number of cationic lipids in ; the number of on the inner leaflets of the two bilayers, NCL in DNA counterions, NDNA-counterion; and the cationic lipid counin terions, NCL-counterion , trapped inside the complex as a function of the total number of lipids at a CL/DNA charge ratio of RCL/DNA ) 1.5. Data are the average of five independent samples at each given lipid number. When there are more lipids in the system, the CLs attracted to the inner leaflet of the bilayers increase slightly because at larger DNA spacing, the DNA counterion screening effects become weaker. There are also weak fluctuations for CL and DNA counterions, as shown in Figure 4b and c. But the total charges, Nincharge, as given in Figure 4d, are very close to zero at any Nlipid. The interaxial distance, dDNA, as a function of the total number of lipid molecules simulated in a NPT ensemble is given in Figure 5a for RCL/DNA of values of 1.3, 1.5, and 1.7. We found the dDNA increases almost linearly as a function of the number
Figure 4. The number of (a) cationic lipids, (b) CL counterions, (c) DNA counterions, and (d) net charge trapped inside the CL-DNA complexes as a function of the number of lipids for a CL/DNA charge ratio of RCL/DNA ) 1.5.
Figure 5. (a) DNA spacing versus the total number of lipids. (b) DNA in spacing as a fucntion of 1/φin c with φc , the molar ratio of cationic lipids to neutral lipids plus cationic lipids on the inner monolayers.
of lipids when Nlipid > 2800, and is independent of RCL/DNA. An isoelectric complex model,15 which assumes complete counterion release and neutralization of anionic DNA and CL, predicted that dDNA is related to the number ratio of the charged lipids to the total lipids in the inner monolayer, φc, and the area per lipid, alipid, by
dDNA )
alipidλDNA 1 2e φc
(1)
Although the DNA and CL counterions are not completely released in our simulations, the neutralization in the interior of the complexes keeps such linear relation still approximately valid, as shown in Figure 5b. At low φc, dDNA increases linearly with 1/φc which means the area per lipid, alipid, is a constant. The linear fits of the three curves give alipid as 0.69, 0.71, and 0.71 nm2 for RCL/DNA ) 1.3, 1.5, and 1.7, respectively. The nonlinear behavior of dDNA versus 1/φc at higher φc means that the area per lipid varies with φc, which is given in Figure 6. The repulsion between CLs, the attraction between CLS and DNA, and the repulsions between DNA strands have effects on the lipid area, but the first effect is very weak. We already found that for a freely standing bilayer, the area per lipid as a function of φc is given by alipid ) 0.65 + 0.008eφc/0.482 nm2 (see Figure 6, the solid line), so the effect is negligible. The attraction between CLs and DNA strands tends both to stretch and to
7264
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Gao et al. can be applied to many biological systems and soft matters in which electrostatics and hydrodynamic modes play important roles. Acknowledgment. We thank L. Golubovic for useful discussions. This work is supported by the National Science Foundation of China (Grant No. 20873007) and the Major State Basic Research Development Programs (Grant No. 2004CB719903). References and Notes
Figure 6. The averaged area per lipid for lipids on the inner monolayers as a function of φc.
undulate the membrane. The repulsion between DNA strands tends to stretch the membrane, too. We note that the existence of the DNA chains, especially when they are highly compacted at high CL concentration, suppresses the membrane undulation, so the electrostatic stress on the membrane is mainly balanced by the surface tension, which results in a larger alipid. When φc decreases, the electrostatic interactions becomes weak; thus, the stretching stress as well as the lipid area, alipid, decrease very quickly. We found the asymptotic values of alipid for the three sets of CL/DNA charge ratios are all close to 0.69 nm2, in good agreement with the fitting results. From these discussions, it is also understandable that DNA spacing as a function of Nlipid is independent of charge ratio RCL/DNA: At low RCL/DNA, fewer CLs are attracted to the inner monolayers so that 1/φc is larger, which should lead to a larger dDNA. At the same time, the electrostatic attraction between the DNA and CLs is weak so that the bilayers are less stretched, which should result in a smaller dDNA. The combination of these two effects gives rise to the same dDNA for a given Nlipid at various RCL/DNA’s. Conclusion In summary, the dissipative particle dynamics simulations with explicit solvent and counterions successfully describes the self-assembly of lamellar cationic lipid-DNA complexes. In a certain range of CL concentration, a well-ordered sandwiched structure can form. The formation of the complexes is associated with the releasing of DNA and lipid counterions. In the complex phase, the bilayers are tensed to balance the electrostatic stress induced by DNA-CL attraction and DNA-DNA repulsions so that the area per lipid as well as the interaxial DNA spacing are dependent on the CL concentration. This mesoscale model
(1) Espanol, P.; Warren, P. B. Europhys. Lett. 1995, 30, 191. (2) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423. (3) Shillcock, J. C.; Lipowsky, R. J. Chem. Phys. 2002, 117, 5048. (4) Kranenburg, M.; Nicolas, J.; Smit, B. Phys. Chem. Chem. Phys. 2004, 6, 4142. (5) Gao, L.; Shillcock, J. C.; Lipowsky, R. J. Chem. Phys. 2007, 126, 015101. (6) Ortiz, V.; Nielsen, S. O.; Discher, D. E.; Klein, M. L.; Lipowsky, R.; Shillcock, J. J. Phys. Chem. B 2005, 109, 17708. (7) Shillcock, J. C.; Lipowsky, R. Nat. Mater. 2005, 4, 225. (8) Grafmu¨ller, A.; Shillcock, J.; Lipowsky, R. Phys. ReV. Lett. 2007, 98, 218101. (9) Gao, L.; Lipowsky, R.; Shillcock, J. C. Soft Matter 2008, 4, 1208. (10) Venturoli, M.; Smit, B.; Sperotto, M. M. Biophys. J. 2005, 88, 1778. (11) Gao, L.; Fang, W. Soft Matter 2009, 8, 3312. (12) Groot, R. D. J. Chem. Phys. 2003, 118, 11265. (13) Gao, L.; Fang, W. J. Chem. Phys. 2010, 132, 031102. (14) Ra¨dler, J. O.; Koltover, I.; Salditt, T.; Safinya, C. R. Science 1997, 275, 810. (15) Koltover, I.; Salditt, T.; Ra¨dler, J. O.; Safinya, C. R. Biophys. J. 1999, 77, 915. (16) Koltover, I.; Salditt, T.; Ra¨dler, J. O.; Safinya, C. R. Science 1998, 281, 78. (17) Golubovic, L.; Golubovic, M. Phys. ReV. Lett. 1998, 80, 4341. (18) O’Hern, C. S.; Lubensky, T. C. Phys. ReV. Lett. 1998, 80, 4345. (19) Bruinsma, R. F. Euro. Phys. J. B 1998, 4, 75. (20) Harries, D.; May, S.; Gelbert, W. M.; Ben-Shaul, A. Biophys. J. 1998, 75, 159. (21) Fleck, C.; Netz, R. R.; von Gruberg, H. H. Biophys. J. 2002, 76, 76. (22) Bandyopadhyay, S.; Tarek, M.; Klein, M. L. J. Phys. Chem. B 1999, 103, 10075. (23) Farago, O.; Gronbech-Jensen, N.; Pincus, P. Phys. ReV. Lett. 2006, 96, 018102. (24) Farago, O.; Gronbech-Jensen, N. Biophys. J. 2007, 92, 3228. (25) Farago, O.; Gronbech-Jensen, N. J. Am. Chem. Soc. 2009, 131, 2875. (26) Khalid, S.; Bond, P. J.; Holyoake, J.; Hawtin, R. W.; Sansom, M. S. P. J. R. Soc. Interface 2008, 5, s241. (27) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177. (28) Jakobsen, A. F. J. Chem. Phys. 2005, 122, 124901. (29) Hu¨nenberger, P. H. J. Chem. Phys. 2002, 116, 6880.
JP102115M
