{ "cells": [ { "cell_type": "markdown", "id": "f080d3b5", "metadata": {}, "source": [ "# POPC:POPE:CHOL bilayer (Coarse-Grained).\n", "\n", "**Last executed:** Jan 3rd, 2024 with MDAnalysis 2.7.0 and Python 3.11\n", "\n" ] }, { "cell_type": "markdown", "id": "047539da", "metadata": {}, "source": [ "**Packages required:** MDAnalysis, NumPy.\n", "\n", "**Packages for visualization:** Matplotlib, SciPy.\n", "\n", "In this tutorial, we are going to use MembraneCurvature to derive surfaces and calculate curvature of a lipid bilayer of lipid composition POPC:POPE:CHOL and a 5:4:1 ratio from a Molecular Dynamics (MD) simulation carried out using the [Martini](http://cgmartini.nl/) force field." ] }, { "cell_type": "code", "execution_count": 1, "id": "148160aa", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "MDAnalysis : INFO MDAnalysis 2.7.0 STARTED logging to 'MDAnalysis.log'\n", "MDAnalysis : INFO MDAnalysis 2.7.0 STARTED logging to 'MDAnalysis.log'\n", "MDAnalysis : INFO MDAnalysis 2.7.0 STARTED logging to 'MDAnalysis.log'\n" ] } ], "source": [ "import MDAnalysis as mda\n", "from membrane_curvature.base import MembraneCurvature\n", "from membrane_curvature.tests.datafiles import MEMB_GRO, MEMB_XTC\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "b42e011f", "metadata": {}, "source": [ "Because the system included in this tutorial uses the [Martini](http://cgmartini.nl/) coarse grain representation, we prefer to suppress the warnings from MDAnalysis. Please note that this shouldn't be the standard approach for your code as warnings may give you some heads-up about issues with your system. " ] }, { "cell_type": "code", "execution_count": 2, "id": "a31001a0", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "id": "06cd8436", "metadata": {}, "source": [ "---\n", "\n", "For a better understanding of how MembraneCurvature help us to calculate membrane curvature, let's split the process into five main steps:\n", "\n", "[1. Load MDA Universe](#1.-Load-MDA-Universe)\n", "\n", "[2. Select Atoms of Reference](#2.-Select-Atoms-of-Reference)\n", "\n", "[3. Run MembraneCurvature](#3.-Run-MembraneCurvature)\n", "\n", "[4. Extract Results](#4.-Extract-Results)\n", "\n", "[5. Visualize Results](#5.-Visualize-Results)\n", "\n", "---\n" ] }, { "cell_type": "markdown", "id": "a5e1e45e", "metadata": {}, "source": [ "## 1. Load MDA Universe" ] }, { "cell_type": "markdown", "id": "e3c217e5", "metadata": {}, "source": [ "We start by loading our molecular dynamics output files (coordinates and trajectory) into an MDAnalysis Universe: " ] }, { "cell_type": "code", "execution_count": 3, "id": "b3beb4b6", "metadata": {}, "outputs": [], "source": [ "universe = mda.Universe(MEMB_GRO, MEMB_XTC)" ] }, { "cell_type": "markdown", "id": "8971f3f4", "metadata": {}, "source": [ "A quick check of our Universe gives us a better idea of the system used in this tutorial. To check the dimensions of the simulation box, number of lipid residues, total number of beads, as well as the number of frames in the trajectory contained in our Universe, we can run: " ] }, { "cell_type": "code", "execution_count": 4, "id": "4f2f1f1f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Universe info:\n", "\n", "Box dimensions[x,y,z]: [241.12502 241.12502 233.69261].\n", "\n", "2046 lipid residues and 23736 beads.\n", "\n", "Lipid types: \n", "-CHOL\n", "-POPE\n", "-POPC\n", "\n", "The trajectory includes 11 frames.\n" ] } ], "source": [ "print(\"Universe info:\")\n", "print(\"\\nBox dimensions[x,y,z]: {}.\" \\\n", " \"\".format(universe.dimensions[:3]))\n", "\n", "print(\"\\n{} lipid residues and {} beads.\" \\\n", " \"\".format(universe.residues.n_residues, universe.residues.n_atoms))\n", "\n", "print(\"\\nLipid types: \")\n", "[print(\"-\"+lipid) for lipid in set(universe.residues.resnames)]\n", "\n", "print(\"\\nThe trajectory includes {} frames.\"\\\n", " \"\".format(universe.trajectory.n_frames))" ] }, { "cell_type": "markdown", "id": "bcca59be", "metadata": {}, "source": [ "It probably gets easier if we can have a quick visual inspection of the system. To visualize our system, we can make use of the [NGL Viewer](https://github.com/nglviewer/nglview)." ] }, { "cell_type": "code", "execution_count": 5, "id": "cde8d33e", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1fcbdd9c59284841aa4671ee5a7e8010", "version_major": 2, "version_minor": 0 }, "text/plain": [] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import nglview as nv\n", "\n", "def color_by_resname(system, zoom, color_, box=False):\n", "\n", " \n", " view = nv.show_mdanalysis(system)\n", " view.add_representation('ball+stick', selection='not CHOL', radius=0.5, color=\"grey\")\n", " view.add_representation('ball+stick', selection='.PO4', radius=2.00, color=color_)\n", " view.add_representation('ball+stick', selection='not POPC and not POPE', radius=0.5, color=\"yellow\")\n", " view.add_representation('ball+stick', selection='.ROH', radius=2.00, color=\"orange\")\n", " \n", " if box == True:\n", " view.add_unitcell()\n", " \n", " view.camera='orthographic'\n", " view.control.zoom(zoom)\n", " view.control.rotate(\n", " mda.lib.transformations.quaternion_from_euler(\n", " -np.pi/2, np.pi/3, np.pi/12, 'rzyz').tolist())\n", " \n", " \n", " return view" ] }, { "cell_type": "code", "execution_count": 6, "id": "4593c2ba", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2b18d79aa79c4bd3a5780d9715d09088", "version_major": 2, "version_minor": 0 }, "text/plain": [ "NGLWidget(max_frame=10)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "color_by_resname(universe, 0.5, \"brown\")" ] }, { "cell_type": "markdown", "id": "d1543026", "metadata": {}, "source": [ "The NGL Viewer enables us to identify two main types of lipids: phospholipids, POPC and POPE lipids, shown in grey colour; and sterols, CHOL lipids, shown in yellow. All phospholipid head groups are shown in red, and sterol head groups are shown in orange." ] }, { "cell_type": "markdown", "id": "60ca1249", "metadata": {}, "source": [ "Since membraneCurvature uses the dimensions of the box to set the dimensions of the grid (See [Algorithm](https://membrane-curvature.readthedocs.io/en/latest/source/pages/Algorithm.html#set-grid)), it is important to run a preliminary check to confirm that our system does not display significant box size fluctuations over time. \n", "\n", "We can quickly check the fluctuations of our simulation box for the `x` and `y` dimensions with:" ] }, { "cell_type": "code", "execution_count": 7, "id": "ea7c4ab8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 [241.12502 241.12502]\n", "1 [241.33353 241.33353]\n", "2 [240.43678 240.43678]\n", "3 [240.13823 240.13823]\n", "4 [240.14366 240.14366]\n", "5 [240.20067 240.20067]\n", "6 [239.7728 239.7728]\n", "7 [240.65321 240.65321]\n", "8 [239.9748 239.9748]\n", "9 [240.36978 240.36978]\n", "10 [240.5304 240.5304]\n" ] } ], "source": [ "for ts in universe.trajectory:\n", " print(ts.frame, ts.dimensions[:2])" ] }, { "cell_type": "markdown", "id": "dbc17150", "metadata": {}, "source": [ "With a fluctuation of $~1\\%$ in each dimension, we know that the grid set by MembraneCurvature using the `u.dimensions` will contain all the elements in our system.\n" ] }, { "cell_type": "markdown", "id": "a8e069ff", "metadata": { "tags": [] }, "source": [ "Now that we have a better idea of what the system used in this tutorial contains and we have verified the fluctuations of the box are not significant, we can go to step 2 and select our atoms of reference." ] }, { "cell_type": "markdown", "id": "d44185f5", "metadata": {}, "source": [ "## 2. Select Atoms of Reference" ] }, { "cell_type": "markdown", "id": "30530fd7", "metadata": {}, "source": [ "The [MembraneCurvature algorithm](https://membrane-curvature.readthedocs.io/en/latest/source/pages/Algorithm.html) uses an `AtomGroup` as a reference to derive surfaces. Therefore, it's important to identify what constitutes a sensible selection for an `AtomGroup` of reference. Keep in mind that choosing atoms that display flip-flop during your simulation is not a good idea. The atoms in your `AtomGroup` of reference should remain in the same leaflet over the simulation time.\n", "\n", "Typically in lipid bilayers, phospholipid headgroups are the most straightforward `AtomGroup` from which to derive surfaces. In this tutorial, we are going to choose phospholipid headgroups of phosphate beads. In the Martini representation, these beads are called `PO4`. Then, we can select the phospholipid head groups (`PO4` beads) in our Universe with:" ] }, { "cell_type": "code", "execution_count": 8, "id": "4f32b1ae", "metadata": {}, "outputs": [], "source": [ "PO4_headgroups = universe.select_atoms('name PO4')" ] }, { "cell_type": "markdown", "id": "30fb6345", "metadata": {}, "source": [ "