10.7: Simulations

Introduction

Computer simulations are often used to understand and predict the behavior and structure of molecules. With modern computers it is practical to perform simulations involving  macromolecules. This section provides brief descriptions of molecular dynamics simulations (MD) for predicting behavior and machine learning (ML/AI) for the prediction of macromolecular structures. 

Molecular dynamics

Classical molecular dynamics (MD) simulations attempt to track the time evolution of the positions and velocities of the atoms and thereby the molecules in a system. In theory with the exact interaction potentials and adequate computing power, one could simulate the motion of all the species involved in a chemical reaction, diffusion across a membrane, folding of a protein, etc. These simulations essentially calculate the physical trajectory of all the simulated particles versus time. This is essentially the same problem as numerical simulation of chemical kinetics except for the scale of the problem. The position vector \(\vec{q}\) and rate of change of \(\vec{q}\) along the three coordinates (velocity) must be kept track of and updated at each time step. This requires being able to calculate the force on each particle based on the position of all the other particles. This is done by taking the partial derivative of the potential V (e.g. \(\partial{V}/dx\), \(\partial{V}/dx\), \(\partial{V}/dx\), hereafter represented as \(\partial{V}/\partial {q}\)). So at each time step one calculates a new position \(q\) and velocity (dq/dt) using Newton's equations:

\[q(t+\delta t) = q_0 + \left(\dfrac{dq}{\delta t}\right)_0 dt\]

\[\dfrac{dq}{dt} (t+\delta t) = \left(\dfrac{dq}{dt}\right)_0 - \delta t \left[\left( \dfrac{∂V}{∂q} \right)_0\dfrac{1}{m_q}\right].\]

Here \(m_q\) is the mass of the particle. Just as with numerical integration of chemical concentrations versus time there are sophisticated algorithms for efficiently calculating the values after the next time step.

Force fields (the potential)

Let’s interrupt our discussion of MD propagation of coordinates and velocities to examine the ingredients that usually appear in the force fields mentioned above. In Figure \(\PageIndex{1}\), we see a molecule in which various intramolecular and intermolecular interactions are introduced.

Figure \(\PageIndex{1}\): Depiction of a molecule in which bond-stretching, bond-bending, intramolecular van der Waals, and intermolecular solvation potentials are illustrated.

The total potential of a system containing one or more such molecules in the presence of a solvent (e.g., water) it typically written as a sum of intramolecular potentials (one for each molecule in the system) and intermolecular potentials. The total potential experienced by a particle thus has the following generic form:

\[V = V_{stretch} + V_{bend} + V_{torsion} + V_{Coulomb} + V_{vanderWaals}\]

The covalent interactions describing how the energy varies with bond stretching, bond bending, and dihedral angle distortion are depicted are:

\[V_{stretch} = \sum_{bonds}{(1/2)k_{f,stretch}(R-R_e)^2}\]

\[V_{bend} = \sum_{angles}{(1/2)k_{f,bend}(\theta-\theta_e)^2}\]

\[V_{torsion} = \sum_{dihedrals}{\sum_\eta{\frac{V_\eta}{2}\left(1 + cos(\eta\phi - \gamma)\right)}}\]

where \(\eta\) is the periodicity of the dihedral rotation (e.g 2 = cis-trans), \(V_\eta\) is the barrier height for torsion, \(\gamma\) is the phase.

The non-covalent interactions describing electrostatic and van der Waals interactions among the atoms in the molecule are

\[V_{\rm noncovalent}=\sum_{i<j}^{\rm atoms} \left[\dfrac{A_{i,j}}{r_{i,j}^{12}}-\dfrac{B_{i,j}}{r_{i,j}^{6}}+\dfrac{q_iq_j}{\varepsilon r_{i,j}}\right].\]

These functional forms would be used to describe how the energy \(V(q)\) changes with the bond lengths (\(r\)) and angles (\(\theta,\phi\)) within, for example, each of the molecules shown in figure \(\PageIndex{1}\) (let’s call them solute molecules) as well as for any water molecules that may be present (if these molecules are explicitly included in the MD simulation).

