Ammm:Free energy simulations

Free Energy Calculations

• Equilibrium Methods (FEP, TI)
• Nonequilibrium Methods

Density of States

• $\Omega$ = density of states
• $1/h^{3N}$ = factor of proportionality (necessary to retrieve correct correspondence with high-temperature quantum-mechanical prediction)

• Calculations

• In canonical ensemble, $A=-\beta \ln {Q(N,V,T)}$

Ergodicity

• The ergodic hypothesis allows generating trajectories (i.e. MD, MC) to sample from phase-space to compute ensemble averages.
• The ensemble average of the property $f$ of a system is:

$\left\langle f\right\rangle =\int fP(\mathbf {x} ,\mathbf {p} _{x})d\mathbf {x} d\mathbf {p} _{x}$

• The ergodic hypothesis states that

$\left\langle f\right\rangle =\lim _{T\rightarrow \infty }{\frac {1}{T}}\int _{0}^{T}f[\mathbf {q} (t),\mathbf {p} (t)]dt$ for any $\mathbf {q_{0}} =\mathbf {q} (0),\mathbf {p_{0}} =\mathbf {p} (0)$

• Quasi-nonergodic - a system that does not explore phase space over time.
  • high energy barriers separating volumes of phase space
  • volume in phase space diffuses very slowly
• Proof of Quasi-ergodic hypothesis
• Proof of ergodic theorem

Free Energy Perturbation

Perturbation Formalism

• We want to calculate the free energy difference between state 0 and state 1 of Hamiltonians ${\mathcal {H}}_{0}(\mathbf {x} ,\mathbf {p} _{x})$ and ${\mathcal {H}}_{1}(\mathbf {x} ,\mathbf {p} _{x})$

• The difference in Helmholtz free energy between is $\Delta A=-{\frac {1}{\beta }}\ln {\frac {Q_{1}}{Q_{0}}}$ where $\beta =(k_{B}T)^{-1}$

• $Q_{0}$ and $Q_{1}$ are their corresponding partition functions $Q={\frac {1}{h^{3N}N!}}\iint \exp[-\beta {\mathcal {H}}]d\mathbf {x} d\mathbf {p} _{x}$

• $\Delta A=-{\frac {1}{\beta }}\ln {\frac {\iint \exp[-\beta {\mathcal {H}}_{1}(\mathbf {x} ,\mathbf {p} _{x})]d\mathbf {x} d\mathbf {p} _{x}}{\iint \exp[-\beta {\mathcal {H}}_{0}(\mathbf {x} ,\mathbf {p} _{x})]d\mathbf {x} d\mathbf {p} _{x}}}$

• $\Delta A=-{\frac {1}{\beta }}\ln {\frac {\iint \exp[-\beta \Delta {\mathcal {H}}(\mathbf {x} ,\mathbf {p} _{x})]\exp[-\beta {\mathcal {H}}_{0}(\mathbf {x} ,\mathbf {p} _{x})]d\mathbf {x} d\mathbf {p} _{x}}{\iint \exp[-\beta {\mathcal {H_{0}}}]d\mathbf {x} d\mathbf {p} _{x}}}$

• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A=-{\frac {1}{\beta }}\ln \iint \exp[-\beta \Delta {\mathcal {H}}(\mathbf {x} ,\mathbf {p} _{x})]P_{0}(\mathbf {x} ,\mathbf {p} _{x})d\mathbf {x} d\mathbf {p} _{x}}
  • since $P_{0}(\mathbf {x} ,\mathbf {p} _{x})={\frac {\exp[-\beta \Delta {\mathcal {H}}_{0}(\mathbf {x} ,\mathbf {p} _{x})]}{\iint \exp[-\beta {\mathcal {H}}_{0}(\mathbf {x} ,\mathbf {p} _{x})]d\mathbf {x} d\mathbf {p} _{x}}}$

• Finally, Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A=-{\frac {1}{\beta }}\ln \left\langle \exp[-\beta \Delta {\mathcal {H}}(\mathbf {x} ,\mathbf {p} _{x})]\right\rangle _{0}}

