Support to AMBER engine for ensemble equilibrations

[silviamenin]

Hi!
It would already be very useful to have the possibility to produce and customize the equilibration protocol, as well as to choose the MD engine for ensemble equilibration. In particular, AMBER could be a great option here, especially for membrane systems where recent improvements make it a strong candidate, also because it is already supported in Biosimspace.
Thanks for considering!

[Roy-Haolin-Du]

Thanks silvia, thanks for your valuable thinking.

As suggested by Julien,
not sure whether AMBER support the restraints selection
https://github.com/michellab/a3fe/blob/main/a3fe/run/leg.py#L682.
And there are a few things that would need updating to support interoperability between MD engines, such as the call above to _BSS.Trajectory.Trajectory.
I will start identifying what would need to be updated in the code, and whether code updates are needed in BioSimSpace itself

[Roy-Haolin-Du]

Introduction and Goal
This document outlines the proposed modifications for supporting for the AMBER engine within the A3FE workflow. The primary goal is to enable A3FE to utilize AMBER for running ensemble equilibration simulations, with the trajectories generated by AMBER subsequently used for Boresch restraint selection. The restraints generated will be formatted for use with the SOMD engine for the subsequent free energy calculation.

1. A3FE Codebase Modifications

• System Preparation (a3fe/run/system_prep.py): Modify run_ensemble_equilibration and run_process function: Implement conditional logic to instantiate and replace process = _BSS.Process.Gromacs(system, protocol, work_dir=work_dir) with _BSS.Process.Amber. And output File Handling: Adapt trajectory files (e.g., .nc) and final coordinate files (.rst7).

• Trajectory Loading (a3fe/run/leg.py): Modify traj = _BSS.Trajectory.Trajectory(trajectory=f"{outdir}/gromacs.xtc" ) to be trajectory=f"{outdir}/amber_trajectory.nc")# Example AMBER trajectory filename.

2. BioSimSpace Verifications and Dependencies.

• Trajectory Reading (_trajectory.py): Confirmation of AMBER .nc file reading by BioSimSpace/Trajectory/_trajectory.py. While it's likely supported. https://github.com/OpenBioSim/biosimspace/blob/devel/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py.

• Restraint Analysis (_restraint_search.py): Verification that the static RestraintSearch.analyse method in (which relies on MDAnalysis internally) works correctly produces restraints with AMBER trajectories. As noted, this is expected to work. https://github.com/OpenBioSim/biosimspace/blob/devel/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py.

• Restraint Formatting (_restraint.py): Just check restraint.toString(engine="SOMD") work well. Beacuse AMBER is only used for the initial ensemble equilibration and restraint generation, with FEP stages run in SOMD. https://github.com/OpenBioSim/biosimspace/blob/devel/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py.

Make sure Amber installation well and path, specifically the location of the AMBER executables, be correctly configured and accessible to A3FE env.

[fjclark]

Just commenting that this all looks great and will be a nice addition to a3fe. Thanks for the detailed plan @Roy-Haolin-Du! As you say, in theory this should be a minor change and all the functionality in BioSimSpace should work as is (although not sure I've ever tried it). The reason I went with GROMACS is that the setup stages seemed to be a bit more stable. Also, a minor terminology comment - constraints are a bit different to restraints in that constraints are absolute restrictions which must be satisfied ~ exactly (e.g. on bond lengths to H) whereas restraints are energetic biases.

[Roy-Haolin-Du]

In A3FE, Gromacs is used for minimization, NVT, NPT, and ensemble equilibration.

For membrane protein systems, special .mdp parameters are required.
(e.g. Semi-isotropic pressure coupling required (XY ≠ Z directions); Different compressibility for membrane; Gentle pressure control to avoid membrane destabilization and so on....(sorry not familiar with these)

Currently, these parameters are not automatically included. But allow to custom gmx .mdp file for membrane systems in configuration generation step.

The code works as follows:

In a3fe, the run_heat_and_preequil_input and run_ensemble_equilibration functions all use run_process in system_prep.py, which internally runs GROMACS.

a3fe/a3fe/run/system_prep.py

Lines 8 to 9 in f8e1948

	"heat_and_preequil_input",
	"run_ensemble_equilibration",

In BioSimSpace, _gromacs.py, there's an extra_options field that allows you to pass additional parameters, which is handled inside the internal method _generate_config by using generateGromacsConfig function. https://github.com/OpenBioSim/biosimspace/blob/8ba863f612baee764c8b5935df9c40dbc2daef1e/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py#L425.
The detailed gromacs.mdp file is generated via generateGromacsConfig function in _config.py, where the extra_options are incorporated for membrane protein-specific control parameters. https://github.com/OpenBioSim/biosimspace/blob/8ba863f612baee764c8b5935df9c40dbc2daef1e/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py#L421.

Temporary a3fe modification.
System Preparation (a3fe/run/system_prep.py):

class SystemPreparationConfig(_BaseModel):
    membrane_protein: bool = Fieldxxx  ### add mem pro config para
    def get_membrane_equilibration_config(self) -> dict: ### membrane pro config para
    def heat_and_preequil_input():
    def run_ensemble_equilibration():
    system = run_process(equilx, protocol, extra_options=cfg.get_membrane_equilibration_config()) ###NPT and equilibration with membrane support
    def run_process(extra_options: _Optional[dict] = None) ### add a new to let extra_options passed.
    process = _BSS.Process.Gromacs(system, protocol, work_dir=work_dir, extra_options=extra_options) #### BioSimSpace automatically merges extra_options into mdp via generateGromacsConfig

Then run:

import a3fe as a3
cfg = a3.SystemPreparationConfig(membrane_protein = True)
calc = a3.Calculation()
calc.setup(bound_leg_sysprep_config = cfg, free_leg_sysprep_config = cfg)
calc.run and analyse as usual...

Also, for the a3fe somd run part. SOMD support barostat_membrane.
https://github.com/OpenBioSim/sire/blob/67673e439ae17cbaca494c82b9d2db2d5fbcbdb3/wrapper/Tools/OpenMMMD.py#L271
Just change the template_config.cfg file

### Barostat ###
barostat = False
barostat_membrane = True

Future Enhancements：

1. Automatic membrane detection and temperature coupling groups (protein/membrane/solvent separation), also stable equilibration pressure protocols (Berendsen → Parrinello-Rahman).
2. BioSimSpace enhancements for native membrane config support for gromacs.

[Roy-Haolin-Du]

@silviamenin
Hi Silvia, I've made a simple modification in a3fe to support membrane protein, via gromacs for minimization, NVT, NPT, and ensemble equilibration.

The only change in the input required is that the use of barostat_membrane, rather than barostat, should be specified in template_config.cfg,

    ### Barostat ###
    barostat = False
    barostat_membrane = True

You can use:

    import a3fe as a3
    cfg = a3.SystemPreparationConfig(membrane_protein=True)
    calc = a3.Calculation()
    calc.setup(bound_leg_sysprep_config = cfg, free_leg_sysprep_config = cfg)

However, need provide the bound_solv and free_solv files externally, as automated parameterization and solvation are not yet supported.

Additionally, 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.

You can git clone my GitHub branch from:
https://github.com/Roy-Haolin-Du/a3fe/tree/feature-membrane-protein-eq. and give it a test.
Julien suggested merging it into michellab only after testing. I had created a PR #49. You’re also welcome to take a look.

Let me know if you run into any issues. And thanks a lot!

Best,