Volume 137, Issue 23, 21 December 2012

The atomistic characterization of the transition state (TS) is a fundamental step to improve the understanding of the folding mechanism and the function of proteins. From a computational point of view, the identification of the conformations that build out the transition state is particularly cumbersome, mainly because of the large computational cost of generating a statistically sound set of folding trajectories. Here we show that a biasing algorithm, based on the physics of the ratchetandpawl, can be used to approximate efficiently the transition state. The basic idea is that the algorithmic ratchet exerts a force on the protein when it is climbing the freeenergy barrier, while it is inactive when it is descending. The transition state can be identified as the point of the trajectory where the ratchet changes regime. Besides discussing this strategy in general terms, we test it within a protein model whose transition state can be studied independently by plain molecular dynamics simulations. Finally, we show its power in explicitsolvent simulations, obtaining and characterizing a set of transitionstate conformations for AcylCoenzyme ABinding Protein (ACBP) and Chymotrypsin Inhibitor 2 (CI2).
 PERSPECTIVES


Perspective: Alchemical free energy calculations for drug discovery
View Description Hide DescriptionComputational techniques see widespread use in pharmaceutical drug discovery, but typically prove unreliable in predicting trends in proteinligand binding. Alchemical free energy calculations seek to change that by providing rigorous binding free energies from molecular simulations. Given adequate sampling and an accurate enough force field, these techniques yield accurate free energy estimates. Recent innovations in alchemical techniques have sparked a resurgence of interest in these calculations. Still, many obstacles stand in the way of their routine application in a drug discovery context, including the one we focus on here, sampling. Sampling of binding modes poses a particular challenge as binding modes are often separated by large energy barriers, leading to slow transitions. Binding modes are difficult to predict, and in some cases multiple binding modes may contribute to binding. In view of these hurdles, we present a framework for dealing carefully with uncertainty in binding mode or conformation in the context of free energy calculations. With careful sampling, free energy techniques show considerable promise for aiding drug discovery.
 Top

 COMMUNICATIONS


Communication: Nanoscale electrostatic theory of epistructural fields at the proteinwater interface
View Description Hide DescriptionNanoscale solvent confinement at the proteinwater interface promotes dipole orientations that are not aligned with the internal electrostatic field of a protein, yielding what we term epistructural polarization. To quantify this effect, an equation is derived from first principles relating epistructural polarization with the magnitude of local distortions in water coordination causative of interfacial tension. The equation defines a nanoscale electrostatic model of water and enables an estimation of protein denaturation free energies and the inference of hot spots for protein associations. The theoretical results are validated visàvis calorimetric data, revealing the destabilizing effect of epistructural polarization and its molecular origin.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

The tension of a curved surface from simulation
View Description Hide DescriptionThis paper demonstrates a method for calculating the tension of a system with a curved interface from a molecular dynamics simulation. To do so, the pressure of a subset of the system is determined by applying a local (virtual) mechanical deformation, fitting the response to that of a bulk fluid, and then using the YoungLaplace equation to infer the tension of the interface. The accuracy of the method is tested by calculating the local pressure of a series of water simulations at various external pressures. The tension of a simulated curved octanewater interface is computed with the method and compares well with the planar tension (≈ 46.7 dyn/cm). Finally, an ambiguity is resolved between the Harasima and IrvingKirkwood methods of calculating the local pressure as a means for computing the tension.

Generalized KleinKramers equations
View Description Hide DescriptionA generalized KleinKramers equation for a particle interacting with an external field is proposed. The equation generalizes the fractional KleinKramers equation introduced by Barkai and Silbey [J. Phys. Chem. B104, 3866 (2000)10.1021/jp993491m]. Besides, the generalized KleinKramers equation can also recover the integrodifferential KleinKramers equation for continuoustime random walk; this means that it can describe the subdiffusive and superdiffusive regimes in the longtime limit. Moreover, analytic solutions for first two moments both in velocity and displacement (for forcefree case) are obtained, and their dynamic behaviors are investigated.

