Water Exclusion Problem of OA-G3 Host-guest Binding Complex

# Water Exclusion Problem of OA-G3 Host-guest Binding Complex ###### tags: `ResearchNote` This research note documents the simulation details about the investigation in the water exclusion problem of OA-G3 host-guest binding complex using computational methods. The goal of this project, is to figure out a way to overcome the slow kinetics of water going inside and out of the binding cavity / host cavity of OA-G3 binding complex. ## Section 1: Preparation of the input files ### Section 1-1. Files provided by SAMPL6 challenge The simulation input files (`.gro` and `.top` files) of 5 different configurations of the host-guest binding complex (from `OA-G3-0` to `OA-G3-4`) can be found in [SAMPL6](https://github.com/samplchallenges/SAMPL6.git) repository, or specifically, in the directory`SAMPL6/host_guest/SAMPLing`. Note that these files are all equilibrated, so there is no need to proceed energy mnimization, NVT equilibration and NPT equilibration. (However, carrying out these steps should still provide a reasonable starting structures. They are just not necessary.) However, since expanded ensemble and replica exchange simulation can only be performed in an NVT ensemble, in this case, we will conduct all the other advanced sampling simulations (if any) in an NVT ensemble. Therefore, here we first run an MD simulation in an NPT ensemble and extract the configuration from the output trajectory which has the volume the same as the average volume of the simulation, which can serve as the input strucutre for an NVT simulation. To start with the preparation, create the directory `Prep` and enter it, create another directory `files_from_SMPLing` and download all the relevant files provided by SAMPL6. ### Section 1-2. Input files for NPT simulations As mentioned above, the files provided by SAMPL6 challenge are all equilibrated. To prepare a reasonable strucutre for NPT simulations. Here conduct an NPT simulation for 5 ns at conditions of interest (like 1 bar and 298 K) so that the output strucutre can serve as the input strucutre for any NPT simulations (even with advanced sampling methods, if any). Accordingly, in `Prep`, create folder `MD` and execute the following command to generate a .tpr file and use it as the input file to `mdrun`: (Note: GROMACS 2020.1 was used.) ``` gmx grompp -f md.mdp -c ../files_from_SAMPLing/complex.gro -p ../files_from_SAMPLing/complex.top -o complex_md.tpr gmx mdrun -s complex_md.tpr -o complex_md.trr -c complex_md.gro -g md.log -e md.edr ``` The content of the parameter file `md.mdp` is shown below: ``` ; Run control integrator = md tinit = 0 dt = 0.002 nsteps = 2500000 ; 5 ns comm-mode = Linear nstcomm = 1 ; Output control nstlog = 1000 nstcalcenergy = 1 nstenergy = 1000 nstxout-compressed = 1000 ; Neighborsearching and short-range nonbonded interactions nstlist = 10 ns_type = grid pbc = xyz rlist = 1.0 ; Electrostatics cutoff-scheme = verlet coulombtype = PME coulomb-modifier = Potential-shift-Verlet rcoulomb-switch = 0.89 rcoulomb = 0.9 ; van der Waals vdw-type = Cut-off vdw-modifier = Potential-switch rvdw-switch = 0.85 rvdw = 0.9 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = AllEnerPres ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.10 ; EWALD/PME/PPPM parameters fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-05 ewald_geometry = 3d epsilon_surface = 0 ; Temperature coupling tcoupl = v-rescale nsttcouple = 1 tc_grps = System tau_t = 0.5 ref_t = 298 ; Pressure coupling is on for NPT pcoupl = Parrinello-Rahman pcoupltype = isotropic tau-p = 5.0 ref-p = 1.0 ; bar compressibility = 4.5e-5 ; For water at 1 atm and 300 K the compressibiilty is 4.5e-5 bar^(-1) ; refcoord_scaling should do nothing since there are no position restraints. gen_vel = yes gen-temp = 298 gen-seed = 6722267; need to randomize the seed each time. ; options for bonds constraints = h-bonds ; we only have C-H bonds here ; Type of constraint algorithm constraint-algorithm = lincs continuation = no ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 lincs-iter = 2 ``` As a result, it took about 86 minutes to finish the short MD simulation. Approximately, 1 ns of standard MD simulation takes 18 minutes. To enable easier copying of the files in the future, in `Prep`, create `Final_inputs` and execute the following command in `Final_inputs`: ``` mkdir NPT NVT && cd NPT cp ../../MD/complex_md.gro complex.gro cp ../../files_from_SAMPLing/complex.top . ``` Finally, `complex.gro` and `complex.top` can serve as the input files for any NPT simulations at 1 bar and 298 K. ### Section 1-3. Input files for NVT simulations As mentioned above, in this project, we will run all the simulations with advanced sampling methods in NVT ensembles. Note that since we don't know under what conditions were the files provided by SAMPL6 equilibrated. Therefore, instead of directly running an NVT simulation using the provided files, from the output trajectory generated by the NPT simulation we just conducted, we extract the configuration that has the volume the same as the average volume of the simulation, which is actually the volume corresponding to the target pressure (1 bar) specified in the NPT ensemble. Only by doing so, we can ensure that the system is under the correct pressure (1 bar) and the correct volume that corresponds to it. To start with, execute the following command to output the data of system volume as a function of time: ``` gmx energy -f md.edr -o volume.xvg ``` As a result, a message like below should show: ``` Last energy frame read 2500 time 5000.000 Statistics over 2500001 steps [ 0.0000 through 5000.0000 ps ], 1 data sets All statistics are over 2500001 points Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Volume 80.2288 0.016 0.444893 0.029049 (nm^3) ``` To plot the system volume as a function of time and find the configuration corresponding to the average volume, we can use the following command enabled by the installation of [MolSci_analysis](https://github.com/wehs7661/MolSci_analysis.git) to extract the configuration: ``` plot_2d -f volume.xvg -x "Time (ns)" -y "Volume ($ nm^{3} $)" -t "Volume as a function of time" -n "volume" ``` The STDOUT message will be something like: ``` Data analysis of the file: volume.xvg ===================================== Analyzing the file ... Plotting and saving figure ... Note that the first 1% of the data is truncated, which is the data that the following statistics is based on. The average of volume: 80.229 nm^{3} (RMSF: 0.005 nm^{3} max: 80.229 nm^{3} , min: 80.229 nm^{3} )The maximum occurs at 3.5 ns, while the minimum occurs at 4.024 ns. The configuration at 0.898 ns has the volume (80.228859 nm^{3} ) that is cloest to the average volume. ``` As shown above, the configuration at 0.898 ns has the volume that is cloest to the average value. To extract this configuration from the trajectory, execute (select `0` in the prompt): ``` gmx trjconv -s complex_md.tpr -f traj_comp.xtc -o complex_nvt.gro -dump 898 ``` Finally execute the command in `Final_inputs/NVT` to group up the input files: ``` cp ../../MD/complex_nvt.gro complex.gro cp ../../files_from_SAMPLing/complex.top . ``` As such, `complex.gro` and `complex.top` in `NVT` can serve as the input files for any NVT simulations. ## Section 2: System of interest ### Section 2-1. Important structural entities of the binding complex First of all, let's first visualize the strucutre to get a general sense about our system of interest. However, note that the coordinates in `.gro` files provided in SAMPL6 directory have 8 decimal places, while the standard format of a `.gro` file only allows 3 decimal places for coordinates (and 4 decimal places for velocities). Therefore, VMD will not be able to directly visualize the strucutre for the `.gro` files downloaded from SAMPL6. If you want to visualize the provided structure, just run `gmx editconf` and the precision of the coordinates will be automatically changed to 3 decimal places. (For example, use `gmx editconf -f xxx.gro`.) Alternatively, you can just visualize the structures we just stored in `Final_inputs`, which should have correct formats. Using VMD, we can visualize the important structural entities of the binding complex. The following are some summarized notes. - Note that the indices in `.gro` file and PLUMED all start from 1, while the indices in VMD start from 0. - There are 1 host molecule, 1 guest molecule, 12 sodium ions, 3 chloride ions and 2586 water molecules in `complex.gro` (in `NVT` folder). - **Host molecule**: atom 1 to atom 184 (in VMD: `index 0 to 183`) - The smaller ring (the red ring): atom 1 to 40 (in VMD: `index 0 to 39`) - The bigger ring (the green ring): atom 57 to 112 and atom 116 to 121 and atom 123 to 128 (in VMD: `index 56 to 111 or index 115 to 120 or index 122 to 127`) - In this project, we define the space between the bigger ring and the small ring as the binding cavity. - **Guest molecule**: atom 185 to atom 201 (in VMD: `index 184 to 200`) - Using atom 185 (`index 184` in VMD) and atom 187 (`index 186` in VMD), one can define the guest axis. Atom 185 and 187 are shown in yellow as below. <img src=https://i.imgur.com/hcdA1o3.jpg width=350><img src=https://i.imgur.com/Ok0rajL.jpg width=345> ### Section 2-2. The problem of water exclusion of OA-G3 Generally, host-guest binding complexes are known for the problem of slow kinetics of the water molecules going inside or out of the binding cavity. Specifically, for example, the water molecules that become trapped in the binding site and hinder the entry of the ligand can severely affect the convergence rate of a vanilla simulations by causing hysteresis effect. To tackle this problem, instead of performing a standard MD simulation, we will first try to use metadynamics to accelerate the process. In addition, rather than simply bias the distance between the host molecule and the guest molecule or the angle between the host axis and the guest axis (which has been proved to be inefficient), it would be better to directly adopt the number of water molecules in the binding cavity as the collective variable (CV). However, calculating the the number of water molecules in the binding cavity, or $N_w$, is not an easy task. One potential method we can try here is to set up two virtual atoms in metadynamics, which are the center of mass of each of the 2 rings of the host molecule (colored in green and red). Then, $N_w$ would be close (but not exactly equal) to the sum of the coordination number with water molecules of these two virtual sites ($C_w$), if the follwing two problems can be solved or avoided: - **Problem 1**: The water molecules close to the rings of the host molecule but not in the binding cavity are counted as the water molecules in the cavity. - **Problem 2**: The two virtual sites are not far enough from each other such that some of the water molecules coordinate with both the virtual atoms at the same time, which leads to overcounting of $N_w$. While it is hard to gurantee that both the assumptions will be satisfied if we adopt CV as the total coordination number of the two virtual sites, the violation of both assumptions will not influence the results if our purpose is to just accelerate the sampling isntead of getting the correct number of water molecules in the cavity. (Also, $N_w$ and $C_w$ should at least be close to each other.) As long as the CV defined above accounts for the kinetics of water molecules, it doesn't matter if the CV is exactly the number of water molecules inside the cavity. (After all, the bias still accounts for the number of water molecules in the cavity and that is enough to accelerate the sampling.) Therefore, we will just adopt this CV for now and assess the results to see if the CV should be improved. ## Section 3: Collective variable to use in metadynamics ### Section 3-1. Decoupling the interactions of the guest molecule in a vanilla MD simulation To examine if the problems of using the sum of the coordination numbers of two vitual sites will have the two problems decribe above and to get a sense about what range we should bias the CV in metadynamics if none of the problems occur, in this section, we are going to do some tests with vanilla MD simulations before we proceed to metadynamics. Specifically, what we want to do here is to run a vanilla MD simulation in which all the interactions of the guest molecule are turned off such that the guest molecule will decouple and unbind from the host molecule and the water molecules will energy the binding cavity. After that, we then need to write a PLUMED parameter file to define the collective variable as the sum of the coordination numbers of the two virtual sites ($C_{w}$)and use `plumed driver` to analyze the output trjactory of the MD simulation, calculate $C_{w}$ as a function of time, and calculate its maximum, minimum, and average, etc. To perform such an MD simulation in which all the interactions are turned off, we need to add the following parameters to the `.mdp` file: ``` ; Free energy parameters free-energy = yes calc-lambda-neighbors = -1 sc-alpha = 0.5 couple-moltype = GST couple-lambda0 = vdw-q couple-lambda1 = none couple-intramol = yes init-lambda-state = 0 ; index starts from 0 nstdhdl = 1000 dhdl-print-energy = total ; lambda-states = 1 coul-lambdas = 1.0 vdw-lambdas = 1.0 ``` Since `free-energy` is specified as `yes` and `init-lambda-state`, the system will just stay at state 0, whose `coul-lambdas` and `vdw-lambdas` are all 0, corresponding to 0 interactions. The simulation length we specified is specified as 100 ns, but actually 5 to 10 ns should be sufficient. ### Section 3-2. Calculate the number of water molecules in the binding cavity of OA-G3 using `plumed driver` #### Section 3-2-1 Method 1: Sum of the coordination numbers of two virtual sites In this section, our goal is to estimate the number of water molecules by summing up the water coordination numbers of the geomectric centers of the rings with the radius of the each ring as the cutoff distance to determine whether a water molecule is inside or outside the binding cavity. However, before we determine this parameters for setting up the collective variable, it is recommended to first take a look at some important dimensions of the binding complex as follows: <center> | $d$ |Distance (nm) | Average | Maximum | Minimum | RMSF | | :------: |-------- | :------: | :-----: | :---------: | :--: | | $d_{1}$| distance between the rings | 0.473 | 0.843 | 0.440 | 0.033 | | $d_{2}$ |bigger radius of the bigger ring| 0.857 | 0.955 | 0.789 | 0.026| | $d_{3}$ |smaller radius of the bigger ring | 0.684 | 0.738 | 0.619 | 0.025 | | $d_{4}$ | radius of the smaller ring | 0.574 |1.026 | 0.531 | 0.033 </center> <img src=https://i.imgur.com/XBabnaJ.png width=345><img src=https://i.imgur.com/TZYpvGf.png width=345> As shown above, $d_{1}$ is calculated by measuring the distance between the geometric centers of the big ring (colored in blue) and the small ring (colored in yellow). For the bigger ring, we define the larger radius as the distance between the geometric center of the whole ring and the geometric center of the white hexagon in CPK representation, while the smaller radius is defined as the distnace between the centers of the whole ring and the black hexagon in CPK representation. Lastly, the radius of the small ring is defined as the distance between the centers of the small ring and the center of the 8 atoms colored in white in the small ring. Specifically, these distances are calculated by `plumed driver` given a PLUMED parameter file as follows: ``` com1: CENTER ATOMS=1-40 # center of the smaller ring com2: CENTER ATOMS=57-112,116-121,123-128 # center of the bigger ring com_r1: CENTER ATOMS=88-93 # endpoint of the bigger radius of the bigger ring com_r2: CENTER ATOMS=57-62 # endpoint of the smaller radius of the bigger ring com_r3: CENTER ATOMS=5-11,1 # endpoint of the smaller ring d1: DISTANCE ATOMS=com1,com2 # distance between geometric centers d2: DISTANCE ATOMS=com1,com_r1 # bigger radius of the bigger ring d3: DISTANCE ATOMS=com1,com_r2 # smaller radius of the bigger ring d4: DISTANCE ATOMS=com2,com_r3 # radius of the smaller ring PRINT ARG=d1,d2,d3,d4, STRIDE=10 ``` And the command of using `plumed driver` in this case is ``` plumed driver --plumed plumed_dist.dat --mf_xtc traj_comp.xtc --pdb complex.pdb > dist.dat ``` where `complex.pdb` can be generated by ``` gmx editconf -f complex_decouple.gro -o complex.pdb ``` After the results are printed to `dist.dat`, execute the following command enabled by the installation of [MolSci_analysis](https://github.com/wehs7661/MolSci_analysis.git): ``` plumed_driver -i dist.dat -x "Time (ns)" -y "Distance (nm)" -t "Distance as a function of time" -n "distance" -ts 2 ``` As shown in the table above, apparently, using the radius of the each ring as the cutoff distance will lead to **Problem 2**, not to mention that **Problem 1** is fundamentally hard to avoid. Therefore, this collective variable might not be the best choice to estimate the number of water molecules in the cavity. To validate this, one could use `plumed driver` to analyze the trajectory with the parameter file below and one should find that the estimated number of water molecules does not match the number observd from the visualization of the structure. Therefore, we have to resort to other approach, such as **Method 2** introduced in the next section. ``` com1: COM ATOMS=1-40 # COM of the smaller ring com2: COM ATOMS=57-112,116-121,123-128 # COM of the bigger ring water_group: GROUP ATOMS=217-7974:3 # oxygen atom of the water molecules x1: COORDINATION GROUPA=com1 GROUPB=water_group R_0=0.3 x2: COORDINATION GROUPA=com2 GROUPB=water_group R_0=0.5 y: MATHEVAL ARG=x1,x2 VAR=a,b FUNC=a+b PERIODIC=NO PRINT ARG=x2, STRIDE=10 ``` #### 3-2-2. Method 2: Apprixmate the binding cavity with a sphere In this section, we try to approximate the binding cavity with a sphere and count the number of water molecules in the sphere to provide an estiamte. To start with, we first use VMD to visualize the binding complex and use the following simple script in Tk Console of VMD to determine a reasonable radius and center of the sphere. ```Tcl # Step 1: set up the geometric centers of the rings set small_ring [atomselect 0 "index 0 to 39"] set com1 [measure center $small_ring] set big_ring [atomselect 0 "index 56 to 111 or index 115 to 120 or index 122 to 127"] set com2 [measure center $big_ring] # Step 2: Write a procedure (proc) - inside_point proc inside_point {coef1 lst1 coef2 lst2 expr} { foreach item1 $lst1 item2 $lst2 { lappend r [expr [string map {%x $item1 %y $item2 %a $coef1 %b $coef2} $expr]] } return $r } # Step 3: Determine the radius and the center of the sphere inside_point 0.4 $com1 0.6 $com2 {%a * %x + %b * %y} draw color cyan draw material Transparent draw sphere $com radius 5 resolution 24 ``` Here are some comments about the Tcl script above: - In Step 1, we set up the geometric centers of the rings. - To set a variable, use `set var_name [var_value]`. - To print out the value of a variable, use `puts var_name`. - `atomselect` and `measure center` are the functions/procedures in VMD, which are for selecting and determining the geometric centers of a group of atoms, respectively. The atom selection language in Tk Console is the same as the one we used in VMD representation. - In Step 2, we write a function/procedure(`proc`) `inside_point` to calculate the coordinates of an inside point given the end points and their weights - The data type of variables such as `com1` and `com2` are lists. Not like many other language, to perform a element-wise operation like, we cannot code something like `expr com = 0.4 * com1 + 0.6 * com2`, but have to write a function to loop over the whole list and append the calculated elements to a new list. - The first line of the procedure `proc inside_point {coef1 lst1 coef2 lst2 expr}` indicates the procedure name (`inside_point`) and input arguments in order (`coef1`, `lst1`, `coef2`, `lst2`, and `expr`). - `foreach item1 $lst1 item2 $lst2` loops over every element (`item1` and `item2`) in `lst1` and `lst2`. An analog of `foreach item $lst1` in Python is `for i in lst1:`. - `lappend` appends list element onto a variable. Its synopsis is `lappend var_name val1 val2 ...`. An example is shown as below: ``` set x 1 # x = 1 lappend var 2 # now x = 1 2 lappend var 3 4 5 # now x = 1 2 3 4 5 ``` - To evaluate a mathematical expression, use `expr`. - `string map` replaces characters in string based on the key-value pairs in `charMap`. `charMap` is a list of key value key value ... as in the form returned by `array get`. Each instance of a key in the string will be replaced with its corresponding value. Both key and value may be multiple characters. Replacement is done in an ordered manner, so the key appearing first in the list will be checked first, and so on. `string` is only iterated over once, so earlier key replacements will have no effect for later key matches. For example,`string map {abc 1 ab 2 a 3 1 0} "1abcaababcabababc"` will return the string `01321221`. - Therefore, given the input argument `$expr`, `[expr [string map {%x $item1 %y $item2 %a $coef1 %b $coef2} $expr]]` basically substitutes `%x` with `item1` in `lst1`, `%y` with `item2` in `lst2` and so on. Then, `string map` returns a string of expression which will be performed by `expr`. The whole bracket here will be appended to the varlue `r`, which will be returned by the function in the end. - In Step 3, we determine a reasonable radius and center of the sphere by visualizing the sphere. - By trying different coefficients used in `inside_point 0.4 $com1 0.6 $com2 {%a * %x + %b * %y}`, we decide the center of the sphere. - To specify the settings of the drawing method, use `draw color` and `draw material`. - `sphere $com radius 5`: After severl trials, have the coefficients as 0.4 and 0.6 and a radius as 5 should be a good choice. - As shown below, with this center and radius, the sphere should approximate the binding cavity reasonably well. Specificall, the sphere does not include any space outside the binding cavity and the space in the cavity that was not included by the sphere seems not enough for any water molecules to be adopted. Therefore, we will adopt this set of parameters in the PLUMED parameter file. <center><img src=https://i.imgur.com/bmBO5Bp.jpg width=600></center> </br> To analyze the trajectory, we use the following PLUMED file adopting the radius and the center of the sphere decided above: ``` center1: CENTER ATOMS=1-40 # geometric center of the smaller ring center2: CENTER ATOMS=57-112,116-121,123-128 # geometric center of the bigger ring center: CENTER ATOMS=center1,center2 WEIGHTS=2,3 # center of the sphere water_group: GROUP ATOMS=217-7974:3 # oxygen atom of the water molecules n: COORDINATION GROUPA=center GROUPB=water_group R_0=0.5 # unit of length: nm PRINT ARG=n, STRIDE=10 ``` Similarly, here some comments about the first half of the code above: - Using `center` method, we assign weights to `center1` and `center2` to set up the inside point such that the ratio of its distance to `center1` and `center2` is 2 to 3. - `water_group: GROUP ATOMS=217-7974:3`: Instead of looping over all the atoms of water molecules, we can only consider the oxygen atoms of the water molecules to decrease the computational cost. Note that `COORDINATION` in the code above calculates the number of contacts between two groups of atoms and is defined as $\sum_{i \in A} \sum_{j \in B} s_{ij}$ where $s_{ij}$ is 1 if the contact between atoms $i$ and $j$ is formed, zero otherwise. In actuality, $s_{ij}$ is replaced with a switching function so as to ensure that the calculated CV has continuous derivatives. The default switching function is: $$s(r) = \frac{(\frac{r-d_{0}}{r_{0}})^{n}}{(\frac{r-d_{0}}{r_{0}})^{m}}$$ where $r_0$ is the cutoff distance, $d_0$ is a parameter that we can use to shift the center of the sphere, and the exponent $n$ and $m$ can be used to adjust the shape of the switching function (The default value of $n$ and $m$ are 6 and 12, respectively.) For more information about this method in PLUMED, please refer to the [documentation of PLUMED](https://www.plumed.org/doc-master/user-doc/html/_c_o_o_r_d_i_n_a_t_i_o_n.html). Accordingly to the equation above, apparently, $s(r)$ approaches to $n/m$ as $r$ approach to $r_0$. With the default values of $n$ ($n=6$) and $m$ ($m=12$), $s(r)=0.5$ at $r_0$, which means that the water molecules near the entry of the binding cavity (but not in the cavity) will also be “partially counted”. Therefore, to find a reasonable set of parameters $n$ and $m$, we should try different values of $m$ (while fixing the val of $n$) to analyze the trajectory as below. (Note that $n=1$ and $d_0=0$ for all the cases tabulated below.) | Value of $m$ | Average | Maximum | Minimum | $s(r)$ at $t=0$ | $s(r)$ at $1.01 r_{0}$ | | :------: | :------: | :------: | :-: | :-: | :-: | | 12 | 2.703 | 3.945 | 0.148 | 0.2367 | 0.0788 | | 25 | 2.463 | 3.678 | 0.002 | 0.05535 | 0.0354 | | 30 | 2.451 | 3.664 | 0.000 | 0.04824 | 0.0287 | | 50 | 2.433 | 3.648 | 0.000 | 0.03599 | 0.0155 | | 100| 2.426 | 3.645 | 0.000 | 0.02860 | 0.0059 | | 200| 2.425 | 3.645 | 0.000 | 0.02681 | 0.0016 | Based on the results above, adopting $m=30$ (and $n=1$, $d_{0}=0$) should be a reasonable choice based on the following observations: - The number of the water molecules in the sphere is nearly 0 (0.04824) at $t=0$. - The minimum of the number of water molecules in the sphere is 0, which occurs at 33.29 ns. Visualizing the structure at this time frame, what is observed is that there are two or three water molecules near the entry. There was only one molecule in the cavity, but it was relatively close to the entry of the cavity. Therefore, the method might not be perfect, but at least provide a reasonable estimate. - The maximum occurs at 26.25 ns and the value also matches the visualization of the structure. (There are 3 to 4 water molecules in the cavity.) - Note that with a such high exponent ($m=30$), we should keep an eye on whether the simulation diverges due to numerical instability. However, if the simulation does not diverge, such a $m$ value should be fine. As a supplementary note, here are some points about how the observation above are made and another way to calculate exactly the same CV: - The data tabulated above can be obtained by using `plumed driver` and `plumed_driver`. - Note that the molecule might cross the periodic boundary and looks broken. To fix it when visualizing the structure at 33.29 ns, use `gmx trjconv` as below: ``` gmx trjconv -s complex_decouple.tpr -f traj_comp.xtc -o complex_nojump.xtc -center -pbc nojump gmx trjconv -s complex_decouple.tpr -f complex_nojump.xtc -o complex_center_min.gro -center -pbc mol -ur compact -dump 33290 ``` Also, use the following command in Tk Console of VMD as needed: ``` pbc wrap -centersel "index 0 to 200" -center com -compound residue -all ``` - In addition to using `COORDINATION` as above, we can also use `DENSITY` and `INSPHERE` to calcuate exactly the same thing. Specifically, with the code below, we can compare the number of water molecules in the binding cavity as a function of time obtained by these two approaches in the figure below. - The content of the code: ``` center1: CENTER ATOMS=1-40 # geometric center of the smaller ring center2: CENTER ATOMS=57-112,116-121,123-128 # geometric center of the bigger ring center: CENTER ATOMS=center1,center2 WEIGHTS=2,3 # center of the sphere water_group: GROUP ATOMS=217-7974:3 # oxygen atom of the water molecules n1: COORDINATION GROUPA=center GROUPB=water_group R_0=0.5 NN=1 MM=30 # unit of length: nm d: DENSITY SPECIES=217-7974:3 n2: INSPHERE ATOM=center DATA=d RADIUS={RATIONAL R_0=0.5 NN=1 MM=30} # number of water molecules in the sphere PRINT ARG=n1,n2 STRIDE=10 ``` - Output message of `plumed_driver`: ``` Data analysis of the file compare.dat: ====================================== The average of n1: 2.451 (RMSF: 0.194 max: 3.664 , min: 0.000 ) The average of n2: 2.451 (RMSF: 0.194 max: 3.664 , min: 0.001 ) ``` - Taking the data of $n_{1}$ as the reference, the RMSD between $n_1$ and $n_2$ is only about $5.36 \times 10^{-5}$. <center><img src=https://i.imgur.com/MYn7Pw3.png width=500></center> Based on the data analysis above, in the metadynamics below, we will use `COORDINATION` to calculate $n_{1}$ and use it as the CV. ## Section 4. Metadynamics accelerating the water exclusion/insertion process ### 4-1. 1-dimensional metadynamics #### 1. Set up and run an 1D metadynamics As mentioned above, we decided to use `COORDINATION` to estimate the number of water molecules in the binding cavity and bias it in our 1-dimensional metadynamics. Specifically, we use the following PLUMED parameter file (`plumed.dat`): ``` center1: CENTER ATOMS=1-40 # geometric center of the smaller ring center2: CENTER ATOMS=57-112,116-121,123-128 # geometric center of the bigger ring center: CENTER ATOMS=center1,center2 WEIGHTS=2,3 # center of the sphere water_group: GROUP ATOMS=217-7974:3 # oxygen atom of the water molecules n: COORDINATION GROUPA=center GROUPB=water_group R_0=0.5 NN=1 MM=25 # unit of length: nm METAD ... ARG=n SIGMA=0.02 HEIGHT=0.3 PACE=500 BIASFACTOR=10 TEMP=298 LABEL=metad GRID_MIN=0.0 GRID_MAX=6.0 GRID_BIN=600 FILE=HILLS_COORDN ... METAD PRINT STRIDE=10 ARG=n,metad.bias FILE=COLVAR ``` With the parameter file above, we can run a metadynamics using PLUMED by simply adding `-plumed` flag to the `gmx mdrun` command if PLUMED has been successfully patched with GROMACS: ``` gmx grompp -f complex.mdp -c complex.gro -p complex.top -o complex.tpr gmx mdrun -deffnm complex -plumed plumed.dat ``` Note that here we run the metadynamics in an NVT ensemble, so the only differences between `complex.mdp` and `md.mdp` used in Section 1.2 are: - `nsteps = 250000000` (500 ns) instead of `nsteps = 2500000` (5 ns) - No pressure coupling in `complex.mdp` #### 2. Data analysis of the 1D metadynamics Although we set `nsteps = 250000000` in `complex.mdp`, we don't have to complete the whole simulation of 500 ns. Here, we just want to see how the CV varies in a long enough simulation, so we stopped the simulation at about 104 ns and there were two files output by PLUMED (with the name specified in `plumed.dat`), including `COLVAR` and `HILLS_COORDN`. - The CV as a fucntion of time To examine the number of water molecules in the binding cavity as a function of time, execute: ``` plot_2d -f COLVAR -x 'Time (ns)' -y 'Number of water molecules in the cavity' -t "Number of water molecules in the cavity as a function of time" -n "n_waters" -cx "ps to ns" ``` As a result, the following STDOUT is shown: ``` Data analysis of the file: COLVAR ================================= Analyzing the file ... Plotting and saving figure ... The average of number of water molecules in the cavity: 1.581 (RMSF: 0.746 max: 4.637, min: 0.000) The maximum occurs at 60.5654 ns, while the minimum occurs at 30.2923 ns. The configuration at 2.28124 ns has the number of water molecules in the cavity (1.580708) that is cloest to the average volume. ``` In addition, a figure `n_waters.png` will be generated, as shown below. According to the figure shown below, it took about 25 ns for the water molecules to be inserted in or excluded from the binding cavity. <center><img src=https://i.imgur.com/D3gJY7d.png width=500></center> - Visuaization of the trajectory To see what really happened during the simulation, we can visualize the trajectory, which often involves the following two steps to preprocess the trajectory: - **Step 1**: Center on the binding complex and remove jumps by `-pbc nojump`: ``` gmx trjconv -s complex.tpr -f complex.xtc -o complex_nojump.xtc -center -pbc nojump ``` - **Step 2**: Center on the binding complex and remove pbc by ``-pbc mol -ur compact`: ``` gmx trjconv -s complex.tpr -f complex_nojump.xtc -o complex_center.pdb -center -pbc mol -ur compact -dt 250 ``` - **Step 3**: If the molecule still seems broken, try the following command in the Tk Console in VMD: ``` pbc wrap -centersel "index 0 to 200" -center com -compound residue -all ``` Using VMD to visualize `complex_center.pdb`, we can see that (like if the water molecules were indeede inserted to and excluded from the bindin cavity or if the guest molecule bind to and unbind from the host molecule several times) - The free energy surface To get the free energy surface, execute: ``` plumed sum_hills --hills HILLS_COORDN ``` then `fes.dat` will be output by PLUMED. To visualize, execute: ``` plot_2d -f fes.dat -x "Number of water molecules in the binding cavity" -y "Free energy (kcal/mol)" -t "Free energy surface of OA-G3 binding" -n "OAG3_fes" -cy "kJ/mol to kcal/mol" ``` ``` plumed sum_hills --hills HILLS_COORDN --stride 20000 ``` ![](https://i.imgur.com/Vm0lI6c.png) ``` plot_2d -f "fes_0.dat" "fes_1.dat" "fes_2.dat" "fes_3.dat" "fes_4.dat" -l "20 ns" "40 ns" "60 ns" "80 ns" "100ns" -x "Number of water molecules in the binding cavity" -y "Free energy surface (kcal/mol)" -cy "kJ/mol to kcal/mol" -t "Free energy surface of OA-G3 binding" -n "OAG3_fes" ``` ![](https://i.imgur.com/DuFjaWM.png)

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.