Introduction

Molecular dynamics (MD) has grown to be a valuable and powerful tool in studying a variety of systems in molecular detail. Advances in force fields and computer hardware have enabled the use of MD in increasingly complex systems, exemplified by recent simulations of, e.g. realistic cell membranes^1,2, virus particles^2,3, and even complete aerosol droplets⁴. However, there is a growing interest in studying systems of even greater complexity, culminating in molecularly detailed simulations of whole organelles^5,6 and the set goal of simulating entire cells^7–9. Moreover, the growing demand for computer-aided rational design relies on high-throughput (HT) simulations with millions of systems simulated in parallel.^10–12 Currently, the computational demand of MD methods representing all atoms explicitly severely limits the access to the spatial-temporal resolution needed to simulate the aforementioned systems. Coarse-grained (CG) MD methods overcome this challenge by grouping several atoms into one effective interaction site called bead and thus reducing the number of degrees of freedom that have to be simulated.

Among the most popular CG methods is the Martini force field.^13,14 Within the scope of the Martini force field about 2–5 non-hydrogen atoms are grouped into one bead. Nonbonded interactions between beads are defined in a discrete interaction table calibrated to reproduce thermodynamic data, whereas bonded interactions are matched to underlying atomistic reference simulations. Molecule parameters created following this approach are transferable between different systems and chemical contexts.^13,14 This transferability-based approach allows Martini simulations to easily reach the aforementioned complexity scale. However, to really prepare Martini for the HT and whole cell scale simulation era, automated workflows that enable fast and efficient setup of complex systems are of fundamental importance.

The Martini community has a long-standing history of easy-to-use and freely accessible scripts and programs, which helped researchers to set up, run, analyze, and backmap simulations. A non-exhaustive overview can be found in our recent review of the 20-year history of Martini.¹⁵ However, the codes and scripts developed so far present a collection of separate scripts that share no common framework or backend even though they share many common operations such as resolution transformation or mapping of coordinates. In addition, input files which define molecule parameters or fragments thereof, are not transferable between the tools, with each one of them often defining their own input file formats. We consider that unifying the operations as well as input streams into a single framework will speed up program development and also the robustness of code design to bugs. In addition, it will allow the implementation of modern software techniques such as code review, continuous integration (CI) testing, and version control which generally improve code quality and resilience.¹⁶

Designing and coding a unified framework to support general Martini software development is a massive undertaking with many facets as the original scripts and programs deal with different stages of MD simulations. To start the development, we focused the design of the framework on topology generation. A topology lies at the heart of each simulation and defines the starting coordinates as well as input parameters for the simulation. For example, to run protein simulations within Martini, a script called martinize¹⁷ takes atomistic protein coordinates, maps them to the coarse-grained resolution, and generates the protein molecule definitions from building blocks. This workflow is quite classic and underlies many scripts and programs for topology generation both at the coarse-grained but also at the all-atom level.^17–36 With the latest release of version 3 of Martini (M3), proteins have been thoroughly reparametrized.¹³ The new capabilities of M3 proteins are exemplified by their use of HT drug binding assays^11,37, which are an essential step in computer-aided drug design (CADD). Part of the improved protein properties come from the redefined Martini interaction table. However, another part of the improvement is the result of protein-specific methods such as the use of structure-biased dihedrals³⁸ (often referred to as side-chain corrections), specific elastic networks³⁹, or integration of Go-like models.⁴⁰ All features are additional specific biasing steps applied after the generation of the original topology file for the protein and are not part of the capabilities of the previous martinize script. Hence, we choose to co-develop a unified framework for topology generation together with a new martinize version, martinize2.

In this paper, we present the VERsatile MOdular Universal Transformation Helper (vermouth) library, a general python framework aiding in the design of programs that can create topologies for complex systems at all-atom (AA), united-atom (UA), and CG resolution. On top of vermouth, we built the martinize2 program, as the successor of the martinize script^17,34. The goal of martinize2 is to encompass all functionality required to generate Martini protein parameters (supporting the older versions Martini 2^17,39,41 as well as the latest Martini 3) and be compatible with high-throughput workflows as needed in CADD approaches based on Martini. In addition, both vermouth and martinize2 are designed to have sufficient flexibility and robustness to ready Martini for the era of high-throughput high complexity simulations.

