Applying forces from three-dimensional densities#

In density-guided simulations, additional forces are applied to atoms that depend on the gradient of similarity between a simulated density and a reference density.

By applying these forces protein structures can be made to “fit” densities from, e.g., cryo electron-microscopy. The implemented approach extends the ones described in 182, and 183.

Overview#

The forces that are applied depend on:

  • The forward model that describes how atom positions are translated into a simulated density, \(\rho^{\mathrm{sim}}\!(\mathbf{r})\).

  • The similarity measure that describes how close the simulated density is to the reference density, \(\rho^{\mathrm{ref}}\), \(S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]\).

  • The scaling of these forces by a force constant, \(k\).

The resulting effective potential energy is

(422)#\[U = U_{\mathrm{forcefield}}(\mathbf{r}) - k S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]\,\mathrm{.}\]

The corresponding density based forces that are added during the simulation are

(423)#\[\mathbf{F}_{\mathrm{density}} = k \nabla_{\mathbf{r}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]\,\mathrm{.}\]

This derivative decomposes into a similarity measure derivative and a simulated density model derivative, summed over all density voxels \(\mathbf{v}\)

(424)#\[\mathbf{F}_{\mathrm{density}} = k \sum_{\mathbf{v}}\partial_{\rho_{\mathbf{v}}^{\mathrm{sim}}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] \cdot \nabla_{\mathbf{r}} \rho_{\mathbf{v}}^{\mathrm{sim}}\!(\mathbf{r})\,\mathrm{.}\]

Thus density-guided simulation force calculations are based on computing a simulated density and its derivative with respect to the atom positions, as well as a density-density derivative between the simulated and the reference density.

Usage#

Density-guided simulations are controlled by setting .mdp options and providing a reference density map as a file additional to the .tpr.

All options that are related to density-guided simulations are prefixed with density-guided-simulation.

Setting density-guided-simulation-active = yes will trigger density-guided simulations with default parameters that will cause atoms to move into the reference density.

The simulated density and its force contribution#

Atoms are spread onto the regular three-dimensional lattice of the reference density. For spreading the atoms onto the grid, the discrete Gauss transform is used. The simulated density from atoms at positions \(\mathbf{r_i}\) at a voxel with coordinates \(\mathbf{v}\) is

(425)#\[\rho_{\mathbf{v}} = \sum_i A_i \frac{1}{\sqrt{2\pi}^3\sigma^3} \exp[-\frac{(\mathbf{r_i}-\mathbf{v})^2}{2 \sigma^2}]\,\mathrm{.}\]

Where \(A_i\) is an amplitude that is determined per atom type and may be the atom mass, partial charge, or unity for all atoms.

The width of the Gaussian spreading function is determined by \(\sigma\). It is not recommended to use a spreading width that is smaller than the grid spacing of the reference density.

The factor for the density force is then

(426)#\[\nabla_{r} \rho_{\mathbf{v}}^{\mathrm{sim}}\!(\mathbf{r}) = \sum_{i} - A_i \frac{(\mathbf{r_i}-\mathbf{v})}{\sigma} \frac{1}{\sqrt{2\pi}^3\sigma^3} \exp[-\frac{(\mathbf{r_i}-\mathbf{v})^2}{2 \sigma^2}]\,\mathrm{.}\]

The density similarity measure and its force contribution#

There are multiple valid similarity measures between the reference density and the simulated density, each motivated by the experimental source of the reference density data. For the density-guided simulations in GROMACS, the following measures are provided:

The inner product of the simulated density,

(427)#\[S_{\mathrm{inner-product}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] = \frac{1}{N_\mathrm{voxel}}\sum_{v=1}^{N_\mathrm{voxel}} \rho^{\mathrm{ref}}_v \rho^{\mathrm{sim}}_v\,\mathrm{.}\]

The negative relative entropy between two densities,

(428)#\[S_{\mathrm{relative-entropy}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] = \sum_{v=1, \rho^{\mathrm{ref}}>0, \rho^{\mathrm{sim}}>0}^{N_\mathrm{voxel}} \rho^\mathrm{ref} [\log(\rho^\mathrm{sim}_v)-\log(\rho^\mathrm{ref}_v)]\,\mathrm{.}\]

The cross correlation between two densities,

(429)#\[S_{\mathrm{cross-correlation}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] = \frac{\sum_{v}\left((\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})\right)} {\sqrt{\sum_v(\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})^2 \sum_v(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})^2}}\mathrm{.}\]

Declaring regions to fit#

A subset of atoms may be chosen when pre-processing the simulation to which the density-guided simulation forces are applied. Only these atoms generate the simulated density that is compared to the reference density.

Performance#

