Perspective on the Martini model

The Martini model, a coarse-grained force field for biomolecular simulations, has found a broad range of applications since its release a decade ago. Based on a building block principle, the model combines speed and versatility while maintaining chemical specificity. Here we review the current state of the model. We describe recent highlights as well as shortcomings, and our ideas on the further development of the model.


Introduction
The use of coarse-grained (CG) models in a variety of simulation techniques has proven to be a powerful tool to probe the spatial and temporal evolution of systems on the microscale, beyond what is feasible with traditional all-atom (AA) models. A large diversity of coarse-graining approaches is available, ranging from qualitative, often solvent-free, models to models including chemical specificity. Models within the latter category are typically parameterized based on comparison to atomistic simulations, a so-called bottom-up or systematic multi-scale approach. They are designed to match a set of specific target distributions from the atomistic simulations, using iterative Boltzmann, 1,2 force matching, 3,4 minimization of relative entropy, 5 or conditional reversible work 6 approaches. For an extended overview of these methods, the reader is referred to recent reviews. [7][8][9][10][11] The Martini force field, [12][13][14] developed by the groups of Marrink and Tieleman, has also been derived in close connection with atomistic models, especially for bonded interactions. However, the philosophy of our coarse-graining approach is different. Instead of focusing on an accurate reproduction of structural details at a particular state point for a specific system, we aim for a broader range of applications without the need to reparameterize the model each time. We follow a top-down approach by extensive calibration of the non-bonded interactions of the chemical building blocks against experimental data, in particular thermodynamic data such as oil/water partitioning coefficients. Processes such as lipid self-assembly, peptide-membrane binding, and protein-protein recognition depend critically on the degree to which the constituents partition between polar and non-polar environments.
The overall aim of our coarse-graining approach is to provide a simple model that is computationally fast and easy to use, yet flexible enough to be applicable to a large range of biomolecular systems. Example input files for many systems can be downloaded from http://cgmartini.nl. The first version of the CG force field, with parameters for lipids only, was published by the Marrink group in 2004. The name 'Martini' for the force field was coined in 2007 with the release of version 2.0 for lipids. 13 The subsequent extension to peptides and proteins 14 was released as version 2.1, with recent improvements to version 2.2, the current version. 15 Although the Martini model was originally developed for use with the GROMACS software suite, 16 the general form of the potential energy functions has allowed other groups to implement the Martini model into other major simulation packages such as NAMD, 17 GROMOS 18 and Desmond. 19 Note that the groups of Schulten 17 and Sansom 20 have developed CG protein force fields compatible with the Martini lipid force field, but different from the Martini protein force field.
In this perspective we describe the current state of the Martini model. We highlight its applications, carefully consider its shortcomings, and sketch the road ahead. The rest of this perspective is organized as follows. Below we describe the design principles of the Martini model, and provide a comprehensive overview of the field of applications. Subsequently, the limitations of the model are described. An extensive outlook section ends this paper.
chosen as an optimum between computational efficiency on the one hand and chemical representability on the other hand. Mapping of water is consistent with this choice, as four real water molecules are mapped to a CG water bead. Ions are represented by a single CG bead, which represents both the ion and its first hydration shell. To represent the geometry of small ring-like fragments or molecules (e.g., benzene, cholesterol, and several of the amino acids), the general four-to-one mapping approach is too coarse. Ring-like molecules are therefore mapped with a higher resolution of up to two non-hydrogen atoms to one Martini particle. The Martini mapping of major classes of (bio)molecules is shown in Fig. 1.
Based on the chemical nature of the underlying structure, the CG beads are assigned a specific particle type with more or less polar character. The Martini model has four main types of particle: polar (P), non-polar (N), apolar (C), and charged (Q). Within each type, subtypes are distinguished either by a letter denoting the hydrogen-bonding capabilities (d = donor, a = acceptor, da = both, 0 = none) or by a number indicating the degree of polarity (from 1 = low polarity to 5 = high polarity), giving a total of 18 particle types or 'building blocks'.

Non-bonded interactions
Non-bonded interactions are described by a Lennard-Jones (LJ) 12-6 potential. The strength of the interaction, determined by the value of the LJ well-depth e ij , depends on the interacting particle types i,j. The value of e ranges from e ij = 5.6 kJ mol À1 for interactions between strongly polar groups to e ij = 2.0 kJ mol À1 for interactions between polar and apolar groups, mimicking the hydrophobic effect. The effective size of the particles is governed by the LJ parameter s = 0.47 nm for all normal particle types. For the special class of particles in ring-like molecules, slightly reduced parameters are defined to model ring-ring interactions: s = 0.43 nm, and e ij is scaled to 75% of the standard value. The full interaction matrix can be found in the original publication. 13 In addition to the LJ interaction, charged groups (type Q) bear a charge AEe and interact via a Coulombic energy function. Coulombic interactions are screened implicitly with a relative dielectric constant e rel = 15 to account for the reduced set of partial charges and resulting dipoles that occur in an atomistic force field. Note that the non-bonded potential energy functions are used in a shifted form. The non-bonded interactions are cut off at a distance r cut = 1.2 nm. The LJ potential is shifted from r shift = 0.9 nm to r cut . The electrostatic potential is shifted from r shift = 0.0 nm to r cut . Shifting of the electrostatic potential in this manner mimics the effect of a distance-dependent screening. To alleviate some of the limitations of the implicit screening of electrostatic interactions caused primarily by a water model without charges or dipole, polarizable water models were also recently introduced for use with Martini. 21,22 With these models a much lower dielectric constant is used because of the increased explicit screening (e rel = 2.5 for the Martini polarizable water model, 21 e rel = 1.3 for the Big Multipole Water (BMW) model 22 ).
The non-bonded interactions of the Martini model have been parameterized based on a systematic comparison to experimental thermodynamic data. Specifically, the free energy of hydration, the free energy of vaporization, and the partitioning free energies between water and a number of organic phases were calculated for each of the 18 different CG particle types. Martini reproduces the correct trend in free energies of hydration and vaporization. However, the actual values are systematically too high, implying that the CG condensed phase is not as stable with respect to the vapour phase as it should be. The same is true with respect to the solid phase. This is a known consequence of using a LJ 12-6 interaction potential, which has a limited fluid range (see 'Limitation' section below). As long as its applications are aimed at studying the condensed phase and not at reproducing gas/fluid or solid/fluid coexistence regions, the most important thermodynamic property is the partitioning free energy. Importantly, the water/oil partitioning behaviour of a wide variety of compounds can be accurately reproduced with the current version of Martini.

Bonded interactions
Bonded interactions are described by a standard set of potential energy functions common in classical force fields, including harmonic bond and angle potentials, and multimodal dihedral potentials. Proper dihedrals are primarily used to impose secondary structure on the peptide backbone. Improper dihedrals are mainly used to prevent out-of-plane distortions of planar groups. LJ interactions between nearest neighbours are excluded.
To parameterize the bonded interactions, we use structural data either directly derived from the underlying atomistic geometry (such as bond lengths of rigid structures) or obtained from comparison to atomistic simulations. In the latter procedure, the higher resolution simulations are first converted into a ''mapped'' CG (MCG) simulation by identifying the centre of mass of the corresponding atoms as the MCG bead. Second, the distribution functions are calculated for the mapped simulation and compared to those obtained from a true CG simulation. Subsequently the CG parameters are systematically changed in an iterative way until satisfactory overlap of the distribution functions is obtained.
For proteins, the PDB databank has been used as an additional source of atomistic reference geometries. Based on MCG distributions from the PDB, bonds, angles, and dihedral interactions were optimized. Importantly, the bonded parameters depend on the sequence, and are used to stabilize the secondary structure elements of the protein; the lack of directional hydrogen bonds prevents realistic folding at the Martini level of coarse-graining. On top of this, an elastic network can be used to further constrain the protein close to a particular, e.g. native, state. The elastic network approach, named ElNeDyn, 23 has been optimized with respect to atomistic reference simulations.