Finally, we note that much of the progress of Martini has resulted from an active community of researchers contributing scripts, programs, and parameters. However, as is the case for most research software in the field they often fail to adhere to the principles of FAIR: findability, accessibility, interoperability, and reusability.^42–44 The FAIR principles⁴³, originally designed to improve data management and reproducibility in science, have recently been extended to research software in a more general sense. This extension is aimed at fostering more sustainable software development in science.⁴² To meet these standards the software tools we present here are distributed under the permissive open-source Apache 2.0 license on GitHub and are developed using contemporary software development practices, such as continuous integration testing. To make adoption as easy as possible, they have few dependencies, are distributed through the Python Package Index, and can be installed using the Python package manager pip. Other researchers are encouraged and welcome to contribute parameters and code as outlined in our contribution workflow.

Results

In this section, we first outline the design and API of the vermouth library. Then we discuss how the vermouth library is used to construct a pipeline for generating protein input parameters for the Martini force field. This pipeline constitutes the new martinize2 program. Finally, we present some benchmarks and selected test cases to demonstrate the capabilities of martinize2 and assess its fitness for generating complex system topologies and high-throughput workflows, surpassing the capabilities of the previous martinize script.

The Vermouth Library

In the literature, many scripts and programs have been described that can create topologies for linear molecules and some specific software exists that also handles branched molecules such as carbohydrates²⁴, or dendritic polymers²⁵. However, to the best of our knowledge, there is at present not a general program that can create topologies from atomistic structures for any type of system, and at any resolution, presenting an extendable and stable API. Based on the existing software, we can, however, define a number of required and desirable features for such a general program and library to have: 1) it must be force field and resolution agnostic; 2) it must be MD engine agnostic; 3) it must use data files that can be checked, made and modified by users, and 4) it must be able to process any type of molecule or polymer, be it linear, cyclic, branched, or dendrimeric, and mixtures thereof.

To start designing a library that can fulfill the above requirements we note that most workflows used for topology generation can be decomposed into six fundamental stages (Figure 1): First, reading input data, usually an atomistic coordinate file (e.g. from the protein data bank); second, identifying the parsed atoms, to find how they correspond to the atoms in the data files describing the building blocks; third, optionally a resolution transformation step; fourth, the generation of the actual topology and assigning particle types and bonded interactions; fifth, any type of post-processing; and finally, sixth, writing the required output files. Even though these stages are generally shared for topology generation pipelines, they also apply to other workflows commonly encountered in Martini programs. Especially, stages 1, 3, 5, and 6 can be found in almost all Martini programs, which generate simulation input files in the broader sense.^17,45–47 Separating these stages, therefore, helps to define an API with data structures and independent processes, which optimally support such workflows. In addition, the clear distinction in stages helps to externalize any data files, which can be edited by the user or force field developers. Vermouth is built on the idea and definition of processors, which are tasks arranged in a pipeline. This design was inspired by the ubiquitous workflow managers available in the field.⁴⁸ We formalize the idea of processors by introducing an abstract base class the Processor. New pipeline stages can be created as subclasses of this base class. All Processors operate on the central data structure class System, which contains any number of Molecule data structures (see Figure 2). A Molecule is defined as the graph of a molecule or assembly of molecules, which are connected by bonded interactions. The nodes of a Molecule usually correspond to atoms or coarse-grain beads but can be any form of particle as defined by the force field.

Figure 1

Nodes can have attributes that describe additional information such as a residue name or charge. However, only the atom name, residue name, and residue number are required as attributes. In addition, the edges of the Molecule follow the connectivity as defined by bonds, angles, or other bonded interactions. For example, two protein chains connected by a disulfide bridge would be considered a single Molecule. In contrast, a cofactor, which is only interacting via non-bonded interactions, would be a separate Molecule. Operations on Molecules usually add or remove bonded interactions or change node attributes. For convenience, Processors can also operate on a collection of molecules, which are defined by the System class (see Figure 2). A list of all available processors is given in the documentation.

Figure 2

Processors operate on Molecules. However, often additional data is required to perform the pipeline as defined by the Processor. The additional data can be provided in the form of one of the four other main data structures (Blocks, Links, Modifications, Mappings) or arguments of the processors that can be set in a script or via the command line interface. These four other data structures contain all molecular level information required to fully define and/or modify a topology for any type of MD code (e.g., atom types, bonded interactions, and positions) as well as enable transformations between topologies. For example, a Mapping consists of two molecular fragments at different resolutions and a correspondence between their particles. In contrast, Blocks, Links, and Modifications are graphs, which describe these molecular fragments, the links between those, and possible changes to fragments respectively. They are all subclasses of a Molecule and an extension of the graph class from the networkx library⁴⁹ (see Figure 2).

