Protein Modeling with Wolfram Language

Robert B Nachbar, Wolfram Solutions

#WolframTechConf

ABSTRACT

Molecular modeling has been supported in Wolfram Language for many years, and initial support for biomolecules (proteins and nucleic acids) was recently added. The current capabilities include importing from PDB, sequence folding with ESMAtlas, visualization, and some property computation. In last year’s presentation we described the computation of protein-ligand binding energy, and the prerequisite biomolecule editing need to obtain an accurate result. That editing and computation were done on the Molecule representation of the protein. In this talk we will describe the direct editing of the BioMolecule representation of proteins to support the computation of binding energy, molecular dynamics, and accurate properties. Specific operations covered include supplying experimentally unobserved hydrogen atoms; adjusting the amino acid residue protonation state for acidic, neutral, and basic pH; and 3D point mutation of residues. Examples and discussion of the prototype code will be included.

Protein modeling

Aren’t proteins and nucleic acids just big molecules?

Yes and No

Biopolymer molecules are generally modeled differently than organic and inorganic small molecules

Proteins and nucleic acids (DNA, RNA) are polymeric

There is a limited number of allowed monomers, or residues

For proteins: Ala, Arg, Asn, Asp, Cys, Gln, Glu, Gly, His, Ile, Leu, Lys, Met, Phe, Pro, Ser, Thr, Trp, Tyr, Val

There is a standard numbering scheme for the residues

There are standard labels for the atoms in each residue

Protein structure data

There are three main sources of 3-dimensional protein structures:

Protein Data Bank (PDB)

Experimental: x-ray and neutron diffraction, NMR, cryo-EM

AlphaFold

Theoretical: deep learning NN sequence folding

ESMFold

Theoretical: deep learning NN sequence folding

All 3-dimensional structures are not created equal

Hydrogen atoms are often not observed by experimental methods, either completely or partially

Short stretches of residues are sometimes not observed experimentally due to disorder

Occasionally C, N, and/or O atoms are  not observed experimentally due to disorder

The NN methods omit all hydrogens

Thus the need for editing biomolecules to support the computation of:

binding energy

molecular dynamics

accurate properties

But wait, there’s more!

What about changing one residue into a different one, for example a phenylalanine into a tryptophane?

Row[{MoleculePlot[StandardProteinResidue["Phe"],ImageSize->Small],Style["  ⟶  ",28],MoleculePlot[StandardProteinResidue["Trp"],ImageSize->Small]}]

And, we want to do it without perturbing the rest of the protein structure

This is a very common task in protein design and analysis

In fact, it was a user request to change one amino acid residue into another and study the effect on the geometry of the protein with molecular dynamics (MD) that initiated this code development

The workflow for this task will provide the outline for this presentation

So, where do we begin?

Our approach has been to treat the protein as a set of interchangeable parts at the residue level

We can:

Extract residues as Molecule objects

Replace residues with Molecule objects (provided they are one of the 20 canonical amino acids)

We can also conform the shape of one residue to that of another

This allows the replacement with minimal interference from spatially nearby side chains

Even adding missing hydrogen atoms and changing the ionization state is handled residue by residue

All of the above functionality is made possible using a dataset of standard protein residues with 3-dimensional geometries

The result is a new ProteinEditing paclet

It is still in initial development, and available is here

Standard protein residues

The data for the 20 canonical (proteinogenic) amino acids came from the Protein Data Bank (PDB)

wwPDB Protonation Variants Companion Dictionary

These are the building blocks for doing point mutation of the 3-dimensional geometry of proteins

They have also proven to be useful for

Assessing the completeness of an experimentally determined structure

Adding and removing hydrogen atoms

There are the various backbone linkage forms:

N-terminal

C-terminal

Linkage

Free amino acid

Zwitter ion