Simulation parameters
The force field terms are coupled to the simulation algorithms Martini was parameterized for, although there is some flexibility in this regard. Key parameters include the timestep, neighbour list update frequency, cut-off radii for the nonbonded potentials, and the exact form of the switch function. Martini simulations generally are stable with timesteps of up to 40 fs or 20 fs in the presence of rings such as in cholesterol or proteins. Lower values than 20 fs we consider an inefficient use of computer time as errors in the model will strongly dominate errors due to numerical integration time steps. 24 Energy conservation also depends strongly on the use of neighbour lists, but the implementation of these lists differs in different software and for GROMACS in different versions. More frequent updates or the use of large buffer zones removes errors associated with neighbour lists but come at a computational cost as well. Martini is parameterized with a switch function that reduces both LJ and Coulomb interactions to zero at 1.2 nm. The form of this function together with the accurate calculation of temperature turned out to be important in implementing Martini in other software as errors may lead to more serious temperature artefacts. 25 Validation Currently, the Martini force field provides parameters for a variety of biomolecules, including many different lipids, 12,26-28 sterols, 13,29 peptides and proteins, 14,20,30 sugars, 31,32 polymers, [33][34][35][36] nanoparticles, 37-39 dendrimers, 40,41 and more. This family of Martini molecules is constructed using the building block principle, i.e. mapping the underlying chemical structure to the corresponding flavor of CG particle types, and stringing them together to reproduce the overall topology of the target molecule. The basic assumption underlying Martini is that the carefully parameterized properties of the individual beads are transferable to the molecule as a whole. This basic assumption requires validation, which may come either from comparing to more detailed atomistic simulations or to experimental data.
Recent validation examples of the first category are: partitioning of amino acid side chain analogues (SCAs) across the membrane, which shows profiles matched to within 1-2 kT for most SCAs; 42 dimerization of SCAs in solvents of different polarities, reproducing most atomistic dimerization free energies to within 1 kT; 43 and the potential of mean force (PMF) between a pair of fullerenes in different solvents, reproducing the atomistic profiles. 39 Examples of validation with respect to experimental data are: the area per lipid, typically reproduced to within 0.1-0.2 nm 2 or experimental accuracy for many types of lipid membranes; 12,13 the ternary phase behaviour of lipid mixtures, showing phase diagrams in semi-quantitative agreement; 44 the relative binding free energy of a systematic series of pentapeptides to the water/ lipid interface, 45 consistent with the experimental hydrophobicity scale derived from these peptides; dimerization free energies of transmembrane (TM) helices within the range of experimental values from FRET data; 46 2 H-NMR quadrupolar splittings of WALP peptides, 47 reproduced to a better extent than atomistic simulations; and the structure of the glycophorin A dimer, 48 compared to NMR data.
Many more critical tests have been performed in several of the papers listed below in the Application section. Sometimes it turns out that the set of standard Martini beads is not sufficient to reproduce the desired accuracy, and further optimizations are made. One of the advantages of Martini is its limited set of parameters, which makes it relatively easy to adjust or optimize the interactions. We quote the original characterization of the Martini model: 13 ''how a few simple ingredients (read: chemical building blocks) can be endlessly varied to create a complex palette of taste''. Cheers.

Martini based applications
The first applications of the Martini model, dating back 10 years, concerned the self-assembly 49 and fusogenicity 50 of small lipid vesicles. Since these first studies, the list of applications has grown dramatically, reflecting the flexibility and transferability underlying our coarse-graining protocol. Currently, applications can be grouped as: characterization of lipid membrane properties; lipid polymorphism; protein-lipid interplay; membrane protein oligomerization; self-assembly of soluble peptides and proteins; protein conformational changes; binding and pore-formation in membranes by membrane active peptides; design of drug and gene delivery systems; structure and dynamics of lipoprotein particles; membrane fusion; compression and expansion of monolayers; self-assembly of surfactants; characterization of carbohydrate based systems; structure and dynamics of polymers; and interaction of nanoparticles with membranes. From this overview it is apparent that lipid membranes are central in many applications. Given the origin of the Martini force field for simulations of lipid systems this is not surprising. Below we discuss the main application areas to date.

Lipid membrane characterization
Although simple lipid membranes can be modelled at the fully atomistic level, Martini offers three significant advantages given the current state of atomistic simulations: (i) the use of a CG model allows many independent simulations in which state conditions are systematically varied, e.g. the in silico design of robust membranes, 51 tethered or supported membranes, [52][53][54] or correlating lipid type and membrane properties; 55,56 (ii) direct, unbiased, sampling of rare processes becomes possible, such as lipid flipflop and lipid desorption, 57,58 deformation of asymmetric membranes, [59][60][61][62] or the membrane binding and permeation of drugs and amphiphiles; [63][64][65][66][67] and (iii) the Martini model allows simulations requiring large system sizes, e.g. modelling of lipid vesicles 49,68-73 and calculation of membrane bending moduli 74-77 and viscosity. 78 As an example of the latter category, Baoukina et al. 79 simulated the formation of membrane tethers ( Fig. 2A), mimicking micropipette aspiration experiments. The largest systems studied contained close to 4 000 000 CG particles, making it one of the biggest Martini-based simulations so far.

Lipid polymorphism
Another important area of applications of the Martini model is lipid phase behaviour, requiring both large systems due to the collective nature of phase transitions and long simulation times to observe the critical nucleation events. Examples include the formation of gel phases [80][81][82][83][84][85][86] and inverted hexagonal 26,87-89 and cubic phases, 90,91 or transitions between micelles, bicelles, and vesicles. 28 A snapshot of an inverted cubic phase stabilized by fusion peptides, obtained by Fuhrmans and Marrink 91 using a self-assembly approach, is shown in Fig. 2B. An important breakthrough in our ability to model realistic, heterogeneous, membranes was reported by Risselada and Marrink, 44 who simulated the spontaneous formation of liquid-ordered (Lo) and liquid-disordered (Ld) domains in ternary mixtures of saturated and unsaturated lipids together with cholesterol. Probing the structural and dynamical properties of these fluid domains has received a lot of attention, as it is presumably linked to the formation of lipid nano-domains (''rafts'') in real cells. Follow up studies have further explored the properties of these domains in a number of ways, including: systematic evaluation of composition, 29,92-96 studying partitioning of lipids and amphiphiles at Ld/Lo domain boundaries, 97,98 probing conditions for domain registration, 99 observing sorting anomalies of domains under stress, 100 looking into the effect of immobilization on domain formation, 101 and analysis of raft dynamics. [102][103][104][105] A recent review by Bennett and Tieleman covers in more detail simulations of membrane domains and membrane asymmetry, including many Martini simulations. 106

