Dual Resolution Membrane Simulations Using Virtual Sites.

All-atomistic (AA) and coarse-grain (CG) simulations have been successfully applied to investigate a broad range of biomolecular processes. However, the accessible time and length scales of AA simulation are limited and the specific molecular details of CG simulation are simplified. Here, we propose a virtual site (VS) based hybrid scheme that can concurrently couple AA and CG resolutions in a single membrane simulation, mitigating the shortcomings of either representation. With some adjustments to make the AA and CG force fields compatible, we demonstrate that lipid bilayer properties are well kept in our hybrid approach. Our VS hybrid method was also applied to simulate a small lipid vesicle, with the inner leaflet and interior solvent represented in AA, and the outer leaflet together with exterior solvent at the CG level. Our multiscale method opens the way to investigate biomembrane properties at increased computational efficiency, in particular applications involving large solvent filled regions.

Introduction

Molecular dynamics simulations have been successfully applied to investigate configurations and dynamics of biomolecular processes for many years.[1−3] The accessible time and length scales are, however, limited for all-atom (AA) simulations. The computational efficiency has been accelerated by orders of magnitude by combining atoms into single effective interaction sites in coarse-grain (CG) simulation.[4] However, the simplifications inherent to a CG model, i.e., the loss of some atomic detail comes at the price of a reduced accuracy. This limitation can be overcome by combining AA and CG resolution in a multiscale model. Several multiscale schemes have been put forward, falling into four classes: sequential, adaptive, mixed, and resolution exchange methods. In the sequential approach, one starts the system from the CG resolution and some key states or interesting configurations found in this representation are transformed to the AA resolution.[5,6] In this way, more efficient sampling is achieved at more coarse representations, while more detailed information is obtained at less coarse representations. The adaptive approach allows atoms or particles to adapt their resolutions on the fly and transfer freely between spatially localized resolution domains without any barrier.[7,8] The mixed approach enables molecules with AA/CG resolutions to coexist in the same system, interacting with each other, while the resolution of the particles cannot be changed, analogous to the well-known QM/MM approach that partitions a system in nuclei and electrons, on the one hand, described at the quantum-chemical level, and an embedding described at a force field level.[9−11] The resolution exchange approach simulates a batch of replicas with different resolutions in parallel and the configurations can exchange between replicas when the detailed balance condition is satisfied, e.g., Hamiltonian replica exchange.[12−14] Further comprehensive understanding of multiscale modeling approaches can be found in a number of reviews.[15−17]

The different multiscaling strategies are associated with different technical and theoretical challenges as well as different computational costs. For the sequential multiscale method, the determination of when to change resolution (frequency and choice of states) is not trivial and the change of resolution can move a conformation away from its initial state before resolution exchange during relaxation under the new force field (especially from CG to AA). For the adaptive resolution method, the simulation is expensive and not compatible with any standard molecular dynamics package (e.g., Gromacs). This is because an added thermodynamic force[18] is required to compensate for the difference in chemical potentials and equations of state at the AA and CG resolutions. This force must be computed first and then added as an external force. For the mixed method, the accuracy of the model in the resolution interface cannot be guaranteed, especially for AA/CG hybrid simulation.[19] Furthermore, the types of potential energy functions (e.g., Coulombic, Lennard-Jones) may not be the same for both resolutions, especially in force-matching type potentials, requiring the use of tabulated potentials.[20,21] For the resolution exchange method, the simulation is computationally demanding, e.g., in Hamiltonian replica exchange, several replicas are needed to simulate the same system in parallel while only one of them can be used in the final analysis.[12]

In this paper, the aim is to combine molecules described in a particle-based manner at two different resolutions in a single simulation using a single set of standard potentials, in particular atomistic (AA) and coarse-grain (CG), but the method should be applicable with other particle-based combinations as long as the potentials are available. Such an approach is expected to be technically most straightforward and computationally most efficient, because it uses only the core functionalities of existing molecular modeling code. Here, we demonstrate this method combining the Martini CG force field[22] with the atomistic GROMOS[23] and CHARMM[24] force fields, employing the Gromacs code,[25] and apply the method to lipid systems. The method builds on the experience of embedding of atomistic particles in a CG environment with the help of virtual sites (VS), described earlier by Rzepiela et al.[10] The main idea of this VS hybrid scheme is that AA particles interact with AA particles according to the AA force field, and CG particles interact with CG particles according to the CG force field. AA particles do not interact directly with CG particles, but the interaction between molecules at AA resolution with molecules at CG resolution is achieved through VS. These VS carry interactions at the CG level; i.e., the AA molecule “sees” the embedding surroundings at the level of the CG model. Conversely, the CG molecules “see” the AA molecules as CG molecules. A complication arises when dealing with the electrostatic interaction between AA and CG subsystems. The AA–CG electrostatic interaction is either ignored or it is acting at the level of the charged particles, be they part of the AA or CG subsystem. Both approaches are possible with straight cutoff or smoothed electrostatic potentials (shift[26] or reaction field[27]), but the latter is the only option when using full electrostatics (such as particle mesh Ewald (PME)[28]). Even in this way, the polar (charged) particles of different resolutions interacting in close proximity to each other still cause problems, as shown by Wassenaar et al.[19] In their simulations, VS hybrid solutes, e.g., amino acids, were solvated in CG solvent. Compared to fully atomistic results, the hybrid AA/CG model was successful for apolar molecules. However, their procedure could not reproduce correct potentials of mean force (PMFs) between pairs of amino acid side chain analogues in water and partitioning free enthalpies of amino acid side chain analogues of charged and polar molecules, despite adjustments in the relative dielectric permittivity to couple the AA–CG electrostatic interactions.[19] Therefore, to avoid polar resolution interfaces, we propose a dual resolution membrane structure with the AA–CG interface at the interleaflet region, constituting an apolar environment.

