Tutorials designed by Martin Calvelo (martin.calvelo@gmail.com) and Manuel N. Melo (m.n.melo@itqb.unl.pt). Check out their repository.

These tutorials were designed to run in Collab. Alternatively you can follow along in your own system.

Coarse-grain MD tutorial

Lipidomics study: Mag2 in a POPC and POPC-POPG 3:1 membrane

Tutorial notes

In this tutorial you’ll find two directories. Directory worked has all the intermediate steps and results, whereas directory minimal has only the essential files to get you going, and you’ll have to do most of the work yourself. Feel free to follow the tutorial using either approach, or even a mix of the two.
This tutorial should be run in a T4 GPU Colab runtime (you can check and change this in the upper right ‘Connect’ menu). Click on “Change runtime type” and select Python 3 (default) and T4 GPU.
Alternatively, you can download it and run it in your own system. You will need to have gromacs installed, as well as the Python libraries vermouth, mdanalysis, and nglview.
Much of the tutorial involves running shell commands. In Colab, this is achieved by prepending an exclamation mark to the command (!some shell_command). You can see this exemplified below in the installation. Note that if you need to change directories, the exclamation mark method won’t work; you need to use a percent sign, and run the command in its own cell (%cd target_dir).
To run any of the following cells of code, you can type Shift + Enter to excute the code in a cell.

Step 0 - Downloading and installing of the tutorial software and materials.

The text output of these steps was minimized so as not to overcrowd the notebook. Feel free to take out the -q flags or add the -v flag to tar if you want full output or for debugging purposes.

# We install `gromacs` and `dssp` binaries, pre-built for the colab runtime. Moreover, we install some Python packages, and import the ones we'll need:
# 1. The `vermouth` (provides the `martinize2` script) and `insane` packages from the Martini team
# 2. The `MDAnalysis` package
# 3. The `nglview` molecular viewer (requires enabling custom notebook widgets, which we also do below)

import os
if os.getenv("COLAB_RELEASE_TAG"):
  !wget https://www.itqb.unl.pt/labs/multiscale-modeling/downloads/gromacs-2023.3-gcolab_build.tar.gz -q -O - | tar --absolute-names -xzf -
  !pip install -q vermouth mdanalysis
  !pip install -q ipywidgets==7.7.2 nglview
  # Needed for nglview to work in colab
  from google.colab import output
  output.enable_custom_widget_manager()

import nglview as nv
import MDAnalysis as mda

Finally, we download and unpack the tutorial’s material.

!wget https://zenodo.org/records/13890148/files/2-Tutorial_Lipidomics.zip
!unzip 2-Tutorial_Lipidomics.zip

We will start in subdirectory minimal, but feel free to change to worked.

initial_path=!pwd
path=initial_path[0]
print(path)

%cd {path}/minimal

Step 1 - Introduction

Antimicrobial peptides (AMPs) play a vital role in the innate immune system of all living organisms.\(^{1,2}\) These peptides can identify and disrupt lipid patterns in bacterial membranes without harming healthy cells.\(^{3}\) Notably, AMPs often possess shared properties: they typically consist of 10–40 residues, adopt an α-helix secondary structure when they interact with membranes, maintain a positive net charge, and contain a significant number of non-polar residues. These features, along with their pronounced transversal hydrophobic moment in their folded state, facilitate interactions with lipid bilayers. Thus, the use of these molecules as therapeutic agents is an obvious lesson from Nature that can be exploited to design new antibiotics, anticancer or antiaging drugs.\(^{4}\)

In this tutorial, we will investigate the potential antimicrobial properties of Magainin-II (Mag2), one of the most well-known natural AMPs.\(^{5}\) We will simulate it in the presence of two different membrane models: one representing a healthy membrane (POPC), and the other a bacterial bilayer (POPC-POPG 3:1). We’ll place several AMPs in the aqueous phase of our systems and run simulations, analysing potential spontaneous peptide-bilayer interactions.

Lancet Infect Dis. 9, e216; doi: 10.1016/S1473-3099(20)30327-3 (2020)
Front. Microbiol. 11; doi: 10.3389/fmicb.2020.582779 (2020)
Colloids Surf. B Biointerfaces 196, 111349; doi: 10.1016/j.colsurfb.2020.111349 (2020)
Front. Immunol. 15, doi: 10.3389/fimmu.2024.1320779 (2024)
BBA. 1327, 119; doi: 10.1016/S0005-2736(97)00051-5 (1997)

