Start the workshop on Rosi

Connect to Rosi in an SSH terminal and initialize your private workshop copy:

. /home/tut04/pore-workshop/start

The launcher verifies release 0.9.0, installs it as ~/pore-workshop-0.9.0, enters that directory, and runs the cluster preflight. Leave this terminal open for the simulation commands below.

Introduction and Background: transmembrane pore formation

A membrane pore is a topological transition: an intact bilayer becomes connected by a water-filled polar defect. The chain coordinate ξ_ch turns that connectivity into a continuous value that can be restrained and used to calculate a potential of mean force (PMF).

The formation of membrane pores plays key roles in numerous processes of biological or medical relevance: Pores are formed during apoptosis, by antimicrobial peptides, or during membrane fusion (which in turn is involved in signaling, viral infection, and fertilization). During electroporation, pores are formed with an external electric pulse with the aim delivering cargos into cells. In MD simulations, pores do not form spontaneously because pore formation is typically associated with a high free energy cost. Therefore, to observe pores formation in MD simulations and to study the roles of lipids, membrane-active pepties, or electric fieds on pore formation, pore formation must be simulated by pulling along reaction coordinate (or collective variable): the chain coordinate.

A fromation of a transmembrane pore is a topological transition: an intact bilayer becomes protruded by a water-filled continuous polar defect. The chain coordinate ξ_ch turns the degree of connectivity of the transmembrane defect into a continuous value that can be restrained and used to calculate the potential of mean force (PMF) of pore formation.

Warmup: watch an equilibrium membrane

Before opening a pore, look at a normal Martini 3 bilayer moving at equilibrium. VMD runs on your own computer; the simulations remain on the cluster.

1. Download and install the appropriate VMD build from the official VMD alpha download.
2. Add the directory containing the vmd executable to your local shell’s PATH. For example, replace the example directory below with the actual installation directory:

export PATH="/opt/vmd/bin:$PATH"
command -v vmd

To keep this setting, add the same export PATH=... line to ~/.bashrc for Bash or ~/.zshrc for Zsh, then open a new terminal. The installation directory depends on your operating system and chosen VMD package.

Open another terminal on your own computer and copy the prepared example from Rosi. Set ROSI_LOGIN to the hostname or SSH alias you normally use for Rosi:

ROSI_LOGIN=rosi
LOCAL_EXAMPLE="$HOME/workshop-visualization/equilibrium-membrane"

mkdir -p "$LOCAL_EXAMPLE"
scp \
    "${ROSI_LOGIN}:pore-workshop-0.9.0/reference/equilibrium-vmd/equilibrium-conect.pdb" \
    "${ROSI_LOGIN}:pore-workshop-0.9.0/reference/equilibrium-vmd/equilibrium-whole.xtc" \
    "$LOCAL_EXAMPLE/"

vmd \
    "$LOCAL_EXAMPLE/equilibrium-conect.pdb" \
    "$LOCAL_EXAMPLE/equilibrium-whole.xtc"

If you already mount your Rosi home directory with SSHFS or an SFTP filesystem client, you may open the same two files directly from pore-workshop-0.9.0/reference/equilibrium-vmd/. Run all simulation commands in the Rosi SSH terminal, not on your local computer.

VMD now shows and example of a previously carried out equilibriums simulation of a membrane of POPC, without any pore-opening restraint. Play or scrub the trajectory, rotate the bilayer, and identify water, headgroups, and the hydrophobic core.

The example includes equilibrium.gro, a one-frame GROMACS structure, and equilibrium-whole.xtc, the trajectory. VMD loads equilibrium-conect.pdb first because a GRO or XTC does not carry the Martini bond topology.

./scripts/14_prepare_vmd.sh topol.tpr traj.xtc index.ndx visualization
vmd visualization/structure-conect.pdb visualization/trajectory-whole.xtc

From polar connectivity to pore size

The chain coordinate ξ_ch measures the degree of connectivity of the transmembrane defect. To quantify “connectivity”, it uses the fraction of polar-occupied slices in a membrane-spanning cylinder (see left image). In case you also want to simulate pore expansion, an extension of the chain coordinate may be used to switch from pore nucleation (left image) to pore expansion (measured as as the pore radius), together implemented as a joint coordinate ξ_p for pore nucleation and expansion. Adapted from Hub, JCTC 2021.