In this paper, the VS hybrid scheme is introduced first, followed by a number of validations of embedded lipid systems by comparing their conformational and dynamic behavior against corresponding reference systems. Finally, some applications are proposed and the computational efficiency is compared against the corresponding pure AA simulations.

Methods

VS Hybrid Model

The virtual site (VS) hybrid scheme employed here is based on the combination of Martini and GROMOS force fields described by Rzepiela et al.,[10] but the idea can be applied to other combinations as well. The system is partitioned into two subsystems, i.e., all-atom (AA) and coarse-grain (CG). The crux of the method is that each subsystem interacts with itself at its corresponding level of description and that the interaction between the subsystems is at the CG level, thus avoiding the need to define and refine direct interactions between two models at different resolutions and with possibly quite different interaction potentials. This approach requires a mapping of the particles in the AA model to the CG beads, defining interaction sites through which the AA system interacts with the CG system. This hybrid technique is achieved through virtual sites located at the center of mass of the corresponding AA atoms:M and R are the total mass of the particles assigned to VS bead k and the position of the VS bead, respectively. The kth VS bead is constructed from N AA atoms with mass m and position r for the ith atom of the kth VS bead. Depending on the mass of the AA particles, the force acting on the VS is distributed over its corresponding AA atoms asf and F are forces acting on AA atom k belonging to bead k and on VS bead k. The VS hybrid scheme naturally combines the CG and AA models, where the AA model is a more detailed description of the molecule of interest. CG particles feel forces only at the CG level, originating either at particles in the CG subsystem or at the VS. Direct interactions between AA atoms and CG beads are not included, because these would constitute a double counting of interactions. The interaction is schematically shown in Figure for lipid molecules.

Figure 1

Multiscale system encompasses both AA/CG molecule representations. Each subsystem (AA and CG, respectively) interacts within itself according to its own resolution. Molecules in the AA region interact at CG level with the molecules in the CG region through virtual sites that represent the AA region. When employing full electrostatics (PME), long-range electrostatic interactions between AA and CG charges are included (green dashed line). The AA region is shown as bonds, the virtual sites representing the AA system as well as the CG region is shown as spheres. Lipid tail, linker, and head are represented by cyan, green, and orange, respectively.

Matching AA and CG Force Fields

When comparing different force fields at AA and CG resolutions, they usually differ in choices made in the interaction functions and/or simulation settings. For example, the Martini and the GROMOS force fields use cutoffs of 1.1 and 1.4 nm, respectively. Also, the shape of the nonbonded potentials is slightly different in these force fields, because they are either smoothed (and/or shifted) or straight cutoffs are applied. A simulation engine such as Gromacs allows only one choice for the potential form and cutoff. One way to deal with the different cutoffs and potential shapes in a hybrid simulation is to use tabulated potentials.[10] However, the tabulated potential can largely limit the computational efficiency and does not support GPU acceleration when using Gromacs.[29] In order to avoid having to use a tabulated Lennard-Jones (LJ) potential, a cutoff of 1.4 nm was used in the hybrid model to guarantee the accuracy of the AA components in the model. The change in cutoff for the CG force field means that CG interactions are compromised, and a solution is proposed in the Results of this paper. Another issue is the interaction between charges. Martini and GROMOS use a different dielectric constant because the screening in the AA force field is done by explicit solvent charges (SPC water for example), whereas standard Martini solvent does not have partial charges and screening must be made explicit by using a dielectric constant not equal to 1 in the Coulombic potential. In order to avoid having to use a Coulombic tabulated potential, the charge of CG beads was therefore scaled down with a factor of 0.258, which is computed by , where ϵrCG = 15 represents the screening factor in standard Martini setup. Thus, we can apply a uniform dielectric screening factor (ϵr = 1) for the whole system. Full details regarding the force fields can be found in the Supporting Information.

