Adding new modifications¶
In Adding new residues and links we added a whole new Block in order to describe a phosphoserine residue. This had the (dis)advantage that we had to redefine all the default serine interactions and parameters as well as the inter-residue links. There must be a better way!
Fortunately there is. Rather than describing a whole new SEP Block and all that entails we can instead describe just the way the phosphorylation modified the normal SER residue. This is exactly what modifications are for. As before, we need to add the new modification in 3 places: input force field, output force field, and mapping between the two. Note that the parameters presented here are for demonstration purposes only and not fit for actual science or simulations!
All input files for this tutorial can be found on github.
The input force field¶
During repair the regular SER atoms will be repaired, the missing hydrogen (HG) will be added (!), and the phosphate atoms will be annotated as being “nonstandard”. During Identify modifications the Processor will try to identify these tagged atoms by finding a minimal set of modifications that describe all relevant atoms. For a modification to apply here it must be subgraph isomorphic to the input structure.
If we run martinize2 -f ala-sep-ala.pdb -o topol.top -x AJA.pdb
we get, as
expected, the warning that not all modifications could be identified:
WARNING - unknown-input - Could not identify the modifications for residues ['SER3'], involving atoms ['21-O1', '22-O2', '23-O3', '24-P']
So let’s define the modification in force_fields/charmm/modification.ff
:
; THESE PARAMETERS ARE FOR DEMONSTRATION PURPOSES ONLY. DO NOT USE.
[ modification ]
SER-phos
[ atoms ]
O1 {"element": "O", "PTM_atom": true}
O2 {"element": "O", "PTM_atom": true}
O3 {"element": "O", "PTM_atom": true}
P {"element": "P", "PTM_atom": true}
OG {"element": "O"}
HG {"element": "H", "replace": {"atomname": null}}
[ edges ]
OG P
P O1
P O2
P O3
As before, the input force field does not define any MD parameters or
interactions. This modification contains nodes and edges. The edges are not
very interesting, and just define the connections between nodes. Nodes on the
other hand define 2 things: 1) the atom name as is should be (first column),
and 2) any constraints the node must satisfy during the subgraph isomorphism
as a JSON formatted mapping. The constraints should define at least 2
properties: the element
, and PTM_atom
. The element property is
self explanatory, but the PTM_atom
needs more explanation.
Modifications contain 2 types of nodes:
Nodes that are already described by the parent block (
PTM_atom
is false, this is the default). We call these nodes “anchors”.Nodes that are not yet described by the parent block (
PTM_atom
is true).
In addition, it’s worth noting that repair
reconstructed the HG atom (see -write-repair
) since it’s not in the input
PDB. We use the “replace” property to describe all node attributes that need
to change because of this modification. In this case we indicate that this atom
should be removed again, by setting its atomname to “null”. You can use the
“replace” property to change any node property, including e.g. atom type
and charge.
The output force field¶
We have to add a similar modification for the output force field in
force_fields/martini3001/modification.ff
:
[ modification ]
; THESE PARAMETERS ARE FOR DEMONSTRATION PURPOSES ONLY. DO NOT USE.
SER-PO4
[ atoms ]
BB {"PTM_atom": false}
SC1 {"PTM_atom": false, "resname": "SER", "replace": {"atype": "Q5n", "charge": -1}}
[ bonds ]
BB SC1 1 0.33 5000
Nothing new here compared to the modification for the input force field. Note that here we do define the simulation parameters, and we define a bond.
The mapping¶
Finally, we need to add the mapping describing how to get from charmm to
martini3001 in mappings/SEP.mapping
:
; THESE PARAMETERS ARE FOR DEMONSTRATION PURPOSES ONLY. DO NOT USE.
[ modification ]
[ from ]
charmm
[ to ]
martini3001
[ from blocks ]
SER-phos
[ to blocks ]
SER-PO4
[ from nodes ]
N
HN
CA
HA
C
O
CB
HB1
HB2
[ from edges ]
N HN
N CA
CA HA
CA C
C O
CA CB
CB HB1
CB HB2
CB OG
[ mapping ]
N BB
HN BB
CA BB
HA BB 0
C BB
O BB
CB BB
CB SC1
HB1 SC1 0
HB2 SC1 0
OG SC1
P SC1
O1 SC1
O2 SC1
O3 SC1
Firstly, notice that this is a different file format than the backwards format
we used before. In this case we have to define between which force fields we’re
going to define a mapping (charmm
and martini3001
), and between which
modifications (or blocks) (SER-phos
and SER-PO4
). This mapping has to
define how to map the phosphate moiety (at least). This moiety will be mapped
to the SC1 bead, so we will need to describe the complete mapping for that bead.
In addition this mapping affects the mapping of the BB bead (since CB will now
also contribute in part to it).
The charmm modification already define some nodes (see above), but not all the nodes required to describe the complete mapping for the BB and SC1 nodes, so these need to be described under from nodes and from edges. Finally, the actual mapping section should be self explanatory.
Now if we run martinize2 -f ala-sep-ala.pdb -x AJA.pdb -o topol.top -ff-dir force_fields/ -map-dir mappings/
we see INFO - general - Applying modification mapping ('SER-phos',)
Now we need to check the produced itp file:
; THESE PARAMETERS ARE FOR DEMONSTRATION PURPOSES ONLY. DO NOT USE.
[ atoms ]
1 Q5 1 ALA BB 1 1
2 TC3 1 ALA SC1 2 0.0
3 P2 2 SER BB 3 0.0
4 Q5n 2 SER SC1 4 -1
5 Q5 3 ALA BB 5 -1
6 TC3 3 ALA SC1 6 0.0
[ bonds ]
3 4 1 0.33 5000
; Backbone bonds
1 3 1 0.350 4000
3 5 1 0.350 4000
#ifdef FLEXIBLE
; Side chain bonds
1 2 1 0.270 1000000
5 6 1 0.270 1000000
#endif
[ constraints ]
#ifndef FLEXIBLE
; Side chain bonds
1 2 1 0.270
5 6 1 0.270
#endif
[ angles ]
; BBB angles
1 3 5 10 127 20
; BBS angles regular martini
1 3 4 2 100 25
3 5 6 2 100 25
; First SBB regular martini
2 1 3 2 100 25
What we see here is that the atom type and bond we specified in the
modification have been applied, and we can also no longer see the BB-SC1 bond
that comes with the normal serine residue (BB SC1 1 0.287 7500
) is no
longer present. In addition, we find the usual backbone/protein interactions.