StandardProteinResidue["Dataset"][Select[MatchQ[#Tuple,{"ALA",___}]&],{"Molecule","Tuple"}]//Normal//MapApply[MoleculePlot[#1,PlotLabel->#2,ImageSize->75,BaseStyle->6]&]//Partition[#,UpTo[5]]&//Grid[#,Frame->All,FrameStyle->LightGray]&

There are various protonated forms depending on the specific amino acid

StandardProteinResidue["Dataset"][Select[MatchQ[#Tuple,{"HIS",___}]&],{"Molecule","Tuple"}]//Normal//MapApply[MoleculePlot[#1,PlotLabel->#2,ImageSize->75,BaseStyle->6]&]//Partition[#,UpTo[5]]&//Grid[#,Frame->All,FrameStyle->LightGray]&

Note that the N-terminal (LSN3), linking (LL), and C-terminal (LEO2, LEO2H) forms have a negative formal charge on the N and/or C atoms which represent the connection points to adjoining residues

StandardProteinResidue

To get started, let’s first load the paclet

Needs["RobertNachbar`ProteinEditing`"]

The function StandardProteinResidue gives the user access to these structures

?StandardProteinResidue

Here is an N-terminal histidine

resmol=StandardProteinResidue[{"HIS","LSN3"}]

Row[{MoleculePlot[resmol],MoleculePlot3D[resmol]}]

Here is a linking histidine with a deprotonated imidazole

resmol=StandardProteinResidue[{"HIS","LL","DHD1"}]

Row[{MoleculePlot[resmol],MoleculePlot3D[resmol]}]

All the atoms have standard labels:

Row@moleculePlotWithLabels[{StandardProteinResidue[{"HIS","LL"}],StandardProteinResidue[{"ASP","LL","DHD2"}]}]

They are used internally by the code, and by the user in ProteinResidueConform

We plan to use them to allow the user to, for example, compute the distance from HE2 of histidine residue 3 to OD2 of aspartic acid residue 34

The details need to be worked out

Extracting a residue as a Molecule expression

To replace a residue with minimal disturbance to the rest of the structure, one would like the new residue to have roughly the same shape as the one being replaced

This seemed most easily accomplished by working with the residues in isolation, thus the need to extract a residue

?ProteinResidueExtract

First, let’s get a protein, the avian pancreatic polypeptide:

biomol=BioMolecule[(ExternalIdentifierURL[ExternalIdentifier["PDBStructureID","1PPT"], "PDB structure", "1PPT", https://rdf.wwpdb.org/pdb/1PPT])]

BioMoleculePlot3D[biomol,ColorRules->{{"A",16}->Black,Splice@blueToRed[biomol,"A"]}]

The residues are colored blue at the N-terminal end of the chain through orange, yellow, and green to red at the C-terminal end

Residue 16 is near the beginning of the α-helix, and has been colored black

We can get the type of residue from the "ResidueIDs" and "Residues" properties of the biomolecule:

MapThread[Thread[#1->#2]&,{biomol["ResidueIDs"],biomol["Residues"]}]

So, residue 16 is an aspartic acid (Asp, D)

resAsp=ProteinResidueExtract[biomol,{"A",16}]

Row[{MoleculePlot[resAsp],MoleculePlot3D[resAsp]}]

Note that there are no hydrogens in the 3-dimensional structure

We can see how the extracted residue fits in the full protein

The atom coordinates of the extracted residue are exactly the same as those of the in situ residue

Show[{BioMoleculePlot3D[biomol,ColorRules->blueToRed[biomol,"A"]],MoleculePlot3D[resAsp,PlotTheme->"Spacefilling"]}]

Conforming one residue to another

Let’s suppose we want to replace the Asp with a phenylalanine (Phe, F)

resPhe=StandardProteinResidue[{"Phe","LL"}]

Row[{MoleculePlot[resPhe],MoleculePlot3D[resPhe]}]

The process of conforming a residue has two steps:

Align the backbone atoms (N, CA, C, and O) of the Phe to those of the Asp

Rotate the flexible bonds of the Phe to maximize the overlap of the backbone and side chain atoms with those of the Asp

?ProteinResidueConform

resPheConformed=ProteinResidueConform[resAsp,resPhe]

MapThread[MoleculePlot3D[#1,IncludeHydrogens->"ExplicitOnly",ColorRules->{"C"->#2}]&,{{resAsp,resPheConformed},{StandardGreen,StandardYellow}}]//Show//Legended[#,SwatchLegend[{StandardGreen,StandardYellow},{"Asp","Phe"}]]&

We can make the phenyl ring of the Phe align with the carboxylate oxygen atoms of the Asp by manually specifying the atom mapping

resPheConformed=ProteinResidueConform[resAsp,resPhe,{"CG"->"CG","OD1"->"CD1","OD2"->"CD2"}]

MapThread[MoleculePlot3D[#1,IncludeHydrogens->"ExplicitOnly",ColorRules->{"C"->#2}]&,{{resAsp,resPheConformed},{StandardGreen,StandardYellow}}]//Show//Legended[#,SwatchLegend[{StandardGreen,StandardYellow},{"Asp","Phe"}]]&

Replacing a residue in a BioMolecule

Residue replacement is done with ProteinResidueReplace

?ProteinResidueReplace

D16F=ProteinResidueReplace[biomol,{"A",16}->resPheConformed]

Show the new biomolecule D16F with it’s residue 16 rendered as space-filling:

Show[{BioMoleculePlot3D[D16F,ColorRules->blueToRed[D16F,"A"]],MoleculePlot3D[ProteinResidueExtract[D16F,{"A",16}],PlotTheme->"Spacefilling"]}]

Adding unmodeled hydrogens

Recall that the experimentally determined structure that we retrieved from PDB did not have any hydrogens

The function ProteinSetHydrogens can add them

Moreover, it can add or remove hydrogen atoms on ionizable moieties for a given pH

?ProteinSetHydrogens

So, here is D16F

BioMoleculePlot3D[D16F,PlotTheme->"Spacefilling"]

Note the lone residue with hydrogen atoms on the right, it’s the Phe that we used to replace the Asp

Add the missing hydrogens using the default pH 7.4:

D16F7=ProteinSetHydrogens[D16F];
BioMoleculePlot3D[%,PlotTheme->"Spacefilling"]

All the Asp and Glu residues are ionized at pH 7.4

At pH 1 they will be protonated:

D16F1=ProteinSetHydrogens[D16F7,1];
BioMoleculePlot3D[%,PlotTheme->"Spacefilling"]

At pH 13 both the acidic and the basic moieties will be deprotonated:

D16F13=ProteinSetHydrogens[D16F7,13];
BioMoleculePlot3D[%,PlotTheme->"Spacefilling"]

Here is a census of the atoms in the various versions of the protein:

biomol["AtomicSymbols"]//Values//Flatten//Counts//KeySort

D16F7["AtomicSymbols"]//Values//Flatten//Counts//KeySort

D16F1["AtomicSymbols"]//Values//Flatten//Counts//KeySort

D16F13["AtomicSymbols"]//Values//Flatten//Counts//KeySort

Odds and ends

We mentioned at the beginning of the presentation that the 3-dimensional protein structures from the usual sources are often less than adequate for modeling

We’ve already seen how the "AtomicSymbols" property can be used to discover if hydrogens atoms are included

Here is the structure of the protein crambin from AlphaFold:

crambin=Import[FileNameJoin[{NotebookDirectory[],"AF-P01542-F1-model_v6.pdb"}]]

crambin["AtomicSymbols"] // Values // Flatten // Counts // KeySort

BioMoleculePlot3D[#,PlotTheme->"Spacefilling"]&/@{crambin,ProteinSetHydrogens[crambin]}//Row

We added three new properties to BioMolecule to make it easier to identify other types of problems

"NoncanonicalResidues"

"UnmodeledResidues"

"PartiallyModeledResidues"

The protein 5V63 contains the noncanonical residues ornithine, sarcosine, and iodophenylalanine:

BioMolecule[(ExternalIdentifierURL[ExternalIdentifier["PDBStructureID","5V63"], "PDB structure", "5V63", https://rdf.wwpdb.org/pdb/5V63])]["NoncanonicalResidues"]

We plan to expand the dataset of standard residues

Several N- and C-terminal residues of the protein 4LHF were not observed:

BioMolecule["4LHF"]["UnmodeledResidues"]

We are experimenting folded sequences

Several side chain atoms were not observed in the first two residues of 1FNA:

BioMolecule["1FNA"]["PartiallyModeledResidues"]

These residues can be replaced with the corresponding standard residues

Thank You

Questions

Initialization

Quit[]

PacletInstall[ResourceObject["https://wolfr.am/1xZtE0ewy"]]

Get["RobertNachbar`ProteinEditing`"]

SetOptions[MoleculePlot3D,IncludeHydrogens->"ExplicitOnly"];

$standardResidues={"Ala","Arg","Asn","Asp","Cys","Gln","Glu","Gly","His","Ile","Leu","Lys","Met","Phe","Pro","Ser","Thr","Trp","Tyr","Val"};

blueToRed[bm_BioMolecule,chain_String]:=With[{n=bm["ResidueCounts"][chain]},Thread[Thread[{chain,Range@n}]->Table[Hue[2/3((-i+n)/(-1+n))],{i,n}]]];

ClearAll[moleculePlotWithLabels]
Options[moleculePlotWithLabels]=Options[MoleculePlot];
moleculePlotWithLabels[mols:{__?MoleculeQ},opts___]:=moleculePlotWithLabels[#,opts]&/@mols
moleculePlotWithLabels[mol_Molecule,opts:OptionsPattern[]]:=MoleculePlot[mol,IncludeHydrogens->"ExplicitOnly",AtomLabels->Thread[Range@AtomCount[mol,IncludeHydrogens->"ExplicitOnly"]->mol["MetaInformation"]["AtomLabels"]],opts,AtomLabelStyle->Directive[StandardBlue,Italic,10]]