Considerations for the Gromacs Implementation

In our implementation, we use the standard Gromacs 2016 simulation package.[30] In technical terms, the absence of direct interactions between AA and VS, VS and VS, and AA and CG particles may be achieved through setting nonbonded interaction parameters (in particular Lennard-Jones parameters) between these two subsystems to zero in the Verlet update scheme[31] or by explicitly excluding the interactions between the AA and CG, AA and VS, and VS and VS particle groups in the Group cutoff scheme. VS beads do not carry charge. If the system of interest is such that AA and CG charged components are so distant that they are beyond the cutoff for Coulombic interactions, the former approach is suggested, since it has better computational efficiency. The latter approach is more general and functional for nonbonded short-range potentials, but it does lead to a complication regarding the long-range electrostatic interactions when a full electrostatics method, such as PME[28] is used, since PME automatically considers all the long-range interactions between charges in the system.

VS Hybrid Scheme Applied to a Lipid Membrane

An intuitive way to build a dual resolution membrane system is a combination of AA membrane and CG water, in the same way as the earlier model for soluble protein,[19] and as illustrated in Figure S1a. The AA dipalmitoylphosphatidylcholine (DPPC) membrane is embedded in solvent water at the level of the CG model. However, this setup leads to major artifacts: compared to the standard AA simulations (reported in Table ), the hybrid system has a much smaller area per lipid (0.521 ± 0.004 nm2) and bigger membrane thickness (5.718 ± 0.031 nm). The lipid head groups aggregate too much and apolar acyl chains are exposed to the CG water, as shown in Figure S1b. These phenomena are likely due to the lack of screening for AA Coulombic interactions in the lipid head components from the CG solvent (without partial charge). In previous works,[32,33] this problem was solved for solutes and proteins by restraining a layer of AA water around the AA solute solvated in the CG water environment. However, it is hard to guarantee the correct thermodynamic and dielectric properties of water close to the resolution interface, which can still affect the solvation of the solute. It was also shown that the VS hybrid GROMOS/Martini model can properly reproduce the potentials of mean force (PMFs) between pairs of apolar amino acid side chain analogues, yet failed to reproduce correct PMFs for the polar and charged AA solutes in CG solvent.[19] Therefore, in this work, the VS hybrid membrane systems are constructed as shown in Figure , to keep the resolution interface at the apolar lipid tail region.

Table 1

Membrane Properties of DPPC Membrane at 323 K

resolution	setting	diffusion constant (10–6 cm2/s)	area per lipid (nm2)	membrane thickness (nm)
AA	standard	0.033 ± 0.002a	0.611 ± 0.004b	3.92 ± 0.03b
AA	hybrid	0.032 ± 0.003	0.626 ± 0.012	3.82 ± 0.05
CG	standard	0.85 ± 0.06	0.627 ± 0.010	4.08 ± 0.04
CG	hybrid	0.25 ± 0.02	0.572 ± 0.004	4.20 ± 0.02
CG	hybrid with scaled ε	0.61 ± 0.06	0.627 ± 0.003	4.04 ± 0.03
VS hybrid	hybrid with scaled ε	AA: 0.052 ± 0.004	0.628 ± 0.003	3.89 ± 0.03
		CG: 0.40 ± 0.04

The standard error of diffusion constants considering each lipid independently (Figure S3).

Standard deviation obtained after the value reached equilibrium. Standard errors are <0.1%.

Figure 2

Virtual site hybrid membrane setups: (a) planar DPPC membrane; (b) DPPC vesicle. Lipids in the AA region carry virtual sites. Lipids tail, linker, and head parts are represented by cyan, green, and orange, respectively.

Simulation Setup

There are several GROMOS models for the DPPC lipids. Poger lipids[34] developed with the GROMOS 53a6 force field can only reproduce the reported membrane properties (i.e., area per lipid value of 0.63 nm2) with Gromacs versions 4.0.X and earlier.[35] The Gromacs developers were contacted and details regarding this problem can be found at http://redmine.gromacs.org/issues/1400. Thus, we used a modified version of the GROMOS 53a6 force field and lipid topology, which are more compatible with the newer version of Gromacs 2016 software. The corresponding files are included in the Supporting Information. We simulated planar membranes of 162 lipids (81 in each leaflet), as illustrated in Figure a. Water molecules could permeate through the biological membrane, so flat-bottom position restraints were added to AA water molecules to prevent them from permeating into other resolution regions (or crossing the resolution boundary), as illustrated in Figure S2. A force constant of 100 kJ/mol–1 was used and the distance between resolution interface and the beginning of the flat bottom potential is 0.5 nm. It is reported that, in the Martini model, the magnitude of the unidirectional flux is 1 CG water (4 water molecules) per 100 ns for a 256 DPPC bilayer patch.[36] Considering the time scale of our simulation, a flat bottom constraint for CG water is not necessary. However, this constraint on CG water is suggested for longer simulations. In our simulations, even though the unrestrained CG water may permeate through the membrane and go to the AA region, the CG water cannot interact with the AA water and acts as water vapor.

