GōMartini 3: From large conformational changes in proteins to environmental bias corrections

Abstract

Coarse-grained modeling has become an important tool to supplement experimental measurements, allowing access to spatio-temporal scales beyond all-atom based approaches. The GōMartini model combines structure- and physics-based coarse-grained approaches, balancing computational efficiency and accurate representation of protein dynamics with the capabilities of studying proteins in different biological environments. This paper introduces an enhanced GōMartini model, which combines a virtual-site implementation of Gō models with Martini 3. The implementation has been extensively tested by the community since the release of the reparametrized version of Martini. This work demonstrates the capabilities of the model in diverse case studies, ranging from protein-membrane binding to protein-ligand interactions and AFM force profile calculations. The model is also versatile, as it can address recent inaccuracies reported in the Martini protein model. Lastly, the paper discusses the advantages, limitations, and future perspectives of the Martini 3 protein model and its combination with Gō models.

Fig. 1. Mapping and bead chemical types of the Martini 3 protein model.

The colors indicate the main classes of bead chemical types: P (polar, in red), N (intermediately polar, in blue), C (nonpolar, in gray), and Q (charged, in green). Different bead sizes are also indicated, ranging from the bead with the largest radii (regular, no symbol) to small (S) and tiny (T) beads.

Fig. 2. PI(4,5)P₂ binding of the PLCδ1 PH domain studied at the CG Martini level.

A Setup of the simulation box containing a POPC bilayer (gray) with one PI(4,5)P₂ (orange) in the binding pocket of the PLCδ1 PH domain (red) solvated in water (light blue transparent surface). B Magnified view of PLCδ1 PH domain with the green arrows indicating the vectors used to determine the orientation: membrane normal z (left), α-helix_115,129 (middle), and α-helix_15,24 (right). C Potential of mean force for the PI(4,5)P₂ binding of the PLCδ1 PH domain. The protein is modeled with the GōMartini model (red) as well as two different elastic network models of type 6, force constant of 500 kJ/(mol nm²), cutoff 0.8 nm (violet), and of type 1, force constant 700 kJ/(mol nm²), cutoff 0.8 nm (blue). Solid lines represented the mean values; the shaded area the standard deviation from bootstrapping with N = 100. D Probability distribution of the protein-membrane distance evaluated for ten replicas of 2 µs for each protein model. The distance is measured between the PI(4,5)P₂ head group and the center of mass of the protein. Solid lines represent the mean value; the shaded area the standard error of the mean calculated for N = 10 replicas. E Orientation of the PLCδ1 PH domain measured by evaluating the probability distributions of the angles between α-helix_15,24/α-helix_115,129 and the membrane normal z. The colors in (D, E) are the same as in (C). Source data are provided as a Source Data file.

Fig. 3. Unbiased simulations of ligand binding to L99A T4 lysozyme with elastic network and GōMartini models.

A Distribution of the average protein backbone RMSF for the EN (blue) and GōMartini (red) models. B Average RMSF per residue of the protein backbone bead (BB) for simulations performed with EN (left panel) and GōMartini (right panel) models. C Radial ligand-receptor PMFs obtained with benzene using EN (blue) and GōMartini (red) models. D Benzene density around L99A T4 lysozyme obtained from averaging 0.9 ms of CG simulations for EN (left panel) and GōMartini (right panel) models. The blue, cyan, and red isosurfaces can be translated to the free energy values shown at the color map. Source data are provided as a Source Data file.

Fig. 4. Comparison of the effect of elastic network and Gō models in the allosteric pathways of SOD1.

A, B Flexibility of the protein backbone of the WT with the elastic network (A) and GōMartini (B) model. Snapshots were taken every 1 μs. The color scale represents the average backbone RMSF per residue. C, D Change in RMSF between the WT and G93A with the elastic network (C) and GōMartini (D) model. Blue indicates rigidification in G93A. Red indicates increased flexibility. E, F Matrix representation of the integrated absolute difference in the distance distributions between all backbone beads. Results are presented for the EN (E) and GōMartini (F) model. The bottom left triangle represents the full data set. In the top right triangle, only values of > 0.3 are depicted; all other values are colored blue. Source data are provided as a Source Data file.

Fig. 5. Nanomechanics of the RBD:H11-H4 complex studied by GōMartini simulations at different pulling speeds.

