Issues with python plotting graphs

[cgseitz]

Hello,

After running chap and selecting water to analyze, I created the free energy profile, hydrophobicity profile, radius profile and solvent number density profile graphs from here:
https://www.channotation.org/docs/plotting_python/

The free energy profile and radius profile graphs look realistic.

1. The hydrophobicity profile graph is a straight line at 0. There are no dots on the graph to show the hydrophobicity profile. Do you know why the hydrophobicities are not being read?

2. The solvent number density profile graph for water contains consistent density around ~5-10 nm^-3. However, the pore being analyzed is huge (3-4 nm wide) so I would assume there is a free passage of water through the pore, and that the density should be closer to free water density (~33 nm^-3). Do you know why the solvent density would be lower than free water? Is this realistic?

Thanks for your help!

[Inniag]

I can think of two reasons for why hydrophobicity is not computed properly:

• Your pore does not consist of amino acids and there is not hydrophobicity information available in CHAP. You will need to provide your own scale as documented here: http://www.channotation.org/docs/hydrophobicity_scales/
• For larger than typical pores, you might need to adjust the parameters used for identifying pore-lining residues. In particular, check -pm-pl-margin and -pm-pf-sel.

As to whether or not a below bulk density of water is sensible, I could not tell without having much more information on the system you study. It really depends on your use case, the simulation data you have, and the underlying physical system.

[cgseitz]

I tried setting the -pm-pl-margin flag, at values between 0.5 and 5 it gives this error:

Standard library runtime error (possible bug):
(exception type: St13runtime_error)
Could not solve tridiagonal system in 1D interpolation. LAPACK error code: -7

This system consists of pores inside the proteins, and pores in between the proteins. For a short simulation looking at the pore inside one protein, I can generate hydrophobicity data. I also wanted to examine a pore in between proteins. I sliced my trajectory so the center of mass would be on the inter-protein pore, but with short trajectories I cannot generate hydrophobicity data.

The system consists of proteins in a waterbox. The proteins are all standard amino acids, and the waterbox is water and ions. There are no large conformational changes happening in the simulation that would affect the pore size. The simulations are all atom, normal MD. What other information could help determine what is a reasonable density of water in this system?

[Inniag]

The error above suggests that Chap can not interpolate the hydrophobicity values at the individual residue positions into a continuous profile. The reason for this is hard to diagnose, it could be that no pore-lining residues are found due to the choice of the margin parameter.

Generally speaking, I would suggest to start with just a single frame to identify good parameters. Then you can start going through trajectories to see if maybe some frames are particularly problematic. If need be, you can skip some pathological frames, so long as you have enough frames for statistical sampling. For hydrophobicity in particular, I would not expect much variation over time, if, as you say, there are no large conformational changes.

Regarding water density: I assume visual inspection in e.g. VMD confirms the pore is water-filled? For radii > 1nm I would certainly expect that. One potential pitfall here is that CHAP does not pick up the waters as belonging to the right group of particles. It uses the libgromacs library for particle group parsing, which in my experience works best when sticking to Gromacs naming conventions. If you are using trajectories from other MD engines or use different residue names from what Gromacs expects, this might throw libgromacs off. Water molecules in particular can be named variously, according to force field and setup pipeline, as e.g. WAT, SOL, TIP, HOH etc. Can you confirm that CHAP picks up solvent particles? You might want to look at the -sel-solvent flag as well.

[cgseitz]

Thanks for your help! After a lot of work, a coworker discovered the likely problem. For my system CHAP seemed to find multiple pores, and would crash when it tried to get results from distinct pores at the same time. I fixed this by manipulating my trajectories so that each one contained a single protein pore, but all of the solvent present in the system. Adjusting the -max-free-dist and -pf-cutoff only really helped once I eliminated the extra pores. This fixed most of the issues. However, I found something quite odd and unfortunate. If I added -sel-solvent 'Ions' to my CHAP run command, this would select all of the available solvent (water and ions). However, if I added -sel-solvent alone to my CHAP run command, then it would prompt me to select what I wanted the solvent to be. I can select ions, and then it properly selects ions. This was very confusing and I didn't notice it for a while, but should be rectified for future users.