Pitfalls of the Martini Model.

The computational and conceptual simplifications realized by coarse-grain (CG) models make them a ubiquitous tool in the current computational modeling landscape. Building block based CG models, such as the Martini model, possess the key advantage of allowing for a broad range of applications without the need to reparametrize the force field each time. However, there are certain inherent limitations to this approach, which we investigate in detail in this work. We first study the consequences of the absence of specific cross Lennard-Jones parameters between different particle sizes. We show that this lack may lead to artificially high free energy barriers in dimerization profiles. We then look at the effect of deviating too far from the standard bonded parameters, both in terms of solute partitioning behavior and solvent properties. Moreover, we show that too weak bonded force constants entail the risk of artificially inducing clustering, which has to be taken into account when designing elastic network models for proteins. These results have implications for the current use of the Martini CG model and provide clear directions for the reparametrization of the Martini model. Moreover, our findings are generally relevant for the parametrization of any other building block based force field.

Introduction

Coarse-grain (CG) models play an increasingly important role in computational science and are nowadays a tool as important as atomically detailed models.[1−6] By grouping atoms into effective interaction sites, often called beads, CG models focus on essential features, while averaging over less vital details. This provides significant computational and conceptual advantages compared to more detailed models, allowing for the probing of the temporal and spatial evolution of systems on the mesoscale.

Among the philosophies of CG modeling, we find both systematic (also known as hierarchical) and building block.[2,3,5] CG models developed on the basis of the former, purely “bottom-up” principle focus on the accurate reproduction of the underlying atomistic structural details at a particular state point for a specific system but require reparametrization whenever any condition changes. This translates into a more time-consuming parametrization procedure. Moreover, complex potential forms are often required, which can result in lower performance and thus less sampling. On the other side, building block approaches usually rely more heavily on a “top-down” approach, where macroscopic properties (e.g., thermodynamic data) are used as the main target of their parametrization. Top-down CG models are often cheaper—due to simpler potential forms and only partial parametrization required—and transferable, as the parametrization of the building blocks allows for reusing them for similar moieties in different molecules. However, the structural accuracy of top-down models is limited as the representation of the atomistic detail is suboptimal. The line that separates these two methodological philosophies is, however, thin. Many successful force fields have been developed combining top-down and bottom-up approaches.[3,5]

One important example of the building block philosophy applied to CG modeling is the Martini force field.[7−9] Designed as a model for simulations of lipids and surfactants,[10] this force field has become the most widely used CG model for simulations of biomolecules,[8,11,12] and it is increasingly popular in soft materials science.[13−20] The Martini model mainly relies on a four-to-one mapping scheme, where on average four non-hydrogen atoms are mapped into a CG regular (R) bead. Finer mappings of up to two-to-one non-hydrogen atoms per CG site are employed when the symmetry of the molecules requires it or for ringlike structures. In the latter cases, small (S) or tiny (T) CG beads are employed.[7,21] There exist four main types of particles: polar (P), nonpolar (N), apolar (C), and charged (Q). These types are in turn divided into subtypes based on their hydrogen-bonding capabilities (with a letter denoting the following: d = donor, a = acceptor, da = donor and acceptor, 0 = not involved in hydrogen bonds) or their degree of polarity (with a number from 1 = low polarity to 5 = high polarity). This gives a total of 18 particle types: the Martini building blocks. The “flavor” of each building block is determined by the nonbonded interactions, which are described by a Lennard-Jones (LJ) 12-6 potential:The LJ σ parameter, determining the effective size of the beads, is 0.47 nm for regular interactions. For the smaller sizes, it is reduced to 0.43 and 0.32 nm in the case of S–S and T–T interactions, respectively. The LJ well-depth ε parameters, determining the strength of the interactions between bead pairs, can vary from 5.6 kJ mol–1 to 2.0 kJ mol–1. These values are scaled down by 75% in the case of S–S or S–T interactions. Together, these LJ parameters determine how the building blocks interact with each other, giving rise to the Martini interaction matrix.[7] Nonbonded interactions were parametrized based on thermodynamic data describing the different affinities of chemical groups toward different solvent phases, namely, free energies of transfer between water and a number of organic solvents—octanol, chloroform, ether, and hexadecane–in a top-down approach.[7] Bonded interactions, described by a standard set of potential energy functions common in classical force fields, are parametrized from the underlying atomistic geometry, usually comparing to experimental data or atomistic simulations in a bottom-up approach. This seemingly simple approach, based on a thorough parametrization of the hydrophilicity/hydrophobicity of the building blocks of the model, resulted in a wide range of successful applications in the modeling of (bio)molecular processes.[8]

The parametrization of any force field is performed under a set of specific conditions. Parametrizations are necessarily carried out on a limited set of small systems described by a number of standard parameters such as the range of LJ parameters and bond lengths and assuming a number of simulation settings which specify how the simulations are carried out, such as the treatment of interactions between particles and temperature and pressure coupling schemes. In the case of the Martini force field, the parametrization was mainly carried out using isolated regular beads or linear molecules composed of such beads and a number of standard parameters for the models, such as bond lengths and angles.[7] Moreover, specific settings were employed for the treatment of the interactions between particles, such as the cutoff treatment. The latter settings will not be discussed in the present work, and the interested reader is referred to ref (22) for a recent work discussing these choices in Martini. Overall, the conditions employed during the parametrization allow for more or less freedom, but there are always boundaries.

Here, we investigate cases of pushing the limits of the parametrization of a building block based CG model, with focus on the Martini force field. Given its very wide use, its more modest initial boundaries of parametrization have been pushed to their limits. In particular, section discusses problems arising from the lack of size-dependent Lennard-Jones interaction parameters. We then explore how going too far from the original bonded parameters affects the behavior of the force field, both in terms of solute (section ) and solvent phases (section ). In section , we then demonstrate how bond length distributions can be affected by weak bond force constants, with consequences for the behavior of the model. Finally, section concludes discussing the implications for the use of the current version of the Martini force field and directions for reparametrizations.

Methods

All-Atom and Coarse-Grained Models