• For systems of with the same masses, $\Delta A=-{\frac {1}{\beta }}\ln \left\langle \exp[-\beta \Delta U]\right\rangle _{0}$

• We can also swap states 0 and 1 and use Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A=-{\frac {1}{\beta }}\ln \left\langle \exp[-\beta \Delta U]\right\rangle _{1}}
  • Although forward and backward calculations are formally equilvalent, they may have different convergence properties

• Forward and backward calculations can be combined (BAR).

Stratification

• P($\Omega$) is a probability distribution over the domain $\Omega$
• Partition $\Omega$ in to disjoint regions Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Omega _{i}\subset \Omega } so that $\bigcup _{i\epsilon I}\Omega _{i}=\Omega$
• Example
  • Consider a transformation of state 0 to state 1 described by parameter Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda }
  • Suppose states 0 and 1 are separated by a high-energy barrier
  • Partition Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda } in to a number of smaller intervals
  • P(Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda } ) can be recovered with great savings in computer time
• Decreases variance of ensemble average
• Provides solution for quasi-nonergodic systems

Importance sampling

• Sample regions in phase space that are important for estimating quantity of interest more frequently
• Results of the simulation need to be weighted to ensure that estimator is unbiased
• Example
  • Transformation from state 0 to state 1 along Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda }
  • Define P(Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda } ) to be the true distribution
  • Let $P'(\lambda )=P(\lambda )\exp[\beta \eta (\lambda )]$ where Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \eta (\lambda )} is the weighting factor
  • Then Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A=-\beta ^{-1}\ln {\frac {P(\lambda _{1})}{P(\lambda _{0})}}=-\beta ^{-1}\ln {\frac {P'(\lambda _{1})}{P'(\lambda _{0})}}+\eta (\lambda _{1})-\eta (\lambda _{0})}

Order Parameter

• Coupling parameters, Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda } , that are used to describe transformations from the reference system to the target one.
• Reaction coordinate/path - an order parameters that corresponds to the path along which the transformation takes place in nature
• Examples
  • Simple/coupled torsional degrees of freedom
  • End-to-end reaction coordinate (for protein folding)
  • Mutagenesis of atoms (relative free energy)
  • Distance separating chemical species (PMF)
  • Annihilation of nonbonded interactions (solvation)
• Choice of Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda } significantly affects efficiency, accuracy of calculations
  • Smooth energy landscapes require less intermediate states than rough ones

Topology

• Single Topology paradigm - a common topology is used for initial and final states
  • Missing atoms are treated as vanishing particles by setting nonbonded parameters to 0.
  • Bonded terms
    • If chemical bonds are not strongly deformed between initial and target states, bonded terms are often not considered
    • However, FEP and TI can simulations can be used as well.
  • Examples of systems that can use single topology

• Dual Topology paradigm -
  • States 0 and 1 contain different topologies
  • Atoms not common to 0 and 1 never interact during the simulation.
  • Bonded parameters are not perturbed

• • 

Thermodynamic Cycle

• Hydration free energy of a small solute
  • Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A_{hydration}=\Delta A_{0}-\Delta A_{1}}

• Binding free energy of ligand to protein
  • $\Delta A_{binding}=\Delta A_{0}-\Delta A_{1}$

• Relative binding free energy of ligand to protein
  • Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta \Delta A_{AB}=\Delta A_{bind\!-\!B}-\Delta A_{bind\!-\!A}=\Delta A_{1}-\Delta A_{0}}

Application of FEP

• Enthalpy and entropy can be computed based on their thermodynamic relationship using finite difference.
• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta S=\partial \Delta A/\partial T}
• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta U=\partial \beta \Delta A/\partial \beta }
• 

Thermodynamic Integration