As shown in Figure 2, to make the data structures that are force field specific (Blocks, Links, Modifications) easier to use, vermouth offers a second collection class called a ForceField. Every molecule must have a ForceField associated with it. Additional information on the data structures is given in the documentation.

Finally, the vermouth library also contains a number of parsers that return instances of the data structures from common input file formats. For example, the in-house ff format defines Blocks, Links, and Modifications, while the backwards style mapping format can be read to return an instance of the Mapping class. Table S1 in the Supporting Information summarizes all input parsers as well as the format and data structure they return. We note that vermouth is also able to read content associated with the ‘[molecules]’ directive of the GROMACS topology file, which is colloquially referred to as itp file. This allows to directly manipulate GROMACS molecule files within vermouth. We note that as neither parsing nor the Molecule itself depends on GROMACS code, the library can easily be extended to other MD engines.

Martinize2

Martinize2 is a pipeline constructed of vermouth Processors with a command line interface (CLI), with the purpose of transforming atomistic structure data to a CG Martini topology including both coordinates and simulation parameters. Martinize2 is the successor of the martinize script, which was used for generating input parameters for Martini version 2 proteins, DNA, or RNA. However, different branches had to be used for proteins and DNA (martinize.py¹⁷, martinize-dna.py³⁴) or RNA³⁶. In contrast, martinize2 is designed to generate topologies for the Martini force field for proteins, DNA, and in principle any other arbitrarily complex molecule.

Martinize2 consists of different Processors which fulfill the basic stages of topology generation as shown in Figure 1. We note that the design of martinize2 is general and applies to arbitrarily complex polymers consisting of arbitrary monomeric repeat units (MRUs). However, to increase the readability of the following sections the layout of the program is described in terms of residues in proteins.

The martinize2 pipeline starts by reading an atomistic structure, which describes a single molecule (e.g., protein) or assembly of any size. Subsequently, bonds between the atoms are inferred either by distance calculation, atom names within residues, or using CONECT records of the PDB file. All atoms that are connected by bonds form a Molecule. Thus, martinize2 creates a System of Molecules at the atomistic resolution at the end of the input reading stage. In stage 2, Identify and Repair, each residue of each molecule is compared against its canonical definition. Canonical definitions are selected by residue name from the library files. This comparison identifies missing or additional atoms on a residue and fixes all atom names (Figure 3a). To efficiently do these comparison operations, martinize2 relies on a number of algorithms coming from graph theory (e.g. subgraph isomorphism), which reduces the dependence on accurate atom names, since these occasionally differ based on the source of the input structure. Which algorithms are used in the code is described in more detail in the Supporting Information. Once the residues have their canonical atom names, martinize2 checks if the missing or additional atoms are described by any of the modification files (Figure 3b). Modifications describe changes in residues from their canonical form, e.g. different protonation states, termini, or post-translational modifications (PTMs), and the effect these have on the topology.

Figure 3

After completing the repair stage everything is in place to perform the mapping to coarsegrained resolution. The mapping descriptions are read from the mapping input files in the library and tie together residue definitions at the all-atom and CG level and the correspondence between them (Figure 3c). Mapping to CG level in martinize2 is done with a multistep subgraph isomorphism procedure, which is general enough to cover edge cases such as when mappings span multiple atomistic residues. A detailed description is provided in the Supporting Information. The mapping provides a System of Molecules at the CG level. These molecules already define all bonded interactions within the residues as well as the coordinates of the CG system. To generate the interactions that link the residues, a simple graph matching with library link definitions is done in the create-topology stage (Figure 3d). Finally, after that, we end up with the full CG topology, which is ready for post-processing steps. Post-processing summarizes all biases and modifications that have to be done on the CG molecule and its CG coordinates. For example, an elastic network is needed to keep the tertiary structure of the protein and is applied in the post-processing stage. Finally, martinize2 writes out the CG coordinates and the CG topology file that are production-ready.

Custom Protonation States and PTMs