Step 2 - System preparation

In your directories you’ll find Mag2.pdb file, containing a α-helix structure of the Magainin II peptide.

Before coarse graining anything we visualize our starting protein structure using nglview. Other viewers, such as pymol or VMD are also great, if you run this workflow outside Colab.

view = nv.NGLWidget()
at = view.add_component(nv.FileStructure('Mag2.pdb'))
view

Coarse-graining with `martinize2`

For CG MD simulation of a protein one needs the CG coordinates and the CG topology. Both can be obtained with the martinize2 tool (despite its name, it works with Martini 3).

You’ll need to: * specify Mag2.pdb as the input structure; * specify an output CG structure (let’s call it Mag2_cg.pdb); * specify an output .top topology master file (let’s call it topol.top); * specify the target forcefield. Martini 3 has code martini3001; * ask for an elastic network to be set, restraining secondary/tertiary structure; * ask for fixing some issues with the excess of flexible side-chain behaviour; * we don’t need DSSP in this case, since we will assume directly that the peptide will adopt an alpha-helix structure. You can provide directly the secondary structure with the flag -ss HHHHHHHHHHHHHHHHHHHHHHH;

#Try yourself!
#Uncomment to run

#!martinize2 -f ...

# @title Solution

!martinize2 -f Mag2.pdb \
 -x Mag2_cg.pdb \
 -o topol.top \
 -ff martini3001 \
 -elastic \
 -scfix \
 -ss 'HHHHHHHHHHHHHHHHHHHHHHH'

Have a look at the produced files with and editor or a molecular structure visualizer. Observe the atomic structure and the CG beads overlaid (in white and yellow spheres).

view.remove_component(at)

aa = nv.FileStructure('Mag2.pdb')
cg = nv.FileStructure('Mag2_cg.pdb')

view.add_component(aa, default=False)
view.add_component(cg, default=False)

view.clear_representations(component=0)
view.clear_representations(component=1)

view.add_licorice(selection='Protein', component=0)
view.add_ball_and_stick(aspectRatio=7.5, opacity=0.45, component=1)

view.center()
view

POPC membrane

Let’s continue with the tutorial. The steps are the same for both membranes. Let’s start with the POPC one:

%cd {path}/minimal/POPC

Insert Mag2 peptides

The first step will be to insert some Mag2 peptides into our system (let’s add 10, but feel free to choose a different number if you prefer). For this, we’ll use the gmx insert-molecules command from GROMACS. Execute it with the flag -h for having more info.

!gmx insert-molecules -h

We have a POPC membrane in water box already equilibrated in a file called POPC.gro. Let’s place the peptides in the aqueous phase. In this way, we’ll see if they spontaneously interact with the membrane. We can do this by instructing gmx insert-molecules to remove water molecules (W) to make room for our peptide (use the option -replace). Take into account that the CG coordinates of Mag2 are in the previous directory, so we need to write ../ before the file name (../Mag2_cg.pdb). Let’s call the output system.gro. Try yourself!

#Try yourself!
#Uncomment to run

#!gmx insert-molecules -f ... -ci .../... -nmol ... -o system.gro -replace W -try 100

# @title Solution

!gmx insert-molecules -f POPC.gro \
                      -ci ../Mag2_cg.pdb \
                      -nmol 10 \
                      -o system.gro \
                      -replace W \
                      -try 100

Open the output system.gro with a molecular visualizator. Here, you can see the peptides in the solvent, with the POPC lipids in grey. Note that, for clarity, water and ions were not displayed.

view_membrane_prot = nv.NGLWidget()
view_membrane_prot.add_component(nv.FileStructure('system.gro'))
view_membrane_prot.add_representation('tube', selection='protein', radius=1)
view_membrane_prot.add_representation('surface', selection='POPC', color='grey')
view_membrane_prot

Update topology (1)

The file initial_top.top is a good start, but we need to complete it. Specifically, we need to add the number of W, CL, NA, and peptides in the system. The first step is to create a copy of the file initial_top.top (let’s call this new file top.top).

!cp initial_top.top top.top

We can automatically update the number of W, CL, NA, and peptides molecules by using the python script update_topology.py. It will count the number of those molecules in a input file (system.gro) and add it directly to a topology file (top.top). Run the script like this:

!python3 ../update_topology.py system.gro top.top

Add also a line with the number of peptides inserted:

%%bash
echo "molecule_0  10" >> top.top

