In this session we will deal with the setup of a Molecular Dynamics simulation with a cell membrane. Membranes will not likely be the focus of the course, but we deem it useful to give you the opportunity to get your hands dirty with some lipid bilayers ;).

Membranes in a nutshell

As you may imagine, the cell membrane separates the inner part of the cell, the cytoplasm, from the extracellular space, the so-called extracellular matrix (ECM). As usual, biology is complicated, so there are many other membranes in your body and in the cell. One example is the nuclear membrane, which separates the chromosomes from the rest of the cell.

The cell membrane is made mainly by a lipid bilayer, that is two arrays of lipids stacked one on the other. With such an arrangement the hydrophilic lipidic heads are exposed to the solvent on both sides of the membrane, while the hydrophobic tails are separated from water in the inner region. Lipid molecules can diffuse in the membrane and the membrane itself can be deformed.

Membranes are formed principally by glycerophospholipids, but other type of lipids are present as well. Different cells and/or organisms can have a different concentration of those lipids in their membranes, resulting in different physical and chemical properties (e.g. more rigidity).

As usual, biology is complicated (2.0).

The cell membrane controls the inward and outward flux of substances, including ions (concept of selective permeability). This is possible thanks to the presence of proteins. Indeed, there are two main categories of proteins interacting with the membranes: the integral proteins, which are completely embedded in the membrane. Among them, transmembrane proteins cover the whole length of the lipid bilayer. Ion channels,like the potassium channels we saw in tutorial1, are an example of integral protein. They allow the flux of ions from and to the cell. Peripheral proteins play a different role, being only loosely attached to the membrane. These structures are usually involved in maintaining the shape of the bilayer or in signalling processes. Proteins can also diffuse on the membrane.

What membrane motions, if any, should we consider if we are interested in modelling a membrane protein such as a pore?

suspense…

In molecular modelling of membrane proteins we usually restrict our analysis to a limited portion of the lipid bilayer (typically from 10 to 30 nanometers) in order to save computational time. In simulations we tend to neglect the curvature of the membrane, therefore modelling it as a perfect square.

An artistic rendering of recoverin attached to a membrane patch (Crespi & Mazza, comp. biophys report 2019-2020).

The force field

Nowadays, standard force fields model the lipid-lipid interactions and the potentials between lipids and proteins in a pretty accurate way. Nonetheless, there can be differences, in general. Even without the membrane, different ff can make some regions of a folded protein more or less stable.

Different force fields can give different results

The image was taken by this nice article by Lindorff-Larsen and coworkers.

VMD Membrane builder extension

VMD has an extension that allows the user to build a custom squared lipid bilayer. If you have it (it is indeed an extension) try to build a 50x50 membrane with POPC lipids.

Launch VMD. Go to Extensions->Modeling->Membrane builder. There you can define the chemical composition and the size of your brand new membrane. The most important parameter here is the topology, which asks for which topology you are building this object. As you already know, it is better not to mess up with topologies and force fields so you should not mix them randomly.

If you are not able to locate the extension don’t worry, there are better ways to build membranes.

The G Protein Coupled Receptor

In this tutorial we will analyse a G Protein Coupled Receptor (GPCR) a member of this class of transmembrane protein that is involved in a huge variety of cellular processes. In normal conditions this receptor is in the inactive form, but when it is activated by external agents in the extracellular space it becomes active and triggers a huge variety of cellular responses and signalling cascades in the cytoplasm.

“108 G-protein-coupled receptors (GPCRs) are the targets of 475 (∼34%) Food and Drug Administration (FDA)-approved drugs and account for a global sales volume of over 180 billion US dollars annually” (source).

From an experimental perspective, obtaining structures of transmembrane proteins is an authentic nightmare. The pdb we are about to download was the first example of GPCR protein that was crystallised by structural biologists, thus paving the way to a new era of experimental structure determination.

Download the pdb 4GBR and load it in VMD

If we look at the protein structure (and at the pdb Id, never forget to do so) we can immediately see that the pdb was obtained with the help of a ligand and of an extracellular protein. Specifically, the experimentalists made a lisozyme bind to the extracellular domain of GPCR, while a small ligand CAU is attached to the core region of the protein.

