Two decades of Martini: Better beads, broader scope

The Martini model, a coarse‐grained force field for molecular dynamics simulations, has been around for nearly two decades. Originally developed for lipid‐based systems by the groups of Marrink and Tieleman, the Martini model has over the years been extended as a community effort to the current level of a general‐purpose force field. Apart from the obvious benefit of a reduction in computational cost, the popularity of the model is largely due to the systematic yet intuitive building‐block approach that underlies the model, as well as the open nature of the development and its continuous validation. The easy implementation in the widely used Gromacs software suite has also been instrumental. Since its conception in 2002, the Martini model underwent a gradual refinement of the bead interactions and a widening scope of applications. In this review, we look back at this development, culminating with the release of the Martini 3 version in 2021. The power of the model is illustrated with key examples of recent important findings in biological and material sciences enabled with Martini, as well as examples from areas where coarse‐grained resolution is essential, namely high‐throughput applications, systems with large complexity, and simulations approaching the scale of whole cells.

70s represented each residue with two particles located on Cα and the center of the side chain, a simplification necessary to enable pioneering simulations of protein dynamics. 1 Similarly, simulations of polymers relied on CG representations, going back to the freely jointed chain models of Binder and coworkers in the early 80s. 2 With increasing computational power and software optimization, more detailed all-atom simulations became the standard in the field of biomolecular simulations during the 80s and 90s. However, the hope that, with a continuing increase in capacity, all processes could be captured with such detailed models turned out to be idle. Even today, with access to exascale computing power, all-atom simulations are still limited to spatiotemporal scales less than 100 nm and typically covering few microseconds. [3][4][5] Considering the gap with processes in real life, the need for a CG description is still as urgent as in the early days. [5][6][7][8][9][10][11] The requirement of more sampling speed became (re)apparent in the late 90s, and caused a surge in the development of numerous enhanced sampling methods as well as the re-valuation of the principle of coarsegraining. In this context, the Martini model was conceptualized and evolved over the years to the current "Martinidome" (Figure 1).
The major assumption underlying the Martini model is that small chemical fragments can be represented by specific CG bead types (the building blocks), irrespective of the molecule to which the fragment belongs, and independent of the state conditions. The interactions between the different bead types are modeled with Lennard-Jones (LJ) potentials, calibrated with experimental thermodynamic data, such as partitioning free energies, as major target. Effective bonded interactions are obtained through optimization of bond distributions with respect to all-atom simulations. The Martini  26 model thus combines top-down and bottom-up parameterization strategies. It is important to realize this coarsegraining philosophy of the Martini model is shared with other CG force fields (e.g., SPICA 12,13 and SIRAH 14,15 ), but is distinctly different from bottom-up models that aim at reproducing structural features of higher resolution models with high accuracy. 16,17 The Martini model, in being a generic force field, also sets itself apart from CG models that are dedicated to specific (bio)molecular classes (e.g., oxDNA, 18 UNRES 19 ), or to coarser (supra-CG) models that do not aim to capture chemical specificity. Despite the more limited structural accuracy of Martini models as compared to structurebased CG models, the model compatibility provided by Martini allows for complex systems to be set up more easily. Moreover, structure-based CG models often require complex potential forms for nonbonded interactions, which may result in lower computational performance than LJ-based CG models. With respect to supra-CG models, the fact that they represent effectively all chemical species as beads connected by springs makes them suitable to describe a wide range of physics, but much less suitable to describe specific chemical systems. Hence, Martini has an edge with respect to coarser CG models when aiming to simulate systems that require chemically specific predictive CG modeling.
Looking back over the past two decades, we have seen many exciting developments around the Martini force field, which has slowly matured during these years. Remarkably, different communities appreciate the model for different reasons: Biologists use it to study cellular processes at length and time scales inaccessible with all-atom models 20,21 ; Chemists appreciate the building-block principle and use Martini to cover the vast chemical space in high-throughput applications in material science 22 ; Physicists take the opportunity to test predictions based on mean-field assumptions and other mesoscale models with a molecularly detailed model, 23 and finally, nonscientists use Martini for illustrative purposes in education 24 and even as a basis for creation of art objects. 25 The still increasing citation numbers (Fig 1b) are further proof of the high demand for a CG model such as Martini in modern integrative research.
The remainder of this review is structured as follows. In the next section, we provide a historical overview of the major developments underlying the Martinidome. Subsequently, we review the latest Martini 3 model, highlighting the major strategic choices underlying the model as well as current best practices of using it. Then, we provide key examples from the recent literature to illustrate the use of Martini, divided into six categories that benefit from the power of Martini: (i) applications in the field of biophysics and structural biology, (ii) applications in the bio-nano interface, (iii) applications in materials sciences, (iv) high-throughput applications, (v) systems with high complexity, and (vi) applications connecting to the full-cell scale.

| A BRIEF HISTORY OF MARTINI
In this section, we will guide the reader through a historical perspective on Martini, from the very first lipid models to the current Martinidome, depicted in Figure 1.

| Lipids only
Originally, the Martini force field was developed to simulate lipid systems only. A prototype model (version 1.0) was ready in 2002, now 20 years ago, and applied to simulate lipid vesicle formation and fusion, 27,28 something that was not feasible with atomistic resolution at the time. The completed version of the first Martini model, version 1.4, was published in 2004, 29 building on pioneering work from Smit et al. 30,31 and the group of Lipowsky, 32,33 who already showed that simplistic CG models of surfactants and lipids can be used to simulate their self-assembly into a variety of aggregation states. However, contrary to these models which are physics based, expressed in reduced units, the Martini model is a chemistry-based model with CG beads representing chemical fragments with realistic units of size and energy.
In this first model, four different bead types were considered to distinguish between apolar, intermediately polar, polar, and charged fragments. The basic idea, which was also adopted by the Klein group around the same time to model surfactant systems, 34 was to use Lennard-Jones potentials to capture the effective nonbonded interactions between CG beads, calibrated to reproduce thermodynamic data, together with standard harmonic bonded potentials parametrized to match the conformational flexibility of reference atomistic simulations. This basic underlying philosophy combining top-down and bottom-up approaches has not changed. Using this concept, the Martini model mapped to real lipids, and was able to distinguish, for example, saturated lipids from unsaturated ones, and a phosphatidylcholine (PC) from a phosphatidyl-ethanolamine (PE) headgroup, with structural and thermodynamic membrane properties "semi-quantitatively" matching experimental data. Note that, from the very beginning, the Martini model used a screened Coulombic potential to capture the interaction between charged particles, including those of zwitterionic lipid heads.
In 2007, the name Martini was coined for version 2.0. 35 The model was still limited to lipid-based systems, but with the increased number of bead types (18) and interaction levels, many more lipids could be realistically represented. Labels were introduced to differentiate between hydrogen bond donor and acceptor groups. This more sophisticated version allowed further pioneering studies, for instance on membranes showing liquid-ordered/liquid-disordered coexistence. 36 More recently, the lipid model has also been reparametrized as an implicit solvent model, appropriately called Dry Martini, 37 and coupled to either stochastic rotational dynamics 38 or lattice Boltzmann 39 to capture hydrodynamics. In the meantime, the Martini 2 lipid library has been steadily increasing, with CG models generated for different classes of more specialized lipids, including glycolipids, 40,41 cardiolipins, 42 polymer-modified lipids, 43,44 lipopolysaccharides, [45][46][47] sterols, and hopanoids. 48

| Extensions to other (bio)molecules
With the apparent success of Martini to simulate lipid dynamics, the natural next step was to extend the force field toward (membrane) proteins, resulting in version 2.1 in 2008, 49 which was later improved to version 2.2 in 2013. 50 To model peptides and proteins we faced an obvious limitation of the use of isotropic beads: capturing the directionality of H-bonds, and hence protein folding, would be out of scope. To keep the secondary structural elements in place, specific bonded parameters were calibrated for helices, β-strands, coil, turns, and bends. 49,50 Additionally, an elastic network approach, coined "ElNeDyn," was introduced to keep tertiary structures folded and stable in relation to a specific reference structure. 51 Early successes of the Martini protein model included the oligomerization of G-protein-coupled receptors (GPCRs), showing a membrane-thickness depending aggregation propensity in accordance with experimental data, 52 and the tension-induced gating of the mechanosensitive channel of large conductance (MscL). 53 Note that the groups of Sansom and Schulten developed their own initial protein models compatible with Martini. 54,55 Extension of the model to other major classes of biomolecules followed over the years: carbohydrates in 2009, 56 DNA in 2015, 57 and RNA in 2017. 58 More specific biomolecules were also added to the Martini library, such as protein posttranslational modifications 59,60 and a variety of cofactors and metabolites. 61,62 In the meantime, Martini had found its first applications outside the field of biomolecular sciences, with models for polymers (PEG) dating back to 2009, 63 and nanoparticles (fullerenes) in 2008. 64 One of the strengths of Martini is that, due to the building block approach, Martini models can be combined allowing researchers from the very beginning to study how different classes of biomolecules interact with each other, and with other materials. A pioneering example is the solvation of fullerenes in a lipid membrane environment. 64

| Polarizable versions and other developments
With the increasing number of applications, also shortcomings of the model became apparent. A number of those are inherent to the process of coarse-graining, such as entropy/enthalpy compensation which affects the temperature dependency, the lack of directionality of H-bonds limiting simulations of protein conformational changes or nucleotide hybridization, and the interpretation of the timescale which is nontrivial as a consequence of the nonhomogeneously smoothened potential energy surface on which the CG dynamics take place. 65,66 Another limiting issue is the accuracy of the electrostatic interactions. Although Martini features charged bead types, and may even be used together with common approaches to include the long-range forces such as Particle Mesh Ewald summation, the major problem is due to the lack of nuclear polarizability of Martini water. In standard Martini, water is modeled as an uncharged bead representing four real water molecules, incapable of performing the electrostatic screening of the real solvent. The screening in Martini is represented implicitly, in a distant-dependent manner using a combination of a dielectric constant of 15 and a shifted electrostatic potential.
To improve on this aspect, multipolar and polarizable versions of Martini water 67,68 and polarizable ions 69 have been developed. In addition, a polarized version of the protein force field was parameterized (version 2.2p) to better capture the interactions between charged amino acids. 50 Polarized Martini beads also found their way into other models; in particular, they are proving helpful to stabilize supramolecular polymers. 70,71 Cation-pi interactions were further added to the mix in version 2.3, released in 2020. 72 Improvements were also achieved in terms of computational speed, through optimization of the use of cutoffs and GPUs. 73 Most recently, a constant pH version of Martini, "Martini-sour", was developed, featuring explicit proton-mimicking beads. 74,75 2.4 | The road toward Martini 3 Despite the gradual maturation of the Martini force field, and the increasing awareness of the intrinsic limitations to CGing in general, 17,76,77 it became clear that another, more specific problem was haunting the model. Several research groups reported on the apparent overstabilization of biomolecular interactions, noted for soluble 78 and membrane proteins, 79,80 as well as carbohydrates. 41,60,81 This shortcoming became known as the stickiness problem. Ad-hoc remedies included downscaling of solute-solute or upscaling solute-solvent interactions, but the underlying cause of the stickiness problem could be traced back to the imbalance between the small and regular bead types, in combination with the use of weak bond force constants or bonds between beads that were too short in the description of biomolecules. 82 The stickiness problem, together with the demand for higher chemical resolution has prompted the development of the recently released Martini 3 version, 83 featuring more than 800 bead types and a thoroughly recalibrated interaction matrix. It has now reached a level of what can be called a "generic" force field, applicable to a wide range of problems, not only in biomolecular simulation but also increasingly in materials science, 22 at par with well-known generic allatom force fields such as CHARMM 84 and Amber. 85 The expanded set of bead types of Martini 3 already has enabled a further extension of the model into the areas of, for example, green solvents (ionic liquids, 86,87 deep eutectic solvents 88 ), coacervates, 89 and intrinsically disordered proteins (IDPs). 90

