Martini 2.2P generating unstable model

[TobyEdwardKing]

Hi everyone,

I recently have been trying to Martinise insulin using the martinize2 script with the following command:

!martinize2 -f 3i40.pdb -ff martini22p -x cg.pdb -elastic -p backbone -o topol.top -cys auto -ignore HOH

However, when trying to minimise this protein, I get an Fmax of the order of 10^7. This is found on residue A21:ASN.
I have simulated in polarisable water and sodium chloride, tried multiple minimisations in vacuum and with different algorithms, multiple steps of equilibration, etc. All lead to the same result, which is either pressure scaling by >1% followed by a segmentation fault, or a "cannot fix PBC" error and box vectors of -nan.

I have done a series of point mutations of ASN21 with every natural amino acid. The following mutations cause a similar error, or the .pdb produced by martinize2.py produces nan coordinates for the sidechain of residue 21 (and thus cannot be simulated):

GLN, ARG, SER, THR, LYS, HIS, PHE, TRP, GLY, TYR

The following point mutations resolved the high force and are stable enough to simulate fully:

ASP, GLU, ALA, VAL, ILE, LEU, MET, CYS, SEC, PRO

I have tried with pretty much every single human insulin .pdb file I could find, and this happens invariably with all of them. Therefore, it is probably not a bad initial structure that is the problem. ASN21 is adjacent to a disulfide bridge, but none of the other residues near bridges have this issue.

What can I do to check and fix the initial structure after the fact? I would rather not work with a variant of insulin if I can avoid it.

Thank you to @fgrunewald for his help and correspondance so far, and thanks in advance for any help.

[pckroon]

Thanks for raising this. It could be there's a mistake in the 2.2P datafiles, or the way we create the charge dummies that you stumbled across...
Does the martinize2 output give anything intelligent with that command line? Does martinize 1 produce a stable model? Does EM work, and when you equilibrate, do you use berendsen for pressure and temperature coupling?

I have done a series of point mutations of ASN21 with every natural amino acid. The following mutations cause a similar error, or the .pdb produced by martinize2.py produces nan coordinates for the sidechain of residue 21 (and thus cannot be simulated):

We can't generate coordinates yet. So the mutations will only work out-of-the-box for residues that are smaller or the same size as the original. By extension, if the resulting AA is too different, the structure will be... suboptimal. The mutation flag is more a way to conveniently generate the topologies. As you've noticed martinize2 is not very clever about coordinates.

martinize2 provides the flags -write-graph, -write-repair and -write-canon to write out PDB files of intermediate structures at various stages of the process. I think you want the one produced by write-canon since that's the one that gets written right before mapping.

[TobyEdwardKing]

Does the martinize2 output give anything intelligent with that command line? Does martinize 1 produce a stable model?

I get the following output from that command line:

INFO - general - Merging two molecules/chains because there is a CONECT record between atoms 49A-CYS7:SG and 231B-CYS7:SG
INFO - general - Read 1 molecules from PDB file 3i40.pdb
INFO - step - Guessing the bonds.
INFO - general - 1 molecules after guessing bonds
INFO - step - Repairing the graph.
INFO - general - Applying modification N-ter to residue A-GLY1
INFO - general - Applying modification C-ter to residue A-ASN21
INFO - general - Applying modification N-ter to residue B-PHE1
INFO - general - Applying modification C-ter to residue B-ALA30
INFO - step - Dealing with modifications.
INFO - general - Identified the modifications ['N-ter', 'N-ter'] on residues ['GLY1', 'GLY1', 'GLY1', 'GLY1']
INFO - general - Identified the modifications ['C-ter'] on residues ['ASN21', 'ASN21', 'ASN21']
INFO - general - Identified the modifications ['C-ter'] on residues ['ALA30', 'ALA30', 'ALA30']
INFO - step - Read input.
INFO - step - Creating the graph at the target resolution.
INFO - general - Applying modification mapping ('C-ter',)
INFO - general - Applying modification mapping ('C-ter',)
INFO - general - Applying modification mapping ('N-ter',)
INFO - general - Applying modification mapping ('N-ter',)
INFO - step - Averaging the coordinates.
INFO - step - Applying the links.
INFO - step - Placing the charge dummies.
INFO - step - Applying position restraints.
INFO - step - Setting the rubber bands.
INFO - step - Writing output.
INFO - general - Please cite: Yesylevskyy, S O; Schaefer, L V; Sengupta, D; Marrink, S J; PLoS Comput Biol 2010
INFO - general - Please cite: Monticelli, L; Kandasamy, S K; Periole, X; Larson, R G; Tieleman, D P; Marrink, S J; Journal of chemical theory and computation 2008
INFO - general - Please cite: de Jong, D H; Singh, G; Bennett, W D; Arnarez, C; Wassenaar, T A; Schaefer, L V; Periole, X; Tieleman, D P; Marrink, S J; Journal of chemical theory and computation 2013
INFO - general - Please cite: Marrink, S J; Risselada, H J; Yefimov, S; Tieleman, D P; De Vries, A H; The journal of physical chemistry B 2007
INFO - general - 'Twas a woman who drove me to drink, and I never had the courtesy to thank her for it -- W.C. Fields

So there doesn't seem to be anything particularly obvious from that.

Does EM work, and when you equilibrate, do you use berendsen for pressure and temperature coupling?

EM doesn't provide a low enough energy to be stable psat equilibration. I use Berendsen for pressure and V-rescale for temperature during equilibration, but it doesn't progress very far before blowing up, usually getting pressure scale warnings at steps 0-50 of equilibration.

