Many-Body Dissipative Particle Dynamics with the MARTINI “Lego” Approach

MARTINI is a popular coarse-grained (CG) force-field that is used in molecular dynamics (MD) simulations. It is based on the “Lego” approach where nonbonded interactions between CG beads representing chemical units of different polarity are obtained through water–octanol partition coefficients. This enables the simulation of a wide range of molecules by only using a finite number of parametrized CG beads, similar to the Lego game, where a finite number of brick types is used to create larger structures. Moreover, the MAR-TINI force-field is based on the Lennard-Jones potential with the shortest possible cutoff including attractions, thus rendering it very efficient for MD simulations. However, MD simulation is in general a computationally expensive method. Here, we demonstrate that using the MARTINI “Lego” approach is suitable for many-body dissipative particle (MDPD) dynamics, a method that can simulate multi-component and multi-phase soft matter systems in a much faster time than MD. In this study, a DPPC lipid bilayer is chosen to provide evidence for the validity of this approach and various properties are compared to highlight the potential of the method, which can be further extended by introducing new CG bead types.


Introduction
Molecular dynamics [1] (MD) has been established as an indispensable tool for scientific research in a wide range of areas, such as solid state physics, soft matter, a e-mail: carnevale@ifpan.edu.plb e-mail: panos@ifpan.edu.plfluid dynamics, and biophysics.The method has also become popular due to the existence of well-documented, open-source, free, parallel software, such as LAMMPS [2,3] and GROMACS [4], which are used and supported by vast and vibrant scientific communities.In the heart of a successful MD simulation lies the forcefield (model), which describes the interactions between the particles (or group of particles/atoms) in a system.The force-field should be able to faithfully reproduce specific properties (e.g., surface tension) of the system under specific conditions (e.g., a certain temperature or pressure).For example, in the case of systems with surfactants, the viscous and surface-tension forces are important properties that the force-field should capture well, in order to describe phenomena, such as droplet coalescence [5].Therefore, it is natural that much of research in MD (MD force-fields are also directly suitable for the Monte Carlo method [6]) concerns the development of high-quality force-fields that renders a simulation reliable to reproduce system's properties.
However, MD simulations are generally computationally expensive, especially in the case of all-atom force-fields.This has led to the quest for more computationally efficient force-fields, where a single interaction site would correspond to a group of atoms, namely coarse-grained (CG) models [7][8][9][10].In the case of CG force-fields the number of interacting sites that represent the system is therefore reduced, which translates into less computations, and in turn importantly a much smaller carbon footprint of simulations.To this end, one of the most popular CG force-fields is MARTINI [11][12][13][14], which has been proposed almost two decades ago and is still under development with the MARTINI 3.0 [15] recently released.The MARTINI force-field can be used to simulate a wide spectrum of systems, such as proteins [16][17][18][19], polymers [20], carbohydrates [21], glycolipids [22], glycans [23], DNA [24], RNA [25], water [26] and various solvents [27], while various extensions include simulations for specific pH [28] and more recently chemical reactions (reactive MARTINI) [29].Without going into the nitty-gritty details of the MAR-TINI force-field here (some of them will be discussed in the next Section), the reasons for its popularity are manifold.Most importantly, it is versatile and transferable, thus enabling the simulation of a wide range of systems, as mentioned above.This is mainly due to the "Lego" approach, where beads with different polarity can be used to represent various chemical units, which in turn are used to build the various molecules of the system, that is complex soft matter systems of different components.Such an approach might hinder a detailed representation of a chemical unit by a CG bead, but this might not always be required for the CG modeling of a system.MARTINI also uses the shortest possible cutoff of Lennard-Jones interactions including attractions, which renders MARTINI computationally efficient for MD.Moreover, it is actively developed and tested by a large community with information, examples and applications being abundant in the literature and online.MARTINI has therefore been an established forcefield for MD simulation.However, another computational method beyond MD could probably offer further possibilities for faster simulations.Importantly, this method should include attractive interactions in order to enable the simulation of multiphase systems.Here, we show that the many-body dissipative particle dynamics (MDPD) [30][31][32][33][34][35] method can be used to substitute MD simulations in certain areas, and, when combined with the MARTINI "Lego" approach, can provide accurate descriptions of a system at a much lower computational cost.As in the case of the MARTINI approach, the intermolecular interactions between CG beads are obtained via a top-down approach based on the water-octanol partition coefficients [36].In general, MDPD is suitable for simulating multi-phase and multicomponent systems that are much larger than in MD [33,37,38].Moreover, in view of the rapid development of MDPD with applications in soft matter and fluids physics, we anticipate that our study opens new possibilities for this method in these areas by taking advantage of the MARTINI force-field approach.
In the next section, our methods and system descriptions are provided.Then, we juxtapose the results of the MD MARTINI with the MDPD MARTINI model for a lipid system, demonstrating the reliability of the MDPD MARTINI model with an accompanied lower computational cost.Then, in the last section, we draw our conclusions and outline the prospects for further development.
2 Theoretical Methods