| Associated tools
In connection to the actual Martini force field, we, and many others, have developed a wide range of tools to facilitate setting up, running, and analyzing MD simulations with this model (Figure 1). Examples include the widely used tools insane 91 and CHARMM-GUI Martini-maker 92,93 to build multicomponent lipid-based systems, tools to facilitate highthroughput sampling of protein-protein interactions such as DAFT 94 or PANEL, 95 the MERMAID tool 96,97 to help setting up and analyzing membrane proteins, the martinize tool to build Martini topologies from biomolecular structural data, 50 polyply, a python script to build complex starting structures involving arbitrary (bio)polymers, 98 and the Nanodisc Builder, to construct starting structures for protein-lipid nanodiscs. 99 Another category of auxiliary tools consists of automatic topology builders, ideally converting a simple smiles string into a complete Martini topology file. Some of these focus on generating bonded potentials, such as the pyCGtool 100 and Swarm-CG, 101 whereas others target the challenge of mapping an underlying chemical structure to its CG Martini representation like autoMartini 102 and the graph-based cg_param algorithm. 103 Recently, the multi-objective based Swarm-CG algorithm has also been used as a tool to improve the nonbonded interactions of Martini lipids. 104 Dedicated analysis tools include Prolint 105 and PyLipID 106 for identifying and analyzing protein-lipid interactions, MemSurfer 107 to define and analyze membrane surfaces, MDVoxelSegmentation 108 for automated lipid leaflet detection, cgHeliParm 109 to analyze CG DNA, a tool based on machine learning to identify membrane domains called NMFk, 110 a tool to directly calculate SAXS spectra from Martini simulations, 111 a tool to analyze the aggregation state of self-assembling systems (MorphoScanner 112 ), and a tool to analyze diffusion on curved surfaces (CurD 113 ).
Furthermore, many enhanced sampling strategies have been optimized for Martini, including metadynamics, 114,115 temperature replica-exchange, 116 replica-exchange umbrella sampling, 117 ensemble-based CGMD, 118 Monte-Carlo enhanced MD, 119 iMapD, 120 AI-driven transition path sampling, 121 and Markov state models (MSM). 122,123 Also worth mentioning is the recent use of Martini in guiding genetic algorithms for inverse design problems (evo-MD). 124

| Multiscale options
The Martini force field has been interfaced with higher or lower resolution models in a variety of multiscale setups ( Figure 2). In sequential multiscaling, CG systems are back mapped to all-atom resolution, often used as a validation test assessing whether or not certain CG configurations remain stable when simulated with a higher resolution force field, or, to obtain more detailed local structural information. Sometimes CG simulations with Martini are only used as an intermediate step during system relaxation, which can be more difficult when directly performed at the all-atom level (e.g., to equilibrate the packing of lipids around a membrane protein). Also, in these cases, a backmapping step is required after the initial equilibration phase. Backmapping can be done in different ways, either using a geometric approach (e.g. with the backward tool 125 ), or a fragment-based approach (using CG2AT 126,127 ), or more recently with a machine learning-based approach (GLIMPS 128 ). At the other end, Martini has been recently connected to mesoscale membrane models with the TS2CG tool, allowing backmapping from triangulated surfaces to the CG level. 129,130 In concurrent multiscaling approaches, the system is partly described at the all-atom level and partially with Martini. In a typical setup, the region of interest is described in atomistic detail, and the surrounding solvent at the CG level, potentially combining accuracy with speed although often resulting in being neither accurate nor fast. The challenge lies in how the two force fields are combined. This can be achieved with the use of virtual sites, [131][132][133] or making use of explicit cross-parameterization as is done in the PACE force field [134][135][136][137][138] or in the work of the Jorgensen group, mixing all atom solutes and polarizable Martini water. 139 In principle such a concurrent multiscale setup can be extended to more resolution levels, as was nicely demonstrated in a triple layer QM/MM/Martini simulation. 140 Coupling to a continuum solvent has also been realized. 141 A third multiscale approach is the use of adaptive resolution schemes, in which the resolution of compounds can change on the fly. Here the main challenge is to define a transition zone that facilitates the change in resolution in a thermodynamically consistent way. Pioneering efforts in this direction, featuring Martini as the CG force field, have been reported for a number of model systems. 142,143 Using such an adaptive resolution scheme, Martini can also be efficiently interfaced with continuum solvent models allowing volume changes on the fly. 144 A connected multiscale approach is that of resolution exchange, a version of Hamiltonian replica exchange coupling parallel simulations F I G U R E 2 Selected Martini multiscale options. (a) Transformation from triangulated surfaces with TS2CG, adapted with permission from Pezeshkian et al. 130 Copyright 2020 Springer Nature. (b) Backmapping from supra CG and to all-atom with insane, adapted with permission from Wassenaar et al. 125 Copyright 2014 ACS. (c) Adaptive resolution setup for water using SWINGER, adapted with permission from Zavadlav et al. 143 Copyright 2019 RSC. (d) Adaptive resolution setup for an all-atom protein in a CG environment, reprinted with permission from Zavadlav et al. 142 Copyright 2014 AIP. (e) Concurrent CG/AA membrane protein systems using the PACE crossinteractions, taken from Qi et al. 138 (f) Concurrent CG/AA lipid vesicle using virtual sites, taken from Liu et al. 133 (g) Resolution exchange protocol for coupling AA and CG resolution, adapted from Liu et al. 145 running at different resolutions. Until now, we are aware of only few attempts to use the Martini model in such a setup. [145][146][147]

