Skip to content
Snippets Groups Projects

Edwardimport

Merged Edwardimport
edwardimport into main
Merged Edward Kim requested to merge edwardimport into main
Viewing commit 4e7864d4
Prev Next
Show latest version
1 file
+ 19
− 282
Compare changes
  • Side-by-side
  • Inline
wiki/pages/drylab.html
+ 19
− 282
@@ -5,8 +5,6 @@
{% block page_content %}
<script type="text/javascript" async="" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"></script>
<head>
</head>
@@ -662,24 +660,7 @@
<div id="program" class="mt-4">
<h2>Program Flow</h2>
<pre class="mermaid">
---
title: Bind Coefficient Calculation
---
graph LR;
pdb[("test.pdb")]
win["windows"]
mt["metafile.txt & cv data"]
pmf["pmf.txt"]
traj["smd_traj.dcd"]
wl["wham_logs.txt"]
kd["$$Kd$$"]
pdb -- "SMD Pulling" --> traj & win
win -- "Umbrella Sampling" --> mt
mt -- "WHAM" --> pmf & wl
pmf -- "Binding Coeff. Algo" --> kd
win -- "Binding Coeff. Algo" --> kd
</pre>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/mermaid-diagram-2024-09-28-193939.png"></img>
</div>
<div id="troubleshooting" class="mt-4">
@@ -693,139 +674,17 @@
<div id="definitions">
<h3>Definitions</h3>
<p>Important terms used:</p>
<ul>
<li><span class="math inline"><em>F</em></span>: Free energy</li>
<li><span class="math inline"><em>k</em><sub><em>B</em></sub></span>:
Boltzmann constant</li>
<li><span class="math inline"><em>T</em></span>: Temperature</li>
<li><span class="math inline"><em>Z</em></span>: Canonical partition
function</li>
<li><span class="math inline"><em>ξ</em></span>: Collective Variable
(CV) <span class="math inline">→</span> represents the progress of a
reaction into a single value of a <strong>state</strong></li>
<li><span class="math inline"><strong>r</strong></span>: Atomic
coordinates</li>
<li><span class="math inline"><em>P</em>(<em>ξ</em>)</span>: Probability
distribution of the system as a function of CV</li>
<li><span class="math inline"><em>β</em></span>: Inverse temperature ($=
1 / k_B T $)</li>
<li><span class="math inline"><em>U</em>(<strong>r</strong>)</span>:
Potential energy of the system</li>
<li><span
class="math inline"><em>U</em><sub><em>i</em></sub><sup>bias</sup></span>:
Biasing potential for window $i $</li>
<li><span class="math inline"><em>k</em></span>: Force constant for the
biasing potential</li>
<li><span class="math inline"><em>ξ</em><sub><em>i</em></sub></span>:
Target value of the CV in window $i $</li>
<li><span
class="math inline"><em>P</em><sub><em>i</em></sub><sup>bias</sup>(<em>ξ</em>)</span>:
Biased probability distribution for window $i $</li>
<li><span
class="math inline"><em>n</em><sub><em>i</em></sub>(<em>ξ</em>)</span>:
Number of samples in bin <span class="math inline"><em>ξ</em></span> for
window $i $</li>
<li><span
class="math inline"><em>F</em><sub><em>i</em></sub><sup>bias</sup>(<em>ξ</em>)</span>:
Biased free energy in window $i $</li>
<li><span class="math inline"><em>F</em>(<em>ξ</em>)</span>: Unbiased
free energy (Potential of Mean Force, PMF)</li>
<li><span
class="math inline"><em>δ</em><sub><em>a</em></sub>(<em>x</em>)</span>:
Dirac delta function</li>
</ul>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/definitions.png"></img>
</div>
<div id="equations">
<h3>Equations</h3>
<ol type="1">
<li>Reaction Coordinate</li>
</ol>
<p><span class="math display"><em>ξ</em>(<strong>r</strong>) = the
progress of a reaction as function of atomic
configuration/coordinates</span></p>
<p>Note: different atomic configurations <span
class="math inline"><strong>r</strong></span> can have the same reaction
coordinate due to <span class="math inline"><em>ξ</em></span>’s limited
scope, especially in reversible reactions/processes. However, this is a
necessary sacrifice because we need to reduce complexity and dimension
of the problem. Essentially, <span class="math inline"><em>ξ</em></span>
is a simplification of <span
class="math inline"><strong>r</strong></span>.</p>
<ol start="2" type="1">
<li>Free energy in the canonical ensemble:</li>
</ol>
<p><span
class="math display"><em>F</em> = −<em>k</em><sub><em>B</em></sub><em>T</em>ln <em>Z</em></span></p>
<p>In context, <span class="math inline"><em>F</em></span> can be
expressed in terms of a specific state <span
class="math inline"><em>ξ</em></span>, so <span
class="math inline"><em>F</em>(<em>ξ</em>) = −<em>k</em><sub><em>B</em></sub><em>T</em>ln <em>P</em>(<em>ξ</em>) + <em>C</em></span></p>
<ol start="3" type="1">
<li>Dirac Delta Function: ( <span
class="math inline"><em>δ</em>(<em>x</em>) = <em>δ</em><sub>1</sub>(<em>x</em>)</span>
):</li>
</ol>
<p><span class="math display">$$
{\displaystyle \delta _{a}(x)={\frac {1}{\left|a\right|{\sqrt {\pi
}}}}e^{-(x/a)^{2}}}
$$</span></p>
<p>In the case of this theory, it projects the space of <span
class="math inline"><strong>r</strong></span> to the space of <span
class="math inline"><em>ξ</em></span>. It picks out the relevant values
of <span class="math inline"><em>ξ</em>(<strong>r</strong>)</span> close
to the value <span class="math inline"><em>ξ</em></span> in question and
concentrates/amplifies those specific values of <span
class="math inline"><strong>r</strong></span>.</p>
<p><span
class="math display"><em>P</em>(<em>ξ</em>) = ∫<em>P</em>(<strong>r</strong>)<em>δ</em>(<em>ξ</em> − <em>ξ</em>(<strong>r</strong>))<em>d</em><strong>r</strong></span></p>
<ol start="4" type="1">
<li>Partition Function (discrete):</li>
</ol>
<p><span class="math display"><em>Z</em> = ∑energies of i<sup>th</sup>
state weighed by boltzmann distribution
factor = ∑<sub><em>i</em></sub><em>e</em><sup>−<em>β</em><em>E</em><sub><em>i</em></sub></sup></span></p>
<p>Encapsulates the total distribution of energy, which can later be
used to normalize probability distributions over boltzmann factors.</p>
<ol start="5" type="1">
<li>Probability Distribution in <span
class="math inline"><strong>r</strong></span></li>
</ol>
<p><span class="math display">$$
P(\textbf{r}) = \frac{\exp(-\beta E(r))}{Z} : \beta = \frac{1}{k_B T}
$$</span></p>
<p>Here higher values of <span class="math inline"><em>β</em></span>, or
coldness, makes states with lower temperatures more probable than
otherwise. While lower values of <span
class="math inline"><em>β</em></span>, leading to higher temperatures,
makes states with higher temperatures more probable than otherwise.
Regardles of positive <span
class="math inline"><em>β</em> &gt; 0</span>, lower temperature states
are more probable than higher temperature states. Note that <span
class="math inline"><em>λ</em> ≡ <em>β</em></span>.</p>
<ol start="6" type="1">
<li>Probability Distribution (continous):</li>
</ol>
<p><span class="math display">$$
P(\xi) = \int P(\textbf{r})\delta(\xi−\xi(\textbf{r}))d\textbf{r} =
\frac{1}{Z} \int e^{-\beta
E(\textbf{r})}\delta(\xi−\xi(\textbf{r}))d\textbf{r}
$$</span></p>
<p>Which is the integration over all possible degrees of freedom of a
reaction coordinate <span class="math inline"><em>ξ</em></span></p>
<p>In all, its basically multiplying all the possible factors that
affect the probability of having state <span
class="math inline"><em>ξ</em></span> and integrating them because
again, probability is a statistical quantity that is approximated.</p>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/equations1-6.png"></img>
</div>
<div id="Adding Biasing">
<h3>Adding Biasing</h3>
<div class = "LaTeX">
<p><span class="math display">
$$
w_i(\xi) = \frac12 k(\xi - \xi_i)^2
$$
</span></p>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/equations7.png"></img>
</div>
<p>
It acts a spring force with strength $k $ that restrains the system
@@ -852,152 +711,30 @@
with biasing potentials, we subtract the biased results by the biasing
potential and retrive the true energy curve. </strong>
</p>
<ol start="8" type="1">
<li>
Energy at window <span class="math inline"><em>i</em></span> to unbiased version:
</li>
</ol>
<div class="LaTeX">
<p><span class="math display">
$$
E_i^\text{bias}(\textbf{r}) = E_i(\textbf{r}) + w_i((\xi(\textbf{r}))) \quad [8]
$$
</span></p>
</div>
<ol start="9" type="1">
<li>
Probability distribution at window <span class="math inline"><em>i</em></span> to unbiased version:
</li>
</ol>
<div class="LaTeX">
<p><span class="math display">
$$
P_i^\text{bias}(\xi) \propto \int e^{-\beta (E_i^\text{bias}(\textbf{r}))} \delta (\xi - \xi(\textbf{r})) dr \propto P_i(\xi) e^{-\beta w_i(\xi)}, \quad [9]
$$
</span></p>
</div>
<ol start="10" type="1">
<li>
Free Energy at window <span class="math inline"><em>i</em></span> to unbiased version:
</li>
</ol>
<div class="LaTeX">
<p><span class="math display">
$$
F^\text{bias}(\xi) = -k_B T \ln P_i(\xi) + w_i(\xi) + C_i = F_i(\xi) + w_i(\xi) + C_i. \quad [10]
$$
</span></p>
<p>Combining equations <span class="math inline">[8], [9], [10]</span>,</p>
<p><span class="math display">
$$
F_i(\xi) = F_i^\text{bias}(\xi)-w_i(\xi) + C_i = -k_B T \text{log}(P_i^\text{bias}(\xi)) - w_i + C_i.
$$
</span></p>
</div>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/equations8-10.png"></img>
</div>
<div id="Calcualting Binding Energy">
<p>Calculating Binding Constant of a Ligand-Receptor System from the Free Energy Curve</p>
<p>Theoretically, <p><span class="math display">$$K_d = \frac{[L][R]}{[LR]}$$</span></p> where <span
class="math inline">[<em>L</em>], [<em>R</em>], [<em>L</em><em>R</em>]</span>
are the cocentrations of the ligand, receptor, and ligand-receptor
complexes respectively. However, since this is computationally
intractible to calculate, a more convenient equation would be a
relationship between binding free energy and the <span
class="math inline"><em>K</em><sub><em>d</em></sub></span> constant.</p>
<ol type="1"><li>Absolute Binding Free Energy</li></ol>
<p><span class="math display">
$$
\Delta F^\circ = k_B T\ln K_\mathrm{d} \Rightarrow K_\mathrm{d} =
e^{\frac{\Delta F^\circ}{k_B T}}
$$
</span></p>
<p>Where <span class="math inline">1 M</span> is 1 molar concentration
and <span class="math inline"><em>R</em></span> is the universal gas
constant. We determine <span
class="math inline"><em>Δ</em><em>G</em><sup>∘</sup></span> with the PMF
graph we have obtanied. In order to do that, we can use this
equation:</p>
<ol start="2" type="1">
<li></li>
</ol>
<p><span class="math display">$$
\Delta G^\circ = -k_B T \ln \frac{\int_{\mathrm{pocket}}
e^{-\frac{\Delta F(\xi)}{k_B T}} \, \mathrm{d}\xi}{\int_{\mathrm{bulk}}
e^{-\frac{\Delta F(\xi)}{k_B T}} \, \mathrm{d}\xi} = -k_B T
\frac{Z^\circ_{\mathrm{pocket}}}{Z^\circ_{\mathrm{bulk}}}
$$</span>
</p>
<ul>
<li><span class="math inline">pocket := {<em>ξ</em> : ligand at
state <em>ξ</em> is bound to the receptor}</span></li>
<li><span class="math inline">bulk := {<em>ξ</em> : ligand at
state <em>ξ</em> is not interacting with the receptor}</span></li>
</ul>
<p>So, what we have right now is the <span
class="math inline"><em>F</em>(<em>ξ</em>)</span> curve. Recognizing
these as the continous definitions of the partition function, we
simplify it down to:</p>
<ol start="3" type="1">
<li></li>
</ol>
<p><span
class="math display"><em>Δ</em><em>G</em><sup>∘</sup> = <em>F</em><sub>pocket</sub> − <em>F</em><sub>bulk</sub></span>
</p>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/bindingcoeff.png"></img>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/bindingcoeff1-3.png"></img>
</div>
</div>
<div id="parameters" class="mt-4">
<h2>Selecting Parameters</h2>
<p>The most important parameters are:
<ol type="1">
<li><p><span class="math inline"><em>L</em><sub><em>i</em></sub></span>
and <span class="math inline"><em>L</em><sub><em>f</em></sub></span>:
the starting and ending lengths of the complex you are trying to
stretch</p></li>
<li><p><span class="math inline">index1</span> and <span
class="math inline">index2</span>: the indexes of the two atoms we want
to pivot as endpoints and stretch accordingly.</p></li>
<li><p><span class="math inline">fc<sub>pull</sub></span>: the strength
of bias</p></li>
<li><p><span class="math inline">v<sub>pulling</sub></span>: the speed
of exeucting the pulling loop to streth the two pivot atoms (the
simulation will thus adjust the total time of the pulling loop based on
<span class="math inline"><em>L</em><sub><em>i</em></sub></span> and
<span class="math inline"><em>L</em><sub><em>f</em></sub></span> as
well)</p></li>
</ol>
<p>Virtually, <span
class="math inline"><em>L</em><sub><em>i</em></sub></span> can be easily
determined by looking at the original PDB file through a MD visualizer,
and <span class="math inline"><em>L</em><sub><em>f</em></sub></span> can
be any value <span
class="math inline"><em>L</em><sub><em>f</em></sub> ≫ <em>L</em><sub><em>i</em></sub></span>
(by maybe a factor less than <span
class="math inline">10<sup>2</sup></span> ish) if you are
stretching.</p>
<p>Additionally, since the simulation is purely discretized, <span
class="math inline">v<sub>pulling</sub></span> should be picked at a low
value to ensure it captures energy landscape adequetly.</p>
<p>However, the pivot indexes can be a little tricky to determine. But,
the best way to start off is by identifying your objective and looking
at the starting PDB to get a glipmse at what you’re working with.</p>
</p>
<img src ="https://static.igem.wiki/teams/5114/umbrellasampling/params.png"></img>
<p>However, the pivot indexes can be a little tricky to determine. But,
the best way to start off is by identifying your objective and looking
at the starting PDB to get a glipmse at what you’re working with.</p>
<p>
<p><strong>Generally Pick the strongest parts of the Complex</strong></p><pre><code>
<span class="hljs-number">1.</span> <span class="hljs-keyword">For</span> proteins: C-alpha atoms are picked because they<span class="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>
<span class="hljs-number">2.</span> <span class="hljs-keyword">For</span> nucleic acids: Phosphate atoms <span class="hljs-keyword">or</span> specific base atoms (because again, they<span class="hljs-comment">'re the strongest backbones).</span>
<span class="hljs-number">3.</span> <span class="hljs-keyword">For</span> small molecules <span class="hljs-keyword">or</span> ligands: <span class="hljs-keyword">Select</span> atoms that are involved <span class="hljs-keyword">in</span> <span class="hljs-keyword">key</span> interactions <span class="hljs-keyword">or</span> are structurally significant.
</code></pre>
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.
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><strong>Generally Pick the strongest parts of the Complex</strong></p>
<span class="hljs-number">1.</span> <span class="hljs-keyword">For</span> proteins: C-alpha atoms are picked because they<span class="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>
<span class="hljs-number">2.</span> <span class="hljs-keyword">For</span> nucleic acids: Phosphate atoms <span class="hljs-keyword">or</span> specific base atoms (because again, they<span class="hljs-comment">'re the strongest backbones).</span> <br>
<span class="hljs-number">3.</span> <span class="hljs-keyword">For</span> small molecules <span class="hljs-keyword">or</span> ligands: <span class="hljs-keyword">Select</span> atoms that are involved <span class="hljs-keyword">in</span> <span class="hljs-keyword">key</span> interactions <span class="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>