docs/reference-manual/special/density-guided-simulation.rst

   1 .. _density-guided-simulation:
   2
   3 Applying forces from three-dimensional densities
   4 ------------------------------------------------
   5
   6 In density-guided simulations, additional forces are applied to atoms that depend
   7 on the gradient of similarity between a simulated density and a reference density.
   8
   9 By applying these forces protein structures can be made to "fit" densities
  10 from, e.g., cryo electron-microscopy. The implemented approach extends the ones
  11 described in \ :ref:`182 <refOrzechowski2008>`, and \ :ref:`183 <refIgaev2019>`.
  12
  13 Overview
  14 ^^^^^^^^
  15
  16 The forces that are applied depend on:
  17
  18  * The forward model that describes how atom positions are translated into a
  19    simulated density, :math:`\rho^{\mathrm{sim}}\!(\mathbf{r})`.
  20  * The similarity measure that describes how close the simulated density is to
  21    the reference density, :math:`\rho^{\mathrm{ref}}`, :math:`S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]`.
  22  * The scaling of these forces by a force constant, :math:`k`.
  23
  24 The resulting effective potential energy is
  25
  26 .. math:: U = U_{\mathrm{forcefield}}(\mathbf{r}) - k S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]\,\mathrm{.}
  27           :label: eqndensone
  28
  29 The corresponding density based forces that are added during the simulation are
  30
  31 .. math:: \mathbf{F}_{\mathrm{density}} = k \nabla_{\mathbf{r}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r})]\,\mathrm{.}
  32           :label: eqndenstwo
  33
  34 This derivative decomposes into a similarity measure derivative and a simulated
  35 density model derivative, summed over all density voxels :math:`\mathbf{v}`
  36
  37 .. math:: \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{.}
  38           :label: eqndensthree
  39
  40 Thus density-guided simulation force calculations are based on computing a
  41 simulated density and its derivative with respect to the atom positions, as
  42 well as a density-density derivative between the simulated and the reference
  43 density.
  44
  45 Usage
  46 ^^^^^
  47
  48 Density-guided simulations are controlled by setting ``.mdp`` options and
  49 providing a reference density map as a file additional to the ``.tpr``.
  50
  51 All options that are related to density-guided simulations are prefixed with
  52 ``density-guided-simulation``.
  53
  54 Setting ``density-guided-simulation-active = yes`` will trigger density-guided
  55 simulations with default parameters that will cause atoms to move into the
  56 reference density.
  57
  58 The simulated density and its force contribution
  59 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  60
  61 Atoms are spread onto the regular three-dimensional lattice of the reference
  62 density. For spreading the atoms onto the grid, the discrete Gauss transform is
  63 used. The simulated density from atoms at positions :math:`\mathbf{r_i}` at a
  64 voxel with coordinates :math:`\mathbf{v}` is
  65
  66 .. math:: \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{.}
  67           :label: eqndensfour
  68
  69 Where :math:`A_i` is an amplitude that is determined per atom type and may be
  70 the atom mass, partial charge, or unity for all atoms.
  71
  72 The width of the Gaussian spreading function is determined by :math:`\sigma`.
  73 It is not recommended to use a spreading width that is smaller than the
  74 grid spacing of the reference density.
  75
  76 The factor for the density force is then
  77
  78 .. math:: \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{.}
  79           :label: eqndensfive
  80
  81 The density similarity measure and its force contribution
  82 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  83
  84 There are multiple valid similarity measures between the reference density and
  85 the simulated density, each motivated by the experimental source of the
  86 reference density data. For the density-guided simulations in |Gromacs|, the following
  87 measures are provided:
  88
  89 The inner product of the simulated density,
  90
  91 .. math:: S_{\mathrm{inner-product}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
  92                 \frac{1}{N_\mathrm{voxel}}\sum_{v=1}^{N_\mathrm{voxel}} \rho^{\mathrm{ref}}_v \rho^{\mathrm{sim}}_v\,\mathrm{.}
  93         :label: eqndenssix
  94
  95 The negative relative entropy between two densities,
  96
  97 .. math:: S_{\mathrm{relative-entropy}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
  98            \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{.}
  99         :label: eqndensseven
 100
 101 The cross correlation between two densities,
 102
 103 .. math:: S_{\mathrm{cross-correlation}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
 104            \frac{\sum_{v}\left((\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})\right)}
 105            {\sqrt{\sum_v(\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})^2 \sum_v(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})^2}}\mathrm{.}
 106         :label: eqndenscrosscorr
 107
 108
 109 Declaring regions to fit
 110 ^^^^^^^^^^^^^^^^^^^^^^^^
 111
 112 A subset of atoms may be chosen when pre-processing the simulation to which the
 113 density-guided simulation forces are applied. Only these atoms generate the
 114 simulated density that is compared to the reference density.
 115
 116 Performance
 117 ^^^^^^^^^^^
 118
 119 The following factors affect the performance of density-guided simulations
 120
 121  * Number of atoms in the density-guided simulation group, :math:`N_{\mathrm{atoms}}`.
 122  * Spreading range in multiples of Gaussian width, :math:`N_{\mathrm{\sigma}}`.
 123  * The ratio of spreading width to the input density grid spacing, :math:`r_{\mathrm{\sigma}}`.
 124  * The number of voxels of the input density, :math:`N_{\mathrm{voxel}}`.
 125  * Frequency of force calculations, :math:`N_{\mathrm{force}}`.
 126  * The communication cost when using multiple ranks, that is reflected in a constant :math:`c_{\mathrm{comm}}`.
 127
 128 The overall cost of the density-guided simulation is approximately proportional to
 129
 130 .. math:: \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{.}
 131           :label: eqndenseight
 132
 133 Applying force every N-th step
 134 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 135
 136 The cost of applying forces every integration step is reduced when applying the
 137 density-guided simulation forces only every :math:`N` steps. The applied force
 138 is scaled by :math:`N` to approximate the same effective Hamiltonian as when
 139 applying the forces every step, while maintaining time-reversibility and energy
 140 conservation. Note that for this setting, the energy output frequency must be a
 141 multiple of :math:`N`.
 142
 143 The maximal time-step should not exceed the fastest oscillation period of any
 144 atom within the map potential divided by :math:`\pi`. This oscillation period
 145 depends on the choice of reference density, the similarity measure and the force
 146 constant and is thus hard to estimate directly. It has been observed to be
 147 in the order of picoseconds for typical cryo electron-microscopy data, resulting
 148 in a `density-guided-simulation-nst` setting in the order of 100.
 149
 150 Combining density-guided simulations with pressure coupling
 151 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 152
 153 Note that the contribution of forces from density-guided simulations to the
 154 system virial are not accounted for. The size of the effect on the
 155 pressure-coupling algorithm grows with the total summed density-guided simulation
 156 force, as well as the angular momentum introduced by forces from density-guided
 157 simulations. To minimize this effect, align your structure to the density before
 158 running a pressure-coupled simulation.
 159
 160 Additionally, applying force every N-th steps does not work with the current
 161 implementation of infrequent evaluation of pressure coupling and the constraint
 162 virial.
 163
 164 Periodic boundary condition treatment
 165 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 166
 167 Of all periodic images only the one closest to the center of the density map
 168 is considered.
 169
 170 The reference density map format
 171 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 172
 173 Reference input for the densities are given in mrc format according to the
 174 "EMDB Map Distribution Format Description Version 1.01 (c) emdatabank.org 2014".
 175 Closely related formats like ``ccp4`` and ``map`` might work.
 176
 177 Be aware that different visualization software handles map formats differently.
 178 During simulations, reference densities are interpreted as visualised by ``VMD``.
 179
 180 Output
 181 ^^^^^^
 182
 183 The energy output file will contain an additional "Density-fitting" term.
 184 This is the energy that is added to the system from the density-guided simulations.
 185 The lower the energy, the higher the similarity between simulated and reference
 186 density.
 187
 188 Adaptive force constant scaling
 189 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 190
 191 To enable a steady increase in similarity between reference and simulated
 192 density while using as little force as possible, adaptive force scaling
 193 decreases the force constant when similarity increases and vice versa. To avoid
 194 large fluctuations in the force constant, change in similarity is measured
 195 with an exponential moving average that smoothens the time series of similarity
 196 measures with a time constant :math:`tau` that is given in ps. If the exponential
 197 moving average similarity increases, the force constant is scaled down by
 198 dividing by :math:`1+\delta t_{\mathrm{density}}/tau`, where
 199 :math:`\delta t_{\mathrm{density}}` is the time between density guided simulation steps.
 200 Conversely, if similarity between reference and simulated density is decreasing,
 201 the force constant is increased by multiplying by :math:`1+2\delta t_{\mathrm{density}}/tau`.
 202 Note that adaptive force scaling does not conserve energy and will ultimately lead to very high
 203 forces when similarity cannot be increased further.
 204
 205 Future developments
 206 ^^^^^^^^^^^^^^^^^^^
 207
 208 Further similarity measures might be added in the future, along with different
 209 methods to determine atom amplitudes. More automation in choosing a force constant
 210 as well as alignment of the input density map to the structure might be provided.