Water inside protein

[jacopo382]

I need to calculate water inside an aquaporin during an md simulation . What can i do to determe the number of water in protein using MD analisis?

[orbeckst]

You can look into the waterdynamics MDAKit.

You can use some of the geometric selections in MDAnalysis

• cyzone: select all water molecules inside a cylinder
• cylayer: select all waters inside a cylindrical sphere

With dynamics selections you can get the changing number of waters:

pore_waters = u.select_atoms("name OW and cyzone 15 4 -8 protein and resid 42", updating=True)
n_pore_waters = []
for ts in u.trajectory:
   n_pore_waters.append([ts.frame, ts.time, pore_waters.n_atoms])

# make data structure a numpy array
n_pore_waters = np.array(n_pore_waters)

# plot with matplotlib
plt.plot(n_pore_waters[:, 1], n_pore_waters[:, 2])
plt.xlabel("time (ps)")
plt.ylabel("number of waters in pore")

This code selects all "OW" atoms (water oxygen — you may need a different selection) that are also inside a cylinder that is defined as follows with cyzone: the center of geometry of protein and resid 42, and creates a cylinder of external radius 15 Å centered on the center of geometry (COG). In z, the cylinder extends from 4 Å above the COG to 8 Å below. Positive values for zMin, or negative ones for zMax, are allowed. This defines your Aqp's pore region. We then iterate over the trajectory and count the number in the atomgroup pore_waters and store them in a list.

For post processing, I turn the list into a numpy array and then do a simple plot of the timeseries of the number of pore waters.