The benzene all-atom (AA) models used in Figure are standard GROMOS (53A6)[23] (retrieved from the ATB server[24]) and OPLS[25] models. Standard Martini 2.2 models[7,21,26] (available on the Martini portal http://cgmartini.nl) were used for the solvents considered in Figure .

Figure 1

Effect of lack of size-dependent Lennard-Jones parameters between particles of different sizes on the dimerization. (a) Potentials of mean force for the dimerization of three-particle Martini ring molecules (see inset) described by regular (red), small (green), or tiny (blue) beads and for atomistic GROMOS (cyan) and OPLS (black) benzene models in water. (b) Schematic of the T-solute in R-solvent dissociation (side view): because the solute is seen by the solvent molecules as described by regular beads, a solvent molecule cannot insert between the two solute molecules until the distance between the positions of the beads is at least twice the diameter of a regular bead.

Figure 4

Behavior of enthalpies of vaporization and the cost of creating a cavity in a solvent in the Martini model. (a) Comparison of computed and experimental enthalpies of vaporization. (b) Computed free energy cost of creating a cavity for a P5 Martini particle type in a Martini solvent vs the experimental enthalpy of vaporization of the solvent divided by the number of non-hydrogen atoms present in the molecule. In both figures, rings (described by short bond length models) do not follow the Martini trend. Note that the cost of placing a bead of type P5 is plotted in this case, but results are qualitatively the same for any Martini bead type (see, for example, Figure S8b). Benzene and water are highlighted.

Simulation Settings

A unique set of GROMACS atomistic run parameters was used for the AA simulations. The Verlet neighbor search algorithm was employed to update the neighbor list, and a 1.4 nm cutoff for LJ and for Coulomb (reaction-field) interactions was employed. The Nosé–Hoover[27,28] thermostat (coupling parameter of 1.0 ps) and the Parrinello–Rahman barostat[29] (coupling parameter of 5.0 ps) were employed to maintain temperature (298.15 K) and pressure (1 bar), respectively. Settings for the CG simulations follow the “new” Martini set of run parameters.[22] Specifically, the Verlet neighbor search algorithm is used to update the neighbor list, with a straight cutoff of 1.1 nm. The velocity-rescaling thermostat[30] (coupling parameter of 1.0 ps) and the Parrinello–Rahman barostat[29] (coupling parameter of 12.0 ps) were employed to maintain temperature (298.15 K) and pressure (1 bar), respectively. CG simulation setting files are available on the Martini portal http://cgmartini.nl. GROMACS[31] 2016.x was employed to run the simulations.

Potential of Mean Force Calculations

The Potential of Mean Force (PMF) profiles were obtained from umbrella sampling simulations.[32] The two solute molecules, either an atomistic or a Martini benzene model or single Martini beads, were placed in a box of at least 5 × 5 × 5 nm3 and solvated in water using the SPC[33] and the TIP3P[34] water models in the GROMOS and OPLS case, respectively. The standard Martini water model was used at the CG level.[7] Umbrella windows were spaced 0.1 nm apart along the reaction coordinate, this being the distance between the centers of mass of the solute molecules. In each window, the distance was restrained by applying a harmonic potential with a force constant of 1500 kJ mol–1 nm–2. Each window was equilibrated for 3 ns (1 ns) and then simulated for 500 ns (150 ns) in the CG (AA) case. A stochastic integrator was employed for both AA and CG simulations, while other settings were the same as the general settings described above. The free energy profiles were calculated using the weighted histogram analysis method (WHAM)[35] as implemented in the GROMACS tool gmx wham.

Free Energies of Transfer Calculations

Alchemical transformations were used to compute free energies of solvation ΔG in a solvent S. The solute was solvated in a pre-equilibrated solvent box of size of at least 5 × 5 × 5 nm3. A series of 11 simulations with equally spaced λ points going from 0 to 1 were performed where solute–solvent interactions were scaled by (1 – λ)—from full solute–solvent interactions at λ = 0 to the disappearance of such interactions at λ = 0. A stochastic integrator was employed; simulations were equilibrated for 2 ns, and each λ point was run for 10 ns. A soft-core potential was employed to avoid singularities due to solute–solvent particle overlaps as interactions were switched off.[36] The soft-core potential Vsc, as implemented in GROMACS,[31] has the following formwhere VA and VB are the full LJ potentials in state A (λ = 0) and state B (λ = 1), respectively, α is the soft-core parameter (set to 0.5 by setting sc_alpha in the .mdp file), p is the soft-core λ power (set to 1 with sc_power in the .mdp file), and σA and σB are the LJ radii of interaction. The free energies and corresponding errors were finally computed using the Multistate Bennett Acceptance Ratio (MBAR).[37] The free energy associated with transferring a solute from a solvent S1 to a solvent S2 (ΔG) is then computed as the difference ΔG – ΔG. In the repulsive TI calculations (e.g., Figure c) the solvent–solute attractive interactions are switched off, that is, the solute–solvent Lennard-Jones dispersion constant C6 is set to zero. The free energies obtained for placing such a purely repulsive LJ particle in a solvent phase capture the free energy cost of creating a cavity in that solvent.

Figure 3

Experimental and Martini partitioning behavior of molecules upon removal of aliphatic carbon atoms for the same chemical functional group: short bond length solvent (benzene) to water case. (a) The benzene → water free energy of transfer for a molecule is plotted against the difference between the same free energy and the one for the corresponding molecule (same functional group) where 1 (green circles), 2 (blue squares), and 4 (purple triangles) aliphatic carbon atoms have been removed. Fits are also shown for the various data sets. Experimental data points (filled) are complemented with COSMO-RS predicted (empty) free energies from ref (43)—see Table S3. (b) The benzene → water free energy of transfer for a Martini two-bead “standard” molecule (i.e., 4 atoms-to-1 CG site mapped molecules with a bond length of 0.5 nm) is plotted against the changes of the free energy of transfer upon bond length reduction to 0.4 nm (green circles), 0.3 nm (blue squares), and 0.2 nm (purple triangles), corresponding to the removal of 1, 2, and 4 aliphatic carbon atoms. Fits are also shown for the various data sets (solid lines). (c) Same plot as (b) but considering only the repulsive component of the LJ interaction between solvent and solute (i.e., only the LJ repulsive constant C12 is nonzero, see also the Methods section); this is approximated as the cost of creating a cavity in the solvent.

Enthalpies of Vaporization Calculations

The enthalpy of vaporization (ΔHvap) has been computed according towhere Ugas and Uliq are the total energies (per mole) of the gas and liquid phase, respectively. The gas phase is approximated as one molecule in a large (7 × 7 × 7 nm3) empty simulation box, and the liquid phase is approximated as an equilibrated box of dimensions of about 5 × 5 × 5 nm3. Gas (liquid) phase simulations were performed in the NVT (NPT) ensemble at 298 K (and 1 bar).

1:1 Mixture and Icosahedron System Setup

The 1:1 mixture of dodecane and dodeca-2,5-diene was set up using the gmx insert-molecules tool of GROMACS. 350 molecules each were added to a rectangular box of 3.5 × 3.5 × 20 nm3. The starting configuration of the eight icosahedrons in a cubic box (8.5 × 8.5 × 8.5 nm3) was generated using the gmx insert-molecules tool to add the eight icosahedrons and the gmx solvate tool to solvate the system with 4634 CG water beads. Both systems were energy minimized (steepest descent, 500 steps), equilibrated for 2 ns (time step of 20 fs), and simulated for 500 ns (time step of 20 fs). A leapfrog integrator was employed.

Polyleucine System Setup

To embed the nine polyleucine peptides in the 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) bilayer, the program insane.py was employed.[38] The peptides were placed on a cubic grid at a distance of 5 nm resulting in a bilayer patch of 15 × 15 nm2 and an overall box size of 15 × 15 × 10 nm3. The POPC bilayer consisting of 567 lipids was solvated with 14211 CG water beads, and 0.13 M NaCl was added after neutralizing the system. The system was then energy-minimized (steepest descent, 500 steps), equilibrated for 500 ps (time step of 10 fs), and simulated for 15 μs (time step of 20 fs). A leapfrog integrator was employed, and reaction-field was used for Coulomb interaction (cutoff 1.1 nm)—as this system contains charged particles, following the “new-rf” set of Martini run parameters.[22] Two different polyleucine (LYS2-LEU26-LYS2) protein models were used: a standard Martini 2.2 model without any elastic network and a Martini 2.2 model with an elastic network of GROMACS bond type 1. Bond type 1 in GROMACS means that the nonbonded interactions between the connected beads are excluded. In both protein models, identical regular bonds, angles, and dihedral angles are applied; their only difference is the elastic network. This allows us to exclude any changes due to different bonded parameters. However, this is a setting for the elastic network which is not commonly used when simulating proteins with Martini and an elastic network. There are two elastic network options commonly used: Martini 2.2 with an elastic network of GROMACS bond type 6, which does not exclude the nonbonded interactions, or the ElNeDyn[39] model, which uses GROMACS type 1 bonds in the elastic network but entails in addition different definitions of bonded interactions. The error in Figure c was estimated as the standard error of the mean when the last 10 μs of the simulations were split into blocks of 1 μs and analyzed separately.