The interactions among the solute and solvent molecules are also often expressed in a form involving electrostatic and van der Waals interactions between pairs of atoms- one on one molecule (solute or solvent) and the other on another molecule (solute or solvent). Usually hydrogen bonding is included in this as a primarily electrostatic interaction.

\[V_{\rm intermolecular} = \sum_{i<j}^{\rm atoms} \left[\dfrac{A_{i,j}}{r_{i,j}^{12}}-\dfrac{B_{i,j}}{r_{i,j}^{6}}+\dfrac{q_iq_j}{\varepsilon r_{i,j}}\right].\]

The Cartesian forces on any atom within a solute or solvent molecule are then computed for use in the MD simulation by using the chain rule to relate derivatives with respect to Cartesian coordinates to derivatives of the above intramolecular and intermolecular potentials with respect to the interatomic distances and the angles appearing in them.

Because water is such a ubiquitous component in condensed-phase chemistry, much effort has been devoted to generating highly accurate intermolecular potentials to describe the interactions among water molecules. In the popular TIP3P and TIP4P models, the water-water interaction is given by

\[V = \dfrac{A}{r_{OO}^{12}}-\dfrac{B}{r_{OO}^{6}}+\sum_{i,j}\dfrac{kq_iq_j}{r_{i,j}}.\]

where r_OO is the distance between the oxygen atoms of the two water molecules in Å, and indices \(i\) and \(j\) run over 3 or 4 sites, respectively, for TIP3P or TIP4P, with \(i\) labeling sites on one water molecule and \(j\) labeling sites on the second water molecule. The parameter \(k\) is 332.1 Å kcal mol^-1. A and B are conventional Lennard-Jones parameters for oxygen atoms and q_i is the magnitude of the partial charge on the ith site. In Figure \(\PageIndex{2}\), we show how the 3 or 4 sites are defined for these two models.

Figure \(\PageIndex{2}\): Location of the 3 or 4 sites used in the TIP3P and TIP4P models.

Typical values for the parameters are given in the table below.

	rOH(Å)	HOH angle degrees	rOM(Å)	A (Å^12kcal/mol)	B (Å^6kcal/mol)	qOor qM	qH
TIP3P	0.9572	104.52		582 x10^3	595	-0.834	0.417
TIP4P	0.9672	104.52	0.15	600 x10^3	610	-1.04	0.52

In the TIP3P model, the three sites reside on the oxygen and two hydrogen centers. For TIP4P, the fourth site is called the M-site and it resides off the oxygen center a distance of 0.15 along the bisector of the two O-H bonds as shown in Figure \(\PageIndex{2}\). In using either the TIP3P or TIP4P model, the intramolecular bond lengths and angles are often constrained to remain fixed; when doing so, one is said to be using a rigid water model.

There are variants to these two 3-site and 4-site models that, for example, include van der Waals interactions between \(H\) atoms on different water molecules, and there are models including more than 4 sites, and models that allow for the polarization of each water molecule induced by the dipole fields (as represented by the partial charges) of the other water molecules and of solute molecules. The more detail and complexity one introduces, the more computational effort is needed to perform MD simulations. In particular, water molecules that allow for polarization are considerably more computationally demanding because they often involve solving self-consistently for the polarization of each molecule by the charge and dipole potentials of all the other molecules, with each dipole potential including both the permanent and induced dipoles of that molecule. This web page provides links to numerous software packages that use these kinds of force fields to carry out MD simulations. These links also offer more detailed information about the performance of various force fields as well as giving values for the parameters used in those force fields.

The parameter values are usually obtained by