Molecular Dynamics
We first simulate the formation of a DPPC lipid bilayer in water with MD simulation.This system has been used as a benchmark in the MARTINI online tutorials [39] by using GROMACS software [4].Here, we have used LAMMPS [2,3] software to carry out both the MD and the MDPD simulations.Regarding the MD simulations, initial files for the simulation of the DPPC bilayer can be found as a moltemplate example [40].
In particular, we initiated simulations for systems comprising 128, 256, and 512 lipids.These lipid-containing systems were initially placed within a simulation box characterized by dimensions L x = L y = L z = 70, 90, and 120, respectively, in conjunction with a water component maintaining a 10:1 proportion of water beads.The MARTINI DPPC model, employed in our simulations, encompasses a representation consisting of distinct beads: Q 0 to signify the choline headgroup, Q a beads representing the phosphate headgroup, two N a beads embodying the glycerol ester moiety, and two chains, each comprised of four C 1 beads, to emulate the alkane tails.Water entities were simulated using P 4 beads.
To initiate our simulations, we adapted the example files provided by moltemplate [40].We commenced the simulations by initially minimizing the system's energy, followed by conducting a simulation in the NPT ensemble employing the Nosé-Hoover thermostat at 300K and 1 bar.This step aimed to relax the system and attain an equilibrium box volume.Subsequently, we executed simulations under NVT conditions, elevating the temperature to 450K to facilitate the removal of any undesired structures that may have formed, paving the way for the bilayer formation.Further refinement of the system was achieved by running simulations within the NPT ensemble, employing a zero surface tension condition while gradually lowering the temperature to 300 K.This phase culminated in the attainment of the final bilayer configuration, characterized by the absence of defects.Lastly, we kept the NPT run at 300K to systematically measure key bilayer properties.

Many-body dissipative particle dynamics
The MDPD method has been applied for fluids with different properties [31-35, 41, 42], such as density and surface tension [37], as well as a range of soft matter multi-component systems, such as systems with surfactant molecules [43].In the case of the MDPD model, the Langevin equation (Equation 1) of motion is solved as is done in MD Langevin dynamics [44], but forces are directly defined in the MDPD method rather than derived from a potential form as in the case of MD.Equation 1 describes the motion of each particle, i, interacting with its neighbors, j, through a conservative force F C , namely, The random force, F R , and the dissipative force, F D , act on each particle, i, and function as system's thermostat, since they are related via the fluctuationdissipation theorem.The mass m is kept the same for all particles and set to unity.A main difference between the MDPD and the DPD [45][46][47][48][49] approaches is the expression of the conservative force, which also includes attractive interactions, and reads A < 0 and B > 0 are parameters of the force-field, which correspond to attractive and repulsive interactions, respectively.r ij is the distance between particles, e ij the unit vector connecting particles i and j, while ω C (r ij ) and ω d (r ij ) are linear weight functions defined as Here, r c is a cutoff distance for the interactions, which is usually set to unity.Moreover, ω d (r ij ) is of the same form, however, its cutoff distance r d = 0.75, smaller than r c .The many-body contributions to the force-field are included in the repulsive term and depends on the local neighborhood densities, ρi and ρj , via the following expression: The random and dissipative forces read with γ being the dissipative strength, ξ the strength of the random force, v ij the relative velocity between particles, and θ ij a random variable from a Gaussian distribution with unit variance.In this case, the fluctuationdissipation theorem require that γ and ξ are related to each other by The temperature of the system in our simulations is set to T = 1 (MDPD units) while the weight functions for the forces are Finally, the equations of motion (Equation 1) are integrated by using a modified velocity-Verlet algorithm as implemented in LAMMPS software [2,3] with a time step ∆t = 0.01.Describing the various interactions for the system amounts to determining the attractive parameter, A ij , of the conservative force between pair of particles, while the repulsive part of the interactions expressed through B will remain constant [50].To build molecules, one also needs to use bond and angle potentials.In particular, harmonic interactions, mathematically expressed as are used to bind consecutive beads in order to be able to form a macromolecular chain (e.g., lipid, surfactant).
Here, as is usually done in the case of CG models [51] and for the sake of simplicity, the same parameters are used for all bonds, namely k = 150 and r 0 = 0.5.Moreover, a harmonic angle potential is applied for a triad of sequential beads along the molecular chains, where k A = 5 and θ 0 , which depends on the specific molecule.In the following section, we discuss the forcefield parametrization of the intermolecular interactions.

