May 18 - 22, 2009
In this talk I shall describe two kinds of soft matter system in which
the interfacial behaviour is of interest. Firstly, the adsorption of
polymers and peptides on solid surfaces is of growing interest because
of the possible applications in biotechnology. We have recently improved
a Monte Carlo simulation method for studying lattice polymers near to
surfaces, and I shall present some recent results. Secondly, I shall
discuss the simulation of confined liquid crystals, with special
emphasis on hard platelets. These idealized systems may be related to
experiments on colloidal suspensions.
We will present results on the application of the treecode method to the computation of the real space part of the Ewald sum in molecular dynamics simulation. We show that for reasonable parameters, the treecode speeds up the Ewald sum for large systems and reproduces the structural and dynamical properties of the liquid within error bars.
All-atom explicit solvent molecular dynamics simulations of (un)folding and other conformational changes in proteins remain a great challenge because of the long time scales involved. These long times are due to high free energy barriers between metastable states. Rare event simulation methodology, such as Replica Exchange (REMD) and Metadynamics can overcome these free energy barriers, and give thermodynamic equilibrium properties. However, they do not preserve the dynamical pathways needed to gain insight into kinetics and mechanism of conformational change. Transition Path Sampling (TPS) is a technique designed for harvesting an ensemble of unbiased dynamical transition paths, and allows computation of the rates and extraction of the reaction coordinates. In this presentation I will discuss the various rare event methods their applicability and limitations. I will illustrate these methods on several examples: the folding of the 20-residue Trp-cage, the folding of a 37-residue WW-domain and the unfolding step in the photocycle of the 125-residue bacterial sensor Photoactive Yellow protein (PYP). For each of these systems we have studied the mechanism, reaction coordinates, folding rates, and/or barriers heights. This work shows that by employing a toolbox of computational rare event techniques we can overcome the timescale gap inherent in Molecular Dynamics of rare events and gain more insight in complex transitions such as protein folding and unfolding.
Metropolized integrators for ergodic
stochastic differential equations (SDE) are exhibited which (i) are
ergodic with respect to the (known) equilibrium distribution of the
SDE and (ii) approximate pathwise the solutions of the SDE on
finite time intervals. Both these properties are demonstrated
including precise strong error estimates. It is also
shown that the Metropolized integrator retains these properties
even in situations where the drift in the SDE is nonglobally
Lipschitz, and vanilla explicit integrators for SDEs typically
become unstable and fail to be ergodic.
Invariant theory concerns an algebra K[V]G: the polynomials on a
vector space V over a field K that are invariant under the action
of a group G that is represented on V. In the case of a finite
group acting on a finite-dimensional vector space over a field of
characteristic zero, and in other cases not of present concern, this
algebra is finitely generated in the Cohen-Macaulay form of a
finite-dimensional free module over a polynomial ring. Generators of
that ring (primary invariants) and a basis for the module (secondary
invariants) may be found by computer algebra. Molecular modelling
encompasses many things: molecular dynamics for the study of reaction
and isomerization processes, diffusion Monte Carlo for ground state
properties, path-integral Monte Carlo for thermal properties, and
direct diagonalization of a configuration-interaction matrix for
calculation of a ro-vibrational spectrum. A key tool in molecular
modelling for any specific system is an analytical, fitted potential
energy surface: the potential energy of the system as a function of
the spatial configuration of nuclei, as is obtained from first
principles by solution of the electronic Schroedinger equation. In
the work to be described here molecular modelling meets invariant
theory and computer algebra in the construction of potential energy
surfaces. A polynomial model in suitable variables is employed and
the model is explicitly invariant under the full molecular permutation
symmetry group. The MAGMA computer algebra system is employed to
obtain the generators. Using such a representation we are able to
construct quite routinely highly accurate potential energy surfaces
for systems of N=7 or more atoms in full (3N-6) dimensionality,
and with use of a many-body expansion also for larger systems such as
This project is devoted to the theoretical analysis of the
zero-temperature string method to identify minimum
energy paths (MEP's) on a given energy landscape. The MEP's are curves
connecting local minima on the energy landscape which are everywhere
tangent to the gradient of the potential except possibly at critical
points. In practice, MEP's are mountain pass curves that play a special role
e.g. in the context of rare reactive events that occur when one
considers a steepest descent dynamics on the potential perturbed by a
small random noise. The string method aims at identifying MEP's by
moving each point of the curve by steepest descent on the energy landscape.
Here we address the question of whether such a curve evolution
necessarily converges to a MEP. Surprisingly, the answer is no, for
an interesting reason: MEP's may not be isolated, in the sense that
there may be families of them that can be continuously deformed into
one another. This degeneracy is related to the presence of critical
Morse index 2 or higher along the MEP. In this paper, we elucidate
this issue and characterize completely the limit set of a curve
evolving by the string method. We establish rigorously that the limit
set of such a curve
is again a curve when the MEP's are isolated. We also show under the
that the string evolution
to an MEP. However, we also identify and classify situations where the limit set
a curve and may contain higher-dimensional
parts. We present a collection of examples where the limit set of a
path contains a 2D region, a 2D surface, or a 3D body. In some of
our examples the evolving path wanders around without converging to its limit set.
In other examples it fills a region, converging to its limit set
which is not an MEP.
We present the need (especially critical for multiscale simulations) of extending standard equilibrium Molecular Dynamics to nonequilibrium situations. We show how this is possible in full rigor by use of the nonlinear Kubo-Onsager relation, connecting dynamical nonequilibrium averages to equilibrium (more genarally, initial distribution) averages over suitable time-evolved variables. To show the power of the method we apply it to study the onset of the convective circulation in liquids establishing as stationary asymptotic state the convective rolls know also as Benard cells.
We have used the ring polymer molecular dynamics method to study the
Azzouz-Borgis model for proton transfer between phenol (AH) and trimethylamine
(B) in liquid methyl chloride. When the A-H distance is used as the reaction
coordinate, the ring polymer trajectories are found to exhibit multiple
recrossings of the transition state dividing surface and to give a rate
coefficient that is smaller than the quantum transition state theory value by an
order of magnitude. This is to be expected on kinematic grounds for a
heavy-light-heavy reaction when the light atom transfer coordinate is used as
the reaction coordinate, and it clearly precludes the use of transition state
theory with this reaction coordinate. As has been shown previously for this
problem, a solvent polarization coordinate defined in terms of the expectation
value of the proton transfer distance in the ground adiabatic quantum state
provides a better reaction coordinate with less recrossing. These results are
discussed in light of the wide body of earlier theoretical work on the
Azzouz-Borgis model and the considerable range of previously reported values for
its proton and deuteron transfer rate coefficients.
Joint work with Amirali Kia and Jose Solomon.
The Mori-Zwanzig formalism is an effective tool to derive differential equations describing the evolution of a small number of resolved variables. We will present its application to the derivation of generalized Langevin equations and generalized non-Markovian Fokker-Planck equations. We show how long time scales rates and meta-stable basins can be extracted from these equations. Numerical algorithms are proposed to discretize these equations. An important aspect is the numerical solution of the orthogonal dynamics equation which is a partial differential equation in a high dimensional space. We propose efficient numerical methods to solve this orthogonal dynamics equation. In addition, we present a projection formalism of the Mori-Zwanzig type that is applicable to discrete maps. Numerical applications are presented from the field of protein modeling.
The study of the interaction between (large) molecules and inorganic (metal) surfaces is a typical multiscale problem. In fact while the specific chemistry emerges
at the quantum scale with the interaction of the electrons of the surface with the valence
electrons of the molecule, the molecular conformations are instead the expression of a large scale statistical behavior. The delicate interplay
between these two aspects, intimately linked during the adsorption process, gives rise to
the interface properties. I will show the application of a sequential modeling procedure which links the quantum to the classical atomistic scale in a consistent way.
Within this framework the ideal approach would be to identify and use only those degrees of freedom
strictly required by the problem. To achieve this goal I will illustrate the first step of the development of a spatially adaptive simulation scheme (AdResS) which allows to switch from the atomistic to a coarse-grained description on demand.
Anton, a recently completed special-purpose machine for molecular
dynamics simulations, achieves a speedup of several orders of magnitude
over state-of-the-art general-purpose systems by mapping novel
algorithms to massively parallel application-specific hardware. This
talk will describe algorithms and numerical techniques used by Anton,
including some developed especially for Anton and others adapted to
exploit its architecture. These algorithms underlie Anton's ability to
avoid communication bottlenecks even at high levels of parallelism, to
effectively map computation onto fast, hardwired arithmetic pipelines,
and to produce accurate results efficiently and deterministically.
Milestoning is a theory for stochastic coarse graining in space and time and an algorithm to compute kinetics and reaction dynamics of complex molecular systems. I will discuss the theory and algorithm properties, implementation, and an application to the folding of a helix and to the recovery stroke of the protein myosin. I will review the difficulties of simulating functions of biological molecules on relevant time scales and suggest that Milestoning addresses some of these problems. Folding pathways of a helical peptide with an energy landscape that is expected to be a funnel is illustrated. The myosin example shows coupling between local molecular events and global transitions which extends to the millisecond time scale.
A long-standing challenge has been to develop realistic atomistic-level
models of key surface reactions on single-crystal metal catalyst surfaces.
We have developed realistic multi-site lattice gas models for these
systems incorporating DFT energetics, and thus accurately describing
the complex nanoscale ordering of reactants . Model behavior is
determined from KMC simulation. A second challenge is to incorporate
these atomistic models into a multiscale treatment of spatio-temporal
behavior on a much larger scale (of microns). This is achieved using
our Heterogeneous Coupled Lattice Gas (HCLG) simulation strategy ,
together with a precise description of collective diffusion in these
mixed reactant adlayers. The result is a "first-principles" description
of reaction-diffusion front propagation .
 e.g., Liu, Evans, J. Chem. Phys. 124 (2006) 154705.
 Tammaro, Sabella, Evans, J. Chem. Phys. 103 (1995) 10277.
 Liu, Evans, Surf. Sci. (2009) doi:10.1016/j.susc.2008.10.058
The maintenance of the genetic information in a cell is essential for its survival. DNA polymerases play a central role in this process since they replicate and repair DNA. These enzymes must be versatile and specific to accommodate the four combinations of the Watson-Crick base pairs while preventing errors such as incorrect nucleotide insertion.
DNA polymerase (pol) lambda is a mammalian enzyme involved in the base excision repair (BER) and non-homologous end-joining (NHEJ) DNA repair pathways. It is from the X-family of polymerases and is shaped like a hand with finger, palm, and thumb subdomains. It also has 8-kDa and BRCT domains that aid its participation in the repair pathways. Pol lambda has an unusual fidelity profile since it produces many deletion errors and relatively few base substitution errors.
Substantial experimental and computational work on a related enzyme, pol beta, has elucidated the conformational changes induced upon its binding the correct and incorrect nucleotide. These motions involve a large-scale rearrangement of the thumb and subtle motions of active-site protein residues that serve as "gate-keepers" in regulating pol beta's transition to an active state. Binding of the incorrect nucleotide hampers these transitions.
Here we describe dynamics simulations of pol lambda/DNA complexes to uncover the important motions and rearrangements prior to chemistry that account for pol lambda's specific error profile. As in pol beta, a series of active-site protein residues are hypothesized to modulate the assembly of the active site. However, DNA motion is also involved. Arg517 emerges as a key residue for discriminating against incorrect nucleotide incorporation. This residue also plays an important role in the stabilization of the active position of the DNA even when the DNA is misaligned and poised for a deletion error.
Implicit solvent models, which treat the solvent as a macroscopic
continuum while admitting a microscopic atomic description for the
biomolecule, are efficient multiscale approaches to complex, large scale
biological systems. We summarize our recent advances in mathematical
methods for the Poisson-Boltzmann (PB) Model.
These include rigorous mathematical treatments of interface jump
conditions, surface geometric singularities, charge singularities, and
associated electrostatic force computation for molecular dynamics.
Large scale conformational rearrangements represent both a
challenge and an opportunity for rational drug design.
Exploring the conformational space of a target protein with sufficient detail is computationally very demanding and most current docking programs are very limited in this respect.
If it were possible, however, it could open the avenue to the design of
more selective drug candidates. Here we show how atomistic molecular
dynamics together with methods developed to accelerate rare events can
be used to study both kinase conformational transitions and ligand binding.
Using a new sampling method which is able to find the low free energy
channel between an initial and final state we determined the atomistic dynamics of the open-to-closed movement of the cyclin dependent kinase 5 (CDK5). 
The inactivation movement has a two-step mechanism. First
alpha-C-helix rotates about 45 allowing the interaction between Glu51
and Arg149, then T-loop has to refold to assume the closed conformation.
We found out the key role of Arg149 in the closure movement as the
link residue between the N-lobe and the C-lobe, allowing a concerted movement
of the two lobes. We do not find evidence for the DFG out flip. A complementary
method, Metadynamics,[2,3] was used to study the undocking path of a
congeneric series of ligands to CDK2. Also in this case the large scale dynamics
of the kinase plays a fundamental role.
1. D. Branduardi, F. L. Gervasio and M. Parrinello. Rare events
in molecular systems: From A to B in Free Energy Space
J. Chem. Phys., 126, 054103, 2007.
2. F. L. Gervasio*, A. Laio and M. Parrinello. Flexible Docking in Solution Using Metadynamics J. Am. Chem. Soc., 127(8): 2600-2607, 2005.
3. A. Laio and F. L. Gervasio* Metadynamics: a method to simulate rare
events and reconstruct the free energy in biophysics, chemistry and material science
Rep. Prog. Phys. , 71, 126601 (22p), 2008.
4. A. Berteotti, A. Cavalli, D. Branduardi, F. L. Gervasio*, M. Recanatini and M. Parrinello, Protein Conformational Transitions: The Closure Mechanism of a Kinase Explored by Atomistic Simulations nism of CDK5. J. Am. Chem. Soc.,131, 244-250, 2009.
Numerous simulation studies comparing the classical and quantum molecular
dynamics of simple empirical liquid water models have suggested that quantum
mechanical zero point energy and tunnelling increase the rates of translational
and orientational dynamics by a factor of about 1.4. In this work, we suggest
that the magnitude of this effect has been overestimated, principally as a
result of the use of either rigid water models or models in which
intramolecular flexibility is described by simple harmonic functions.
We have developed a new simple point charge model for liquid water, q-TIP4P/F,
in which O-H stretches are described by Morse-type functions. This model
correctly describes the liquid structure, translational diffusion and dipole
absorption frequencies in quantum (path integral-based) simulations, and also
gives a good description of the temperature-dependence of the liquid density.
By comparing classical and quantum simulations, we find that quantum effects
increase the rates of translation diffusion and orientational relaxation by a
factor of 1.15, an effect which is smaller than observed in all previous
simulations, regardless of quantum simulation method or water model. We
attribute the small quantum effect in our water model to two competing
phenomena. First, zero point energy weakens the intermolecular hydrogen-bonding
network resulting in a less viscous liquid and faster diffusion. Second,
intramolecular zero point energy changes the average water monomer geometry,
resulting in a larger molecular dipole moment and stronger intermolecular
interactions. By comparison to other simulations of liquid water models, we
suggest that the small quantum effect in our model is associated with its
ability to reproduce a single broad O-H stretching band in the infra-red dipole
I will present a computational method for simulating the dynamics of atomic systems on time scales much longer than can be accessed with classical dynamics. Possible reaction mechanisms available to the system are found by exploring the potential energy surface from minima to find nearby saddle points. Reaction rates are then calculated using harmonic transition state theory, and the system is propagated in time according to the kinetic Monte Carlo algorithm. The method can be run in parallel, and a distributed computing system has allowed us to simulate dynamics on metal surfaces and in grain boundaries over experimental time scales using an embedded atom potential. The algorithm is efficient enough to model the evolution of systems with ab-initio forces as well, for which I will show a few examples, including metal cluster formation on oxides and catalytic reactions on metal surfaces.
Replication of chromosomal DNA is accomplished by a dynamic protein assembly termed the replisome. Within the replisome, the function of replicative DNA polymerases is to rapidly and faithfully duplicate the cell’s genetic material prior to cell division. During replication proliferating cell nuclear antigen (PCNA, sliding clamp) serves as an accessory protein whose role is to topologically tether DNA polymerase to DNA. PCNA also acts as a scaffold in the recruitment of proteins involved in cell-cycle control and DNA repair. Mounting recent experimental evidence has pointed to the involvement of sliding clamps in almost every aspect of DNA metabolism. Yet many mechanistic details of how these fascinating proteins are loaded onto DNA and function within the machinery of the replisome have remained unknown.
To perform its vital functions the clamp has to be opened and resealed at DNA primer-template junctions by the action of a clamp loader ATPase – replication factor C (RFC). Recent computational work on the RFC/PCNA complex will be presented. It established that upon clamp opening the complex undergoes a large conformational rearrangement, leading to the formation of an extended interface between the surface of the sliding clamp and RFC. An interface complementary to ring-open PCNA transforms the free energy landscape underlying the closed- to open state transition, effectively trapping PCNA in an open conformation. Thus, careful comparison of free energy profiles for clamp opening in the presence and absence of RFC has allowed us to substantiate the role of the clamp loader in the initial stage of the clamp-loading cycle.
Joint work with Gernot J. Pauschenwein.
The success of experiments to design colloidal particles that interact
like classical hard spheres reintroduced lively work on systems
interacting via a hard core and various attached potential tails. We
investigated the equilibrium crystal structures at zero temperature of
hard--core systems interacting additionally either via soft shoulder
(polymer--grafted colloids) or Yukawa interactions (equally charged
colloids in salty solvent). Our investigations are based on
optimisation strategies that use ideas of genetic algorithms. The
problem of excluding a priori structures with overlapping cores
from search space was solved by introducing a particular
parametrisation of all crystal lattices. This parametrisation
drastically improved the performance of the genetic algorithm.
The application of this approach to the soft shoulder potential
revealed a large number of equilibrium structures in dependence of
ambient pressure, ranging from clusters over columns and lamellae to
finally compact structures. The close--packed structure with the
lowest energy per particle can, in dependence of the shoulder width,
be one with a different stacking of hexagonal layers than both fcc
(ABC) or hcp (AB), as the genetic algorithm revealed and was proven
theoretically. For hard core Yukawa systems the genetic algorithm
gives, as expected, fcc or bcc depending on pressure and interaction
range over the vast majority of parameter space. However our
investigation gives evidence that the phase transformation from bcc to
fcc at high pressure is happening through a continuous Bain
transformation, giving rise to a small regime of centered tetragonal
Joint work with Julia Fornleitner, Federica LoVerso, and Christos N. Likos.
Genetic algorithms (GAs) are powerful tools, applied in many fields of
physics as they offer an efficient strategy in high-dimensional
optimisation problems. We have developed GA-techniques that allow us
to study equilibrium crystal structures in an unbiased and
parameter-free way. Recently we have refined our approach and extended
it to binary mixtures of dipolar particles in a two-dimensional
geometry, realisable in experiments via, e.g. magnetic colloids on a
pendant water droplet and exposed to an external field or,
equivalently, via polystyrene particles floating at an
oil-water-interface. Several non-trivial and exotic structures are
discovered and their dependence on concentration and size ratio of the
particle species is discussed.
Protein association plays an essential role in many biological
processes. Understanding the general principles governing the
protein-protein interactions at the molecular level is of great
importance. In this study, we present a computational based alanine
scanning method that is used to elucidate the interactions between
alpha-lactamase and the alpha-lactamase inhibitor protein (BLIP). The
advantage of such approach is that it provides a direct
quantification of the specific per residue interactions between the
two proteins which can be directly validated with experimental
studies. Such analysis provides insight behind the structural and
energetic basis behind the molecular recognition process. Here, we
show the association between BLIP and TEM-1 is governed by clusters
of "hot spots" residues, which co-operatively modulate the formation
of the BLIP-TEM-1 complex. We also show that mutations in each
cluster do not influence the binding energy contribution of the
residues comprising other "hot spots" clusters, in agreement with
experimental results. The method presented here have broad range
application to study effect of mutation in drug binding as well as
de-novel protein design.
A treecode algorithm is presented for evaluating the electrostatic potential in charged particle systems undergoing
screened Coulomb interactions in 3D. The method uses a far-field Taylor expansion in Cartesian coordinates to compute particle-cluster interactions. The Taylor coefficients are evaluated using new recurrence relations which permit efficient computation of high order approximations. Two types of clusters are considered, uniform cubes and adapted rectangular boxes. The treecode error, CPU time, and memory usage are reported and compared with direct summation for randomly distributed particles inside a cube, on the surface of a sphere, and on an 8-sphere configuration. For a given order of Taylor approximation, the treecode CPU time scales as O(NlogN) and the memory usage scales as O(N), where N is the number of particles in the system. Results show that the treecode is well suited for non-homogeneous particle distributions as in the sphere and 8-sphere test cases. This is joint work with Peijun Li (Purdue) and Hans Johnston (U. Mass - Amherst), and has appeared in J. Comput. Phys., vol. 228 (2009) 3858-3868.
Molecular dynamics and replica-exchange simulations are performed for a model
peptide, NAcetyl-ALA5-Amide, using several popular force fields. We obtain helix populations,
melting curves, folding and nucleation times. The distributions of conformer populations are used to measure folding cooperativity. Finally, a statistical analysis of the sample of helix-coil transition paths is performed, uncovering interesting microscopic properties of the folding process.
In this work we present a mathematical approach to the study of diffusing particles confined to narrow geometries, such as water and ion channels in cell membranes. In this model, a strong potential (such as an electric field) drives the particle to remain close to the center of the channel, so that the dynamics in the direction parallel to the channel may be described by an averaged equation with an effective Fixman potential driving the motion of the particle. The Fixman potential depends on the particular geometry of the channel, accounting for entropic barriers (bottlenecks) whose effective heights are proportional to temperature. We compare transport properties such as flux rate and effective diffusivity when the typical barrier is either entropic or energetic.
A Kinetic Monte Carlo (KMC) model is developed to simulate non-symmetrically the cathode side of a Yttria Stabilized Zirconia (YSZ) fuel cell, in order to translate experimental, and ultimately theoretical rates into an atomistic model of
the fuel cell. The KMC model consists of a set of several electrochemical reaction rates, adopted from experiments and
first-principles calculations. The KMC simulations are used to model these simultaneously occurring events, and to determine potential limitations in cathode/YSZ performance. The focus of this work is ionic current density (J), studied as a function of various physical parameters: oxygen partial pressure, external applied bias voltage, operating temperature, dopant concentration, relative permittivity of YSZ, and geometrical features of the YSZ electrolyte. This simple model can be used as a baseline to translate elementary electrochemical reaction rates into atomistic simulations of a working solid oxide fuel cell cathodes, pertinent to the complete set of experimental operating conditions.
Computing canonical averages is a standard task in molecular dynamics. The main difficulty comes from the existence of metastable states. In order to facilitate sampling, and to obtain a better understanding of the system, low-dimensional reaction coordinates are often introduced. The statistics of these reaction coordinates are completely described by the associated free energy. However, the link between free energy barriers and dynamical informations, such as transition rates from one well to another, is less clear.
In this work, we design an effective, low-dimensional, dynamics on the reaction coordinate, which is a good approximation of the reference evolution, computed from the dynamics of the complete system. The accuracy of the effective dynamics is supported by error estimators. We also show, on a simple system, that our dynamics correctly reproduces the residence times in the wells.
For molecular dynamics to be useful for recovering macroscopically-relevant information (e.g. thermodynamic averages), it must be regulated by auxiliary control laws. Typically these would represent thermal or pressure controls. In this talk I will discuss a general family of stochastic-dynamic methods, including gentle thermostats which mildly perturb averaged dynamics and some very aggressive thermostats intended to enhance sampling of the canonical ensemble.
I will discuss a coarse-grained model of molecular dynamics in
crystalline solids. Such system usually has an underlying lattice
structure, and the empirical potentials are of quite simple forms.
Thus it provides a nice setting in which coarse-graining (CG) methods
can be considered. The goal of the CG method is to reduce the
atomic degrees of freedom, and arrive at an effective model in
which only a small number of atoms are explicitly involved.
Joint work with
Andrij Baumketnera, Wei
Shaozhong Dengb, Donald
Jacobsa, and Zhenli
A new solvation model is proposed for simulations of
biomolecules in aqueous solutions that combines the strengths
of explicit and implicit solvent representations. Solute
molecules are placed in spherical cavities filled with explicit
water, thus providing microscopic detail where it is needed.
Solvent outside of the cavities is replaced with a dielectric
continuum whose effect on the solute is modeled through the
reaction field corrections. With this explicit/implicit model,
the electrostatic potential represents a solute molecule in an
infinite bath of solvent, thus avoiding unphysical interactions
between periodic images of the solute commonly used in the
lattice-sum explicit solvent simulations. For improved
computational efficiency, our model employs an accurate and
efficient multiple-image method to compute reaction fields
together with the fast multipole-expansion technique for the
direct Coulomb interactions. To minimize the surface effects,
periodic boundary conditions are employed for non-electrostatic
interactions. The proposed model is applied to study liquid
water. The effect of main geometric parameters of the model,
which include the size of the cavity, the number of charge
images used to compute reaction field and the thickness of the
buffer layer, is investigated in comparison with the
particle-mesh Ewald simulations as a reference. An optimal set
of parameters is obtained that allows for a faithful
representation of many structural, dielectric and dynamic
properties of the simulated water, while maintaining manageable
a computational cost. With predicable increasing accuracy of
the multiple-image charge representation of the reaction field,
it is concluded that the proposed model achieves convergence
with only one image charge in the case of pure water. Future
applications to pKa calculations, conformational sampling of
solvated biomolecules and electrolyte solutions are briefly
a: Department of Physics and Optical Science, University of
North Carolina at Charlotte, Charlotte, NC 28223, United
b: Department of Mathematics and Statistics, University of
North Carolina at Charlotte, Charlotte, NC 28223, United
A prominent limitation of the commonly-used QM/MM methods is that no partial transfer is allowed between the QM and MM subsystems. A challenging example is the Eigen cation (Fig. 1), where one models the inner moiety H3O+ by QM and the surrounding three H2O molecules by MM. Both mutual polarization and partial charge transfer are expected between the QM and MM subsystems. This work aims to go beyond the above limit by incorporating partial charge transfer between the QM and MM subsystems, which (we hope) could lead to more accurate description for the QM/MM electrostatic interactions.
In this talk I will discuss some of the problems arising in connection with the dynamics of polymers under confinement (such as proteins traversing transmembrane pores). Confinement often partitions a polymer's configuration space into distinct basins of attraction separated by large free energy barriers. Because of the predominantly entropic nature of such barriers, transitions between different basins often cannot be characterized in terms of dominant pathways. I will particularly focus on the problem of a polymer reversal inside a long pore and compare exact rate calculations with predictions from polymer scaling theory and with various approximate theories.
Joint work with E. Laine (1), J.D.
Yoneda (2), and A. Blondel
Among the toxins secreted by Bacillus anthracis the edema factor EF, an adenylate
cyclase, provokes severe cellular dysfunction by accumulating cAMP from ATP. EF is
activated by calmodulin (CaM), involved in many calcium signaling pathways. The
stability of the EF-CaM complex depends on the level of calcium bound to CaM while
the architecture of the complex loaded with 2, 3 or 4 Ca2+ ions remains practically
unchanged. That is why modeling the electrostatic effect of Calcium through EF-
CaM structure is challenging.
Here, we aim at describing the calcium-induced changes in EF-CaM dynamics and
energetics through a consensual view of its residue network organization.
The analysis of molecular dynamics (MD) simulations of EF-CaM with 0,2 and 4 Ca2+
ions helped characterize CaM conformational plasticity and led to a model of the EF-
CaM interaction, in which CaM acts as a spring that maintains EF in an open active
conformation (Laine et al., 2008).
The complex was then analyzed (Laine et al.,2009) in order to extract simplified features describing its energetics and dynamics. The dynamics were analyzed using the LFA (Local Feature Analysis) (Zhang and Wrigers, 2006) and the generalized correlations (Lange and Grubmuller, 2006). These two approaches permitted to determine communication paths along the network of EF residues, from the helical domain to the domain CB. Beside, a new definition of the complex domains arose from the LFA analysis. This new definition, along with an approach derived from the energetic dependency maps (Hamacher et al, 2006), permitted to draw a simplified model of the energetic influences inside the complex. The variations of dynamics and energetics according to the level of Calcium complexation provide phenomenological reasons for the stability of the complex EF-CaM, by showing a disruption of the communication path or an imbalance of the energetic influences.
The computation of various dynamical covariances and energetic dependency maps
from the MD trajectories further raised the concept of residue network
connectedness. This connectedness quality provides a frame for unifying the
dynamics and energetics of the complex and a criterion for assessing its stability. .
Hamacher, K., J. Trylska, and J. McCammon. 2006. Dependency map
of proteins in the small ribosomal subunit. PLoS Comput. Biol. 2:e10.
Laine E, Blondel A, Malliavin TE (2009) Dynamics and energetics: a consensus analysis of the impact of calcium on EF-CaM protein complex.
Biophys J 96:1249-1263.
Laine E, Yoneda JD, Blondel A, Malliavin TE (2008) The conformational plasticity of calmodulin upon calcium complexation gives a model of its interaction with the oedema factor of Bacillus anthracis.
Lange OF, Grubmuller, H (2006) Generalized correlation for biomolecular dynamics. Proteins 62:1053-1061.
Zhang Z, Wriggers, W (2006) Local feature analysis: a statistical theory for reproducible essential dynamics of large macromolecules. Proteins 64:391-403.
(1) Unite de Bioinformatique Structurale,
Institut Pasteur, 28,
rue du Dr. Roux, F-75015 Paris, France
Programa de Posgraduac¸ao em Quimica Organica, Universidade
Federal Fluminense, Outeiro de S. Joao Batista s/n, 24020150,
Niteroi, RJ, Brazil.
A quantum simulation of an imaginary time path integral typically requires around n times more computational effort than the corresponding classical simulation, where n is the number of ring polymer beads (or imaginary time slices) used in the calculation. However, this estimate neglects the fact that the potential energies of many systems can be decomposed into a sum of rapidly-varying short-range and slowly-varying long-range contributions. For such systems, the computational effort of the path integral simulation can be reduced considerably by evaluating the long-range forces on a contracted ring polymer with fewer beads than are needed to evaluate the short-range forces. Using this approach it is shown that, for a typical application to a flexible model of liquid water, near classical simulation times can be achieved while still obtaining the exact path integral result.
A density operator treatment has been developed for the
dissipative quantum dynamics of molecular systems in condensed matter.
It incorporates instantaneous and delayed rates of dissipation
originating in electronic and vibrational motions of the molecular
environment, respectively. The theory is applied here to the optical
response of the adsorbates CO/Cu(001) and Ag/Si(111).
When a system is driven out of equilibrium by a time-dependent protocol that modifies the Hamiltonian, it follows a nonequilibrium path. Samples of these paths can be used in nonequilibrium work theorems to estimate equilibrium quantities, such as free energy differences. Here, we consider analyzing paths generated with one protocol using another one. It is posited that analysis protocols which minimize the lag, the difference between the nonequilibrium and the instantaneous equilibrium densities, will reduce the dissipation of reprocessed trajectories and lead to better free energy estimates. Indeed, when minimal lag analysis protocols based on exactly soluble propagators or relative entropies are applied to several test cases, substantial gains in the accuracy and precision of estimated free energy differences are observed.
The Activation-Relaxation Technique nouveau (ARTn) is an eigenvector following method for systematic search of saddle points and transition pathways on a given potential energy surface. We propose a variation of this method aiming at improving the efficiency of the local convergence close to the saddle point. We prove the convergence and robustness of this new algorithm. The efficiency of the method is tested in the case of point defects in body centered cubic iron.
We consider the mean field kinetic equations describing the relaxation dynamics of a lattice model of a fluid confined in a porous material. The dynamical theory embodied in these equations can be viewed as a mean field approximation to a Kawasaki dynamics Monte Carlo simulation of the system, as a theory of diffusion, or as a dynamical density functional theory. The solutions of the kinetic equations for long times coincide with the solutions of the static mean field equations for the inhomogeneous lattice gas. The approach is applied to a lattice gas model of a fluid confined in a finite length slit pore open at both ends and in contact with the bulk fluid at a temperature where capillary condensation and hysteresis occur. The states emerging dynamically during irreversible changes in the chemical potential are compared with those obtained from the static mean field equations for states associated with a quasistatic progression up and down the adsorption/desorption isotherm. In the capillary transition region the dynamics involves the appearance of undulates (adsorption) and liquid bridges (adsorption and desorption) which are unstable in the static mean field theory in the grand ensemble for the open pore but which are stable in the static mean field theory in the canonical ensemble for an infinite pore.
Molecular dynamics perturbs Hamiltonian dynamics in order to compute
correct ensemble averages. On the other hand computation of dynamical
averages such as autocorrelation functions requires the molecular
dynamics trajectory to be close to a microcanonical dynamics. This
suggest that, ideally we would like to have a small growth of the
perturbation and a fast rate of convergence to the Boltzmann-Gibbs
measure. Here we show how to achieve such a gentle thermostat by
combining a deterministic method (e.g. Nose-Hoover) with Langevin
dynamics. We use the concept of hypoellipticity of the forward
Kolmogorov equation, and Hörmander's condition, to investigate the existence of a unique invariant measure, which
implies ergodicity. We compare different methods in numerical
This is a joint work with Ben Leimkuhler (School of Mathematics, University of
Edinburgh) and Florian Theil (Mathematics Institute, University of Warwick).
 B. Leimkuler, E. Noorizadeh and F. Theil, J. Stat. Phys. 135, 261-277 (2009)
Ion channels gate the passage of ions through cell membranes
in response to external stimuli. In the case of the
archetypal potassium ion channel, KcsA, which opens and closes
upon a change in pH, recent experimental results have
demonstrated the existence of two physical gates: an
intracellular gate and a gate at the selectivity filter.
Lowering the pH opens the intracellular gate allowing ions to
pass. After the intracellular gate opens, however, the
channel can still inactivate by constricting the selectivity
filter and impeding the flow of ions even though the bottom
gate remains open. We use a combination of advanced
simulation techniques such as umbrella sampling, the string
method with swarms of trajectories, free energy perturbation
theory and Markov State models to probe the detailed mechanism
of this physiologically important process.
We present a triple-scale simulation of liquid water by concurrently coupling
the atomic, coarse-grained, and continuum descriptions of the liquid.
The triple-scale scheme, which is shown to correctly describe the
hydrodynamics, successfully sorts out the problem of large molecule
insertion in the hybrid particle-continuum simulations of molecular liquids.
The presented approach allows for efficient grand-canonical molecular dynamics
simulations of open systems involving an explicit solvent, e.g., liquid water.
We discuss theoretical and computational methodologies for quantitatively
describing how cell-membrane topologies are actively mediated and
manipulated by intracellular protein assemblies. Such scenarios are
ubiquitous in intracellular trafficking mechanisms, i.e., active transport
mechanisms characterized by vesicle nucleation and budding of the cell
membrane orchestrated by protein-interaction networks. We will describe the
development and application of the kinetic Monte Carlo time-dependent
Ginzburg Landau (KMC-TDGL) algorithm for unified dynamics of
curvature-inducing proteins and membranes. We will also describe the
surface-evolution method which predicts minimum energy conformations of
highly curved axis-symmetric membrane structures and an extension of the
TDGL method which relaxes the assumption of axis-symmetry. Together, these
methods enable us to quantify thermal effects in membrane mediated events.
We will briefly focus on an application to clathrin-dependent endocytosis.
Carbon nanofibers have attractive properties as support materials for
catalyst nanoparticles. Nonetheless, there is still a need for more
insights into the actual atomic structure of a metal particle interacting
with carbon nanostructures. Detailed information on the local structure of
adhered nanoclusters is still inaccessible to direct experimental imaging.
In contrast, molecular dynamics (MD) simulations can provide atomistic
structural information of a nanoscale system. In addition, the Reax force
field (ReaxFF) aims to achieve a good compromise between accuracy and
computational efficiency. We present novel results of MD simulations of
platinum and nickel clusters adsorbed on several different carbon supports
by using the ReaxFF. This work seeks to give insights into recent
experiments that have found a high correlation between carbon
nanofiber-induced strain and catalytic activity of Ni clusters. We show
how some metal atoms are detached from the surface of the nanoclusters
over time scales unreachable by electronic-structure methods. The adatom
migration is compensated by a rearrangement of the atomic structure of the
cluster. We also show differences in the bond length distribution of the
cluster when adsorbed to carbon nanocones with different radii as well as
differences found between Pt and Ni clusters. In addition, we brief on
preliminary results on diffusion of oxygen, hydrogen and carbon atoms over
the adsorbed cluster.
Liquid silicon has previously been shown to exhibit a liquid-liquid
transition in the supercooled state at zero pressure, using computer
simulations employing the Stillinger-Weber (SW) potential. The
associated liquid-liquid (LL) critical point lies at negative
pressures. Computer simulation evidence of such a negative pressure
critical point is presented, along with characterization of structural
and dynamical properties of the liquid in the vicinity of the critical
point. Preliminary results are presented on crystal nucleation, which
near the liquid-liquid critical point will be governed by the
interplay of critical density fluctuations and nearly tetrahedral,
crystal-like, local order in the liquid that grows with decreasing
Joint work with Marco Sarich (Freie Universität Berlin).
We map a discrete-time Markov Chain onto another Markov Chain describing transitions between some clusters only.
Every such aggregated Markov chain and affiliation function can be lifted again onto the full state space to define the so-called lifted transition matrix.
The optimal aggregated Markov chain and affiliation function can then be determined by minimizing some appropriately defined distance between the lifted transition matrix and the transition matrix of the original chain.
The talk will center around "conformation dynamics" and "rare
Given two metastable states A and B of a biomolecular system,
the problem is to compute a representative transition path by
which the mechanism of reaction can be analyzed. It is often
necessary to find such a path in collective variable space for
large biomolecular system. The maximum flux transition path (MFTP)
is defined as a path in collective variable space
that crosses each isocommittor at a point which (locally) has the
highest crossing rate of distinct reactive trajectories.
(Here the committor is defined to be the probability that a trajectory
at that point will reach B before A in collective variable space.)
Such an algorithm and computational performance will be discussed.
In rigidity theory, specifically combinatorial rigidity, one can simply count
vertices and edges (constraints) in a graph and its subgraphs to determine the
rigidity and flexibility of a corresponding framework. The 6|V|–6 counting
condition for 3-dimensional body-hinge structures (modulo molecular
conjecture), and a fast `pebble game' algorithm which tracks the underlying
count in the multigraph, have led to the development of the program FIRST, for
rapid analysis of the rigidity and flexibility of proteins and other large
macromolecular structures such as viral capsids. An important recent
application of FIRST is the introduction of the program FRODA which is a Monte
Carlo based geometric simulation technique, which uses the predictions of rigid
and flexible regions from FIRST to simulate the motions of proteins. FRODA can
run up to a million times faster than the Molecular Dynamics Simulations. We
will give a brief description and illustration of the 6|V| – 6 pebble game
algorithm. We further extend the pebble game algorithm to quantify the relative
degrees of freedom of a specified region (core) in the multigraph and identify
the regions that are relevant as constraints with respect to the core and those
regions that are irrelevant and that can be ignored both for focus and for speed
of simulations. With these new extensions (algorithms) we have developed a novel
approach in detecting hinge motions between protein domains. We can accurately
locate hinge points using only a single protein structure. Other application of
this work is to detect regions required for allostery; and the rigidity impact
of mutation studies. We present examples of these new extensions along with
The self-guided dynamics algorithms are a pair of techniques to
increase the speed of sampling from the possible conformations for a
macromolecule. Wu and Wang (1998), and later Wu and Brooks (2003),
proposed a numerical scheme that introduced a biasing term based on a
running average of past states of the system into the equations of
motion to increase the rate of barrier-crossings. In the present work,
an analytical form for the guiding force is derived from the numerical
methods. That analytical form is used to study the sampling properties
of the two algorithms. It is shown that, up to a separability
assumption, the systems do not sample from the Boltzmann distribution,
and the magnitudes of the error are estimated. Finally, a new
thermostat for the self-guided algorithms that conserves the Boltzmann
distribution is proposed.
I study from a mathematical viewpoint a nonlinear stochastic dynamics (called the Adaptive Biasing Force method), which allows to adaptively bias a metastable dynamics by computing the free energy associated with a chosen slow direction. Some applications in computational statistical physics as well as Bayesian statistics are presented.
The sequence dependent folding landscapes of nucleic acid
hairpins reflect much of the complexity of biomolecular
folding. Recently, mechanical folding trajectories, generated
using single molecule force clamp experiments by attaching
semiflexible polymers to the ends of hairpins have been used to
infer their folding landscapes. Using simulations and theory,
we study the effect of the dynamics of the attached handles on
the RNA free energy profile F(zm), where
zm is the molecular
extension of the hairpin. Accurate measurements of
requires stiff polymers with small L/lp, where L is the contour
length of the handle, and lp is the persistence length.
Paradoxically, reliable estimates of hopping rates can only be
made using flexible handles. We show that F(zm) at an external
tension fm, the force (f) at which the folded and unfolded
states are equally populated, in conjunction with Kramers'
theory, can provide accurate estimates of the force-dependent
hopping rates in the absence of handles at arbitrary values of
f. Our theoretical framework shows that zm is a good reaction
coordinate for nucleic acid hairpins under tension.
The concept of embedding some fragment in a larger system or reservoir appears everywhere in our scientific deliberations. Examples cover such diverse situations as an atom in a crystal, a functional group in a molecule, a molecular wire between a source and sink, or a protein in aqueous solution. Virtually every atomistic potential energy surface is built on one embedding strategy or another. This situation presents a difficulty for systems where charges are being transferred among fragments. Governance of transfer of charge is embodied in the concept of the chemical potential. When an embedded fragment and a reservoir can be readily identified, transfer of electrons between them can be described by an open ensemble. However, when interactions are strong as happens in crystals and molecules, the ensemble picture is not strictly applicable. Nevertheless, we still think in terms of open systems. Customary atomistic models use a linear model of chemical potential that is highly problematic. Here one strategy is presented to accommodate the open system point of view, while overcoming several limitations. That strategy focuses on decomposition into fragments of the many-body electronic hamiltonian itself, rather than decomposing the ground-state energy. The strategy produces a new, nonlinear model of chemical potential that has numerous ramifications for designing atomic potentials that can account for charge transfer.
Coauthors: J. L. F. Abascal, M. M. Conde and J. L. Aragones.
The performance of several popular water models (TIP3P, TIP4P,
TIP5P and TIP4P/2005) is analysed.
For that purpose the predictions for ten different properties
of water are investigated, namely: 1. vapour-liquid equilibria
(VLE) and critical temperature; 2. surface tension; 3.
densities of the different solid structures of water (ices); 4.
phase diagram; 5. melting point properties; 6. maximum in
density at room pressure and thermal coefficients α and
7. structure of liquid water and ice; 8. equation of state at
high pressures; 9. diffusion coefficient; 10. dielectric
constant. For each property, the performance of each model is
analysed in detail with a critical discussion of the possible
reason of the success or failure of the model. A final
judgement on the quality of these models is provided.
TIP4P/2005 provides the best description of almost all
properties of the list, with the only exception of the
dielectric constant. In the second position, TIP5P and TIP4P
yield an overall similar performance, and the last place with
the poorest description of the water properties is provided by
The ideas leading to the proposal and design of the TIP4P/2005
are also discussed in detail. TIP4P/2005 is probably close to
the best description of water that can be achieved with a non
polarizable model described by a single Lennard-Jones
(LJ) site and three charges.
It is often desirable to reduce the description of an atomic system to the dynamics of only a few relevant, usually, slow variables. Projection methods offer a rigourous framework to perform such reduction of variables. The main outcome of the method is the definition of a memory or friction kernel for the reduced dynamics. The memory kernel is equal to the time correlation of a so-called random force acting on the relevant variables. The evolution of the random force is not governed by the natural microscopic dynamics but by a projected dynamics. We have tried to illustrate this projected dynamics through the design of a method to directly compute the projected correlation functions from Molecular Dynamics simulations. Not only did this allow us to compute the memory kernel for the self-diffusion of Lennard-Jones particles but il also allowed for the decomposition of the memory kernel into short range and long range interactions contributions. Finally, long time behavior of the projected dynamics is investigated.
The fundamental model for diffusion is the Stokes-Einstein approach. It relates the diffusion coefficient to the friction coefficient which is in turn extracted from the hydrodynamic friction felt by a moving sphere in a liquid environment. The hydrodynamic flow around the sphere, commonly called Stokes flow, depends on the boundary conditions on the sphere. Two boundary conditions are usually considered: stick boundary conditions (no velocity at contact), as in the orginal approach, or slip boundary conditions (finite tangential velocity). As these two boundary conditions lead to different diffusion coefficient as a function of particle size, the usual way to investigate these from numerical simulations has been to compute the diffusion coefficient as a function of particle size. This leads for diffusion of Lennard-Jones particles to the suggestion of slip boundary conditions.
Here, we want to present a direct calculation of the microscopic velocity field around a diffusing particle from numerical simulations. This allow for comparison between the atomic model and the hydrodynamics approach. It is first demonstrated that the hydrodynamics flow is well recovered after only a few atomic radius from the tagged particle. However, two effects of the velocity fluctuations of the diffusing particle are evidenced. The boundary conditions are shown to include a finite normal velocity at contact, which would seem at first in contradiction with the non-penetrability of the particles but is an effect of fluctuations. Then, the flux of momentum in the flow is not determined solely by viscosity but also from the diffusion of the tagged particle.
This gives new insight in the diffusion process in liquids and in particular on the role of fluctuations at the atomic level. We will also try to show how this method could be also used in order to provide a more mechanistic description of the diffusion processes, in terms similar to reaction paths and transition states.
Coarse-graining the potential energy surface into the basins of
attraction of local minima, provides a computational framework for
investigating structure, dynamics and thermodynamics in molecular
science. Steps between local minima form the basis for global
optimisation via basin-hopping and for calculating thermodynamic
properties using the superposition approach and basin-sampling. To treat
global dynamics we must include transition states of the potential
energy surface, which link local minima via steepest-descent paths. We
may then apply the discrete path sampling method, which provides access
to rate constants for rare events. In large systems the paths between
minima with unrelated structures may involve hundreds of stationary
points of the potential energy surface. New algorithms have been
developed for both geometry optimisation and making connections between
distant local minima, which allow us to treat such systems. Applications
will be presented for a wide variety of atomic and molecular clusters,
glass-formers and biomolecules. Results for folding, misfolding and
aggregation of peptides and proteins illustrate how experimental time
and length scales can now be addressed for such systems.
The behavior of supercooled polymer melt composed of short
chains with ten beads between rapidly oscillating plates is
simulated by using a hybrid simulation of molecular dynamics
and computational fluid dynamics.
The flow profiles of polymer melt near an oscillating plate are
quite different from those of Newtonian fluid.
The viscous boundary layer of the melt is much thinner than
that of the Newtonian fluid due to the shear thinning of the
Three different rheological regimes, i.e., the viscous fluid,
viscoelastic liquid, and viscoelastic solid regimes, form over
the oscillating plate according to the local Deborah numbers.
The melt behaves as a viscous fluid when ftR < 1, and the
crossover between the liquid-like and solid-like regime takes
place around ftR = 1 (where f is the angular frequency of the
plate and tR and ta are the Rouse and the alpha relaxation
Conformational changes in proteins are intrinsic to almost every function, but sampling such changes remains one of the most difficult bio-computational problems. The challenge lies on two levels. First, algorithms which can efficiently focus computer time on the correct path ensemble are necessary. Second, even if perfect path-sampling algorithms were available, intrinsic "fast" timescales of protein motions are too costly for atomistic models; coarse-grained models are therefore required in most systems of interest. The group's work on both these issues will be described. Algorithmically, the "weighted ensemble" method is generalized and set on a firm statistical footing. Toward the goal of determining the most chemically realistic models which can be fully sampled, we describe a novel class of semi-atomistic models, based on pre-calculated residue libraries.