Markovian dissipative coarse grained molecular dynamics for a simple 2D graphene model
View Description Hide DescriptionBased upon a finiteelement “coarsegrained molecular dynamics” (CGMD) procedure, as applied to a simple atomistic 2D model of graphene, we formulate a new coarsegrained model for graphene mechanics explicitly accounting for dissipative effects. It is shown that, within the Moriprojection operator formalism, the reversible part of the dynamics is equivalent to the finite temperature CGMDequations of motion, and that dissipative contributions to CGMD can also be included within the Mori formalism. The CGMD nodal momenta in the present graphene model display clear nonMarkovian behavior, a property that can be ascribed to the fact that the CGMDweighting function suppresses highfrequency modes more effectively than, e.g., a simple center of mass (COM) based CG procedure. The present coarsegrained graphene model is also shown to reproduce the short time behavior of the momentum correlation functions more accurately than COMvariables and it is less dissipative than COMCG. Finally, we find that, while the intermediate time scale represented directly by the CGMD variables shows a clear nonMarkovian dynamics, the macroscopic dynamics of normal modes can be approximated by a Markovian dissipation, with friction coefficients scaling like the square of the wave vector. This opens the way to the development of a CGMD model capable of describing the correct long time behavior of such macroscopic normal modes.

Hybrid models of molecular machines and the nopumping theorem
View Description Hide DescriptionSynthetic nanoscale complexes capable of mechanical movement are often studied theoretically using discretestate models that involve instantaneous transitions between metastable states. A number of general results have been derived within this framework, including a “nopumping theorem” that restricts the possibility of generating directed motion by the periodic variation of external parameters. Motivated by recent experiments using timeresolved vibrational spectroscopy[Panman et al. , Science328, 1255 (2010)], we introduce a more detailed and realistic class of models in which transitions between metastable states occur by finitetime, diffusive processes rather than sudden jumps. We show that the nopumping theorem remains valid within this framework.

An algorithm for the efficient evaluation of twoelectron repulsion integrals over contracted Gaussiantype basis functions
View Description Hide DescriptionA new algorithm for the evaluation of twoelectron repulsion integrals optimized for high contraction degrees is derived. Both the segmented and general contraction versions of the algorithm show significant theoretical performance gains over the asymptotically fastest algorithms published in the literature so far. A preliminary implementation of the algorithm shows good agreement with the theoretical results and demonstrates substantial average speedups in the evaluation of twoelectron repulsion integrals over commonly used basis sets with varying degrees of contraction with respect to a mature, highly optimized quantum chemical code.

The RARE model: A generalized approach to random relaxation processes in disordered systems
View Description Hide DescriptionThis paper introduces and analyses a general statistical model, termed the RAndom RElaxations (RARE) model, of random relaxation processes in disordered systems. The model considers excitations that are randomly scattered around a reaction center in a general embedding space. The model's input quantities are the spatial scattering statistics of the excitations around the reaction center, and the chemical reaction rates between the excitations and the reaction center as a function of their mutual distance. The framework of the RARE model is versatile and a detailed stochastic analysis of the random relaxation processes is established. Analytic results regarding the duration and the range of the random relaxation processes, as well as the model's thermodynamic limit, are obtained in closed form. In particular, the case of powerlaw inputs, which turn out to yield stretched exponential relaxation patterns and asymptotically Paretian relaxation ranges, is addressed in detail.

Efficient and automatic calculation of optical band shapes and resonance Raman spectra for larger molecules within the independent mode displaced harmonic oscillator model
View Description Hide DescriptionIn this work, an improved method for the efficient automatic simulation of optical band shapes and resonance Raman (rR) intensities within the “independent mode displaced harmonic oscillator” is described. Despite the relative simplicity of this model, it is able to account for the intensity distribution in absorption (ABS), fluorescence, and rR spectra corresponding to strongly dipole allowed electronic transitions with high accuracy. In order to include temperatureinduced effects, we propose a simple extension of the time dependent wavepacket formalism developed by Heller which enables one to derive analytical expressions for the intensities of hot bands in ABS and rR spectra from the dependence of the wavepacket evolution on its initial coordinate. We have also greatly optimized the computational procedures for numerical integration of complicated oscillating integrals. This is important for efficient simulations of higherorder rR spectra and excitation profiles, as well as for the fitting of experimental spectra of large molecules. In particular, the multimode damping mechanism is taken into account for efficient reduction of the upper time limit in the numerical integration. Excited stateenergy gradient as well as excited state geometry optimization calculations are employed in order to determine excited state dimensionless normal coordinate displacements. The gradient techniques are highly costeffective provided that analytical excited state derivatives with respect to nuclear displacements are available. Through comparison with experimental spectra of some representative molecules, we illustrate that the gradient techniques can even outperform the geometry optimization method if the harmonic approximation becomes inadequate.

