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 to study systems of even greater complexity, culminating in molecularly detailed simulations of whole organelles^5,6 and the set goal of simulating entire cells^7,8. Moreover, the growing demand for computer aided rational design relies on high-throughput simulations with millions of systems simulated in parallel.^9–11 Currently, the computational demand of MD methods representing all-atoms explicitly severely limits the access to 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 coarse-grained methods is the Martini force field.^12,13 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.^12,13 This transferability-based approach allows Martini simulations to easily reach the aforementioned complexity scale. However, to really prepare Martini for the high-throughput 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 setup, run, analyze, and backmap simulations. A non-exhaustive overview can be found in our recent review on the 20-year history of Martini.¹⁴ However, the codes and scripts developed so far present a collection of separate scripts which 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 robustness of code design to bugs. In addition, it will allow 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, in order 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.^16–35 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 high-throughput drug binding assays^10,36, 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 Gō-like models.³⁹ All features are additional specific biasing steps applied after 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 which can create topologies for complex systems at all-atom (AA), united-atom (UA), and coarse-grain (CG) resolution. On top of vermouth, we built the martinize2 program, as the successor of the martinize script^16,33. The goal of martinize2 is to encompass all functionality required to generate Martini protein parameters (supporting the older versions Martini 2^16,38,40 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.^41–43 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 to foster 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 which can create topologies for linear molecules and some specific software exists that also handle 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, 6 can be found in almost all Martini programs, which generate simulation input files in the broader sense.^16,44–46 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 to 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

Figure 2

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.

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 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 to transform atomistic structure data to a coarse-grained 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. To efficiently do these comparison operations, Martinize2 relies on a number of algorithms coming from graph-theory (e.g. subgraph isomorphism), which reduces the dependance 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. 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.

After completing the repair stage everything is in place to perform the mapping to coarse-grained 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. 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 coordinates of the CG system. In order to generate the interactions that link the residues, a simple graph matching with library link definitions is done in the create-topology stage. 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 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^50–52 – also at the Martini level^53,54 – 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 3). 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 3). Subsequently, martinize2 checks if there are any modifications which 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 ε-tautomer (3 residues) and the δ-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 3

The same procedure used for setting protonation states also applies to identifying any other common post-translational modification (PTMs). 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 example (PDB 2GS2). Residue TYR845 (see Figure 4), 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 4

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 postprocessing step.^38,40 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 high-throughput 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 force-constant, 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 Martini2 DNA³³, for example. Martini2 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 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 5a and Figure 5b. 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 5b). 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 5b whereas there are in Figure 5a.

Figure 5

Furthermore, martinize2 allows to define 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 5c and Figure 5d show the FtsZ protein of E-coli as predicted by alpha-fold.^60,61 FtsZ possesses a structural unit and two disordered tail domains. With the region option martinize2 allows to generate 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 Data Base (MAD), which offers a further utility to remove certain elastic bonds selectively.⁶² We note that elastic networks can only be applied within protein molecules at the moment.

Incorporating a Ligand small-molecule Database

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. 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, which maps and parameterizes the system. Having a single step for topology generation greatly facilitates high-throughput workflows such as protein-ligand binding, one of the cornerstones of CADD.

We test this on two protein complexes. The first test case concerns Flavin Reductase (see Figure 6a), which consists of two chains that have flavinmononucleotide 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 data-base. 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 6b). 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 6

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 data-base for the M3 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 data-base, 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. 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 include 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 7). 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 coarse-grained protein in vacuum to validate that the generated structures and topology could be processed by GROMACS 2022.3.

Figure 7

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^66,67. 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 Cβ 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 clashed or extended bonds) lead to high energies in the coarse-grained systems, which causes a failure in the energy minimization routine. Likely these starting structures are also not numerically stable in a coarse-grained simulation.

As a second test case to assess the robustness of martinize2 we processed a subset of the AlphaFold Protein Structure Database.^60,61 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 8). 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 8

Discussion