Neutralize the system

MD simulations must be executed with a total net charge of 0. Since Mag2 has a charge of +2, adding the peptides has resulted in an imbalance in the system’s charge. We need to fix this. We can do this using the GROMACS tool gmx genion. To start, we need to create a .tpr file. Generate it by running the following command:

!gmx grompp -f ../em.mdp -c system.gro -p top.top -o ions.tpr -maxwarn 1

Now let’s run the gmx genion command. We’ll achieve a neutral charge with the -neutral flag and will need to specify which molecules it will replace to add extra ions to balance the charge (water, in this case).

%%bash
echo "W" | gmx genion -s ions.tpr -neutral -o system_neutralized.gro

This command will generate a file with the coordinates of the neutralized system called system_neutralized.gro.

Update topology (2)

Since we’ve changed the number of water molecules and ions when neutralizing the system, we need to update the topology again. Let’s start again from the initial_top.top file:

!cp initial_top.top top.top

And let’s run the update_topology.py script again, but this time using system_neutralized.gro as the input file.

!python3 ../update_topology.py system_neutralized.gro top.top

Remember to add also a line with the number of peptides:

%%bash
echo "molecule_0  10" >> top.top

Take a look to the top.top file:

!cat top.top

Create and index file

Let’s create an index file defining the beads that are part of the protein and membrane and those that are part of the solvent. Run first the gmx make_ndx command of GROMACS using as input the system_neutralized.gro file:

%%bash
echo 'q' | gmx make_ndx -f system_neutralized.gro

Using the create_index.py script, we’ll create 3 groups: * Bilayer: beads of the membrane. * Bilayer_Protein: beads of the membrane and the AMPs. We’ll use this for analysis. * Solvent: Water and ions.

Run the script:

!python3 ../create_index.py system_neutralized.gro

POPC-POPG 3:1 membrane

As we said, in this tutorial we are going to simulate two membranes. Now that you have seen how to prepare the system for one of them, try preparing the other one yourself!

%cd {path}/minimal/POPC-POPG_31/

Insert Mag2 peptides

Take into account that now our membrane file is called POPC-POPG_31.gro.

# Complete it

!gmx insert-molecules

# @title Solution

!gmx insert-molecules -f POPC-POPG_31.gro \
                      -ci ../Mag2_cg.pdb \
                      -nmol 10 \
                      -o system.gro \
                      -replace W \
                      -try 100

Update topology (1)

!cp initial_top.top top.top

!python3 ../update_topology.py system.gro top.top

%%bash
echo "molecule_0  10" >> top.top

Neutralize the system

# Complete it

!gmx grompp -f ../em.mdp

# @title Solution

!gmx grompp -f ../em.mdp -c system.gro -p top.top -o ions.tpr -maxwarn 1

%%bash
# echo ... | !gmx genion

# @title Solution

!echo "W" | gmx genion -s ions.tpr -neutral -o system_neutralized.gro

Update topology (2)

!cp initial_top.top top.top

!python3 ../update_topology.py system_neutralized.gro top.top

%%bash
echo "molecule_0  10" >> top.top

Create index file

%%bash
echo 'q' | gmx make_ndx -f system_neutralized.gro

!python3 ../create_index.py system_neutralized.gro

Step 3 - Simulation

You can now run the simulations. In CG it is enough to perform a short energy minimization and a single pressure/temperature equilibration step before production. In GROMACS, we need to use the command gmx grompp to create a processed binary file (with the extension .tpr) for running minimization/MD calculations. The instructions for the calculation are in .mdp files. Again, let’s start with the POPC membrane:

POPC membrane

%cd {path}/minimal/POPC/

Energy minimization

Let’s create the tpr file (em.tpr) to perform the energy minimization. You will need a mdp file, available in the previous directory (../em.mdp). Open it with a text editor to check it.

Use the grompp tool of GROMACS. You can run gmx grompp with the flag -h for listing all the possible options.

#Try yourself!
#Uncomment to run

#!gmx grompp ...

# @title Solution

!gmx grompp -f ../em.mdp -c system_neutralized.gro -p top.top -n index.ndx -o em.tpr -maxwarn 2

Run now the minimization with:

!gmx mdrun -deffnm em

Equilibration

We will now equilibrate the system using the file ../eq.mdp. You can donwload it and check it with a text editor. Note that we wiil use a semiisotropic pressure coupling (isotropic in the x and y direction, but different in the z direction), something typical when we work with membranes.