Particle lifetime in cylindrical cavity with absorbing spot on the wall: Going beyond the narrow escape problem
View Description Hide DescriptionThe mean lifetime of a particle diffusing in a cylindrical cavity with a circular absorbing spot on the cavity wall is studied analytically as a function of the spot radius, its location on the wall, the particle initial position, and the cavity shape determined by its length and radius. When the spot radius tends to zero our formulas for the mean lifetime reduce to the result given by the solution of the narrow escape problem, according to which the mean lifetime is proportional to the ratio of the cavity volume to the spot radius and is independent of the cavity shape, the spot location on the cavity wall, and the particle starting point, assuming that this point is not too close to the spot. When the spot radius is not small enough, the asymptotic narrow escape formula for the mean lifetime fails, and one should use more general formulas derived in the present study. To check the accuracy and to establish the range of applicability of the formulas, we compare our theoretical predictions with the results of Brownian dynamics simulations.

Stochastic model reduction using a modified Hilltype kinetic rate law
View Description Hide DescriptionIn the present work, we address a major challenge facing the modeling of biochemical reaction networks: when using stochastic simulations, the computational load and number of unknown parameters may dramatically increase with system size and complexity. A proposed solution to this challenge is the reduction of models by utilizing nonlinear reaction rate laws in place of a complex multireaction mechanism. This type of model reduction in stochastic systems often fails when applied outside of the context in which it was initially conceived. We hypothesize that the use of nonlinear rate laws fails because a single reaction is inherently Poisson distributed and cannot match higher order statistics. In this study we explore the use of Hilltype rate laws as an approximation for gene regulation, specifically transcription repression. We matched output data for several simple gene networks to determine Hilltype parameters. We show that the models exhibit inaccuracies when placed into a simple feedback repression model. By adding an additional abstract reaction to the models we account for secondorder statistics. This split Hill rate law matches higher order statistics and demonstrates that the new model is able to more accurately describe the mean protein output. Finally, the modified Hill model is shown to be modular and models retain accuracy when placed into a larger multigene network. The work as presented may be used in gene regulatory or cellsignaling networks, where multiple binding events can be captured by Hill kinetics. The added benefit of the proposed splitHill kinetics is the improved accuracy in modeling stochastic effects. We demonstrate these benefits with a few specific reaction network examples

Variable timestepping in the pathwise numerical solution of the chemical Langevin equation
View Description Hide DescriptionStochastic modeling is essential for an accurate description of the biochemical network dynamics at the level of a single cell. Biochemically reacting systems often evolve on multiple timescales, thus their stochastic mathematical models manifest stiffness. Stochastic models which, in addition, are stiff and computationally very challenging, therefore the need for developing effective and accurate numerical methods for approximating their solution. An important stochastic model of wellstirred biochemical systems is the chemical Langevin Equation. The chemical Langevin equation is a system of stochastic differential equation with multidimensional noncommutative noise. This model is valid in the regime of large molecular populations, far from the thermodynamic limit. In this paper, we propose a variable timestepping strategy for the numerical solution of a general chemical Langevin equation, which applies for any level of randomness in the system. Our variable stepsize method allows arbitrary values of the timestep. Numerical results on several models arising in applications show significant improvement in accuracy and efficiency of the proposed adaptive scheme over the existing methods, the strategies based on halving/doubling of the stepsize and the fixed stepsize ones.

Optical studies of random disorder of colloidal photonic crystals and its evolution in evaporation induced selfassembly
View Description Hide DescriptionSelfassembled photonic structures have been under theoretical and experimental study for decades, whereas previous theories on optical properties were mainly concerned with perfect structure or some certain limited kinds of disordered photonic crystals (PCs), making them unsuitable for characterizing the real selfassembled PCs. In order to improve our understanding of the mechanism of selfassembly and provide more crucial clues to further grow perfect crystals, we extended previous widely used scalar wave approximation (SWA), making it be able to characterize longrange disorder (β) and shortrange disorder (α) in PCs synthetically in a simple and effective way. Excellent agreement with in situ observed reflectance of evaporation induced selfassembled colloidal photonic crystals (CPCs) was obtained, demonstrating that the introduction of the parameters α and β in SWA can successfully characterize the disorder in selfassembled CPCs. Furthermore, extended SWA was further used to study the disorder formation in selfassembly, and it was found that during growing stage both β and α drop down, whereas in drying stage β stays nearly unchanged while α increases significantly. It turned out that the growing stage of selfassembly is a stage when the structure transforms from disordered to ordered one, and growth induced disorder mainly arises in drying stage. The results obtained provide an insight into the growth mechanisms of selfassembly and theoretical basis for characterizing optical properties of disordered PCs.

