Free Energy

Creating a Alchemical Free Energy Simulator Project

Let’s walk through how to create a Alchemical Free Energy project. Before running the Alchemical Free Energy Simulator, you must first build the following components:

  • PDB Reader

  • Solution Builder or Membrane Builder

  • Ligand Reader & Modeler

Both the Solution/Membrane Builder and the Ligand Reader and Modeler must use the same PDB Reader, which must include the reference ligand. All builders must employ the same protein force field. Be especially careful with the ligand force field(s) applied in the Ligand Reader and Modeler—you can use multiple ligand force fields if needed.

create_project options:

  • title (str): project title

  • referenc_project_id (str): Solution Builder or Membrane Builder project id

  • ff (ForceField): project force field ‘amberff’ or ‘charmff’ (default: ForceField.AMBER)
    • PROGRAM: Available ForceField

      • AMBER: AMBERFF

      • GROMACS: AMBERFF, CHARMMFF

      • OpenMM: AMBERFF, CHARMMFF

      • BLaDE: AMBERFF, CHARMMFF

  • calculate_program (Program): (default: Program.AMBER)

  • calculate_application (Application): (default: Application.RELATIVE_BINDING)

  • amber_options (dict[str, Enum]): The default selections for amberff are as follows:

    {
    "protein": Protein.FF19SB,
    "water": Water.OPC
}

NOTE:

Choose Program and Application according to your free energy calculation goals:

  • Program: AMBER, GROMACS, OpenMM, BLaDE

  • Application:
    • Relative Binding F.E. (RBFE)

      • AMBER-TI, GROMACS Noneq-TI and OpenMM FEP: compute relative free energy between ligand A–B as defined by each morph

    • Absolute Binding F.E. (ABFE)

      • AMBER-TI and GROMACS Noneq-TI: perform absolute comparison of protein–ligand interaction energies

    • Mutation F.E.

      • AMBER-TI and GROMACS Noneq-TI: compute relative change in binding affinity resulting from a specific amino acid substitution

Here is an example workflow for AMBER-TI Relative Binding:

from molcube.free_energy.enum import *

# alt: solution_project_id = sbuilder.projectId
solution_project_id = '1910f751-a91d-450b-8634-ba3c6349dca5'

# Create a Alchemical Free Energy Simulator project
afes = molcube.create_afes_project()
afes.create_project(
    title='my_afes_project',
    reference_project_id=solution_project_id,
)

After creating the project, you can check the reference ligands using the check_reference_ligands method.

>>> available_ligands = afes.check_reference_ligands()
Chain Index: HETE_C, Residue Name: 53U

You can read a Ligand Reader & Modeler project using the read_ligand_modeler_project method.

>>> afes.read_ligand_modeler_project(lrmbuilder.projectId)
Available ligands for openffFast:
  - 1a
  - 1b

AMBER-TI Relative Binding

  • BASE WORKFLOW (after create_project())

    1. afes.generate_morph_set() - Based on the imported LRM project, recommend and generate both closed and radial morph sets:

    • Closed morph: Closed minimal perturbation path based on the clustering. (Fingerprint-based)

    • Radial morph: Radial shape of perturbation path from the ligand 1 (no clustering). (Fingerprint-based)

    • Minimum Spanning Tree morph (MST): Minimum spanning tree of the perturbation path. (MCS-based)

    • LOMAP morph: LOMAP-based perturbation path generation. (MCS-based)

    • Custom morph: Manually define perturbation path in the next step.

    • opitons:

      • reference_ligand (dict[str]): e.g. {“chainIndex”: ‘YOUR_LIGAND_CHAIN_INDEX’, “residueName”: “YOUR_LIGAND_RESIDUE_NAME”}

      • ligand_modeler_project_id (str): Ligand Reader & Modeler Project ID. If you have already executed the read_ligand_modeler_project method, you do not need to provide this argument.

      • ligand_ff (list[LigandForceField]): Select the ligand force fields generated in the LRM project that you want to use for Free Energy Calculation.

      • selected_ligands (dict[LigandForceField, list[str]]): Specify which ligands generated in the LRM project will be used for Free Energy Calculation, grouped by ligand force field. If you use an empty dictionary ({}), all selectable ligands will be included. (default: {})

        e.g. {LigandForceField.GAFF_FAST: [“LIG1”, “LIG2”], LigandForceField.GAFF: [“LIG3”]}

      • generate_morphset_type (GenerateMorphsetType): morphset type (default: GenerateMorphsetType.CLOSED)

      • threshold (int): threshold for clustering (default: 1)

    2. afes.processing_ligand() - Processes ligand information and prepares data for free energy calculations.

    3. afes.generate_system() - Generates a molecular system for free energy calculations based on user-specified parameters and preprocessed results. - options:

    • margin (float): Water box margin for the system. (default: 22.5)

    • ion_type (str): Type of ions to use. (default: “kcl”)

    • ion_conc (float): Ion concentration. (default: 0.15)

    • use_hmr (bool): Whether to use hydrogen mass repartitioning. (default: False)

    • temperature (float): Temperature for the system in Kelvin. (default: 310.0)

  • Optional Methods

    a. After ``generate_morph_set()`` - afes.delete_morph_set()

    • Removes a morph set from the generated morph sets.

    • options:

      • lig_ff (LigandForceField): ligand force field to delete among the generated morph sets

      • morph_set (dict[str]): {‘start’: INDEX or Name, ‘end’: INDEX or Name}

        e.g. {‘start’: 1, ‘end’: 2}

    • afes.add_morph_set():

      • Adds a morph set to the generated morph sets.

      • options:

        • lig_ff (LigandForceField): ligand force field to add among the generated morph sets

        • morph_set (dict[str]): {‘start’: INDEX or Name, ‘end’: INDEX or Name}

          e.g. {‘start’: 1, ‘end’: 2}

    • afes.get_similarity_score():

      • Returns the similarity score between ligands for the specified ligand force field.

      • options:

        • lig_ff (LigandForceField): ligand force field to check similarity score for.

    • afes.check_ligand_2d_structure():

      • Displays the 2D structures of two ligands for visual comparison.

      • options:

        • file_name_1 (str): The file name of the ligand you want to check e.g. ‘1a.sdf’

        • file_name_2 (str): The file name of the ligand you want to check e.g. ‘1b.sdf’

    b. After ``processing_ligand`` - afes.get_core_indices()

    • Returns the core atom indices for each ligand in the morph set or ligand group.

    • options:

      • lig_ff (LigandForceField): ligand force field to check the core indices for.

      • set_num (int): set number to check core indices (defualt: 1)

    • afes.delete_core_indices()

      • Removes core atom indices from a ligand in the morph set or ligand group.

      • options:

        • lig_ff (LigandForceField): ligand force field from which to remove core indices.

        • set_num (int): set number from which to remove core indices. (default: 1)

        • indices (list[int]): list of core indices to remove

        • ligand_name (str): Name of the ligand from which to remove core indices

    • afes.add_core_indices()

      • Adds core atom indices to a ligand in the morph set or ligand group.

      • options:

        • lig_ff (LigandForceField): ligand force field to which to add core indices.

        • set_num (int): set number to which to add core indices. (default: 1)

        • indices (list[int]): list of core indices to add

        • ligand_name (str): Name of the ligand from which to add core indices

