Ammm:Umbrella sampling and WHAM

Reaction Coordinate

• A reaction coordinate is an abstract one-dimensional coordinate which represents progress along a reaction pathway.
  • Real coordinate system: bond length, bond angle, torsion
  • Non geometric parameters: bond order, Hydrogen bonds, RMSD
• Reaction coordinates are often plotted against free energy to demonstrate in some schematic form the potential energy profile associated to the reaction.

Umbrella Sampling (US)

Umbrella sampling is a method in computational physics and chemistry, which can enhance the sampling of different systems when it is hard to realize ergodicity due to the energy landscape. It involves the importance sampling in statistical mechanics and was first recommended by Torrie and Valleau in 1977.In many cases, the configuration space may have multiple high-energy barriers which are hard to conquer in the limited simulation time with common molecular dynamics (MD). Therefore, it is of high probability that the trajectory is trapped in the local minimum (local state) and leaves other inaccessible important states unsampled. In order to remove the impact of the energy barriers (activation case), the researchers can chosose the umbrella sampling (US) to bridge the gap between different configuration states by adding an extra bias potential. This bias potential should cancel out the influence of these energy barriers so that the energy landscape will become more smooth, making it easier to reach other states.

Why do we need biased potential?

Features of method

• Pre-determined collective variables are required to describe the configuration space related to reaction.
• The target of this method is to handle activation problems in sampling, where high energy barriers can be smoothened by adding extra bias potential. (Convert activation into diffusion)
• The extra bias potential causes the loss of temporal properties.
• This method supports parallel MD simulations with different bias potentials at the same time so that computational cost is reduced.
• This method demands post-processing to remove the effects of bias potential.

Weighted ensemble

The original ensemble can be modified by using a weight function. The weight function should be positive here.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle <exp(-\beta \Delta U)>_{0}= \frac{\int\, exp(-\beta \Delta U) exp(-\beta \Delta U_{0}) dx} {\int\, exp(-\beta \Delta U_{0}) dx} }

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle <exp(-\beta \Delta U)>_{0}= \frac{\int\, exp(-\beta \Delta U) w w^{-1} exp(-\beta \Delta U_{0}) dx} {\int\, exp(-\beta \Delta U_{0}) w w^{-1} dx}=\frac{<w^{-1} exp(-\beta \Delta U)>_{w}} {<w^{-1}>_{w}}}

In the equation above, the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle <...>_{w}} denotes the weighted ensemble average and the probability distribution function can be shown as:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \rho_{w}(x)=\frac{w exp(-\beta \Delta U_{0})} {\int\, w exp(-\beta \Delta U_{0}) dx}}

To reflect the connection between bias potential and the weight function, we can represent the weight function as an exponential form: $w=exp(-\beta V)$, where $V$ means the bias potential.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle <exp(-\beta \Delta U)>_{0}=\frac{<exp(\beta V) exp(-\beta \Delta U)>_{w}} {<exp(\beta V)>_{w}}}

Free energy calculation

After defining a collective variable (CV) $\xi$, we can calculate the free energy at any point along the direction of CV by using the probability distribution function. (The reference point is chosen at $\xi _{0}$)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta A(\xi)=-k_{B} T ln \frac {\rho_{0}(\xi)} {\rho_{0}(\xi_{0})}}

When the bias potential is added, the free energy change under the weighted ensemble can be given by:

$\Delta A_{w}(\xi )=-k_{B}Tln{\frac {\rho _{w}(\xi )}{\rho _{w}(\xi _{0})}}$

In this case, the original free energy change should be recovered by combining these two equations above: (the denominators have been offset in the ratio)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta A(\xi)=-k_{B} T ln \frac {\rho_{w}(\xi) w(\xi_{0})} {\rho_{w}(\xi_{0}) w(\xi)}}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta A(\xi)=\Delta A_{w}(\xi)+k_{B} T (ln w(\xi)-ln w(\xi_{0}))=\Delta A_{w}(\xi)-[V(\xi)-V(\xi_{0})]}

The best bias potential can help remove all the possible barriers on the energy landscape, which means the free energy change under this bias potential should be 0 everywhere. So, it is obvious to know that the best bias potential must counterbalance the true free energy.

$V(\xi )-V(\xi _{0})=-\Delta A(\xi )$

Challenges of Umbrella Sampling

• Pick up the best bias potential along the CVs to realize the sufficiently flat energy landscape
• Ensure complete scanning of all the configuration space along the CVs under bias potential
• Perform discretization of the probability distribution along CVs

WHAM

Historty

• Initially developed in 1989 by Ferrenberg and Swenden: Ferrenberg-Swenden reweighting
• Generalized by Kumar et al.: Weighted Histogram Analysis Method

Basic Idea

