Purpose
GMCT and
GCEM
form a pair of programs for studying the binding thermodynamics of
macromolecular receptors.
GCEM
is a program for the automated preparation of the necessary input
for
GMCT
from a continuum electrostatics / molecular mechanics model.
The software allows a detailed modeling of complex
ligand-receptor systems in structure based calculations.
The primary targets of the software are biomolecular
receptors like proteins that bind or transfer
protons, electrons or other small-molecule ligands.
The software might, however, also be useful for
studying polyelectrolytes in a broader sense or
other systems like Ising or Potts models.
The properties of large systems can be studied using a
variety of modern simulation methods.
This makes the software ideally suited to study the thermodynamics
of ligand binding
and charge transfer processes in bioenergetic
complexes and other complex biomolecular systems.
The receptor model
The image on the right shows a schematic view of the
receptor model of
GMCT
with all possible features at the example of a
membrane-intrinsic receptor.
The image gives explanations on mouse-over if
your browser supports
SVG.
-
The receptor with the bound ligands is described in explicit detail
while the receptor environment with the unbound ligands is modeled implicitly.
-
Within the continuum electrostatics model of
GCEM,
the receptor and its environment can be further subdivided into regions
of differing dielectric constants and accessibility to mobile ions.
Details about the continuum electrostatics model are given
→ here.
-
There can be multiple types of ligands, where the unbound ligands are
represented implicitly by their electrochemical potentials.
-
The receptor can possess multiple global conformations.
-
The receptor can possess multiple sites.
-
Each site can possess multiple instances, where an instance is a certain form of a site
which is characterized by binding topology, charge distribution, conformation and the
number of bound ligands of each type. With this description of the sites one can represent,
e.g., multiple redox and/or protonation forms, tautomeric forms, sidechain rotamers and spin states.
-
For a realistic modeling of trans-membrane proteins, it is possible to account
for an electrochemical trans-membrane potential that consist of a chemical
and an electrical component. The chemical component arises from concentration differences
of compounds between the two phases, while the electrical component arises either from charge
separation across the membrane or from an externally applied electrostatic potential, as for
example in some electrophysiological experiments.
|
|
Theoretical basis
Our software rests on a general formulation of binding theory in terms of electrochemical potentials,
instead of not directly comparable quantities like pH value and reduction potential.
This formulation increases the transparency of calculation results by making them
particularly easy to interpret.
A detailed derivation of the formalism can be found in
→ this thesis.
The essentials of the formalism can be displayed below.
show more detail
The electrochemical potential of a compound measures its thermodynamic stability and how well
the compound is accommodated within its environment.
A low (more negative) electrochemical potential indicates a thermodynamically more favorable
state, while a high (more positive) electrochemical potential indicates a thermodynamically less favorable
state.
The electrochemical potential of a compound \(i\) is designated by the symbol \(\bar{\mu}\)
and the chemical potential is designated by the symbol \(\mu\).
The electrochemical potential is given by
\[
\bar{\mu}_{i} = \mu^{\circ}_{i} + \mathrm{z}_{i} \mathrm{F} \phi_{i} + \beta^{-1} \ln a_{i}
\]
where \(\mu^{\circ}\) is the standard chemical potential and
\(\bar{\mu}^{\circ} = \mu^{\circ}_{i} + \mathrm{z}_{i} \mathrm{F} \phi_{i}\)
is the standard electrochemical potential.
The second term on the right hand side describes the influence of an electrostatic potential
\(\phi_{i}\) at the location of the compound on \(\bar{\mu}\),
where \(\mathrm{z}_{i}\) is the net formal charge of compound \(i\)
and \(\mathrm{F}\) is the Faraday constant.
The last term on the right-hand side describes the dependence of the
electrochemical potential on the activity and thus on the concentration of the compound, where
\(\beta^{-1} = \mathrm{k}_{\mathrm{B}} \mathrm{T}\) is the thermal energy.
The logarithmic activity \(\ln a_{i}\) is equal to zero under standard conditions.
An integrated form of the Gibbs-Duhem equation relates the Gibbs energy
of the system in a certain microstate \(\ms{n}\)
to the electrochemical potentials
(partial molar free energies) \(\bar{\mu}_i\) and the
stoichiometric coefficients \(\nu_{i}\) of all components
\(i\) present in the system
\[
\Gmicro{n} = \sum\limits_{i} \nu_{i} \bar{\mu}_{i}
\]
Here, the electrochemical potentials of the
compounds are mutually interdependent, that is,
adding a compound to the system or altering the
configuration of the system can shift the
electrochemical potentials of all compounds.
The symbol \(E\) was chosen to distinguish the microstate energy
from a fully qualified free energy which is, in general, given by
the thermodynamic average over an ensemble of microstates.
Such an ensemble can be referred to as macrostate or substate of a system.
The receptor behavior is determined by relative electrochemical potentials
of the receptor species.
We can define the intrinsic energy of a receptor microstate as its
relative standard electrochemical potential with respect to a freely chosen
reference microstate
\[
\GintS{\ms{n}} = \bar{\mu}^{\circ}_{\ms{n}} - \bar{\mu}^{\circ,\mathrm{ref}}
\]
This definition can also be used if the receptor is subdivided into
receptor constituents.
In our formulation, the receptor constituents comprise the sites and
the so-called background that is formed by all parts of the receptor
that do not belong to any site.
Within the continuum electrostatics model, the microstate energy can be written as
a sum of contributions of the receptor constituents and pairwise interaction energies
between them
\[
\Gmicro{n} = \Gconf{c}
+ \sum\limits_{i=1}^{\Nsites}
\left(
\Gint{c}{i}{k}
-\sum\limits_{m}^{\Nlig} \nu_{c,i,k,m} \bar{\mu}_{m}
\right)
+ \sum\limits_{i=1}^{\Nsites} \sum\limits_{j=1}^{j \lt i}
\Ginter{c}{i}{k}{j}{l}
\]
The first term on the right-hand side is the global conformational energy, where
\(c\) is the global conformation of the receptor.
The second term on the right-hand side describes the contributions of the
individual sites \(i\) in their respective instances \(k\left(i\right)\)
and the energetic cost of removing the bound ligands from the surrounding solution.
The outer sum runs over all \(\Nsites\) sites and the inner sum runs over all
\(\Nlig\) ligand types and their stoichiometric coefficients \(\nu_{c,i,k,m}\)
bound to the instance of the site within the global receptor conformation.
The third term on the right-hand side sums over the interaction energies \(W\)
between all pairs of sites \(i, j\) in their respective instances
\(k\left(i\right),l\left(j\right)\). Like the intrinsic energies, the interaction
energies do also depend on the global conformation \(c\) of the receptor.
Observables and thermodynamic functions of the system are determined
by the partition functions of involved substates of the system or the overall
system.
The partition function of the system is given by
\[
\partitionf = \sum\limits_{\ms{n}}^{\Nstates}
\exp\left[-\beta \Gmicro{\ms{n}} \right]
\]
where the sum runs over all microstates \(\ms{n}\) of the system.
The partition function of a substate \(a\) is given by the sum
over all microstates of the system that are compatible with the definition
of the substate
\[
\partitionf_{a} = \sum\limits_{\ms{n}}^{\Nstates} \delta_{\ms{n},a}
\exp\left[-\beta \Gmicro{\ms{n}} \right]
\]
where \(\delta_{\ms{n},a}\) is equal to \(1\) if the microstate is
compatible with the definition of the substate \(a\) and equal to
\(0\) otherwise.
The definition of substate \(a\) could, for example, state that a
certain site binds two protons.
The probability of observing substate \(a\) in equilibrium
is given by
\[
p_{a} = \frac{\partitionf_{a}}{\partitionf_{\phantom{a}}}
\]
The free energy difference between two substates \(a\) and \(b\)
of a system is determined by their partition functions
\[
\Delta G_{a \rightarrow b} = -\beta^{-1}\ln\left[\frac{\partitionf_{b}}{\partitionf_{a}}\right]
\]
One could, for example, be interested in the free energy cost
of increasing the number of protons bound to a certain site from
one to two.
An analytical calculation of thermodynamic properties and observables
is in practice impracticable already for systems of moderate size,
because the number of possible microstates grows exponentially
with the number of sites.
Simulation methods make it possible to accomplish such calculations
even for large systems, like the protein complexes occurring in bioenergetics,
in acceptable time.
Monte Carlo (MC) simulation methods are often particularly efficient.
Simulation methods
Currently,
GMCT
offers two basic MC methods: Metropolis MC and Wang-Landau MC.
Unique features of
GMCT
are accurate and efficient
free energy calculation methods
that can
be used to calculate free energy differences
for freely definable transformations and
for the calculation of free energy measures
of cooperativity.
Namely, the free energy perturbation method,
thermodynamic integration,
the non-equilibrium work method and the
Bennett acceptance ratio method have been implemented.
The coupling between events in molecular systems can be
quantified with covariances or free energy measures of cooperativity.
Approximate semi-analytic methods
Since version 1.1,
GMCT
supports also calculations with the hybrid method
that combines exact statistical mechanics within clusters
of strongly interacting residues with a mean-field approach
for the description of inter-cluster interactions.
Details can be found in the user manual.
Applications
A variety of modern Monte Carlo simulation methods
can be used to study overall properties of the
receptor as well as properties of individual sites.
The description of the system in terms of discrete
microstates of the receptor and chemical potentials
of the ligands renders the simulations computationally
very inexpensive relative to all-atom simulations.
This computational efficiency enables very accurate
calculations of receptor properties with low
statistical uncertainty.
Properties of binding processes that can be calculated
are for example binding probabilities
(titration curves), binding free energies and binding constants.
These properties can be computed from a
microscopic viewpoint for studying the behavior of separate
sites or groups of sites in the receptor or from a macroscopic
viewpoint for studying the overall behavior of the receptor.
Midpoint reduction potentials \(\mathcal{E}_{1/2}\) and
\(\mathrm{p}K_{1/2}\) values can be derived from computed titration curves.
Binding free energies can be expressed in terms of
thermodynamically defined reduction potentials and
\(\mathrm{p}K_{\mathrm{a}}\) values.
The free energy calculation methods can also be used
to study charge transfer reactions, conformational transitions
and any other process that can be described within
the receptor model of
GMCT
.
A particularly interesting feature of GMCT is the possibility
to calculate free energy measures of cooperativity that
can be used to study the coupling of different processes
in the receptor.
An example of special interest in our lab is
the coupling between binding and transfer processes of
charged ligands in bioenergetic protein complexes.
GMCT
can also be helpful in setting up and complementing
molecular dynamics (MD) simulations.
The preparation of a protein structure for MD
simulations does often require the specification of protonation
states and tautomeric states occupied by titratable residues.
This information can be obtained from Metropolis MC calculations
with
GMCT
.
In addition, protonation state calculations
can be used to assess whether
the modeling of a protein or other polyelectrolyte
could require a constant-pH
MD method.
Documentation
GMCT
and the continuum electrostatics / molecular mechanics model of
GCEM
are described in →this paper.
See the subdirectory doc of the
GMCT
distribution for detailed documentation,
including the theoretical basis of both programs
(basically comprises the above mentioned paper with
more detail at some points plus a detailed description of all programs).
Information about the extended
MEAD
library utilized by
GCEM
and usage of the program can be found
→here.
Examples for the usage of
GMCT
can be found in the directory examples of the
GMCT distribution.
Examples for the usage of
GCEM
can be found in the directory examples/gcem of the
MEAD
distribution.
Downloads
GMCT
Version 1.2.5 provides a cleaner, more robust Makefile.
Set your favorite C++ compiler via the environment variable CXX,
otherwise the alias c++ is used or g++.
Version 1.2.4 fixed compilation issues with newer compilers
(GCC >= 7, Clang >= 8, ICC >= 2018).
The MEAD distribution was updated to version 2.3.5 which fixes
compilation issues with newer compilers.
Version 2.3.2 fixed a critical bug that caused
segfaults in calculations with very big and/or fine grids.
Details can be found on the
web page of the distribution.
The software is free software
in the sense of the definition given by the Free Software Foundation.
GMCT
and
GCEM
are distributed under the terms of the
GNU Affero General Public License as published by the
Free Software Foundation; either version 3, or (at your option) any later version.
For more details, see the documentation of
GMCT
found in the directory doc of the distribution.
MEAD
is distributed under the terms of the
GNU General Public License as published by the
Free Software Foundation; either version 1, or (at your option) any later version.
For more details, see the files README, INSTALLATION and COPYING of the distribution.
Feedback and bug reports
Your opinion and hints are very welcome.
Please provide detailed information about the program input and output when reporting a bug.
Often, running a program with a higher verbosity level (set the parameter blab to 3)
helps to clarify the source of an error.
Related software
There is a variety of software that utilizes a continuum electrostatics
approach for the description of binding equilibria.
In principle, program input for
GMCT
can also be generated with
alternative continuum electrostatics software some of which is included
in the list below (some reformatting and unit conversion could be necessary).
Likewise, the energy terms computed by
GCEM
should also be usable with other simulation software, as far as
the features are supported.
The following list contains the corresponding software known to us,
grouped according to the authoring research groups.
For details about the software, please refer to the respective
linked websites.
-
Donald Bashford
wrote the original
MEAD
software package.
Besides
GCEM,
MEAD
contains an older program called Multiflex that can also be used for the
calculation of continuum electrostatics energies, but does not provide all
of
GCEM's
,
features.
The program Redti can be used for the analytic computation of titration curves
using an approximate approach.
-
XMCTI is a program from
Paul Beroza
that can be used for Metropolis MC simulations of protonation equilibria.
-
HYBRID is a program from
Michael Gilson
that can be used for the computation of titration curves with an approximate method.
-
The group of Ernst-Walther Knapp
provides a pair of software for continuum electrostatics calculations
on protonation and reduction equilibria.
Karlsberg
allows Metropolis MC titrations as function of pH value and reduction potential.
Replica exchange parallel tempering and biased MC are also possible.
TAPBS allows the
calculation of the electrostatic energies used by
Karlsberg
by interfacing to the software
APBS.
-
Jens Erik Nielsen is the author of
various related software.
The WhatIf pKa package
allows protonation state calculations within a continuum-electrostatics model,
where the continuum electrostatics part utilizes the
Delphi
software.
A newer variant of the package named
pdb2pka
utilizes
APBS for the continuum electrostatics
part.
-
MCCE is a program for
continuum electrostatics calculations on protonation equilibria from
the group of Marilyn Gunner.
The program allows Metropolis Monte Carlo titrations of protonation and reduction equilibria.
MCCE
uses the programs DELPHI
or APBS
for the continuum electrostatics calculations.
-
The group of Antonio Baptista
provides several programs for continuum electrostatics calculations on protonation
and reduction equilibria.
The programs
Petite
and
MCRP
perform Metropolis MC simulations of protonation and reduction equilibria.
Correlations between pairs of sites can also be computed.
The electrostatics part of their software utilizes the
MEAD
software package.
-
The group of Alexei Stuchebrukhov
is also developing a pair of software for continuum electrostatics calculations
on protonation and reduction equilibria.
Monte
allows Metropolis MC titrations and
replica exchange parallel tempering.
pKip
allows the
calculation of the electrostatic energies used by
Monte
by interfacing to the software
APBS.
Currently, the programs do not seem to be publicly available.
-
DOPS is a program for
continuum electrostatics calculations on protonation equilibria from the
Antosiewiczs group.
The name appeared in an older publication, but the program seems not to be publicly available.
-
Zap TK is a
commercial programming library which can be used for
continuum electrostatics calculations on protonation equilibria and
also seems to contain functionality for Metropolis MC titrations.
The library can be licensed from
OpenEye Scientific Software.