In the previous section we have presented the vermouth python library for facilitating topology generation and manipulation. In order 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^70–72, 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 form 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^41,42 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 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 build 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⁵³, 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. 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 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. Especially, martinize2 is capable of detecting and using protonation states and PTMs, and capping groups automatically. In addition, martinize2 allows to fine-tune the elastic-network 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 high-throughput 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^66,67 to assign secondary structures and some specific missing atoms can prevent DSSP to work 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 is not corresponding 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 data base (~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 data-base 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 high-throughput and high complexity tasks.

Ultimately, the robustness comes at a price. Martinize2 uses a subgraph isomorphism in order to identify atoms based on their connectivity, and then issue a warning or repair 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 call GROMACS functions, so support for other MD engines can be added in 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.^34,76 The HIS-tag of 1MJ5 was removed.

All-atom simulations. For 2GS2 and 1L84 CHARMM all-atom parameters⁵⁵ were created using the CHARMM-GUI solution builder^34,76 and a small equilibration simulation (20ns) was run before the structures were converted with martinize2. The all-atom simulation used the recommended non-bonded 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 (τ = 12 ps) after initial equilibration with the Berendsen⁸⁰ barostat.

Coarse-grained simulations. All coarse-grained 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 reaction-field (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 (τ = 6 ps) and separate coupling groups for solvent and proteins. Pressure was maintained at 1 bar using the Berendsen barostat for equilibrations (τ=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.

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 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 by 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.

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.

Data availability

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

Author contributions

PCK and SJM conceived the project; PCK, JB, 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 accuracy of the models; 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.

Supplementary Information

1 – Input Parsers & Output Writers

Table 1

Table 2

2 – Related Tools

Table 1:

3 – Martinize2 Pipeline

In this section, we describe the pipeline underlying the martinize2 program in more depth highlighting the algorithms used.

Step 1 Parse input. Reading different input file formats is trivial, and all that is needed is to select the correct parser based on the file name provided. At time of writing parsers are available for pdb and gro files (coordinate files in Gromacs format). The input is commonly a list of atoms with associated properties such as atom names, coordinates and MRU (monomeric repeat unit) names. Sometimes the input also provides information about bonds in the system, such as PDB ‘CONECT’ records. These will be used if available. Otherwise, bonds will be added between the atoms based on simple geometric criteria. At the very least we require MRU names and numbers, elements, and either coordinates or bonds. In the end, the input has been parsed and transformed into an undirected graph with atoms as nodes and bonds as edges.

Step 2 Identify and Repair. To identify the parsed atoms the current generation of tools takes the combination of atom name and MRU name as leading, even though this is the most variable between models. For instance, the atom names assigned in the experimental data often do not match the atom names expected by the force field causing existing tools to either throw an error, or even produce incorrect output. We identify atoms based on their MRU names, connectivity, and their elements by overlaying the MRU with its canonical form (Figure 1a).

Figure 1:

Doing so allows us to identify deviations from the canonical structure such as PTMs, different protonation states and capping groups. In addition, this method reveals which atoms are missing in the input data, allowing us to reconstruct them. We rely on graph theory to perform the overlaying of input and reference structures (see the dedicated section on graph algorithms below).

In order to do this, every MRU in the input molecule is overlaid with its canonical reference structure with the constraint that the elements of corresponding atoms must be the same. To get the relevant canonical structure it is assumed the MRU names in the input molecule are correct and that for each MRU a corresponding block can be found in the library. If the corresponding block cannot be found an error is raised and execution is terminated. Since the library files are designed to be human readable and writable, users can add any data to the library they need.

In the best case finding the overlay is an induced subgraph isomorphism problem where M_r ⫇ R_r with M_r an MRU of the input molecule and R_r the corresponding canonical form. However, this is treated as a largest common induced subgraph problem (see below) since M_r can contain “unexpected” atoms not described by R_r such as PTMs or capping groups. If there are multiple solutions, the solution where most atom names correspond is taken. Either way, a correspondence between the input molecule and its canonical form is obtained. This correspondence is used to a) identify and add missing atoms, b) correct the atom names for the atoms that are there, and c) find which atoms are not described by the canonical MRUs. It should be noted that in this paradigm PTMs, non-standard protonation states, termini, and capping groups are all considered unexpected atoms and treated the same way.

Next, we try to identify all these unexpected atoms by overlaying them with modification template graphs from the library (Figure 1b). This is a graph covering problem where we aim to find a minimal combination of templates that covers all unexpected atoms (see below). This does mean that unless there is clear additional metadata there can be no missing atoms in the found modifications since it is not known what they should look like beforehand. The found correspondences are then used to correct the atom names. The MRUs these atoms are part of are labelled so that the correct mappings and interactions can be applied later on. In the end, the input molecule is complete, has correct atom names, and MRUs that deviate from the reference are labelled. At this point all information contained in the atom definitions in the input file and their connectivity has been used. Any atoms that could not be recognized will be removed. A warning is issued to the user if this is the case.

Step 3 Resolution Transformation. The resolution transformation step maps the input molecule to the desired output resolution (Figure 1c). We must assume that these mappings are many-to-many correspondences, and that in a mapping from e.g. AA to CG a single AA atom can be mapped to multiple CG beads. Unfortunately, this generalization prevents the use of methods developed in graph theory for this problem so far^1,2. Instead, we perform the transformation using the same type of overlay we used to identify atoms in the input molecule. This requires a ‘Mapping’ object, which consists of two molecular fragments at different resolutions, and a correspondence between their particles. These Mapping objects are taken from a library. Including this resolution transformation step in the pipeline makes VerMoUTH resolution agnostic, capable of also generating CG topologies.

