Calculation of true water permeation through lipid bilayer membrane

[Eisha Khilji]

Hello I want to calculate only those water molecules which must pass through the protein (in z coordinate direction). I have a MD system of a lipid bilayer membrane with a water channel protein embedded in it. I am using the following code for analysis, but i am unable to add more than one residues in atomselection command. Is there any way out by using this? 

Or if anyone can suggest something more helpful to solve the problem.

import numpy as np
import MDAnalysis as mda

u = mda.Universe("step5_input.psf","step7_production.dcd")

water_oxygens = u.select_atoms("around 3 protein and (resid 163 and resid 187 and resid 24)", updating=True)

fw=open("watermolecules.txt","w")

for ts in u.trajectory:
    i=0
    while i <len(water_oxygens):
        if "TIP3" in str(water_oxygens[i]):
            frame=str(ts.frame)
            tip=str(water_oxygens[i]).split()[1][:-1]
            fw.write("TIP3: {0} at Frame {1}\n".format(tip, frame))

        i+=1
fw.close()

[Hugo Macdermott-Opeskin]

Hey Eisha, 

I think you might be using AND  where you mean OR 

For example, u.select_atoms("resid 163 and resid 164") selects atoms that belong to BOTH  resid 163 and 164. To get atoms that belong to either try and use the OR keyword. 

To simplify your workflow a bit, try using something like the following selection string "around 3 protein and resname TIP3P and name O" to select water oxygens within 3 of the protein. 

Let me know if you ahve any other questions. 

Cheers

Hugo

[Hugo Macdermott-Opeskin]

Sorry accidentally sent a private mesasge as well, copying here:

Hugo: Secondly, while I am no expert on your specific research problem, you may need some additional complexity to capture water molecules that pass through a protein. An approach I have used it is to define two zones, and track whether water molecules flip betweent he two zones without  passing through the periodic boundary. 

Eisha: Yes you just mentioned "define two zones, and track whether water molecules flip between the two zones without  passing through the periodic boundary" .. This is the real thing I need, but I do not know what to do. I have used cyzone , and prop z , but all in vain, they also calculates false events ( molecules whose coordinates are changing by passing through periodic boundary condition). I would be very thankful if you can help me with this.

if i remember correctly from when I did this in the past, I used four zones, the  upper and lower parts of a cylinder, and zones for solvent. I am happy to share my code with you as this took a long time for me to get right. One thing to note is that you should try and save your trajectory frequently for this to work optimally. 

See this piece of analysis code https://github.com/hmacdope/TmEnc_Flux/blob/main/Analysis/pore_flux_2x_cl.py, and corresponding discussion in the paper here

https://www.science.org/doi/10.1126/sciadv.abl7346. If you don't have SciAdv access I can send you a copy, just let me know. 

Cheers

Hugo

[Eisha Khilji]

Thanks I got the paper and the script as well, I will let you know if I have more queries, after reading it. 
Thanking in advance.

--
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/QUk_fcS2aoA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/580a5a30-4a93-49d6-b863-72003763149en%40googlegroups.com.

[Eisha Khilji]

 I have some queries regarding the code. As I am new to python so I request you to help me. I mentioned the queries below:

1. u = mda.Universe("no_water.M29_ions_1264.prmtop", "M29_ions_1264_combined_imaged_full_no_water.nc")

 It requires only  two files as input, right? a .prmtop file and .nc file.  I run simulations on NAMD with input files generated through CHARMM force field and I have files with different formats, so here I have to use a .psf file (contains coordinates information), and a trajectory file (.dcd).. please correct me if I am wrong. 

2.     def check_in_cylinder_upper_or_lower(pos_array, cylinder_center, cylinder_rad, cylinder_height): Here, you defined a function,, what do you mean by pos array?

3. center_pore_ag_1 = u.select_atoms("(resid 25 97 506 99) and name CA")

print(center_pore_ag_1.resnames) 

center_pore_ag_2 = u.select_atoms("(resid 796 868 870 249) and name CA")

print(center_pore_ag_2.resnames)

center_pore_ag_3 = u.select_atoms("(resid 1125 1127 1020 1053) and name CA")

print(center_pore_ag_3.resnames) .......................

In this part of your code, which type of residues you mentioned here? If they are the residues present in pore lining? 

I am working on a water channel protein, named aquaporin. It is a tetramer, having 4 water channel pores. Depending on 4 pores, I should mention the pore lining residues of all four pores here?

4.

# backup constants
	cyl_rad = 15 #angstrom
	cyl_height = 10 #angstrom
	twiddle_z_pore = -0.0 #shift pore down slightly

Here, it build one cylinder as a whole, or it builds cylinder for each pore?

I am attaching the images for reference. One is a view with protein and water box (I hide the lipids), and the other picture is from the above of protein, showing four pores.

Thanks.

Respects,

[Yogesh Sharma]

It goes quite tricky to calculate permeation from each monomer. The script seems to be working for a single channel. You can try simulating aqp as monomer if that helps. 

You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAGPgGqkhYxJ%3DPH5C-qRTr7-eRcsm6SZoGv_Biq7b2TZs2mkhkQ%40mail.gmail.com.

[Hugo Macdermott-Opeskin]

Hi Eisha,

1. >> Your NAMD and CHARMM simulations will work! This is one of the great things about MDAnalysis as your analysis scripts are largely independent of your file format.

Just replace the names with those of your files.

2. >> pos_array is the numpy array of ion positions returned by ion_ag.positions It will have a shape of (n_ions, 3).