Membrane protein-lipid interplay
Martini has proven to be very useful to probe protein-lipid interactions of membrane embedded proteins. Recent applications can be broadly classified into six categories: (i) predicting binding modes of proteins to membranes as well as membrane adaptation around proteins, as exemplified by many simulations of a variety of peptides 46,[107][108][109][110][111][112] and proteins. [113][114][115][116][117][118] These kinds of simulations can be efficiently done in high-throughput mode using a self-assembly approach pioneered by Sansom and coworkers. 119,120 3D stress profiles around membrane proteins [121][122][123] have been resolved using Martini as well; (ii) simulating protein sorting, notably the propensity of TM peptides and membrane proteins to partition in either Ld or Lo domains. [124][125][126][127][128][129] Domanski et al., 129 for instance, modelled the sorting of TM model peptides (the synthetic WALP peptides) under crowding conditions. They found that, at lipid/protein ratios characteristic of real membranes, the peptides induced Ld/Lo domain segregation (Fig. 3C). Unravelling the complicated dynamics of membrane proteins in homogeneous membranes is a related area of interest; [129][130][131][132][133][134] (iii) identifying specific lipid binding sites on membrane proteins is a rapidly growing area. Current examples include the enrichment of short-tail lipids around OmpA, 135 phosphatidylinositolphosphate (PIP) binding to the pleckstrin homology domain 136 and to inwardly rectifying potassium (Kir) channels, 137,138 cardiolipin binding to mitochondrial creatine kinase (MtCK), 139 cytochrome c oxidase, 140 and cytochrome bc1, 141 cholesterol binding to the serotonin (1A) receptor, 142 enrichment of anionic lipids at the gap junction hemichannel connexion-26, 143 at potassium channels KcsA and chimeric KcsA-Kv1.3, 144 and at TM and juxtamembrane  79 Lipid head groups are coloured grey, lipid tails green. The tether is partly cut to view the inside. (B) Stabilization of an inverted cubic phase by fusion peptides. 91 Fusion peptides are depicted as red cylinders, the water/lipid interface is rendered as a green surface. The space occupied by the lipids is coloured black.
domains of cytokines, 145 specific interactions of phosphatidic acid (PA) and PIP2 with the actin capping protein (CP), 146 and binding of cholesteryl esters to cholesteryl ester transfer protein; 147 (iv) simulating protein mediated membrane remodelling, i.e. the deformation of membranes under the action of proteins, is also suitable for 'Martinization' as proven by an increasing number of applications. A fine example is the large-scale simulation of Davies et al., 148 probing the role of F1F0-ATP synthase dimers in shaping the mitochondrial cristae. Based on these simulations, the authors propose that the assembly of ATP synthase dimer rows is driven by the reduction in the membrane elastic energy, rather than by direct protein contacts, and that the dimer rows enable the formation of highly curved ridges in mitochondrial cristae (Fig. 3A). Other examples in the category of membrane remodelling are simulations of the curvature field induced by a-synuclein, 149 and the shaping of membranes by BAR domains. 150 Studies on (v) membrane-mediated protein clustering and (vi) gating are described in sections below.

Membrane protein oligomerization
Many integral membrane proteins assemble in oligomeric structures in biological membranes. The Martini force field makes it possible to study these self-assembly processes at near-atomic detail over time scales of micro-to milliseconds. The first study on large scale protein assembly, by Periole et al., 151 considered rhodopsin, a prototypical G-protein coupled receptor (GPCR). Starting from random initial positions, the proteins formed linear aggregates due to competition of nonspecific lipid mediated forces and specific sidechain-sidechain interactions at the protein surface. This explained in detail previous experimentally observed linear aggregates. An extensive systematic study on aggregation of model membrane proteins, modulated by hydrophobic mismatch, membrane curvature, and protein class, was performed by Parton et al. 152 An example from this study is presented in Fig. 3B, showing clustering of generalized b-barrel proteins under mild mismatch conditions. Other examples of protein-protein self-assembly simulated with Martini include the dimerization studies of glycophorin A, 48,153 the repressor of primer protein, 23 the c0-subunit of the ATP synthase complex, 154 TM WALP peptides, 46,47,124,129 gramicidin channels, 155 the delta opioid receptor, 156 TM domains of ErbB2 receptors 157 and Fukutin-I, 158 integrin TM helix heterodimers, 159,160 the integrin-talin complex, 161,162 the cohesion-dockerin system, 163 the TM domain of T cell receptor complex, 164 the bacterial stressresponse peptide TisB 165 and the serine receptor Tsr, 166 as well as tetramerization of the TM domain of the M2 channel 167 and modelling of caseinolytic peptidase B (ClpB) hexameric ring structures. 168 The fast sampling speed of Martini also allows computation of the dimerization free energy of membrane proteins; recent studies of GPCRs 169,170 and OmpF 171 reveal specific, favorable, association interfaces.

Self-assembly of soluble peptides and proteins
In addition to membrane proteins, aggregation of soluble proteins is a potentially large field of Martini applications. So far, however, self-assembly simulations have been restricted mainly to small peptides. A number of these studies deal with dipeptides. In a very nice high-throughput study by Frederix et al., 172 the Martini force field was used to screen all 400 dipeptide combinations and predict their ability to aggregate as a potential precursor to their self-assembly. Systems that showed strong aggregation tendencies were selected for longer simulations in which supramolecular structures were formed consistent with known aggregation states of dipeptides reported in the literature (Fig. 4A). Related studies show self-assembly of diphenylalanine peptides into a range of morphologies, 173 and of N-(fluorenyl-9-methoxycarbonyl)-dialanine peptides into hydrogels. 174 Aggregation of larger peptides has also been the subject of a number of Martini studies. Lee et al. 175 addressed the selfassembly of peptide amphiphiles, composed of a 13-residue hydrophilic peptide connected to a hydrophobic tail. In line with experimental data, they observed spontaneous formation of fibres (Fig. 4B). Likewise, self-assembly of amphiphilic peptides into peptide-vesicles (''peptosomes'') 176 or micelles 177 has been simulated. A particular challenging case for CG models such as Martini, which lack directional hydrogen bonds, is the formation of amyloid fibrils. For a realistic simulation of amyloid formation, a more detailed description of the peptide backbone is probably required; 178 more qualitative studies, however, can be performed. For instance, Sorensen et al. 179 modelled the assembly of amylin (20)(21)(22)(23)(24)(25)(26)(27)(28)(29) peptides, preassembled into protofibril fragments. The protofibril fragments were kept together with an elastic network using the ElNeDyn 23 approach. The simulations pointed to an elongation growth mechanism of the protofibrils into fibres (Fig. 4C). Another example is the study on the effect of lipid concentration on the aggregation propensity of the amyloidogenic peptide apolipoprotein C-II (60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70). 180 Martini studies of interactions between soluble proteins are still limited, but do exist, for example on the stability of the spider silk's N-terminal protein domain 181 and of keratin filaments, 182 the stability of engineered nanofibres, 183 the mechanical properties of protein filaments, 184,185 (20)(21)(22)(23)(24)(25)(26)(27)(28)(29) peptide fragments. 179 Phenylalanine side chains, purple; other side chains, yellow; backbone, cyan. and dissociation of gas phase protein complexes. 186,187 A related application of the Martini model is protein adsorption on solid supports 188,189 or at fluid interfaces. 190 Martini could also be used to screen for protein-ligand interactions in principle. Thus far, we are aware of one example, in which peptide binding to the OppA transporter was studied. 191

Conformational changes in proteins
Although secondary structure changes of proteins cannot be modelled in Martini yet, changes in tertiary structure are unrestricted and in principle realistic within the general approximations underlying our CG model. A fine example is the gating process of the mechano-sensitive channel of large conductance (MscL). When the membrane is put under tension, MscL channels undergo significant conformational changes in accordance with an iris-like expansion mechanism and reach a conducting state on a microsecond timescale. 192,193 The gating of an MscL channel embedded in a liposome has been simulated in a large-scale study by Louhivuori et al. 194 After putting the liposome under stress by adding excess water, mimicking hypoosmotic shock conditions, the MscL channel opened to release excess solvent (Fig. 5A). The energetics of the gating process have been analysed in more detail elsewhere. 195 Full closure of the channel has also been achieved in simulations of the reversible gating of MscL by Ingolfsson et al. 196 ( Fig. 5B). Structural transitions of other channel proteins that have been simulated with Martini include voltage-gated potassium channels, 197 the SecY channel, 198 and Ca-ATPase. 199 Furthermore, acyl carrier protein (ACP) substrate-shuttling within the fatty acid synthase reaction chamber has been modelled, 200 providing valuable insight into the role of linkers and crowding in the shuttling process. The dynamics of membrane-bound and soluble cytochrome P450 has also been studied, revealing correlations between opening and closing of different tunnels from the enzyme's buried active site. 201 Another study addresses collective motions of RNA polymerases that might contribute significantly to the conformational transition between the open and closed states. 202 A final example is a Martini study on the dynamics of TM helices of bacterial chemoreceptors, supporting a piston model of signalling. 203