Now, let’s return to the AMBER-TI Relative Binding example.

The following method, based on the imported LRM project, recommends and generates both closed and radial morph sets.

Note:

  • Closed morph: Closed minimal perturbation path based on the clustering. (Fingerprint-based)

  • Radial morph: Radial shape of perturbation path from the ligand 1 (no clustering). (Fingerprint-based)

  • Minimum Spanning Tree morph (MST): Minimum spanning tree of the perturbation path. (MCS-based)

  • LOMAP morph: LOMAP-based perturbation path generation. (MCS-based)

  • Custom morph: Manually define perturbation path in the next step.

>>> selected_ligands = {
...     # or ['1a', '1b', '1c']
...     LigandForceField.OPENFF_FAST: ['1a.sdf', '1b.sdf'],
... }

>>> afes.generate_morph_set(
...     ligand_ff=[LigandForceField.OPENFF_FAST],
...     selected_ligands=selected_ligands,
...     generate_morphset_type=GenerateMorphsetType.CLOSED,
... )

Force Field: openff-fast
 Morph 1: 1a(1) -> 1b(2)

Morph set generated with threshold 1 for 1 force fields.
To update the morph set, use 'delete_morph_set()' or 'add_morph_set()' method.

The following method processes ligand information and prepares data for free energy calculations.

>>> afes.processing_ligand()
1a.sdf (51)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

1b.sdf (51)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

This produces an image highlighting the maximum common subgraph (MCS). If no alterations are made, these atoms will be used as the basis for alignment. If you want to add or remove atoms from the alignment group, you can use the add_core_indices() and delete_core_indices() methods.

>>> afes.delete_core_indices(
...     lig_ff=LigandForceField.OPENFF_FAST,
...     set_num=1,
...     ligand_name='1a', indices=[14])
First name: 1a, End name: 1b
1a.sdf (50)
1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

1b.sdf (50)
1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

>>> afes.add_core_indices(
...     lig_ff=LigandForceField.OPENFF_FAST,
...     set_num=1,
...     ligand_name='1a', indices=[52])
Core index 52 does not exist in the original indices.
1a.sdf (50)
1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

1b.sdf (50)
1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

After each change, an image is generated highlighting the new alignment group atoms.

Finally, we will generate a molecular system for free energy calculations based on user-specified parameters and preprocessed results. The following arguments are available for the method:

  • margin (float): Water box margin for the system. (default: 22.5)

  • ion_type (str): Type of ions to use. (default: “kcl”)

  • ion_conc (float): Ion concentration. (default: 0.15)

  • use_hmr (bool): Whether to use hydrogen mass repartitioning. (default: False)

  • temperature (float): Temperature for the system in Kelvin. (default: 310.0)

The available options for ions are as follows:

  • ions: [‘kcl’, ‘nacl’, ‘cacl2’]

afes.generate_system()
afes.download_project(filename="afes.tar.gz")