You can convince yourself of this by creating the following representations in VMD (Graphics->Representations->Create New):

with only the lysozime protein, selecting the chain B and New Cartoon as Drawing method. Select the Coloring Method that you think it’s more suitable to this part of the structure;
with the GPCR protein, that is chain A and not resname CAU. I recommend using New Cartoon and Secondary structure as Drawing and Coloring methods, respectively;
the small ligand CAU can be visualized in Licorice. I think you know how to see only this pretty little guy.
let’s highlight now the atoms that have chain A and not resname CAU and resid > 301 and draw them in red (see Color id) selecting again New Cartoon as Drawing method. This alpha helix is the main actor of the cellular response: when the receptor receives a signal from the extracellular space the GPCR undergoes a conformational rearrangement to the active form that causes the exposition of this small helix. The alpha helix, once exposed, can bind to the G protein and then trigger the overall cellular response, which is specific to the type of signal received. As usual, biology is complicated but also amazing (3.0).

Our aim now is to save only the pdb with the channel, neglecting the two ligands. In order to do so, let’s go to the VMD main window. Go to File->Save Coordinates. Here you can decide to save one of the current selections to a new file. It is not necessary to have a pdb in output but you can select the format you like the most. Here we will choose the pdb format for the output and the second representation. Let’s call the output file 4GBR_alone.pdb.

Load the new pdb in VMD and check that you actually saved the right coordinates

Now let’s inspect a bit the differences between the file. Open both of them with your favourite text editor. You will see that saving the visualization with VMD is not optimal, since it makes you lose a lot of experimental information about the system.

For example, its experimentally determined disulfide bond network. Let’s retrieve the data about S-S bonds from our original pdb:

grep -rn 4GBR.pdb -e SSBOND

Save the Cysteine residues that are involved in disulfide bonds in a separate text file.

Charmm-gui

Charmm-gui is not a guy with a lot of charm, rather it is a platform that is appositely designed to perform interactive setups of complex Molecular Dynamics simulations. It is written in order to produce input parameters for all the main MD engines. Although it is based on Charmm (and on the charmm force field), it also includes the possibility of generating input files that employ AMBER potentials. In this part of the tutorial (for the last time in the course) we will stick to the Charmm force field.

From the main page of Charmm-gui go to the section called Input Generator: here you can see all the possible Charmm-gui workflows to generate input files for a huge variety of computational experiments. Among them, it is worth mentioning the following:

Glycolipid Modeler, which deals with glycolipids;
Ligand Reader & Modeler to read ligand structures and parametrise them with the Charmm force field;
Membrane builder (all the possible flavors of membrane available);
Martini input generators and backmapper, to perform Coarse-grained simulations with the Martini force field.

We are going to use the Membrane/Bilayer Builder to create the correct setup for our GPCR protein.

It will be a six-step process, as descibed at the bottom of the page:

Read protein coordinates: here you are supposed to provide your input file to Charmm-gui, either in form of a pdb or using an already oriented OPM file;
Protein-membrane orientation: in this crucial step Charmm-gui should orient the membrane and the protein in the correct way, depending on user’s input;
system size selection: how many molecules do we want in our system? the bigger the membrane the better, except for the computational overhead!
system build, where the fundamental constituents of the system (protein, membrane bilayer, water and ions) are created;
system assembly: everything is assembled together, like in LEGOs;
Equilibration: Charmm-gui provides six files for progressive equilibration of the membrane-protein-water system. We will analyse them more in detail later.

Let’s start!

Retrieve the pdb file 4gbr_alone.pdb that you generated before

If you lost it or if you were not able to generate it from VMD you can find a copy inside the tutorial folder available here.

1.Read coordinates

Upload the pdb file in the Membrane builder input cell, making sure you are selecting Protein/Membrane system as system type. Check the pdb box. Press the red arrow on the bottom right of the page to proceed to the next step.

In the next page you can select the model you want to use in your setup. This may be useful if you’re submitting a pdb extracted with NMR method, which can contain up to tens of different input structures, or if you have multiple snapshots of your protein obtained, for instance, with Molecular Dynamics.

Leave everything with default values and go to the next page

We don’t have multiple models in our pdb and we don’t want to cut it here (we already did it!).

Check the job id on the top right side of the page

