LipidEnrichment and LeafletFinder class possible error

[lunamorrow]

Hi all! I have just started using MDAnalysis as an RA. I have been modifying/implementing a script to calculate the lipid depletion and enrichment around a protein which is only embedded into the cytosolic leaflet. Unfortunately it has been giving me a lot of grief and any help or ideas would be greatly appreciated. Essentially, I run a trajectory through the script, which uses LeafletFinder and LipidEnrichment, then converts the data from LipidEnrichment into a pandas dataframe for csv export, but the data in the csv indicates that at random timesteps in the trajectory all or some of the lipids from one leaflet are being detected on the other leaflet. I have visualised all of these trajectories in VMD and this is certainly not happening! I have also tried just about every processing option and combination possible through gromacs trjconv to varying benefit. The best fit for the script was from creating a larger pbc box, shifting the system by 3 on the z-axis to prevent it crossing the pbc and specifying pbc mol. However, the issue is certainly still there (about 1 in every 10-20 frames), and I cannot present meaningful analysis from this data as these flipflops of the whole membrane is producing insane standard errors. I tried some progressive fits to the protein and these were much.

I know at least one of my peers has heard of this issue but was unsure how/if it was resolved, and I am in contact with someone much more familiar with the MDAnalysis world to give me a hand too. But, has anyone else experienced this phenomena? Aside from making a copy of these classes and modifying them so I can attempt to find the root of this issue, I am out of ideas.

As an aside, I have tested the script with a plain membrane and it performed ok. I have yet to test it with a more 'typical' protein that spans the whole membrane. Also, the protein I am working with induces quite a lot of movement in the membrane, which makes it very curvy, which could be adding to my issues?

[hmacdope]

Hey @lunamorrow, this sounds like an issue with lipyds: https://github.com/lilyminium/lipyds and thus the issue is probably better discussed on the issue tracker there. I will also tag the amazing @lilyminium, the original author of lipyds who may be able to transfer the issue there.

Just quickly here though, in my experience the clustering and graph methods used to assign to the leaflets often need manual tuning, see the options here https://lipyds.readthedocs.io/en/latest/grouping.html#lipyds.leafletfinder.grouping.SpectralClusteringMethod etc. You may need to scan some hyperparameters using a script for example in cutoff, delta and cosine space to get the correct results. You can do this by something like the below:

select = "resname PO4"
deltas = [10,20,30]
cutoffs = [25, 30, 35]
for d in deltas:
    for c in cutoffs # if nesting is getting annoying use itertools.product 
        lf = lipyds.leafletfinder.leafletfinder.LeafletFinder(u, select,  method="spectralclustering", cutoff=c, delta=d)
        # blah blah 
        # visualise your output and decide if its working

@lilyminium may be able to give you more info, hope this helps.

[lunamorrow]

Hey Hugo, fancy seeing you here! Thank you for taking the time to reply.

I have my environment setup from Lily's Github (https://github.com/lilyminium/2023-10-10_lipyd-example/tree/main). I can post on lipyds too, but I have already contacted Lily directly and she is going to look into it for me when she has time.

I will check out the manual and give this fine-tuning a go until I hear back from her, and see if it helps. It'd be great to get this bug sorted!

Update: those adjustments have fixed it!