The following factors affect the performance of density-guided simulations

  • Number of atoms in the density-guided simulation group, \(N_{\mathrm{atoms}}\).

  • Spreading range in multiples of Gaussian width, \(N_{\mathrm{\sigma}}\).

  • The ratio of spreading width to the input density grid spacing, \(r_{\mathrm{\sigma}}\).

  • The number of voxels of the input density, \(N_{\mathrm{voxel}}\).

  • Frequency of force calculations, \(N_{\mathrm{force}}\).

  • The communication cost when using multiple ranks, that is reflected in a constant \(c_{\mathrm{comm}}\).

The overall cost of the density-guided simulation is approximately proportional to

(430)#\[\frac{1}{N_{\mathrm{force}}} \left[N_{\mathrm{atoms}}\left(N_{\mathrm{\sigma}}r_{\mathrm{\sigma}}\right)^3 + c_{\mathrm{comm}}N_{\mathrm{voxel}}\right]\,\mathrm{.}\]

Applying force every N-th step#

The cost of applying forces every integration step is reduced when applying the density-guided simulation forces only every \(N\) steps. The applied force is scaled by \(N\) to approximate the same effective Hamiltonian as when applying the forces every step, while maintaining time-reversibility and energy conservation. Note that for this setting, the energy output frequency must be a multiple of \(N\).

The maximal time-step should not exceed the fastest oscillation period of any atom within the map potential divided by \(\pi\). This oscillation period depends on the choice of reference density, the similarity measure and the force constant and is thus hard to estimate directly. It has been observed to be in the order of picoseconds for typical cryo electron-microscopy data, resulting in a density-guided-simulation-nst setting in the order of 100.

Combining density-guided simulations with pressure coupling#

Note that the contribution of forces from density-guided simulations to the system virial are not accounted for. The size of the effect on the pressure-coupling algorithm grows with the total summed density-guided simulation force, as well as the angular momentum introduced by forces from density-guided simulations. To minimize this effect, align your structure to the density before running a pressure-coupled simulation.

Additionally, applying force every N-th steps does not work with the current implementation of infrequent evaluation of pressure coupling and the constraint virial.

Periodic boundary condition treatment#

Of all periodic images only the one closest to the center of the density map is considered.

The reference density map format#

Reference input for the densities are given in mrc format according to the “EMDB Map Distribution Format Description Version 1.01 (c) emdatabank.org 2014”. Closely related formats like ccp4 and map might work.

Be aware that different visualization software handles map formats differently. During simulations, reference densities are interpreted as visualised by VMD.

Output#

The energy output file will contain an additional “Density-fitting” term. This is the energy that is added to the system from the density-guided simulations. The lower the energy, the higher the similarity between simulated and reference density.

Adaptive force constant scaling#

To enable a steady increase in similarity between reference and simulated density while using as little force as possible, adaptive force scaling decreases the force constant when similarity increases and vice versa. To avoid large fluctuations in the force constant, change in similarity is measured with an exponential moving average that smoothens the time series of similarity measures with a time constant \(tau\) that is given in ps. If the exponential moving average similarity increases, the force constant is scaled down by dividing by \(1+\delta t_{\mathrm{density}}/tau\), where \(\delta t_{\mathrm{density}}\) is the time between density guided simulation steps. Conversely, if similarity between reference and simulated density is decreasing, the force constant is increased by multiplying by \(1+2\delta t_{\mathrm{density}}/tau\). Note that adaptive force scaling does not conserve energy and will ultimately lead to very high forces when similarity cannot be increased further.

Mapping input structure to density data with affine transformations#

To align input structure and density data, a transformation matrix \(\mathbf{A}\) and shift vector \(\mathbf{v_{\mathrm{shift}}}\) may be defined that transform the input structure atom coordinates before evaluating density-guided-simulation energies and forces, so that

(431)#\[U = U_{\mathrm{forcefield}}(\mathbf{r}) - k S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{A} \mathbf{r}+\mathbf{v}_{\mathrm{shift}})]\,\mathrm{.}\]
(432)#\[\mathbf{F}_{\mathrm{density}} = k \nabla_{\mathbf{r}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{A} \mathbf{r}+\mathbf{v}_{\mathrm{shift}})]\,\mathrm{.}\]

Affine transformations may be used, amongst other things, to perform

  • rotations, e.g., around the z-axis by an angle \(\theta\) by using \(A=\begin{pmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix}\).

  • projection, e.g., onto the z-axis by using \(A=\begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\). This allows density-guided simulations to be steered by a density-profile along this axis.

  • scaling the structure against the density by a factor \(s\) by using \(A=\begin{pmatrix} s & 0 & 0 \\ 0 & s & 0 \\ 0 & 0 & s \end{pmatrix}\). This proves useful when, e.g., voxel-sizes in cryo-EM densities have to be adjusted.

  • and arbitrary combinations of these by matrix multiplications (note that matrix multiplications are not commutative).

Future developments#

Further similarity measures might be added in the future, along with different methods to determine atom amplitudes. More automation in choosing a force constant as well as alignment of the input density map to the structure might be provided.