Of the 20 common amino acids, there are four (GLU, ASP, LYS, HIS), which can readily change their protonation state as a function of pH or environment. Whereas commonly those amino acids are still considered to be in their pH7 protonation state, it is more appropriate to determine their local pKa from for example continuum electrostatics.⁵⁰ Subsequently the appropriate charge of the amino-acid can be determined from that pKa and set for the simulation. Even though recently more advanced methods became available for dynamically treating protonation states^51–53 - also at the Martini level^54,55 - the fixed charge approach is still the most common and for Martini most computationally efficient. However, the previous martinize version lacked the functionality to treat protonation states for all amino acids. Only histidine protonation states could be set interactively but only for two of three possible protonation states.

Other protonation states as defined by the atomistic structure coordinates or residue names were ignored without warning. In addition, the interactive setting of protonation states becomes very cumbersome for large protein complexes.

To overcome this problem and make protonation state handling easier and more robust, we utilize a dual strategy in martinize2 to identify and correctly set the protonation states (see Figure 4). In route a) the user provides atomistic structure coordinates with AA residue names including those of non-default protonation states corresponding to the naming conventions used in CHARMM⁵⁶ or AMBER⁵⁷. Protonation states can be obtained from online servers such as H++⁵⁸ or propKa⁵⁹, for example. If the residue names are correctly given, they can be matched against the parameters in the library and the CG residue obtains the correct protonation state. In the alternative route b), the residue name is simply that of the default pH 7 amino acid, however, the structure file contains an additional hydrogen. In the repair and identify step the chemical graph of the amino acid is compared to the building blocks in the library and any unexpected atoms are flagged. For example, in the case of protonated histidine, the additional hydrogen is labeled (see Figure 4). Subsequently, martinize2 checks if there are any modifications that would match the complete input graph if added to the original block. In the Histidine example, the modification contains the additional hydrogen which together with the original histidine block make up a protonated histidine. The modification also changes the mapping such that the correct protonation state is set at the CG level. This route is more appropriate for example when processing crystal structure files, which are not necessarily named according to any force field convention. We have tested this feature on two protein structures taken from the PDB (1MJ5, 3LZT) and processed as described in the Methods section. In 1MJ5 there are six Histidine residues of which one is predicted to be charged at pH 7. The others are neutral. However, they are divided between the E-tautomer (3 residues) and the 5-tautomer (2 residues). Martini 3 parameters are different for the two tautomers in contrast to Martini2, which is accordingly recognized by martinize2. In addition, for lysozyme, we have considered residue GLU35 protonated, which would be appropriate at a pH of 6 or less. For both examples, the appropriate protonation states and tautomers are generated at the CG level.

Figure 4

The same procedure used for setting protonation states also applies to identifying any other common PTM. Using this procedure, lipidation, phosphorylation, amination or acetylation can be taken into account automatically. To demonstrate that martinize2 can handle PTMs, we have implemented dummy parameters for testing of Tyrosine phosphorylations in the M2 force field and generated a Martini topology for the EGFR kinase as an example (PDB 2GS2). Residue TYR845 (see Figure 5), which is located in the activation loop of the EGFR kinase, is phosphorylated when the kinase is activated.⁵⁶ Martinize2 was able to convert the structure in one go to M2 resolution. We note that at the time of writing the M3 force field is lacking parameters for these PTMs and they are therefore not implemented in martinize2 yet. In this case, a warning is issued by the program.

Figure 5

Expanding the Options of Elastic Network Fine-Tuning

Due to the limitations in most coarse-grained protein models (e.g. lack of explicit hydrogen bonding directionality), the tertiary structure has to be enforced with a structural bias called elastic network (EN).⁶⁰ An EN for Martini proteins consists of weak harmonic bonds between backbone beads of residues (within a chosen cut-off distance) and is generated after the resolution transformation as a postprocessing step.^39,41 Martinize offered only two types of EN options, the regular model and the Elnedyn³⁹ approach, both of which are also implemented in martinize2. However, as the EN fixes the tertiary structure, changes in the structure upon, e.g., ligand binding are not captured. To improve protein models in this aspect recently Go-like models have been applied to Martini.⁴⁰ In a Go-like model the harmonic bonds are substituted by custom Lennard-Jones (LJ) interactions that can dissociate, thereby allowing for some tertiary structure changes. Within the scope of Martini, a workflow is available to replace the elastic network with a Go model that is generated from a contact map.