Peptide induced membrane permeabilization
Antimicrobial peptides (AMPs) are still in the spotlight as potential new antibiotics essential to overcome bacterial resistance to conventional antibiotics, while similar peptides are studied as anti-cancer and anti-fungal agents as well as vehicles for drug delivery. At a certain threshold concentration, AMPs permeabilize the membrane, either by forming a discrete pore or by disrupting the bilayer structure. Using the Martini model it is in principle possible to simulate AMP induced membrane poration. A first example was shown as validation example of the protein force field. 14 Magainin peptides stabilized a toroidal pore with a structure similar to those seen with atomistic models. In Martini simulations, spontaneous pore formation typically has been observed only under rather favourable conditions, such as in very thin bilayers, 204 at very high peptide concentrations, 205,206 or with highly asymmetric bilayers. 207 Other attempts to simulate AMP-induced membrane poration revealed potential problems with Martini, such as the formation of dehydrated pores 208,209 or membrane buckling, 210,211 see the section on 'Limitations' below. Stabilization of pores by human islet amyloid polypeptides has also been investigated, 212 as well as pores stabilized by fusion peptides. 90,91,213 Likewise, Martini allows for systematic calculations of the binding and translocation free energy of AMPs 214,215 and other membrane active peptides such as caveolin, 216 but also in these cases the shortcoming of Martini in modelling lipid-lined membrane pores has to be kept in mind. The ability of AMPs in a surface bound state to cluster charged lipids does not suffer from this problem, and is therefore expected to be a real phenomenon well suited for Martini. 217,218 Drug and gene delivery systems Membranes pose a barrier for typical drug molecules and gene fragments and hinder their use in practice. Therefore, design of delivery vehicles that help these molecules to cross the cell membrane is an active area of research in which Martini has found a growing number of applications. A series of simulation studies of polyamidoamine (PAMAM) dendrimers show these nanoparticles are able to porate membranes. 40,41 PAMAMs can be functionalized to optimize drug loading. 219 Other polymerbased delivery systems that have been investigated using the Martini model are poly(gamma-glutamyl-glutamate) paclitaxel (PGG-PTX), a series of semiflexible polymer-drug constructs varying in PTX loading fraction [220][221][222] and ABA type triblock copolymers. 223 The design of self-assembling nanofibres composed of amphiphilic peptide derivatives for delivery of therapeutic peptides has also been addressed. 224 Other examples of Martini-based studies of drug-delivery systems include simulations of encapsulation of propofol in quaternary ammonium palmitoyl glycol chitosan micelles, 225 investigation of the permeability of lysolipid-incorporated liposomes for enhanced permeation of anticancer drugs, 226 and fusion of perfluorocarbon-based nanoemulsions to target liposomes. 227 Considering gene delivery, the structure of lipoplexes (DNA-lipid complexes) has been modelled. 228 A recent large scale Martini application by Leung et al. 229 involves modelling of lipid nanoparticles containing ionizable cationic lipids, enabling therapeutic applications of siRNA (Fig. 6A).

Lipoprotein particles
Low-density and high-density lipoproteins (LDL/HDL) transport cholesterol in the bloodstream and play an important role in the development of atherosclerosis. They consist of a mixture of components including phospholipids, triglycerides, cholesterol, and associated proteins. The Martini model has proven useful to study the structure and dynamics of these particles, starting with the self-assembly of model HDL by Shih et al. 230 using the Martini variant developed by the Schulten group. Subsequently, the structure and dynamics of HDL 231,232 with realistic size and lipid composition have been modelled with standard Martini. The interfacial tension and surface pressure of various lipid droplets have also been investigated. 233 A model for LDL was also developed, by Murtola et al., 234 and consists of an ApoB-100 protein wrapped around a 20 nm thick lipid droplet (Fig. 6B).
The simulations reveal a complex distribution of the lipid constituents in the particle, and point to biologically relevant contacts between the surrounding protein and some of the core molecules. Another type of lipoprotein particle is the complex formed by Saposin A, a protein involved in lipid transfer. Simulation studies show that the protein can solubilize various phospholipids into small discs. 235 Martini based simulations have been helpful to reveal the molecular packing of artificial lipoprotein particles, so-called nanodiscs, which are increasingly used by experimentalists to characterize membrane proteins. Examples are nanodiscs stabilized by amphipathic polymers (amphipols) 236 or engineered proteins. 237

Membrane fusion
Membrane fusion was one of the very first applications of the Martini model 50 and is still a hot topic where simulation studies can resolve details that are not easily probed experimentally. Possible fusion pathways between lamellar membranes and between vesicles are now quite well established, 50,87,238,239 and current efforts are directed to calculate the energetics of the various intermediates [240][241][242][243][244] and to investigate the role of fusion promoting molecules such as PEG. 245 An important new direction is peptide and protein mediated fusion. Baoukina and Tieleman 246 simulated the fusion of small unilamellar vesicles mediated by lung surfactant protein B (SP-B). They found SP-B monomers capable of triggering fusion events, consistent with the experimentally observed fusogenic effect of SP-B, by anchoring two vesicles and facilitating the formation of a lipid bridge between the proximal leaflets. Another recent example is the Martini simulation of neuronal SNARE-mediated membrane fusion by Risselada et al. 247 (Fig. 7A). The simulations reveal that SNARE complexes mediate membrane fusion in a cooperative and synchronized way, requiring at least one single SNARE complex consistent with singlemolecule fluorescence studies. After fusion, the zipping of the SNAREs extends into the membrane region, in agreement with the recently resolved X-ray structure of the fully assembled state. Force transduction of the SNARE protein synaptobrevin has  also been looked into with Martini. 248 Furthermore, wild-type hemagglutinin (HA) influenza fusion peptides were shown to stabilize a highly fusogenic pre-fusion structure, i.e. a peptide bundle formed by four or more TM fusion peptides. 213 Studies on mutant peptides provided a rationale for experimentally observed incomplete fusion pathways. Related studies also point to the ability of HA fusion peptides to stabilize fusion intermediates. 90,91 For a recent review of simulations on fusion, including many Martinibased works, see Markvoort and Marrink. 249

Lipid monolayers
At polar-apolar interfaces, lipids form monolayers that reduce the surface tension. Besides being of fundamental interest for surface science, lipid monolayer collapse is crucial for maintaining low surface tension at the gas-exchange interface in the lungs during breathing. The Martini force field can be used to study the behaviour of lipid monolayers in detail, both for pure lipid systems and in combination with peptides and proteins. Although quantitative reproduction of the pressure-area isotherms at low surface pressure remains challenging, 250 the model is well suited to study properties of monolayers in the liquid-expanded and liquid-condensed phases, 250,251 and to study for instance the collapse of lipid monolayers at high surface pressure. 252 Mixed monolayers have also been studied, e.g. to model the tear fluid lipid layer, [253][254][255] or to study phase transitions 256 and phase coexistence, 103,104 and the effect of the presence of lung surfactant proteins. 257 In a large scale study Baoukina and Tieleman simulated fusion of lipid reservoirs to monolayers mediated by lung surfactant protein SP-B. 258 They found SP-B induced establishment of a lipid-lined connection similar to the hemifusion stalk (Fig. 7B) allowing lipid flow in a surface-tension-dependent manner. These findings are in line with existing hypotheses on SP-B activity in lung surfactant and explain its molecular mechanism.