• If Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \xi } is a reaction coordinate, the free energy is Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle A(\xi )=-k_{B}T\ln P(\xi )}
• For a sufficiently smooth function, free energy difference can be written as $A(\xi _{1})-A(\xi _{0})=\int _{\xi _{0}}^{\xi _{1}}dA/d\xi \,d\xi$ where Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\frac {dA}{d\xi }}=\left\langle {\frac {\partial {\mathcal {H}}}{\partial \xi }}\right\rangle _{\xi }}
• Order parameters can be torsion angle, bond lengths/angles between atoms or groups of atoms, hydration number, gyration radius

Constrained Simulations

• The Hamiltonian is supplemented by a Langrange multiplier

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\mathcal {H}}^{\lambda }\equiv {\mathcal {H}}+\lambda (\xi -\xi (\mathbf {x} ))=K+U+\lambda (\xi -\xi (\mathbf {x} ))}

• Since Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \xi (\mathbf {x} )} is constant, Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\ddot {\xi }}=0} we have

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \lambda =Z_{\xi }^{-1}\left(\sum _{i}{\frac {1}{m_{i}}}{\frac {\partial \xi }{\partial x_{i}}}{\frac {\partial U}{\partial x_{i}}}-{\dot {\mathbf {x} }}^{t}\cdot \mathbf {H} \cdot {\dot {\mathbf {x} }}\right)} where $Z_{\xi }=\sum _{i}{\frac {1}{m_{i}}}\left({\frac {\partial \xi }{\partial x_{i}}}\right)^{2}$ and Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \mathbf {H} _{ij}={\frac {\partial ^{2}\xi }{\partial x_{i}\partial x_{j}}}} is the hessian of Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \xi }

• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle A(\xi )=\int \langle \lambda \rangle _{\xi ,{\dot {\xi }}}d\xi -k_{b}T\ln \langle Z_{\xi }^{-1/2}\rangle }

Adaptive Biasing Force Method

• 
• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\frac {dA}{d\xi }}=-\left\langle {\frac {d}{dt}}\left(m_{\xi }{\frac {d\xi }{dt}}\right)\right\rangle } where Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle m_{\xi }=Z_{\xi }^{-1}}
• Using velocity verlet integrator where:
  • $v(t+dt/2)=v(t-dt/2)+dt\;a(t)$
  • Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle x(t+dt/2)=x(t)+dt/2\;v(t+dt/2)}
• Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\frac {d}{dt}}\left(m_{\xi }{\frac {d\xi }{dt}}\right)={\frac {1}{2}}\left({\frac {p_{\xi }^{+}(t+dt)-p_{\xi }^{+}(t)}{dt}}+{\frac {p_{\xi }^{-}(t+dt)-p_{\xi }^{-}(t)}{dt}}\right)}

where Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle p_{\xi }^{+}(t)\equiv m_{\xi }(t)grad(\xi (t))\cdot \left[v(t+dt/2)-{\frac {dt}{6}}a(t+dt)\right]} and Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle p_{\xi }^{-}(t)\equiv m_{\xi }(t)grad(\xi (t))\cdot \left[v(t-dt/2)-{\frac {dt}{6}}a(t-dt)\right]}

• Ramping function is used to prevent nonequilibrium effects and systematic bias in the beginning of the simulation.

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle R(n)={\begin{cases}n/N_{0}&{\textrm {if}}\;n\leq N_{0}\\1&{\textrm {if}}\;n>N_{0}\\{\textrm {with}}\;N_{0}\cong 100\end{cases}}}

Examples

Deca-L-alanine

• ABF was used to probe the reversible unfolding of deca-L-alanine in vacuo
• The reaction coordinate is the distance separating first and last alpha carbon from 12-32 Angstrom
  • 5ns MD trajectory
• 
• Minimum at around 14 Angstrom
• 
• Free energy profile essentially identical to a reversible, 20-ns steered MD simulation
• Park et al.

Glycophorin A

• ABF was used to model the helix-helix association of Glycophorin A
• 
• Difficulty of this simulation is that once the two helices are close to each other, it is very difficult to slide one with respect to the other or change relative orientation due to steric clashes.
• Constrained sampling would keep positions in local minimum in phase-space.
•