Support membrane protein in a3fe via Gromacs-based System Preparation

[Roy-Haolin-Du]

Description

Ensure ABFE calculations with membrane proteins are supported via gromacs for minimization, NVT, NPT, and ensemble equilibration.

• The input required is that the use of barostat_membrane, rather than barostat, should be specified in template_config.cfg.
• Modified the SystemPreparationConfig object via cfg = a3.SystemPreparationConfig(membrane_protein=True).
• Users must provide the bound_solv and free_solv files externally, as automated parameterization and solvation are not yet supported.
• The configuration files (gromacs.mdp) for these preparation steps specified in get_membrane_equilibration_config function of system_prep.py file in which is general and not highly specialized; Thus, it is recommended that users perform these pre-equilibration steps externally. Future updates will address these limitations.

Close #48.

Todos

• More robust code logic.

Automatic detection of membrane protein systems.

Determine right barostat usage

• Write user guide documentations.
• Design the test logic carefully.

Status

• Ready to go
• Waiting for real-case testing.

[fjclark]

Thanks for this Roy - looks good! Let me know if/ when you'd like a review and I'm happy to take a look!

[Roy-Haolin-Du]

@fjclark Thanks, Finlay, I definitely appreciate your advice!
(I was worried you might be too busy to bother, 😄)

The current idea is to first modify the system preparation in A3FE to support membrane-specific configurations, and later on, BioSimSpace will natively support membrane configuration for gromacs.

The change I made is pretty low-level, I simply took advantage of _BSS.Process.Gromacs(extra_options) to pass in membrane-related parameters. However, the current gmx.mdp is somewhat crude and need further refinement. I’ll discuss this part in more detail with Julien.

So could please review it? I’d be sincerely grateful for any suggestions or comments you might have.

In addition, Julien suggested having Silvia test it before merging. I’ll also listen to her feedback.

Thanks.

[fjclark]

Looks good, thanks Roy! Just a couple of minor/ pedantic comments.

[fjclark]

Thanks, docs look great! Only thing (which you were probably planning anyway) is to update the changelog.

[Roy-Haolin-Du]

Mention for myself.

For temperature coupling groups (protein/membrane/solvent separation) in .mdp file, tc-groups = SOLU MEMB SOLV

We can create a custom wrapper for gmx make_ndx in a3fe to generate membrane protein temperature coupling groups.

SOLU: Protein + Ligand
MEMB: Membrane lipids (PA + PC + OL)
SOLV: Solvent (ions + water)

However, the limitation discovered
BioSimSpace's _generate_binary_run_file method lacks -n index.ndx support
The gmx grompp command generated by BioSimSpace is:
gmx grompp -f gromacs.mdp -po gromacs.out.mdp -c gromacs.gro -p gromacs.top -r gromacs_ref.gro -o gromacs.tpr
from https://github.com/OpenBioSim/biosimspace/blob/baad60690aff9fc09da3a3a7edfdee11c74639b8/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py#L628-L647.

Practical Workaround: Use gmx predefined groups
tc-grps = "Protein non-Protein" # Instead of "SOLU MEMB SOLV"
,which is less granular but functional within BioSimSpace's current limitations.

[Roy-Haolin-Du]

Add membrane protein temperature coupling with Protein non-Protein groups.

• Support separation of Protein and Non-Protein for temperature coupling configuration in NVT, NPT and ensem equil.
• Support bound/free leg configuration differences in equilibration steps.

As BSS seems not support custom index files for these groups.(-n index.ndx) .OpenBioSim/biosimspace#405.
We use gmx predefined groups
tc-grps = "Protein non-Protein" # Instead of "SOLU MEMB SOLV".

Hi Finlay @fjclark
Would you mind having a quick look at these two small changes，please?
Thanks a lot!

[Roy-Haolin-Du]

Hi Finlay,
I believe you must be very busy, but that's totally fine.
I have already completed a full test workflow run and analysis, got ΔG.
I will go ahead and merge it into the michellab feature-membrane-protein-eq branch, as we have a meeting with Silvia’s membrane team this Friday.
Of course, you can review it whenever you have time afterwards.
Thank you!

[silviamenin]

@Roy-Haolin-Du I don't know if you mentioned it already, but when simulating the FREE leg in the context of membrane protein systems, the barostat should be the normal isotropic one in the somd.cfg.

In my personal version I put a switch in the leg.py, after line 832:

if self.leg_type == _LegType.FREE:
    options_to_write['barostat'] = True
    options_to_write['barostat_membrane'] = False

In this way, the bound leg with the membrane is simualted with the semi-isotropic barostat, while the free leg it is always simulated using the isotropic one, since the system is always a simple water box with the ligand

[Roy-Haolin-Du]

Hi @silviamenin Thanks a lot for pointing this out, I didn’t know!
I’ll make default barostat and barostat_membrane to False in the template, and rewrite barostat_membrane to True only if a membrane is detected from the system.
Great suggestion and much appreciated!