We have also applied the VS hybrid scheme to a vesicle with a diameter of 20 nm, as illustrated in Figure b. We chose 20 nm because it is the minimal size of DPPC vesicles observed in experiment.[37] To build the VS dual resolution vesicle, we first built a CG vesicle using the CHARMM-GUI Martini maker.[38] The number of lipids within each leaflet of the vesicle was estimated on the basis of the area per lipid. During the equilibration processes, six artificial pores in the vesicle were kept open with constraints to enhance lipid flip-flop and to equilibrate the number of water inside and outside the vesicle.[39] After equilibration, the vesicle was closed by healing the water pores by switching off the restraints. CG water inside the vesicle and the inner leaflets were then mapped into AA resolution with the backmapping software tool backward.[40] The final configuration contained 996 lipids in the inner leaflet and 1748 lipids in the outer leaflet. Similar to the planar membrane, spherical flat-bottomed position restraints were added to AA water in the production run (Figure S2).

We have used an integration time step of 2 fs in standard AA (both GROMOS and CHARMM) simulations and CHARMM/Martini VS hybrid simulation, and 20 fs in standard Martini simulations. However, we have used 4 fs for GROMOS/Martini VS hybrid simulations and GROMOS related simulations in hybrid settings, because the conformational and dynamical properties of united atom DPPC lipids are correctly reproduced in simulations with an integration time step up to 5 fs.[41] The Verlet cutoff scheme was used for the system, and the allowed energy error due to the Verlet buffer is 0.005 (kJ/mol)/ps per atom. The potential modifiers were used to shift the complete LJ and Coulombic potential value to zero at the cutoff. Electrostatic interactions were treated using PME with a 0.12 nm Fourier grid spacing. The temperature was maintained at 323 K by integrating the equations of motion with the leapfrog stochastic dynamics integrator[42] with an inverse friction constant of 1 ps to CG (or VS) particles and 0.1 ps to AA particles. Noted that even though an inverse friction constant is given to the temperature coupling of VS components, VS particles cannot really be a part of the thermostat, since they have no mass. AA lipids, AA water, CG lipids, CG water, and VS beads were coupled separately to a thermostat. The pressure was kept at 1 bar independently in the lateral and normal directions with the Berendsen barostat[43] in a semi-isotropic pressure bath (τp = 0.5 ps) in the planar membrane, while an isotropic pressure bath was used in the vesicle. The compressibility was 4.6 × 10–6 bar–1. The AA bonds were constrained with the LINCS algorithm, and simple point charge (SPC) water[44] was constrained using SETTLE.[45] We refer to this setup as the hybrid setting. It is applied to a DPPC lipid membrane and validated against standard AA and CG simulations in the Results. Note that, compared with this hybrid setting, the only difference is that the GROMOS force field in a standard setting uses a time step of 2 fs and applies the reaction field to calculate the long-range Coulombic interactions.

Free Energy Calculations

Alchemical transformations were used to compute free energies of solvation ΔGsolv, hydration ΔGhydr, and partitioning ΔGpart. The hydration and solvation free energies of target bead types were calculated by the thermodynamic integration (TI) method,[46] decoupling the solute bead from its surrounding solvent (water or hexadecane) molecules by turning all solute–solvent interactions off. The partitioning free energy was obtained from ΔGpart = ΔGhydr – ΔGsolv. A series of 11 simulations with equally spaced λ points going from 0 to 1 were performed at 300 K for each selected CG bead and each λ point lasted for 10 ns. A soft core potential[47] was used to avoid singularities due to solute–solvent particle overlaps as interactions were switched off, where α = 0.5, σ = 0.47, and the soft-core λ power is set to 1. The derivative of the free energy with respect to λ was integrated through the trapezoidal method. Error estimations were obtained using the Multistate Bennett Acceptance Ratio.[48]

Analysis

The behavior of the lipid assemblies were analyzed at an aggregate and molecular level. Aggregate level analysis included calculation of the projected area per lipid, bilayer thickness, temperature, etc. At the molecular level, second rank order parameters were calculated according toWhen the order parameter is computed at the CG level, θ represents the angle between the CG bond of the lipid and the bilayer normal in planar membrane or between the CG bond and the line connecting the center of mass of vesicle to the center of the CG bond in the vesicle. When the deuterium order parameter is computed in the AA level, θ represents the angle between the orientation of the C–H bond vector with respect to the bilayer normal in planar membrane or between the C–H bond vector with respect to the vector from the center of mass of the vesicle to the target carbon atom. Angular brackets denote the average over the ensemble of bonds and over time.