Force-field parametrization
To obtain the water-octanol partition coefficients, we have measured the molar concentrations of the solute molecules in the water, [S] w , and 1-octanol phases, [S] o , following a previous study [36].Then, the water-octanol partition coefficient is defined as The advantage of this approach based on the wateroctanol partition coefficient lies in the fact that experimental data are abundant for various solute molecules [52,53].Moreover, the coefficient can be readily calculated in particle-based simulations by directly using the above equation and determining the quantities  11).
[S] w and [S] o , or via umbrella sampling simulations, as has previously been done in the case of the MAR-TINI model [13,15].However, umbrella sampling simulations are not expected to provide any advantage in the case of MDPD simulations, since MDPD uses soft interactions that provide better sampling and a more efficient path to equilibrium.Moreover, umbrella simulations would require an expression for the underlying potential, which may further complicate things.Here, a limitation in measuring log P with the current approach comes from the requirement that at least one bead is present in each phase.In this case, the maximum log P value that can be obtained is log N , when the number of solute beads is N .Our system contains 1500 solute beads (around 3% of the concentration), which implies a limit log P o/w = ±3.17able to describe both apolar and polar solutes.This is in general satisfactory, since the usual range in experiments is between −4 and 7.
Here, we have chosen to parametrize a composite system consisting of water and DPPC lipids [54], which serves as a benchmark in the simulation of lipid membranes.According to the MARTINI model, four heavy atoms are usually represented by one CG bead.While we follow here the MARTINI recipe, a different CG level might in general be chosen, for example as in the case of models based on the Statistical Associating Fluid Theory [51].The parametrization process outlined in this study involves several crucial steps to characterize and define the interactions between different types of beads in a complex molecular system.Firstly, the system is represented by three primary bead types in the reference system, namely the water (W) bead, the octanol (C) bead used for the alkane chain, and the G bead representing the polar alcohol portion of the chain.Notably, the W bead corresponds to three water molecules, while the C bead is employed to model a four carbon apolar chain, and the G bead captures the polar components.
Secondly, the self-interaction parameters, denoted as A ii , are determined for each bead type.These parameters are derived from surface tension measurements and are specific to different bead types.For instance, A W W pertains to water, A CC corresponds to octane (comprised of two bonded C beads), and A GG is associated with Diethylene Glycol [43] (composed of two bonded G beads).The surface tension can be measured by placing a slab of the desired component in a large simulation box and as usual computed from the pressure tensor (mechanical approach [55]) as where P zz is the pressure in the direction normal to the surface.The cross-interaction parameter between C and G beads, specifically A CG , is also calculated based on surface tension measurements involving octanol.Thirdly, additional cross-interaction parameters, namely A W C and A W G , are determined through a comparison with the partition coefficient (log P o/w ) for water beads in the octanol phase and octanol molecules in the water phase.This comparison aligns with experimental data, with water in octanol exhibiting a log(P ) of 1.3 (experimental) and 1.4 (simulated), while octanol in water shows a log(P ) of 3.1 (experimental) and 2.9 (simulated).The accuracy of the model is further validated by comparing the calculated Diethylene Glycol log(P ) values (-1.5), which closely match the experimental value -1.4.
Lastly, the parametrization process extends to other bead types by categorizing them as either polar or apolar based on their behavior and measuring their partition coefficients (log(P )) to capture their respective characteristics.An illustrative example is provided in Figure 1, where the H bead (solute) in the lipid model is identified as highly polar with a log(P ) value of -1.7.All the attractive parameters are presented in the interaction matrix of Table 1.