Universal programmable quantum circuit schemes to emulate an operator
View Description Hide DescriptionUnlike fixed designs, programmable circuit designs support an infinite number of operators. The functionality of a programmable circuit can be altered by simply changing the angle values of the rotation gates in the circuit. Here, we present a new quantum circuit design technique resulting in two general programmable circuit schemes. The circuit schemes can be used to simulate any given operator by setting the angle values in the circuit. This provides a fixed circuit design whose angles are determined from the elements of the given matrix–which can be nonunitary–in an efficient way. We also give both the classical and quantum complexity analysis for these circuits and show that the circuits require a few classical computations. For the electronic structure simulation on a quantum computer, one has to perform the following steps: prepare the initial wave function of the system; present the evolution operator U = e ^{−iHt } for a given atomic and molecular Hamiltonian H in terms of quantum gates array and apply the phase estimation algorithm to find the energy eigenvalues. Thus, in the circuit model of quantum computing for quantum chemistry, a crucial step is presenting the evolution operator for the atomic and molecular Hamiltonians in terms of quantum gate arrays. Since the presented circuit designs are independent from the matrix decomposition techniques and the global optimization processes used to find quantum circuits for a given operator, high accuracy simulations can be done for the unitary propagators of molecular Hamiltonians on quantum computers. As an example, we show how to build the circuit design for the hydrogen molecule.

On transition rates in surface hopping
View Description Hide DescriptionTrajectory surface hopping (TSH) is one of the most widely used quantumclassical algorithms for nonadiabaticmolecular dynamics. Despite its empirical effectiveness and popularity, a rigorous derivation of TSH as the classical limit of a combined quantum electronnuclear dynamics is still missing. In this work, we aim to elucidate the theoretical basis for the widely used hopping rules. Naturally, we concentrate thereby on the formal aspects of the TSH. Using a Gaussian wave packet limit, we derive the transition rates governing the hopping process at a simple avoided level crossing. In this derivation, which gives insight into the physics underlying the hopping process, some essential features of the standard TSH algorithm are retrieved, namely (i) nonzero electronic transition rate (“hopping probability”) at avoided crossings; (ii) rescaling of the nuclear velocities to conserve total energy; (iii) electronic transition rates linear in the nonadiabatic coupling vectors. The wellknown LandauZener model is then used for illustration.
 Advanced Experimental Techniques

Cavity ring down spectroscopy with 5 × 10^{−13} cm^{−1} sensitivity
View Description Hide DescriptionThe ultimate sensitivity performances obtained with a continuous wavecavity ring down spectroscopy setup in the near infrared are investigated. At fixed frequency, the noise of the photodetector is found to be the main limitation and the best limit of detection (about 10^{−11} cm^{−1}) is reached after a 10 s averaging. We show that long term baseline fluctuations can be efficiently averaged over several days allowing us to reach a detection limit as low as 5 × 10^{−13} cm^{−1}. The achieved sensitivity is illustrated on narrow spectral intervals where the weakest lines detected so far by absorption spectroscopy are observed: (i) ultraweak transitions of the a ^{1}Δ_{g}(0)−X ^{3}Σ_{g} ^{−}(1) hot band of ^{16}O_{2} near 1.58 μm and (ii) first detection of an electric quadrupole transition in the second overtone band of nitrogen (^{14}N_{2}) near 1.44 μm.

A light scattering study of non equilibrium fluctuations in liquid mixtures to measure the Soret and mass diffusion coefficient
View Description Hide DescriptionWe use dynamic near field scattering to measure the dynamics of concentration non equilibrium fluctuations at the steadystate of Soret separation. The analysis reveals that above a threshold wave vector q _{ c }, the dynamics is governed by diffusion while at smaller wave vectors, gravity dominates. From the measurements, we extract both the mass diffusion and the Soret coefficients. Comparing our results with literature data, we find good agreement confirming that the proposed experimental technique can be considered a sound approach for the study of thermodiffusion processes.
 Atoms, Molecules, and Clusters