For a comparison between hybrid and CG models, atoms in the AA simulation are first mapped into CG beads in the analysis processes and then the analysis is performed at the CG level. Membrane thickness and area per lipid were computed using software called FATSLIM.[49] PO4 beads of the PC lipids were used as references to compute membrane thickness between two leaflets. Since the fluctuations or curvature of the membrane may introduce noise, the thickness and area per lipid are computed on the basis of neighborhood-averaged coordinates to smooth the fluctuations. The specific explanation and software are freely available on Web site http://fatslim.github.io/.

Partial density distributions were obtained from a simulation of a bilayer, where two leaflets of membrane and water are considered separately. The radii of the vesicle were computed by the average distance between the PO4 beads, either in the inner leaflet or in the outer leaflet, and the center of mass of the vesicle.

The diffusion constant was obtained by fitting the linear part of the mean square displacement (MSD) of lipids in the lateral (xy) plane. The error estimation is obtained from the standard error of the distribution of diffusion constants of each individual lipid.

Results

Consequences of Using Single Set of Parameters for Two Different Resolution Models

The essence of the VS hybrid model is that the highest computational efficiency with the Gromacs simulation engine can be achieved when atomistic (AA) and coarse-grained (CG) resolution models are treated with a single set of standard force field functions and run settings. Our proposed hybrid setup in general therefore may necessitate that the AA or CG model or both models are run with different parameters (functional form of the interactions, cutoff, electrostatics scheme, time step) from the ones they were developed with. To gain insight into the effect of changing these input parameters to the ones used in the hybrid setting described in the Methods, we simulated planar DPPC membranes in the AA and CG models and compared the conformational and dynamic properties of the membranes to their corresponding references in the standard settings. Note that the only difference between hybrid setting and standard setting for AA simulation is the way the long-range Coulombic interaction is calculated (PME versus RF) and the time step (4 fs versus 2 fs). It can be seen in Table and Figure that, compared to the AA model in standard setting, the AA model in a hybrid setting has a slightly lower membrane thickness and deuterium order parameter and a slightly higher area per lipid. The position of the water/lipid interface and the extent of hydrophobic tail interdigitation between two leaflets were estimated from the partial density profiles across the membrane, shown in Figure e,f. The AA water/lipid interfaces and the interdigitation in the hybrid setting agree well with their corresponding references in standard setting. The lateral diffusion constants of the lipids in the AA models in the two settings are the same. MSD for individual lipids are computed and shown in Figure S3. The distributions of the MSD are also close for simulations in hybrid and standard setting. Therefore, the membrane properties of the AA simulation in a hybrid setting are, even though not the same, very close to their AA references and also agree well with the standard CG simulation after mapping to the CG model, as shown in Table .

Figure 3

Membrane properties in different settings. (a) and (b) represent the order parameter in planar membrane and (c) represents the order parameter of outer leaflet in vesicle. (a) represents the deuterium order parameter (−SCD) of two DPPC tails. (b) and (c) represent the CG order parameter, which are computed on the basis of the bonds connecting CG beads. (d) and (e) represent the partial density of membranes in CG and AA resolution, respectively. (f) represents the number density of VS and CG beads in VS hybrid membrane. The scaled Martini potential was applied in the hybrid setting and VS hybrid model.

CG simulations in a hybrid setting were also compared to those in the standard CG setting. It can be seen in Table that the CG model in hybrid setting has a lower area per lipid, a lower diffusion constant, and a higher membrane thickness than in the standard setting and the agreements are less ideal than the correspondence between the AA model in the two settings. To further investigate the cause of the disagreement, we simulated several CG systems in a hybrid setting with nonbonded cutoffs ranging from 1.1 nm (the standard cutoff for the CG model) to 1.4 nm (the standard cutoff for the AA model and the chosen cutoff for the hybrid setting). The rest of the input parameters are the same as in the hybrid setting. In Table S1, we report that membrane thickness and order parameter increase with the cutoff, while area per lipid and diffusion constant show an opposite trend. This is because a longer cutoff has stronger nonbonded interactions, which can order the lipid arrangement. The slower diffusion rate can be either caused by more tightly packed lipid arrangements or by the stronger interactions due to the longer cutoff. To answer this question, we ran CG simulations at constant area, but with different LJ cutoffs. In Table S2, it can be seen that the diffusion constant is similar for simulations with the same cutoff, even though the area per lipid is different. Therefore, a different cutoff (or interactions) is the reason behind the different diffusion rates.