Even though Go models offer better flexibility, they are currently limited to single monomeric protein units and require some fine-tuning to get the optimal performance.⁴⁰ Especially, for HT workflows the EN approach is therefore the preferred option. To further improve upon the elastic networks generated by the old martinize, martinize2 offers several options to fine-tune the EN and get better behavior within the constraints of the EN approach. Besides the cut-off and forceconstant, martinize2 now implements a residue minimum distance (RMD). The RMD is defined as a graph distance and dictates how far residues need to be apart in order to participate in elastic bonds. Defining the RMD as a graph distance means that no bonds are generated between residues that are for example bound by a disulfide bridge. It thus presents a more rigorous implementation than in the previous version. Usually, the residue minimum distance is 3 in order to avoid the EN competing with the bonds, angles, and dihedrals between the backbone beads.

We note that this is part of the Martini protein model and should not be changed. Additionally, martinize2 allows to select which beads to generate the EN between. This option is needed for Martini 2 DNA³⁴, for example. Martini 2 DNA offers a stiff EN version, where also sidechain beads are included. Furthermore, martinize2 allows to define where in the protein to apply the elastic network. This is done with the EN unit option. The EN unit can be a molecule, chain, all, or ranges of residue indices. The most trivial option is all in which case an EN is applied between all protein molecules in the system. The option molecule and chain yield the same network, if distinct molecules are also distinct chains. However, when two chains are connected by a disulfide bridge, for example, they would be one Molecule in the martinize sense. On the other hand, if the interface is not very well defined or more flexible, biasing the two chains separately could improve the EN. In that case, the chain option can be used. This use-case is shown for the human insulin dimer in Figure 6a and Figure 6b. The human insulin dimer consists of two chains, which are connected by two disulfide bridges. If the molecule or all option is used an EN is generated within the chains and between the chains (Figure 6b). However, to avoid generating the EN between the two chains the chain option can be supplied in which case the EN is only generated within chains. As the zoom-in on the tail part shows there are no more bonds between the two chains in Figure 6b whereas there are in Figure 6a.

Figure 6

Furthermore, martinize2 allows the definition of regions of residue IDs where an EN should be generated. This feature gives maximum flexibility and allows to bias structural regions of proteins whereas an EN in intrinsically disordered regions (IDRs) can be avoided. For example, Figure 6c and Figure 6d show the FtsZ protein of E-coli as predicted by alpha-fold.^61,62 FtsZ possesses a structural unit and two disordered tail domains. With the region option, martinize2 allows the generation of an EN only for the structural domain. Within the old martinize superfluous bonds needed to be removed manually.

Finally, we note that martinize2 is now implemented in the Martini Database (MAD) GUI, which offers a further utility to remove certain elastic bonds selectively.⁶³ We note that this option is only available for protein molecules at the moment.

Beyond Proteins: Incorporating other Molecules in Martinize2

Legacy martinize is only applicable to one category of molecule (i.e. proteins or DNA), which is one of its biggest drawbacks even for setting up simple protein simulations. Martinize2 allows the inclusion of new classes of molecules without adjusting the codebase. For instance, proteins frequently have other molecular units associated such as ligands, cofactors, metal ions, or lipids. The general workflow of martinize2 allows us to convert these molecules in one go provided that the library files are present. In this way, no post-hoc step is needed, that maps and parameterizes the system. Having a single step for topology generation greatly facilitates HT workflows such as protein-ligand binding, one of the cornerstones of CADD.

We test this feature on two protein complexes. The first test case concerns Flavin Reductase (see Figure 7a), which consists of two chains that have flavin mononucleotide ligands (FMN) and one NAD cofactor bound (2BKJ). M2 parameters and mappings from the GROMOS force field were previously published.⁶⁴ Parameters and mappings have been added to the vermouth database. Subsequently, the system could be converted in one step. During a short simulation, the cofactors remain well bound, indicating that no inappropriate parameters or faulty geometries were generated. Next, we created topologies and starting structures for Lysozyme with a benzene molecule bound (1L84), using the M3 force field (Figure 7b). The protein and ligand were again converted in one step and then simulated for a short period. As previously the ligand stays bound, showing that the protocol generates reasonable starting structures and correct parameters.

Figure 7

To fully leverage this new feature, ligand data files are required to be present. Thus, we implemented mappings and parameters from a previously published small molecule database for the Martini 3 force field.⁶⁵ The set comprises 43 small molecules, which are often part of drugs or drug precursors. All small molecules have corresponding parameters in the CHARMM ligand database, which allows users to directly convert atomistic CHARMM simulations to Martini. Mapping directly from crystal structures as present for example in the PDB or other force fields is also possible. In these cases, the residue names may have to be adjusted to be the same as in the CHARMM naming convention.

