[Solvent] Accessible Surface Area in LOOS

[acnash]

Good morning LOOS members,

Thank you for your earlier support. I decided this needed a new discussion despite all my questions being part of the same research project. Before asking I spent several hours reading through the class description and each of the prebuilt tools, however, I wasn't able to find any code immediately related to calculating a protein's solvent accessible surface area. Did I miss something? There are several packages out there that can do this, however, I've put in considerable effort building the tools I need so far using LOOS that I would like to keep my code as succinct as possible.

Many thanks.

[agrossfield]

Ah, you've touched a bit of a sore point. :-) We chose not to implement solvent accessible surface area, because doing it well is quite hard (I don't know if it's still true, but for years you got wildly different answers depending on what package you used).

In general, we've found there are generally other ways to get at what you're looking for that aren't as brittle. For example, for the problem you've been discussing you could easily count contacts to water, lipid, and peptide to see how occluded the protein structures are. The tool fcontacts can do this, for example, if you treat your protein as the probe and lipid and water as the targets.

[acnash]

Yep! A very good point. When I couldn't find a SASA class/etc., I looked for an algorithmic breakdown so I could implement it myself. As you said, I found several methods but implementing one is well beyond the time scope of this project. That's when I found myself going through MDAnalysis (no luck, they're proposing an implementation this year as part of a Google Code internship) and MDTraj, and it was at this point that I felt like I was fragmenting my work across too many environments and APIs. On that note, I'm really glad that you've got a workaround. I'll look into the contact tool. Many thanks.

[agrossfield]

If you need to implement something yourself (eg because you want to subset your data by oligomerization state), I'd look at the hardContact() and logisticContact() methods to AtomicGroup. They're both good for counting contacts between groups (I'd recommend using just heavy atoms to make it faster), with the only difference being the function used (logisticContact gives a continuously varying value, while hardContact is either 0 or 1 for each pair of atoms).

[acnash]

Brilliant. Just what I was after. Thank you.

[acnash]

Hi @agrossfield. Just an update on logisticContact() and I'm hoping I could get your thoughts. I'm sorry for the length, questions are towards the end, the rest is context and methodology.

As you know I'm after a SASA function to calculate the surface area of TM protein domains (helical) in a model bilayer (POPC:DPPC mix). It's encoded using the latest CG Martini 3 release and I'm using Gromacs 2021.2 (May 2021). The system has four helical peptides. The simulation begins with the peptides arranged 5 nm apart and is then left to evolve. With the peptides being CG helical, they won't unravel and although each system of 4 peptides is unique, each helix is always 26 residues long.

I'm trying to design a metric to represent the oligomeric state. I could use a distance rule-based system but as the number of peptides increases the coding gets to be a real problem. As a smart alternative, I'm using the logisticContact() score between the protein mass and the bilayer/water environment. Of the possible oligomeric states in a 4 peptide system, I imagine the score order would be something similar to 4_monomers > 2_monomer_1_dimer > 2_dimers > trimer_monomer > 4_mer. Ideally, four monomers will have the greatest score as more protein is exposed to the lipid/water environment.

I've ran 50 simulations collecting the logisticContact() scores for each respective oligomeric state; this had to be done by watching 50 simulations and noting down the frame numbers associated with the above oligomeric states. It certainly took a while. Sadly, the results were a bit of a mess so I reduced the analysis down to the differences in two states over four randomly picked simulations.

I trialled the sigma (1, 2) and the radius (10, 14), and found that on the first 4 monomers vs 2 dimers I got very good median separation in logisticContact score when sigma was 2 and the radius 4. But this was not the case when I tried these values on a new set of 4 monomers vs 2 dimers. The median scores looked like this:

System A - 4 monomers - Median: 0.76, IQR: 0.21
System A - 2 dimers - Median: 0.44, IQR: 0.16
System B - 4 monomers - Median: 0.48, IQR: 0.19
System B - 2 dimers - Median: 0.45, IQR: 0.17

The System A shows a great separation in scores, however, scores from System B are proving to be a problem.

I was hoping I could get some insight from your experience building LOOS. The documentation states:

"Compute the number of contacts between the centroid of this AG and the centroid of another AG, using a smooth logistic function S = 1/(1 + dist/radius)**sigma"

1. Ideally, for a vdW potential short-range cut off between 1 nm and 1.4 nm, what in your experience should the radius be? I'm seeing a greater score as the radius increases. I presume that's because the point of the centroid will sit approximately in the middle of the helix.

2. Is the centroid calculated accounting for the mass of the molecule?

3. What is meant by "Compute the number of contacts"?

4. As I write this, I've just spotted something. In the documentation, it states: "between the centroid of this AG and the centroid of another AG". Could you help give me some insight as to what the centroid of a water box or a water-bilayer system would be? I'm finding this very hard to conceptualise.

[agrossfield]

Don't apologize for the lengthy question -- I think your bug reports/questions are exemplary, because they give me everything I need to figure out what's going on with little to nothing extraneous.

Answering your questions:

1. I'd pick the radius based on what the centroid-centroid distance is when the molecules you're interested in are in contact vs. not -- the vdw cutoff is a bit of a red herring, since that's atom-atom, and this is centroid-centroid. For lipid-lipid contact, we've typically used roughly 12 Ang, so I suspect it would need to be larger for a helix-lipid contact.
2. centroid() is not mass-weighted, while centerOfMass() is. In practice, the two will give virtually indistinguishable answers when you're averaging over something as large as a helix. Moreover, since this is MARTINI, all particles have the same mass, so they're literally the same. We got in the habit of using centroid() most of the time, because depending on what you're using to define your system you might not have masses (IIRC, we set default values of 1 so that codes don't break).
3. You've uncovered an error in the documentation, left over from an earlier version of the method that took a vector instead of a single group as the argument. I will fix that and update it as soon as I finish writing this note.The current version operates on the centroid of the calling group and one other.
4. You are totally correct that working with the centroid of the lipids as a whole would be nonsensical (and the water even worse, since the centroid would be roughly at the origin of the box, where there is no water!) The way to implement this would be to break your target up first, eg something like this:

protein_groups = protein.splitByMolecule()  # list of AtomicGroups, each containing 1 helix, assuming you already selected the proteins as a whole into protein
lipid_groups = lipids.splitByMolecule() # list of lipid molecules, assuming you already selected all lipids

total_contacts = 0.0
for helix in protein_groups:
    contacts = 0.0
    for lipid in lipid_groups:
        contacts += helix.logisticContact(helix, radius, sigma, box)
   total_contacts += contacts 

I wouldn't bother doing the water -- I strongly suspect water exposure won't change much with oligomerization of TM helices, at least not compared to lipid contacts, and obviously there are a lot of waters so it'll be slower. Also, you might consider using logisticContact2D(), because (assuming your membranes are small enough they're pretty flat) computing the distance in the xy plane might be closer to what you're thinking of, since otherwise there will be a roughly constant z component (the centroid of the helix will be at roughly z=0, while each lipid will be at either +8 ang or -8 ang (or something like that).