Results and Discussion
We initiated MDPD and MD MARTINI simulations by randomly distributing lipid molecules and water beads (Figure 2) within the simulation box.The system was  let to evolve until a defect-free bilayer structure emerged.
Remarkably, the behavior of MD MARTINI lipids in our simulations closely mirrored results that have been previously reported [11].In Figure 2, we present the results of an MDPD simulation featuring 512 lipids.Initially, the lipids exhibited a tendency to aggregate into small clusters, which subsequently coalesced to form a more extensive molecular assembly.Over time, this structure underwent further transformations, ultimately giving rise to a bilayer configuration.Throughout these dynamic transitions, we observed the formation of transient defects, such as water pores and lipid bridges, analogous to the observations made in MD MARTINI simulations.These defects vanished over time, leading to the attainment of the final defect-free bilayer structure.For comparison, Table 2 presents a time-based assessment of MDPD and MD MARTINI simulations for various system sizes under consideration.Notably, our findings indicate that MDPD is capable of generating defect-free bilayers akin to those produced by MD MARTINI simulations.However, a striking acceleration in computational efficiency was achieved through the utilization of soft, short-range interactions in MDPD.
System size 128 256 512 MDPD 20s 1.5min 5min MD 70s 5min 35min Table 2 Time comparison of simulation runs on eight cores Intel Core i7-10700 @ 2.9Ghz for the formation of a bilayer with no defects.Successful formation highly depends on initial configuration but MDPD is always four to seven times faster than MD.Three system sizes were considered, with 128, 256 and 512 lipid molecules; all with a 10:1 proportion of water beads for both MDPD and MD.
The area per lipid (APL) quantifies the surface area occupied by an individual lipid molecule.This measure is subject to influence by various factors, including the lipid type, temperature, and the presence of other molecules within the system.In our investigations, we adopted a method wherein the area of plane aligned tangentially to the bilayer surface in the simulation box, was divided by half the total number of lipid molecules in order to determine the APL.On the one hand, for the MDPD DPPC model, the resulting APL was found to be 47.9 Å2 , a value that aligns closely with experimental data and recent MD studies [56].In the MD MARTINI model, on the other hand, the calculated APL was 61.9 Å2 , a measure comparable to that of a DPPC bilayer at a temperature of 323K [57].
Another fundamental property we assessed was the bilayer thickness, a parameter essential for characterizing the structural aspects of lipid bilayers.This property was determined by identifying the peaks in the density distribution profile of the choline headgroup bead (H1 in the MDPD model and Q 0 in the MAR-TINI model), as illustrated in Figure 3.The thickness of the MDPD bilayer was computed to be 45 Å, while that of the MD bilayer was 42 Å.Notably, both values closely mirror experimental measurements [58], affirming the accuracy and reliability of our simulation results in capturing the essential structural characteristics of DPPC lipid bilayers.
Lipid bilayers, possessing elastic properties represented by the bending modulus, k c , play a crucial role in membrane dynamics.The bending modulus quantifies the intensity of undulations within the lipid membrane and can be derived from the fluctuation spectrum [59].In our analysis, equilibrated bilayer systems, each comprising 512 lipids, are replicated four times in both x and y directions, resulting in a large membrane consisting of 8192 lipids.Undulations are described as a surface function u(x, y), representing the average position of two monolayers defined by the positions of phosphate head beads.The undulation spectrum ûis obtained by performing a Fourier transform on the surface function using the expression: From a continuum model of the bilayer, the undulation intensities in Fourier space are given by where A = L x L y is the lipid membrane area, q = |q| = q 2 x + q 2 y the wavevector length in the radial direction Fig. 4 Spectrum intensities of the undulation modes on the bilayer.Long wavelengths obey the q −4 behavior and transition around a 1 nm wavelength to either a constant value in the MD simulations or to a q −2 regime in MDPD simulation.The data is averaged for 100 time steps of a system consisting of 8192 lipids.
and γ is surface tension.The analysis of the undulation spectrum, as illustrated in Figure 4, reveals that both MD and MDPD simulations exhibit the expected q −4 scaling for the long wavelength modes.The dashed line represents a reference from the experimental bending modulus value of k c = 5 × 10 −20 J [60].For wavelengths around 1 nm and below, a transition occurs, with MD simulations showing a constant intensity and MDPD scaling exhibiting a q −2 regime.This deviation may stem from differing surface tension conditions between simulations, with MDPD performed on an NVT ensemble and MD on an NPT ensemble employing a Nosé-Hoover barostat to maintain a zero surface tension condition.The q −2 line shown is the figure is obtained