In the published atomistic definition, water oxygen and lipid phosphate oxygen atoms contribute when they lie inside the cylinder. An intact membrane leaves central slices empty; water entry occupies more slices; ξ_ch approaches unity when the defect spans the entire membrane. The chain coordinate then saturates and cannot describe further pore expansion. The joint coordinate ξ_p solves that problem by switching smoothly to the expansion term.

Explore the full published pathway

Use the slider to follow atomistic structures from an intact membrane through pore nucleation and expansion. The values at the bar are ξ_p; in this workshop, however, we only use the chain coordinate ξ_ch to focus on the pore nucleation phase.

What you will do

The workshop has two successive applications of the chain coordinate.

First, you will study pore formation in a Martini 3 membrane. You will choose the cylinder resolution and opening rate, generate 27 short umbrella windows, and reconstruct a pore-formation PMF.

After completing that workflow, you will apply the same ideas to stalk formation between two facing Martini 2 membranes. You will calibrate the coordinate for a double membrane, follow a 50 ns stalk-opening trajectory, generate 19 umbrella windows, and reconstruct a stalk-formation PMF.

The workflows use separate run directories and scripts, so their simulation files remain independent. Conceptually, the stalk exercise builds on what you learned during the pore exercise.

Everything else uses tested defaults. Optional explanations are available in expandable Advanced panels.

Run the showcase →

Run the showcase

Keep this page open beside a terminal. Commands beginning with ./scripts/ must be run from the copied workshop folder.

1. Choose a membrane and create your run

SYSTEM is the exact membrane ID. RUN is the new private directory that will hold your calculation. For example:

export SYSTEM=popc-small
export RUN="$(realpath -m "$PWD/../${USER}-${SYSTEM}")"
./scripts/01_create_run.sh "$SYSTEM" "$RUN"
./scripts/02_inspect_system.sh "$RUN"

Replace only popc-small. Your run is independent and is never combined with another participant’s files.

GROMACS file types

structure + topology + groups + settings  →  grompp  →  .tpr
                                                       ↓ mdrun
                                                .xtc + .xvg

printf 'SYSTEM=%s\nRUN=%s\n' "$SYSTEM" "$RUN"
ls -lh "$RUN/base"
tail -1 "$RUN/base/equilibrated.gro"
less "$RUN/base/topol.top"
less "$RUN/base/index.ndx"

2. Choose the number of slices (or length of the cylinder)

The membrane-spanning cylinder is divided into N slices with thickness of 2Å (1Å for all-atom simulations). Hence, the definition of ξ_ch requires choosing N. Choose N (and thereby the length) such that, for a planar intact membrane (with no pore), ξ_ch takes a value of ~0.2 meaning 20% of the cylinder slices are occupied by polar atoms, namely the slices in the two head group regions. This ensures that the cylinder spans the entire membrane including some parts of the polar headgroup reagion. If the cylinder is too short such that no slice contains a polar atom, you get an error message. In practice, you run a series of so-called reruns from a given equilibrium simulation trajectory, compute ξ_ch with various values of N form the frames. Then choose a good value. We did this for you in advance:

Inspect the results from these reruns, which we did in advance for you:

./scripts/03_scan_slices.sh "$RUN" --reference
column -t -s $'\t' "$RUN/calibration/scan-status.tsv" | less

Look for three things:

1. reject values marked FAIL - cylinder is too short, N too small;
2. observe how the computed ξ_ch changes with the choice of N;
3. choose a passing N with an mean close to 0.2–0.25, making sure that the cylinder spans the membrane.

Make a provisional choice and replace CHOOSE_FROM_SLICE_SCAN:

${EDITOR:-nano} "$RUN/input/coordinate.env"
./scripts/03_scan_slices.sh "$RUN" --check-choice

Try another plausible value if you want to compare. Your first choice does not have to match the tested one. Decide from the evidence before opening the solution, then reconcile your value before selecting a fallback opening.

./scripts/03_scan_slices.sh "$RUN" --solution
cat "$RUN/calibration/selected-slices.tsv"

3. Choose the opening rate

The pore-opening pulling simulation moves the harmonic-restraint centre from ξ_ch=0.1 to 1.0. Inspect the timestep and number of steps in the MDP file:

grep -E '^[[:space:]]*(dt|nsteps|pull-coord1-init)' \
    "$RUN/input/opening.mdp"

First calculate the duration dt × nsteps in ps. Then calculate:

pull rate = (desired final value − initial value) / duration in ps