Surfactant self-assembly
The Martini force field is ideally suited to study formation of micelles, as it enables simulations that are long enough to equilibrate many surfactant solutions. However, the ability to obtain an equilibrated micellar size distribution depends on the type of surfactant; exchange of surfactants between micelles, and fusion and fission of small micelles are only observed for surfactants with a relatively high CMC. For instance, in the first study on surfactants using the Martini model, DPC (dodecyl-phosphatidylcholine) could be observed to self-assemble into micelles, but convergence of the size distribution was very slow and not achieved on the sub-microsecond time scale. 12 Recent Martini simulation studies, involving acyltrimethylammonium chloride, 259,260 dodecyltrimethylammonium bromide (DTAB), 73 cetyltrimethylammonium bromide (CTAB), 261 sodium dodecyl sulphate (SDS), 262 dihexanoyl-PC (DHPC), 263 and pentaethylene glycol monododecyl ether (C12E5), 264 however, achieve longer (multi-microsecond) simulation times and approach equilibrium distributions. Sampling of interesting phenomena, such as the sphere-to-rod transition, 259,261,264 becomes possible. In a study comparing a number of surfactants, 265 the Martini model reproduced experimental CMCs and aggregation numbers for non-ionic surfactants reasonably well. The temperature dependence of these properties, however, is not correct. Reproducing the correct temperature dependence is a known challenge for CG models in general (see 'Limitation' section below). Martini simulations of micelles can also be used to provide starting structures for fine-grained models, as shown in a recent study on lyso-phospholipids, 266 and to provide details on the absorption and desorption process of non-ionic surfactants. 267 Bicelles, composed of double tail and single or short tail lipids, have also been the topic of a number of studies. 28,268,269 Carbohydrates Carbohydrates (saccharides) constitute a fundamental class of biomolecules and are present in a variety of emerging classes of biomimetic materials. The large size of most polysaccharides warrants the use of a coarse-grained model, yet the complexity of carbohydrate physico-chemical properties makes this a very challenging undertaking. Recently, common mono-and disaccharides have been parameterized for Martini, 31 providing a basis for further carbohydrate modelling. Simulations of glucose and trehalose interacting with a DPPC membrane show suppression of gel phase formation, 31 in agreement with the cryoprotective properties of these sugars. Based on the Martini carbohydrate model, oligosaccharides such as amylose and curdlan, 31 and cellulose 32,270 have also been parameterized. Martini parameters for the important class of glyco-lipids were developed as well, 27 including biologically relevant lipids for the thylakoid membrane, signalling lipids such as PI and its phosphorylated analogues PIPn, as well as gangliosides. A pioneering study on multi-component membranes containing the ganglioside monosialotetrahexosylganglioside (GM1) shows formation of GM1 enriched nano-domains that act as protein shuttles between membrane regions with different degrees of order. 127

Polymers
Coarse-graining has been of fundamental importance in the field of polymer modelling. In order to study the large size of polymer chains and the associated slow relaxation processes, a reduction of the degrees of freedom has proven absolutely necessary. Although the Martini model was designed primarily with biomolecular applications in mind, there is no reason why the same philosophy could not be extended to soft matter in general. Indeed, there is a growing number of basic polymer systems for which Martini parameters have been derived, currently including polyethyleneglycol (PEG), 33,271 polystyrene, 34 triblock copolymers polyethyleneoxide-polypropyleneoxidepolyethyleneoxide (PEO-PPO-PEO), 223,272 polyurea, 273 Nafion ionomers, 274 polyester coatings composed of two dicarboxylic acids and a diol, esterified to neopentyl glycol monomers and crosslinked by hexa(methoxymethyl)melamine, 275 polymer nanofibres composed of nylon-6 (polycaprolactam), 276 PAMAM dendrimers 40,41,219 and PEG-conjugated PAMAM dendrimers, 277 and amphipols. 236 Two recent examples are illustrative of the  278 show the effect of nanoparticle (NP) concentration, size, and polydispersity on the behaviour of the polymer matrix. For instance, the presence of NPs leads to elongation of the polymer chains as the radius of gyration is increased (Fig. 8A). Simulations of grafted polystyrene brushes in good solvent (benzene), by Rossi et al., 279 explore the structural and dynamical configuration of the brushes as a function of their grafting density in a high grafting density regime (0.1-0.3 chains nm À2 ) (Fig. 8B). The presence of metastable collapsed states, with free chain ends trapped close to the substrate, is revealed at all grafting densities. These collapsed states are shown to be long-lived, surviving over a time scale of several microseconds.

Nanoparticles
The two-dimensional carbon sheet graphene and its three dimensional cousins fullerene and carbon nanotubes (CNTs) enjoy a huge scientific interest due to their interesting mechanical and electrical properties, with many realised and potential applications. The interaction of these nanoparticles with biological material is of specific importance with respect to their toxicity. To address these concerns, a number of groups have started developing Martini models for fullerene, 37,39,280 graphene, 281,282 and CNTs, 38,283-287 and explore how they interact with lipids and surfactants. For instance, in a breakthrough study by Wong-Ekkabut et al., 37 the lipid membrane was shown to be a very good solvent for fullerenes (Fig. 9A). Clustering in the aqueous phase, the fullerene aggregates spontaneously enter the membrane and disperse on a microsecond time scale. Likewise, as demonstrated in a pioneering study by Titov et al., 281 graphene nanoflakes dissolve into a lipid membrane and adopt a flat orientation in between the membrane leaflets (Fig. 9B).
The affinity of graphene-based particles for lipids is also elegantly shown in a study by Wallace and Sansom 38 on CNTs.
Dispersing lipids or surfactants randomly in the aqueous solution, self-assembly of a CNT-lipid complex is observed (Fig. 9C). In addition, the interaction of other nanoparticles such as gold clusters, 288-292 nanocrystals, 293 and coated nanoparticles 294,295 with lipid membranes has been studied using Martini. As an example, in the work of Lin et al., 288 the interaction of functionalized gold particles with DPPC/DPPG membranes was systematically investigated, showing penetration, disruption, pore formation or wrapping (Fig. 9D) depending on the surface charge of both the nanoparticle and the membrane. In a related study, 296 the aggregation of monolayer-protected gold nanoparticles in various solvents was investigated. We expect many more of such studies in the near future, for instance along the lines of a very recent study of the effect of PEG functionalized phospholipids on nanotube bundling 297 and on morphological characterization of polymer-functionalized gold nanoparticles. 298 Finally, we mention three other non-biological applications, namely nanopore inhibition, 299,300 surface wetting, 301 and physisorption of organic molecules on graphite. 302

Limitations of Martini
Martini, as any other model, has a number of limitations. It is obviously important to know these limitations, both to make sure the model is used appropriately and to further improve the model. Some of Martini's limitations are shared with coarsegrained models at a fundamental level, such as the chemical and spatial resolution, which are both limited compared to atomistic models; a shifted balance between entropy and  enthalpy due to the reduced number of degrees of freedom; and kinetics that are modified in unpredictable ways. Other limitations are more tightly linked to Martini itself and could be improved in the future at the level of coarse-graining inherent in Martini, or by accepting more detailed descriptions in parts of the system of interest. This presents an interesting philosophical problem, as the essence of a coarse-grained model is that it contains less detail than a particular reference model (atomistic models, in Martini's case). While it is clear that adding back details in such a model can improve its accuracy, this also negates (some of) its computational advantages. Below we discuss a number of known limitations of Martini.