In addition to creating topologies for linear biopolymers, martinize2 is now also able to handle topologically more complex molecules. For example, crown ether (Figure 7c) consists of six polyethylene glycol (PEG) repeat units and is cyclic. To test whether martinize2 can handle cyclic molecules of multiple repeat units it was converted to Martini2 resolution applying the latest PEG parameters.⁶⁶ The second example is branched polyethylene (PE), where we chose a sequence that begins with two linear units followed by three branched ones and two linear units after. Also, this molecule is converted to Martini2 resolution⁶⁷ by martinize2 (Figure 7d).

Finally, we have set up instructions on how researchers can submit parameters to the database allowing it to grow and support other researchers. In addition, martinize2 facilitates dynamic linking of citations to parameters. With this mechanism, citations are printed at the end of the run that dynamically includes citations to all parameters used in the final topology. Such a system also allows researchers to easily receive credit for contributed parameters.

Complexity Benchmark

To assess the robustness of martinize2 in a high-throughput use case, we processed the template library used by the I-TASSER⁶⁸ protein prediction software (Figure 8). At the time of download (26 March 2021), the dataset contained 87084 protein structures. We processed each of these structures with martinize2 to get M2.2 models with elastic networks. We then minimized the energy of the CG protein in a vacuum to validate that the generated structures and topology could be processed by GROMACS 2022.3.

Figure 8

Of the 87084 structures in the dataset, 63680 (73%) could be processed through the whole workflow without error. The main cause of failure (25% of the structures) was missing coordinates in the input structures. When all the atoms that compose a bead are missing from the input, martinize2 can generate a topology but it cannot generate coordinates for the bead. Note that if only some atoms are missing, then vermouth does estimate the position of the bead. 876 structures (1%) had missing coordinates in the backbone that prevented the use of DSSP^69,70. Finally, 802 input structures (1%) had at least one residue that was inconsistent with the library. Upon further inspection, most of these structures contain malformed glycine residues with an unexpected Cp atom. Martinize2 detected these inconsistencies and emitted a warning for each of them; warnings can be explicitly and selectively ignored, if they are not no output is written to avoid subsequent workflow steps working with corrupted files.

All the 63680 input structures that were successfully processed by martinize2 could be processed by the GROMACS pre-processor (grompp). However, 13 structures failed the energy minimization. A visual inspection of some of these failing inputs shows the input atomistic structures can be problematic. Erroneous interatomic distances (steric clashes or extended bonds) lead to high energies in the CG systems, which causes a failure in the energy minimization routine. Likely these starting structures are also not numerically stable in a subsequent simulation.

As a second test case to assess the robustness of martinize2, we processed a subset of the AlphaFold Protein Structure Database.^61,62 200,000 randomly chosen unique protein structures (see Supporting Information) were given to martinize2 and subsequently an energy minimization was performed, if the structure could successfully be converted to coarse-grained representation. Of the 200.000 structures in the dataset, 7 structures (see Supporting Information) raised an error during the conversion step. Upon further (visual) inspection of the problematic structures, we concluded that all errors were caused by inaccurate initial atomistic coordinates. These inaccurate atomic positions caused bonds to not be identified or additional superfluous bonds to be detected (Figure 9). In these cases, the unrecognizable residues were detected and caused martinize2 to emit a warning. The remaining 199,993 successfully converted structures could be processed by the GROMACS pre-processor (grompp) and it was possible to perform an energy minimization.

Figure 9

Discussion

In the previous section, we have presented the vermouth python library for facilitating topology generation and manipulation. For researchers to use vermouth as a framework for software development it presents a clear API separated into data structures, parsers, and processors. In addition, the library relies on only three permissibly licensed open-software projects namely numpy⁷¹, scipy⁷², and networkx⁴⁹. This allows researchers more freedom in licensing their code and reduces the potential for bugs introduced by dependency changes. Furthermore, the central data structure represents molecules as graphs. Representing molecules as graphs allows to leverage algorithms from graph theory. Using graph theory for many of the workflows underlying the Processors makes them faster and more robust towards edge cases. Even though applying graph theory to molecules is not a new idea^73–75, vermouth is specifically designed to also handle coarse-grained level molecule transformation focusing on the Martini force field.