Figure 5

Effect of weak bond force constants on the behavior of Martini systems of increasing complexity. (a) Relative number of DOD-DODE contacts (number of DOD-DODE contacts over the total number of contacts made by DODE molecules) in a 1:1 DOD:DODE mixture as a function of the force constant used in the two bonds of a three-bead DODE molecule. The corresponding effective bond length distributions for such bonds are shown in Figure S9a. The insets show renderings of the DOD(cyan):DODE(gray) mixture for two data points. (b) Cluster size distribution for a simulation of eight icosahedrons (inset) described by P4 Martini beads in water (also described by P4 beads) as a function of the force constant used for the six bonds on the surface of the icosahedrons (red bonds in the inset)—and (d) corresponding effective bond length distributions for such bonds. The force constant is varied from 1250 to 500 kJ mol–1 nm–2. (c) Cluster size distribution for a simulation of nine polyleucine transmembrane α-helical peptides embedded in a POPC bilayer modeled without (black line) and with (red line) an elastic network using a common force constant of 500 kJ mol–1 nm–2; results with the ElNeDyn model (blue line), using the standard force constant of 500 kJ mol–1 nm–2, are also shown (see also the Methods section). The distribution is computed over the last 10 μs of a 15 μs long simulation.

Results and Discussion

Differences in Bead Sizes: The Desolvation Problem

Different Bead Sizes Can Lead to Artificial Free Energy Barriers

Mixing different particle sizes without introducing also mixed resolution LJ parameters can lead to artificial free energy barriers. This is the case in Martini when interactions between small or tiny and regular beads take place. As described above, the LJ σ parameter in the case of S–S interactions is reduced to 0.43 nm, from the 0.47 nm used for R–R interactions. At the same time, the LJ ε between two S-particles, εS–S, is also reduced by a quarter as compared to the interaction between two regular Martini particles, i.e., εS–S = 0.75 εR–R. In the case of T-beads an even smaller σT–T of 0.32 nm is used, while no εT–T scaling is applied. However, for simplicity, the LJ ε and σ for R–S (R–T) interactions are kept the same as the ones for R–R interactions. Therefore, S- or T-beads are seen by R-beads as regular particles, while they interact with other S- or T-beads with the reduced σ values. We find that this leads to the formation of artificial free energy barriers. This can be observed in, for example, potentials of mean force (PMFs) of dimerization of molecules described by S- or T-beads solvated in an R-bead solvent. Such a case is shown in Figure , where PMFs of dimerization of two Martini three-particle ring molecules in water are shown for the three different bead sizes available in the force field. The PMFs have been computed by umbrella sampling, as described in the Methods section.

