- Posts: 4
DPPC/DUPC/CHOL: Lo phase has gel-like properties
- da294
-
Topic Author
- Offline
- Fresh Boarder
Less
More
10 years 1 month ago #3544
by da294
DPPC/DUPC/CHOL: Lo phase has gel-like properties was created by da294
Hello Everyone,
I have simulated a three component lipid bilayer of DPPC/DUPC/CHOL = 40/40/20 using the martini_v2.1.itp forcefield as updated in 2007. The simulation contains 4608 lipids and was run for 25 microseconds at 295 K. The bilayer does indeed phase separate into a DPPC/CHOL rich "Lo" phase and a DUPC rich "Ld" phase. However, the Lo phase seems to have gel-like properties.
When measuring the correlation functions of, for instance, cholesterol to cholesterol in this Lo phase, I see persistent correlations that last upwards of 10 nm. Additionally, the diffusion in the Lo phase is extremely slow. When I look only at Lo lipids >3 nm from the Lo-Ld phase interface, and I correct for their average COM motion, I get diffusion coefficients on the order of 5E-10 cm^2/s. This is nearly 100 times slower than reported elsewhere. If I only correct for the overall leaflet COM motion, I get larger diffusion coefficients. However, I believe this might be misleading as the lipids deep within the Lo phase move in a concerted way, which when subtracted, reveals their almost non-existent diffusion with respect to one another.
Preliminary results on shorter simulations of smaller systems do not seem to have this slow diffusion, but I am still looking into that.
I was just wondering if anyone else has come across this, and whether or not this is an artifact.
I have copied my .mdp file below for clarity.
Thank you very much for your time,
David
________________________________________________________________________
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require
20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp
file.
integrator = md
dt = 0.02
nsteps = 1250000000
nstcomm = 10
comm-grps =
nstxout = 250000
nstvout = 250000
nstfout = 0
nstlog = 250000
nstenergy = 250000
nstxtcout = 0
xtc_precision = 0
xtc-grps =
energygrps = DPPC DUPC CHOL W
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy
conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit
screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a
lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with
Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf =
infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale
(V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = DPPC DUPC CHOL W
tau_t = 1.0 1.0 1.0 1.0
ref_t = 295 295 295 295
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 4.0 4.0 ;parrinello-rahman is more stable with
larger tau-p, DdJ, 20130422
compressibility = 5e-5 5e-5
ref_p = 1.0 1.0
gen_vel = yes
gen_temp = 295
gen_seed = -1
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
I have simulated a three component lipid bilayer of DPPC/DUPC/CHOL = 40/40/20 using the martini_v2.1.itp forcefield as updated in 2007. The simulation contains 4608 lipids and was run for 25 microseconds at 295 K. The bilayer does indeed phase separate into a DPPC/CHOL rich "Lo" phase and a DUPC rich "Ld" phase. However, the Lo phase seems to have gel-like properties.
When measuring the correlation functions of, for instance, cholesterol to cholesterol in this Lo phase, I see persistent correlations that last upwards of 10 nm. Additionally, the diffusion in the Lo phase is extremely slow. When I look only at Lo lipids >3 nm from the Lo-Ld phase interface, and I correct for their average COM motion, I get diffusion coefficients on the order of 5E-10 cm^2/s. This is nearly 100 times slower than reported elsewhere. If I only correct for the overall leaflet COM motion, I get larger diffusion coefficients. However, I believe this might be misleading as the lipids deep within the Lo phase move in a concerted way, which when subtracted, reveals their almost non-existent diffusion with respect to one another.
Preliminary results on shorter simulations of smaller systems do not seem to have this slow diffusion, but I am still looking into that.
I was just wondering if anyone else has come across this, and whether or not this is an artifact.
I have copied my .mdp file below for clarity.
Thank you very much for your time,
David
________________________________________________________________________
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require
20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp
file.
integrator = md
dt = 0.02
nsteps = 1250000000
nstcomm = 10
comm-grps =
nstxout = 250000
nstvout = 250000
nstfout = 0
nstlog = 250000
nstenergy = 250000
nstxtcout = 0
xtc_precision = 0
xtc-grps =
energygrps = DPPC DUPC CHOL W
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy
conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit
screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a
lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with
Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf =
infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale
(V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = DPPC DUPC CHOL W
tau_t = 1.0 1.0 1.0 1.0
ref_t = 295 295 295 295
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 4.0 4.0 ;parrinello-rahman is more stable with
larger tau-p, DdJ, 20130422
compressibility = 5e-5 5e-5
ref_p = 1.0 1.0
gen_vel = yes
gen_temp = 295
gen_seed = -1
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
Please Log in or Create an account to join the conversation.
- mnmelo
-
- Offline
- Admin
Less
More
- Posts: 89
10 years 1 month ago #3554
by mnmelo
Replied by mnmelo on topic DPPC/DUPC/CHOL: Lo phase has gel-like properties
Hi David,
Thank you very much for the report. Others have also been letting us know of similar observations, and I too got an apparently over-ordered lo phase when rerunning some of the mixtures from the Risselada 2008 PNAS paper.
Your .mdp looks very similar to the one I used. The only difference is that I tried to emulate older GROMACS' behavior by setting nsttcouple and nstpcouple to 1 (which didn't really affect lo phase behavior).
We're currently looking into the 2008 trajectory files to try and determine whether this over-ordering has always been present or whether it is somehow dependent on the GROMACS version or particular starting conditions. So far it seems Martini cholesterol fails to efficiently disorder DPPC packing, becoming itself part of the gel lattice. This is in contrast to what was recently observed in microseconds-long all-atom simulations (Sodt et al. "The Molecular Structure of the Liquid Ordered Phase of Lipid Bilayers." JACS, 2013).
We'll keep people posted of our findings, including if we need improved topologies to better reproduce the properties of lo phases.
Cheers, and thanks again for reporting,
Manuel Melo
Thank you very much for the report. Others have also been letting us know of similar observations, and I too got an apparently over-ordered lo phase when rerunning some of the mixtures from the Risselada 2008 PNAS paper.
Your .mdp looks very similar to the one I used. The only difference is that I tried to emulate older GROMACS' behavior by setting nsttcouple and nstpcouple to 1 (which didn't really affect lo phase behavior).
We're currently looking into the 2008 trajectory files to try and determine whether this over-ordering has always been present or whether it is somehow dependent on the GROMACS version or particular starting conditions. So far it seems Martini cholesterol fails to efficiently disorder DPPC packing, becoming itself part of the gel lattice. This is in contrast to what was recently observed in microseconds-long all-atom simulations (Sodt et al. "The Molecular Structure of the Liquid Ordered Phase of Lipid Bilayers." JACS, 2013).
We'll keep people posted of our findings, including if we need improved topologies to better reproduce the properties of lo phases.
Cheers, and thanks again for reporting,
Manuel Melo
Please Log in or Create an account to join the conversation.
- da294
-
Topic Author
- Offline
- Fresh Boarder
Less
More
- Posts: 4
10 years 1 month ago #3556
by da294
Replied by da294 on topic DPPC/DUPC/CHOL: Lo phase has gel-like properties
Thank you for your prompt response.
In case it helps, I have also been testing different systems for these effects. I have seen that the strong and persistent correlations exist in the Lo phase of DPPC/DUPC/Chol bilayers regardless of using martini_v2.0.itp, martini_v2.1.itp (updated as of 2007), or martini_v2.1.itp (updated in 2011). The correlations also exist in pure Lo phases of just DPPC/Chol bilayers at molar ratios of 2:1. Again, for large system sizes, these correlations can extend upwards of 10 nm.
The extremely slow diffusion, on the other hand, is not always as easily apparent. However, discrepancies appear in larger simulations containing ~4600 lipids. In such simulations, the Lo lipids in the bulk of the Lo phase (more than several nm from the Lo-Ld interface) move in a concerted way. Over 1 microsecond, they have a mean squared displacement of < 0.5 nm^2, when correcting for their average COM motion. This leads to a diffusion coefficient nearly 100 times slower than what it should be. I also see similarly slow diffusions in pure Lo bilayers.
I believe that in the smaller simulation sizes, the Lo-Ld interface makes up a significant amount of the Lo phase and thereby enables faster diffusion of Lo lipids into the Ld phase. In the larger simulations, Lo lipids can be several nm from an interface, meaning they do not readily exchange with their surrounding over even microsecond timescales. This is probably also coupled to the overly condensed and gel-like state of the Lo phase.
Since the Martini Lo phase seems more gel-like than Lo-like, I was just wondering what part of the Martini simulation results should be taken as accurate versus an artifact of the parameterization?
I really appreciate your help.
Thanks,
David
In case it helps, I have also been testing different systems for these effects. I have seen that the strong and persistent correlations exist in the Lo phase of DPPC/DUPC/Chol bilayers regardless of using martini_v2.0.itp, martini_v2.1.itp (updated as of 2007), or martini_v2.1.itp (updated in 2011). The correlations also exist in pure Lo phases of just DPPC/Chol bilayers at molar ratios of 2:1. Again, for large system sizes, these correlations can extend upwards of 10 nm.
The extremely slow diffusion, on the other hand, is not always as easily apparent. However, discrepancies appear in larger simulations containing ~4600 lipids. In such simulations, the Lo lipids in the bulk of the Lo phase (more than several nm from the Lo-Ld interface) move in a concerted way. Over 1 microsecond, they have a mean squared displacement of < 0.5 nm^2, when correcting for their average COM motion. This leads to a diffusion coefficient nearly 100 times slower than what it should be. I also see similarly slow diffusions in pure Lo bilayers.
I believe that in the smaller simulation sizes, the Lo-Ld interface makes up a significant amount of the Lo phase and thereby enables faster diffusion of Lo lipids into the Ld phase. In the larger simulations, Lo lipids can be several nm from an interface, meaning they do not readily exchange with their surrounding over even microsecond timescales. This is probably also coupled to the overly condensed and gel-like state of the Lo phase.
Since the Martini Lo phase seems more gel-like than Lo-like, I was just wondering what part of the Martini simulation results should be taken as accurate versus an artifact of the parameterization?
I really appreciate your help.
Thanks,
David
Please Log in or Create an account to join the conversation.
- davhak
-
- Offline
- Fresh Boarder
Less
More
- Posts: 1
10 years 1 month ago #3557
by davhak
Replied by davhak on topic DPPC/DUPC/CHOL: Lo phase has gel-like properties
Dear All,
We are also simulating DPPC/DUPC/CHOL Martini systems and so far observed a regular Lo/Ld phases at, particularly, 15% of CHOL concentration [ pubs.acs.org/doi/abs/10.1021/jp312245y or www.plosone.org/article/info%3Adoi%2F10....journal.pone.0087369 ], (up to 30% of CHOL we were still observing and Lo phase for DPPC/CHOL domains) for a system size of ~20x20 nm.
We are also simulating DPPC/DUPC/CHOL Martini systems and so far observed a regular Lo/Ld phases at, particularly, 15% of CHOL concentration [ pubs.acs.org/doi/abs/10.1021/jp312245y or www.plosone.org/article/info%3Adoi%2F10....journal.pone.0087369 ], (up to 30% of CHOL we were still observing and Lo phase for DPPC/CHOL domains) for a system size of ~20x20 nm.
Please Log in or Create an account to join the conversation.
Time to create page: 0.088 seconds