Computational Biochemistry · Membrane Protein Simulation
Membrane Protein: KALP₁₅ in DPPC Bilayer
Transmembrane peptide · Berger lipids · GROMOS96 53A6 · InflateGRO packing · Membrane-specific analysis
GROMACS
GROMOS96 53A6
Berger Lipids
DPPC Bilayer
InflateGRO · GridMAT-MD
System
Prep
Prep
Prepare Topology
Step 1
pdb2gmx
Modify Force Field
Step 2
Berger lipids
Build & Pack Bilayer
Step 3
InflateGRO
Solvate & Add Ions
Step 4
solvate · genion
Equil. ·
Prod. · Analysis
Prod. · Analysis
Energy Minimization
Step 5
mdrun EM
NVT + NPT Equil.
Step 6
323 K · 1 ns
Production MD
Step 7
1 ns
Membrane Analysis
Step 8
SCD · density · APL
1
System Preparation
▼
System Overview
KALP₁₅ is a model transmembrane peptide with sequence Ac-GKK(LA)₄LKKA-NH₂, oriented with its helix axis along the z-axis. The bilayer is 128-lipid DPPC, lying in the x-y plane. This tutorial is based on Kandasamy & Larson's study of hydrophobic mismatch.
Orient the Peptide
The KALP peptide is built with ideal α-helix geometry (φ = −60°, ψ = −40°) and oriented along the z-axis using editconf -princ. The pre-oriented structure is in KALP-15_princ.pdb.
Generate Topology
gmx pdb2gmx -f KALP-15_princ.pdb -o KALP-15_processed.gro -ignh -ter -water spc
# Select: GROMOS96 53A6 · "None" for both termini (acetyl/amide caps already present)
# Select: GROMOS96 53A6 · "None" for both termini (acetyl/amide caps already present)
- -ignh — ignore H atoms in input; the xLeap (AMBER) H naming differs from GROMOS96 — pdb2gmx will add only what it needs
- -ter — interactively assign termini; choose "None" for both N- and C-termini so pdb2gmx connects to the acetyl and amide capping groups
- GROMOS96 53A6 — compatible with Berger lipid parameters (which will be added in Step 2)
Note (2023): The original Berger lipid parameters are no longer available online. This tutorial serves as a methodological reference. For modern membrane simulations, CHARMM36 with CHARMM-GUI is recommended.
▼
Why a Custom Force Field?
The Berger lipids are a hybrid of GROMOS atom types and OPLS partial charges. They cannot be simply #included alongside forcefield.itp — they must be merged into the force field files themselves. The Ryckaert-Bellemans (RB) dihedral potential is used for alkane chains, which differs from standard GROMOS periodic dihedrals.
Step 1 — Copy the Force Field Directory
cp -r /usr/local/gromacs/share/gromacs/top/gromos53a6.ff/ ./gromos53a6_lipid.ff
Step 2 — Merge Parameters (6 required changes)
# 1. Copy [atomtypes] from lipid.itp → ffnonbonded.itp (add atomic number column)
# 2. Copy [nonbond_params] from lipid.itp → ffnonbonded.itp
# 3. DELETE ";; parameters for lipid-GROMOS interactions" and ALL subsequent lines
# 4. Remove/rename all HW lines in [nonbond_params] → will cause fatal error if left
# 5. Copy [pairtypes] from lipid.itp → ffnonbonded.itp
# 6. Append [dihedraltypes] from lipid.itp → ffbonded.itp
New atom type entries in ffnonbonded.itp (with atomic numbers added):
LO 8 15.9994 0.000 A 2.364e-03 1.590e-06 ;carbonyl O, OPLS
LP2 6 14.0270 0.000 A 5.874e-03 2.265e-05 ;RB CH2, Berger LJ
LP3 6 15.0350 0.000 A 8.777e-03 3.385e-05 ;RB CH3, Berger LJ
Step 3 — Update topol.top
#include "gromos53a6_lipid.ff/forcefield.itp" ← changed from gromos53a6.ff
...
#ifdef POSRES
#include "posre.itp"
#endif
#include "dppc.itp" ← ADD: DPPC molecule topology
#include "gromos53a6_lipid.ff/spc.itp"
Scaling factor: Berger lipids use a 1-4 LJ scaling factor of 0.125 (vs 1.0 in standard GROMOS). This is handled internally by the RB dihedral parameters and is critical for correct acyl chain behavior.
▼
1 — Orient the Peptide in the Bilayer Frame
# Remove periodicity from DPPC bilayer
gmx grompp -f minim.mdp -c dppc128.pdb -p topol_dppc.top -o dppc.tpr
gmx trjconv -s dppc.tpr -f dppc128.pdb -o dppc128_whole.gro -pbc mol -ur compact
# Place KALP at bilayer center using the DPPC box vectors
gmx editconf -f KALP-15_processed.gro -o KALP_newbox.gro -c -box 6.41840 6.44350 6.59650
# Concatenate protein + bilayer
cat KALP_newbox.gro dppc128_whole.gro > system.gro
# Remove duplicate headers/box vectors, update atom count on line 2
2 — Pack Lipids with InflateGRO
InflateGRO first inflates the bilayer laterally (scale = 4) to create space around the protein, removes overlapping lipids, then shrinks iteratively (scale = 0.95) with energy minimization at each step until the area per lipid reaches the target value.
# Generate strong protein position restraints
gmx genrestr -f KALP_newbox.gro -o strong_posre.itp -fc 100000 100000 100000
# Inflate lipids (scale factor = 4, delete lipids within 14 Å of protein)
perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat
# EM on inflated system
gmx grompp -f minim_inflategro.mdp -c system_inflated.gro -p topol.top -r system_inflated.gro -o system_inflated_em.tpr
gmx mdrun -deffnm system_inflated_em
# Shrink iteratively: run ~26 times with scale=0.95 until APL ≈ 71 Ų
perl inflategro.pl system_inflated_em.gro 0.95 DPPC 0 system_shrink1.gro 5 area_shrink1.dat
Inflate factor
4.0
Shrink factor
0.95 × 26
Target APL
~71 Ų
Exp. APL DPPC
62 Ų
3 — Solvate & Remove Core Water
gmx solvate -cp system_shrink26_em.gro -cs spc216.gro -o system_solv.gro -p topol.top
# Remove water molecules that filled the hydrophobic bilayer core
perl water_deletor.pl -in system_solv.gro -out system_solv_fix.gro -ref O33 -middle C50 -nwater 3
# Update SOL count in topol.top with number reported by script
Water in acyl chains: gmx solvate fills ALL voids including the hydrophobic bilayer interior. These water molecules MUST be removed before continuing —
water_deletor.pl uses upper/lower leaflet ester atoms (O33) as reference and the bilayer midpoint (C50) to define the deletion zone.▼
gmx grompp -f ions.mdp -c system_solv_fix.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o system_solv_ions.gro -p topol.top -pname NA -nname CL -neutral
The KALP₁₅ peptide contains 4 lysine residues bearing a net charge of +4e at neutral pH. Four Cl⁻ counterions are added to neutralize the system.
Peptide charge
+4e (4 Lys)
Cl⁻ added
4
Net charge
0
DPPC molecules
128 − deleted
2
Equilibration
▼
gmx grompp -f minim.mdp -c system_solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
Membrane-Specific EM Issues
Membrane systems are far more challenging to minimize than simple aqueous proteins. Common problems and solutions:
- Intra-headgroup H-bonding collapse (PE/PG) — headgroups fold into solvent voids. Fix: use position restraints or freeze groups; temporarily reduce H-atom charges (restore before production!)
- Acyl chain overlap — occurs during InflateGRO packing. Use appropriate cutoff in initial inflation step; do not over-pack
- Protein-lipid overlap — check the InflateGRO cutoff radius (14 Å used here)
- Ion-headgroup overlap — a Cl⁻ placed next to a phosphate can cause explosive forces. genion random placement can be problematic near lipid headgroups
Potential Energy — Energy Minimization
KALP₁₅/DPPC Complex · Steepest Descent · GROMOS96 53A6 + Berger
Convergence confirms system is ready for dynamics — verify Epot is negative and Fmax < 1000 kJ mol⁻¹ nm⁻¹
▼
Create Special Index Groups
gmx make_ndx -f em.gro -o index.ndx
> 1 | 13 ← merge Protein (1) + DPPC (13) → Protein_DPPC
> q
NVT Equilibration
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt
Key NVT .mdp Differences from Aqueous Protein
- T = 323 K — must be above DPPC phase transition temperature (315 K); using 308 K would freeze the bilayer into a gel phase
- tc-grps = Protein_DPPC Water_and_ions — couple protein+lipids together; DO NOT couple ions separately (insufficient DoF)
- COM motion removal — comm-grps = Protein_DPPC Water_and_ions separately, preventing lateral drift of bilayer vs water phases in opposite directions
- rcoulomb = rvdw = 1.2 nm — compromise between GROMOS96 (0.9 nm) and Berger original (1.0 nm) parametrization
Temperature — NVT Equilibration
KALP₁₅/DPPC · 323 K target · V-rescale thermostat
Temperature stabilizes at 323 K — above DPPC gel–liquid phase transition (315 K)
NPT Equilibration — Membrane-Specific Settings
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -deffnm npt
- tcoupl = Nosé-Hoover — produces correct kinetic ensemble for NPT; NOT appropriate for NVT (too wide fluctuations)
- pcoupltype = semiisotropic — bilayer deforms independently in x-y vs z; isotropic would suppress membrane area fluctuations
- Two compressibility values — x-y plane (lateral) and z (bilayer normal) treated independently
Pressure
KALP₁₅/DPPC · NPT
Running avg stabilizes near 1 bar
Box Vectors (X & Y)
Membrane lateral dimensions
Stable Box-X/Y confirm membrane lateral area equilibrated
3
Production MD & Membrane Analysis
▼
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_1.tpr
gmx mdrun -deffnm md_0_1
All position restraints are absent from the production .mdp file. The system continues with Nosé-Hoover thermostat and semiisotropic pressure coupling established during NPT equilibration.
1 ns is insufficient for lipid equilibration. Lipid lateral diffusion and membrane reorganization require 50–100 ns of simulation. This tutorial provides instructional value only — results at 1 ns have no rigorous physical meaning. Additionally, GROMOS96 53A6 under-stabilizes α-helices; consider 54A7 for long simulations.
▼
1. Deuterium Order Parameters (−SCD)
Order parameters measure the orientational order of acyl chain C-H bonds relative to the bilayer normal. A gel-phase bilayer shows high order; a liquid-crystalline bilayer shows lower, more uniform values. Each chain must be analyzed separately using a dedicated index group.
# Create sn-1 chain index (C34, C36-C50 by name)
gmx make_ndx -f md_0_1.tpr -o sn1.ndx
> a C34
> a C36
... a C50
> del 0-21
> q
# Calculate order parameters along z-axis
gmx order -s md_0_1.tpr -f md_0_1.xtc -n sn1.ndx -d z -od deuter_sn1.xvg
# Repeat for sn-2 chain (C15, C17-C31)
Deuterium Order Parameters (−SCD)
KALP₁₅/DPPC · sn-1 and sn-2 palmitoyl chains
Lower −SCD at chain ends (terminal methyl) is expected. High plateau values confirm liquid-crystalline (non-gel) phase at 323 K
2. Membrane Density Profile
Density is analyzed by splitting DPPC into three groups: headgroups, glycerol/ester region, and acyl chains. Together with protein and water, this gives a complete picture of bilayer structure along z.
gmx make_ndx -f md_0_1.tpr -o density_groups.ndx
> 13 & a C1 | a C2 | a C3 | a N4 | ... | a O11
> name 22 Headgroups
> 13 & a C12 | a C13 | a O14 | a C15 | a O16 | a C32 | a O33 | a C34 | a O35
> name 23 Glycerol_Ester
> 13 & ! 22 & ! 23
> name 24 Acyl_Chains
gmx density -s md_0_1.tpr -f md_0_1.xtc -n density_groups.ndx -o dens_headgroups.xvg -d Z
Membrane Electron Density Profile
KALP₁₅/DPPC · Along z-axis · headgroups · glycerol · acyl · water · protein
Characteristic bilayer profile: water peaks flank the headgroup regions; low-density acyl core in the center
3. Area per Lipid & Bilayer Thickness — GridMAT-MD
Standard GROMACS tools cannot compute area per lipid in the presence of an embedded protein (the protein displaces lipids). GridMAT-MD solves this by projecting a 2D grid onto the bilayer plane and counting headgroup occupancy per grid cell, excluding protein-occupied regions.
4. Lateral Diffusion of Lipids
gmx make_ndx -f md_0_1.tpr -o p8.ndx
> a P8 ← one reference atom per DPPC headgroup
> q
gmx msd -s md_0_1.tpr -f md_0_1.xtc -n p8.ndx -lateral z
Mean-Square Displacement — Lateral Diffusion
DPPC headgroup P8 atoms · lateral (x-y plane)
D(P8) ≈ 0.027 × 10⁻⁵ cm²/s at 1 ns — 1 ns is insufficient for accurate D; longer simulations (50+ ns) required for convergence
🧾
Common Lipid Parameters
| Lipid | Area per Lipid (Ų) | Phase Transition (K) | Reference |
|---|---|---|---|
| DPPC | 62.9 – 64 | 315 | Nagle (1993) Biophys. J. 64: 1476 |
| DMPC | 60.6 | 297 | Wohlert & Edholm (2006) J. Chem. Phys. 125: 204703 |
| POPC | 65.8 | 271 | Tieleman et al. (1998) Biochem. 37: 17554 |
| POPE | 56 | 298 | Tieleman et al. (1998) Biochem. 37: 17554 |
| POPG | 53 | 269 | Dickey & Faller (2008) Biophys. J. 95: 2636 |
| POPA | 51–52 | 301 | Dickey & Faller (2008) Biophys. J. 95: 2636 |
| POPS | 55 | 300 | Mukhopadhyay et al. (2004) Biophys. J. 86: 1601 |
| DMTAP | 71 | 310 | Gurtovenko et al. (2004) Biophys. J. 86: 3461 |
⚙
Simulation Parameters
Protein FF
GROMOS96 53A6
Lipid FF
Berger (UA)
Water model
SPC
Temperature
323 K
Thermostat
Nosé-Hoover
Pressure type
Semiisotropic
rcoulomb/rvdw
1.2 nm
DPPC lipids
128 − deleted
📈
Membrane Analysis Metrics
−SCD
Deuterium order params
sn-1 & sn-2 chains
sn-1 & sn-2 chains
Density
Electron density profile
along bilayer normal
along bilayer normal
APL
Area per lipid
GridMAT-MD
GridMAT-MD
Thickness
Bilayer thickness
z-projection
z-projection
Dlateral
Lipid diffusion
gmx msd -lateral z
gmx msd -lateral z
RMSD
Peptide stability
P-N vector, tilt
P-N vector, tilt