# Protein Molecular Dynamics in water

$$\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?

• And you have to add other non-contact term like van der Waals forces between protein and water atoms.
– Camps
Oct 20, 2022 at 13:14

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).