VERY IMPORTANT POINT relevant to all below points: The code is validated for my system only,  it is unlikely to work for your system out of the box without some modification. It is really really important to be skeptical when doing science and to try and validate any workflow for your system. Never trust a tool straight out of the box just because you want it to work.  For example there are a number of hardcoded constants relevant to my system only such as the cylinder height etc etc.

3. >> Yes you are right, these are residues that line the pore in my system that has only one pore

3.I "I am working on a water channel protein, named aquaporin. It is a tetramer, having 4 water channel pores. Depending on 4 pores, I should mention the pore lining residues of all four pores here?" 

>> The setup in my script will only work for a single pore. As Yogesh says you would need to modify it to work for all four pores as once. 

alternately you could work on one pore at a time and collate the results later. 

4. >> It only builds a single cylinder for ONE  pore. You would have to repeat this process or modify the script significantly to get multiple pores.

Additional note: Are the aquaporin pores formed between monomers? If so simulating as a monomer is not going to work.

Cheers

Hugo

[Eisha Khilji]

Yes, I am modifying this script for my system. The protein has its pore channels inside the monomers, so yes I can trim each monomer and it will work as an individual water channel.

But, I'm still confused at this point. Following is a portion from your script.

center_pore_ag_1 = u.select_atoms("(resid 25 97 506 99) and name CA")

print(center_pore_ag_1.resnames)

center_pore_ag_2 = u.select_atoms("(resid 796 868 870 249) and name CA")

print(center_pore_ag_2.resnames)

center_pore_ag_3 = u.select_atoms("(resid 1125 1127 1020 1053) and name CA")

print(center_pore_ag_3.resnames)

center_pore_ag_4 = u.select_atoms("(resid 611 613 1277 1193) and name CA")

print(center_pore_ag_4.resnames)

center_pore_ag_5 = u.select_atoms("(resid 354 356 763 282) and name CA")

print(center_pore_ag_5.resnames)

1.       If I edit it in the way mentioned below, it builds a single cylinder around all of these residues, OR it takes the residues as a centre and builds each cylinder around each chain separately (4 cylinder) ???

2. If it only builds a single cylinder then how should I edit this part to work it on my system ? I mean what type of residues I can store here in pore_ag_1 and 2 and 3 and 4 and 5, Should I mention the pore lining residues for one single channel and separate them differently in pore_ag_1 and pore_ag_2 and so on ?

center_pore_ag_1 = u.select_atoms("index 346 2301 2690")     #pore lining residues for 1st pore present inside chainA
print(center_pore_ag_1.resnames)
center_pore_ag_2 = u.select_atoms("index 3736 5691 6080")    #pore lining residues for 2nd pore present inside chainB
print(center_pore_ag_2.resnames)
center_pore_ag_3 = u.select_atoms("index 7126 9081 9470")    #pore lining residues for 3rd pore present inside chainC
print(center_pore_ag_3.resnames)
center_pore_ag_4 = u.select_atoms("index 10516 12471 12860")    #pore lining residues for 4th pore present inside chainD
print(center_pore_ag_4.resnames)

( I use index values for CA of residues instead of resnames in order to avoid the confusion, because chains are identical and have same residues at same positions)

[Eisha Khilji]

I got these results:

## RESULTS ##

transitions with oscillations 11369
transitions discounting oscillations 585
total pore interaction events 68141

transitions with oscillations 8390
transitions discounting oscillations 440
total pore interaction events 70039

transitions with oscillations 8599
transitions discounting oscillations 419
total pore interaction events 66401

transitions with oscillations 10998
transitions discounting oscillations 372
total pore interaction events 70823

[[11369   585 68141]
 [ 8390   440 70039]
 [ 8599   419 66401]

 [10998   372 70823]]

But I also need residue ids of water molecules. What should I do?

[Hugo Macdermott-Opeskin]

Hi Eisha,

Sorry for the confusion, re-reading the script, it does indeed create 5 cylinders for 5 separate pores, which i think is what you want. The version for a single pore is at the same github link.

You will need to visualise your cylinder to see if it is correct. The `cylinder` command in VMD is likely your best option. See: https://www.ks.uiuc.edu/Research/vmd/vmd-1.7.1/ug/node169.html 

The center of the pore is printed by the line:

print(f"Pore COM is {[pore_ag_com[0], pore_ag_com[1], pore_ag_com[2] + twiddle_z_pore]}")

To get the resids, you will need to modify the script yourself.

Cheers

Hugo

[Eisha Khilji]

Thank you for your response. It's all clear now. I have just one last question. 

What is the main difference between transitions with oscillations and

transitions discounting oscillations. 

--

You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/QUk_fcS2aoA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/95d845f1-bf9f-4d40-9f22-49827e6c4486n%40googlegroups.com.

[Hugo Macdermott-Opeskin]

Hi Eisha,

Transitions with oscillations count flipping between the upper and lower portion of the cylinder as an event, eg movement of water molecules back and forth in the pore. This is not as useful. 

Transitions without oscillations only count an event once a flip has occured and the ion has left the cylinder without flipping back. So therefore there will be fewer of them. 

Just one thing, remember to cite MDAnalysis and the paper that the script comes from :)

[Eisha Khilji]

I will definitely cite this when my work go for publishing. And thank you for this much cooperation and guidance. I am really grateful.

Thanks,

[Hugo Macdermott-Opeskin]

No worries, let me know if you have any more questions. :)

[Eisha Khilji]

Yes, thank you.

--

You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/QUk_fcS2aoA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/f164d095-1759-47fb-b8e9-e6f9418c63ecn%40googlegroups.com.