We are now almost ready to simulate the system of SBCG BAR domains with SBCG membrane. In this section we will discuss first how to write a NAMD configuration file for an SBCG system. We will then perform the simulation, and discuss the file outputs and simulation result.
To perform exersises, navigate to the folder 5_simulation/.
Since SBCG was designed to be compatible for NAMD, an SBCG configuration file looks similar to a normal, all-atom, NAMD configuration file that you might have used before.
1. A sample configuration file sim.conf has been prepared for you in the directory example-output/. Copy it to the folder where you want to run the simulation, and open it with a text editor. Remember to create in the same folder subfolder output and input, and add you CG files to the folder input analogously to how it is done in example-output/. Note that you will need here an SBCG parameter file for lipids, which we have not used before. For this purpose, copy the file example-output/input/lipid-ion.par to your input folder.
2. The configuration file contains many options (entries in the first column), followed by their parameters (entries in the second column) specifically chosen for the simulated system. Assuming readers already have experience with NAMD simulations, here we will only go through those options that require special adjustments for an SBCG system. New NAMD users are encouraged to consult the NAMD Tutorial and NAMD User's Guide.
3. In the text editor displaying the content of sim.conf, scroll down to the section # Force-Field Parameters. Note all lines beginning with # are comments ignored by NAMD.
Under # Force-Field Parameters, you will find six simulation options that might need different parameters than those of an all-atom simulation. These options define how you want the interactions between beads to be computed. Since each SBCG bead represents a cluster of atoms, parameters such as cutoff, switchdist, and pairlistdist would typically have bigger values than those for an all-atom simulation. The values listed here have been tested to work well for the SBCG BAR domain system, and can be used as starting values for other SBCG systems. Note, however, the cutoff parameter should always be bigger than the longest bond length in your simulation system, otherwise your simulation will crash.
4. Scroll down to the section # Integrator Parameters.
The parameter timestep has a value of 100.0, implying that the integration timestep of the simulation is 100 fs/step. A typical all-atom simulation uses 1 or 2 fs/step, hence the SBCG gives a speedup of 100 from the choice of integration timestep alone. The choice of the timestep depends on how fast beads are moving in the simulation, and, thus, the maximal timestep possible (so that the simulation does not crash) is determined by the strength of interactions, e.g., stiffness of bonds, as mentioned above. If your simulation crashes with a timestep of 100 fs/step, starting the simulation with a shorter timestep might fix the problem. Then timestep can be increased when the system becomes stable later in the simulation.
5. In the cellBasisVector1, cellBasisVector2, cellBasisVector3, and cellOrigin parameters, one specifies the periodic boundary of the system. Compare these values to the system size you have measured in the previous section.
You should see that only cellBasisVector2 matches the y-dimension of your system, while cellBasisVector1 and cellBasisVector3 are both much larger than the x- and z-dimension of the system. This is because we only want the system to be continuous in the y-direction, while in the other two directions we want to allow the membrane to bend freely.
6. Constant temperature is maintained in this SBCG simulation using Langevin dynamics, as usually done in all-atom simulations. You can take a look at these parameters under # Constant Temperature Control. In addition to the temperature control, the Langevin dynamics provides means to simulate the solvent effect implicitly. Namely, the Langevin dynamics introduces viscous drag and random forces acting on each CG bead, which can be used to mimic the viscosity of the solvent and the Brownian motion due to random hits from the molecules of the solvent. A single parameter, langevinDamping, is used to account for these effects. Here, langevinDamping is set to 2ps, as this value was found to reproduce the viscosity of water for the coarse-graining level used (150 atoms per CG bead).
7. In the last few lines of sim.conf, you can see that the simulation is designed to be minimized by 1000 steps, and then to run for 50,000,000 steps. This corresponds to 50,000,000 steps 100 fs/step = 5 s simulation time.
8. Close the text editor displaying sim.conf. Run the simulation by typing namd2 sim.conf > sim.log in a terminal window.
On a one-processor machine, this simulation will take about eight days to complete, but actually we do not need to run the full 5s. The general trend is obvious already after about the first 10ns, which can be achieved within half an hour or so. If you do not wish to run the simulation yourself, you can use the files provided in example-output/ for the following discussion on file outputs and results.
1. Open the logfile of the simulation, sim.log, with a text editor. If you did not run the simulation, use example-output/sim.log.
The logfile of a simulation contains useful information. When your simulation crashes, checking the logfile for the error message is the first step of fixing the problem. The logfile can also give you an estimate on how long a simulation will run. Find the words ``Benchmark time'' in sim.log, here you can find the speed of the simulation. Now let's examine the system via VMD.
2. Close the text editor. Open VMD, and load the psf file of the system, 6bar.psf. If you did not run the simulation, make sure you use the provided example-output/input/6bar.psf.
3. Load the output dcd file, sim.dcd. If you did not run the simulation, use example-output/output/sim.dcd. In this case, you will have 150 frames loaded in VMD, one frame for each nanosecond of the simulation.
Use VMD to take a look at the simulation result. The BAR domain system originally sits on a flat membrane, and after 150ns of simulation time, the membrane is curved into an arch. A simple SBCG simulation used here demonstrates very well how BAR domain proteins perform their job of sculpting membrane shape. For more information on BAR domain-induced membrane curvature, please see Arkhipov et al., Biophys. J., 95:2806, 2008.
4. In your VMD session, use the VMD animation tool to display the last frame, which should show a highly curved BAR domain/membrane system. Open the Graphical Representations window via Graphics Representations, and choose the Periodic tab in the lower half of the GUI window. Check on the box of +X.
This step should create an additional periodic image of the system in the +x-direction. In the VMD OpenGL window, you can see that you now have two systems separated by a large distance. This is the result of giving cellBasisVector1 a value much larger than the actual x-dimension of the system. This was done to cut off the system from its x-direction periodic neighbors, such that the system is effectively isolated in this direction, making simulation of curvature formation feasible. Same thing was done in the z-direction.
5. Uncheck the box of +X, and check +Y, and increase the value in ``Number of images" on the bottom of the Graphical Representations window to, for example, 10.
The y-direction is continuous, hence you can see in the VMD OpenGL window a long membrane patch, curved as to form a membrane tube.
6. This is the end of the tutorial. You are now ready to use SBCG!