As computational biologists have generated and studied ever larger and more complex biological systems in silico, modeling large-scale lipid membranes with carefully considered curvature/geometry has become increasingly critical. We here present a computer program called LipidWrapper, written by Jacob Durrant in the lab of Rommie E. Amaro, capable of creating curved membrane models with geometries derived from various possible sources, both experimental and theoretical. We are hopeful that this utility will be a useful tool for the computational-biology community. A copy can be downloaded free of charge from SourceForge.
1. Defining the Mesh
The LipidWrapper algorithm begins by considering mesh points that define a bilayer surface (Figure 1A). These mesh points can be generated in four ways: the user can provide a PDB file with atom coordinates that define the points; the points can be defined by a gray-scale image, where black is topologically "low" and white is topologically "high"; the user can provide an equation that defines z when given (x, y) coordinates; and the user can provide a collada DAE file with both mesh points and surface polygons (triangles) specified.
For the first three of these options, a tessellated/triangulated surface must be generated from the mesh points. All mesh points are projected onto the x-y plane, and the resulting set of coplanar points is subjected to a 2D Delaunay triangulation. This triangulation is then projected back onto the 3D mesh (Figure 1B). As the triangulation requires an intermediate 2D projection, bilayer models than cover concave or enclosed surfaces are not permitted. For complex shapes, including concave and enclosed surfaces, the user can provide a collada DAE file with the surface tesselation/triangulation already specified. LipidWrapper has been specifically tested with collada files exported from the free, open-source program Blender, a full-featured 3D modeling package that has been used in diverse projects ranging from animation to video-game development to scientific imaging. Once a three-dimensional model has been loaded into Blender, the user need only switch to the Edit Mode, triangulate the surface mesh of the desired object (Control-T after selection), and export the modified object as a collada DAE file. When this file is passed to LipidWrapper, the algorithm will skip its own triangulation step entirely, relying instead of the surface triangulation specified by the DAE file.
2. Positioning Bilayer Models onto the Tessellated Triangles
With surface triangulation complete, the program next considers a user-specified PDB model of a planar lipid bilayer that is oriented parallel to the x-y plane, with the extracellular side facing the positive z direction (Figure 2A). The plane that bisects this bilayer is determined by averaging the z values of all lipid-residue headgroups, as judged by the position of the lipid phosphorus and cholesterol hydroxyl oxygen atoms (by default). Copies of the previously tessellated surface triangles are then randomly placed on this lipid-bisecting plane such that they are entirely contained within the region defined by the headgroups (Figure 2A). The transformation matrices required to rotate and translate the triangle copies back onto their original 3D surface triangles are then calculated. These same matrices are then applied to the lipid molecules that fall within each placed triangle, effectively positioning them onto the three-dimensional surface (Figure 2B and 2C).
3. Resolving Steric Clashes
Next, the user can optionally request that LipidWrapper delete lipid residues as needed to resolve steric clashes. As identifying steric clashes ultimately requires computationally expensive pairwise distance comparisons, the program goes to extraordinary measures to keep the number of distance calculations to a minimum. First, Lipidwrapper assumes that lipids belonging to the same tessellated triangle do not clash, as they are all derived from a single bilayer model that was presumably carefully prepared. Second, as lipids belonging to distant tessellated triangles are unlikely to clash, LipidWrapper next seeks to identify proximate triangles that share corners and/or have centers that come within a user-defined distance of each other (50.0 Å by default, Figure 3). Third, as only lipids near triangle edges can clash, lipids with headgroups that, when projected onto their respective triangle planes, do not come within a user-specified distance of the triangle edge (25.0 Å by default), are discarded (Figure 3A, shown in purple/pink). Finally, as lipids with distant headgroups are unlikely to clash, lipids from proximate triangles with headgroups that are greater than a user-specified distance apart (50.0 Å by default) are likewise discarded.
Once the above pruning techniques have been completed, the remaining lipids, judged to be potential clashers, are next considered. First, the program checks to see if the minimum coordinate of one lipid is so much larger than the maximum coordinate of its potential clash partner (in any of the three dimensions) so as to render a clash impossible. If not, all atoms in the overlapping region of bounding boxes encompassing each lipid are identified. A pairwise distance comparison is then performed on these atoms, and if any two of them come within a user-specified distance (2.0 Å by default), the two lipids are judged to clash. To resolve all lipid clashes, lipids are removed one by one, starting with the lipid that clashes with the greatest number of its neighbors, until all steric clashes have been resolved (Figure 3B).
4. Filling in the Resulting Bilayer Holes
Lipid removal invariably produces “holes” in the bilayer. The user can optionally instruct LipidWrapper to fill these holes with additional lipid molecules. First, the region of each tessellated triangle near the triangle edges (where holes resulting from lipid deletions are likely to occur) is flooded with equidistant points. For each of these coplanar points, two new points are generated both above and below the triangle plane, positioned roughly at the level of the lipid headgroups. After points too close to the surrounding lipids are deleted, the remaining points are located in bilayer holes (Figure 4A, in green). LipidWrapper then iterates through each of these points a user-specified number of times (10 by default). A lipid molecule is placed at each point. If the molecule does not clash with any of its neighbors, it is incorporated into the bilayer (Figure 4B, in blue).
LipidWrapper accepts a number of commandline parameters:
Methods for creating a surface mesh =================================== --surface_equation: Generate a surface mesh from a python-formatted equation defining z, given x and y. The --min_x, --max_x, --min_y, and --max_y parameters are used to specify the region over which the function should be evaluated. The --step_x and --step_y parameters define the x-y distance between adjacent points. Python functions from the math, numpy, and scipy modules can be used. Example: --surface_equation "z = 250*numpy.sin(x*x/60000 +y*y/60000)" --surface_filename: If this parameter specifies a file with the PDB extension, a surface mesh is generated from the coordinates of the PDB atoms. Example: --surface_filename mymesh.pdb --surface_filename: If this parameter specifies a file that does not have the PDB extension, the file is assumed to be a gray-scale image, where black represents regions that are topologically low, and white represents regions that are topologically high. The --min_x, --max_x, --min_y, and --max_y parameters are used to specify the region where the mesh should be generated. The --step_x and --step_y parameters define the x-y distance between adjacent points. The --max_height parameter determines the height of the bilayer model at those locations where the image is white; black regions are assigned a height of 0. This feature is only available if the python PIL module has been installed on your system. Example: --surface_filename mymesh.png The initial lipid model ======================= --lipid_pdb_filename: This parameter specifies a PDB file containing an all-atom model of a planar lipid bilayer. LipidWrapper will wrap this lipid around the user-generated mesh. Example: --lipid_pdb_filename lipid.pdb --lipid_headgroup_marker: A unique atom representing the headgroup of each lipid residue must be specified. The --lipid_headgroup_marker accepts a comma-separated lists of atom specifications (RESNAME_ATOMNAME). If either RESNAME or ATOMNAME is omitted, any value will be accepted. By default, LipidWrapper identifies lipid headgroups by looking for any atom named "P" (_P) or any atom named "O3" belonging to a cholesterol molecule (CHL1_O3). Example: --lipid_headgroup_marker "_P,CHL1_O3" Methods for resolving lipid clashes =================================== --delete_clashing_lipids: It's common for lipids to sterically clash at the interface of two adjacent surface-mesh tessellated triangles. If this parameter is set to TRUE, any clashing lipids are deleted. Example: --delete_clashing_lipids TRUE --clash_cutoff: If you do choose to delete clashing lipids, this parameter determines how close two atoms must be (in Angstroms) to constitute a steric clash. Example: --clash_cutoff 2.0 --fill_holes: Deleting lipids often leaves holes in the membrane. If this parameter is set to TRUE, LipidWrapper tries to fill the hole. Example: --fill_holes TRUE --fill_hole_exhaustiveness: Essentially, how long LipidWrapper should try to fill the holes. Example: --fill_hole_exhaustiveness 10 --clashing_potential_margin: Lipid clashes occur at the edges of adjacent tessellated triangles. If these triangles are very large, it's faster to only check for clashes and holes near the triangle edges. This variable specifies how far from the edges, in Angstroms, that LipidWrapper should look for clashes and holes. Example: --clashing_potential_margin 25.0 --very_distant_lipids_cutoff: LipidWrapper determines if two lipids clash by comparing the distance between every atom in the first lipid with every atom in the second lipid. This can be computationally expensive. However, sometimes two lipids are so distant from each other, that it's obvious there are no clashes, making the pair-wise comparison unnecessary. Before performing this expensive pair-wise comparison, LipidWrapper calculates the distance between one atom of each lipid. If this distance is greater than this user-specified cutoff, the program will simply assume there are no clashes. WARNING: Remember to consider the width of your lipid bilayer when choosing this value. Adjacent lipids on opposite sides of the bilayer can seem distant when considering the distance between their headgroups, for example. Example: --very_distant_lipids_cutoff 50.0 --triangle_center_proximity_cutoff_distance: Lipid steric clashes/holes typically occur between lipids that belong to adjacent tessellated triangles. However, if tessellated triangles are small enough, clashes are possible between lipids that belong to non-adjacent triangles as well. Consequently, in addition to checking for adjacency, LipidWrapper also checks the distance between the triangle centers, using this user-specified value as a cutoff. Example: --triangle_center_proximity_cutoff_distance 50.0 --memory_optimization_factor: When the tessellated triangles are very large and consequently contain many individual lipids, the extensive pairwise distance comparisons required can result in memory errors. This parameter tells lipid Wrapper to divide the list of atoms being compared into smaller chunks. The pairwise distance comparison is performed piecewise on each chunk-chunk pair and so uses less memory, albeit at the expensive of speed. Only increase the value of this parameter if you run into memory errors. Example: --memory_optimization_factor 1 Additional options ================== --number_of_processors: Using multiple processors can significantly increase the speed of the LipidWrapper algorithm. Example: --number_of_processors 8 --show_grid_points: Aside from producing PDB coordinates for lipid atoms, additional coordinates will be appended to the bottom of the output containing "atoms" named "X" that specify the location of the surface mesh points. Example: --show_grid_points TRUE --create_triangle_tcl_file: A separate file named "triangles.tcl" will be generated containing a tcl script that can be run in VMD to visualize the mesh surface. Example: --create_triangle_tcl_file TRUE --output_directory: If an output directory is specified, all LipidWrapper output files, as well as additional files representing the intermediate steps required to build the final bilayer, will be saved in that directory. Example: --output_directory ./my_output/ --use_disk_instead_of_memory: For very large systems, storing the growing model in memory can be problematic. If this parameter is set to TRUE, the growing model will be stored on the hard disk instead. However, expect longer execution times if this parameter is set to TRUE. Example: --use_disk_instead_of_memory TRUE --compress_output: Depending on the user options selected, LipidWrapper output can require a lot of disk space. If this parameter is set to TRUE, the output will be automatically compressed using the gzip algorithm (Lempel-Ziv coding LZ77). The files can be uncompressed with the UNIX gunzip utility, or similar Windows-based packages. Example: --compress_output TRUE
Generate a bilayer from an equation, using 24 processors. Store the growing model on the hard disk rather than in memory. Resolve steric clashes and fill the resulting bilayer holes.python lipidwrapper.py --surface_equation "z = 250*numpy.sin(x*x/60000 +y*y/60000) * (-numpy.sqrt(x*x+y*y)/(560 * numpy.sqrt(2)) + 1)" --min_x -500 --max_x 500 --min_y -500 --max_y 500 --step_x 25 --step_y 25 --lipid_pdb_filename lipid_example.pdb --delete_clashing_lipids TRUE --fill_holes TRUE --clash_cutoff 1.0 --output_directory example1_output --number_of_processors 24 --use_disk_instead_of_memory TRUE --show_grid_points TRUE --create_triangle_tcl_file TRUE --clashing_potential_margin 15.0 --memory_optimization_factor 1
Generate a bilayer from a series of points in a PDB file, using 24 processors. Store the growing model in memory, not on the hard disk (faster). Resolve steric clashes and fill the resulting bilayer holes.python lipidwrapper.py --surface_filename influenza_mesh.pdb --lipid_pdb_filename lipid_example.pdb --delete_clashing_lipids TRUE --fill_holes TRUE --clash_cutoff 1.0 --output_directory example2_output --number_of_processors 24 --use_disk_instead_of_memory FALSE --show_grid_points TRUE --create_triangle_tcl_file TRUE --clashing_potential_margin 15.0 --memory_optimization_factor 1
Generate a bilayer from an image, using 24 processors. Store the growing model in memory, not on the hard disk (faster). Resolve steric clashes and fill the resulting bilayer holes.python lipidwrapper.py --surface_filename face.png --min_x -500 --max_x 500 --min_y -500 --max_y 500 --step_x 25 --step_y 25 --lipid_pdb_filename lipid_example.pdb --max_height 150 --delete_clashing_lipids TRUE --fill_holes TRUE --clash_cutoff 1.0 --output_directory example3_output --number_of_processors 24 --use_disk_instead_of_memory FALSE --show_grid_points TRUE --create_triangle_tcl_file TRUE --clashing_potential_margin 15.0 --memory_optimization_factor 1
Generate a bilayer from a blender-generated DAE file, using 24 processors. Store the growing model in memory, not on the hard disk (faster). Resolve steric clashes and fill the resulting bilayer holes.python lipidwrapper.py --surface_filename concave_model.dae --lipid_pdb_filename lipid_example.pdb --delete_clashing_lipids TRUE --fill_holes TRUE --clash_cutoff 1.0 --output_directory example4_output --number_of_processors 24 --use_disk_instead_of_memory FALSE --show_grid_points TRUE --create_triangle_tcl_file TRUE --clashing_potential_margin 15.0 --memory_optimization_factor 1
With some versions of python, running LipidWrapper on multiple processors produces the following warning:Exception AssertionError: AssertionError() in <Finalize object, dead> ignored
This warning is caused by a known python bug. Fortunately, it does not affect the LipidWrapper output and so can be ignored.