1. fitting the intramolecular or intermolecular functional form (e.g., as shown above) to energies obtained in electronic structure calculations at a large number of geometries, or
2. adjusting them to cause MD simulations employing the force field to reproduce certain thermodynamic properties (e.g., radial distribution functions, solvation energies, vaporization energies, diffusion constants), or some combination of both. It is important to observe that the kind of force fields discussed above have limitations beyond issues of accuracy. In particular, they are not designed to allow for bond breaking and bond forming, and they represent the Born-Oppenheimer energy of one (most often the ground) electronic state. There are force fields explicitly designed to include chemical bonding changes, but most MD packages do not include them. When one is interested in treating a problem that involves transitions from one electronic state to another (e.g., in spectroscopy or when the system undergoes a surface hop near a conical intersection), it is most common to use a combined QM-MM approach. A QM treatment of the portion of the system that undergoes the electronic transition is combined with a force-field (MM) treatment of the rest of the system to carry out the MD simulation.

Practical considerations

By applying one of the time-propagation algorithms to all of the coordinates and momenta of the \(N\) molecules at time t, one generates a set of new coordinates \(q(t+\delta t)\) and new velocities \(dq/dt(t+\delta t)\) appropriate to the system at time \(t+dt\). Using these new coordinates and momenta as \(q_0\) and \((dq/\delta t)_0\) and evaluating the forces \(–(∂V/∂q)_0\) at these new coordinates, one can again use the propagation equations to generate another finite-time-step set of new coordinates and velocities. Through the sequential application of this process, one generates a sequence of coordinates and velocities that simulate the system’s behavior. 

These MD simulations are difficult to carry out in a parallel manner. One can certainly execute many different classical trajectories on many different computer nodes, which when averaged can provide values for average (thermodynamic) outcomes. However, to distribute one trajectory over many nodes is difficult. The primary difficulty is that, for each time step, all \(N\) of the molecules undergo movement to new coordinates and momenta. To compute the forces on all \(N\) molecules requires on the order of \(N^2\) calculations (e.g., when pairwise additive potentials are used).

Another factor that complicates MD simulations has to do with the wide range of times scales that may be involved. For example, for one to use a time step dt short enough to follow high-frequency motions (e.g., O-H stretching) in a simulation of an ion or polymer in water solvent, dt must be of the order of 10^-15 s. To then simulate the diffusion of an ion or the folding of a polymer in the liquid state, which might require 10^-4 s or longer, one would have to carry out 10¹¹ MD steps. In the table below we illustrate the wide range of time scales that characterize various events that one might want to simulate using some form of MD, and we give a sense of what was practical using MD simulations in the year 2010. Since that time things have not improved much. This may seem counter to the massive improvement in computer usability, but because these calculations do not parallelize well the explosion in number of cpu cores with only an increase in cpu speed of 2 to 3 has not helped much.

Examples of dynamical processes taking place over timescales ranging from 10^-15 s through hundreds of seconds, each of which one may wish to simulate using MD.

10^-15-10^-14 s	10^-12s	10^-9s	10^-6s	10^-3 s	110 s
C-H, N-H, O-H bond vibration	Rotation of small molecule	Routinely accessible time duration for atomistic MD simulation	Time duration for heroic atomistic MD simulation	Time duration achievable using coarse-graining techniques	Time needed for protein folding

Because one can not afford to carry out simulations covering 10^-3 -100 s using time steps needed to follow bond vibrations 10^-15 s, it is necessary to devise strategies to focus on motions whose time frame is of primary interest while ignoring or approximating faster motions. For example, when carrying out long-time MD simulations, one can ignore the high-frequency intramolecular motions by simply not including these coordinates and momenta in the Netwonian dynamics (e.g., as one does when using a rigid-water model discussed earlier). In other words, one simply freezes certain bond lengths and angles. Of course, this is an approximation whose consequences must be tested and justified, and would certainly not be a wise step to take if those coordinates played a key role in the dynamical process being simulated. Another approach, called coarse graining involves replacing the fully atomistic description of selected components of the system by a much-simplified description involving significantly fewer spatial coordinates and momenta.

Quite impressive simulations can be made using some approximations. For example see this simulation of a COVID virus particle in a water drop.

Machine learning (AI) prediction 3D protein structure