The strong decrease in area per lipid of the CG model in a hybrid setting leads to an undesirable mismatch of lipid bilayer properties between AA and CG models, which can lead to artifacts in the hybrid setup. We chose to fix the CG disagreement in a hybrid setting against the standard setting, by applying a scaling factor of 0.898 on ε in all the LJ potentials for the Martini subsystem in the hybrid setting. This scaling factor was arrived at by investigating pure hexadecane and water systems. In the newly parametrized Martini potential, the densities of hexadecane and water are similar to the simulations in a standard setting; see Table S3. In addition to achieving correct densities, a correct representation of the mutual solubilities of hexadecane and water phases is also important in the Martini force field since it is based on partitioning. Randomly mixed systems containing hexadecane and water molecules were prepared and were observed to quickly phase separate to form an aqueous slab and an oil slab. The density profiles of the two components were compared between Martini simulations in both settings and can be seen to agree well in Figure S4. The reproduction of correct free energy trends is the hallmark of the Martini force field.[50] With the newly scaled interaction table, the interactions of building blocks have been changed with respect to that of the standard Martini beads. Consequently, the free energy of solvation in different solvents of the CG particle types needs to be re-evaluated. In Figure , the results of the free energy calculations in a hybrid setting are presented and compared to standard Martini values for representative chemical building blocks (or Martini beads). The scaled CG model reproduces almost identical values for hydration free energy, hexadecane solvation free energy, and partition free energy when compared to the standard Martini CG model. Thus, the scaled LJ potential can ideally reproduce the Martini force field with a cutoff of 1.4 nm in hybrid setting. The Coulombic interaction cutoff does not play an important role in the Martini force field, since Martini systems are very sparsely charged and the Coulombic interaction beyond the cutoff can be taken into account in reaction field or PME approaches. We have computed properties of the zwitterionic lipids (DPPC and DLiPC) and negatively charged lipids (DPPG) using Martini simulations with the scaled potential in a hybrid setting. The long-range interactions were calculated using either reaction field or PME approaches. In Table S5, it can be seen that the structural properties, like area per lipid, membrane thickness, and order parameter, are very close to the standard Martini simulations. The order parameter profile and the partial density profile, which indicates water/lipid interface and interdigitation level between two leaflets, of the DPPC membrane are also similar to standard Martini simulations; see Figure b,d. In Table , it is shown that the diffusion constant is higher than the simulation in a hybrid setting without scaling and lower than the simulation in a standard setting, which agrees with our previous finding that the interactions rather than membrane structure are the main reason for the different diffusion rates.

Figure 4

Free energy comparison between Martini simulations in hybrid and standard settings. The data points in the plot represent C1, C4, Na, P4, Q0, and Qa beads, respectively. Hydration free energy, hexadecane solvation free energy, and partition free energy are compared in the plot and details can be found in Table S4.

Validation of the VS Hybrid Scheme for Lipid Bilayer Systems

If we apply the unscaled LJ parameters to the CG part of the VS hybrid approach on a membrane, a DPPC membrane forms a ripple-like phase, as illustrated in Figure S1c. This is caused by the too strong CG interactions with a cutoff of 1.4 nm. Therefore, in the rest of the paper, the scaled LJ potential is included in the hybrid setting.

The conformational properties of VS hybrid models were investigated for the double lipid bilayer system shown in Figure a. In Table and Figure , the values of area per lipid, membrane thickness, and order parameter (both AA and CG components) of the VS hybrid model are shown to agree well with their corresponding references in standard setting. As shown in Figure d,e, the water/lipid interface and interdigitation level between two leaflets estimated from partial density are reasonable compared to their corresponding references. Therefore, the conformational properties in the VS hybrid model are promising.

Thermodynamic and dynamic properties of the VS hybrid model are also tested and compared with pure AA simulations. Temperatures of different components in the VS hybrid model are under good control by the thermostat (Table S6), regardless of the hybrid structure. In addition, the diffusion constant of lipids in the AA leaflet in the hybrid model is slightly higher than that of the pure AA simulations in hybrid setting and the diffusion constant of the lipids in the CG leaflet in the hybrid model is slightly lower than that of the pure CG model in the hybrid setting (Table and Figure S3). This is because in the VS hybrid model, the AA lipids experience a lower friction due to the presence of a CG leaflet.

Applying the VS Hybrid Scheme to a Vesicle

The VS hybrid scheme was also applied to a vesicle as illustrated in Figure b. Since it is very expensive to simulate an equilibrated vesicle in AA simulation, the VS hybrid vesicle was mostly compared against CG simulation. Area per lipid and radius of both inner and outer leaflet of the vesicle, membrane thickness (Table ), and CG order parameter of outer leaflet (Figure c) in the VS hybrid model fit well with CG vesicle references. In general, the conformational properties of the VS hybrid vesicle are very reasonable.

