Computational Modeling and Simulations of Biomolecular Systems

The following sections are included:

• Dedication
• Preface
• Acknowledgments
• Contents

Let us envision a large biological macromolecular system. For example, the system may comprise a protein embedded in a phospholipid bilayer membrane, surrounded by a large number of water molecules and ions. A typical system is schematically illustrated in Figure 1.1. We seek to describe the evolution and behavior of such biological macromolecular systems at the most fundamental level based on the laws of physics. The most realistic and detailed description would involve the quantum mechanical (QM) nature of the microscopic world. This means that the evolution of the atom nuclei and electrons of the biological macromolecular system must be represented through the time-dependent nuclear–electronic wave function Ψ(R₁, …, R_N, r₁, …, r_n, t), which is obtained by solving the time-dependent Schrödinger equation [Merzbacher (1998)],

To provide a complete classical description of an atomic system, we need to specify the position r_i and momenta p_i of all the atoms. But in practice, when we probe a molecular system experimentally, we do not have access to all this information. In fact, even in the case of computer simulations where we have access to all this information, we may still wish to account for the behavior of the system via a reduced set of key variables. Statistical mechanics is a powerful theoretical framework that allows us to achieve a rational management of the limited information about these systems. For this reason, it is the language of choice when trying to communicate about complex molecular systems.

Computer simulations in which a large number of solvent molecules are treated explicitly represent one of the most detailed approaches to study complex biomolecules. Nevertheless, a significant computational cost is associated with the large number of solvent molecules required to model a bulk solution. In practice, a large fraction of computer time is spent calculating the trajectory of the solvent molecules, even though it is primarily the behavior of the solutes that is of interest. For this reason, it is desirable to develop a language and a set of conceptual tools to describe the influence of solvent on biomolecules. One of these conceptual tools is the solvation free energy, which plays an important role in understanding simulations generated with explicit solvent molecules.

We have seen that the charging free energy of a simple monoatomic solute can be approximated by representing the surrounding solvent as a dielectric continuum. It can be shown that this type of treatment leads to the Born model [Born (1920)], which was established long ago. More generally, continuum electrostatics approximations in which the solvent is represented as a featureless dielectric medium are very popular approaches to represent the influence of solvent on biological molecules. An additional factor that greatly affects the solvation forces acting on biomolecules is the influence of salt in the solution.

Issues of molecular recognition, involving the noncovalent association of small ligands to large macromolecules with high affinity and specificity, plays a crucial role in biology and medicinal chemistry. Computational studies can help elucidate the fundamental principles governing those issues at the molecular level. Moreover, improving our ability to screen large databases of compounds in silico to identify potential lead drug molecules with an accurate prediction of binding affinities could have a great impact on structure-based drug design. The main quantity of interest is the standard equilibrium binding free energy, $\text{Δ}{G}_{\text{eq}}^{(\circ )}$. To calculate $\text{Δ}{G}_{\text{eq}}^{(\circ )}$ with computer simulations, a mathematical relation has to be established between macroscopic observables and microscopic variables. The standard binding free energy is defined from the equilibrium constant by $\text{Δ}{G}_{\text{eq}}^{(\circ )}=-k_{\text{B}}T\ln \left(C^{(\circ )}{K}_{\text{b}}\right)$, where C^(°) is standard concentration, k_B is the Boltzmann constant, and T is the absolute temperature.

So far we have mainly considered averages like ⟨x⟩ or ⟨x²⟩, ⟨v⟩ or ⟨v²⟩. Equilibrium properties of complex molecular systems can be analyzed and understood on the basis of the Maxwell–Boltzmann distribution. But different tools and concepts must be introduced to tackle time-dependent kinetic and dynamical properties, where timescales become important. One of the most important concepts in the theoretical treatment of dynamical properties is that of time-correlation functions. For example, the time-correlation function C(t) = ⟨v(t)v(0)⟩, represents the average of the velocity v at time 0 and the velocity at another time t. This type of joint average at two different times reports information about the system that is inherently dynamical. We will see that time-correlation functions play a critical role in the theoretical formulation of dynamical properties.

The most striking drawback of all-atom molecular dynamics (MD) is that the trajectory of a large number of atoms must be calculated, although only a small number among them might be truly relevant to understanding the function that one wishes to investigate. An attractive strategy would be to focus only on the dynamics of the most relevant variables, and yet retain the ability to represent all their microscopic properties (structural and dynamical) correctly. The problem of representing the dynamics of a reduced set of degrees of freedom that is part of a large complex system, is a subject that has been extensively explored in statistical mechanical theories over the last 50 years. At the most rigorous level are the ideas introduced by Zwanzig (1965) and Mori (1965), which consists of “projecting out” uninteresting degrees of freedom, in order to develop an effective but physically correct dynamical scheme representing the time evolution of the most “relevant” degrees of freedom realistically. According to this analysis, the reduced set of coordinates treated explicitly evolves according to an effective nondeterministic dynamics that displays an unpredictable random and stochastic character. In contrast, because it accounts for all the degrees of freedom explicitly, the effective classical MD trajectory of the all-atom system, calculated according to Newton’s classical laws of motion, is completely deterministic. The effective dynamics of reduced models provides important conceptual frameworks for representing a wide range of diffusional and relaxation processes in complex molecular systems.