With the job id you can retrieve your job easily using the Job Retriever.

Proceed to the next page. Here we can see many different PDB Manipulation Options. We could mutate the protein, add phosphoryl groups to some residues and so on and so forth. We leave the N-terminal patching as it is, so that the N-terminal and C-terminal regions are modelled with Charmm default parameters.

Select the disulfide bonds box and fill the network with the bonds present in 4gbr.pdb

Let’s proceed to the next step. We are almost at the end of the first part, we have only to decide how to orient our protein.

Here you will have few files wrapping up all what you’ve done so far.

In the psf file you can see the charmm topology of the protein, while step1_pdbreader.pdb contains the preprocessed pdb file.

Let’s check that the disulfide bonds are correctly included in the topology

In order to do so, retrieve the atom ID of the sulphur atoms involved in S-S bonds. Then go at the end of the bonds section to see if they have been included correctly.

If everything looks fine to you, it’s time to select the orientation of our protein. We would like to have our GPCR axis parallel to the z axis, so that the lipid bilayer can be built in the xy plane.

Check the Align the First Principal Axis Along Z box and go to the next page

2.Orientation

Here you can visualize the properly aligned protein in step2_orient.pdb file. In the page you should see the calculated Cross Sectional Area of the protein, with two bumps corresponding to the two pores on the sides of the membrane. You can even download all the information about cross sectional areas from the files box.

It’s time to choose the chemical composition of our membrane!

In my attempt I build a homogeneous lipid bilayer using phosphatidylcholine (POPC) lipids, but as you can see you have the possibility to design your membrane with the ratio of lipids that mimics the most the environment you’re simulating.

Off-topic: last year project was about recoverine, a membrane protein that is important for the human vision. Therefore the membrane was modelled using the estimated lipids ratio in the human retina cell membranes, namely 20% phosphatidylserine (POPS), 40% Phosphatidylethanolamine (POPE) and 40% of POPC.

Take home message: always talk to experimentalists or read a lot of articles before choosing the composition of your cell membrane

Select 90 (Angstroms) as xy size of the membrane and locate the lipids you want to include, specifying the desired ratios for both sides of the bilayer.

As an example try with 100% POPC.

Once you have specified the desired chemical composition click on the Show System Info panel: it will show you some parameters of your membrane. You should have around 100 POPC lipid molecules in both leaflets.

Let’s go to the next page: Determine system size.

3.Size determination

Download the file step3_packing.pdb

Don’t get scared!

This is only the pdb of the protein with the positions of the lipids, that are going to be added in the following step. At the four sides of the structure you can see small pieces of protein. These represents (I think) the Periodic Boundary Conditions of your system. You won’t have such artifacts in the output files.

Choose the replacement method to embed the protein in the membrane. Then leave 0.15 M as ion concentration.

Calculate the number of ions that will be present in the box.

Let’s proceed: Build Component!

4.Building components

Here you can see the assembled protein-mebrane system you just built in the previous step (step4_lipid.pdb).

Download the pdb file and load it in VMD.

Check that the average distance between lipid molecules is around 90 Angstroms.

Hint: use the 2 keyword!

At this point of the setup Charmm-gui automatically checks for Protein surface penetration and Lipid ring penetration. In the first case the engine checks for the presence of lipid tails that enter deeply inside the protein structure, while in the second one Charmm-gui monitors if any tail penetrates through a chemical ring structure. This ring can either be located on the protein (aromatic amino acids such as tryptophan) or on other lipids (cholesterol).

Charmm-gui developers say: Protein surface penetration and Lipid ring penetration are unphysical and will cause unstable simulations, and it is also very difficult to identify such penetrations solely by visual inspection.

If you do not have such artifacts you can proceed to the next step!

Assemble components!

5.Equilibration

Load the step5_assembly.pdb file.

Play around with it, checking that all the elements have been correctly included and patched together.

It’s time to select the desired force-field with which you want to simulate your system. As we said before you can choose between charmm36 and AMBER.

Check the Charmm36m force field.

Then you’re asked to specify the desired format (i.e. MD engine) for your simulation files.

Check the GROMACS format.

Leave the other options unchanged: we will do a NPT equilibration.

Since we really care about the health of your laptops, we won’t make you run any of these steps on your PC.