Table 2

Membrane Properties of the DPPC Vesicle

resolution		area per lipid (nm2)	radius (nm)	membrane thickness (nm)
CG	inner leaflet	0.558 ± 0.001a	6.9 ± 0.9a	3.79 ± 0.01a
	outer leaflet	0.804 ± 0.001	10.6 ± 0.9
VS hybrid	inner leaflet	0.579 ± 0.001	6.8 ± 0.5	3.65 ± 0.01
	outer leaflet	0.781 ± 0.002	10.4 ± 0.5

Standard deviation obtained after the value reached equilibrium. Standard errors are <0.1%.

Note that the standard deviation of membrane thickness and area per lipid in the vesicle is smaller than those of the planar membrane (cf. Table ). This is because, in FATSLIM software, the noise introduced from the fluctuations or curvature of the membrane is smoothed out, while the fluctuations of the simulation box in xy directions in the planar membrane are included in the standard deviation.

Computational Efficiency

To assess the performance of the VS hybrid scheme, its computational efficiency was compared to a pure AA DPPC model with the same compositions. Table shows that the number of atoms or beads in the VS hybrid model is about 2 times lower than a pure AA model for the planar membrane model and about 6 times lower for the vesicle model. Thus, in the VS hybrid model, less computer memory is needed and fewer interaction pairs need to be taken into account. The computational speed is about 2 times faster for the planar membrane and about 5 times faster for the vesicle. The improved performance will be larger in a bigger system, when the size of the CG region is further increased.

Table 3

Performance for DPPC Systems Using Different Simulation Schemesa

system	resolution	number of particles	performance (ns/day)b	performance (step/day × 105)
vesicle	AA	2 402 906	0.14	0.67
vesicle	VS hybrid	394 631	0.66	1.7
planar membranec	AA	83 624	3.6	18
planar membranec	VS hybrid	46 663	8.6	23

AMD 2.3 GHz, 12 CPU (12 MPI processes).

The AA resolution used a time step of 2 fs and VS hybrid used 4 fs.

The planar membrane systems include two membranes as shown in Figure a.

Applying the VS Hybrid Scheme to the CHARMM and Martini Force Field Combination

We also apply the VS hybrid model to combine CHARMM36[24,51] with Martini force fields. The standard Martini interaction table was used in CHARMM/Martini VS hybrid model, since the CHARMM force field uses the same nonbonded cutoff as the Martini force field in the old setting[22] (1.2 nm). Even though the nonbonded interaction potentials are different between the two force fields (details in Supporting Information), the Martini model still reproduces the lipid membrane structure reasonably well in the CHARMM setting, as shown in Table . Therefore, the hybrid setting uses the standard AA setting to guarantee the accuracy in AA resolution. However, in the CHARMM/Martini VS hybrid model, the area per lipid in the VS hybrid model is too high and the membrane underwent a strong interdigitation between the two leaflets, shown in Figure S5a and Table . This is because CHARMM DPPC[52] has partial charges in the acyl tails that can interact with the Martini DPPC full charge head components without screening, which means the attractions between the two leaflets are overestimated. Note that this problem can be partially solved by excluding the short-range interactions between CG and AA components with the group cutoff scheme. However, this cutoff scheme sacrifices the computational efficiency. Therefore, we choose to decrease the LJ interactions between C1 and VC1 beads (interactions between tail parts of CG and AA leaflets) in the VS hybrid model from ε = 3.5 kJ/mol to ε = 3.25 kJ/mol. By doing so, the interaction between AA–AA and CG–CG subsystems are not affected, only the coupling (or attraction) between two resolutions is slightly decreased. The membrane structure properties (area per lipid, membrane thickness, partial density, and deuterium order parameter) are more reasonable, as shown in Table and Figure S5.

Table 4

Membrane Properties of the CHARMM DPPC Membrane at 323 K

resolution	setting	area per lipid (nm2)	membrane thickness (nm)
AA	standard	0.622 ± 0.010a	3.82 ± 0.07a
CG	hybrid	0.633 ± 0.089	4.05 ± 0.04
VS hybrid	hybrid	0.681 ± 0.014	3.66 ± 0.07
VS hybrid	hybrid with C1-VC1 LJ ε = 3.25 kJ/mol	0.621 ± 0.008	3.85 ± 0.05

Standard deviation obtained after the value reached equilibrium. Standard errors are <0.1%.

Discussion

In this study, we propose a VS hybrid model to bridge AA and CG resolutions in biomolecular membrane simulations. The accuracy of the hybrid scheme compared with the standard GROMOS/CHARMM or Martini models is satisfactory, with a significant increase in the computational efficiency.