Let’s create the tpr file (eq.tpr), and remember to include the index.ndx file:

#Try yourself!
#Uncomment to run

#!gmx grompp -f ../eq.mdp ...

# @title Solution

!gmx grompp -f ../eq.mdp -c em.gro -p top.top -n index.ndx -o eq.tpr -maxwarn 2

Run now the equilibration:

!gmx mdrun -v -deffnm eq

Production

Finally, we can create the .tpr and run the production using the md.mdp file.

#Try yourself!
#Uncomment to run

#!gmx grompp ...

# @title Solution

!gmx grompp -f ../md.mdp -p top.top -c eq.gro -n index.ndx -o md.tpr

Since this step is quite computationally expensive, we’ll take the trajectory file (md.xtc) from the worked directory (also, some additional files to avoid incompatibilities).

!cp {path}/worked/POPC/md.xtc .
!cp {path}/worked/POPC/index.ndx .
!cp {path}/worked/POPC/md.gro .
!cp {path}/worked/POPC/md.tpr .

POPC-POPG 3:1 membrane

Now, try to do it yourself with the other membrane!

%cd {path}/minimal/POPC-POPG_31/

Energy minimization

# Create the em.tpr file

!gmx

# @title Solution

!gmx grompp -f ../em.mdp -c system_neutralized.gro -p top.top -n index.ndx -o em.tpr -maxwarn 2

# Run the minimization

!gmx

# @title Solution

!gmx mdrun -deffnm em

Equilibration

#Create the eq.tpr file

!gmx

# @title Solution

!gmx grompp -f ../eq.mdp -c em.gro -p top.top -n index.ndx -o eq.tpr -maxwarn 2

#Run the equilibration

!gmx

# @title Solution

!gmx mdrun -v -deffnm eq

Production

#Create the md.tpr file

!gmx

# @title Solution

!gmx grompp -f ../md.mdp -p top.top -c eq.gro -n index.ndx -o md.tpr

Don’t run the production, let’s copy the files from the worked directory:

!cp {path}/worked/POPC-POPG_31/md.xtc .
!cp {path}/worked/POPC-POPG_31/index.ndx .
!cp {path}/worked/POPC-POPG_31/md.gro .
!cp {path}/worked/POPC-POPG_31/md.tpr .

Step 4 - Analysis

We’re going to make some adjustments to the trajectory and perform some analyses. After that, we will visualize the results.

POPC membrane

Let’s start this section analysing the POPC membrane trajectory:

%cd {path}/minimal/POPC/

Trajectory processing

Let’s remove the water and ions and make all molecules whole:

%%bash
echo 'Bilayer_Protein' | gmx trjconv -f md.xtc -s md.tpr -pbc whole -o md_pbc.xtc -n index.ndx

Let’s do the same, but using only the last frame of the simulation:

%%bash
echo 'Bilayer_Protein' | gmx trjconv -f md.gro -s md.tpr -pbc whole -o md_pbc.gro -n index.ndx

Let’s convert the production tpr file into another with fewer atoms compatibl with our centered trajectory (since we now ignore the solvent).

%%bash
echo 'Bilayer_Protein' | gmx convert-tpr -s md.tpr -n index.ndx -o md_pbc.tpr

Lipid contacts counting

To count contacts we first get a list of closest distances between any lipid and any AMP over time, using gmx mindist. Display first the help options with the -h flag.

!gmx mindist -h

Try now to calculate the contacts using gmx mindist. You will need to choose the options Protein and non-Protein (since there are not water or ions). Moreover, use the flag -group (check in the info what it does). Name the output contacts_amp-bilayer.xvg using the flag -o.

#Try yourself!
#Uncomment to run

#%%bash
#echo 'Protein' 'non-Protein' | gmx mindist ...

# @title Solution

!echo 'Protein' 'non-Protein' | gmx mindist -s md_pbc.tpr -f md_pbc.xtc -on contacts_amp-bilayer.xvg -group

POPC-POPG 3:1 membrane

Now, let’s do it with the other membrane.

%cd {path}/minimal/POPC-POPG_31/

Trajectory processing

Let’s remove the waters and ions and make molecules whole. We will also do the same for just one frame and create a new .tpr file compatible with the treated trajectory.