The LJ ε value for the self-interaction of the solute molecules is kept constant in the three cases (we chose the scaled down intermediate level IV, i.e., 2.625 kJ mol–1), so as to exclude its effect. Bond lengths between the ring particles are also constrained in the three cases and are set to 0.27 nm, the bond distance used in the standard Martini benzene.[7] Atomistic PMFs, computed employing the GROMOS (53A6)[23] and OPLS[25] force fields, are also shown in Figure a for the dimerization of benzene molecules, which the S-ring model may be taken to represent. It is very evident from the plot that, going from the R- to the T-ring, an energy barrier arises at around 0.8 nm, while no such barrier is present in the atomistic PMFs. The barrier increases as the difference in size between the solvent and solute beads increases.

Lack of Size-Dependent Cross Interactions

We rationalize the appearance and increase of the barrier by looking at a simple picture representing the system, a schematic of which is reported in Figure b. Note that, despite referring to the T-rings solvated in water, the following description applies generally to all (T-, S-, R- and atomistic) systems, as a barrier is present in all cases. However, using the R–R LJ parameters for the R–S and R–T LJ cross interactions artificially increases the height of the barrier in the case of the S- and T-rings. When the two T-rings are in close contact—which happens at about 0.36 nm, i.e., the diameter of a T-bead—the interaction is favorable, and the free energy is at its absolute minimum on the profile shown in Figure a. As the two T-rings are pulled apart, a cavity starts to form between them. This cavity cannot be filled by any solvent molecule until the distance between the two solute molecules becomes equal or larger than the diameter of interaction dictated by the solute–solvent LJ σ parameter (i.e., the σR–T parameter in this case). This translates into an energy barrier, which is generally observed in conjunction with dimerization.[40] However, for a solvent R-bead to make its way between the two T-systems the distance between the positions of the beads must be not only 0.89 nm (that is, the diameter of a T-bead, 0.36 nm, plus the diameter of a R-bead, 0.53 nm) but instead 1.06 nm (that is, twice the diameter of a R-bead), as the T-systems are seen as composed by regular beads by the solvent particles. This translates into the formation of a cavity which is larger than what it would be if the σR–T were tailored for R–T interactions (e.g., one could take the arithmetic or geometric average between the σR–R and σT–T; such a choice is found to remove the artificial free energy barrier, as can be seen in Figure S1). This larger cavity thus has an associated increased cavity cost which leads to the artificially higher free energy barrier observed for the S- and T-systems. The barrier increases the larger the mismatch between the solute–solute and solute–solvent σ parameters is.

Generality and Consequences

The formation of an artificial energy barrier in dimerization has been shown for the case of a prototypical CG ring molecule. The effect is most obvious in Martini for ring systems that use S- and T-beads. However, the effect is not limited to ring geometries. It also plays a role in the simplest case, i.e., the dimerization of single beads. The PMFs obtained in this case are shown in Figure S2a and demonstrate again the appearance of an energy barrier as the difference in size between solute and solvent increases—only smaller, as fewer particles are involved. While free energy barriers are usually observed in conjunction with dimerization,[40] the increase of such free energy barriers in Martini is a direct consequence of the lack of size-dependent cross interactions. It is thus a consequence of the design of the model itself one should be aware of. In the T-case, the situation is evidently problematic, and this effect is very noticeable in the stacking and base-pair PMFs of the Martini DNA bases.[21] It should be noted that in case of small S-bead solutes, up to a few particles, the effect is relatively mild (compare the green curve to the atomistic case in Figure a). However, as the number of particles in the ring structure increases, the effect increases, as can be seen from the PMF of the polycyclic aromatic molecule pyrene (Figure S2b).

Solutes with Short Bond Lengths: Effects on Oil/Water Partitioning

We now look at the effects of bonded parameters on the behavior of the Martini CG force field. In particular, we investigate the robustness of the Martini model upon changes in the bond lengths which connect the various building blocks. To this end, we systematically study the effect that varying the bond lengths has on the reproduction of the main experimental parametrization target of the Martini model, i.e., the partitioning behavior of molecules. This is done by comparing changes in experimental and computed free energies of transfer (ΔGtransfer) upon changes in bond lengths. It is useful to define the change in the free energy of transfer (ΔΔGtransfer) as followsthat is, we can divide the change in the free energy of transfer in contributions due to solute–solvent and solvent–solvent interactions. In this work, we discuss uncharged systems, so the interactions involved in eq are controlled by the σ and ε LJ parameters (eq ) associated with the solute and solvent particles. Bond lengths, along with the LJ σ, determine the density of interaction sites that will be found in the simulation. In turn, the density of interaction sites affects the strength of the interactions between molecules and therefore their thermodynamic properties. The LJ ε parameters between the building blocks of the Martini model were parametrized mostly based on single R-beads or molecules composed of linear R-bead chains employing a standard bond length of 0.47 nm. In this section, we vary the bond lengths of the solute molecule, while using Martini solvents either consisting of single R-beads (like water) or described by models composed of linear R-bead chains with standard bond lengths of 0.47 nm (like hexadecane). Thus, we first look at the impact of shortening bond lengths on the ΔΔGsolute–solvent of eq , while using standard Martini solvents and thus well-calibrated solvent–solvent interactions.

Experimental Behavior Corresponding to Shortening Bond Lengths

Before describing the results, it is instructive to consider what changing a bond length in a CG model means in terms of the actual molecules and what behavior(s) should be thus captured by the model. Shorter CG bond lengths arise when the number of atoms mapped with the same number of beads is lower (e.g., when a two-bead model describing octane is adapted to represent heptane), or when the molecule is branched or cyclic (e.g., going from octane to tetramethylbutane or dimethylcyclohexane). Here, we focus on studying the partitioning behavior of molecules upon removal of aliphatic carbon atoms. The behavior of fully branched and cyclic organic compounds with respect to their corresponding linear isomers is shown in Figure S6. We gathered a large set of partitioning data[41] from which to extract experimental trends. The data are plotted in Figure a and show how the hexadecane → water free energy of transfer (ΔGHD→W), for the same chemical functional group, changes upon removal of aliphatic carbons (for an example, see Figure b; more examples can be found in the Supporting Information, Table S1). ΔGHD→W has been chosen because it comprises the two prototypical extremes of hydrophobicity and hydrophilicity, and because there are numerous experimental data points available. It is evident that the hydrophilicity of molecules increases upon reduction of their size, i.e., when removing carbon atoms. This is to be expected, given the higher cost of creating a cavity in water as compared to hexadecane[40] which translates into a higher hydrophilicity of the molecule upon size reduction, as the free energy gain in creating a smaller cavity in water outweighs the one in hexadecane. Branched and cyclic molecules in comparison with their linear versions show a similar trend, also related with the size reduction upon branching (Figure S6); in particular, going from a linear to a fully branched or cyclic isomer is equivalent to the removal on one aliphatic carbon atom in terms of ΔGHD→W (Figure S6).