Replace CHOOSE_OPENING_RATE_XI_PER_PS in the MDP file with your result and test it. An incorrect value is safe here: the validator reports the endpoint and refuses to prepare a simulation until it is close to 1.0.

${EDITOR:-nano} "$RUN/input/opening.mdp"
./scripts/04_validate_configuration.sh "$RUN"

Revise the rate if validation fails. Reveal the solution only after making a first calculation.

(1.0 − 0.1) / 20,000 ps = 4.5 × 10⁻⁵ ps⁻¹

Once both choices validate, prepare and submit. Here the prepare and submit scripts essentially run the Gromacs commands gmx grompp and gmx mdrun, but help you with getting files into the right folders:

./scripts/05_prepare_opening.sh "$RUN"
./scripts/06_submit_opening.sh "$RUN"
./scripts/07_check_opening.sh "$RUN"

If the status says your opening is complete, select it:

./scripts/07_check_opening.sh "$RUN" --select-own

If it is still queued or running, select the validated opening simulation for the same membrane instead:

./scripts/07_check_opening.sh "$RUN" --use-fallback

grep -E 'nsteps|pull-nstfout|pull-coord1-(geometry|init|rate|k)' \
    "$RUN/opening/opening.mdp"

GROMACS writes the coordinate to opening_pullx.xvg and the applied force to opening_pullf.xvg. Take a look, e.g. with xmgrace:

xmgrace opening_pullx.xvg

and inspect the simulation with VMD. Do you see the pore?

4. Generate and submit the umbrella sampling simulations

Umbrella windows restrain the system at neighboring ξ_ch values. From each window, a histogram along ξ_ch is collected. From all histograms the PMF is computed with the so-called weighted histogram analysis method (WHAM).

./scripts/08_generate_window_grid.sh "$RUN"
./scripts/10_submit_windows.sh "$RUN"
./scripts/11_status.sh "$RUN"

This creates 27 independent umbrella windows, each simulated for 2 ns. Slurm starts each one when a GPU is available.