Mappings from the input force field to the required output force field are taken from the library. However, since these mappings can cross MRU boundaries this is a graph covering problem. This is a variant of the exact cover problem and therefore an NP-hard problem^3,4. Because in this case it applies to the full polymer, this is intractable. We sidestep this problem by approaching it as if it were an induced subgraph isomorphism problem where all possible places a mapping fits on the input graph are found, respecting the constraints that atom and MRU names must match. In addition, the mapping may only cross MRU boundaries where it is explicitly allowed by the mapping. If mappings overlap an error is raised. For every mapping that is applied interactions from the corresponding Block are added to the output graph.

Once done, the found modifications can be mapped. First, the modifications are grouped together by connectivity with their MRUs. This is done because with multiple modifications for a single MRU their interactions may influence each other, e.g. (partial) charges in zwitterionic amino acids. Based on these groups the modification mappings that apply to most of those modifications at once are found by solving the exact set covering problem. The found modifications are then applied by finding the corresponding subgraph isomorphisms. Warnings are issued if multiple modification mappings affect the same particle or interaction.

Step 4 Create Topology. Left then is generating the topology. Generating the inter-MRU interactions by applying the appropriate Links is a series of induced subgraph isomorphism problems where all possible ways a link fits on the produced output graph are found. A link can be applied multiple times, and can overlap with other links. Whenever a link is applied the interactions it defines are added to the output graph. In addition to adding interactions, links can also change interactions already set by blocks. For example, the particle type or partial charge may depend on neighboring MRUs. Because of this, links are non-commutative and the order in which they are applied matters. To resolve this, we solve the subgraph isomorphism problems in the order the links are defined in the library (Figure 1d). At this point the output graph represents a molecule at the desired resolution with most interactions defined and coordinates can be generated. Usually these can be trivially taken from the input coordinates. However, in case atoms were missing in the input this might not be possible. In those cases, we generate coordinates based on the coordinates of the neighboring atoms.

Step 5 Post-Processing. Post-processing can consist of any number of steps, and can perform all sorts of force field specific dress-up. For example, it can create an elastic network⁵, or generate the parameters required for Gō interactions^6,7. These steps have access to the complete molecule with coordinates and canonical atom names, even if they were missing in the input, and they have access to the full topology with all associated interactions. Separating these steps out into separate Processors helps to keep them independent of each other, and allowing for any type of post-processing helps in making the program force field agnostic. There can be any number of this kind of processors depending on what was requested by the user.

Step 6 Write Output. Lastly, the output topology and coordinate files have to be written. Since this is just a matter of file formatting, this is trivial. Separating it out from the rest of the pipeline makes the program agnostic of the MD engine used. At the time of writing VerMoUTH is capable of writing Gromacs compatible topologies.

4 – Graph algorithms

Steps 2-4, which form the core of VerMoUTH rely heavily on graph algorithms, because molecules and polymers can be very naturally described as undirected graphs^8–11. In our case nodes correspond to atoms, and edges to bonds between atoms. In addition, polymers can also be described as a coarser graph, where nodes correspond to MRUs and edges to bonds between MRUs. Graph theory is a subfield of mathematics that deals with graphs, making it a particularly powerful tool in the context of this work. We primarily use methods from graph theory to identify atoms. First when curating the provided input data (Step 2), but also when performing the resolution transformation (Step 3) and when applying links (Step 4). Our primary tools for this are algorithms for finding induced subgraph isomorphisms^12–15, and for finding largest common induced subgraphs^16,17.

Largest Common Induced Subgraph. When repairing the provided molecule correspondences between the MRUs in the input molecule (M_r) and their canonical forms (R_r) are needed. In the case where M_r is not a subgraph of R_r and contains atoms that are not described by R_r, this is a largest common induced subgraph (LCIS) problem. The solution to this problem is the largest graph G that is an induced subgraph of both M_r and R_r, and the correspondences between the nodes in G and M_r; and between the nodes in G and R_r. This problem belongs to the class of NP hard problems^3,4. A possible solution to the LCIS problem is to approach it as a repeating subgraph isomorphism problem where initially G = M_r, and nodes are removed from G in a breadth-first manner until an induced subgraph isomorphism G ⫇ R_r is found¹⁶. Once a subgraph isomorphism between G and R_r is found the subgraph is not shrunk further since that would always result in a smaller common subgraph. We have based our implementation on the ISMAGS subgraph isomorphism algorithm^13,18 since, generally, molecules can be described as very sparse and (locally) symmetric graphs. The ISMAGS algorithm exploits these properties and produces only symmetrically distinct answers which reduces the runtime significantly compared to both other subgraph isomorphism algorithms, such as VF2¹³ and other LCIS algorithms, such as Koch’s¹⁷. Since our implementation of the ISMAGS is more generally applicable than just in the context of VerMoUTH we have collaborated with the authors of the popular Python graph library NetworkX¹⁹ to include our implementation.