Model resolution and accuracy
With its 4-1 mapping and its particular range of interactions, Martini can reproduce the thermodynamics of a large number of different organic compounds. This mapping clearly limits the chemical resolution. Phosphatidylcholine lipids can be used to illustrate this in a biological context. Key properties of such lipids, including melting temperature and bilayer stability and thickness depend strongly on the length of the acyl chains. The lipid dilauroylphosphatidylcholine (DLPC) has two acyl chains with 12 carbon atoms, represented in principle by 3 Martini beads, while the lipid dipalmitoylphosphatidylcholine (DPPC) has two acyl chains with 16 carbon atoms, represented by 4 Martini beads. The common lipid dimyristoylphosphatidylcholine (DMPC) with 14 carbon atoms falls somewhere in the middle. Thus the mapping to specific lipids is not unique. Martini reproduces properties such as lipid bilayer density profiles, but within the limitations of this mapping and the spatial resolution afforded by the size of Martini particles. Similar limitations occur in common biomolecules such as sterols, where small changes in ring structure or tail can have large effects on the thermodynamics of membranes that contain these sterols; carbohydrates, where stereochemistry and other details that are too subtle to easily represent in a coarse-grained model occur; and proteins, where all amino acids have a different representation but it is unlikely subtle differences between similar side chains are accurately represented.
Martini in many cases is comparable in accuracy to atomistic models, particularly in thermodynamics, but it also amplifies a number of limitations in atomistic models, where lack of electronic polarizability in the standard force field is one of the main limitations. Important interactions such as cation-pi interactions in proteins or the strong electrostatically driven interactions between benzene molecules are difficult to include in a coarse-grained force field to reproduce both the strength of these interactions and, more challenging, the peculiar geometries resulting from these interactions. However, the wide use of Martini and extensive testing in the original papers clearly shows the degree of agreement with experiments and atomistic simulations and allows an assessment of whether Martini is accurate enough for a particular application. Recent progress in linking Martini more closely to atomistic simulations through backmapping and through hybrid simulations appears promising in terms of extending the use of Martini simulations to problems that may be outside the current resolution (see 'Outlook' section).

Effective time scale
In atomistic simulations the time scale is well-defined and properties such as water diffusion coefficients are used to judge the quality of the underlying model. In coarse-grained simulations this is not normally the case. Coarse-graining involves modifying the energy landscape to become smoother, which effectively results in more sampling of the energy landscape in a given time period, speeding up the kinetics of the system. This is one of the main advantages of coarse-grained models, but the speed-up is not easily predictable and is not likely to be the same for all degrees of freedom. 303 In Martini the current best estimate of a semi-universal factor of speed-up compared to atomistic simulations is about 4, based on lateral diffusion coefficients of lipids 13 and TM peptides 130 in membranes. Some papers therefore report an effective simulation length, which is 4 times the formal simulation length. However, the use of a global speed-up factor has to be considered with care. In a systematic study by de Jong 304 on a range of compounds, speedup factors based on self-diffusion ranged between 1.2 and 22 comparing Martini to experiments and 0.79 and 17 comparing Martini to atomistic simulations. The speed-up is dependent on the type of molecule: Martini alkanes have a small speed-up compared to experiments or even a reduced speed compared to atomistic simulations, whereas alcohols show large speed-up factors. This can be rationalized by the lack of explicit hydrogen bonds in the CG model. Water, although highly polar, shows a relatively small speed-up (E2.5), but this is due to the bundling of four molecules together. Lipids show a large speed-up, comparable to polar fluids suggesting that their diffusion is mainly controlled by the polar head groups.
These estimates of speed-up factors are based on singlemolecule diffusion, but the underlying mechanism has not been carefully investigated and may be different for tracer diffusion, viscosity, collective diffusion, etc. For instance, the viscosity of standard Martini water has been determined 78,305 at Z = 7 Â 10 À4 Pa s (at 323 K) in close agreement with the experimental value of 5.5 Â 10 À4 Pa s and actually closer to the experimental value than the widely used SPC water model. For lipid bilayers, apparent surface shear viscosity has been calculated by Baoukina et al. 79 and is similar to experimental estimates; for lipid monolayers, the apparent surface shear viscosity was similar to previously reported atomistic simulations but much smaller than measured experimentally for monolayers. 252 However, a direct comparison for surface viscosities is difficult because experimental shear rates are much lower and simulation results are expected to depend on the system size. Baoukina et al. 103 also calculated the effect of the presence of domains on monolayer surface viscosity and found a clear increase in viscosity in the presence of domains or at higher compression levels, in qualitative agreement with experiment. Thus although there is some data on viscosities, more thorough investigations would be welcome. In addition, density fluctuations (collective diffusion) have not been carefully analyzed yet. Another example of comparing effective time scales is the recent study on the dynamics of the onset of phase separation in membranes composed of ternary lipid mixtures. 105 Comparing atomistic and Martini simulations, the authors conclude that, while sharing the structural features of phase separation in the CG model, the onset of demixing for the atomistic model is 40 times slower.
In addition to the effect of a simplified energy landscape, there is also a much smaller source of uncertainty in kinetics due to the choice of masses of the Martini particles. The actual time step used in molecular dynamics depends on the masses of the particles involved, which correspond to real atomic masses in atomistic simulations but to effective masses in coarse-grained simulations. All Martini particles have a mass of 72 amu for full particles or 45 amu for scaled particles, regardless of the underlying atoms. Effectively this scales the simulation time differently for groups with different underlying masses. This could be corrected by using more accurate masses but since the masses do not affect the ensemble obtained from a molecular dynamics simulation this is not critical. As atomistic simulations increase in reach we expect more detailed comparisons of kinetics between atomistic and Martini simulations, which would help calibrate the effective timescale of Martini simulations. One useful application of Martini is to assess the amount of sampling required to study a particular phenomenon, such as lipid mixing, before attempting this with atomistic simulations.

Free energies, enthalpies, entropies
A third category of properties that are fundamentally affected in CG models are thermodynamic properties, particularly the balance between enthalpy and entropy. Martini is parameterized to reproduce accurate free energies, with a reduced number of degrees of freedom compared to atomistic simulations. Reducing the number of degrees of freedom affects the entropy of the simulation system, which is compensated for by reduced enthalpic terms in the model. However, this means that the temperature dependence of coarse-grained models is a priori not correct. It also means that although free energy differences are accurate, a break-down of free energies into enthalpies and entropies may not be accurate. The classical example of this is the potential of the mean force (PMF) between two hydrophobic molecules in solution. MacCallum et al. performed detailed all-atom calculations of the PMF between two hydrophobic helices, showing a deep well in the free energy, despite an enthalpic barrier at a slighter higher separation than the equilibrium conformation. 306 Aggregation was driven by a favorable entropy, related to the structure of water around the hydrophobic helices and the release of water trapped between the two helices at near-contact. Corresponding Martini simulations with the standard water model showed the right free energy profile, but the basis for this profile is almost entirely enthalpic. 307 The polarizable Martini water model shows a bit more water-ordering effects, but still did not reproduce the atomistic simulations. The BMW water model did at least reproduce the correct sign although the features of the entropy/enthalpy decomposition still differed from the atomistic simulations. This water model is not a universal solution due to its interactions with other parts of the Martini force field, but does point out the importance of water structure if one wants to reproduce more accurately the underlying thermodynamics. A recent study compared atomistic simulations with a different CG model in their ability to reproduce different components of the thermodynamic driving force for a ligand binding to a protein. 308 The authors concluded that an anisotropic component to the water model interactions combined with translational entropy was sufficient to retain a reasonable enthalpy/entropy balance and temperature dependence of ligand binding.
Thus in solution, PMF decompositions based on CG models should be used with great caution and in many cases results with both the standard and polarizable water models in Martini will be unrealistic. In hydrophobic solvents or membranes, the story may be different as water ordering is probably one of the least accurate components of Martini, while lipid ordering or hydrophobic solvents are likely to be represented much more accurately. PMFs for dimerization of two fullerene molecules in octane gave good agreement between atomistic and coarsegrained simulations, including decomposition in entropy and enthalpy profiles (Monticelli, personal communication). A study of the aggregation of two TM helices gave thermodynamic parameters, including enthalpy and entropy, in reasonable agreement with experiment, 47 although obtaining these parameters from atomistic simulations is out of reach at the moment. Thus it seems entropy/enthalpy decompositions in apolar environments give at least qualitatively useful results and warrant further testing.
In a number of studies, Baron et al. compared thermodynamic properties of alkanes and lipids in atomistic and coarse-grained simulations, including configurational entropies of alkanes 309 and lipids, 310 as well as a larger range of thermodynamic properties in liquid alkanes. 18 In general, both models showed similar behaviour, although the temperature dependence of alkanes and lipids modelled with Martini was notably weaker than with the atomistic force field or compared to experiment.
Other properties related to strong water interactions and water ordering are also challenging to reproduce in coarsegrained models like Martini. Free energies of hydration and transfer of groups with full charges are so large that they are difficult to fit in the range of interactions Martini is capable of reproducing based on the current LJ and limited Coulomb interactions. Charge-charge interactions in apolar solvents are too weak because the implicit screening in Martini and the size of the Martini particles reduce interactions too much. One way to improve this may be to place charges off-centre to allow a closer interaction, analogous to recent changes in the protein force field. 15 The functional form of the non-bonded potential The use of a LJ 12-6 potential to describe the non-bonded interactions in Martini is, in hindsight, not the best choice.
The steep repulsion leads to overstructuring of fluids compared to atomistic models, as evidenced for instance in the radial distribution functions (RDFs) for simple alkanes. 309 Furthermore, the LJ 12-6 potential has only a limited fluid range. Consequently, Martini water is prone to freezing even at room temperature, although the switch function reduces attraction and increases the melting temperature compared to a pure LJ function. The melting temperature of Martini water is 290 AE 5 K. 13 In confined geometries freezing is more pronounced. A practical, if chemically difficult to justify, solution is the use of 'anti-freeze' particles, water particles with a slightly larger radius. 13 Switching to a different non-bonded interaction potential could, in principle, improve the relative stability of the fluid phase (see 'Outlook' section). A related property that is not accurate in current versions of Martini is the surface tension of the air/water interface, which is much too low. This results in unrealistic behaviour in monolayers at the air/water interface, where lipids do not spread homogeneously at high surface tensions but rather form pores in tighter packed monolayers. 250 It is also apparent in water/oil interfaces, which show properties of consisting of two interfaces (water/air and air/oil) close to each other, rather than a single interface. 251 Using the polarizable water model, 21 these properties improve but only to a limited amount.