In Chapter 7, the concept of reduced models comprising a small subset of degrees of freedom was introduced. We saw that these reduced models evolve according to some effective nondeterministic stochastic dynamics to mimic the influence of the degrees of freedom that are left out. The reduced model serves to generate an ensemble of random trajectories that, in turn, can be used to calculate different types of averages. An alternative representation of outcome is through the evolution of a time-dependent probability distribution. In this case, the averaging process is exploited to construct a differential equation to calculate the evolution of the time-dependent probability distribution. While the two representations are formally equivalent, the latter often provides a more convenient route to derive closed-form expressions for important properties of the reduced model. This is the focus of this chapter.

A situation where long timescales are important is when there are rare transitions between two stable states. While the fluctuations within the stable states may occur on a short timescale, the actual transitions between the stable states may take place on a far longer timescale. A typical case is illustrated in Figure 9.1 where two stable states, A and B, are separated by a large energy barrier. An important concept illustrated in Figure 9.1 is the identification of a “reaction coordinate,” x, that can be used to monitor the progress of the transition. Theoretical analysis to determine the transition rate between stable states can vary greatly in sophistication and complexity. One of the simplest treatments goes back to transition state theory (TST).

In Chapter 9, we saw that the long timescale transitions between two stable states take the form of a simple rate dependence. This simplified picture can be generalized to represent the dynamical evolution of a complex system in terms of kinetic models with multiple discrete states. When the transitions between the different states do not depend on the past history, the system is then said to evolve dynamically with no memory. Mathematically, such a kinetic model is called a discrete-state continuous-time Markov chain (named after the Russian mathematician Andrey Andreyevich Markov, 1856–1922). A discrete-state continuous-time Markov chain is a stochastic model in which the probability of the transition from state i to state j depends only on the current state i, and not the previous history of the system. A typical Markov chain is depicted in Figure 10.1. The dynamics of this type of model is governed by the state-to-state transition rates, which are the input parameters. While this simplified representation may look somewhat artificial, it is a physically valid approximation when long-timescale properties are dominated by the rare transitions between a set of long-living metastable states. Furthermore, it is important to recall that the transition rates arise from the well-defined statistical mechanical features of detailed atomic models (see Chapter 9). Therefore, kinetic models with discrete states may be constructed in such a way that remains faithfully representative of the microscopic system of interest with respect to long-timescale properties.

Many of the effective dynamical models that were introduced incorporate elements of stochasticity and randomness. For example, the random walk is a process that evolves through discrete steps in unpredictable directions, the dynamics from the Langevin equation obeys Newton’s equation of motion but with a Gaussian random force, and the dynamics of a kinetic Markov model is a stochastic trajectory that comprises periods of staying in a given state for a random period of time followed by a discrete jump to another state. To simulate these formal models on a computer requires random numbers with specific probability distributions. In everyday life, familiar mechanical devices are used to produce randomness (Figure 11.1). For example, throwing a coin yielding head or tail can be used to generate random numbers with probability one-half. Or throwing a dice with six faces can be used to generate random numbers between 1 and 6. One may combine various combinations of the dice’s faces together to generate random numbers with a probability of 1/6, 2/6, 3/6, 4/6, and 5/6. The latter example shows that it is possible to manipulate the outcome from a given device to generate random numbers with different probabilities. Nevertheless, such simple mechanical devices have obvious limitations when a very large number of repetitive steps must be taken. There are some electronic devices that appear to genuinely be random (e.g., involving particle decay), which could be linked to a computer to simulate stochastic processes. But commonly, it is more convenient to use lists of “pseudo-random” numbers that are generated via some deterministic algorithms. One advantage of pseudo-random numbers is that one can control and repeat the suite of numbers whenever desired, which is useful for debugging and analyzing code performance.

The previous chapters provided a broad overview of the theoretical foundations that underlie many important concepts such as marginal distribution, potential of mean force, free energy, continuum electrostatics, Poisson– Boltzmann theory, ionic screening, binding constants, effective dynamics of reduced models, diffusion processes, transition rates, and Markov chains. These various conceptual elements, can be creatively combined together as “building blocks” to construct models representative of complex macromolecular “machines.” Calling them in this manner is not merely a metaphor: they are very much like real machines, using energy to change their shape, drive motion, and do useful work…

The following sections are included:

• Bibliography
• Further reading
• Index