We extended our implementation of the ISMAGS algorithm to also solve the LCIS problem in order to further exploit the symmetry breaking constraints used in the subgraph isomorphism problem. The symmetry breaking constraints are used when finding subgraph isomorphisms (see ¹³) and when shrinking the subgraph: when nodes are equivalent the node with the highest index is removed from G preferentially. In addition, to ensure common subgraphs are preferentially found using nodes with a lower index (analogous to the ISMAGS algorithm), the candidate subgraphs are sorted by their node indices. In this way we obtain good performance because in our case it is generally true that: a) there are only a few nodes not part of the reference, and b) those nodes have the highest node index. Because of this we can terminate the algorithm after the first common subgraph is found.

To demonstrate how this works we consider an example where we will try to find all LCISs between graph X and subgraph Y. The example is illustrated in Figure 2. Note that at this point the distinction between “graph” and “subgraph” is arbitrary, except for symmetry detection and performance. Nodes are represented by a letter that reflects their underlying attributes (e.g. atom type). We will consider nodes compatible if they have the same letter, and we distinguish nodes with the same letter by subscripts. First all symmetries in subgraph Y are found. This reveals A₁ to be equivalent to A₂. In the first iteration we try to find a subgraph isomorphism between X and Y (Iteration 1). Since none can be found, subgraph Y is shrunk. This yields the subgraphs in box “Iteration 2”. Since the subgraph made from the nodes {A₁, B, E, F} is symmetry equivalent to the subgraph made from nodes {A₂, B, E, F}, only the first is taken into consideration. Because no subgraph isomorphism can be found between X and any of these four subgraphs for this iteration, they are shrunk further, resulting in seven subgraphs with unique symmetries consisting of three nodes each. These are depicted in box “Iteration 3”. Of these seven subgraphs, at least one is isomorphic to a subgraph of X ({A₁, A₂, B}), therefore all subgraph isomorphisms between X and these seven subgraphs are exported in order and the algorithm is terminated.

Figure 2:

The algorithm presented is not without faults however: symmetry of X is not taken into consideration, which could reduce the search space dramatically depending on the graphs in question. In addition, some operations are performed multiple times. As an example, many of the subgraphs in Figure 2 contain the motif {A₁, B} (in bold). This results in the subgraph isomorphism algorithm reaching the conclusion that {A₁, B} is isomorphic to {A₁, B} and {A₂, B} multiple times. This can be avoided by starting the algorithm using small subgraphs, and growing them as the algorithm progresses. The results of the smaller isomorphism problems can be used to restrict the search space of the larger ones. Since in most of our cases M_r contains only a few nodes that are not isomorphic to nodes in the R_r we do not expect a (large) performance gain. It may be worthwhile to implement an adaptive algorithm that switches strategy after a few iterations of either strategy however.

Graph Covering. To identify unexpected atoms, we need to cover all those atoms with known fragments describing e.g. PTMs. We aim to find the solution where all unexpected atoms are covered exactly once, preferentially using fragments with a lower index. In VerMoUTH we sort the fragments by size so that larger fragments are used preferentially. This is a variant of the exact cover problem, making it NP hard^3,4. We solve this problem by a recursive backtracking algorithm: in order, we try to fit the fragments on the unexpected atoms until all are covered. If applying a fragment result in atoms that can no longer be covered, the solution is rejected, and the next fit is tried.

4 – AlphaFold Benchmark

The following 7 structures from the AlphaFold benchmarks produced an error, which lead martinzie2 to abort the input file generation:

AF-O80995-F1-model_v3.pdb

AF-Q58295-F1-model_v3.pdb

AF-B1GZ76-F1-model_v3.pdb

AF-A1ZA47-F1-model_v3.pdb

AF-J9VQ06-F1-model_v3.pdb

AF-F1QWK4-F1-model_v3.pdb

AF-P64653-F1-model_v3.pdb

A list of all surveyed models is available at https://github.com/marrink-lab/martinize-examples/blob/master/AlphaFoldBenchmark/surveyed_models.txt.