• The Weighted Histogram Analysis Method (WHAM) is a practical technique for the calculation of potentials of mean force (PMFs) from a group of umbrella sampling simulations. It can be applied to analyze the energy data based on not only real collective variables (have relations to the atom coordinates) but also virtual variables (temperature or interactions in the system). This technique can give out a systematic way to improve the choice of the bias potential from a set of different bias potentials so that it yields uniform sampling along the collective variables. Meanwhile, the multiple simulations can start from different points along collective variables, then the full sampling can be realized. Also, the underlying theory of this method provides us a practical way to do histogram analysis and calculate free energy under discretized language.

Discretization of configuration space along collective variable

Assume we have run $R$ simulations at the same time, with each using a different or same bias potential and temperature. Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle }

$U_{k}=U_{0}+V_{k}(\xi ),k=1,...,R$

The whole configuration space can be divided into Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle L} bins. The middle bins share the same width Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle d} while the beginning and end bins have the width $d/2$. After the division, the configuration space becomes a histogram graph, and the energy value for each bin $i$ in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} should be given by the following equation: (energy for middle bins is set to be the value of the middle point of this bin, while the energy of the beginning bin is set to be the value at the starting point and the energy of the end bin is set to be the value at the endpoint)

$V_{k,i}={\begin{cases}V_{k}(\xi _{i}),&(i=1,L+1)\\V_{k}({\frac {\xi _{i}+\xi _{i+1}}{2}}),&(i=2,...,L)\end{cases}}$

Then, the bias potential energy of the simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} along collective variable can be expressed in discretized language:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k}(\xi)=\sum_{i=1}^L V_{k,i} \chi_{[\xi_{i},\xi_{i+1}]}(\xi)}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \chi_{[\xi_{i},\xi_{i+1}]}} is characteristic function: if Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi} is in the range Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [\xi_{i},\xi_{i+1}]} then it takes value Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1} , otherwise it is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0} .

From the operation of dividing the whole configuration space into discrete bins, it is obvious that more bins indicate higher resolution, but more snapshots from the simulations are required. To make the expression more clear, we directly use Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k,i}(\xi)} to denote the function Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k,i} \chi_{[\xi_{i},\xi_{i+1}]}(\xi)} in the following.

State density & probability distribution function

The most important property in WHAM is state density (or energy state density) Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega} . In an NVT ensemble, the generalized state density can be shown as:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(N,V,U)=\frac{1} {h^{3N} N!} \iint\limits_{V^{3N}} \delta (H(\mathbf{q},\mathbf{p})-U)\, d\mathbf{q}\,d\mathbf{p}}

then Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(N,V,U)dU} represents the volume of microstates in Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [U-dU/2,U+dU/2]} . The partition function can be illustrated by using state density:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q(N,V,T)=\int\limits exp(-\beta U) \Omega(N,V,U)\, dU}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q(N,V,T)=\frac{1} {h^{3N} N!} \iint\limits_{V^{3N}} exp(-\beta U(\mathbf{q},\mathbf{p}))\, d\mathbf{q}\,d\mathbf{p}=\frac{1} {\Lambda^{3N} N!} \int\limits exp(-\beta U(\mathbf{q}))\, d\mathbf{q}=\frac{1} {\Lambda^{3N} N!} Z(N,V,T)}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Z(N,V,T)} is the configurational integral. From these equations above, the probability distribution function can be expressed as:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \rho(U;N,V,T)=\frac{exp(-\beta U) \Omega(N,V,U)} {Z(N,V,T)}}

In WHAM, the probability and state density can be calculated by using statistical data: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_{k}(\xi)} is the number of counts at Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi} in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_{k}} is the total number of snapshots in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} .

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \rho_{k}(U(\xi),T_{k})=\frac{N_{k}} {n_{k}}=\frac{exp(-\beta_{k} U) \Omega_{k}(U)} {Z_{k}}}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{k}(U)=a_{k} \frac{N_{k}} {n_{k}} exp(-\beta_{k} U)}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_{k}} denotes the proportionality constant unique for k-th simulation (related to configurational integral), which can be modified further to make the following equations more pretty: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_{k}=exp(-\beta_{k} f_{k})} .

Derivation of WHAM

When bias potential is shown explicitly, the state density becomes:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{k}(\xi)=\frac{N_{k}} {n_{k}} exp(\beta_{k} U_{0}(\xi)+\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)-\beta_{k} f_{k})}

Because the original energy distribution is the same for all the simulations and will be reduced in the following deduction, then we can remove that in state density.

$\Omega _{k}(\xi )={\frac {N_{k}}{n_{k}}}exp(\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )-\beta _{k}f_{k})$

The state density with the best bias potential can be given by a linear combination of state density functions from all $R$ simulations, with each one being assigned a weight $\omega _{k}$. The sum of all the weights should be $1$.

$\Omega (\xi )=\sum _{k=1}^{R}\omega _{k}\Omega _{k},\sum _{k=1}^{R}\omega _{k}=1$

To get the best state density, we need to minimize the mean square deviation (MSD) $\delta ^{2}\Omega$

$\delta ^{2}\Omega =\sum _{k=1}^{R}\omega _{k}^{2}\delta ^{2}\Omega _{k}$