Figure 2

Experimental—(a), (b)—and Martini—(c), (d), (e)—partitioning behavior of molecules upon removal of aliphatic carbon atoms for the same chemical functional group. (a) The hexadecane → water free energy of transfer (ΔGHD→W) for a molecule (e.g., octane) is plotted against the difference between the same free energy and the one for the corresponding molecule (same functional group) where 1 (green circles, e.g., heptane), 2 (blue squares, e.g., hexane), and 4 (purple triangles, e.g., butane) aliphatic carbon atoms have been removed. Fits are also shown for the various data sets. Data points are available in the Supporting Information, Table S1. (b) Summary of how molecular properties change upon molecular size reduction, including the example of the change of ΔGHD→W for octane upon removal of 1, 2, and 4 carbon atoms. (c) The ΔGHD→W for a Martini two-bead “standard” molecule (i.e., 4 atoms-to-1 CG site mapped molecules with a bond length of 0.5 nm) is plotted against the changes of the free energy of transfer upon bond length reduction to 0.4 nm (green circles), 0.3 nm (blue squares), and 0.2 nm (purple triangles), corresponding to the removal of 1, 2, and 4 aliphatic carbon atoms (that is, to 3.5-to-1, 3-to-1, and 2-to-1 atoms-to-CG sites mapping schemes, respectively). Fits are also shown for the various data sets (solid lines). (d) Solute–solvent radial distribution functions (RDFs) for a Martini two-bead molecule (in this specific example using CG particles of C3 type, but other bead types lead to the same result) solvated in water, hexadecane, and benzene (from left to right) having two different bond lengths: 0.50 nm (gray) and 0.20 nm (purple). (e) Schematic of how too short a bond length affects the solvation shells and thus the interaction between the solute molecule and its environment.

Overall, the main effects of reducing the size of a molecule can be summarized as shown in Figure b: smaller molecules interact less with the environment and possess reduced solvent accessible surface area; due to their smaller size, their solvation comes with a lower ΔGcavity in any solvent; due to the high cost of creating a cavity in water, a greater discount is obtained on the ΔGcavity in water upon size reduction, which makes smaller molecules more hydrophilic. A quantitative empirical observation can also be extracted: hydrophobic molecules get from 3.0 to 3.5 kJ mol–1 more hydrophilic for each aliphatic carbon atom removed, while hydrophilic molecules get a somewhat smaller free energy gain (2–2.5 kJ mol–1). That the proximity of an aliphatic carbon atom to a polar group alters its hydrophobicity is in line with the proximity-based correcting factors often applied within prediction schemes for partition coefficients (logPs).[42]

Effect of Bond Lengths on the Partitioning of Martini Molecules

We now turn to the CG model, to investigate whether it succeeds in capturing the experimental partitioning behavior of molecules upon size reduction. Computed free energies of transfer (via alchemical transformations as described in the Methods section) are plotted in Figure c, in a similar way to what was done in Figure a. In this case, the removal of 1, 2, and 4 aliphatic carbon atoms corresponds to model bond lengths of 0.4, 0.3, and 0.2 nm, respectively. The free energies on the horizontal axis are the ones for all possible Martini two-bead “standard” molecules, i.e., 4 atoms-to-1 CG site mapped molecules with a bond length of 0.5 nm. Note that the standard Martini bond length is 0.47 nm, but there is no significant difference between using 0.47 or 0.5 nm. The removal of 4 aliphatic carbon atoms should be considered an extreme case: it is not realistic to remove 4 carbon atoms and stick to a two-bead model. The Martini approach would dictate that one bead gets removed at that point. However, it is instructive to push this far to extract trends in the behavior of the force field.

Overall, with respect to the experimental behavior, all the molecules with shorter bond lengths become less hydrophilic than what they should (compare to Figure a). For a discussion on the direct comparison of Figure a and Figure c, see the last paragraph of this section. More importantly, the effect is not constant across the whole hydrophilicity scale spanned by the horizontal axis, a trend most evident when looking at the 0.2 nm case. In particular, short bond length hydrophilic molecules (left-hand side of Figure c) gain some hydrophilicity, while hydrophobic molecules (right-hand side of Figure c) eventually get even more hydrophobic than their corresponding “full-size” molecule, the latter case being in qualitative disagreement with the experimental behavior. The same effect is qualitatively observed in all the other pairs of solvents taken into account in the original Martini parametrization,[7] as shown in Figure S3, with variations which correlate with the relative polarity of the two solvents. We will first rationalize our observations and then come back to the comparison with the experimental data.

The “Bond Length Effect”: Increased Solute–Solvent Interactions

Our observations can be rationalized by analyzing the pair correlation functions between solute and solvent molecules and comparing the standard and short bond length cases. This is done in Figure d, where radial distribution functions (RDFs) are shown for a Martini two-bead molecule solvated in three different solvents (hexadecane, water, and benzene). In all cases, an extra solvation shell appears as the molecule reduces in size (see also the schematic representation shown in Figure e). As the extra shell appears, each particle of the solute effectively sees more solvent molecules, and the overall solute–solvent interaction therefore increases. We dub this the “bond length effect”. First of all, this contradicts the conclusion that shortening of aliphatic chains leads to less solute–solvent interactions drawn from the experimental data set (conclusion 1), Figure b). Because different beads interact with the various solvents with different interaction strengths, an imbalance across the horizontal axis arises, which is observed most clearly for the shortest bond length of 0.2 nm (Figure c). Indeed, due to the “bond length effect”, a very hydrophobic solute molecule interacts more strongly with both water and hexadecane upon shrinking, but, due to a stronger interaction with hexadecane than with water, the resulting free energy of transfer contains some added hydrophobicity. The same line of reasoning holds for very hydrophilic solutes, for which the “bond length effect” provides some added hydrophilicity upon shrinking.