martinize2 provides the flags -write-graph, -write-repair and -write-canon to write out PDB files of intermediate structures at various stages of the process. I think you want the one produced by write-canon since that's the one that gets written right before mapping.

Looking at the -write-canon output, I'm afriad I don't particularly understand what I'm looking at. It seems to be in two halves, one of heavy atoms and the other of hydrogens. The hydrogen half is completely nan. The other half doesn't look particularly problematic.
I'll attach it here as a .txt file for ease of attachment. I could try running with this as a UA forcefield system but I'm not sure what forcefield I would specify to use it.
WRITE_CANON.txt

[pckroon]

Thanks, and sorry for the radio silence (busybusy).

I'm a bit surprised by the following:

INFO - general - Applying modification N-ter to residue A-GLY1
INFO - general - Applying modification N-ter to residue B-PHE1
INFO - general - Identified the modifications ['N-ter', 'N-ter'] on residues ['GLY1', 'GLY1', 'GLY1', 'GLY1']
INFO - general - Applying modification mapping ('N-ter',)
INFO - general - Applying modification mapping ('N-ter',)

I think there's something fishy with the N-terminus of the B chain.
As for the -write-canon output, it's the graph Martinize2 uses right before mapping it to CG. The hydrogens for example were not present in your input structure, and as such M2 assigned them nan coordinates. The most interesting bit are the conect records ;) But, it doesn't show anything suspicious around the N of B-PHE1... (at a glance at least)

Could you run martinize2 once more, this time with the -vv flag and attach its output?
And could you generate the martini2.2p model for this protein using martinize1? Is that stable? If yes, could you attach the files it generates as well?

[TobyEdwardKing]

martini2_insulin-vv.txt
I have attached the output of martini2 with the -vv flag as requested.
Unfortunately, I cannot get the original martinize.py (taken from the cgmartini.nl tutorials) to run in the Google colab notebook I'm using. I suspect this might be due to Python 2 vs Python 3 changes. Is there a way to call it in polyply/ vermouth?
Thanks again for your help.

[pckroon]

hmmmn.
This output is for when you also mutate residue 21 to ALA, which might confound the error. Besides that nothing really jumps out as being suspicious.
I'll spend a little effort on this side to get martinize1 running (i.e. install python2). In the meantime, could you run martinize2 one more time with -vv and -write-canon, without -mutate, and provide the following output files?

• canonical-graph.pdb
• ITP
• CG.pdb
• Screen output from martinize2-vv

[TobyEdwardKing]

Hi,

I've attached the requested as an archive.
Let me know if you need anything else.
vermouth-insulin.zip

[pckroon]

Thanks, I'll dig in to this. In the meantime I managed to get martinize1 running here, and I attached the itp it produced.

python3 martinize1.py -f 3i40.pdb -o m1_cg.top -x m1_cg.pdb -ff martini22p -elastic -p backbone -cys auto -ss dssp_out.ssd 
INFO       MARTINIZE, script version 2.6_3
INFO       If you use this script please cite:
INFO       de Jong et al., J. Chem. Theory Comput., 2013, DOI:10.1021/ct300646g
INFO       Chain termini will be charged
INFO       Residues at chain brakes will not be charged
INFO       The martini22p forcefield will be used.
INFO       Local elastic bonds will be used for extended regions.
INFO       Position restraints will be generated.
WARNING    Position restraints are only enabled if -DPOSRES is set in the MDP file
INFO       Read input structure from file.
INFO       Input structure is a PDB file.
WARNING    Several chains have identical chain identifiers in the PDB file.
INFO       Found 4 chains:
INFO          1:   A (), 171 atoms in 21 residues.
INFO          2:   B (), 240 atoms in 30 residues.
INFO          3:   A (), 12 atoms in 12 residues.
INFO          4:   B (), 23 atoms in 23 residues.
INFO       Removing 12 water molecules (chain A).
INFO       Removing 23 water molecules (chain B).
INFO       Found SS contact linking chains 1 and 2 (0.203893 nm)
WARNING    Merging chains.
WARNING    This may change the order of atoms and will change the number of topology files.
INFO       Merges: [1, 2]
INFO       All chains will be merged in a single moleculetype
INFO       Total size of the system: 51 residues.
INFO       Will read secondary structure from file (assuming DSSP output).
INFO       Writing coarse grained structure.
INFO       (Average) Secondary structure has been determined (see head of .itp-file).
INFO       Checking for cystine bridges, based on sulphur (SG) atoms lying closer than 0.2200 nm
INFO       Detected SS bridge between ('SG', 'CYS', 6, 'A') and ('SG', 'CYS', 11, 'A') (0.202703 nm)
INFO       Detected SS bridge between ('SG', 'CYS', 7, 'A') and ('SG', 'CYS', 7, 'B') (0.203893 nm)
INFO       Detected SS bridge between ('SG', 'CYS', 20, 'A') and ('SG', 'CYS', 19, 'B') (0.203647 nm)
INFO       Created coarsegrained topology
INFO       Written 1 ITP file
INFO       Output contains 1 molecules:
INFO          1->  Protein_A+Protein_B (chains A B)
INFO       Written topology files
WARNING    Bead names of charges in sidechains differ between .top/.itp and .pdb.
WARNING    Using names in topology, as Gromacs does, gives the correct result.
INFO       Note: Cysteine bonds are 0.24 nm constraints, instead of the published 0.39nm/5000kJ/mol.

        There you are. One MARTINI. Shaken, not stirred.


 Why don't you get out of that wet coat and into a dry martini? 
                                                               --Robert Benchley

Protein_A+Protein_B.itp.txt