Proceed to the next page. You should have reached the end of your interactive setup! Here you can download a zip file with all the pieces of your system.

Some fundamental concepts and a bit of analysis

Open all the mdp files you retrieved from charmm-gui

Their structure is more or less the same, but they contain some differences. Here we are going to highlight the most important ones

diff step6.1_equilibration.mdp step6.2_equilibration.mdp

This command shows you the difference between the first two mdp files. Both have position restraints (see the first line), but with different constants, which are higher in the first equilibration.

Do you have an explanation for this choice?

The other difference is that in the first file the velocities are generated (gen_vel = yes), while the second simulation is supposed to be a continuation of the first (continuation = yes) and there’s no gen_vel.

Both simulations do not have pressure coupling but only temperature coupling (NVT).

diff step6.2_equilibration.mdp step6.3_equilibration.mdp

Again, the position restraints are weaker. pcoupl = berendsen means that we are coupling the system to the Berendsen barostat, with a reference pressure (ref_p) of 1.0 bar.

The most important parameter for pressure coupling in membrane is pcoupltype = semiisotropic. This means that the pressure coupling is not completely isotropic, rather we impose an uniform scaling of x-y box vectors, while the z component is independent. We deal with a squared membrane, so we don’t want to obtain a rectangular one at the end of our setup.

In most simulations of membranes the pressure coupling is limited to the z component, thus keeping fixed the x-y shape of the system.

diff step6.3_equilibration.mdp step6.4_equilibration.mdp

The constraints are weakened and we have a change of the timestep (dt = 0.002 instead of dt = 0.001). The system is stable enough to switch to the standard (2 femtoseconds)timestep. The number of steps is also doubled.

diff step6.4_equilibration.mdp step6.5_equilibration.mdp

diff step6.5_equilibration.mdp step6.6_equilibration.mdp

Only the constraints are weakened.

Let’s see the differences between the last equilibration file and the production mdp:

diff step6.6_equilibration.mdp step7_production.mdp

Here the situation changed:

No constraints anymore (the ones on hydrogen atoms are always present);
Weak coupling Berendsen thermostat is substituted by Nose-Hoover temperature coupling, which correctly samples the canonical ensemble;
Weak coupling Berendsen barostat is replaced by Parrinello-Rhaman pressure coupling, which correctly samples the isothermal-isobaric ensemble.

Never use the Berendsen barostat for actual data production!

Analysis

In the tutorial folder you can find equilibrated the system for you (step6.6_equilibration.gro) and 2 nanoseconds of simulation of the equilibrated system. You can download it from (step7_production_noPBC.xtc).

Load the gro file in VMD

You can see the equilibrated system. You can make a visual comparison with the output file you got from charmm-gui. This file has been through different steps of minimisation, NVT equilibration and NPT, so it is more realistic than the first. Nonetheless, with such big systems it is common practice to consider the first nanoseconds of the production run as a further equilibration phase, thus starting actual data collection later on.

Load the trajecory on the gro file in VMD

As usual you can do it going to File/Load Data Into Molecule and selecting the file step7_production_noPBC.xtc.

Unfortunately, it’s only a few frames since the supercomputer of the university crashed :(.

Let’s create a nice representation of our system by removing the water and highlighting the protein and lipids in a different way. Create a new representation (Graphics/Representations/Create Rep):

where we select the atoms of the protein (protein) using New Cartoon and Secondary Structure as Drawing and Coloring method respectively;
where we draw the lipid bilayer using Polyedra as Drawing method;
where we draw the ions (ions) using VdW as Drawing method.

Now we will see how to save to disk only a fraction of the atoms of the system directly from VMD. In this case let’s say that we want to save only the atoms of our protein.

Make sure you selected your trajectory;
Go to File/Save Coordinates;
save the first frame of your trajectory (First = Last = 0) writing only the protein atoms. Save it in a gro format (you can call it step6.6_equilibration_protein.gro);
write the binary trajectory keeping the same selection but changing the Last parameter to the number of frames and selecting dcd as Output format. Call it step7_production_protein.dcd

Analyse the behavior of the RMSD with respect to frame 0 using MDAnalysis or VMD

Further Notes

Charmm-gui 2016 article