Potential energy surface and rovibrational energy levels of the H_{2}CS van der Waals complex
View Description Hide DescriptionOwing to its large dipole, astrophysicists use carbon monosulfide (CS) as a tracer of molecular gas in the interstellar medium, often in regions where H_{2} is the most abundant collider. Predictions of the rovibrational energy levels of the weakly bound complex CSH_{2} (not yet observed) and also of rate coefficients for rotational transitions of CS in collision with H_{2} should help to interpret the observed spectra. This paper deals with the first goal, i.e., the calculation of the rovibrational energy levels. A new fourdimensional intermolecular potential energy surface for the H_{2}CS complex is presented. Ab initiopotential energy calculations were carried out at the coupledcluster level with single and double excitations and a perturbative treatment of triple excitations, using a quadruplezeta basis set and midbond functions. The potential energy surface was obtained by an analytic fit of the ab initio data. The equilibrium structure of the H_{2}CS complex is found to be linear with the carbon pointing toward H_{2} at the intermolecular separation of 8.6 a _{o}. The corresponding well depth is −173 cm^{−1}. The potential was used to calculate the rovibrational energy levels of the paraH_{2}CS and orthoH_{2}CS complexes. The present work provides the first theoretical predictions of these levels. The calculated dissociation energies are found to be 35.9 cm^{−1} and 49.9 cm^{−1}, respectively, for the para and ortho complexes. The second virial coefficient for the H_{2}CS pair has also been calculated for a large range of temperature. These results could be used to assign future experimental spectra and to check the accuracy of the potential energy surface.

A random rotor molecule: Vibrational analysis and molecular dynamics simulations
View Description Hide DescriptionMolecular structures that permit intramolecular rotational motion have the potential to function as molecular rotors. We have employed density functional theory and vibrational frequency analysis to study the characteristic structure and vibrational behavior of the molecule (4^{′},4^{″″}(bicyclo[2,2,2]octane1,4diyldi4,1phenylene)bis2,2^{′}:6^{′},2^{″}terpyridine. IR active vibrational modes were found that favor intramolecular rotation. To demonstrate the rotor behavior of the isolated single molecule, ab initiomolecular dynamics simulations at various temperatures were carried out. This molecular rotor is expected to be thermally triggered via excitation of specific vibrational modes, which implies randomness in its direction of rotation.

Taming the lowlying electronic states of FeH
View Description Hide DescriptionThe lowlying electronic states (X ^{4}Δ, A ^{4}Π, a ^{6}Δ, b ^{6}Π) of the iron monohydride radical, which are especially troublesome for electronic structure theory, have been successfully described using a focal point analysis (FPA) approach that conjoined a correlationconsistent family of basis sets up to augccpwCV5ZDK with highorder coupled cluster theory through hextuple (CCSDTQPH) excitations. Adiabatic excitation energies (T _{0}) and spectroscopic constants (r _{e}, r _{0}, B _{e}, B _{0}, _{e}, ω _{e}, v _{0}, α _{e}, ω _{e} x _{e}) were extrapolated to the valence complete basis set DouglasKroll (DK) augccpwCV∞ZDK CCSDT level of theory, and additional treatments accounted for higherorder valence electron correlation, core correlation, spinorbit coupling, and the diagonal BornOppenheimer correction. The purely ab initio FPA approach yields the following T _{0} results (in eV) for the lowest spinorbit components of each electronic state: 0 (X ^{4}Δ) < 0.132 (A ^{4}Π) < 0.190 (a ^{6}Δ) < 0.444 (b ^{6}Π). The computed anharmonic fundamental vibrational frequencies (v _{0}) for the ^{4,6}Δ electronic states are within 3 cm^{−1} of experiment and provide reliable predictions for the ^{4,6}Π states. With the ccpVDZ basis set, even CCSDTQPH energies give an incorrect ground state of FeH, highlighting the importance of combining highorder electron correlation treatments with robust basis sets when studying transitionmetal radicals. The FPA computations provide D _{0} = 1.86 eV (42.9 kcal mol^{−1}) for the 0 K dissociation energy of FeH and [FeH_{(g)}] = 107.7 kcal mol^{−1} for the enthalpy of formation at room temperature. Despite sizable multireference character in the quartet states, highorder singlereference coupled cluster computations improve the spectroscopic parameters over previous multireference theoretical studies; for example, the X ^{4}Δ → A ^{4}Π and a ^{6}Δ → b ^{6}Π transition energies are reproduced to 0.012 and 0.002 eV, respectively, while the error for the problematic X ^{4}Δ → a ^{6}Δ intercombination excitation is reduced from at least 0.17 eV to about 0.04 eV.