./scripts/09_prepare_representatives.sh "$RUN" --preview
column -t -s $'\t' "$RUN/input/window-regions.tsv"
less "$RUN/state/window-grid.tsv"
grep -H -E 'pull-coord1-(geometry|init|k)' \
    "$RUN"/state/representative-preview/*/umbrella.mdp

See the result →

See the result

Inspect a prepared result

If your simulations are still running, first inspect the prepared result for the same membrane:

cat "reference/results/$SYSTEM/analysis-summary.txt"

Open reference/results/$SYSTEM/report.html in a browser. A useful result has:

1. all 27 windows completed;
2. connected overlap between every neighboring pair;
3. a PMF shown with uncertainty only over the sampled ξ_ch range.

Analyze your result

Restore RUN, check completion, and build your report:

export RUN=/home/alice/alice-popc-small
./scripts/11_status.sh "$RUN"
./scripts/12_analyze.sh "$RUN"

Then open $RUN/analysis/report.html. You now have the complete workflow from an equilibrated membrane to a diagnosed pore-nucleation PMF.

Critically compare PMF of pore formation obtaind with MARTINI with PMF from all-atom simulations

Figure S1 from Ting et al., PRL 2018: atomistic CHARMM36 PMFs for phosphatidylglycerol (PG) membranes.

This figure shows PMF of pore formation obtained with all-atom simulations. Check also the PDF of Ting et al. for PMFs with PC lipids. Note the local minima of the PMFs at ξ_ch =1. What do these minima imply, for instance with respect to the life time of open pores?

Do you see such minimal in the PMFs obtained with MARTINI? Comare with your colleagues who possibly simulated other lipid types.

What can you say about the free energy cost of pore formation in MARTINI compared to all-atom CHARMM36 simulations? Note that $kT \approx 2.5 \mathrm{kJ/mol}$. Note that the probability of pore formation follows:

$$ P_\mathrm{pore} \propto \exp(-\Delta G/k_B T), $$

where $\Delta G$ is the free energy of pore formation. By which factor is the probability of pore formation underestimated by Martini?

cat "$RUN/analysis/analysis-summary.txt"
column -t "$RUN/analysis/overlap-summary.tsv" | less
grep -RE 'Fatal error|Segmentation fault|CUDA error' "$RUN"
find "$RUN/windows" -name 'pullf.xvg' -size +0c | wc -l

Next: form a stalk

Now that you have used the chain coordinate to describe pore formation, apply the same principles to a different topological transition: stalk formation between two facing membranes. The simulations use separate files, but the coordinate-calibration, umbrella-sampling, and PMF concepts carry over.

Continue to Part 2: stalk formation →

Optional: continue with pore expansion

The workshop stops at a connected polar defect. The extended GROMACS Chain Coordinate example continues with the joint reaction coordinate ξ_p for pore nucleation and expansion, including its smooth switch from ξ_ch to pore-radius expansion, and provides a more extensive tutorial.

Optional advanced reference →

Optional advanced reference

This page is optional lookup material.

Essential terms

export SYSTEM=popc-small                 # replace with one listed system ID
export RUN="$(realpath -m "$PWD/../${USER}-${SYSTEM}")"
./scripts/00_preflight.sh
./scripts/01_create_run.sh "$SYSTEM" "$RUN"
./scripts/02_inspect_system.sh "$RUN"
./scripts/03_scan_slices.sh "$RUN" --reference
${EDITOR:-nano} "$RUN/input/coordinate.env"   # enter your provisional N
./scripts/03_scan_slices.sh "$RUN" --check-choice
./scripts/03_scan_slices.sh "$RUN" --solution # reveal only after choosing
${EDITOR:-nano} "$RUN/input/opening.mdp"      # calculate and test the rate
./scripts/04_validate_configuration.sh "$RUN"
./scripts/05_prepare_opening.sh "$RUN"
./scripts/06_submit_opening.sh "$RUN"
./scripts/07_check_opening.sh "$RUN"
# Select exactly one opening source:
./scripts/07_check_opening.sh "$RUN" --select-own
# or: ./scripts/07_check_opening.sh "$RUN" --use-fallback
./scripts/08_generate_window_grid.sh "$RUN"
./scripts/09_prepare_representatives.sh "$RUN" --preview
./scripts/10_submit_windows.sh "$RUN"
./scripts/11_status.sh "$RUN"
./scripts/12_analyze.sh "$RUN"

Sources and further reading →

Sources and further reading

Reaction coordinate

Hub, J. S. and Awasthi, N. (2017). “Probing a Continuous Polar Defect: A Reaction Coordinate for Pore Formation in Lipid Membranes.” Journal of Chemical Theory and Computation. DOI: 10.1021/acs.jctc.7b00106

Hub, J. S. (2021). “Joint Reaction Coordinate for Computing the Free-Energy Landscape of Pore Nucleation and Pore Expansion in Lipid Membranes.” Journal of Chemical Theory and Computation. DOI: 10.1021/acs.jctc.0c01134 · group-hosted PDF

Hub group. “Example system illustrating how to use the joint reaction coordinate for pore formation and expansion.” GROMACS Chain Coordinate tutorial

Applications and interpretation

Ting, C. L., Awasthi, N., Müller, M., and Hub, J. S. (2018). “Metastable Prepores in Tension-Free Lipid Bilayers.” Physical Review Letters 120, 128103. DOI: 10.1103/PhysRevLett.120.128103 · Supporting information

Starke, S., Allolio, C., and Hub, J. S. (2025). “How pore formation in complex biological membranes is governed by lipid composition, mechanics, and lateral sorting.” PNAS Nexus. DOI: 10.1093/pnasnexus/pgaf033

Foundations, convergence, and finite size effects

Awasthi, N. and Hub, J. S. (2016). “Simulations of Pore Formation in Lipid Membranes: Reaction Coordinates, Convergence, Hysteresis, and Finite-Size Effects.” Journal of Chemical Theory and Computation. DOI: 10.1021/acs.jctc.6b00369

The published coordinate and pore-state images are adapted from Hub, JCTC 2021. The atomistic PMF comparison is Figure S1 from the supporting information of Ting et al., PRL 2018. Software versions and system provenance are recorded in PROVENANCE.md.

MD simulations with Gromacs, visualization

University of Illinois. VMD alpha downloads and VMD command-line options. VMD accepts a PDB structure on its command line and can append a trajectory to that molecule.

GROMACS 2021 documentation for gmx trjconv and gmx editconf. The -conect option writes PDB bond records for non-standard molecules such as coarse-grained lipids; -pbc whole makes broken bonded molecules whole.

SSHFS project. SSHFS usage and installation. SSHFS exposes a remote SFTP directory as a local filesystem using sshfs user@host:directory mountpoint.

Continue to stalk formation

You have completed Part 1. Continue with the second chain-coordinate workflow.

Continue to Part 2: stalk formation →