A Force-displacement profiles for RBD:H11-H4 complex at v_pull = 5 × 10^-4nm/ps in SMD simulations using the GōMartini model, with and without position restraints of the BB beads of the RBD. The pulling SMD spring constant was set to 600 kJ/(mol·nm²). Data show the mean, and the whiskers are the standard deviation of the mean value obtained from N  =  50 for each case, respectively. B Same as in (A), but the dissociation of the complex is carried out at v_pull = 10^-5nm/ps, and the pulling SMD spring constant was set to 60 kJ/(mol·nm²). Data show the mean and the whiskers are the standard deviation of the mean value obtained from N  =  50. The inset in A and B shows the reference AA SMD data, note that the x-axis shows the distance (D) between the center of mass of groups pulled in AA SMD protocol, whereas in the GōMartini study, the displacement is associated with increase of z-value along the pulling direction. C Structure of the RBD:H11-H4 complex placed in a box of CG Martini water represented as blue beads in the initial bound state with F = 0 pN. The fixed LYS-528 residue in the RBD and the LYS-128 residue in H11-H4 used for pulling are highlighted by red beads. D Structure of the complex at F_max ~ 434 pN. E Magnified view of the last protein segments in contact before the full dissociation of the protein complex at d ~ 6 nm. The structures in (C–E) are taken from a replica simulated with v_pull = 10^-5nm/ps. Source data are provided as a Source Data file.

Fig. 6. Improving GōMartini to match AA models.

RMSF comparison between original GōMartini (with ε_LJ optimized), modified GōMartini (with the removal of Gō interactions in loops), and AA simulations for three proteins: titin I-band (1TIT), glycoside hydrolase (3W0K), and the transmembrane domain of Ist2. Solid lines represented the RMSF mean values; the shaded area the standard deviation (N = 2 for the AA simulations and N = 3 for the CG ones). The mean absolute error (MAE) for loops and structured regions, calculated as the average of the absolute differences between the RMSF values of GōMartini models and AA simulation, is highlighted in the insets. The bottom-right panel presents the flexibility of the protein backbone beads during simulations using the modified GōMartini model for glycoside hydrolase (3W0K). Source data are provided as a Source Data file.

Fig. 7. Improving IDP global dimensions and condensation using Martini 3 with GōMartini model-based interaction rescaling.

A Radii of gyration of the IDP benchmark set of Thomasen et al. Results are compared between the experimental value (blue) to both the native Martini model (purple) and optimized Martini IDP + Gō model with additional bonded parameters (green). For the experimental data, the error bars indicate the experimental uncertainty taken from ref. . B Illustrations of the increase in ensemble dimensions of ACTR comparing (left) native Martini 3 and (right) the final model for IDPs with additional bonded and non-bonded potentials. C, D Illustrative snapshots of a condensate (with improved IDP parameters) and an aggregate (with default Martini 3 parameters) of an artificial IDP (WT20) known to phase separate. E Snapshots of FLssLF peptide systems with varying increases in the strength of the BB-water interactions. Left; native Martini (0% increase); middle: 6% increase; right: 8% increase of protein BB-water interactions. Source data are provided as a Source Data file.

Fig. 8. Improving transmembrane peptide insertion and beta-sheet aggregation.

Tilt angle distributions from simulations of WALP peptides inserted in DMPC membranes using the GōMartini model with either (A) no additional LJ interaction, or (B) an additional LJ interaction between the virtual Gō sites and water beads of ε = − 1.0 kJ/mol; tilt angle states close to 0° correspond to TM configurations, whereas those close to 90° correspond to peripherally membrane-adsorbed ones. C Representative WALP16 configurations, both fully inserted in its preferred transmembrane configuration and in its peripherally membrane-adsorbed state. WALP16 backbone shown in blue, with side chains in white. Membrane phosphate beads are represented in orange. D Normalized average contacts between two RAD16-I peptide beta strands. Solid lines show the running averages of 500 frames, while the shaded area shows the running standard deviation. Simulations were run with a GōMartini model applied between the two chains (red), with an additional LJ interaction between the virtual Gō sites and water beads with ε = − 0.5 kJ/mol (green), and without any structural or interaction bias (blue). E Representative RAD16-I strand configurations, both aggregated and dissociated. Each backbone chain is colored either brown or green, with the side chains colored in a lighter shade of the same color. Source data are provided as a Source Data file.