# How to see effect of point mutation in PyMOL?

I am trying to model how a specific point mutation would affect the protein structure. So far, I've figured out how to create the point mutation in PyMOL using the mutagenesis wizard. I'm not sure how to take the next step toward actually visualizing the consequences of the mutation though (i.e. would it disrupt a helix, etc). Is there a way to do this in PyMOL, or would I need other software? Thank you!

You could run a molecular dynamics simulation for example with Gromacs or Charm. The problem with pymol is that the structure is static and you don't know how much tension that point mutation creates. First of all you would run an energy minimization to remove bad non-physical overlaps between atoms, then add solvent molecules, and then simulate your system first with position restraints for equilibration and then without for 1ns to 100ns+.

Even then you cannot be sure that what comes out is the 'final' result as sometimes an induced fold change can take much longer. It depends on the system and the energy barriers included and the time scales that you are interested in. If the protein is large it may have many different physiologically relevant conformations and the point mutation might only change the equilibrium distribution of those states or totally prevent a transition.

In case you have more than one crystal structure (in different conformations) you may introduce these point mutations in all of them and simulate them for a while. But it is not that trivial to setup such a simulation you should contact an expert. There are many things that can go wrong. Other software you could use are Modeller and Rosetta, if I am not mistaken.

You could give SAAPpred a try. This tool uses machine learning (in the form of random forest) to predict the local effects of a single point mutation on the protein structure.

I did my PhD in the group and I used SAAPpred to find the most (predicted) damaging mutations of G6PD, to then ran extensive MD simulations to evaluate the final stability of those mutants.

Keep in mind that the tools doesn't return a structure, but just the possible/likely effects of that substitution.

A common metric used to assess if a mutation destabilises the structure is the difference in Gibbs free energy. PyMOL does not give a ∆∆G (difference in Gibbs free energy, kcal/mol), but a very crude steric clash nor does it repack neighbouring sidechains out of the way, or alter the backbone, so unfortunately is not the way forward. However, if one must there is a way (see SIMBAI below).

## Beyond destabilisation

Destabilisation predictors have a ~60% or lower success at predicting if something is pathogenic (ref). This is because the effect a mutation has on a protein is rather complicated as you can have destabilisation, but also aggregation, loss of regulation, mislocalisation, loss of binding to another protein, loss of binding to ligand, loss of activity (catalysis, inability to change conformation) etc etc. Then you can get curious subcases, such as loss of activity that result in the sequestration of a binding partner or homologues in a complex resulting in a worse phenotype than destabilisation or truncation.

## Modelling

The behaviour of a protein can be assess by the energy of the system using forcefields, which are multi-term equations that describes the atomic interactions (attraction, repulsion etc.). The result is the Gibbs free energy mentioned above, which is a potential, so negative is good. Kinetic energy is release when the unfolded protein "rolls down" down the potential energy funnel —ProTherm dataset (an empirical dataset of ∆∆Gs) and the tool SDM give this value inverted as the negative can be confusing.

An uphill change in the Gibbs free energy (∆∆G > 0) generally correlates with a destabilisation of a protein —aggregation and solubility being the caveat in the "generally". The destabilisation a single state (static conformation) can be assessed. gmx energy in Gromacs, FoldX or Rosetta can be used for this.

For altered ability to switch conformation steered MD simulations are required, but generally the ∆∆G at different states can be used.

For the non-destabilisation options, https://www.phosphosite.org/ and http://elm.eu.org/ are very useful.

I have been working on a server that gives a ∆∆G value and note if there are any nearby post-translational modifications, gnomAD or changes in potential linear motifs etc.: Venus (unpublished). However, there are many ∆∆G calculation servers out there and are compared in many papers (example), but it should be said that the dataset used in these comparisons is always ultra-biased despite how much is tried. Missense by the Phyre2 group (Michael Sternberg) give detailed information of the neighbourhood changes the mutation causes. In terms of extra info, Miscast gives Uniprot info.

## SIMBA-I via PyMOL

In terms of PyMOL, the SIMBA-I algorithm can be used by getting the RSA from it —I would not recommend it as it may underestimate the severity of unusual cases. But if one must, in PyMOL one can paste the following into the console with the wild type variant:

# defined by you so change accordingly
from_residue = '👾' # single letter
to_residue = '👾' # single letter
resi = 👾
chain = '👾'

## data
maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0,\
'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0,\
'W': 264.0, 'Y': 255.0, 'V': 165.0}
volumes = {'A': 88.6, 'R': 173.4, 'N': 114.1, 'D': 111.1, 'C': 108.5, 'Q': 143.8, 'E': 138.4, 'G': 60.1,\
'H': 153.2, 'I': 166.7, 'L': 166.7, 'K': 168.6, 'M': 162.9, 'F': 189.9, 'P': 112.7, 'S': 89.0,\
'T': 116.1, 'W': 227.8, 'Y': 193.6, 'V': 140.0}
hydrophobicities = {'A': 8.1, 'R': 10.5, 'N': 11.6, 'D': 13.0, 'C': 5.5, 'Q': 10.5, 'E': 12.3, 'G': 9.0,\
'H': 10.4, 'I': 5.2, 'L': 4.9, 'K': 11.3, 'M': 5.7, 'F': 5.2, 'P': 8.0, 'S': 9.2, 'T': 8.6,\
'W': 5.4, 'Y': 6.2, 'V': 5.9}
cmd.set('dot_solvent', 1)
cmd.set('dot_density', 3)
sasa = cmd.get_area(f'resi {resi} and chain {chain}')
rsa = sasa / maxASA[from_residue]
offset = -1.6
rsa_factor = +2.2
dV_factor = +0.59
dH_factor = -0.29
rsa_dV_factor = -0.54
rsa_dH_factor = +0.49
# calculate
dV = (volumes[to_residue] - volumes[from_residue]) / 100
dH = (hydrophobicities[to_residue] - hydrophobicities[from_residue])
simba = offset + rsa_factor * rsa + dV_factor * dV + dH_factor * dH + rsa_dV_factor * rsa * dV + rsa_dH_factor * rsa * dH
print -simba  # negative stabilised, positive destabilised


## Human reasoning

However, the major problem with all the tools mentioned is that human reasoning is often required to read the literature and make a solid hypothesis.

A resource I cannot stress enough, for human mutations, its utility is gnomAD, which contains phenotypically neutral mutations: for example, one really ought to not claim a homozygous missense is pathogenic because it is causing a loss of function by destabilising the protein, yet there are homozygous truncations roaming around in the human population.

The paper describing the crystal structure used and those of the homologues may point out key features which would missed by automated servers, such as secondary-structure elements that are speculated to perform certain roles. In fact, without doing any calculations seeing in PyMOL if a mutation causes a clash in the core of a protein, alters the active site, or changes a residue to a smaller one in a region known to alter conformation upon a given condition, is a lot more informative even if not quantitative...