Please Cite Source: Wildermuth, K.; Monje-Galvan, V.; Warburton, L. M.; Klauda, J. B., Effect of membrane lipid packing on stable binding of the ALPS peptide. Journal of Chemical Theory and Computation 2019, 15, 1418-1429. (https://pubs.acs.org/doi/10.1021/acs.jctc.8b00945)
This analysis focuses on the study of membrane lipid packing defects – areas of low density lipid packing on the surface of a membrane that expose the nonpolar core of the membrane. Packing defects arise from dispersion of membrane constituents due to external forces such as imposed surface tensions or chemical interactions with other molecules. The analysis of packing defects is useful when quantifying the membrane surface is of interest. For example, I used packing defect quantification to describe the dynamic binding event between the Osh4 ALPS-like peptide and a membrane surface, and to describe the effects of using the HMMM simulation technique.
This process uses image processing techniques to identify packing defects in specially constructed images of the membrane surface.The process is broken into three steps illustrated in Figure 1: (1) constructing images in VMD, (2) identifying packing defects using Python and OpenCV, then (3) analyzing the results. Images produced in step 1 will look something like what is shown in Figure 1 where yellow “blobs” shown in the bottom image are the packing defects. In step 2, each packing defect on the membrane is identified as an individual object and tagged with descriptors such as time, size, and coordinates. The descriptors that are assigned to each defect can be expanded upon by future users. For step 3, I will present some analysis I have developed for my application, but the possibilities reach beyond what I describe here.
It is also worth noting that I developed this process to analyze the binding event of a peptide to a membrane, thus these scripts also incorporate identification and analysis of a protein in the simulation. The process can be simplified to study a membrane-only system (see notes in detailed instructions on where important changes should be made to study a membrane-only system).
Since this code has not been used by anyone but myself as of writing this, feel free to forward questions directly to me at kywildermuth@gmail.com.
Figure 1: Defect Analysis Workflow
Using the tachyon ray tracer in VMD, this script will produce a series of images that expose packing defects on the surface of a membrane from a MD trajectory. The images can be used for quantification and analysis of packing defects in a membrane and their relationship to a binding protein. This script produces a file architecture that is used in accompanied defect quantification and analysis scripts. Most of the user input changes will occur in this step.
The following files should be in the same directory (I will refer to the main directory as home
)
home
using an automated directory architecture created by the script. Two image sets (protein and lipid) for each leaflet of the membrane are created, separated into two directories within home
(top, bottom) for a total of four image sets. An optional side directory will be created if this image set is selected to be written. home
. Contains output lines from VMD operation. Useful for troubleshooting the script. defect_tachyon_scav.tcl
. CREATE REPRESENTATIONS
section of the script to match the species in your system. polar_PS
, polar_PE
). If a lipid in your system is not already listed, it must be defined. Default definition of polar is to include all atoms in the head group down to and including the second carbon in the lipid chain. custom_sel
variable if you want to track a specific selection during the trajectory. For example, if you wanted to track the TRP alpha carbon of a peptide, you could set custom_sel atom selection to resname TRP and name CA
. The associated analysis scripts are constructed to read a custom_sel file, thus if you do not want to track a custom selection, just in the selection with something arbitrary and ignore its output. xsize
, ysize
: pixel dimensions of outputted images (default 600×600), recommended to keep 1:1 ratio.home
: filepath to directory that contains input files pbc wrap -centersel "protein" -center com -compound residue -all
-shiftcenter
moves the box. Increase the magnitude of the z coordinate in the appropriate direction if necessary. render Tachyon $filename "/homes/jbklauda/vmd/vmd-1.9/lib/tachyon/tachyon_LINUXAMD64" $filename -lowshade -res 600 600 -format BMP -o $filename.bmp
-lowshade
option removes shade so that the image produced can be easily processed-res
option can be changed if creating images larger than 600×600.
From command line in home
, call defect_tachyon.csh
. The cshell file calls VMD.
Now we will run an image processing script on the prepared image sets. This step utilizes python and a computer vision library called OpenCV. Both softwares are available on DT1, and these directions will assume use of DT1 for operation of the analysis. Though in theory, this analysis can be run on any system with Python (and numpy) and OpenCV. Furthermore, I assume the use of OpenCV v.3.1.0 and Python v.2.7.8
A visual check of some of the images should be conducted before running the code to make sure things are in order. For example, make sure the green box is visible in the frame and, if there is a protein in the system, make sure the green box is centered on the protein. The protein should not wrap around the edge of the PBC. Protein images should contain only the protein, and lipid images should contain only the lipid.
home
which contains the directories built by the previous VMD script. Note that the files are very similar, though blob_top.py
runs analysis for the top leaflet of the membrane and blot_bot.py
runs analysis for the bottom leaflet of the membrane parent_dir
to the same as home
used in the VMD scriptbox_x
and box_y
into both blob_top.py
and blob_bot.py
. This is very important since this information is used to convert unit pixels to Å^2is_contour_bad
function
From the command line and while in home
, load the two software modules and then simply call the python scripts:
module load python/2.7.8 module load opencv/3.1.0 python blob_top.py
By completing steps 1 and 2, all packing defects have been identified. Now this information can be used to run analysis. Defect analysis will differ case-by-case. For example, if you are measuring defects to analyze a binding event of a protein, you might plot a time series of the local packing defect size relative to the protein area, which you would expect to increase when the peptide binds (Figure 2). If you want to determine how membrane properties differ between two simulations, you might run a statistical test on the set of identified packing defects of each simulation to see if there is a difference in a variable, such as the probability that a packing defect with an area larger than some threshold exists (Figure 3). These analysis can easily be conducted from the output files of the blob_top.py
and blop_bot.py
scripts. More complex analysis can also be created using the information provided from steps 1 and 2. Below is a process that outlines an example.
Figure 2: Local packing defect factor time series
Figure 3: Probability of a defect larger than a threshold existing for three different membrane types
Some proteins bind to membranes at large packing defects. One way to show that this is true is by comparing the probability of the peptide moving towards the membrane to the conditional probability of the peptide moving towards the membrane given that a large packing defect exists near the peptide (theoretical details are outlined in a paper published by our group). Using information from the packing defect analysis, we can run a script to compute the necessary probabilities to conduct this analysis.
The script truDiff.py
can be run in the data_files
directory created by the VMD script. This script will measure the minimum z-distance between the protein and the phosphate layer of a specified leaflet of the membrane. In the script, specify the final binding leaflet
('top' or 'bottom'), and the number of frames in the DCD (num_frame
and prebind_limit
). Execute the script from the command line by using the following code while in the data_files
directory:
module load python/2.7.8 python truDiff.py
There will be a short prompt and the computation will complete quickly. An output file binding_distances.dat
will be made which contains the necessary minimum binding distance time series, as well two other distance time series described in the script's header
Next, the script defect_peptide_frequency.py
will count the number of times the peptide moves towards the membrane by a specified distance (m), how many instances a defect exists larger than a specified area (d), and how many times both of those instances occur simultaneously. These frequencies can be used to compute the probabilities and conditional probabilities described above. Options in the script allow the user to set the range of defect area and movement variable space, distance cutoffs (so you can analyze interactions only in a specific distance range between the membrane and the peptide), and trajectory cutoffs (so you can analyze interactions only in a specific time range, for example while the peptide is unbinded). Choose freq_path
to be the local_defect
data file for the leaflet of interest (typically the binding leaflet) and make sure the dist_path
is set to the truDiff output file. Execute the script from the command line by using the following code while in the data_files
directory:
module load python/2.7.8 python defect_peptide_frequency.py
module load
command does not have to be executed if the python module has already been loaded. freq_stat_total.dat
will contain counts of each event described above in a large table. Probabilities can be calculated from these tallies, and various results plotted (insert example excel sheet)