MDPD MD
Fig. 5 Second rank order parameter calculated for every bond in the lipid molecule.The data is averaged for 100 time steps and over all the bonds in the bilayer consisting of 8192 lipids.
by considering γ = 50 mN/m.The difference could also be due to a possible protrusion regime present in the MDPD model or tilt angle of the lipid molecules as both can give rise to a q −2 dominant regime [61,62].
To further compare the MD and MDPD methods, we compute the second-rank order parameter P 2 defined in Eq. 15 [63,64].This is an important parameter that allows us to identify the alignment of the lipid molecules with respect to the bilayer by measuring the angle θ between each bond and the bilayer normal direction.
Figure 5 shows the P 2 value for each bond in the lipid molecule.A value of P 2 = 1 indicates perfect alignment, P 2 = 0 a random orientation and P 2 = −0.5 antialignment.Our measurements suggest that the MDPD method results in lipids more aligned to the bilayer normal direction, possibly influenced by the presence of surface tension in the MDPD simulation [54].Both methods exhibit a similar trend along the molecule chain, highlighting the robustness of the observed alignment characteristics.

Conclusions
We have examined the efficacy of our proposed MDPD force-field in comparison with the widely recognized MARTINI force-field for CG MD simulations.Both methodologies demonstrated the capability to produce defectfree lipid bilayer structures, accompanied by similar self-aggregation dynamics.Furthermore, our investigation delved into the assessment of fundamental structural attributes of these lipid bilayers, specifically the determination of area per lipid, bilayer thickness, bending modulus, and orientation order parameter.Remarkably, both MDPD and MD MARTINI simulations yielded results that closely approximated experimental observations.A noteworthy advantage of the MDPD method lies in its utilization of soft, short-range interactions, a feature that enables the realisation of faster simulations without comprising on key properties of the system selected here for demonstration purposes.This enhancement in computational efficiency holds paramount significance in the context of large-scale molecular simulations, potentially facilitating the exploration of increasingly intricate and large molecular systems.
Not only our study highlights the effectiveness of the MARTINI "Lego" approach when applied to the MDPD method, allowing for the representation of complex molecular interactions in a CG fashion through bead parameterization based on hydrophobic/hydrophilic characteristics, but also underscores the versatile nature of the MDPD force-field.This adaptability paves the way for the incorporation of additional interaction levels, diverse bead types, and varying bead sizes, thereby extending its applicability to a wider spectrum of molecules and systems.In conclusion, the MDPD approach emerges as a versatile and computationally efficient tool for the simulation of large-scale molecular systems, offering an intriguing alternative to the MD MARTINI framework and holding promise for further advancements in CG simulations.

Fig. 1
Fig.1(a) Typical simulation snapshot of a water-octanol system used for the parametrization of a solute particle (H bead) based on the o/w partitioning (Equation11).

Fig. 2
Fig. 2 MDPD CG model for DPPC lipid and water molecule.Snapshots illustrate the gradual formation of the DPPC bilayer starting from a random distribution of 512 lipid molecules that aggregate into clusters which then coalesce to produce the lipid bilayer.Time τ is in MDPD reduced units.

Fig. 3
Fig. 3 Comparison of Bilayer Thickness in MDPD and MD Simulations.(a) Density profile from MDPD simulation, with a bilayer thickness measured between the peaks of the top H bead (H1) equal to 45 Å.(b) Martini simulation results showing a bilayer thickness of 42 Å.Both simulations consist of 512 lipid molecules.