Dataset DOI: 10.5061/dryad.ttdz08m92

Description of the data and file structure

Overview

This work implements the nucleation collective variable ξ proposed by Hub et. al (10.1021/acs.jctc.7b00106) to expand on our recently developed coarse-grained to atomistic backmapping and umbrella sampling methodology (10.1021/acs.jpcb.4c03276) to rapidly resolve free energy profiles for melittin-lined pores in red blood cell (RBC) lipid membranes with molecular dynamics (MD) simulations. All simulations were run with Gromacs 2021.5 patched with PLUMED 2.8.

Folder Structure

The provided tarball (Supporting_Data.tar.gz) contains the following folders to aid with the reproducibility of the results of this work:

• analysis_output --> Stores all trajectory analysis results from MD simulation trajectories
• backmapping --> Contains all necessary .map files for backmapping CG to AA configurations, including the 20 amino acids, 4 lipids (POPC, POPE, POPS, POSM), and cholesterol (CHOL) 
• ff --> All coarse-grained (martini2.2) and atomistic (charmm36) force field parameters used for simulations in this work
• gro_start --> Initial configurations for all umbrella sampling windows
• idx --> Index files (.ndx) for all systems for use with Gromacs
• long_equil --> All long timescale MARTINI simulations used to initialize MEL pore-lining configurations
• mdp --> Molecular dynamics parameter files (.mdp) for atomistic umbrella sampling simulations
• topology --> Gromacs topology files (.top) for all systems
• umb_samp --> Umbrella sampling trajectories for 23 windows for each system

Long-Timescale CG Equilibration Simulations

All long-timescale CG equilibration simulations (following the methodology from the 'Long-Timescale Pore Equilibration' section in the main text) are provided in the 'long_equil' folder. The following files are provided for each system:

• 'groups.dat' file that stores the atom indices of the membrane (Mem), lipid head phosphates (POxygens), and water W beads (Waters) that are biased as part of the implementation of ξ.
• Nucleation_CV.cpp for implementing ξ in PLUMED
• Files for each simulation replica (T1, T2, T3) for the selected ξ_long for each replica.
  • Includes a .dat file for each replica, which is used to implement a moving restraint that increases ξ over either 5.65 µs (0.8 ≤ ξ_long ≤ 0.875, Figure S7) or 6.2 µs (0.7 ≤ ξ_long ≤ 0.775, Figure 2a) to prepare MEL pore-lining configurations

Backmapping CG to AA Configurations

As described in the manuscript, MEL pore-lined configurations are first achieved with a long timescale pore equilibration step (either 5.65 or 6.2 µs, trajectory results are stored in the long_equil folder), followed by 500 ns of sampling for each window across 23 windows ranging from a flat membrane state (ξ = 0.2) to nucleated pore state (ξ = 1.0). Final CG configurations from each window are then backmapped to the CHARMM36 force field using .map files provided in the backmapping/Mapping directory.

Backmapping from the MARTINI to CHARMM36 force field can be achieved with the provided CG_to_AA.sh script, which has the following user-provided inputs at the top of the script:

• memb --> Lipid membrane type (either POPC or RBC)
• chol --> % cholesterol in membrane (0CHOL or 50CHOL)
• pep --> Whether the membrane contains melittin (8 MEL) or is a bare membrane (Pure_Memb)
• trial --> Which of the 3 independent replicas (Trial_1, Trial_2, Trial_3)
• kink --> If the MARTINI model used for MEL is fully α-helical (no) or contains a kink at T11 (yes)

Running AA Umbrella Sampling Simulations

The run_umbrella_sample.sh script provided can be utilized to set up and run atomistic (AA) umbrella sampling simulations from backmapped CG to AA initial configurations across the 23 windows of increasing pore size.

Trajectory Analysis

The following shell (.sh) and Python (.py) scripts are also provided to aid with the reproducibility of the following figures in the main text:

• Figure 2b --> Calculating number of pore-lining MEL during long-timescale CG equilibration simulations (Figures 2a, S7)
  • First, run mel_pore_line.sh
    • Centers membrane in z direction, calculates MEL density in z direction, and then integrates density profile from z = -0.5 nm to z = 0.5 nm.
  • Second, run mel_pore_line.py
    • Converts MEL pore density to the number of pore-lining MEL by normalizing the density by the density attributed to 1 pore-lining MEL
    • Visualizes Number of Pore-Lining MEL vs. time (Figure S13), which is used to calculate the 'Pore-Lining MEL' metric in Figure 2b.
• Figure 3 --> Calculating POPS flip-flop and lateral enrichment metrics from long-timescale CG equilibration simulations
  • Figure 3a --> run lipid_flip_flop.py
    • Extracts .gro file every 50 ns from simulation trajectory
    • Approximates the membrane center in the z direction by calculating the average z coordinate of all lipid C4A and C4B beads
    • Counts the number of PO4 beads in the upper and lower leaflets for each lipid type (POPC, POPE, POPS, POSM)
    • Visualizes POPS in the upper leaflet vs. time as POPS_flip_flop.svg
  • Figure 3b --> run lipid_enrich.py
    • Extracts .gro file every 50 ns from simulation trajectory
    • Counts the number of PO4 beads within a 2 nm radius of the pore center (x = 5 nm, y = 5 nm) for each lipid type (POPC, POPE, POPS, POSM)
    • Visualizes POPS within 2 nm of pore center vs. time as POPS_enrich.svg
• Figures 4a and 5 --> Average PMF profiles vs ξ across 3 replicas
  • All PMF free energy files are provided as 'freefile.txt' in the umb_samp folder
    • Free energies are resolved with Grossfield WHAM --> http://membrane.urmc.rochester.edu/?page_id=126
• Figure 6 --> Average area per lipid (APL) and bilayer thickness for atomistic umbrella sampling trajectories
  • Figure 6a --> run APL.py
    • Determines the number of lipids + CHOL in the upper leaflet
    • Calculates average membrane area during simulation with gmx energy and normalizes by the number of upper leaflet lipids + CHOL to get APL
  • Figure 6b --> run Thickness.py
    • Creates a density profile of all membrane lipid head phosphorus atoms in the z direction
    • Calculates the z value of the phosphorus peak in both the upper and lower leaflet and subtracts these quantities to get the membrane thickness
• Figure 7a --> Average MEL tilt angle
  • Run pep_helix_angle.py
    • Creates an index file (tilt.ndx) that stores the N-terminal nitrogen and C-terminal carbon for each MEL
    • Calculates the tilt angle for each MEL vs. time
    • Visualizes average tilt angle for each MEL as 'Tilt_Angles.svg'
• Figure 8 --> Calculate Lipid Fraction Near Pore
  • Run lipid_enrich_AA.py
    • 'radius' variable (in nm) in file set to target an average of 30 lipids considered for the fully nucleated pore state (ξ=1) for each RBC system considered:
      • radius = 1.6 for chol = 0CHOL and pep = Pure_Memb
      • radius = 2.0 for chol = 0CHOL and pep = 8MEL
      • radius = 2.0 for chol = 50CHOL and pep = Pure_Memb
      • radius = 2.2 for chol = 50CHOL and pep = 8MEL
    • Script calculates the number of lipid head phosphorus within 'radius' nm of the pore center (x = 5 nm, y = 5 nm) for each lipid and normalizes by the total number of lipid head phosphorus within 'radius'
    • Number of lipids within pore center 'radius' vs time is stored as 'inrange_P.csv'