%%bash
echo 'Bilayer_Protein' | gmx trjconv -f md.xtc -s md.tpr -pbc whole -o md_pbc.xtc -n index.ndx
echo 'Bilayer_Protein' | gmx trjconv -f md.gro -s md.tpr -pbc whole -o md_pbc.gro -n index.ndx
echo 'Bilayer_Protein' | gmx convert-tpr -s md.tpr -n index.ndx -o md_pbc.tpr

Lipid conctacts counting

And now let’s calculate the number of contacts. Do you think you’d be able to do it yourself? Give it a try!

#Try yourself!
#Uncomment to run

#%%bash
#echo 'Protein' 'non-Protein' | gmx mindist ...  -o contacts_amp-bilayer.xvg

# @title Solution

!echo 'Protein' 'non-Protein' | gmx mindist -s md_pbc.tpr -f md_pbc.xtc -on contacts_amp-bilayer.xvg -group

Visualize the last frame

Let’s start by visualizing the last structure of the simulations. Go to the main directory.

%cd {path}/minimal/

First, the one with POPC:

view_membrane_PC = nv.NGLWidget()
view_membrane_PC.add_component(nv.FileStructure('POPC/md_pbc.gro'))
view_membrane_PC.add_representation('tube', selection='protein', radius=1)
view_membrane_PC.add_representation('surface', selection='POPC', color='grey')
view_membrane_PC

Compare it with the POPC-POPG 3:1 membrane.

view_membrane_PCPG = nv.NGLWidget()
view_membrane_PCPG.add_component(nv.FileStructure('POPC-POPG_31/md_pbc.gro'))
view_membrane_PCPG.add_representation('tube', selection='protein', radius=1)
view_membrane_PCPG.add_representation('surface', selection='POPC', color='grey')
view_membrane_PCPG.add_representation('surface', selection='POPG', color='green')
view_membrane_PCPG

Apparently, Mag2 interacts more with the bacterial membrane. However, we are only observing one frame of the trajectory, which might not be representative. Let’s visualize the results of the contacts during the entire trajectory to gain more information.

Print the contacts

After calculating the contacts between the AMPs and the lipids in former sections, we will plot them for visualization. In the main directory we have a Python script called contacts_vs_time.py, which will plot the contacts betwen AMPs and lipids per membrane. Run it.

%matplotlib inline

!python3 contacts_vs_time.py

This script will generate an output in PDF format (contacts_amp-bilayer.pdf) that you can download. To visualize it in this Colab, run the following cell:

import numpy as np
import matplotlib.pyplot as plt

POPC = np.loadtxt('POPC/contacts_amp-bilayer.xvg', comments=('@', '#'))
POPCPG = np.loadtxt('POPC-POPG_31/contacts_amp-bilayer.xvg', comments=('@', '#'))
contacts_POPC = POPC[:,1:]
contacts_POPCPG = POPCPG[:,1:]
time = len(contacts_POPC)
plt.plot(np.arange(1, len(contacts_POPC) + 1), contacts_POPC, label='POPC', color='red')
plt.plot(np.arange(1, len(contacts_POPCPG) + 1), contacts_POPCPG, label='POPC-POPG 3:1', color='blue')
plt.title('Contacs AMP - Bilayer')
plt.legend()
plt.xlabel('Time (ns)')
plt.ylabel('# contacts')
plt.show()

As you can see, the number of contacts is higher in the bacterial membrane model. Do you think Mag2 has antimicrobial potential?

Conclusion

Congratulations! You’ve reached the end of this tutorial! Visit https://cgmartini.nl and continue exploring Martini and Coarse-Grained MD!

Coarse-grain MD tutorial

Lipidomics study: Mag2 in a POPC and POPC-POPG 3:1 membrane

Tutorial notes

Step 0 - Downloading and installing of the tutorial software and materials.

Step 1 - Introduction

Step 2 - System preparation

Coarse-graining with martinize2

POPC membrane

Insert Mag2 peptides

Update topology (1)

Neutralize the system

Update topology (2)

Create and index file

POPC-POPG 3:1 membrane

Insert Mag2 peptides

Update topology (1)

Neutralize the system

Update topology (2)

Create index file

Step 3 - Simulation

POPC membrane

Energy minimization

Equilibration

Production

POPC-POPG 3:1 membrane

Energy minimization

Equilibration

Production

Step 4 - Analysis

POPC membrane

Trajectory processing

Lipid contacts counting

POPC-POPG 3:1 membrane

Trajectory processing

Lipid conctacts counting

Visualize the last frame

Print the contacts

Conclusion

Coarse-graining with `martinize2`