Secondary structure constraints
Realistic modelling of conformational changes in proteins and peptides in biomolecular force fields is essential for many problems, but is a major challenge in CG simulation. Martini is unusual among CG protein force fields in that it has rather sophisticated side chains but a very simple backbone, consistent with the overall design philosophy of Martini. However, that means that conformational changes with the simplified description of Martini cannot be accurately reproduced in general. In the current version of the protein force field, we apply secondary structure constraints to maintain extended and helical secondary structures by strengthening bonded interactions to maintain a preset secondary structure. The secondary structure elements are able to move relative to each other, which in general does not lead to accurate protein structures in solution but enables a number of applications of primarily TM proteins like MscL mentioned above. These proteins are constrained by the membrane environment and consist mainly of helices with short loops where details of the conformations are in some cases not important. Situations where proteins or peptides change structure considerably cannot be reliably handled. A few examples that are hindered by this but otherwise seem very suitable for Martini are membrane-binding peptides that may be unstructured in solution but adopt specific conformations when binding to the membrane, and may adopt yet different structures in aggregates or channels; membrane-binding proteins with loops that may fold into helices or unfold from a protected, perhaps multimeric, structure into a membrane-binding state, and other systems with similar types of conformational changes. These changes can be incorporated as assumptions in Martini simulations, of course, but this is typically less desirable. Perhaps a combination of backmapping and coarse-graining steps with atomistic and coarse-grained models may be useful in some of these cases. Further work in this direction would be interesting. In a more experimental approach, we tried to include backbone pseudo-dihedral potentials that reproduce backbone conformations of atomistic simulations. 178 This method enables some new applications but requires atomistic simulations as input and is not generally transferable between sequences or between chemical environments.
Testing the limits: simulations of pore formation We discussed a range of limitations of Martini, most of which have their origins clearly in the choice of approximations, including the resolution of the model and the functional form of the potential. However, Martini can also break down in cases where there is no obvious reason. One example of this is pore formation in lipid bilayers. This is a particularly interesting example because lipid-lined pores occur in a wide range of problems for which Martini has been used or could be used. Pore formation, or more generally defects, in lipid bilayers occur in processes such as permeabilization of bilayers by antimicrobial peptides, transport of polar molecules across membranes, fusion intermediates, and lipid flip-flop from one membrane leaflet to another. These defects involve water penetrating into the membrane, followed by rearrangements of phospholipids. This is an interesting test case for Martini because of the relatively small characteristic size at the level of water molecules and the complex environment of water and polar head groups in a low-dielectric membrane interior. We have studied membrane defects in great detail by atomistic simulations and tried to reproduce this data with Martini. 311 Fig. 10A shows a snapshot of a DLPC bilayer in which we have restrained a DLPC headgroup in the centre of the bilayer. Atomistically this leads to pore formation and an extensive structural defect, but with Martini (Fig. 10C) the structure is very different although the free energy cost of placing the lipid at the centre is similar. To test whether this is due to the simple water model, we repeated the calculation with the polarizable Martini water, but with the same result (Fig. 10B). To test whether the difference is due to the large size of a Martini water particle compared to an atomistic water molecule, we tethered atomistic water molecules together in tightly-bound groups of four, but this gave nearly the same results as the regular atomistic simulations. Increasing the interaction between water and headgroups did give structures closer to the pore observed in atomistic simulations (Fig. 10D), but worsened other important bilayer properties such as the area/ lipid. We will continue to use a system like this for tests of future refinements of Martini.
on the other hand, there is a lot of work to be done. Below we sketch the main developments we foresee in the near future.

The next generation, Martini 3.x
The use of a LJ 12-6 potential to describe the non-bonded interactions in Martini leads to overstructuring as discussed above. To reproduce the 'softness' of real systems a smaller power for the repulsive term appears more appropriate, as in the recently reparameterized CG model of Klein and co-workers. 312 The other drawback of the LJ 12-6 potential is the limited fluid range. Possible solutions are given in the literature, 313,314 including switching to a power 4 attraction or to a different potential such as the Morse potential. For the next generation of the Martini force field, we aim to change the non-bonded form to a softer potential with longer-ranged attraction. We are currently exploring the use of the Mie potential, 315 a generalized LJ potential with exponents n and m. Preliminary results 304 indicate that key properties can be reproduced with Mie 5-4, 5-3, 4-3, and 4-2 potentials in their shifted form. Another major direction forward is the development of a foldable protein backbone. The lack of structural flexibility of the Martini proteins obviously limits the applicability range of the model. While aiming at a CG model capable of folding a protein from its sequence may be overly ambitious, describing the backbone in a more realistic manner would be an important improvement. One solution is the use of a specific set of bonded parameters optimized to reproduce an ensemble obtained from atomistic simulations or NMR experiments, an approach already used to model amyloids 178 and fusion peptides. 91 For a less biased approach, the directionality of H-bonds needs to be restored. Available coarse-grained models for protein folding represent the backbone H-bonding in quite different manners and generally use a specific potential. [316][317][318][319][320] In contrast we want to keep a simple and generic approach. For instance, we have been exploring the possibility of using a fixed dipole to represent the backbone polarity, see Fig. 11. The dipoles interact with other dipoles and particles through regular Coulomb and LJ potentials. Importantly, this description restores the directionality of backbone interactions, and has a number of added benefits: the overall dipole along the axis of an a-helix is represented, and the charges are natural sensors of the dielectric environment. A similar approach was used to model polar side chains more realistically in Martini 2.2 and may prove to be an effective way to model polar chemical building blocks in general.

More Martini molecules
Nucleotides, RNA and DNA, constitute the major class of biomolecules that has not yet been thoroughly parameterized for Martini. Although there is a Martini model available for a small piece of DNA, 321 a more generic approach is required to make the model applicable to all nucleic acids. We are currently developing parameters for single and double stranded DNA, using a systematic approach analogous to our work on carbohydrates. 31 A preliminary model is shown in Fig. 1F. Another class of biomolecules that requires careful parameterization are glycans (oligosaccharides), e.g. peptidoglycans that are the building block of the bacterial cell wall, or the glycans attached to proteins in the process of glycosylation, or lipopolysaccharides that are the major component of the outer membrane of Gram-negative bacteria. The Martini carbohydrate 31 and glycolipid 27 force fields can serve as starting points. Apart from biomolecules, other molecular systems are in principle amenable to the Martini CG approach. We encourage parameter development for many more polymers and nanoparticles, and expect a growing number of applications  involving hybrid molecules; combining bio-and synthetic molecules into hybrid molecules is a booming research field in the era of synthetic biology. And why not a Martini force field for ionic liquids? Using polarized or polarizable CG particles this might well turn out to be a fruitful approach. Attempts in this direction have recently been undertaken. 322

