Question about protein relative free energy

[obasekigit]

Good day all, I would like to ask if biosimspace can handle charge changes when carrying out protein relative free energy calculations for example valine to lysine

[akalpokas]

Hi @obasekigit, in terms of the creation of the perturbable protein modifications involving charge change, yes BioSimSpace can create these kinds of transformations. However, once you save your system and when it comes to the actual running of the dynamics you then have an issue with a charge-changing transformation. To neutralize the charge change during the transformation, we in general tend to use the co-alchemical water approach, but to the best of my knowledge we don't have a function to call for setting it up easily. That doesn't mean it's difficult to set it up manually, the code that I have used in the past looks something like this:

First we create a solvated protein system based on our WT protein (which you would load anyway during the protein mutation set up):

box_size = 3 # nanometers
box = 3 * [box_size * BSS.Units.Length.nanometer]
solvated_co_alc_system = BSS.Solvent.tip3p(molecule=wt_protein, box=box, ion_conc=0.15)

Since you are asking about valine --> lysine transformation (0 --> +1), you would need to transform a sodium ion to a water molecule to balance out the charge change (+1 --> 0). What I have done in the past is to pick a water and an ion away from the protein (to prevent them coming into a contact with each other as much as possible), create a transformation and then add it back to the original perturbed protein system:

# extract a water molecule and ion from the solvated system
furthest_water = solvated_co_alc_system.search("furthest 1 water from resname GLY") # doesn't have to be resname GLY, can use entire protein for searching
furthest_ions = solvated_co_alc_system.search("furthest 1 (atomname NA) from resname GLY")
furthest_water_mol = furthest_water.molecules()[0]
furthest_ion = furthest_ions.molecules()[0]
furthest_water_mol.translate([-1.5, -1.5, -1.5])
furthest_ion.translate([-1.5, -1.5, -1.5])

# scoring function rmsd_align fails
co_alchemical_water_map = BSS.Align.matchAtoms(furthest_ion, furthest_water_mol, scoring_function="rmsd")

co_alchemical_water = BSS.Align.merge(furthest_ion, furthest_water_mol, co_alchemical_water_map)
co_alchemical_water.toSystem()

The furthest_water_mol.translate([-1.5, -1.5, -1.5]) bit is important. It translates the molecule by -1.5 A in xyz direction (might need to be tweaked in your case). Without it, the co-alchemical water would be placed on top of your previous ion/water molecule that already exists in the system, which can lead to the simulation exploding during the minimization.

You can then simply add the co-alchemical system to your previous protein perturbation system and export everything:

co_alchemical_system = co_alchemical_water + merged_solvated_system

Having said all this, we do have support for automatically the adding co-alchemical water during dynamics run time with our new OpenMM based FEP engine SOMD2, which is BioSimSpace perturbation compatible, although the engine is still in the experimental/validation stage. Hopefully this answers your question.

[lohedges]

The addition of a co-alchemical is also (semi) automatic for SOMD1. I'll look up the configuration option when I'm back at the computer. I think you just need to specify the charge difference between the end states in the config file.

[obasekigit]

Thank you so much for the detailed explanation. I look forward to seeing the co-alchemical configuration file for SOMD1

[lohedges]

Here's the config option:

somd-freenrg -H | grep -A2 'charge difference' 
charge difference = 0
The difference in net charge between the two states

So, you should be able to enable it in somd1 by passing, e.g. extra_options = {"charge difference": 1}, when creating the FreeEnergy.Relative or Process.Somd object.

Note that the charge difference specifies the the charge of the lambda = 0 state minus the charge of the lambda = 1 state, so if the value is -1 then a water in the lambda = 0 state will be perturbed to a chlorine ion at lambda = 1. Note also that the water used as the the alchemical ion is just chosen at random, i.e. there is no constraint on where it is relative to the protein or perturbable region. There are also no restraints on the position of the ion, so it can drift around over time. As such, you might want to check the log to see what water was chosen as the ion, and re-run if it is inappropriately close to other charge groups. In somd2 we always choose the water that is furthest from the ligand's center of geometry and also have support for additional restraints if required.