Protein Molecular Dynamics in water

Question

$$\mathcal{V}(\mathbf{r}^N)=\sum_\text{bonds}\frac{k_i}{2}(l_i-l_{i,0})^2 + \sum_\text{angles}\frac{k_i}{2}(\theta_i-\theta_{i,0})^2 + \sum_\text{torsions}\frac{V_n}{2}(1+\cos(n\omega-\gamma) + \sum^N_{i=1}\sum^N_{j=i+1}\left(4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}} \right)^6 \right] + \frac{q_iq_j}{4\pi\varepsilon_0 r_{ij}} \right)$$

The above is potential energy expression for a molecule in vacuum. This part is clear to me!

Now lets put this molecule (in my case protein) in water and run molecular dynamics. Now I have to take solvent effect in my MD. As far as I know, there are two ways:

** SYSTEM = Protein + Water **

1. Change the $ \epsilon _0 =1 $ to $ \epsilon =80 $ (water has dielectric const. of 80) which will reduce my electrostatic force, and that is all the effect the solvent has on the system.

2. Now second way is to model the waters as rigid 3-sites molecule with formal charge on oxygen and hydrogen such as TIP3P etc. Now in this model we don't change the dielectric constant, that is we keep $ \epsilon =1$ while simulating the system. Now while calculating the potential energy of the molecule we have to add extra terms in the above expression, which can be expressed as $$\sum_{i=0}^N \sum _{j=0}^{3M} {q_i q_j} /{\epsilon _0 r(ij) }$$ where $M$ is the number of water molecules. So this additional electrostatic part is energy due to interaction between all sites of every water molecule in the system and each atom in our protein molecule. Please note that $N$ is the total number of atoms in our protein molecule and $r(ij)$ is the distance between $i^{th}$ site and $j^{th}$ site. Now one more additional part is due to Van der Walls interaction between protein and water, which can be expressed in similar manner as above.

One more overhead in the second way is that I also have to calculate potential energy of each water molecule in similar fashion as above. And during the entire process in second way the dielectric constant is set to 1 (i.e. vacuum).

Is my description of working of the second way is correct?? Can someone please shed some light?

Answer 1

Method 2 is essentially correct

Your description of the method 2, which is usually called explicit-solvent simulation, is essentially correct. In explicit-solvent calculations, atoms of water are in principle treated the same as atoms of the protein. As you suggest, the dielectric constant is taken to be 1. In a simplified picture, there are nonbonded forces between all non-excluded pairs of atoms (atoms covalently bonded to each other or covalently bonded to a third shared atom are usually excluded from nonbonded forces). However, for computational efficiency, Lennard-Jones forces for pairs with separation distances larger than some value (e.g. 12 Å), are ignored, and Coulomb interactions are approximated using the PME method or other algorithms that reduced the computational complexity from $O(N^2)$ to something more manageable ($O(N \log N)$ for PME). Recently, using a PME-like algorithm for Lennard-Jones interactions has become more common as well. As an added complication, the common water models (TIP3P, TIP4P variants, OPC, SPC/E) used are typically held rigid (rigid bonds and rigid angles). Popular protein models use rigid bonds with hydrogen atoms as well (CHARMM and Amber).

Method 1 is very crude and not used with atomistic models

The alternative to explicit-solvent simulation is implicit-solvent simulation. However, the method you described of simply using a dielectric constant of 80 is a very poor model of reality for molecular simulations. Consider the hydrogen bonds inside of a protein alpha helix. No water molecules are present inside the helix and there is no room for water molecules between hydrogen bonding groups of the helix (NH and C=O groups); however, in this method, you are reducing the electrostatic energy by a factor of 80 as if these interactions were fully screened by water. The discrete size of water molecules for molecular interactions plays an important role and cannot by represented by a dielectric continuum. Furthermore, entropic effects must be considered as well: water molecules hydrogen-bonded to protein groups are less free to reorient in response to external electric fields. Implicit solvent methods usually include implicit screening from dissolved ions as well. Hence, implicit-solvent methods are more complicated than what you've described. The most common method used today for implicit-solvent molecular dynamics simulations is probably the generalized Born surface-area method, which is implemented in many MD codes such as Amber, NAMD, and OpenMM. Other more accurate methods may be beginning to supersede it, however.