This effect, shown for the simplest case (two-bead systems, i.e., one bond length systems), is obviously present also in systems containing more beads. In particular, we examined how the effect scales with the number of beads and with the geometry by computing the behavior of three-bead molecules both in a linear and ring geometry. The results are shown in Figure S4. Notably, the effect is stronger in ring geometries than in linear ones. This can be intuitively explained in terms of the larger overlap between beads present in the case of the ring geometry, leading to a higher density of interaction sites for the same number of particles.

Implications: Short Bond Lengths Make Martini Less Intuitive

The direct comparison of Figure a and Figure c clearly shows a systematic underestimation of the hydrophilicity upon molecular size reduction. In order to compensate this underestimation, a standard procedure in the Martini framework is to switch to a more hydrophilic particle type whenever the number of carbon atoms is reduced. For example, butane is described by a C1 particle, while propane is described by a more hydrophilic C2 bead[7]—see the SI for a discussion. Nevertheless, the “bond length effect” makes Martini less intuitive, as the choice of changing the bead type to a more hydrophilic one upon size reduction starts to depend on whether the molecule is hydrophilic or hydrophobic. While for hydrophilic molecules the “bond length effect” already accounts for some added hydrophilicity, and a change in bead type may be too much, for hydrophobic molecules the “bond length effect” accounts for some added hydrophobicity, and a change in bead type may be even too little. This is relevant because short bond lengths occur in several cases when connecting two big fragments, or many small fragments, to build macromolecules. In many cases, no experimental free energy is available for the resulting macromolecules. One example is the case of connecting multiple side chain models of amino acids to build a protein model. In this case, a “naive” (that is, “just attach it and it works”) building block approach would fail, as the particle type cannot be chosen merely on the basis of the chemical group which needs to be represented. The hidden effect on partitioning behavior of shorter bond lengths blurs the intuitive building block procedure.

Solvents with Short Bond Lengths: Excessive Cavity Cost

With the effect of short bond lengths in solute molecules in mind, we now investigate the effect of a solvent phase constituted by molecules with short bond lengths. We will see how short bond lengths affect not only the solute–solvent but also the solvent–solvent term of eq . The ΔGsolvent–solvent is directly linked to the free energy cost of creating a cavity (ΔGcavity) in the solvent: the stronger the interactions between the solvent molecules, the harder it is for them to make some room for accommodating a solute molecule. We selected benzene as a representative solvent phase described with a model containing short bond lengths, both for the (albeit limited) availability of experimental partitioning data and for its importance as prototypical aromatic molecule.

Benzene/Water Partitioning: Discrepancies between Martini and Experiments

Similarly to the previous section, we investigate what happens when computing the free energy of transfer of a two-bead system as a function of the bond length, but this time for the benzene → water (BENZ → W) case. The results, presented in the same format as Figure c, are shown in Figure b. It is evident that smaller solutes are more hydrophobic, i.e., shortening bond lengths favors partitioning to the benzene phase. The question is now whether this behavior corresponds to what is observed experimentally. Figure a shows how experimental BENZ → W free energies of transfer change when shortening alkyl chains by 1, 2, and 4 aliphatic carbon atoms for the same chemical functional group. Due to limited availability of experimental data, we complemented the data with COSMO-RS[43] predicted free energies of transfer. More compact molecules are found to be more hydrophilic in a similar way to what was found in the case of HD → W free energies (Figure a). Given that Figure b shows the opposite trend, we can conclude that the behavior of Martini BENZ → W free energies of transfer upon reduction of the size of the solute molecule does not agree with experimental observations.

Excessive Cavity Cost in Short Bond Length Solvents

In this section we rationalize the main cause of the discrepancies in partitioning data which involve a solvent phase with short bond lengths. With respect to the HD → W case considered before, we now have two different ΔG contributions which may be affected by short bond lengths at the same time: the solute–solvent and solvent–solvent terms of eq . To disentangle these two aspects, we first performed free energy calculations treating the solute beads as purely repulsive particles (see the Methods section). We consider this as an approximation of the free energy cost of creating a cavity in the solvents, which is part of the ΔGsolvent–solvent term, as an increase in solvent–solvent interactions translates into a higher ΔGcavity in that solvent. Figure c depicts the resulting cavity costs. Comparison to Figure b reveals that the cavity cost is the dominant contribution. For a complete picture, Figure S7c shows the attractive contribution, due to solute–solvent interactions, which significantly contributes only at the extreme bond length of 0.2 nm.

The major role of the cavity cost implies that the key issue is caused by the solvent–solvent interactions. It is insightful to compare computed and experimental[44] enthalpies of vaporization (ΔHvap) for various solvents in order to determine whether solvent–solvent interactions deviate from the Martini trend in the case of short bond length solvents. This is done in Figure a, where data points corresponding to various molecular classes have been depicted with different point symbols. We remark that enthalpies of vaporization are not parametrization targets of the Martini force field and are systematically underestimated due to the limited fluid range of the employed 12-6 LJ potential form.[7,8] From F’igure a, it is evident that while most solvents (such as water and the alkanes) follow the Martini trend, models containing short bond lengths—used to describe ring structures such as benzene within the Martini CG model—deviate from the trend and possess systematically higher ΔHvap. This confirms the issue with solvent–solvent interactions. A few other trends are evident, as highlighted in Figure S8 and associated discussion (Supporting Information).