| Open science policy
Following the philosophy of the late Prof. Berendsen to keep things pragmatic and share developments among colleagues, 148 we maintain an active forum and organize regular workshops and tutorials. [149][150][151][152][153][154] The policy of sharing of data, often preliminary models (beta-versions), has helped in popularizing Martini, and provided valuable feedback from the community to the main developers. The Martini web portal (http://cgmartini.nl/) 149 plays a central role in this respect, with currently >6000 registered users. A Martini twitter account (@CG_Martini) was established in 2020 to further expand the communication options. Currently, we are in the process of establishing a Martini database (MAD, https://mad.ibcp.fr/), 155 which will serve the need for a well-documented repository of both validated and developmental Martini topologies across the molecular sciences.
Due to its popularity, Martini has been implemented not only in its traditional host Gromacs, 156 but also in other software packages including NAMD, 157 LAMMPS, 158 ddcMD, 159 openMM, 160 CEMD, 161,162 Materials Studio, 163 and Genesis. 164 In addition, the Martini force field is integrated into the HADDOCK workflow. [165][166][167] Note, however, that for the implementations outside Gromacs it is not always clear which version of Martini is available, and the full range of options might not be supported.

| MARTINI 3
In this section, we review the novel aspects of Martini 3 compared to the previous versions, including strategic choices made for the parametrization of new beads and improvements of the interaction matrix, improved treatment of proteins, and guidelines for the best practices for Martini 3 simulations.

| Strategic choices for the bead parametrization
Given the demand for beads that provide broader and more accurate coverage of chemical space 168 and a clear understanding of the stickiness problem, 82 Martini 3 83 adopted a different strategy for the parametrization of its beads compared to previous versions. Overall, the iterative approach did not use a rigorous separation of calibration and validation tests (Figure 3a). In addition, more experimental targets and a large set of simulation systems of different complexity were used for the parametrization of the Lennard-Jones pair interactions. This strategy was fundamentally different from Martini 2, which mainly used liquid-liquid partitioning and simple systems represented by isolated beads for the main calibration step. 35 As Martini is based on pair interactions, additional target systems, and properties helped version 3 to cover more points in the extensive interaction matrix and to prevent possible overestimation of interaction strength involving systems with two or more connected beads.
The first step of this iterative process was to delineate the "Martini universe," which did not only involve defining the number of generic chemical types and sizes of the model, but also a better definition of mapping choices, interaction levels, and bead assignments. For instance, small and tiny beads were not only considered suited to ring-like compounds as in Martini 2, but were redefined here to be also useful for 3-to-1 and 2-to-1 mapping of linear and branched chemical groups. 83 In addition, both self-and cross-interactions of different bead sizes were included as parameters to be calibrated Another novelty was the introduction of the "size-shape concept," which improves the volume and shape of coarse-grained models of molecules compared to their atomistic references ( Figure 3b). 83,169 This approach requires the definition of consistent mapping rules for the selection of bead sizes and for the calculation of bond lengths between beads. Together with strategic parametrization choices, the new Martini 3 model presents a fully revised interaction matrix and new intermediate interaction levels, added to smoothen the transition between different chemical types. Additional improvements came from the separation of pair interactions into organic, ions, organic-ion, and water blocks, which were optimized independently. In particular, this approach allowed significant improvements for water, which in Martini 3 has its own bead type (W) and specific interaction levels. 83 With a defined Martini universe, the parametrization strategy was divided into "tiers", which represented systems and target properties with different level of complexity. Tier 0 contained simple molecules (with 1 or more connected beads) and many thermodynamic reference data that were used mainly for calibration of Lennard-Jones parameters. In addition to the liquid-liquid partitioning, miscibility was used as one of the main parametrization targets, allowing to probe at the same time self-and cross-interaction between the two molecular species. In addition, the balance of different bead sizes in this tier was intended to be multiresolution, with ideal mixing behavior 170 of pure solvents composed of molecules mapped at different bead sizes. Excess free energies of mixing and trends in liquid-liquid interfacial tension, densities, and solvation-free energies were also used in this tier. Different classes of molecules used complementary targets for the parametrization. For instance, calibration of charged systems also relied on ion transfer properties between ionic liquids and aqueous salt solution, which can be considered analogs to liquid-liquid partitioning used for neutral molecules. Trends according to Hofmeister series 171,172 were also considered ( Figure 3c). Additionally, specific systems such as metal ions, charged amino acids, and simple sugars used osmotic pressures as a reference, which is a target previously used for ad hoc rescaling corrections applied in Martini 2. 60,78,81 In the intermediate tier 1, tests formulated as "yes-or-no" questions were used to evaluate the quality of the model, with parameters not covered in tier 0 also tested and fine-tuned here. Properties of more complex systems such as soluble proteins, sugars, and bilayers were tested to detect possible emerging artifacts at a larger structural scale than in tier 0. Examples of "yes-or-no" questions were simple self-assembly MD simulations, indicating if certain lipids could or not form bilayers or if soluble proteins and sugars mainly stay in the monomeric form at low concentrations instead of forming a unique phase-separated aggregate. Additional systems including peripheral membrane proteins, protein nanopores, transmembrane proteins, ionic liquids, and other soft materials were included in these quality control tests as well. 83 In the final tier 2, quantitative tests involving complex systems were performed, including comparisons with experimental or atomistic simulation data. Here most of the systems were considered for validation, and many comparisons involving dimerization free energies were used, including systems such as aedamers (a biomimetic molecular system useful to study folding and self-assembly), nucleobases, and transmembrane helices. At every point that one of these tests failed, initial choices were reconsidered in an iterative process, including the definition of the initial Martini universe of beads.

| Improved interactions and better coverage of the chemical space
As a result of the parametrization, Martini 3 showed many improvements, in particular for systems which heavily rely on smaller beads, such as aromatic and aliphatic rings or branched fragments. 83,169 Liquid-liquid partitioning involving water/hexadecane, water/chloroform, and water/octanol biphasic systems showed a mean absolute error compared to experimental data of around 2-3 kJ/mol ( Figure 4a). 83,169 At the same time, liquid-liquid mixing behavior was captured accurately, from common solvents 83,169 (Figure 4b) to complex mixtures involving small conjugated molecules, 173 coacervates, 89 ionic liquids, 86 and deep eutectic solvents. 88 As a direct consequence of the improved calibration of the F I G U R E 4 Examples of Martini 3 with improved interactions and use of labels. (a) Water-oil transfer free energies (ΔG) were computed for around 260 data points using Martini 3. The plot includes data used for calibration and validation of the Martini 3 beads. (b) Comparison between Martini and experimental free energy of mixing (ΔG mix ) and excess free energy of mixing (ΔG ex ) for the benzenecyclohexane mixture as a function of the mixture composition. Â1 is the molar fraction of benzene. The analogous Martini 2 system forms a biphasic system at the same condition. (c) The dimerization free energy surface of chlorobenzene is plotted on the 2D coordinate space formed by the distance between the center of geometries of the two molecules and cosθ, where θ is the angle between the two vectors perpendicular to the planes of the two molecules. The profiles are obtained with GROMOS (left panel) and Martini 3 (right panel). (d) Selfassembly of aedamers. The left panel shows the dimerization free energies (ΔG dim ) of pegylated monomers of DAN and NDI. The right panel shows the self-assembled duplex dimer formed by amide-linked tetramers of NDI (green) and DAN (orange). (a) and (d) Adapted with permission from Souza et al. 83 Copyright 2021, Springer Nature. (b) and (c) Adapted with permission from Alessandri et al. 169 cross-interactions between the different bead sizes, large dimerization dissociation barrier artifacts between cyclic molecules previously observed in Martini 2 no longer exist. 83,169 Free energy surfaces of dimerization for many aromatic compounds in water were improved, with Martini 3 showing better agreement than the Martini 2 version in terms of depth and position of minima using atomistic force fields as reference ( Figure 4c). 169 Similar improvements in nucleobase stacking and packing of poly(3-hexyl-thiophene) in bulk heterojunction solar cells were also observed. 83 However, inaccuracy in packing and binding modes may still be in place, given the directionality limitations of the bead interactions. For instance, similarly to H-bonds, T-shaped stacking interactions are also not captured by standard Martini models, including version 3, due to the lack of specific quadripolar electrostatic interactions. 169 In combination with the improvements in the interactions of all bead sizes, Martini 3 showed a substantial increase in the number of chemical bead types, which are separated into seven classes: polar (P), intermediately polar (N), apolar (C), halo-compounds (X), monovalent ions (Q), divalent ions (D), and water (W). 83 W and D are classes composed of a special and dedicated bead type, while the other classes show more diversity of options. The main classes of the organic blocks-P, N, and C-are composed of six bead types with different degree of polarity, enabling a more precise definition of different chemical groups compared to the previous version. In particular, the N class was greatly expanded in relation to Martini 2, allowing differentiation of different chemical groups in the most populated region of the chemical space-as observed by Kanekal and Bereau 168 -with compounds that have octanol-water partitioning free energies between À5 and +5 kJ/mol. In addition to the expansion of the number of levels representing N beads, Martini 3 introduced a new X class dedicated to model halo-compounds. This approach was necessary because the self-interactions of halo-compounds follow the opposite trends with increasing polarity when compared with typical apolar compounds (represented by the C and part of the N bead types). In analogy with the N class, Q beads were also expanded, allowing the separation of soft hydrophobic-like monovalent ions (e.g., tetrafluoroborate, BF 4 À ) from the hard hydrophilic-like ones (e.g., formate, COO À ). Among new applications of the charged beads, the new Martini 3 models of phosphoinositides stand out, as they were able to reproduce experimentally known protein-binding poses as well as phosphoinositide aggregation tendencies in the presence of calcium ions, the latter modeled with D-beads. 174 In addition to the new bead definitions, Martini 3 enhanced the ability to modify bead properties depending on the chemical details of the underlying moieties. These modifications are controlled via labels, added as subscripts to the bead chemical types. This concept already existed in the previous versions of Martini 29,35 but was only used to account for different hydrogen bonding capabilities of N and Q beads. In Martini 3, the label concept was generalized and can be applied to all beads. Labels slightly change the properties of the beads they are attached to, causing small changes in the interaction matrix in a bead-size dependent way. 83 Three of these labels (hydrogen bonding, electron polarizability, and positive/negative) are chemically specific, which means they are only applied to certain chemical types. 83 Hydrogen bonding labels (applied to P and N classes) are sufficient to reproduce trends in nucleobase pairing selectivity, while electronic polarizability labels (used in C and X classes) allow for specific and strong stacking between electron-donor and electron-acceptor aromatic rings in aedamers ( Figure 4d). The other two labels (self-interaction, partial charge) of Martini 3 are more generic, and can be applied on top of all beads of the organic block, including beads that have already been modified with some of the chemically specific labels. Partial charge labels, which are useful for large charged molecules represented by more than one bead, have been used for the modeling of ionic liquids 86,87 and deep eutectic solvents. 88 Together, all combinations of sizes, chemical types, and labels give a total of 843 bead types, which correspond to about 15 times more options than in the previous version of Martini. A recent application of the new bead types by Yu and Wilson 175 to model self-assembling dyes showed their promise: the Martini 3 models outperformed CG models based on a force-matching approach.

| The differences and improvements of protein models in Martini 3
Currently, proteins are possibly the system with the most obvious improvements in Martini 3. One of the key differences in the Martini 3 protein model compared to the previous version is the total reparametrization of the side-chains. In Martini 2.2, the side chains are considered "overmapped", with on average too few atoms incorporated in a Martini bead. In Martini 3, the new mapping avoids the use of regular size beads for fragments containing less than four nonhydrogen atoms. The packing of aromatic residues also greatly improved with the use of "tiny" beads, the smallest Martini 3 bead size. Although most of the remaining bonded parameters were not changed compared to Martini 2, additional improvements were obtained by side-chain corrections, 176 which are dihedral terms involving side-chain and backbone beads based on experimental reference structures. A final improvement was a new definition of backbone bead types, which no longer depend on the secondary structure. They are now represented by P2 beads, except for proline (SP1a), alanine (SP2, with an additional bead for the side chain), and glycine (SP1). These are the beads which best represent the liquid-liquid partitioning of the secondary amides in the protein backbone of these residues. As an alternative to ElneDyn, 51 Martini 3 proteins can also be combined with G o-like models, [177][178][179] which can greatly improve the representation of protein flexibility. In fact, the recent implementation of G oMartini via virtual particles marked a milestone, showing the use of Martini 3 to study long-range structural communication and allosteric changes caused by single-point mutations. 180 After the release of the open-beta version, Martini 3 has been evaluated in many protein systems and properties, including peripheral membrane proteins, protein flexibility, protein-small molecule and protein-protein interactions. The group of Vanni demonstrated that the model is mostly able to accurately characterize the membranebinding behavior of peripheral proteins, identifying key residues found to affect membrane binding. 181,182 As a result of the improvements in packing and interactions in Martini 3, binding of small molecules to protein pockets was successfully reproduced for a variety of systems such as GPCRs, nuclear receptors, and kinases. 183 More of these results are discussed in Section 4.4 of this review. Using parallel tempering metadynamics simulations, the group of Cournia 181 indicated a possible underestimation of the dimerization binding strength for the open-beta release of Martini 3 proteins in solution, while it performs well in describing the association of transmembrane (TM) proteins. Spurious supramolecular protein aggregation, observed in Martini 2.2P, no longer occurs, with near-native dimer complexes identified as minima in the free energy surfaces. 181 As of the time of writing this review, only TM peptides and IDPs have been evaluated with the final release of Martini 3. The dimerization of helical single-pass TM peptides showed high accuracy in predicting binding affinities 83 and dimerization binding modes. 184 The Lindorff-Larsen group tested Martini 3 for IDPs and found that Martini 3 underestimates the global dimensions of IDPs when compared with experimental small-angle x-ray data. 90 IDPs were not considered as part of the calibration and validation tests in Martini 3, explaining this discrepancy. A possible remedy, that would not affect the overall balance of interactions, is to increase the hydration strength of the peptide backbone in disordered regions using virtual sites.

| Best practices and tips
All the improvements incorporated in Martini 3 are only useful if the force field is properly used by the modeling community. Therefore, in this subsection, we summarize key advice and best practices to parametrize new molecules and run simulations with the new version of Martini. Overall, the molecular dynamics settings for Martini 3 simulations follow the "new" Martini set of parameters proposed by De Jong et al. 73 The default time step is 20 fs. A reaction field should be used as standard option for electrostatic interactions. However, the model was also systematically tested with Particle Mesh Ewald (PME). The use of PME showed some improvements for highly charged systems, 86,88 although with significantly reduced computational performance, in particular, due to the reduced computational scaling over processors. With the ongoing development of a polarizable water model compatible with Martini 3, we expect further benefits in modeling accurate electrostatic interactions. Special attention should also be paid to CG models with highly coupled constraints. As an example, these are present in the current cholesterol topology, which contains two obtuse triangles sharing an edge. 48 This configuration can lead to nonconverged constraints, which cause artificial temperature gradients. 185 For such systems, a more conservative set of LINCS settings is recommended. 185 Lastly, the regular size water bead is the standard option to model liquid water in Martini 3. Although small and tiny water beads were calibrated to have near-to-ideal miscibility between all water bead sizes and also to reproduce the hydration free energies of neutral beads solvated in regular water beads, 83 their use and transferability in different systems should be tested and verified by the user with caution. Small and tiny W-bead sizes are currently recommended only for certain applications, such as the modeling of water in protein pockets or channels where regular water beads would not fit. 183 A small fraction of small or tiny water beads could be added to the solution to fill these pockets. Such an approach can potentially improve the small molecule binding to the pockets of certain proteins.
The parametrization of new Martini models has a well-defined set of rules in version 3, which should be followed to build high-quality models and maintain internal consistency with other Martini 3 molecules. A quick guide for unexperienced users is described as follows.
(1) Mapping of the atomistic model to the CG model: Split the molecule in reasonable chemical groups (the building blocks) in accordance with the standard Martini bead sizes: regular (4-atoms-to-1-bead), small (3-atoms-to-1-bead), and tiny (2-atoms-to-1-bead). A higher resolution (smaller bead size) allows for a better reproduction of symmetry and chemical details, while larger bead sizes result in simpler models which are more computationally efficient. In terms of volume and shape criteria, tiny beads are the preferred option for aromatic rings, while small beads are preferred for aliphatic rings. If a molecule contains highly branched groups, halo compounds, and/or charged fragments, consult section C1 of the supplementary information of the Martini 3 paper. 83 (2) Assignment of chemical bead types: According to the chemical group underlying the bead, the user should select a fragment according to the "Martini Bible" (supplementary tables 24 and 25 of the Martini 3 paper 83 or table 1 of the Martini 3 Small Molecule paper 169 ). If no matching chemical group is available in these resources, a suitable bead type can be selected using the liquid-liquid partitioning data of the isolated fragments (see supplementary tables 17, 18, and 19 of the Martini 3 paper 83 ). Sections A4 and C2 of the supplementary information of Souza et al. 83 provide a detailed description of the use of labels.
(3) Bonded parameters, volume, and shape: Bonded parameters can be obtained from atomistic simulation data. Bond distances should be based on the center of geometry of the atoms that are incorporated in a bead, including the underlying hydrogen atoms. This approach showed a good match with molecular volume, surface accessible areas (SASA), and bulk densities. 83,169 Additional fine-tuning can be performed by scaling selected bond lengths as needed to properly match Connolly surfaces and SASA. We recommend the use of Van der Waals radius from the work of Rowald et al. 186 as a reference for the atomistic calculations. In the case of the Martini models, the radius of the CG beads should be used. Appropriate force constants and additional bonded terms such as angles or dihedrals should be calibrated with respect to structural data obtained from either reference atomistic simulations or a structural databank. Note that the bimodality of the distributions of bonds, angles, and asymmetric dihedrals may not be well-captured by Martini given limitations of potentials implemented in Gromacs (or other MD codes), although cosine-based functions such as the Ryckaert-Bellemans potential could be used in certain cases. Section C3 of the supplementary information of the Martini 3 paper provides more details. 83 (4) Improving numerical stability and performance: Two main factors affect the numerical stability of Martini simulations: (i) the number of bond constraints and the use of virtual particles to manage this number; and (ii) the definition of proper dihedrals and collinear 3-body angles. Virtual sites can be conveniently used to model extended stiff molecular structures, which reduces the number of constraints that these systems would need. This approach not only increases the numerical stability, reducing the risk of nonconverged constraints, 185 but it also improves the computational performance. 169 Traditional proper dihedrals used in atomistic MD simulations may also be a source of instability in CG simulations. This is usually the situation if the three-body angles i-j-k and j-k-l in a dihedral i-j-k-l are near to 0 or 180 , as described in the work of Bulacu et al. 187 One possible solution involves the use of restricted bending potentials (function type 10 in Gromacs 188 ). Alternatively, virtual "dummy" sites can be used to redefine the dihedral in a way that collinear three-body angles are avoided. 187 The Martini 3 Small Molecule paper 169 provides additional advanced design strategies involving virtual sites. Additional discussion involving numerical stability issues is also available in section C4 of the supplementary information of the Martini 3 paper. 83 (5) Generating protein models in Martini 3: The accuracy of Martini 3 proteins depends on the reference atomistic structure chosen by the user to generate the model. This choice has an impact on the assumption of the secondary structure of the protein, which directly determines the backbone bonded parameters, as in the previous Martini 2 version. 49,50 Additionally, the reference structure is also used to define the dihedral side-chain corrections 176 and provides a bias toward certain tertiary structures, which is mainly controlled by the elastic network 51 or the G o model. 177 Therefore, the users should be careful with choosing the most appropriate reference structure; detailed knowledge of your protein structure and dynamics is recommended. The possibility of issues in the reference, such as structural crystallographic artifacts or missing regions that need to be computationally modeled, should always be double-checked. As in atomistic simulations, protonation states, the presence of cofactors, or posttranslation modifications cannot be ignored, and maybe need to be manually added. Given a high-quality atomistic reference model, the user still should consider fine-tuning the elastic network and/or G o model, to improve the representation of protein dynamics. For elastic networks combined with Martini 3, initial choices for force constant should be 700 kJ/mol nm 2 , with a cutoff between 0.8 and 0.9 nm. The choice for a higher force constant than Martini 2 is related to the differences of the nonbonded interactions of the models and also to avoid possible pitfalls with low force constants. 82 For G o models, the initial choice of the depth of the potential is 10-12 kJ/mol. Both elastic and G o models can also have their contact map manually refined in MAD 155 to improve certain aspects of the flexibility. In the current implementation, Martini 3 protein models should not be used in applications that involve large conformations changes as occurring during protein folding.

| EXAMPLE APPLICATIONS
In this section, we will highlight the main application areas of Martini, divided into the categories of biological systems (Section 4.1), the bio-nano interface (Section 4.2), and materials science (Section 4.3). We also include dedicated sections illustrating three of the most important reasons to use a CG model such as Martini, namely in high-throughput workflows (Section 4.4), in systems of increasing complexity (Section 4.5), and to bridge detailed models toward the whole cell level (Section 4.6). In each section, we provide key examples of recent work. By no means, we aim to be comprehensive, which would be nearly impossible given the large number of Martini-based publications (Figure 1b). In Table 1, we provide the current "records" in terms of system size, complexity, and time-scale that have been achieved with Martini to date.

| Biological systems
Martini has historically primarily been used for biological systems, driven by its initial development as a lipid model followed by protein parameters. The lower level of detail and its correspondingly lower computational cost has enabled a range of applications that are still difficult to reach atomistically. It is remarkable how many important biochemical and biophysical problems can be addressed when simulations can routinely reach time scales of tens to hundreds of microseconds and length scales of tens to hundreds of nanometers. Below we highlight a number of biological applications to illustrate the breadth of biophysical and biochemical problems that can be addressed using Martini.
Many early applications of Martini were on the properties of lipid systems. As the time-scale of mixing even in simple binary mixtures is microseconds and goes up dramatically in the presence of domains, atomistic simulations until recently were limited in modeling many aspects of lipid biophysics. With Martini, processes such as domain formation, membrane fusion, lipid phase transitions, and monolayer breakdown came into reach. Comprehensive reviews of membrane simulations, including many of these pioneering Martini-based papers, were published in 2019. 20,193 A more recent development is the modeling of more realistic and complex membranes. A landmark study by Ingolfsson et al. in 2014 described an asymmetric membrane with over 60 lipid types on a time scale of tens of microseconds, 194 still one of the most detailed models of a generic plasma membrane. Subsequent models have targeted specific cells, including specific neuronal, 195 skin, 196 and epithelial cells, 197 as well as bacterial membranes. 198,199 With the growth of lipidomics data, 200 such simulations continue to be of major interest in the studies of membrane structure and increasingly biological function when combined with proteins. Martini allows also to model the properties of lipid mixtures across an entire phase diagram. 201,202 This is both of basic interest and a possible route to improving lipid parameters, although it continues to be a challenge to efficiently calculate phase diagrams for models as detailed as Martini, let alone more detailed atomistic models. A recent atomistic study of just a few phase points highlights the sampling challenges. 203 Note that, as an inherent drawback of CGing, the entropy/enthalpy balance is distorted, implying that one has to be careful in interpreting temperature-dependent phenomena. 66 Note: The interpretation of the timescale in Martini is not well-defined. Due to the missing friction from the atomistic degrees of freedom, the potential energy landscape is smoothened, resulting in an effective speedup. The amount of speedup is, however, system dependent and has been estimated as ranging from a factor of 2-10. 66 So, a 1 microsecond Martini simulation would correspond to a 2-10 microsecond simulation with an all-atom model.
A second major area of applications in biological systems consists of membrane proteins, in particular addressing lipid-protein interactions at different scales. The time scales involved in sampling many types of lipid-protein interactions are tens to hundreds of microseconds, 192,204 readily achievable for Martini simulations. The Martini 2.2 lipidome and lipid-protein interactions have been shown by a large body of work to be quite accurate, as reviewed comprehensively by Corradi et al. 21 These applications also typically do not require an accurate description of internal protein dynamics, and thus make use of elastic networks to keep the protein close to the reference state (usually a crystal structure). In many cases, simulations have identified specific binding sites for key lipids. An early example is the binding of the cardiolipin to respiratory proteins, 205,206 where Martini simulations reproduced known binding sites and predicted additional ones that have been verified experimentally. Other examples include the binding of PIP lipids and cholesterol to Kir potassium channels, [207][208][209] which are regulated by these lipids, the binding of ceramides to VDAC, 210 and cholesterol binding to GPCRs 211-215 and various transporters. 216 Martini simulations also lend themselves to identifying specific lipids based on low-resolution cryo-EM or crystallography densities of membrane protein structures. 217,218 As the number of cryo-EM structures continues to explode, this is likely to become increasingly useful to identify specific roles for lipids in membrane protein structure and function.
In addition to specific binding sites, simulations have been used to characterize the local lipid environment around membrane proteins, in a long tradition of research on membrane adaptations to proteins; a tradition going back to theoretical models that introduced the mattress model of hydrophobic matching 219 and some of the earliest atomistic membrane protein simulations. 220 Corradi et al. introduced the concept of lipid fingerprints, the unique local environment near a membrane protein induced by lipid-protein and lipid-lipid interactions in a membrane of a complex composition. 192 This has subsequently been investigated for a range of other membrane proteins, including GPCRs, 221 aquaporins, 222 plant light-harvesting complexes, 223,224 mechanosensitive channels, 225 neurotransmitter transporters, 226 and may play an important role in determining membrane organization. The local membrane environment, as exemplified by a growing number of Martini-based simulations, plays a key role in the organization of membrane protein complexes, 137,227-231 binding of peripheral membrane proteins, 232,233 and the equilibrium distribution of membrane protein conformational states. 191,[234][235][236][237] In some cases, Martini simulations can also directly provide a link to the mechanism and function of membrane proteins. Recent examples include the lipid scrambling by TMEM16K 238 and the modulation of the hERG potassium channel by ceramides (Figure 5a, left panel), where a combination of Martini, atomistic simulations, and electrophysiology enabled a detailed functional model. 239 Another example is the mechanism of partitioning between raft and nonraft-like domains in membranes depending on the preferential interaction with specific lipids (PIP2 240 or polyunsaturated chains 241 ) or on palmitoylation. 240,242 Also, pore formation and translocation of peptides across membranes have been explored in the case of peptides or proteins with stable secondary structures, such as small, highly helical antimicrobial, and cell-penetrating peptides 243,244 or even larger protein toxins. 245 Note that processes such as translocation of polar compounds that require transient pore formation suffers from a significant overestimation of the pore-free energy in the Martini model. 246 The scale of Martini furthermore enables large-scale simulations of the role of proteins and lipids in membrane remodeling (Figure 5a, right panel). For instance, to identify possible roles for SNARE-family proteins [247][248][249][250] in fusion of synaptic vesicles (notably, proteins increase lipid tail protrusion and promote stalk formation), also relevant for viral infection by enveloped viruses like Influenza B. 251,252 Another example is lipid droplet formation, a rapidly growing area in the literature. 253 Early work in this area was published by Khandelia, who showed that triolein phase separates and forms oil lenses in phospholipid bilayer membranes. 254 More recently, the group of Monticelli showed that membrane asymmetry imposes the direction of lipid droplet budding, 255 and suggested that the preference of monotopic membrane proteins for the lipid droplet surface (measured experimentally) may be due to their stronger interaction with triglycerides. 256 Also, Prasanna et al. showed that seipin, an oligomeric ER protein, can trap triglycerides in the ER bilayer via its luminal hydrophobic helices, lowering the nucleation barrier for the initial lens formation (Figure 5b) 257 -a prediction confirmed by simulations with the SPICA coarse-grained force field. 258 In addition, Martini simulations made significant contributions to exploring the roles of proteins and lipids in curvature sensing, 259,260 curvature induction, 261 and membrane deformation. 262,263 For example, the group of Hummer simulated membrane remodeling by FAM134B reticulon homology domain, important in regulating the size, and shape of the endoplasmic reticulum, and by nhTMEM16, 264 a scramblase that dissipates the compositional asymmetry between membrane leaflets, therefore, inducing the conversion of vesicle buds into flat membranes. 265 For more work in this area, see Pezeshkian and Marrink. 266 A growing number of papers use Martini for other protein-based applications besides the membrane proteins mentioned above. In principle, the sampling power of Martini enables a large number of problems to be addressed, but in practice limitations of Martini 2.2 made the interpretation of some simulations challenging. One example is the simulation of protein solutions, one of the systems Martini 3 parameterization has targeted. In this case, a protein solution now gives the right salt-dependent behavior without ad-hoc adjustments to protein-protein interactions. 267,268 A more elaborate parameterization of protein-protein interactions enables more reliable simulations of complex formation and the free energy of protein-protein interactions. 181 Martini 2 had a limited ability to model small molecules, but Martini 3 has sufficient detail to be used for screening simulations for small molecules binding to proteins. 183,269 Simulations of intrinsically disordered proteins have been challenging for atomistic simulations, in part because force fields had not been fine-tuned to give the correct free energies of solvation for disordered chains 270 and in part because sampling disordered proteins are limited and free energy methods are not readily applicable. A recent paper shows that Martini 3 is a viable method to study IDPs (Figure 5c, left panel), although there is room for improvements that should fit well within the available Martini 3 framework. 90 This opens opportunities for a number of different problems in fluid-fluid phase separation in biological systems. A recent example is the work of Tsanai et al. on the formation of a coacervate phase by polylysine and polyglutamate in solvent and the partitioning of solutes between solvent and coacervate (Figure 5c, right panel). 89 A final application we would like to mention is the use of Martini in integrated modeling approaches. In this case, the physics inherent in the Martini model is one ingredient in more complex calculations that may include a variety of experimental data at different levels of resolution. A very successful example is HADDOCK, developed by the Bonvin group and co-workers, which can incorporate Martini simulations as element in modeling workflows for both protein 166 and protein-DNA complexes. 165 In this section, we focused on lipids and protein applications, which make up the bulk of the Martini studies in the literature. There have also been interesting recent applications using the nucleic acid parameters, for example, for modeling aspects of drug delivery 196,271 (Figure 5d), nucleosome formation, 272  DNA-based nanopores, 275,276 and carbohydrates, in modeling cellulose 277 and glycosylated proteins, 60,278 and development of important biopolymers including peptidoglycans 279 for more realistic simulations of Gram-negative bacterial membranes. We expect such applications to become more prevalent as Martini 3 continues to be developed.

| The nano-bio interface
Soon after the Martini model was extended to proteins, Martini models for nonbiological materials were developed by different groups, with a particular focus on nano-sized particles. One of the main goals of such models was to characterize the interaction between nonbiological and biological materials, notably model membranes.
Back in 2008, Wong-Ekkabut et al. published the first Martini model of fullerene 64 and characterized the permeation of the nanoparticle through model membranes, as well as its effect on structural, dynamic, and elastic properties of the membrane. Similar models were used around the same time by the group of Sansom to study the interaction of truncated, extremely short carbon nanotubes (CNTs) with lipid membranes, with the aim of exploring the potential use of nano-injectors for the introduction of drugs into cells, 280 and with short peptides, potentially useful for nanotube solubilization. 281 However, it should be noticed that CNTs could not be parameterized based on the experimental partitioning, as normally done within the Martini framework, due to lack of data and major simplifications in the chemistry of the models. The fullerene model, instead, was developed using experimental data on solubility and solidstate properties. 282 The optimized special bead for fullerene was used in subsequent simulations of carbon nanotubes permeation through lipid membranes. 283 The same optimized fullerene model was also used to answer questions on why lipid membranes are more effective than alkanes at dissolving fullerene aggregates, 284 to explore the interaction of fullerene with monolayer membranes, 285,286 (Figure 6a) and to characterize the size dependence of nanoparticle aggregation in lipid membranes. 287 Models of CNTs were also developed more recently by the group of Hummer, who considered the effect of functional groups at CNT rims, and focused on the interaction of CNTs with lipid membranes in the context of solvent filtration devices. 288 The same group also carried out simulations that enabled an interpretation of cryo-TEM images of CNTs in living cells, suggesting CNT-mediated vesicle fusion (Figure 6c). 289 Another important nano-bio interface is the one between graphene (and its derivatives) and model membranes. Pioneering work by the group of Kral showed that individual, small graphene sheets can, in principle, be incorporated into the central region of lipid bilayer membranes. 290 Simulations also provided an insight into a possible mechanism of formation of such sandwiched graphene-bilayer system, starting from graphene-POPC micelles. 290 A more recent related study by the group of Reigada suggested that, as previously observed for CNTs, graphene sheets may also promote vesicle fusion, which hints at a possible mechanism of biological activity. 291 Martini models have also been used to simulate the interaction of graphene and graphene oxide with supported lipid bilayers. 292 As expected, the hydrophobic (graphene) or hydrophilic (graphene oxide) nature of the solid support has a major effect on the organization of bilayer membranes, and CG simulations contributed to the interpretation of experimental results.
Besides carbon nanoparticles, during the past decade much work has been devoted to the interaction of metal nanoparticles (NPs) with biological systems, and particularly gold (Au) NPs. 293 Au NPs have numerous applications in biomedicine, as recently highlighted by the development of Au NPs for the detection of SARS-CoV2 infection in antigen tests. 294 Pioneering work on Martini models of Au NPs has been carried out by Lin and Alexander-Katz. 295, 296 Lin developed and validated a Martini model of Au NPs coated with alkyl thiol ligands that were hydrophobic, or functionalized with cationic (amino) or anionic (carboxylate) groups. 295 The model was used to study chiefly the interaction of the NPs with model lipid membranes as a function of the chemical nature of the coating ligands. 295 Similar and related questions, regarding the interaction of Au NPs with membranes, were also addressed by a number of other groups (Figure 6b). [296][297][298][299][300][301][302] Recent studies tackled more complex questions regarding the NP-membrane interaction. For instance, Rossi et al. used Martini models for Au NPs coated with hydrophobic and anionic ligands to study the effect of NPs on membrane phase separation, 303 the effect of cholesterol on NP passive uptake, 304 and the induction of both positive and negative curvature in liposomes. 305 The effect of phase separation on Au NP-membrane interactions was studied by Van Lehn and coworkers 306 as well as Liang and coworkers. 307 Alexander-Katz considered the synergistic effects of size, surface charge, and ligand chemistry on nanoparticle translocation across membranes, 308 and the effect of applied voltage on the membrane adsorption of cationic nanoparticles was evaluated by Chiarpotti et al. 309 Overall, in these simulation studies, the advantage of the Martini model compared to continuum and all-atom models was clear: continuum models are efficient but lack the level of detail necessary to answer questions related to the effect of specific spatial arrangement or the specific chemical nature of the ligands on a nanoparticle, while all-atom simulations provide the most accurate representation for the system, but often proved incapable in reaching the required timescales.
Considering metal nanoparticles, the interaction with peptides/proteins and small biological macromolecules also raised significant interest. Rossi and coworkers studied the interaction of Au NPs with human serum albumin, 310 the most abundant protein in human blood plasma and in the protein corona formed around Au NPs in physiological environments. The same group also developed a model for Au NPs with more weakly bound ligands, notably citrate. 311 Citrate capping is particularly relevant since it is often used to increase the colloidal stability of gold NPs and as an intermediate toward the synthesis of Au NPs with other functional groups. The colloidal stability of citrate-capped Au NPs was also studied by De Vivo and coworkers using experimental and theoretical methods, including Martini CG simulations. 312 The same group also characterized the behavior of more complex metal NPs, bound to multiple copies of a polycationic cell-penetrating peptide, in contact with a model membrane. 313 Peptide-coated NPs strongly interacted with membranes without affecting their mechanical stability; however, the fixed secondary structure of the peptides in Martini is a significant limitation in predictions of their interaction with lipid membranes.
A number of polymers were also studied, in the framework of the Martini model, to probe the interface between biological and nonbiological systems. Here again, the motivation for the studies has been twofold: On the one hand, polymeric systems are increasingly used in biomedical applications, such as drug delivery or biosensors; on the other hand, the massive use of polymer plastics in consumer products (and their disposal, not always controlled) raises concerns about the possible noxious effects on organisms. [314][315][316][317] Pioneering work on Martini models of dendrimers interacting with lipid membranes was carried out by Lee and Larson. 318,319 The same authors also pioneered the use of the PME algorithm for long-range electrostatics in Martini simulations, which turned out to be important to reproduce the behavior of polyelectrolytes, in particular the formation of pores in lipid membrane by interaction with charged dendrimers. 318,319 Recent examples related to biomedical applications include the works on poloxamers (i.e., Pluronics) and their interaction with lipid membranes [320][321][322] and studies on the complexation of RNA with polymers. 323,324 Another class of copolymers that attracted significant attention for their interaction with biological membranes is styrene-maleic acid (SMA) copolymers. SMA molecules mixed with lipids form disk-shaped nanoparticles known as SMA lipid particles (SMALPs), and are used to extract membrane proteins from lipid membranes without the use of detergents or membrane scaffold proteins. Martini models for SMA copolymers were developed by several groups and used to characterize the mechanism of the formation of nanodiscs. 325,326 Finally, polymer-membrane interactions were investigated by Rossi and co-workers focusing on industrial plastic polymers. It was found that nano-sized polymeric particles can be dissolved easily in phospholipid membranes and alter most of their properties (Figure 6d), 327,328 as later also observed by all-atom simulations. 329 In addition, Martini PEG polymers have been used in a biological context as crowding agent, 330 to induce membrane fusion, 331 and as hybrid lipid/polymer biomimetics. 332 Overall, simulations of nano-sized materials at the interface with biological systems highlight some of the major assets of the Martini model: its versatility, simplicity, and the benefits of an open science policy. Development of Martini models for various materials has often been carried out by groups not related to the groups of initial Martini developers; this has been made possible by the intuitive nature of the model and the availability on the web of all resources and information necessary to carry out the development. The building block approach, in which beads are parameterized individually, makes the model extremely versatile, as it allows to easily combine of existing (and new) particle types to cover unexplored regions of the chemical space, while retaining the ability to mix them with particles originally developed to reproduce the behavior of biological systems. Such a combination of biological and nonbiological building blocks has also facilitated a growing number of Martini models for "hybrid" molecules, such as PEGylated lipids and proteins, 150 peptide amphiphiles, 333 or drug-polymer conjugates, 334 to name but a few.

| Materials science
Soft functional materials are implicated in many innovative applications and their tailored design is expected to lean increasingly on advanced computational modeling. 336 As discussed, the Martini model has matured into a generic force field; the assumptions underlying the model do not preclude its application in (soft) materials science. Indeed, in the last decade, Martini has been increasingly used to model materials such as (block-co)polymers, 44,337 self-assembled supramolecular systems, 70,71,338 organic semiconductors, 339,340 polymer-nanoparticle composites, 341,342 surfaces, [343][344][345] ion-conducting materials, 346,347 and green solvents. [86][87][88] Here, we highlight a few recent and diverse examples spanning the breadth of studies performed with Martini in materials science. For a comprehensive overview of such materials science studies, we refer to our recent perspective. 22 Self-assembled supramolecular materials, featuring a wide range of nanotechnological applications, 348 have been extensively simulated with Martini to elucidate experimental findings or to predict assembly outcomes solely based on the structure of the molecular building blocks. 338,[349][350][351][352][353] Pioneering studies have used short peptides as building blocks to explore biocompatible supramolecular nanostructures. [350][351][352] More recently, supramolecular materials composed of self-assembling synthetic molecules have been increasingly investigated, in particular, by the group of Pavan. 71,349,354,355 A recent highlight is the work of Sarkar and co-workers, who investigated two-component supramolecular polymerization by both experiments and Martini-based MD. 349 Even with two monomeric components, supramolecular copolymerization can lead to several structural outcomes, including supramolecular homopolymers, random copolymers, alternate copolymers, and block copolymers. Comparable interaction energies between the monomer components resulted in supramolecular random copolymers (Figure 7a). In contrast, supramolecular homopolymers could be achieved by leveraging high kinetic stability of the supramolecular structures or low monomer exchange dynamics. Martini simulations, enhanced with well-tempered metadynamics, delivered mechanistic insights into the strength of the interactions between the different monomer components and the monomer exchange dynamics in the supramolecular structures, and provided to be vital for rationalizing the experimental outcomes.
One of the latest materials science branches where Martini has been used is the field of green solvents, with several recent studies modeling ionic liquids 86,87,356,357 and deep eutectic solvents. 88 Green solvents promise to replace more traditional solvents by providing more environmentally friendly solutions. One of the most intriguing areas of application is separation science, given green solvents' tunable, specific, and performant extraction properties. As is common for organic-based materials, their chemical tunability allows for great customization. However, this flexibility also leads to an extremely large number of potential material candidates, giving rise to a chemical space impossible to explore experimentally. Hence, efficient computational modeling of green solvents is vital to enable the targeted design of green solvents. In a recent work, Vainikka et al. developed Martini 3 models for deep eutectic solvents (DESs). 88 The authors selected some well-known and actively studied DES systems and investigated the ability of the models to capture the properties of these green solvents. In particular, they found the models able to capture known liquid-liquid extraction processes-such as, for example, the extraction of 2-methylthiophene and benzothiophene from octane by the DES tetrabutylammonium chlorideÀacetic acid (Figure 7b)-as well as the morphological changes of the shape of surfactant aggregates as a function of DES concentration. The study opens possibilities for in silico design of DESs aided by Martini simulations.
A growing library of Martini polymer models is available, 22 including simple linear polymers, 44,337 branched and hyperbranched polymers, 358,359 conjugated polymers, 339,340,360 block copolymers, [361][362][363] and polyelectrolytes. 364 Several of these models are included in the library of Polyply. 98 Since the pioneering development of Martini models for polyethylene oxide (PEO) 63 and polystyrene (PS), 337 more recent works in this area have explored different aspects and applications of polymers, such as thoroughly revising the transferability of developed polymer models, 44 combining polymers with nanoparticles, 341,342 or simulating complex multi-component systems (e.g., Figure 7c). 98,365 As an example of complex multi-component systems, Li et al. recently developed a hybrid synthetic structure that mimics mechanical actuation, with Martini simulations providing insights into the microscopic processes taking place during actuation. 365 The hybrid material consists of peptide amphiphile supramolecular fibers covalently linked to polymers which are crosslinked across fibers to form a network. The polymers contain a spiropyran-based switch that can be activated with light: depending on light exposure, either a charged hydrophilic form or a noncharged hydrophobic form of the switch can be realized. Switching to the hydrophobic forms leads to water expulsion and an associated shrinkage of the material. By constructing a minimal Martini model of the hybrid material using two crosslinked fibers, Li et al. found that Martini could reproduce the amount of shrinkage in a quantitative way, while providing a microscopic view of the actuation process. 365 The development of robust relationships connecting molecular structure to material nanomorphology remains a key challenge in organic electronics. 366,367 To this end, Martini-level simulations have started to play an important role in elucidating local molecular packing and morphology features as a function of processing conditions and molecular structure. 173,342,[367][368][369] Computational modeling of systems for organic electronics necessitates being able to simulate mixtures of polymers, nanoparticles, and small molecules, and the widening library of corresponding Martini models has enabled such studies. 173,339,340,369,370 Moreover, synthetic molecules can be mixed with biomolecules, as in the work of Mehandzhiyski and Zozoulenko where the conjugated polymer PEDOT:PSS-polymer poly(3,4-ethylenedioxythiophene) blended with poly(styrenesulfonate) (PEDOT:PSS)-is mixed with cellulose nanofibers to investigate PEDOT:PSS/cellulose composites, a material with prospective applications as conductive paper. 368 The group of Zozoulenko has made extensive use of Martini to model conducting polymer-based systems. 340,369 Besides the PEDOT:PSS/cellulose example, a recent highlight includes the development of a protocol to simulate cyclic voltammetry experiments. 371 The protocol consists of a stepwise change of the oxidation state of PEDOT-compensated by the introduction of appropriate counterions-to which the system is allowed to respond. As such, the protocol allows to study the water intake, swelling, and exchange of ions in the conducting polymer morphologies during cyclic voltammetry. When applied to PEDOT:tosylate morphologies, the authors observed significant changes in the morphology during the redox process, with the loss of mass during reduction being in agreement with spectroscopic findings and revealing how tosylate ions are replaced by chlorine ions as the counterions for PEDOT (Figure 7d). Finally, an added benefit of Martini that is particularly important in the field of organic electronics is the possibility of directly backmapping 125 structures to atomistic resolution given the high degree of molecular specificity of the Martini model. This possibility allows to carry out quantum chemical calculations on structures derived from large-scale self-assembly processes, hence enabling advanced multiscale studies aimed at connecting structural and electronic properties. 342,372

| High-throughput workflows
The ability to perform high-throughput (HT) simulations is key to the rational design process of new materials, improved enzymes, or novel drugs. 373 In addition, HT simulations are important to generate fundamental understanding of a system's behavior by systematically probing many different state conditions. Furthermore, HT workflows can be used to create databases that are subsequently mined by bio-or chem-informatics tools. The inherent speed of simulations with CG resolution, in combination with the ease of generating CG topologies and starting structures, has prompted the use of Martini in a variety of HT pipelines, some of which are illustrated in Figure 8 and further discussed as follows.
An early example of such a HT pipeline is from the folding@home project of the Pande group, utilizing the power of a large number of home computers to simulate lipid membrane fusion. 122 The generated ensemble, based on 10,000 short MD trajectories, was combined with MSM to provide the underlying energy landscape of the fusion process. A more recent example in this area is from Hub and coworkers. 377 In their study, they addressed the propensity of membranes to fuse from a computational lipidomics perspective. In particular, the CG resolution enabled computation of $200 free energies of stalk formation in membranes with different lipid headgroups, tail lengths, tail unsaturations, and sterol content. The results from this dataset were subsequently used to explain the finding that the inner leaflet of a typical plasma membrane is far more fusogenic than the outer leaflet.
The group of Stansfeld has used HT Martini-based simulations to maintain a database of membrane proteins and associated webserver (http://memprotmd.bioch.ox.ac.uk/). 374,378 To fill the database, membrane protein structures are automatically identified in the PDB on a weekly basis and converted to Martini representation (Figure 8b). This is followed by a 100 ns MD simulation in the presence of lipids and water to allow a lipid bilayer to self-assemble around the protein. The simulation is continued for another 900 ns simulation to study the dynamics of the assembled membrane, and subsequently converted to atomistic resolution and analyzed. Currently, the database contains almost 6000 simulations of intrinsic membrane proteins. The database can be explored in many ways, for example, to see how a specific protein is embedded in the membrane, or run statistics on a subset of proteins to explore protein-lipid contacts, for instance. In a recent application, cardiolipin binding sites of a large set of 42 different Escherichia coli inner membrane proteins were analyzed and compared. 379 More than 700 specific binding sites were detected, leading to the identification of the molecular basis of a prototypical cardiolipin binding site. A related database that has just been released, SuPepMem, 380 aims to collect simulations of host defense peptides (HDPs) interacting with lipid membranes, at different resolutions including Martini. A first set of 403 systems resulting from the combination of 52 HDPs and six different membranes models, simulated using the Martini 2.2 force field for 5 μs each, is currently available. Further examples of HT work in this area are studies of the Tieleman group comparing lipid fingerprints of 10 generic plasma membrane proteins 192 as well as of 20 different GPCRs, 221 and work from the Stansfeld group reporting affinities for cholesterol binding to a range of membrane proteins based on the computation of binding saturation curves, designed to mimic experimental protocols. 381 Another example of a HT pipeline is the MuMMI workflow (Figure 8d). This pipeline features HT Martini simulations on RAS-RAF protein complexes interacting with a multi-component plasma membrane model, that are coupled to a mesoscale simulation running in parallel. 191,376 In total, 120,000 Martini simulations were performed on patches consisting of about 3000 lipids and containing one or more proteins, with starting configurations extracted from the mesoscale simulation based on a machine learning algorithm. Each patch was simulated for between 1 and 4 ms, with an aggregated simulation time of almost 0.2 s, a record in terms of aggregated simulation time (Table 1). An extended three-scale version of MuMMI has also been developed, coupling between a continuum model, CG Martini, and AA CHARMM to resolve RAS-RAF protein complex membrane dynamics (featuring over 34,000 CG simulations of nearly 0.1 s aggregated simulation time). 382 These huge efforts required dedicated time on next-generation exascale supercomputers, and while some of the final results are still pending, the establishment of the pipeline by itself is already a major feat. A smaller-scale example of HT simulations involving membrane protein aggregation is provided by the group of Sansom, who performed more than 60 independent simulations of GPCR clustering, combined with MSM to reveal cluster populations and dynamics of oligomerization. 228 The Martini can correctly identify known binding pockets. 183 Moreover, the amount of sampling is so large that accurate statistics on binding free energies and entry and exit pathways are obtained at the same time. A potential pipeline for HT ligand screening is illustrated in Figure 8a, 269 combining the benefits of fast methods such as docking and including protein and solvent dynamics. In a related study, Melo et al. used a high-throughput protocol to explore the effect of bilayer environment on the gating kinetics of the membrane channel MscL. Based on over 500 independent CG MD simulations totaling more than 10 ms of effective sampling time, the study revealed that short-chain alcohols have an impeding effect on the opening threshold of MscL, which was corroborated by an experimental efflux assay. 383 The potential for HT drug lead optimization with Martini has also been explored by the group of Bereau. They make clever use of the fact that CG beads in Martini can actually map to multiple chemical fragments. 168 In a landmark paper, 375 they performed HT simulations to derive a membrane-drug permeability surface in terms of two simple molecular descriptors, namely the bulk partitioning free energy and pKa (Figure 8c). The surface is constructed by exhaustively simulating all CG Martini compounds that are representative of small organic molecules. Connecting back to the atomic resolution, their predictions extrapolate to more than 500,000 compounds, covering a large part of the chemical space and providing a clear connection between specific chemical groups and the resulting permeability coefficient. In a follow-up study, the generated data set was interrogated with an artificial intelligence technique to identify the driving forces governing passive drug transport across membranes. 384 A recent study by Melcr et al. also uses HT lipidomics-based computational screening to unravel the role of membrane properties on permeability. 385 HT screening with the Martini force field is also increasingly used in the field of materials science. A prime example is the simulation by Frederix et al. of all possible 400 dipeptides and 8000 tripeptides in search for new peptide-based hydrogels. 338,350 A classification based on peptide hydrophobicity versus aggregation propensity-as estimated from the simulations-resulted in a selection of promising candidates for hydrogel formation. Some of these peptides were subsequently synthesized, and indeed led to the discovery of novel hydrogels. A very recent example of HT Martini-based simulations, also aimed at finding new hydrogels, is from the work of Cardellini et al. 386 In their study, a protein-polymer compositional phase diagram was constructed based on a large set of simulations, to establish guidelines for proteinpolymer coassembly. Another recent example is from the Tuttle group, screening more than 80,000 tri/tetrapeptides on their ability to adsorb at an interface, aiming to identify peptides suitable to act as emulsifiers. 387 An HT example from a different part of materials science is the work of Alessandri et al. on bulk heterojunction morphologies of polymer/fullerene blends that form the basis of organic photovoltaic devices. 339,342 The experimental setup of solution processing to form the bulk heterojunction was mimicked by stepwise evaporation of the solvent from the system during the simulation. A systematic exploration of polymer weight and drying time on the formed morphologies, as well as the impact of polymer side chains, provided important insight into this process, with the prospect of designing more efficient light-harvesting materials in the future. More HT examples of Martini in soft matter systems can be found in Alessandri et al. 22

| Increasing complexity
The availability of many Martini molecule types and the longer length and time scales afforded by CG have been key enablers of multi-component system simulations. (Number of component copies or total system size can also be evoked, on their own, as criteria for complexity; those examples are addressed over the next sections.) In multi-component systems, a larger size is often a direct consequence because (1) the simulation box must hold all the different components, and (2) for mixtures of different proportions of components, enough molecules must be simulated so that species with the lowest concentration is still present in representative numbers. Together with size requirements, equilibration times of multi-component systems can also be longer than those of similarly sized systems with fewer components. This reflects the need to properly sample the several relevant component-component configurations, at the different length scales of the system. 388,389 Figure 9e shows the time evolution of Martini simulation records, regarding the number of nonsolvent components, since the Martini 2.0 release. All records have been of systems involving lipids-the initial focus of Martini-but not always of bilayer systems. The group of Vattulainen 390,391 has simulated realistic lipoprotein 392 particles using up to eight nonsolvent components (six lipid types and two proteins), in an application example where, incidentally, total system size does not grow with complexity. Monolayer systems such as those of lung 393 and tear fluid surfactant films 394 have also been successfully simulated with up to five protein and lipid components. Figure 9a also highlights the paradigm shift brought by Ing olfsson et al.'s simulation of a 63-component plasma membrane of realistic lipid composition. 194 That effort showcased Martini's maturity and suitability to this type of systems and kick-started an entire sub-field in simulations of realistic lipid bilayers 192,195,395,396 (of which the current complexity record, held by Corradi et al. on lipid interactions of membrane proteins, 192 is an example; Table 1, Figure 9b). These simulations also underscore the simulation size and time requirements that high complexity can entail: Ing olfsson et al. had to simulate a membrane patch of over 70 Â 70 nm 2 so that, at 0.03% mole fraction, the least represented lipid could still be present in multiple copies. 194 In this system, the evolution of heterogeneity in lateral lipid organization approached 10 2 μs timescale. 194,195 Since Risselada and Marrink's initial work 36 on lipid phase coexistence, Martini has been extensively used to study interactions with phase-separated membranes. 397,398 In these examples, complexity also arises from the number of phases that coexist, since systems are of obviously higher complexity than the same lipids in a homogeneous mixture. For these cases, the effective number of components is then perhaps best interpreted scaled by the number of phases, since simulation size and time requirements will typically be those of more complex systems. Highly complex systems can also lead to, or suffer from, subtle artifacts that may be less readily detectable than in simpler set-ups. This was recently demonstrated in the context of phase separation by Thallmair et al., 185 whereby nonconverged constraints of the Martini cholesterol model can result in temperature gradients between different membrane domains if simulation settings are not chosen carefully.
Martini has also been used to simulate systems where complexity arises less from the total number of lipids, and more from the combination of different-natured components, most notably those mixing lipids with proteins and/or  35,192,194,390,391,402,403 nucleic acids. For instance, Van Eerden et al.'s simulations of the photosystem II complex in a realistic thylakoid membrane involved modeling of 19 different protein subunits, six cofactors, five ions, and five lipid types. 399,400 Nonlipid systems have also been studied at some degree of complexity, though nowhere near that of lipids. A biological example is the Martini ribosome, 58 developed with six components (two subunits, each with protein and RNA content, plus mRNA, and tRNA), but that could conceptually be extended to include, among others, nascent chains and chaperones as was done in a mixed resolution study using the PACE force field. 401 Nonbiological applications that stand out, complexity-wise, are the ones involving phase separation and/or morphology. Here, complexity again emerges from heterogeneity rather than from shear component number. One of the earliest examples is Alessandri et al.'s three-component simulation of solvent processing of organic solar cell bulk heterojunctions. 339 More recently, enabled by the development of Martini 3, Vazquez-Salazar et al. simulated the phase behavior and oil extraction capabilities of imidazolium-based ionic liquids 86 -these extraction simulations, at six components, also contend for the most complex nonlipid Martini system to date (Figure 9d). Also leveraging Martini 3, Tsanai et al. have simulated complex coacervate aqueous-aqueous phase coexistence using a polylysine/ polyglutamate mixture, in which additional components (proteins, RNA) exhibited phase selectivity. 89 The successful examples in this section show that complexity can be confidently tackled with Martini and are a testament to the framework's transferability qualities-by which molecular models behave representatively in a range of environments, all coexisting in the same simulation box. With novel tools such as polyply, 98 it becomes increasingly easy to construct starting structures of multi-component systems in a variety of spatial geometries, as further exemplified by a PEGylated lipid vesicle containing a phase-separated polymer mixture (Figure 9c). Given this sound CG basis, we expect also nonlipid simulations to explode in complexity, paralleling the one that occurred for lipid membranes.

| Connecting to the whole cell scale
One of the key features of CG models is the connection between the atomistic detailed simulations and the mesoscale which is traditionally the realm of continuum models. Martini has proven efficient in building this bridge, allowing it to reach spatiotemporal scales of hundreds of nm and hundreds of microseconds (Table 1) while keeping a near-atomic resolution. With increasing computer power, we currently see that Martini-based simulations can truly reach the mesoscale even with modest computational resources, in contrast to all-atom simulations which require massive use of stateof-the-art facilities. [404][405][406] It is not far-fetched to expect that whole cell simulations are within reach, a prospect enabled by the progress made in the experimental characterization of cell organization fuelling integrative research efforts in this direction. 407,408 A first example of a "Martini mesoscale simulation" is given by the work of Vögele et al., 189 who performed simulations of increasingly large lipid bilayers to evaluate the effect of periodic boundary conditions on lipid diffusion in order to test the underlying assumptions of a hydrodynamic theory. They nicely showed how predictions form the theory match those from Martini simulations, converging to the continuum limit for systems of 420 Â 420 Â 100 nm box dimensions, with 132 million CG particles, the largest Martini system simulated to date ( Table 1). Simulations of cellular membranes also steadily increase in size. A first step-up in lateral dimensions was achieved by Arnarez et al., reporting on a mitochondrial membrane patch containing 45 respiratory chain proteins, about 18,500 lipids, and approaching the 0.1 μm length scale. 409 An even larger system was simulated by the Sansom group. 410 They performed a 2.5-μs simulation of a bacterial membrane patch with lateral dimensions close to 0.5 Â 0.5 μm, containing 1152 BtuB + 1152 OmpF proteins for a total of more than 22 million particles. The length scale of the simulated system corresponded to that probed experimentally, but the limited time scale proved still insufficient for a meaningful one-to-one comparison.
An important prospect of Martini is to enable molecularly detailed simulations of entire cells. Although we are not there yet, progress toward this challenging goal is being made by many groups. For instance, simulations of entire viruses have already been reported, [411][412][413][414] and protocols for setting up such simulations are becoming available. 413,415,416 An advantage of using a CG model is that many conditions can be probed, in particular with placing reasonable amounts of lipids in both leaflets and an appropriate amount of interior solvent, including ions. Initial choices often lead to instabilities which only become apparent after many microseconds of simulation. Most recently, the COVID-19 viral envelope has been simulated with Martini 417-419 (Figure 10a). In the work of the Tieleman group, 417 both SARS-CoV and SARS-CoV-2 envelopes were simulated and their dynamics and supramolecular organization explored. Structural proteins were found to form multiple string-like islands in the membrane, and clusters appeared between the heads of spike proteins with noticeable differences between the SARS-CoV and SARS-CoV-2 envelopes. The models of Pezeshkian et al. 418 represent a larger and more accurate SARS-COV-2 envelope, featuring more than 1000 proteins (25 S, 2 E, 1003 M), more than 60,000 lipids (of six different types), and up to 35 million solvent particles. Based on simulations of 1-4 μs, interesting patterns emerged of the abundant M-dimers, in good agreement with electron microscopy data. In the work of Wang et al., 419 a part of the RNA bound to the envelope was also included in their virion model, and the full structure was backmapped to all-atom resolution. Another noteworthy large-scale effort to probe viral envelope dynamics is the simulation of a HIV-I liposome, measuring 150 nm in diameter and containing about 300,000 lipids and 200 million solvent particles, performed by Bryer et al. 420 Beyond viruses, organelles are the next obvious target for detailed simulations. A first example is the autophagosome constructed by Hummer and coworkers 421 (Figure 10b). The autophagosome model features a 60-nm diameter Martini lipid vesicle filled with proteins at stoichiometries as derived from experiments. Although the proteins are represented at different (available) resolutions, and the model is just a static picture, it reveals the crowded nature of such a relatively simple organelle. Another example is the model of a full endosome measuring about 50 nm in diameter with a total of almost 3 million beads from the study of Bruininks et al. 271 Although the endosomal membrane model is lipid only, it allowed the simulation of the translocation of lipoplexed genetic material across the endosomal membrane in a 5 μs simulation. In a comparable study, Khalid and coworkers simulated the interaction of lipid liposomes representing both smooth and rough outer membrane vesicles-particles that are in reality secreted by various bacteria-with realistic models of the host plasma membrane. 422 A primary example of the capability of Martini to simulate entire organelles is the mitochondrion model from Pezeshkian et al. 130 (Figure 10c). The model is lipid only, featuring about 5 million lipids (more than 80 million particles) of different types with distinct compositions for the two leaflets of both inner and outer membranes mimicking real mitochondrial membranes. The model was simulated for a short period (2 ns) with the Dry Martini force field. Although a relatively small mitochondrion, with system dimensions of 341 Â 488 Â 792 nm, 3 it is currently the largest Martini-based system in terms of spatial dimension ( Table 1). The system was constructed from electron microscopy density maps using the TS2CG software. 130 On the way to full cell models, Tajkhorshid and co-workers recently simulated a cell-scale membrane envelope with Martini 423 (Figure 10d probed in terms of lipid and solvent content, eventually resulting in a final protocell model that proved stable during a 500 ns simulation.

| OUTLOOK
In this final section, we discuss some of the ongoing and planned developments, by us and with the help of the larger Martini community, to further improve on the Martini model. With the release of Martini 3 and the availability of many new and more accurate bead types, there is a lot of room to improve the currently available molecular models. This is not a trivial task, given the multi-objective nature of the underlying parameterization strategy. Optimizing parameters to reproduce a single target property is easily feasible, but when aiming to capture many different types of properties simultaneously, the task becomes challenging. To address this challenge, refinement of the different categories of molecules is currently organized into dedicated taskforces.
To improve lipids, the better balance between standard and smaller bead types in Martini 3 enables a finer mapping of the lipid tails, allowing one to distinguish between, for example, C16 and C18 tails, something that was not possible in Martini 2. In addition, the lipid headgroups are being revised, not only considering the new bead types that are available, but also including more sophisticated bonded potentials to better capture the headgroup flexibility in comparison to all-atom simulations. A glimpse of how this can improve lipid behavior is provided by the recently released PIP models. 174 We expect that the forthcoming generation of lipid models will improve the behavior of both pure lipid systems, such as correctly reproducing phase diagrams including liquid-order domain formation, and specific protein-lipid binding interactions.
A main limitation for proteins in Martini is the lack of an accurate description of conformational flexibility, which means most simulations use fixed secondary and tertiary structures with some flexibility for loops or specific domains, requiring specific assumptions by the investigator. Although there are several examples of proteins with relative domain motions, 383,424 including individual helices, 399,425 and proteins can sometimes be simulated in different, but fixed, conformational states, 239,262 we expect a major improvement to be made with the further fine-tuning of G o-type potentials. 177 Such potentials can be efficiently incorporated into Martini models using virtual sites. A similar approach is also being tried to fine-tune protein-solvent interactions, which appears to be required to accurately reproduce the conformational ensembles of disordered proteins. 90 Together with the obvious improvement in Martini 3 in its ability to model elements of a system at higher resolution, we foresee more applications directly aimed at the mechanism of proteins including functional aspects of protein-ligand interactions and transport mechanisms of membrane proteins. Fine-tuning of interactions to more accurately model protein-protein interactions is a less fundamental problem and has been incorporated to a significant degree in Martini 3. 83 The new version already shows a major improvement in its ability to model dimerization of transmembrane helices and the resulting structures are very close to experimental structures in a number of test cases, 83,184 promising further improvements in the treatment of membrane proteins.
For carbohydrates, the restricted choice of particle types and sizes in Martini 2 has limited its possibility to capture the chemical complexity and subtlety of carbohydrates, but with the new Martini 3 bead sizes the available universe of possible structures is much larger. Clever choices of virtual sites also allow for a better mimicking of the conformation flexibilities. Similar considerations apply to the modeling of other types of polymers. Based on the new Martini 3 bead types, we are currently constructing a library of validated topologies for a wide range of common polymers. In addition, the release of the polyply tool 98 enables users to generate topologies and starting structures for both linear and branched copolymers of arbitrary composition. The improvement of nucleic acids involves multiple challenges, including capturing backbone conformations, subtleties in base-base interactions, and accurate treatment of electrostatics, 57 all of which are likely to be improved significantly in Martini 3 with both a broader choice of bead sizes, the option to fine-tune specific interactions, and improved modeling of ions. The ultimate aim is to allow hybridization with the next generation Martini nucleotides, perhaps with the help of Go-potentials in analogy to the proteins.
In addition to the progress made in the dedicated taskforces, we are working on porting the polarizable water model to Martini 3, to account for a more accurate description of long-range electrostatic effects and systems with dielectric gradients. Also, the titratable water model is being integrated into Martini 3, and further expanded to enable constant pH simulations involving proteins. A new implementation of lambda-dynamics into Gromacs offers another route to constant pH simulations with Martini. 426 As a prospect further down the line, we are investigating the possibility to model chemical reactions on the fly with the Martini model, which would open up a large number of new applications. A proof of principle along these lines has recently been demonstrated in the case of silica polymerization reactions. 427 An important additional development is the construction of MAD, the Martini database, 155 allowing deposit and retrieval of all Martini models developed so far, with detailed version tracking, through a simple graphical interface. The MAD sever will also integrate and provide a graphical interface for the most commonly used Martini tools, such as martinize, insane, and polyply, therefore promoting adoption of the force field by new users. Finally, new and accurate automatic topology builders, properly including all the steps of Martini 3 parametrization protocol, should be of the paramount importance for high-throughput applications and for overall progress of the model.
Finally, we would like to stress the importance of the continuous development and support of other CG models. As mentioned in the Introduction, different philosophies to CGing exist, each with their own merits and shortcomings. For certain application areas, specific CG force fields optimized for the task at hand can be a better choice-examples include protein folding studies with the UNRES force field, 19 prediction of nucleotide structure and hybridization with oxDNA, 18 studying molecular fluids with SAFT-γ, 428 or simulation of polymer dynamics either with generic models 429,430 or with bottom-up constructed (structure-based) CG models. [431][432][433] Further advancement of machine learning 434,435 will certainly contribute to this development, either by producing target force fields on demand, or in the process of integrating data from various resources to optimize existing ones, including Martini. 104 Encouraging is also to see the efforts of other CG force field developers to extend their models toward broader application ranges, providing more future opportunities to cross-validate predictions from Martini.

| CONCLUSION
Looking back at two decades of development, the Martini model has taken a firm place among the established force fields for molecular simulations. This is evidenced by the large number of groups that have embraced the model, in fields ranging all the way from cell biophysics to materials chemistry. Due to its inherent speedup compared to all-atom models, and its ease of use, the model has been integrated into many workflows, to aid in the interpretation of experimental data, to unravel molecular driving forces, and in the rational design of novel materials.
With the Martini 3 version recently established, and ongoing efforts to update the basic lipid, protein, saccharide, and nucleotide parameters taking full advantage of the new parameters, we can expect a further increase in accuracy and widening of the application horizons for the coming decade. At the same time, we expect a growing awareness of the limitations of coarse-graining in general (and Martini in particular), implying the ongoing need to validate predictions by experiments, and to benefit from multiscale approaches including higher resolution models. In combination with the ever-expanding number of auxiliary tools, we foresee breakthroughs in the coming years toward whole cell modeling, drug screening, and design of supramolecular materials.

ACKNOWLEDGMENTS
Siewert J. Marrink acknowledges Gerhard Hummer for providing the image of the autophagosome, Chen Song for providing the image of the SARS-CoV2 viral envelope, Helgi Ingolfsson for providing the image of the multiscale MuMMi workflow, and both Helgi and Bart Bruininks for proofreading the manuscript. Riccardo Alessandri acknowledges Charly Empereur-Mot, Petteri Vainikka, and Najmeh Delavari for providing images for Figure 7.