$\delta ^{2}\Omega _{k}={\frac {1}{n_{k}^{2}}}exp(\sum _{i=1}^{L}\beta _{k}V_{k,i}-\beta _{k}f_{k})\delta ^{2}N_{k}$

The MSD of $N_{k}$ can be given by multiplication of a correlation factor and its average: (the correlation factor is roughly equal if the division and running time are the same for all simulations, so it can be canceled)

$\delta ^{2}N_{k}(\xi )=g_{k}{\bar {N_{k}}}(\xi )$

An estimate can be made to obtain the average ${\bar {N}}_{k}$:

$\Omega (\xi )={\frac {{\bar {N_{k}}}(\xi )exp(\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )-\beta _{k}f_{k})}{n_{k}}}$

which can be compared with $\Omega _{k}(\xi )={\frac {N_{k}(\xi )exp(\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )-\beta _{k}f_{k})}{n_{k}}}$.

Then: ${\bar {N_{k}}}(\xi )=n_{k}\Omega (\xi )exp(-\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )+\beta _{k}f_{k})$.

The minimization can be finished by using the method of Lagrangian multipliers:

${\frac {\partial (\delta ^{2}\Omega )}{\partial \omega _{k}}}={\frac {\partial }{\partial \omega _{k}}}(\delta ^{2}\Omega +\lambda (\sum _{i}\omega _{i}-1))$

$2\omega _{k}\delta ^{2}\Omega _{k}+\lambda =0(k=1,...,R)$

$\omega _{1}\delta ^{2}\Omega _{1}=\omega _{2}\delta ^{2}\Omega _{2}=...=\omega _{R}\delta ^{2}\Omega _{R}=C$

Because $\sum {i}\omega _{i}=1$, then $C={\frac {1}{\sum _{k}{\frac {1}{\delta ^{2}\Omega _{k}}}}}$.

$\omega _{k}={\frac {\delta ^{2}\Omega _{k}}{\sum _{k}{\frac {1}{\delta ^{2}\Omega _{k}}}}}={\frac {g_{k}^{-1}n_{k}exp(-\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )+\beta _{k}f_{k})}{\sum _{m=1}^{R}g_{m}^{-1}n_{m}exp(-\sum _{i=1}^{L}\beta _{m}V_{m,i}(\xi )+\beta _{m}f_{m})}}$.

Then:

$\Omega (\xi )=\sum _{k=1}^{R}\omega _{k}\Omega _{k}(\xi )={\frac {\sum _{k=1}^{R}g_{k}^{-1}N_{k}(\xi )}{\sum _{m=1}^{R}g_{m}^{-1}n_{m}exp(-\sum _{i=1}^{L}\beta _{m}V_{m,i}(\xi )+\beta _{m}f_{m})}}$

In bin $i$:

$\Omega _{,i}={\frac {\sum _{k=1}^{R}g_{k}^{-1}N_{k,i}}{\sum _{m=1}^{R}g_{m}^{-1}n_{m}exp(-\sum _{i=1}^{L}\beta _{m}V_{m,i}+\beta _{m}f_{m})}}$

The un-normalized density at bin $i$ in simulation $K$ can be given by:

${\tilde {\rho }}_{K,i}={\frac {exp(-\beta _{K}V_{K,i})\sum _{k=1}^{R}g_{k}^{-1}N_{k,i}}{\sum _{m=1}^{R}g_{m}^{-1}n_{m}exp(-\sum _{i=1}^{L}\beta _{m}V_{m,i}+\beta _{m}f_{m})}}\ \ (*)$

Because $f_{K}$ is related to configurational integral:

$exp(-\beta _{K}f_{K})=\sum _{i=1}^{L}{\tilde {\rho }}_{K,i}\ \ (**)$

Self-consistent algorithm

In equations $(*)$ and $(**)$, the unknown values are $f_{K}$ and ${\tilde {\rho }}_{K,i}$. We can work out their values by using the self-consistent method.

At first, we can assume $f_{K}=0(K=1,...,R)$ and put them into Equation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (*)} to obtain ${\tilde {\rho }}_{K,i}(K=1,..,R;i=1,...,L)$.

Next step: Put all Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}} into Equation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (**)} to get the values of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f_{K}} .

Repeat the two steps above until the deviations of all these values are smaller than pre-set tolerance. The true free energy change between any two bins in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K} without bias potential can be calculated by:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta A_{K,j \to i}=-\beta_{K}^{-1} ln(\frac{\rho_{K,i}} {\rho_{K,j}})}

Applications

• US combined with FEP

• Potential of Mean Force

Examples

Butane

• Protocol
  • 18 independent simulations
    • 500ps
    • Restraint spring constant = 0.02 kcal/mol-deg
  • WHAM
    • 90 bins (40/bin)
    • Enforced periodcity

• Histograms from Individual Trajectories

• Histograms of Combined Trajectories

• Butane PMF

Ion Channel