Both the repulsive contribution (Figure c) and the comparison of experimental and Martini ΔHvap (Figure a) indicate that the discrepancy observed in Figure b is rooted in the ΔGsolvent–solvent. In particular, Figure a shows that solvent–solvent interactions are too strong in short bond length solvents (this despite the 75% reduction in interaction strength between the S-beads—these being the beads used to describe rings in Martini 2). The consequence of this point, as anticipated earlier, is a too high free energy cost of creating a cavity in the short bond length solvent. Experimental data for cavitation free energies are not available, while they can be approximately computed (see the Methods section). However, we expect them to follow trends within Martini that correlate with the ΔHvap data, because both reflect interactions between solvent molecules. We indeed find a strong correlation between experimental ΔHvap data (divided by the number of non-hydrogen atoms) and computed ΔGcavity for a P5 type bead in most Martini solvent models (Figure b). Short bond length models indeed deviate from the trend and present considerably higher ΔGcavity. In particular, we note that the cost of creating a cavity in benzene is higher than in water (Figure b, black data points). Therefore, partitioning to the benzene phase is favored if a solute molecule’s size is reduced, as one gets a larger discount on the cavity cost in benzene than in water. The ΔGsolute–solvent contribution is significant only for the extreme bond length of 0.2 nm (see Figure S7c), and it is responsible for the slope observed in the 0.2 nm data in Figure b. In conclusion, the main reason for the qualitatively wrong behavior of benzene → water partitioning upon shrinking of a solute molecule is the too high cost of creating a cavity in the short bond length solvent benzene.

Short Bond Lengths Caused by Weak Force Constants

We have seen how short bond lengths affect the parametrized behavior of the Martini model. Apart from setting the equilibrium bond distance to a lower value, short bond lengths can also result from weakening the force constant. Previously, a collapse of neighboring backbone beads was observed in the coil secondary structure of Martini 2.1 proteins due to weak force constants.[26] The issue was corrected in the Martini 2.2 protein model by increasing the force constant between directly connected backbone beads in the coil/turn/bend secondary structure from 200/400/500 to 1250 kJ mol–1 nm–2. In this section, we explore the effect of the force constant on the behavior of the model in a more comprehensive way.

Weak Force Constants Impact the Behavior of the Martini Model

To investigate the effect of the force constant, we discuss three systems of increasing complexity. The first one is a 1:1 mixture of two different three-bead molecules, namely dodecane (DOD) and dodeca-2,5-diene (DODE). In the Martini framework they are represented by a C1–C1–C1 and C4–C4–C1 model, respectively. These mimic the lipid tails of dipalmitoyl-phosphatidylcholine (DPPC) and dilinoleyl-phosphatidylcholine (DLiPC), respectively, whose interactions are important for the phase separation in membranes.[45] We decrease the force constant of both bonds, C4–C4 and C4–C1, of the DODE model from 1250 kJ mol–1 nm–2 (the standard Martini force constant for aliphatic chains) to 200 kJ mol–1 nm–2. The system, initially mixed, stays so for force constants above 500 kJ mol–1 nm–2. However, it starts to demix if the force constants become weaker than 500 kJ mol–1 nm–2, as can be seen by the steep decrease in the number of DOD-DODE contacts reported in Figure a as a function of the force constant. Lowering the force constant for the bonds in one of the molecules is thus found to induce phase separation.

For the second test case, we constructed an icosahedron of P4 beads including a central P4 bead and solvated eight of them in water—which is also described by P4 beads in Martini. All 12 outer beads of each icosahedron are connected to the central bead by a harmonic potential with a force constant of 1250 kJ mol–1 nm–2 and a minimum position at 0.47 nm. In addition, six bonds exist on the surface connecting pairs of adjacent corners (red bonds in the inset of Figure b). Thus, each corner is connected to exactly one neighboring corner. In our simulations we varied the force constant of the surface bonds while keeping the force constants to the center constant. The fixed force constant to the central bead ensures that the overall size of the icosahedron stays approximately unaffected, while effective bond lengths on the surface can affect the interaction with the environment. We simulated the association of the eight icosahedrons in water (simulation time 500 ns) with different force constants of the surface bonds. Figure b shows the distribution of solute cluster sizes for the different simulations. When applying a force constant of 563 kJ mol–1 nm–2 or below, the cluster size distribution changes. Our simulations show that short bond length regions on the surface of a given icosahedron interact particularly with such short bond length regions of other icosahedrons. As a result, cluster formation is strongly enhanced.

The third system consists of nine polyleucine transmembrane α-helical peptides embedded in a POPC bilayer (starting as monomers, see Figure S9b). We compare a protein model with an elastic network in the protein backbone using a common force constant of 500 kJ mol–1 nm–2 and a protein model without elastic network. Again, the number of protein clusters is analyzed to investigate if the elastic network impacts the protein–protein interaction. Figure c depicts the distribution of clusters for the last 10 μs of a 15 μs long simulation. See Figure S9c for the number of clusters present as a function of time during the 15 μs simulation. Evidently, the system simulated without elastic network (black line) consists of more and smaller clusters, while in the system with the elastic network (red line) the distribution is shifted toward larger cluster sizes. Comparison to experimental data suggests that the protein model without elastic network is already too sticky. It slightly underestimates the population of the monomeric state.[46−48] For completeness, we also performed simulations of polyleucine with the ElNeDyn[39] model and a standard force constant of 500 kJ mol–1 nm–2 (Figure c, blue line). The clusters appear to be more stable than in the case of Martini 2.2 with elastic network, possibly due to the different bonded parameters used in the Martini 2.2 and ElNeDyn 2.2 models (see also Section 4 of the Supporting Information).

Weak Force Constants Lead to the Formation of Superinteraction Centers

The behavior observed for these three systems, i.e., increased aggregation between models which use weaker force constants, can be rationalized in terms of the effective bond length distributions present in the systems. Such distributions are shown in Figure d for the icosahedron case, but they are qualitatively the same for the two other systems—see Figures S9a and S9d for the DOD-DODE and polyleucine systems, respectively. The equilibrium bond distance of the harmonic bond being fixed to 0.47 nm for all the bonds of Figure d, the effective bond length in the systems differs considerably when force constants lower than 1000 kJ mol–1 nm–2 are used. The weaker the force constant used, the shorter the effective bond length.

