The setup of the test includes the following paraemters:
<divclass="horizontal-scroll">
<tableclass="valuetable">
<thead>
<tr>
<th>Li</th>
<th>Lf</th>
<th>index1</th>
<th>index2</th>
<th>num_win</th>
<th>fc_pull</th>
<th>v_pulling</th>
<th>total_steps</th>
<th>increment_steps</th>
<th>wTotal</th>
<th>wDelta</th>
</tr>
</thead>
<tbody>
<tr>
<td>3.0 nm</td>
<td>3.2 nm</td>
<td>14 (Hydrogen UNL 1)</td>
<td>5439 (Hydrogen SER 349)</td>
<td>100</td>
<td>1000.0 kj/mol/nm^2</td>
<td>0.02 nm/pc</td>
<td>30,000</td>
<td>10</td>
<td>100,000</td>
<td>10</td>
</tr>
</tbody>
</table>
<p>
<p><strong>Generally Pick the strongest parts of the Complex</strong></p>
<spanclass="hljs-number"><strong>1.</strong></span><spanclass="hljs-keyword">For</span> proteins: C-alpha atoms are picked because they<spanclass="hljs-comment">'re usually the backbone atoms that represent the structure very well, so stretching that wont break the system entirely (they are strong enough to stretch).</span><br>
<spanclass="hljs-number"><strong>2.</strong></span><spanclass="hljs-keyword">For</span> nucleic acids: Phosphate atoms <spanclass="hljs-keyword">or</span> specific base atoms (because again, they<spanclass="hljs-comment">'re the strongest backbones).</span><br>
<spanclass="hljs-number"><strong>3.</strong></span><spanclass="hljs-keyword">For</span> small molecules <spanclass="hljs-keyword">or</span> ligands: <spanclass="hljs-keyword">Select</span> atoms that are involved <spanclass="hljs-keyword">in</span><spanclass="hljs-keyword">key</span> interactions <spanclass="hljs-keyword">or</span> are structurally significant.
If you are trying to get the binding affinity of a docked complex (which is what I am trying to do in this case,) I would select the center of mass of the ligand, and the key atom in a resiude in the binding pocket of the receptor. <br><br>
<p>
As stored in the <pre><code>./gh_scoresolver/tests/pdb/ahl_dock_luxr_1.pdb</code></pre> PDB file, we have a ligand protein ahl docked into a receptor protein luxr.
</p>
</p>
</div>
<br><br>
The window 0 pdb file looks as follows with the two atoms that were pulled indicated:
strains. Current Protocols, 1(1). https://doi.org/10.1002/cpz1.20 </li>
</ul>
</div>
<divid="Umbrella-Sampling">
<h2>Umbrella Sampling</h2>
<h3>Overall Purpose</h3>
<p>The Umbrella Sampling technique used in molecular dynamics is especially useful in this project as it manually enhances the sampling of complex energy landscapes of protein data banks, involving free energy calculations and calculations of other thermodynamic variables. It can overcome large energy barries and simulate rare events that other rudimentary simulations would not be able to achieve.</p>
<p></p>
<hr>
<h3>Design</h3>
<p>
The program was designed as a pipeline between various molecular dynamics tools. However, the most important part of the pipeline is the Umbrella Sampling process indicated by the caption on the arrow between the windows rectangle and metafile/cv files rectangle.
Below is the pseudocode outlining the program above including the SMD pulling loop, Umbrella Sampling, wham analysis, and binding coefficient algorithm. This paraphrases the code itno simple words, and the full code in python is in the software gitlab repository.
CREATE directory structure for new test (cv, windows, hist)
PRINT "New folder created"
# Initialize System
LOAD PDB file from pdb_path
CREATE Modeller to check and add missing hydrogens
APPLY forcefield to the system
# Define System and Simulation
SET system with forcefield, constraints, and parameters
INITIALIZE Langevin integrator for temperature and timestep
SET positions and velocities for the system
MINIMIZE system energy
# Parameters for Steered Molecular Dynamics (SMD)
DEFINE starting and ending values for collective variable (CV)
DEFINE pulling force constant, velocity, and number of windows
DEFINE steps and recording intervals for simulation
# Add pulling force to system
CREATE pullingForce as harmonic force on CV
ADD pullingForce to the system
INITIALIZE simulation with pulling force
# Perform Steered Molecular Dynamics (SMD) Loop
FOR each window
PERFORM simulation steps
UPDATE pulling force parameter r0
RECORD collective variable (CV) values
IF CV reaches target value for the current window
SAVE system coordinates
END IF
UPDATE progress bar
END FOR
# Save Windows and CV Data
FOR each saved window
SAVE window coordinates as PDB file
END FOR
# Run Simulations on Each Window
DEFINE function run_window(window_index)
LOAD coordinates from corresponding window
SET system positions to window coordinates
RESET temperature and run simulation for steps
RECORD collective variable (CV) values for each window
SAVE CV values to file
END FUNCTION
CALL run_window() for each window
# Generate Histogram and Save
FOR each window
LOAD CV data
GENERATE histogram of CV values
SAVE histogram plot
CREATE metafile entry for WHAM analysis
END FOR
# WHAM Analysis
PREPARE WHAM input arguments
RUN WHAM on metafile to compute PMF
SAVE WHAM output and log
# Post-WHAM Processing
LOAD PMF data from WHAM output
PARSE coordinates, free energies, and probabilities
PLOT free energy and probability graphs
SAVE graph as PNG file
# Compute Equilibrium Constant (Kd)
COMPUTE delta_G_bound and delta_G_unbound from PMF
CALCULATE standard free energy difference (ΔG)
COMPUTE dissociation constant (Kd) from ΔG
OUTPUT Kd value
END
</code>
</pre>
</div>
</p>
<h3>Test & Results</h3>
<p>
This program was tested with the PDB that had protonated PFOA in FAB-GFP. Although deprotonated PFOA could be simulated in antechamber, we were unable to generated a forcefield file for the program to use to simulate deprotonated PFOA. Thus we had to use protonated pfoa where the forcefield was able to generate with a program called LigParGen, named and stored as: <pre><code>umbrellasampling/tests/forcefields/pfoa_protonated.xml</code></pre><br>
<spanclass="hljs-number"><strong>1.</strong></span><spanclass="hljs-keyword"></span>The protonated PFOA in the PDB was made by manually adding a hydrogen atom to PFOA through ChimeraX and edited into the FAB-GFP again. <br>
<spanclass="hljs-number"><strong>2.</strong></span><spanclass="hljs-keyword"></span>To prevent NaN coordinates, the simulation's energy was minimzied by moving the protiens and PFOA around. The conformal differences of: <br>
<spanclass="hljs-number"><strong>UNL 1</strong></span> Refers to the protonated PFOA, and the conformational changes of the protein fab-gfp in window 0 that follow from the minimized energy process (via openMM) is in red. The added hydrogen according to ChimeraX is in solid green.<br>
<spanclass="hljs-number"><strong>UNL 377</strong></span> Refers to the deprotonated PFOA, and the conformational changes of the protein fab-gfp that follow from the minimized energy process (via antechamber) is in green with reduced opacity.
The Kd and delta G calculated from the PMF graph is:
<divclass="horizontal-scroll">
<tableclass="value-table">
<thead>
<tr>
<th>Dissociation Constant</th>
<th>ΔG°</th>
</tr>
</thead>
<tbody>
<tr>
<td>2.26 M</td>
<td>0.48224 kcal/mol</td>
</tr>
</tbody>
</table>
</div>
</p>
<h3>Learn</h3>
<p> Although we were unable to simulate deprotonated PFOA with the program, the major takeaways with the test was the stark differences in the behaviors of protonated and deprotonated PFOA. <br>
<spanclass="hljs-number"><strong>1.</strong></span><spanclass="hljs-keyword">When the energy was minimized, the PFOA was put outside the binding pocket of FAB-GFP, thus the very high dissociation constant, or low binding strength, was explained by this result. <br>
<spanclass="hljs-number"><strong>2.</strong></span><spanclass="hljs-keyword">The deprotonated PFOA had a much lower dissociation constant, thus much higher bonding strength, explained by being inside the binding pocket of FAB-GFP. <br>
</p>
</div>
<divid="References">
<h2>References</h2>
<ul>
<li>Blinov, M. L., J. C. Schaff, D. Vasilescu, Moraru, II, J. E. Bloom, and L. M. Loew. 2017.
Compartmental and Spatial Rule-Based Modeling with Virtual Cell. Biophysical journal
113:1365-1372. PMC5627391
</li>
<li>Cai, L., Friedman, N., & Xie, X. S. (2006). Stochastic protein expression in individual cells at the single molecule level. Nature, 440(7082), 358–362. https://doi.org/10.1038/nature04599
</li>
<li>Cowan, A. E., Moraru, II, J. C. Schaff, B. M. Slepchenko, and L. M. Loew. 2012. Spatial modeling of cell signaling networks. Methods Cell Biol 110:195-221. PMC3288182
</li>
<li>Schaff, J., C. C. Fink, B. Slepchenko, J. H. Carson, and L. M. Loew. 1997. A general computational framework for modeling cellular structure and function. Biophysical journal 73:1135-1146. PMC1181013
</li>
<li>Ali Azam T, Iwata A, Nishimura A, Ueda S, Ishihama A. Growth phase-dependent variation in protein composition of the Escherichia coli nucleoid. J Bacteriol. 1999 Oct;181(20):6361-70. doi: 10.1128/JB.181.20.6361-6370.1999. PMID: 10515926; PMCID: PMC103771.</li>
<li>Tuttle, A. R., Trahan, N. D., & Son, M. S. (2021). Growth and maintenance of escherichia coli laboratory strains. Current Protocols, 1(1). https://doi.org/10.1002/cpz1.20 </li>