Based on: https://bio.libretexts.org/Bookshelves/Biochemistry/Fundamentals_of_Biochemistry_(Jakubowski_and_Flatt)/01%3A_Unit_I-_Structure_and_Catalysis/04%3A_The_Three-Dimensional_Structure_of_Proteins/4.14%3A_Predicting_Structure_from_Sequence_and_Sequence_from_Structure_Function_(New_10%2F%2F24)#The_Protein_Folding_Problem__Sequence_to_3D_Structure

Protein structure can be experimentally analyzed and the 3D structure of a protein determined using NMR, X-ray crystallography, and cryo-EM. Now, using sequence and structural databases (like the PDB with over 227,000 structures), we can often predict the 3D structure of a protein just from its linear sequence by comparing the sequence of a protein of unknown tertiary structure to homologous proteins (by sequence) whose 3D structures are known. Machine learning and artificial intelligence have extended earlier and simpler "homology modeling" attempts to allow structural predictions for millions of protein sequences using programs such as RoseTTAFold and AlphaFold. RoseTTAFold and AlphaFold produce high-quality structure predictions when trained using the vast sequence information in the Protein Data Bank.  Embedded in those linear sequences is a large amount of hidden (to the human eye) evolutionary information that machine learning and AI can harness to predict 3D structures. The key appears to be that particular sequences lead to the same sequential steps with similar balances between the kinetics and thermodynamics as the proteins fold. These AI tools work less well when very limited sequence comparisons are available.

In general (for smaller proteins), the protein folding problem, the prediction of structure from sequence, appears to have been "solved". The Nobel Prize in Chemistry in 2024 was awarded to Demis Hassabis and John M. Jumper from Google DeepMind for developing AlphaFold and David Baker for developing RoseTTaFold and other powerful techniques described below.   (For more details on how this works see, Chapter 4.13: Predicting Structure and Function of Biomolecules Through Natural Language Processing Tools)

A comparison of protein structures obtained using these programs with known 3D structures obtained through X-ray crystallography or other techniques shows them almost identical. Different metrics can be used to compare predicted structures to the actual ones. The root mean squared deviation (RMSD) is a common one. RoseTTAFold uses a TM-score to assess the topological similarity of protein structures. Compared to RMSD, the TM-score weights smaller distance errors higher than larger ones, making it sensitive to the global fold, not local structural differences. TM values range from 0-100 (100 is a perfect match). Scores below 17 indicate no topology match, while those above 50 suggest a common fold.

AlphaFold uses a "neural network, meaning it simultaneously considers patterns in protein sequences, how a protein’s amino acids interact, and its possible three-dimensional structure. In this architecture, one-, two-, and three-dimensional information flows back and forth, allowing the network to collectively reason about the relationship between a protein’s chemical parts and its folded structure". This is essentially a statistical probability engine that predicts a structure based on what it is most like. Programs of this type might allow the generation of proteins with new therapeutic or commercial potential based on sequences. These include vaccines, sensors, specific immune system suppressors or activators, and antivirals. AlphaFold has now been used to predict the structure of 214 million proteins from more than one million species — essentially all known protein-coding sequences. These AlphaFold iCn3D models are available in the AlphaFold Database - Protein Structure Database.

Figure \(\PageIndex{3}\) shows the backbone tube cartoon of the x-ray pdb structure of the small protein (1xww, cyan) and the structure predicted by both RoseTTAFold program and AlphaFold (magenta) just from its primary sequence. Sulfate, a competitive inhibitor, is shown (spacefill) bound in the active site. The alignment is spectacular, except for the N-terminal 5 amino acids at the bottom of the figure (6 o'clock). This stretch has more disorder even in the x-ray structure as the amino acids have high B-factors, indicating more conformational flexibility. 

Figure \(\PageIndex{3}\): Comparison of the x-ray and computationally predicted structures of human low molecule weight protein tyrosine phosphatase

Left panel: X-ray structure (cyan) of low molecular weight protein tyrosine phosphatase with bound SO₄^2- (1xww) and corresponding structure predicted by the RoseTTAFold (magenta).  Right panel: Same structures using AlphaFold for the structural prediction.