Therefore, vermouth presents additional functionality often lacking from other packages. For example, the handling of virtual-sites, which are ubiquitous in many M3 molecules, is rigorously handled in all Processors. As another more general example, the Processor applying interactions between residues can automatically compute structural biases from the mapped coordinates. Finally, the vermouth library adheres to the FAIR principles^42,43 to allow adoption by non-experts and ensure quality control. In particular, for both the vermouth library and martinize2, continuous integration testing is implemented and code review is required. The software is also semantically versioned, and it is distributed through established channels, most notably the Python Package Index, and hosted openly on GitHub.

We have shown how vermouth was used to shape the martinize2 program. However, martinize2 is not the only program leveraging the power of the vermouth library. The polyply python suite is another library and collection of command line programs built upon vermouth. Polyply enables users to generate both atomistic and coarse-grained simulation input data, i.e. structures and topologies, from sequence information. As such, it allows building system coordinates for arbitrarily complex macromolecular systems and nanomaterials.⁴⁵ Furthermore, the martini-sour package^54,76, which is currently under development, utilizes vermouth to convert topology files from regular Martini to titratable Martini simulations. These examples already illustrate that vermouth has the potential to indeed become the central framework for Martini software development and possibly for other scientific software developments.

Martinize2 enables researchers to prepare simulation input files for arbitrary (bio)polymers, starting from atomistic structure. We have shown in-depth examples focusing on protein-specific applications, given that they are the most important target for martinize2. However, also more complex non-biological molecules such as a cyclic crown ether and branched PE were showcased to demonstrate the capabilities of martinize2. Furthermore, the user has complete control over the data files used. The abstraction of force field data into Blocks, Modifications, and Links allows researchers to reason about model intricacies in a structured manner. This helps the development of optimized models and parameters for complex (polymeric) molecules, as well as clearly defining in which combinations these are validated. The new program uses algorithms from graph theory to identify atoms and assign the appropriate interactions. This makes the program more tolerant towards its input so that the users have to worry less over details such as atom names, or ensuring that all residues are in order and appropriately numbered. In particular, martinize2 is capable of detecting and using protonation states and PTMs and capping groups automatically. In addition, martinize2 allows to fine-tune the EN and as it is not limited to proteins can also generate parameters for ligands, cofactors, or lipids.

In practice, there are decisions a user needs to make when using vermouth and martinize2, especially for HT pipelines. Martinize2 detects but does not reconstruct atoms that are missing from the input structures; these missing atoms can have adverse effects on the result. In the most harmless cases, they only shift the position of a particle in the output structure. When all the atoms for a particle are missing, then the program cannot compute a position for that particle leading to an incomplete output where a particle does not have coordinates. Also, some workflows depend on DSSP^69,70 to assign secondary structures and some specific missing atoms can prevent DSSP from working properly. In those cases, martinize2 issues a warning whenever it cannot automatically take care of pitfalls. Handling of these cases is a central difference between the new and old version. The old version either terminates with an undefined error or, probably worse, runs and gives output that does not correspond to the atomistic structure given as input. To illustrate the robustness of martinize2 towards problematic input, we applied the program to the complete I-TASSER database (~87k structures) as well as a subset of the AlphaFold Protein Structure Database (~200k structures). For the two benchmark cases, martinize2 was able to issue a warning or error for all structures, which contained seriously malformed residues. Of the first database only 13 structures failed in the energy minimization due to problematic starting coordinates but not obviously malformed residues. In the second benchmark set only 7 seriously malformed residues were identified, and all other structures were successfully energy minimized. Thus, we consider martinize2 more robust and fit for HT and high-complexity tasks.

Ultimately, the robustness comes at a price. Martinize2 uses a subgraph isomorphism to identify atoms based on their connectivity, and then issues a warning or repairs the input. However, subgraph isomorphism is an NP-complete problem⁷⁷. As a result, martinize2 is significantly slower than martinize. Nevertheless, considering the flexibility the new program offers, in addition to the fact that it is still fast enough to process all entries in the I-TASSER data bank ⁶⁸, this is deemed to be acceptable. Even though martinize2 will most likely never be as fast as martinize we note that many of the processes can still be optimized to yield further performance increases. Aside from the performance limitations, vermouth, and martinize2 present some other limitations as well. Both martinize2 and vermouth are currently only capable of writing topologies in GROMACS format. However, our library does not use the MD parameters of the produced topologies or calls GROMACS functions, so support for other MD engines can be added in the future. In addition, since vermouth defines an API, it could even be integrated with existing software such as OpenMM.⁷⁸ Furthermore, the processor pipeline underlying martinize2 is currently hardcoded. Future improvements will focus on making the workflow defined by martinize2 more flexible, in order to include the processor pipeline as part of the force field definitions. This would enable the use of different pipelines for different force fields, allowing for easier force field-specific post-processing.