Although our method is straightforward to use, some choices for the setup need careful consideration. Whereas the Martini force field is mostly used in combination with reaction field for the long-range electrostatics, here, we use PME to be compatible with most of the popular AA force fields. Thus, our VS hybrid scheme can be further expanded to combine Martini with other AA force fields, e.g., OPLS,[53] AMBER,[54] etc., in the future. Furthermore, we have used the flat bottom potential to prevent AA water from penetrating into the membrane. In principle, one could also add artificial repulsion between AA water atoms and beads in CG leaflets to avoid the flat bottom potential restraint. To test this idea, we explored a pure repulsive LJ potential (C12 of 10–7 (kJ nm12)/mol and C6 of 0 (kJ nm6)/mol) between the AA water and CG tail beads. When AA water crosses the AA leaflet, the repulsion between the CG leaflet and AA water effectively blocks their further penetration. Note that occasional entering of AA solvent molecules into the CG region is not a problem per se, as they will effectively behave as an ideal gas in the absence of interactions with the CG solvent.

In addition, coupling an AA force field to the Martini model requires some fine-tuning of the interaction parameters to account for the different settings. We showed that, to match the GROMOS force field, a uniform scaling of all Martini LJ interactions by a factor 0.898 suffices to compensate for the use of a longer cutoff. Matching the CHARMM force field, however, required a different approach in which specific cross interactions between the AA and CG models were optimized. Our method works best on membranes for which the area per lipid in both AA and CG resolutions is very similar, because the two leaflets in one membrane are represented in different resolutions. A potential mismatch in the area per lipid could be remedied by introducing additional lipids in one of the leaflets, but this is not ideal.

Concerning the computational efficiency of our method, we showed that the VS hybrid scheme can increase the speed of the simulation by about 5-fold in case of a small vesicle. This speedup will be enhanced if larger vesicles are considered. The efficiency can be further improved by applying mean field force approximation boundary potentials, that replace both the internal and external excess bulk solvent around a membrane[39] or by using the Dry Martini force field[55] to model the outer leaflet. In a planar membrane, the computational efficiency benefit also exists, albeit much smaller (2-fold with the current settings), for any simulations requiring two membrane environments. The computational efficiency can potentially be further increased by using a multiple time step approach for different resolution levels, such as reversible RESPA.[56]

In order to apply our hybrid method, it is essential that the two combined force fields are compatible or can be made compatible, as illustrated in this paper. We have successfully combined both the GROMOS and CHARMM force fields with the Martini force field in the DPPC membrane. We expect that the GROMOS/Martini VS hybrid scheme can be easily applied to simulate any type of lipid or lipid mixtures without further modification, as long as the area per lipid between the two resolutions agrees. However, to combine CHARMM (or other AA force fields) and Martini through the VS hybrid scheme, reparameterization of the cross interactions might be needed for lipid types that feature different atom types in the tails.

In terms of potential applications of the VS hybrid membrane model, we could think of simulations of cytosolic solutions inside a liposome, or simulations of curvature effects on the lipid packing and interaction of other compounds with the inner leaflet. In both examples, the interior region of interest is modeled in full AA detail, whereas the less interesting (but necessarily included) outer region is modeled at the CG level. Simulations that include two planar membranes could also benefit from our hybrid model, e.g., to simulate the preliminary phase (stalk formation) of membrane fusion[57] or the effect of asymmetric ionic concentrations across the membrane,[58] etc. Note that this scheme only gains computational speedup in a double membrane setup. Compared with an AA system, our hybrid scheme replaces one solvent compartment and two membrane leaflets with its CG representations. The amount of computational speedup strongly depends on the size of the CG solvent compartment. In addition, the hybrid scheme could be used to drive the atomistic region across potential barriers. For instance, the Martini force field has successfully captured spontaneous separation of ternary membranes into a liquid-disordered and a liquid-ordered phase[59,60] that can hardly be reached by simulation with all atomistic details.[61,62] One can imagine applying the VS hybrid method on such a ternary membrane and expect the phase separation of the CG leaflets to accelerate and guide this process in the AA leaflets.

Conclusion

In conclusion, we demonstrated a VS based hybrid method that can concurrently couple AA and CG force fields in molecular dynamics simulations of membranes. We have tested the VS hybrid method on both a DPPC planar membrane and a DPPC vesicle and obtained reasonable structural and conformational properties compared to AA reference simulations of the same system. The computational speedup of our VS hybrid method is largely increased compared to AA simulations, in particular for the vesicular geometry. We expect that our new multiscale method finds applications in simulating membrane related processes at large spatiotemporal scales, as well as to be useful in the ongoing development of related multiscale simulation schemes.