Weak force constants let the CG beads within a molecule get too close; when that happens, their LJ interactions with a third bead in the surrounding add up. This increased interaction with the environment results in the creation of a superinteraction center, that is a region with high density of interaction sites—analogous to the situation described in Sections and 3.3. However, for the creation of such a superinteraction center not only the equilibrium position of the bonded potential (as seen in Sections and 3.3) is of importance but also its force constant. If this force constant is too weak, the bonded interaction cannot compete with the nonbonded force. Their imbalance enables the bonded beads to approach closely resulting in a distance distribution which is not centered at the minimum position anymore (Figures d and S9a).

Outlook

We have seen how the lack of size-dependent Lennard-Jones parameters in the Martini model can artificially increase the barrier in free energy profiles of dimerization. This effect is larger, the larger the mismatch between the solute–solute and solute–solvent Lennard-Jones σ parameters and the bigger the solute molecules. We have then investigated in detail the effect that the use of bond lengths shorter than the ones used during the parametrization has on the behavior of the Martini force field. Shortening bond lengths increases the density of CG interaction sites and may thus lead to imbalances. In particular, we have seen that shortening the bond length of a solute molecule increases its interactions with any solvent. Because different beads interact with the various solvents with different interaction strengths, the effect is nonlinear and thus unbalances the carefully parametrized behavior of the Martini force field. We have also shown how the use of a solvent phase constituted of short bond length molecules leads to even bigger discrepancies. The enhanced interactions between solvent molecules increases the cost of creating a cavity in the short bond length solvent disproportionately, disturbing the balance between different solvents. Finally, we have seen that a lower limit for the force constant of bonded interactions described by harmonic potentials exists if they entail exclusions, i.e., nonbonded interactions between the bonded beads are not present. If the lower limit is undercut, the harmonic potential cannot compete with the nonbonded potentials which leads to short bond lengths and thus increased interactions.

Implications for the Current Use of Martini 2

Discrepancies between parametrized and observed behavior may arise when systems are rich in molecules containing short bond lengths. A short bond length phase clearly possesses increased interactions both with solute molecules and between the molecules of the phase themselves. Such short bond lengths arise in Martini models when finer mappings are designed. A perfect match with an atomistic bond distribution should be sacrificed in exchange for more reasonable densities of CG interaction sites. Deviations from the parametrized behavior are mostly expected when mixing standard and short bond length systems. A consistent use of short bond lengths for both solute(s) and solvent(s) may reduce the discrepancies observed in properties such as partitioning or mixing due to consistent shifts in overall behavior. However, properties of models rich in short bond lengths may deviate from the overall behavior of the force field.[14,49−51] Moreover, as soon as both standard and short bond length models are present, short bond length molecules will interact predominantly with other short bond length molecules. This effect may be partly responsible for the need of “custom” beads that emerged when modeling a number of polymers.[52,53] Such systems rely heavily on S-beads, hence contain short bond lengths, and need to behave properly in both S-beads and regular solvents. This effect may also contribute to the observed stickiness of Martini proteins[54−56] or sugars.[57] While this complex multicomponent problem is not straightforward and affects also atomistic force fields,[58] short bond lengths will be part of the problem in the case of Martini, as both sugars and proteins contain short bonds.

More generally, when dealing with models based on a building block approach, not only the calibration of the fragments but also their connection must be considered carefully. Despite careful calibration of the fragments, their connection can introduce artifacts, as it was shown to be the case in the Martini model. The extensive use of a certain model within a wide and various research community can only be beneficial to the improvement of the model, as such nontrivial effects can be spotted more promptly. In a broader view, this can affect also atomistic force fields based on similar building block philosophies. While there is much less variability between bond lengths at atomistic resolution, a similar role to the one played by bonds within Martini may be played by dihedral angles in atomistic force fields. Moreover, the required accuracy of an atomistic force field is typically higher, and small systematic errors may accumulate and lead to significant deviations at the level of macromolecular interactions.

Directions for Reparametrizations

The findings reported in this work lead to clear paths for improvements of the Martini CG model and should be also taken into account in the parametrization of any other building block based force field. Specifically, size-dependent Lennard-Jones parameters are necessary to ensure balanced interactions between CG interaction sites of different sizes and to avoid artifacts such as increased barriers in dimerization profiles. A reader who has some experience with customization of CG force fields might be attracted by the idea to fix this issue in Martini 2 by using the arithmetic or geometric average as the LJ σ of two differently sized beads. However, we discourage these readers to simply do so because the overall balance of interactions will be disturbed. Instead, a full reparametrization is required. The density of interaction sites is a very critical property of the system. If finer mappings are required due to symmetry or necessity of a description with a higher resolution, well calibrated particles with different sizes should be available. Such bead sizes should probably be calibrated in a way that will lead to correct trends for enthalpies of vaporization (and hence cavity costs) for different resolutions. Ideally, models of the same molecules with different resolutions, e.g., a dodecane molecule mapped with three 4-to-1, or four 3-to-1, or six 2-to-1 atoms-to-CG-site, should give the same enthalpy of vaporization and free energy of solvation (hence cavity cost) and hence mix ideally between themselves. The different resolutions are intrinsically coupled to the bond lengths used in the systems. If short bond lengths are necessary, it is because finer mappings or very branched chemical moieties are being represented. Thus, finer mappings imply smaller beads and shorter bond lengths, while coarser ones imply larger beads and longer bond lengths. If this harmony is not maintained, an imbalance in the parametrized behavior of the model is expected. Lastly, the elastic network approach might be replaced by a Go̅-model approach[59,60] to (i) avoid problems with weak bonds and (ii) allow some folding–unfolding at the same time.

These guidelines have been taken into account in the reparametrization of the Martini CG force field which led to the very recent development of Martini 3.0.[61] However, while some choices, like the use of size-dependent LJ parameters, can be taken into account during the parametrization procedure, others, like the coupling between bead sizes and bond lengths, should be kept in mind, as this cannot be guaranteed by the parametrization procedure itself but only by a careful use of the different bead sizes.