Methods

Preparation of protein input files

Crystal structures were obtained from the RCSB for the following proteins (3LZT; 2GS2; 2BKJ; 1L84; 3I40; 3IGM, 1MJ5) or the Alpha Fold Data Bank⁶¹ for FtsZ with the ID A0A7Y6D765. Hydrogens and missing heavy atoms were reconstructed using the PRAS package, if appropriate.⁷⁹ For 3LST and 1MJ5, the pKa and half-way titration point were estimated using the propka package.⁵⁹ For 3LST the GLU35 was protonated using the CHARMM-GUI solution builder.^35,80 The HIS-tag of 1MJ5 was removed.

AA simulations

For 2GS2 and 1L84 CHARMM AA parameters⁵⁶ were created using the CHARMM-GUI solution builder^35,80 and a small equilibration simulation (20ns) was run before the structures were converted with martinize2. The AA simulation used the recommended nonbonded force settings as for CHARMM with GROMACS. ⁸¹ The temperature was maintained using the v-rescale thermostat by Bussi et al.⁸² at 310K and pressure was maintained at 1 bar using the Parrinello-Rahman⁸³ barostat (t = 12 ps) after initial equilibration with the Berendsen⁸⁴ barostat.

CG simulations

All CG MD simulations were run using GROMACS 2021.5⁸⁵ and the recommended mdp parameters for Martini 2⁸⁶ and Martini 3¹³ respectively. In particular, the Lennard-Jones interactions were cut-off at 1.1 nm and electrostatics were treated with reactionfield (cut-off 1.1 nm, dielectric constant 15). The time-step was 20 fs in all cases and the production trajectories were run with the standard leap-frog integrator. Temperature was maintained using the v-rescale thermostat by Bussi et al.⁸² at 310K with (t = 6 ps) and separate coupling groups for solvent and proteins. The pressure was maintained at 1 bar using the Berendsen barostat for equilibrations (t=6ps). The initial systems were solvated using the polyply⁴⁵ package or gmx solvate utility.

Complexity benchmark

The (Swiss-Prot) subset of the AlphaFold protein structure database used for the complexity benchmark contained 542.378 pdb structure files at the time of download (22 December 2022). The testing pipeline we used was written in Python and randomly picked 200.000 structures which were given to martinize2. Possible errors during conversion or the subsequent grompp and energy minimization steps were captured.

Code availability

All code can be found online at https://www.github.com/marrink-lab/vermouth-martinize. In addition, all released versions are also published on the Python Package Index at https://www.pypi.org/project/vermouth. The documentation is available at https://vermouth-martinize.readthedocs.io/en/latest/index.html.

Data availability

Input files and commands required to reproduce the example test cases from this paper are available on GitHub at https://github.com/marrink-lab/martinize-examples. MD trajectories and benchmark data are available upon reasonable request from the corresponding authors.

Acknowledgements

We would like to thank all users that tested the development versions and provided valuable feedback, in particular the members of the SJM group and the participants of the Martini Workshop 2021. We also thank Melanie König for her feedback on the manuscript and figures. Work is supported by an ERC Advanced Grant (“COMP-O-CELLMIC-CROW-MEM”) to SJM. PCTS acknowledges the support of the French National Center for Scientific Research (CNRS) and the research collaboration with PharmCADD. JB acknowledges financial support from the Agencia Estatal de Investigación (Spain), the Xunta de Galicia - Consellería de Cultura, Educación e Universidade (Centro de investigación de Galicia accreditation 2019–2022 ED431G-2019/04 and Reference Competitive Group accreditation 2021–2024, CÓDIGO AXUDA). The European Union (European Regional Development Fund - ERDF) and the European Research Council through consolidator grant NANOVR 866559.

Author contributions

PCK and SJM conceived the project; PCK, JB, and FG implemented the described software; PCK, JB, TAW designed the program structure; PCTS & FG designed the benchmark tests used along the development of the code to guarantee the accuracy of the models; MvT ran the protein benchmark test-case using the alpha-fold data-base, while JB ran and analyzed the iTasser data-base; FG ran and analyzed all other test-cases. PCTS helped to implement the force field files, and managed feedback from beta testers; PCK and FG wrote the manuscript, with contributions from all authors. SJM provided guidance and supervision in the project.

Competing interests

The authors declare no competing interests.