Algorithmic improvements affecting Martini
A number of algorithmic developments are worthwhile to mention as they are likely to improve the accuracy, speed, and stability of Martini simulations in the near future. A first issue is the choice of thermostat. Although Martini was developed with the widely-used Berendsen weak-coupling scheme for temperature and pressure control, recent developments promise efficient temperature control methods that conserve momentum, 323 which may be important in future work on time and length scales where hydrodynamic effects play a role. It remains to be seen how accurately the Martini model can capture hydrodynamic effects; the fact that the solvent viscosity is actually reproduced quite well is promising in that respect (see the discussion on the effective time scale above). Second, to improve speed, the buffered Verlet pairlist scheme recently implemented in GROMACS shows promise. Instead of a switch function, the Verlet scheme uses a straight cut-off in combination with a shift of the potential to zero at the cut-off. Neighbour list searching is much faster with this method. However, as this involves a change in the overall shape of the non-bonded potential, it might affect system properties. Preliminary results (de Jong and Marrink, unpublished) show that key properties of the Martini model remain unaffected when using a slightly reduced cut-off (1.1 nm instead of the current 1.2 nm) in combination with an optimized neighbour list cut-off. These settings result in a speed-up of almost 100%. Third, improving numerical stability is also deemed important, especially in light of high-throughput applications (see below). One common reason for Martini simulations to crash is that dihedrals, formed by four particles, are not defined when three of these particles are co-linear. In atomistic models this does not occur but at the coarser level of Martini this is difficult to avoid. We are working on alternative definitions of local geometry to avoid this numerical problem, making use of combined bendingtorsion potentials or dummy-assisted dihedrals. 324

Mixing Martinis
One of the current challenges in biomolecular simulation is to develop effective multi-scale methods, [325][326][327][328][329] which combine the advantages of atomistic and coarse-grained models. Multiresolution methods can either use a static division as in QM/MM (denoted ''hybrid'' models), or allow particles to change resolution on the fly (''adaptive'' models). Recently we introduced a straightforward scheme to perform hybrid simulations, making use of virtual sites to couple the two levels of resolution. 327 With the help of these virtual sites interactions between molecules at different levels of resolution, i.e. between CG and atomistic molecules, are treated in the same way as the pure CG-CG interactions. As a proof of principle, hybrid GROMOS/Martini simulations of simple liquids like butane and small peptides were shown to behave well. However, in an extension of this work 330 challenges are pointed out in regard to the electrostatic coupling between the two levels of resolution. In particular, a proper description of the interaction between polar AA molecules suffers from the poor short range screening behaviour of the CG solvent. Applications of this multi-scale method dealing with less polar solvents, e.g. a lipid bilayer, are more promising. For instance, the gating of AA mechanosensitive channels in CG bilayers is currently being studied by our group (Fig. 12A). To achieve a quantitatively more accurate method, cross optimization of the interactions between the Martini and the atomistic force field is probably necessary, as has been recently attempted in the PACE force field. 331 A breakthrough implementation of an adaptive multi-scale method is the AdResS (Adaptive Resolution Simulation) approach, developed by Kremer and co-workers. 332 In this method, a transition region allows molecules to pass from atomistic to CG resolution and vice versa as a function of the position of the molecule in the simulation box. The coupling of resolutions is achieved through the use of a thermodynamic force in the transition region that compensates for the chemical potential difference between the two resolutions. We are currently exploring the coupling of Martini to atomistic resolution using the AdResS scheme (Fig. 12B). Within the field of biomolecules, Martini based multi-scale methods appear ideally suited to study e.g. protein-ligand binding, where the active site and ligand are modelled in atomistic detail and the rest of the protein together with the solvent is coarse-grained, or aggregation of membrane proteins where the surface of the proteins is treated at a finer grained level compared to the interior and lipids may increase their resolution in close contact with proteins. Another potentially powerful approach to multi-scale simulations is based on Hamiltonian Replica Exchange MD (HREMD). In HREMD, multiple simulations are run at the same time with systematic differences in the Hamiltonian, in this case the atomistic and CG Hamiltonians. Coupling of CG models to AA models using REMD has been pioneered by several groups, 333,334 but has not yet been used in combination with the Martini model. Not really multi-scale, but worth mentioning as an interesting method that opens a range of new problems that can be studied, is the use of constant pH algorithms with Martini. 335 At the other end of the scale, there is promise in combining Martini with self-consistent field theory, 336 to couple Martini to multi-domain elastic network models, 337,338 and to combine Martini with protein-protein docking approaches. Another scale bridging opportunity is the ability of Martini to yield converged PMFs between proteins. 169 These PMFs may subsequently be used to describe the effective interactions between simplified proteins in a solvent free approach, a method coined Atomistic Resolution Brownian Dynamics (ARBD). 339 In a related line of research, we are currently parameterizing a variant of the Martini force field that is solvent free (''Martini dry'').

Automated high throughput
World-wide the trend remains to invest in massively parallel computing with petaflop clusters currently topping the list of the fastest computers (see http://www.top500.org). Within the foreseeable future, the exascale will be reached, with clusters composed of millions of cores. Even local facilities will approach the petascale soon. This opens the way to simulations of either very huge systems or massively high throughput. The largest systems simulated with Martini to date include systems of a few million particles, e.g. a huge membrane patch from which a tether was pulled (cf. Fig. 2A), totalling almost 4 000 000 particles in a box with dimensions 86 Â 86 Â 65 nm 3 and simulated for 1 ms; 79 a lipid nanoparticle for gene delivery (cf. Fig. 6A) with 1 400 000 particles in total, a box size of 54 Â 54 Â 54 nm 3 and simulated for 10 ms; 229 and a tetrameric row of F1F0 ATPase dimers (cf. Fig. 3A) composed of 4 635 000 particles in a box of 84 Â 336 Â 19 nm 3 and simulated for 100 ns. 148 Examples of high throughput studies are: scoring of the aggregation propensity of all 400 dipeptides based on 400 Â 100 ns simulations of a box of 300 solvated dipeptides; 172 construction of a CG database for the positioning of 91 membrane proteins in a lipid bilayer, based on hundreds of self-assembly simulations; 119 and performing 10 000 separate simulations by using the Folding@Home distributed computing project to model membrane fusion. 239 A recent example of highthroughput simulations is the systematic sampling of proteinprotein interfaces (Fig. 13). In the forthcoming exascale era, and with the aid of the continuous improvement of simulation software, we can expect system sizes to grow to 10 9 particles, timescale reaching seconds, and high throughput studies systematically exploring 1000 s of conditions. To harvest this enormous power, automated workflows should be designed and optimized. Efforts in that direction are already undertaken in a number of groups, for example in the grid-enabled web portal for MD simulations, 340 automated topology builder 341 and automatic embedding of membrane proteins. 342 To facilitate the Martini workflow, we have developed the Martinize script 15 which generates GROMACS topology files for proteins from coordinate input files. We are currently developing scripts to automatically perform and analyse simulations of heterogeneous membranes composed of arbitrary lipids and proteins.

Conclusions
In the ten years since its initial publication, the Martini model has developed from a model for simulations of lipids and surfactants to the most widely-used coarse-grained force field for biomolecular simulation and increasingly in synthetic biology. It has been used in a broad range of applications and been subjected to many critical tests, resulting both in exciting scientific results and clear routes towards further development.