iSoLF is a coarse-grained (CG) model for lipid molecules with the implicit-solvent approximation used in molecular dynamics (MD) simulations of biological membranes. Using the original iSoLF (iSoLFv1), MD simulations of lipid bilayers consisting of either POPC or DPPC and these bilayers, including membrane proteins, can be performed. Here, we improve the original model, explicitly treating the electrostatic interactions between different lipid molecules and adding CG particle types. As a result, the available lipid types increase to 30. To parameterize the potential functions of the new model, we performed all-atom MD simulations of each lipid at three different temperatures using the CHARMM36 force field and the modified TIP3P model. Then, we parameterized both the bonded and non-bonded interactions to fit the area per lipid and the membrane thickness of each lipid bilayer by using the multistate Boltzmann Inversion method. The final model reproduces the area per lipid and the membrane thickness of each lipid bilayer at the three temperatures. We also examined the applicability of the new model, iSoLFv2, to simulate the phase behaviors of mixtures of DOPC and DPPC at different concentrations. The simulation results with iSoLFv2 are consistent with those using Dry Martini and Martini 3, although iSoLFv2 requires much fewer computations. iSoLFv2 has been implemented in the GENESIS MD software and is publicly available.

All-atom (AA) molecular dynamics (MD) is a powerful computational tool for investigating various biological processes in membrane environments.1–3 To simulate moderately large biological systems for relevant time scales, a tremendous computational cost is still necessary despite the recent advancement in both hardware and software.4–7 One popular alternative to scale up in size and time MD simulations of biological systems is using coarse-grained (CG) models. They can reduce the number of particles significantly compared to those with AA models and thereby accelerate MD simulations to reach biologically relevant time scales. Although the atomic details are lost in CG models, essential features of biomolecular structure and dynamics are still obtained in MD simulations.8 It should be noted that there is more than one recipe for coarse-graining the same biomolecules.9–12 The selection of the CG model used in MD simulations, therefore, depends on the scientific question to be addressed.

CG models have been successfully applied to MD simulations of lipid membranes,13–16 allowing the investigation of slow biological processes that would be prohibitively expensive in traditional AA MD.17–19 The resolution of these models can widely vary from one20,21 to a few22–24 and up to several particles per lipid.25,26 The coarser the mapping, the lower the computational cost. However, this comes at the expense of structural information, becoming a trade-off between computational cost and atomistic detail. Therefore, careful consideration is required in the selection of a CG model for each simulation.

The parameters in CG models can be determined following either a bottom–up or a top–down strategy. In the former, the main objective is to reproduce the structural properties of reference data, usually obtained from atomistic MD trajectories, whereas, in the latter, the aim is to reproduce experimental measurements of properties of interest. Therefore, the CG model following the bottom–up approach might lose accuracy when using it at conditions different from their parameterization. Also, it might not be easy to reproduce structural details using the CG model following the top–down strategy. For these reasons, many CG models adopt both strategies, at least partially. The electrostatics based (ELBA) polarized force field27,28 is parameterized following the top–down approach by reproducing with accuracy different properties of lipids, such as density, surface tension, and viscosity. The Southamerican initiative for a rapid and accurate Hamiltonian (SIRAH) model29 was initially developed for proteins and afterward extended to lipids.30 It reproduces the physical properties of lipid bilayers, such as the area per lipid, the membrane thickness, and the order parameters. There are models such as the surface property fitting coarse-graining (SPICA),26 which follow the bottom–up and top–down strategies for the bonded and non-bonded interactions, respectively. This allows the SPICA model to reproduce, for example, the phase separation in lipid bilayers at different cholesterol concentrations. One of the most popular models, Martini,25 also follows top–down and bottom–up strategies for parameterizing the interactions. It has been extended to proteins,31 carbohydrates,32 and other molecules.33 Martini is one of the most mature force fields, with an ecosystem built around it to facilitate the usage for studying a wide variety of systems.

There are also CG models of lipid membranes with implicit solvent approximations where the effect of water molecules is included in the interaction potentials. The main advantage of this implicit representation is the further reduction of the computational cost, allowing the simulation of even larger systems for prolonged periods. The implicit solvent variant in the Martini force field family is called Dry Martini.34 Most notably, it has been used by Pezeshkian et al. to simulate an entire mitochondrial membrane.35 Currently, Dry Martini is parameterized only for lipids, with the protein force field still under development.

One recent alternative is the Implicit Solvent Lipid Force Field (iSoLF).24 This lipid model is designed to have a resolution similar to the Cα protein models, making it suitable to simulate lipid-protein environments in combination with the structure-based CG models.36 One lipid in iSoLF consists of a linear chain containing between 4 and 6 particles, reducing computational costs further compared to Martini or Dry Martini. However, the original version of iSoLF (iSoLFv1) has several limitations: for instance, only two phospholipids (POPC and DPPC) are available together with proteins, and electrostatic interactions are not treated explicitly but are included just as effective potentials. These might cause difficulties in increasing the different types of lipid molecules available in the CG model.

In this work, we take the first step to extend and improve the iSoLF model for the new version (iSoLFv2) by addressing existing shortcomings in iSoLFv1. In the “Methods” section, we present modifications to the force field function: the explicit treatment of electrostatic interactions, the addition of polar interactions, and the introduction of new mappings to represent different lipid molecules. The “Results” section shows the parameterization results of 30 different phospholipids. We simulate the phase behavior of mixed lipid membranes composed of DOPC and DPPC at different concentrations using iSoLFv2 and compare the results with Martini v3 and Dry Martini. Finally, we discuss the status and the limitations of the iSoLFv2 model and future directions for simulating biological membranes with various lipid molecules and membrane proteins.

In both iSoLFv1 and iSoLFv2, two-tailed lipids are represented as single-tailed molecules composed of CG particles classified as hydrophilic head and hydrophobic tail beads (Fig. 1). In iSoLFv1, there are only two different head beads, H1 and H2, corresponding to the phosphatidylcholine and the phosphate bonded to the glycerol group, respectively, for POPC and DPPC. These beads are located at the center of mass of the corresponding groups they represent. The H2 particle is located between the phosphate and glycerol groups. In iSoLFv2, electrostatic interactions are explicitly treated to increase the number of available lipid types. The H2 head particle is separated into two independent CG particles, PHO and MID, representing the phosphate and the glycerol groups, respectively. Additionally, the H1 particle representing the choline group in iSoLFv1 is renamed to CHO in iSoLFv2, and four new particles are added: ETH, GLY, SER, and PHA. The first three particles correspond to the ethanolamine, glycerol, and L-serine head groups, respectively. The last one, PHA, corresponds to the unsubstituted phosphate group of phosphate lipids. In the case of the lipid tails, either two or three beads are used depending on the length of the conforming fatty acids. For 1,2-dilauroyl and 1,2-dimyristoyl lipids, two beads are used, whereas for 1,2-dipalmitoyl, 1,2-dioleoyl, 1-palmitoyl-2-oleoyl, and 1-stearoyl-2-oleoyl lipids, three beads are used. By combining tail and head beads, 30 different lipids in total can be represented with the new mapping (Fig. 1, Table S1).

FIG. 1.

Coarse-grained mapping. The coarse-grained mapping used in iSoLFv2 for different lipid molecules. (a) Representation of atoms contained in each coarse-grained bead for all the particles in the iSoLFv2 model. Gray, red, and blue colored groups represent beads with a net electric charge of 0e, −1e, and 1e, respectively. (b) Comparison between the mapping in iSoLFv1 and iSoLFv2. In iSoLFv2, the phosphate and the glycerol groups in each phospholipid are represented with individual beads. (c) Table summarizing the lipid available in iSoLFv2. Different lipids are represented by combining the corresponding tail and head beads.

FIG. 1.

Coarse-grained mapping. The coarse-grained mapping used in iSoLFv2 for different lipid molecules. (a) Representation of atoms contained in each coarse-grained bead for all the particles in the iSoLFv2 model. Gray, red, and blue colored groups represent beads with a net electric charge of 0e, −1e, and 1e, respectively. (b) Comparison between the mapping in iSoLFv1 and iSoLFv2. In iSoLFv2, the phosphate and the glycerol groups in each phospholipid are represented with individual beads. (c) Table summarizing the lipid available in iSoLFv2. Different lipids are represented by combining the corresponding tail and head beads.

Close modal
In the force field of iSoLFv2, we add two interaction terms, VElectrostatic and VPolar, to the original form in iSoLFv1. We also modify the functional form of the VAttraction potential, resulting in the following six terms:
V=VBond+VAngle+VRepulsion+VAttraction+VElectrostatic+VPolar.
(1)
The first term, VBond, represents the virtual bond interaction between two consecutive CG particles in the same lipid molecule. It is modeled as
VBond=i=1nbondskbond,ibib0,i2,
(2)
where kbond,i is the force constant of the ith virtual bond, bi is the ith virtual bond length, b0,i is the equilibrium distance of the ith virtual bond, and nbonds is the total number of virtual bonds. The second term, VAngle, represents the virtual angle between three consecutive CG particles in the same lipid molecule and is modeled as
VAngle=i=1nangleskangle,iθiθ0,i2,
(3)
where kangle,i is the force constant of the ith virtual angle, θi is the ith virtual angle, θ0,i is the equilibrium angle of the ith virtual angle, and nangles is the total number of virtual angles. The third term, VRepulsion, represents the short-range repulsion between two CG particles, and it is modeled with the Weeks–Chandler–Andersen potential as
VRepulsion=i<jnnlpairs4εijσijrij12σijrij6+14,rij26σij0,rij>26σij,
(4)
where ɛij is the force scaling factor between the ij particle pair, σij is the repulsive range for the ij particle pair, rij is the distance between the ith and the jth particle, and nnlpairs is the total number of non-local pairs of CG beads, i.e., pairs of beads not participating in the same virtual bond or virtual angle interactions. The fourth term, VAttraction, represents the attractive hydrophobic interaction. It is redefined in iSoLFv2 as
VAttraction=i<jnnlpairsεij,rij26σijεij2rij26σijωij33rij26σijωij2+1,26σij<rij26σij+ωij0,26σij+ωij<rij,
(5)
where ɛij and σij have the same values in the VRepulsion potential, ωij is the attraction range for the ij particle pair, rij is the distance between the ith and the jth particle, and nnlpairs is the total number of non-local pairs of CG beads. This function has been previously used to model attractive interactions in liquid–liquid phase separation systems.37 It has a continuous derivative and is computationally less expensive compared to the original potential function that requires the calculation of trigonometric functions. Additionally, it has a tunable width via the ω parameter (Fig. S1), which is important for the lateral fluidity as it has been previously shown by Cooke and Deserno.38 
The fifth term, VElectrostatic, represents the electrostatic interactions between charged particles. It is modeled using the Debye–Hückel theory,
VElectrostatic=i<jnnlpairsqiqj4πε0εrrijerij/λD,
(6)
where qi and qj are the elementary charges of the ith and jth particles, nnlpairs is the total number non-local pairs of CG beads, ɛ0 is the dielectric permittivity of vacuum, ɛr is the relative permittivity of saline water defined as εr=eTa(C), where eT=249.40.788T+7.2×104T2 is a function of the temperature39 and aC=10.2551C+5.151×102C26.889×103C3 is a function of the salt molarity,40 and λD is the Debye screening length defined as λD=kBTε0εr/2NAec2I where kB is the Boltzmann constant, T is the temperature, NA is the Avogadro’s number, ec is the elementary charge, and I is the ionic strength of the solution. The final term, VPolar, represents the interaction between charged and neutral polar particles, and it is modeled as the attractive half of the Lennard–Jones potential,
VPolar=i<jnnlpairsεij,rij26σij,4εijσijrij12σijrij6,rij>26σij,
(7)
where ɛij and σij have the same values in the VRepulsion potential, and nnlpairs is the total number of non-local pairs of CG beads. Combination rules are used to define the parameters in each non-local pair as follows: εij=εiεj, σij = (σi + σj)/2, and ωij = (ωi + ωj)/2.

We perform AA MD simulations using GENESIS41,42 v2.0 and the CHARMM force field.43 We prepare the initial configurations for 30 single-component lipid bilayers using CHARMM-GUI’s membrane builder utility.44–48 Each bilayer is composed of 200 lipids with 100 in each leaflet. Each system contains 20 000 water molecules (100 water molecules per lipid) and is embedded in a 150 mM solution of NaCl. Next, we equilibrate the bilayers using the protocol provided by CHARMM-GUI at three different temperatures of 303.15, 313.15, and 323.15 K and a pressure of 1.0 atm. Then, we perform production runs for 1.2 µs in the NPT ensemble at the same temperature and pressure with a time step of 0.0025 ps, saving the coordinates every 10 ps. Finally, we discard the first 200 ns and use for analysis the remaining 1.0 µs, corresponding to the last 100 000 frames of the trajectories.

We perform CG MD simulations with iSoLFv2 implemented in an in-house version of GENESIS. For the parameter optimization, we prepare the initial configuration for each of the 30 single-component lipid bilayers using GENESIS_CG_tools, a set of scripts written in the Julia language.49 Each system consists of lipids positioned in a grid of 25 × 25 nm2 such that they are close to their corresponding area per lipid at 303.15 K. Table S2 summarizes, for each lipid, the number of molecules in each leaflet. We perform a short equilibration for 6.01 × 105 steps in the NPγT ensemble using Langevin dynamics at three different temperatures (303.15, 313.15, and 323.15 K), a surface tension of 0.0 dyn/cm, a salt concentration of 150 mM, a time step of 10 fs, and with friction parameters for the thermostat and barostat of 0.1 ps−1. During the first 1000 simulation steps, the neighbor list is updated every two steps, and afterward, every ten steps saving coordinates every 500 steps. We use the final 1.0 × 105 steps of the trajectory for calculating properties (Fig. S2).

Three lipid bilayers composed of DOPC/DPPC and DOPC/DOPS mixtures at different relative concentrations (75%/25%, 50%/50%, and 25%/75%) are prepared for CG MD simulations with iSoLFv2 using GENESIS_CG_tools (Table S2). We also compare CG MD simulations of equivalent systems using GROMACS50 v2022 with the Dry Martini34 and GROMACS v2023 with the Martini 325 force fields. We build two bilayer systems of 2048 lipids (1024 per bilayer) using CHARMM-GUI’s martini builder utility.51,52 We equilibrate the Dry Martini and Martini 3 systems using the protocol provided by CHARMM-GUI and then perform production runs using the protocol provided on the Martini webpage. We generate trajectories of 2.5 × 107, 5 × 107, and 1 × 108 steps using a time step of 40, 20, and 10 fs and saving the coordinates every 2.5 × 103, 5 × 103, and 1 × 104 steps for the Dry Martini, Martini 3, and iSoLFv2 systems, respectively. Finally, we use, for analysis, the last 100 frames of each trajectory.

For parameterizing the force field, we calculate the area per lipid (APL) and the membrane thickness in both AA and CG MD simulations. For APL, we use the following formula:
APL=2Axynlipids,
(8)
where Axy is the cross-sectional area of the simulation box in the XY plane, and nlipids is the total number of lipids forming a bilayer. For membrane thickness, we use the distance between the average positions of the beads representing the phosphate group, PHA or PHO, in each bilayer leaflet and refer to it as the phosphate–phosphate distance (PPD) in this work. For the AA MD simulations, we obtain the position of the PHO group by calculating the center of mass of the atoms involved in the PHO group. In contrast, for the CG simulations, we use the position of PHA or PHO for lipids containing the PA head group and lipids without the PA head group, respectively. For the multicomponent CG membrane simulations, we calculate the order parameter, SCC,53 using the following formula:
SCC=123cos2θ1,
(9)
where θ is the angle between the vector connecting two consecutive tail beads and the axis perpendicular to the lipid bilayer (z axis in this work), and angular brackets denote the average over all the consecutive bead pairs in the lipid tails. Finally, to compare iSoLFv2, Martini 3, and Dry Martini, we calculate the distribution of the number of neighbors and the bending modulus. The number of neighbors for each lipid component is calculated by selecting each lipid in the bilayer and counting the neighboring lipids of the same type. We use in-house Python scripts and the MDAnalysis54 and LiPyphilic55 Python libraries to calculate the order parameter and the neighbor count. For the neighbor counting, we use a cutoff distance of 1.2 nm and select GL1 and GL2 for Martini and Dry Martini, and MID for iSoLF for reference. To estimate the bending modulus, we use a relation for the fluctuation modes derived from the Canham–Helfrich theory,
hq2=kBTL2(kcq4+γq2),
(10)
where hq2 is the power spectrum of the height fluctuation, q is the magnitude of the wave vector q=qxqy, kc is the bending modulus, γ is the surface tension, kB is the Boltzmann constant, T is the temperature, and L = Lx = Ly is the lateral length of the membrane.38,56–58 For tensionless membranes, Eq. (10) simplifies to
hq2=kBTkcq4,
(11)
from which kc is obtained by fitting q−4 against hq2 for values of q/π ≤ 0.35 nm−1. Both q and hq are calculated using a nonlocal weighted average with parameters α = 8 and δ = 2.0 nm.59 We use the last 1000 frames of each trajectory in blocks of 100 frames for calculating errors in the bending modulus estimation.

We parameterized the intramolecular and intermolecular interactions separately, similar to the strategy used in iSoLFv1 (Fig. 2). For the intramolecular interactions, VBond and VAngle, we use the multistate iterative Boltzmann Inversion method (MS-IBI).60 This method is an extension of the original single-state IBI method and aims to improve the quality of the fitted potentials by increasing the target data averaged over states. We use AA MD simulations of 30 lipid bilayers as reference datasets, each simulated for 1 µs at three different temperatures. We map the atoms in each lipid to the corresponding CG particles. Table S3 specifies the atom names in the CHARMM topology belonging to each CG particle. After one iteration of the MS-IBI method, we obtain the distributions for each lipid's bond and angle terms. We then fit the intramolecular interaction potentials to their corresponding distributions, obtaining parameters representing the original AA structures in various temperatures. However, for DMPA, DMPE, DPPC, and DPPG, we use only the distributions at 323.15 K, and for DMPS, the distributions at 313.15 and 323.15 K, since these lipids at lower temperatures were not in the liquid-disordered phase.

FIG. 2.

Parameterization strategy. Workflow for the parameterization strategy used in iSoLFv2. Parameters for the intra- and inter-molecular interactions are calculated from all-atom simulations. Intra-molecular interactions are determined with the Boltzmann Inversion method, while inter-interactions are fitted to reproduce the area per lipid and the membrane thickness.

FIG. 2.

Parameterization strategy. Workflow for the parameterization strategy used in iSoLFv2. Parameters for the intra- and inter-molecular interactions are calculated from all-atom simulations. Intra-molecular interactions are determined with the Boltzmann Inversion method, while inter-interactions are fitted to reproduce the area per lipid and the membrane thickness.

Close modal
We also follow a bottom–up approach instead of the top–down approach used in iSoLFv1 for parameterizing the intermolecular interactions. We opt for this approach since there is no experimental data for all 30 target lipids at different temperatures. We select two properties for fitting the interactions, the area-per-lipid, APL, and the membrane thickness, PPD, represented by the average distance between the phosphate groups in each bilayer. To reproduce the properties calculated in the AA simulations, we define a score function for guiding the optimization process given by
S=t30°C,40°C,50°CAPLCG,tAPLAA,tAPLAA,t+PPDCG,tPPDAA,tPPDAA,t,
(12)
S*=S+t30°C,40°C,50°CHt,
(13)
where S* represents the score and is a function of the force field parameters, APL represents the area per lipid, PPD represents the membrane thickness, and H is a regularization function that takes the value of 1 when a pore is formed in the membrane; otherwise, it is equal to 0. The H term is necessary to avoid finding parameter sets that cause the spontaneous formation of pores. We minimize S* by exploring the parameter space in steps that gradually decrease size. First, for each parameter p to be optimized, a range Lp;Up is defined, which might include a global minimum. Then, S* is evaluated at three points: Lp, Mp, and Up, where Mp = (Lp + Up)/2. If the lowest value of S* is obtained at Mp, then the range is shrunk to Lp+Mp2;Up+Mp2; otherwise, if the lowest value occurs at Lp or Up, the range is updated to 2LpMp;Mp or Mp;2UpMp, respectively. This updating procedure continues until each parameter is within the desired range.

The optimization is performed in four phases. First, we start with the optimization of the parameters for 1,2-dipalmitoyl-sn-glycero-3-phosphate (DLPA). In iSoLFv2, the particles used to model this lipid are PHA, MID, DL1, and DL2. However, even with four particles, exploring the entire parameter space is not computationally efficient. Thus, we reduce the dimensionality by setting ɛPHA = ɛMID and σPHA = σMID for the head beads and ɛDL1 = ɛDL2, σDL1 = σDL2, and ωDL1 = ωDL2 for the tail beads. Then, in the second phase, we independently optimize the σ and ɛ parameters for PC, PE, PG, and PS in the DLPC, DLPE, DLPG, and DLPS lipids while keeping all the other parameters constant and by setting the parameters of PHO equal to the ones of PHA. In the third phase, we relax the constraints for the head beads and optimize the parameters for PHO and MID independently while keeping all the other parameters constant. Finally, in the fourth phase, we keep the parameters for the head beads constant and fine-tune the parameters for the tail beads. Since all the DL lipids share the same parameters for the tails, the tuning of ɛ, σ, and ω values is guided by the minimization of the global score function, calculated as the sum of the individual score function for each DL lipid. This ensures that the final parameters for the tail beads best reproduce the properties across different lipids. For the other lipids, we use the tuned parameters for the head beads and start from the fourth phase.

We perform AA MD simulations of single-component lipid bilayers for each one of the 30 lipids considered in this work at three different temperatures. The resulting trajectories are used to parameterize both the bonded and nonbonded interactions of the CG model. In the original iSoLFv1 model, we used experimental values for fitting the nonbonded interactions. However, due to the limited experimental measurements available in the literature at different temperatures, all the reference values used in the parameterization of iSoLFv2 are obtained from atomistic MD simulation trajectories using GENESIS. Tables S4 and S5 summarize the obtained values for the APL and PPD, respectively. These values are consistent with reported experimental measurements and previous calculations with other MD software.43,61–66 Additionally, we observe that the APL increases and the PPD decreases by increasing the temperature from 303.15 to 323.15 K, which is expected for both properties.

Bonded interactions (bond and angle) are parameterized using the MS-IBI method to match the distributions obtained from our AA MD simulations. In the iSoLF model, two-tailed lipids are represented as single-chain molecules. Thus, lipid conformations where the two tails are partially extended along the bilayer plane cannot be represented. This becomes evident when observing the Boltzmann inverted potentials for the bond interactions. Figure S3 shows these potentials for the five bond interactions in the DOPC lipid. From these interactions, only the third, fourth, and fifth bonds involve particles representing the lipid tails (DO1, DO2, or DO3). As a result, the potentials present long tails on the left side, which extend until values close to 0 Å, corresponding to structures where two consecutive beads are almost superimposed. These structures cannot be represented using a harmonic potential and, thus, represent a limitation of the model. A similar explanation follows for the Boltzmann inverted potentials for the angle interaction (Fig. S4).

For tuning the nonbonded interactions, we optimize a cost function by gradually adjusting the parameters for VRepulsion, VAttraction, and VPolar potentials until the APL and PPD at different temperatures are reproduced. In our initial score function, S in Eq. (12), we only considered the L1 distances between the target and the calculated APL and PPD values. This score function guided the optimization toward parameter sets that produced stable lipid bilayers that reproduced the target APL and PPD values. However, we found that for some lipids, the resulting parameter sets favored the formation of stable pores in the membranes. These pores, being comparably small relative to the simulation box size, introduce minor deviations in the APL but have almost no effect on the PPD and cannot be detected by S. To avoid reaching those parameter sets after the optimization process, we include an additional term, H, to our score function, and name the new function S* in Eq. (13). H is defined such that it has a value of 0 when there is no pore in the membrane and a value of 1 otherwise and acts as a regularization factor during the optimization process. With this new function, we could find parameter sets that produce pore-free lipid bilayers (Fig. 3). Interestingly, for 1,2-dioleoyl lipids, the optimization process guided S* toward small values of σDP in the parameter space. Even though pores are not formed during the optimization, parameter sets with σDP<7.1Å do not destabilize pores when used in pore-containing membranes. Therefore, only for 1,2-dioleoyl, we set σDP=7.1Å.

FIG. 3.

Score function modification. (a) and (c) Parameterizing inter-molecular parameters by minimizing a score function (S) consisting of the distance between the calculated and target properties leads to membranes presenting pores (red circles). (b) and (d) A modified score function (S*) obtained by adding a regularization term prevents reaching parameter sets that stabilize pores. The red arrow indicates a point where further modification of the parameters allows the stabilization of pores.

FIG. 3.

Score function modification. (a) and (c) Parameterizing inter-molecular parameters by minimizing a score function (S) consisting of the distance between the calculated and target properties leads to membranes presenting pores (red circles). (b) and (d) A modified score function (S*) obtained by adding a regularization term prevents reaching parameter sets that stabilize pores. The red arrow indicates a point where further modification of the parameters allows the stabilization of pores.

Close modal

Finally, using the function S*, we obtain parameter sets for all 30 lipids. An example of a complete parameterization is shown in Fig. S5. Compared with iSoLFv1, these new parameters reproduce the APL and PPD more faithfully in a broader range of temperatures (Fig. 4). It is interesting to note that lipids presenting a phase transition between 30 and 50 °C reproduce these properties partially, showing that our current model still cannot reproduce accurately transition temperatures (Figs. S6–S10). A list of all the parameters for the bonded and nonbonded interactions can be found in Tables S6–S8, with a summary of the interactions between different particles in Table S9.

FIG. 4.

APL and PPD of 1,2-dioleoyl (DO) lipids. By adjusting iteratively the parameters involved in the non-bonded interaction terms, APL (a) and PPD (b) calculated in all-atom simulations were reproduced at different temperatures.

FIG. 4.

APL and PPD of 1,2-dioleoyl (DO) lipids. By adjusting iteratively the parameters involved in the non-bonded interaction terms, APL (a) and PPD (b) calculated in all-atom simulations were reproduced at different temperatures.

Close modal

To test the quality of the parameters in iSoLFv2, we performed CG MD simulations of multicomponent lipid bilayers composed of DOPC and DPPC at different concentrations. The phase behavior of these multicomponent lipid bilayers has been studied experimentally in the presence of cholesterols by Veatch and Keller.67 In the absence of cholesterols, three zones are identified in the ternary phase diagram. At low concentrations of DPPC, it is completely dissolved in the DOPC phase, giving a single liquid-disordered phase (Ld). By increasing the concentration of DPPC at around 25%, a gel phase starts to form, producing the coexistence of both liquid-disordered and gel phases (Ld/Lβ). This coexistence phase continues until DPPC reaches a threshold value of around 75%, after which a single gel phase (Lβ) remains.

In CG simulations with iSoLFv2, we track the phase behavior of each component in the lipid bilayer by observing the distribution of the order parameter for each lipid type. DOPC in a single component membrane naturally occurs in a Ld phase, presenting a peak of the order parameter around 0.55. On the other hand, DPPC stays in the Lβ phase and presents a peak of around 0.88. In multicomponent DOPC/DPPC simulations, we investigate the phase behaviors as a function of the relative concentrations. At a low DPPC concentration, the peaks of DOPC and DPPC are located around 0.56 and 0.66, respectively [Fig. 5(a)], indicating that both lipids are in the Ld phase [Fig. 5(d)]. In the second simulation, the membrane is in an equimolar concentration of DPPC and DOPC. We observe a shift in the peak of DPPC to the right and a spread of the distribution for DOPC [Fig. 5(b)]. This corresponds to DPPC lipids mainly in the Lβ and DOPC lipids in both Lβ and Ld phases [Fig. 5(e)]. Finally, in our third system, we increase the concentration of DPPC lipids, effectively inverting the concentration in our initial simulation. In this case, we observe the appearance of peaks for DOPC and DPPC around 0.86 and 0.82, respectively [Fig. 5(c)]. The location of these peaks shows that both DOPC and DPPC are in the Lβ phase [Fig. 5(f)]. Figure 6 shows the three DOPC/DPPC systems' top and side views of the last frame in each trajectory.

FIG. 5.

Order parameter for DPPC/DOPC mixtures. Distribution of the order parameter for DPPC and DOPC at different concentrations. (a) At low DPPC concentrations, both lipids are present in a liquid-disordered phase. (b) At equal concentrations, the peak of DPPC is shifted to the right, indicating this lipid's presence in the gel phase coexisting with POPC in the liquid phase. (c) At high DPPC concentrations, both lipids are present in the gel phase. (d)–(f) Overlay conformations of representative lipids aligned in the glycerol particle (MID) for the systems in (a)–(c). DPPC and DOPC molecules are colored red and blue, respectively. Lipids tails in the liquid-disordered phase have greater mobility than in the gel phase.

FIG. 5.

Order parameter for DPPC/DOPC mixtures. Distribution of the order parameter for DPPC and DOPC at different concentrations. (a) At low DPPC concentrations, both lipids are present in a liquid-disordered phase. (b) At equal concentrations, the peak of DPPC is shifted to the right, indicating this lipid's presence in the gel phase coexisting with POPC in the liquid phase. (c) At high DPPC concentrations, both lipids are present in the gel phase. (d)–(f) Overlay conformations of representative lipids aligned in the glycerol particle (MID) for the systems in (a)–(c). DPPC and DOPC molecules are colored red and blue, respectively. Lipids tails in the liquid-disordered phase have greater mobility than in the gel phase.

Close modal
FIG. 6.

Top and side views of DPPC/DOPC systems. (a)–(c) Top and (d)–(f) side view of the DPPC/DOPC systems at different relative concentrations. At low DPPC concentration, the system is in the Ld phase. At equimolar concentration, small domains of DPPC in the Lβ phase and DOPC lipids in the Ld phase coexist. At high DPPC concentration, the system presents a single Lβ phase.

FIG. 6.

Top and side views of DPPC/DOPC systems. (a)–(c) Top and (d)–(f) side view of the DPPC/DOPC systems at different relative concentrations. At low DPPC concentration, the system is in the Ld phase. At equimolar concentration, small domains of DPPC in the Lβ phase and DOPC lipids in the Ld phase coexist. At high DPPC concentration, the system presents a single Lβ phase.

Close modal

We further investigate the behavior of lipid mixtures by comparing the neighborhood of each lipid species across different coarse-grained models. Figure 7 shows the calculated distributions of the number of neighbors for DPPC/DOPC systems at the same three relative concentrations previously described using iSoLFv2, Martini 3, and Dry Martini. These lipids share the same head group (PC) but differ in the tail particle types. In each distribution, the neighbor count is performed over lipids of the same type. Thus, when either DPPC or DOPC presented a low concentration, the peak in each distribution is shifted to the left. For the case where both lipids are present with the same concentration, the peaks are centered around the same neighbor count. A comparison in analogous systems composed of DOPS/DOPC lipids is provided in Fig. S11. In these systems, the lipid components share the same tail particle types (DO1, DO2, and DO3) but differ in the head group and present at the three concentrations a similar trend to the DPPC/DOPC systems. We also evaluate the elastic properties by comparing the bending modulus of the previously described lipid mixtures. Table I and Table S10 show the estimated values from the power spectra (Figs. S12 and S13) for the DPPC/DOPC and DOPS/DOPC systems, respectively, using the three CG models. We observe that both iSoLF and Martini 3 follow an increasing trend in the bending modulus, whereas values obtained with the Dry Martini model remain constant. Finally, we calculated the mixing entropy as described in the supplementary methods for the DPPC/DOPC and DOPS/DOPC systems, observing the phase coexistence in the DPPC/DOPC system (Fig. S14). Overall, we find that the neighborhoods are consistent among different models, and the concentration dependence of the bending modulus in Martini 3 is qualitatively reproduced in iSoLFv2.

FIG. 7.

Neighbor distribution. Neighbor distribution comparison for different concentrations of DPPC and DOPC lipids calculated with iSoLFv2, Dry Martini, and Martini 3. For a given lipid, a surrounding lipid in the same leaflet is considered a neighbor if the distance between any of the MID particles in iSoLFv2 and GL1 or GL2 in Dry Martini and Martini 3 is less than 1.2 nm.

FIG. 7.

Neighbor distribution. Neighbor distribution comparison for different concentrations of DPPC and DOPC lipids calculated with iSoLFv2, Dry Martini, and Martini 3. For a given lipid, a surrounding lipid in the same leaflet is considered a neighbor if the distance between any of the MID particles in iSoLFv2 and GL1 or GL2 in Dry Martini and Martini 3 is less than 1.2 nm.

Close modal
TABLE I.

Bending modulus. Bending modulus of lipid bilayers containing DPPC and DOPC lipids at 25%/75%, 50%/50%, and 75%/25% relative concentrations and at a temperature of 303.15 K.

Bending modulus for DPPC/DOPC mixtures (1020J)
DPPC/DOPC iSoLFv2 Dry martini Martini 3 
25%/75% 7.73±0.03 6.49±0.18 5.21±0.09 
50%/50% 9.82±0.12 6.44±0.10 5.41±0.08 
75%/25% 14.01±0.27 6.42±0.15 6.21±0.07 
Bending modulus for DPPC/DOPC mixtures (1020J)
DPPC/DOPC iSoLFv2 Dry martini Martini 3 
25%/75% 7.73±0.03 6.49±0.18 5.21±0.09 
50%/50% 9.82±0.12 6.44±0.10 5.41±0.08 
75%/25% 14.01±0.27 6.42±0.15 6.21±0.07 

One of the main limitations in iSoLFv1 is the implicit representation of charge interactions. This limitation is especially significant when dealing with protein–lipid interactions and when trying to capture the properties of lipid head groups with different net charges. Thus, one of the key improvements in the iSoLFv2 model is the explicit treatment of electrostatic interactions using the Debye–Hückel model, allowing the representation of charged lipids such as phosphatidyl-L-serine and phosphatidylglycerol. Additionally, the Debye–Hückel model considers the effect of the solvent by screening the electrostatic interactions as a function of the temperature of the solvent and the ionic concentration, opening the possibility to capture the dependency of bilayer properties on the salt concentration by further tuning the parameters. It is important to mention that since the dielectric constant inside lipid bilayers is comparably different from aqueous environments, we might need to refine the electrostatic interactions further, and it represents a limitation of the implicit solvent nature of the model. In the iSoLFv1 model, the center of mass of the H2 bead is localized between the phosphate and glycerol group, adding a non-zero charge to this particle would introduce artifacts due to the localization of the charge. Thus, the iSoLFv2 model represents the phosphate and glycerol with individual beads and assigns integer charges of −1e, 0e, and +1e to the negatively, neutral, and positively charged particles, respectively (Fig. 1). These integer charges are fixed and do not change during the parameterization of the non-bonded interactions. Interestingly, despite presenting charged groups in the atomic structure, the SER head group is represented with a neutral bead in the iSoLFv2 model. Thus, if we exclude the polar interaction, the SER lipid head would only interact with other particles with the repulsive potential. In fact, in the early stages of this work, we did not consider such polar interaction, making the optimization algorithm lead us to a trivial solution for the parameter set of the excluded volume interaction of the SER bead. This can be explained by analyzing the structures, APLs, and PPDs of 1,2-dipalmitoyl-sn-glycero-3-phosphate (DLPA) and 1,2-dipalmitoyl-sn-glycero-3-phospho-L-serine (DLPS). Despite differing structurally only in the presence or absence of the SER lipid head bead, both lipids present similar values for the APL and PPD and keep this trend at different temperatures (Tables S4 and S5). Since the similar property values produce similar scores using Eq. (13), the repulsive parameters for PS always converged toward zero, effectively reducing the lateral repulsion between DLPS lipids and making DLPA and DLPS the same lipid. Therefore, by adding the polar interaction between SER and charged groups, DLPS lipids could approach the target APL and PPD not necessarily by decreasing the repulsion but also by increasing the attraction between polar and charge particles, resulting in a non-trivial solution.

We find that our current parameterization best reproduces the properties of lipids that do not present a phase transition from Lβ to Ld in the temperature range of 30–50 °C, as we observe by the deviations in the APL and PPD when comparing CG and AA simulations. This happens mainly when the phase transition temperature was found in the parameterization range of 30–50 °C since our model probably cannot correctly capture both phases. Additionally, the spontaneous formation of pores in membranes is addressed in iSoLFv2 by including a regularization term in the score function, as opposed to fixing the ratio between the head and tail particles as in iSoLFv1. Replacing this regularization term with explicit measurements of the line tension68,69 will improve the score function's quality for guiding the force field's parameterization.

Finally, by reproducing the properties of lipids used routinely in experiments at different temperatures with the new parameter set, we can apply the iSoLFv2 model to a wider range of lipid environments.70,71 However, key components are still missing in the iSoLFv2 model. The first one is the representation of sterols, specifically cholesterol. In its current state, we can represent only lipid membranes with a low to zero cholesterol content. These molecules are involved in various biological processes72–75 and regulate lipid bilayers' mechanical and chemical properties. The second missing component is the representation of sugars. Sugars are involved in several signaling processes and including them in the model unlocks the representation of glycolipids by combining them with the current parameter set.76 Finally, the third and most important component is the reparameterization of the lipid–protein interactions with the new parameter set since one of the key features of the iSoLF model is compatibility with other implicit solvent CG models. Decoupling the electrostatic interaction from the current, lipid–protein interaction potential is required before iSoLFv2 can be used with other protein models, and it will allow the simulation of proteins with different protonation states, something that is not possible with the current iSoLFv1 model. We are currently working on this last point, adjusting the current parameter set36 to reproduce the orientation and dynamics of a wide variety of membrane proteins, and we expect soon to be able to use iSoLFv2 not only alongside the AICG2+ protein model77 but also with other models for studying intrinsically disordered proteins78 and liquid–liquid phase separation.37,79,80

In this work, we have improved and extended the iSoLF implicit-solvent lipid force field, addressing various shortcomings. We have added an explicit treatment for the electrostatic interactions between charged particles and polar interactions between charged and neutral polar head particles into the force field. Additionally, we have modified the functional form for the attractive interaction between the hydrophobic tails by replacing it with a polynomial switch function, decreasing the computational cost associated with this interaction. To better represent the charge interaction between particles, we have changed the CG mapping from AA structures, representing the phosphate and the glycerol group in each lipid with individual particles instead of the one-particle representation in the original iSoLF model. By fitting the parameters to 30 different lipids, the force field can reproduce the area per lipid and the membrane thickness calculated in AA simulations at three different temperatures. Finally, we have examined the behavior of lipid bilayers composed of DOPC and DPPC at different concentrations and compared it with similar systems using the Martini 3 and Dry Martini force fields. Together, these improvements allow the iSoLF model to explore the behavior of large multicomponent systems. It will become a valuable tool once we update the interaction between lipids and proteins with the iSoLFv2 model and combine it with existing CG models. The iSoLFv2 model is publicly available in the GENESIS MD software.

See the supplementary material for additional information on the iSoLFv2 force fields, parameterizations, and comparisons between all-atom and CG MD simulations of various lipids.

This research was partially supported by RIKEN Pioneering Research Projects (Dynamic Structural Biology/Glycolipidologue Initiative) (to Y.S.), MEXT Program: Data Creation and Utilization-Type Material Research and Development Project Grant No. JPMXP1122714694 (to Y.S.), and MEXT/KAKENHI under Grant Nos. JP19H05645 and JP21H05249 (to Y.S.). We used the computer systems HOKUSAI (Project ID: Q22536, Q23536), provided by the RIKEN Information System Division, and FUGAKU (Project ID: ra00003).

The authors have no conflicts to disclose.

Diego Ugarte La Torre: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Shoji Takada: Conceptualization (supporting); Writing – review & editing (supporting). Yuji Sugita: Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
I. D.
Pogozheva
,
G. A.
Armstrong
,
L.
Kong
,
T. J.
Hartnagel
,
C. A.
Carpino
,
S. E.
Gee
,
D. M.
Picarello
,
A. S.
Rubin
,
J.
Lee
,
S.
Park
,
A. L.
Lomize
, and
W.
Im
,
J. Chem. Inf. Model.
62
,
1036
(
2022
).
2.
J.
Rizo
,
L.
Sari
,
Y.
Qi
,
W.
Im
, and
M. M.
Lin
,
eLife
11
,
e76356
(
2022
).
3.
G.
Shahane
,
W.
Ding
,
M.
Palaiokostas
, and
M.
Orsi
,
J. Mol. Model
25
,
76
(
2019
).
4.
J.
Jung
,
W.
Nishima
,
M.
Daniels
,
G.
Bascom
,
C.
Kobayashi
,
A.
Adedoyin
,
M.
Wall
,
A.
Lappala
,
D.
Phillips
,
W.
Fischer
,
C.
Tung
,
T.
Schlick
,
Y.
Sugita
, and
K. Y.
Sanbonmatsu
,
J. Comput. Chem.
40
,
1919
(
2019
).
5.
S.
Páll
,
A.
Zhmurov
,
P.
Bauer
,
M.
Abraham
,
M.
Lundbor
,
A.
Gray
,
B.
Hess
, and
E.
Lindahl
,
J. Chem. Phys.
153
,
134110
(
2020
).
6.
J. C.
Phillips
,
D. J.
Hardy
,
J. D. C.
Maia
,
J. E.
Stone
,
J. V.
Ribeiro
,
R. C.
Bernardi
,
R.
Buch
,
G.
Fiorin
,
J.
Hénin
,
W.
Jiang
,
R.
McGreevy
,
M. C. R.
Melo
,
B. K.
Radak
,
R. D.
Skeel
,
A.
Singharoy
,
Y.
Wang
,
B.
Roux
,
A.
Aksimentiev
,
Z.
Luthey-Schulten
,
L. V.
Kalé
,
K.
Schulten
,
C.
Chipot
, and
E.
Tajkhorshid
,
J. Chem. Phys.
153
,
044130
(
2020
).
7.
P.
Eastman
,
J.
Swails
,
J. D.
Chodera
,
R. T.
McGibbon
,
Y.
Zhao
,
K. A.
Beauchamp
,
L.
Wang
,
A. C.
Simmonett
,
M. P.
Harrigan
,
C. D.
Stern
,
R. P.
Wiewiora
,
B. R.
Brooks
, and
V. S.
Pande
,
PLoS Comput. Biol.
13
,
e1005659
(
2017
).
8.
W. G.
Noid
,
J.
Chu
,
G. S.
Ayton
,
V.
Krishna
,
S.
Izvekov
,
G. A.
Voth
,
A.
Das
, and
H. C.
Andersen
,
J. Chem. Phys.
128
,
244114
(
2008
).
9.
S. J.
Marrink
,
V.
Corradi
,
P. C. T.
Souza
,
H. I.
Ingolfsson
,
D. P.
Tieleman
, and
M. S. P.
Sansom
,
Chem. Rev.
119
,
6184
(
2019
).
10.
S.
Kmiecik
,
D.
Gront
,
M.
Kolinski
,
L.
Wieteska
,
A. E.
Dawid
, and
A.
Kolinski
,
Chem. Rev.
116
,
7898
(
2016
).
11.
J.
Li
and
S.
Chen
,
Front. Mol. Biosci.
8
,
720937
(
2021
).
12.
T.
Sun
,
V.
Minhas
,
N.
Korolev
,
A.
Mirzoev
,
A. P.
Lyubartsev
, and
L.
Nordenskiöld
,
Front. Mol. Biosci.
8
,
645527
(
2021
).
13.
J.
Parchekani
,
A.
Allahverdi
,
M.
Taghdir
, and
H.
Naderi-Manesh
,
Sci. Rep.
12
,
2371
(
2022
).
14.
T.
Mandal
,
S. E.
Spagnolie
,
A.
Audhya
, and
Q.
Cui
,
Biophys. J.
120
,
3211
(
2021
).
15.
M. I.
Mahmood
,
H.
Noguchi
, and
K.
Okazaki
,
Sci. Rep.
9
,
14557
(
2019
).
16.
A.
Yu
,
A. J.
Pak
,
P.
He
,
V.
Monje-Galvan
,
L.
Casalino
,
Z.
Gaieb
,
A. C.
Dommer
,
R. E.
Amaro
, and
G. A.
Voth
,
Biophys. J.
120
,
1097
(
2021
).
17.
C.
Chng
,
K. J.
Hsia
, and
C.
Huang
,
Soft Matter
18
,
7752
(
2022
).
18.
P.
Appshaw
,
A. M.
Seddon
, and
S.
Hanna
,
Soft Matter
18
,
1747
(
2022
).
19.
I.
Mannelli
,
F.
Sagues
,
V.
Pruneri
, and
R.
Reigada
,
Langmuir
32
,
12632
(
2016
).
20.
H.
Yuan
,
C.
Huang
,
J.
Li
,
G.
Lykotrafitis
, and
S.
Zhang
,
Phys. Rev. E
82
,
011905
(
2010
).
21.
H.
Noguchi
and
M.
Takasu
,
J. Chem. Phys.
115
,
9547
(
2001
).
22.
S.
Foley
and
M.
Deserno
,
J. Chem. Theory Comput.
16
,
7195
(
2020
).
23.
R.
Goetz
and
R.
Lipowsky
,
J. Chem. Phys.
108
,
7397
(
1998
).
24.
D.
Ugarte La Torre
and
S.
Takada
,
J. Chem. Phys.
153
,
205101
(
2020
).
25.
P. C. T.
Souza
,
R.
Alessandri
,
J.
Barnoud
,
S.
Thallmair
,
I.
Faustino
,
F.
Grünewald
,
I.
Patmanidis
,
H.
Abdizadeh
,
B. M. H.
Bruininks
,
T. A.
Wassenaar
,
P. C.
Kroon
,
J.
Melcr
,
V.
Nieto
,
V.
Corradi
,
H. M.
Khan
,
J.
Domański
,
M.
Javanainen
,
H.
Martinez-Seara
,
N.
Reuter
,
R. B.
Best
,
I.
Vattulainen
,
L.
Monticelli
,
X.
Periole
,
D. P.
Tieleman
,
A. H.
de Vries
, and
S. J.
Marrink
,
Nat. Methods
18
,
382
(
2021
).
26.
S.
Seo
and
W.
Shinoda
,
J. Chem. Theory Comput.
15
,
762
(
2018
).
27.
M.
Orsi
and
J. W.
Essex
,
PLoS one
6
,
e28637
(
2011
).
28.
R. M.
de Souza
,
R. H.
Ratochinski
,
M.
Karttunen
, and
L. G.
Dias
,
J. Chem. Inf. Model.
60
,
522
(
2020
).
29.
L.
Darre
,
M. R.
Machado
,
A. F.
Brandner
,
H. C.
González
,
S.
Ferreira
, and
S.
Pantano
,
J. Chem. Theory Comput.
11
,
723
(
2015
).
30.
E. E.
Barrera
,
M. R.
Machado
, and
S.
Pantano
,
J. Chem. Theory Comput.
15
,
5674
(
2019
).
31.
L.
Monticelli
,
S. K.
Kandasamy
,
X.
Periole
,
R. G.
Larson
,
D. P.
Tieleman
, and
S. J.
Marrink
,
J. Chem. Theory Comput.
4
,
819
(
2008
).
32.
V.
Lutsyk
,
P.
Wolski
, and
W.
Plazinski
,
J. Chem. Theory Comput.
18
,
5089
(
2022
).
33.
S. J.
Marrink
,
L.
Monticelli
,
M. N.
Melo
,
R.
Alessandri
,
D.
Peter Tieleman
, and
P. C. T.
Souza
,
Wiley Interdiscip. Rev. Comput. Mol. Sci.
13
,
e1620
(
2023
).
34.
C.
Arnarez
,
J. J.
Uusitalo
,
M. F.
Masman
,
H. I.
Ingólfsson
,
D. H.
de Jong
,
M. N.
Melo
,
X.
Periole
,
A. H.
de Vries
, and
S. J.
Marrink
,
J. Chem. Theory Comput.
11
,
260
(
2015
).
35.
W.
Pezeshkian
,
M.
König
,
T. A.
Wassenaar
, and
S. J.
Marrink
,
Nat. Commun.
11
,
2296
(
2020
).
36.
D.
Ugarte La Torre
and
S.
Takada
,
J. Chem. Phys.
155
,
155101
(
2021
).
37.
Y.
Murata
,
T.
Niina
, and
S.
Takada
,
Biophys. J.
121
,
4382
(
2022
).
38.
I. R.
Cooke
and
M.
Deserno
,
J. Chem. Phys.
123
,
224710
(
2005
).
39.
A.
Catenaccio
,
Y.
Daruich
, and
C.
Magallanes
,
Chem. Phys. Lett.
367
,
669
(
2003
).
40.
A.
Stogryn
,
IEEE Trans. Microwave Theory Tech.
19
,
733
(
1971
).
41.
C.
Kobayashi
,
J.
Jung
,
Y.
Matsunaga
,
T.
Mori
,
T.
Ando
,
K.
Tamura
,
M.
Kamiya
, and
Y.
Sugita
,
J. Comput. Chem.
38
,
2193
(
2017
).
42.
J.
Jung
,
T.
Mori
,
C.
Kobayashi
,
Y.
Matsunaga
,
T.
Yoda
,
M.
Feig
, and
Y.
Sugita
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
5
,
310
(
2015
).
43.
J. B.
Klauda
,
R. M.
Venable
,
J. A.
Freites
,
J. W.
O’Connor
,
D. J.
Tobias
,
C.
Mondragon-Ramirez
,
I.
Vorobyov
,
A. D.
MacKerell
, Jr.
, and
R. W.
Pastor
,
J. Phys. Chem. B
114
,
7830
(
2010
).
44.
S.
Jo
,
T.
Kim
,
V. G.
Iyer
, and
W.
Im
,
J. Comput. Chem.
29
,
1859
(
2008
).
45.
E. L.
Wu
,
X.
Cheng
,
S.
Jo
,
H.
Rui
,
K. C.
Song
,
E. M.
Dávila-Contreras
,
Y.
Qi
,
J.
Lee
,
V.
Monje-Galvan
,
R. M.
Venable
,
J. B.
Klauda
, and
W.
Im
,
J. Comput. Chem.
35
,
1997
(
2014
).
46.
S.
Jo
,
J. B.
Lim
,
J. B.
Klauda
, and
W.
Im
,
Biophys. J.
97
,
50
(
2009
).
47.
S.
Jo
,
T.
Kim
, and
W.
Im
,
PLoS One
2
,
e880
(
2007
).
48.
J.
Lee
,
D. S.
Patel
,
J.
Ståhle
,
S.
Park
,
N. R.
Kern
,
S.
Kim
,
J.
Lee
,
X.
Cheng
,
M. A.
Valvano
,
O.
Holst
,
Y. A.
Knirel
,
Y.
Qi
,
S.
Jo
,
J. B.
Klauda
,
G.
Widmalm
, and
W.
Im
,
J. Chem. Theory Comput.
15
,
775
(
2019
).
49.
C.
Tan
,
J.
Jung
,
C.
Kobayashi
,
D. U. L.
Torre
,
S.
Takada
, and
Y.
Sugita
,
PLoS Comput. Biol.
18
,
e1009578
(
2022
).
50.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
,
SoftwareX
1-2
,
19
(
2015
).
51.
Y.
Qi
,
H. I.
Ingólfsson
,
X.
Cheng
,
J.
Lee
,
S. J.
Marrink
, and
W.
Im
,
J. Chem. Theory Comput.
11
,
4486
(
2015
).
52.
P.
Hsu
,
B. M. H.
Bruininks
,
D.
Jefferies
,
P.
Cesar Telles de Souza
,
J.
Lee
,
D. S.
Patel
,
S. J.
Marrink
,
Y.
Qi
,
S.
Khalid
, and
W.
Im
,
J. Comput. Chem.
38
,
2354
(
2017
).
53.
S.
Seo
,
M.
Murata
, and
W.
Shinoda
,
J. Phys. Chem. Lett.
11
,
5171
(
2020
).
54.
N.
Michaud-Agrawal
,
E. J.
Denning
,
T. B.
Woolf
, and
O.
Beckstein
,
J. Comput. Chem.
32
,
2319
(
2011
).
55.
P.
Smith
and
C. D.
Lorenz
,
J. Chem. Theory Comput.
17
,
5907
(
2021
).
56.
W.
Helfrich
,
Z. Naturforsch. C
28
,
693
(
1973
).
57.
F. L. H.
Brown
,
Annu. Rev. Phys. Chem.
59
,
685
(
2008
).
58.
P. W.
Fowler
,
J.
Hélie
,
A.
Duncan
,
M.
Chavent
,
H.
Koldsø
, and
M. S. P.
Sansom
,
Soft Matter
12
,
7792
(
2016
).
59.
A. K.
Chaurasia
,
A. M.
Rukangu
,
M. K.
Philen
,
G. D.
Seidel
, and
E. C.
Freeman
,
Phys. Rev. E
97
,
032421
(
2018
).
60.
T. C.
Moore
,
C. R.
Iacovella
, and
C.
McCabe
,
J. Chem. Phys.
140
,
224104
(
2014
).
61.
J.
Lee
,
X.
Cheng
,
J. M.
Swails
,
M. S.
Yeom
,
P. K.
Eastman
,
J. A.
Lemkul
,
S.
Wei
,
J.
Buckner
,
J. C.
Jeong
,
Y.
Qi
,
S.
Jo
,
V. S.
Pande
,
D. A.
Case
,
C. L.
Brooks
III
,
A. D.
MacKerell
, Jr.
,
J. B.
Klauda
, and
W.
Im
,
J. Chem. Theory Comput.
12
,
405
(
2016
).
62.
A. P.
Lyubartsev
and
A. L.
Rabinovich
,
Biochim. Biophys. Acta, Biomembr.
1858
,
2483
(
2016
).
63.
X.
Zhuang
,
E. M.
Dávila-Contreras
,
A. H.
Beaven
,
W.
Im
, and
J. B.
Klauda
,
Biochim. Biophys. Acta, Biomembr.
1858
,
3093
(
2016
).
64.
X.
Zhuang
,
J. R.
Makover
,
W.
Im
, and
J. B.
Klauda
,
Biochim. Biophys. Acta, Biomembr.
1838
,
2520
(
2014
).
65.
K.
Pluhackova
,
S. A.
Kirsch
,
J.
Han
,
L.
Sun
,
Z.
Jiang
,
T.
Unruh
, and
R. A.
Böckmann
,
J. Phys. Chem. B
120
,
3888
(
2016
).
66.
N.
Kucerka
,
M.
Nieh
, and
J.
Katsaras
,
Biochim. Biophys. Acta, Biomembr.
1808
,
2761
(
2011
).
67.
S. L.
Veatch
and
S. L.
Keller
,
Biophys. J.
85
,
3074
(
2003
).
68.
T. V.
Tolpekina
,
W. K.
den Otter
, and
W. J.
Briels
,
J. Chem. Phys.
121
,
8014
(
2004
).
69.
H.
Noguchi
,
Soft Matter
8
,
3146
(
2012
).
70.
R. B.
Lira
,
T.
Robinson
,
R.
Dimova
, and
K. A.
Riske
,
Biophys. J.
116
,
79
(
2019
).
71.
I.
Gozen
and
P.
Dommersnes
,
Eur. Phys. J.: Spec. Top.
223
,
1813
(
2014
).
72.
N. K.
Khadka
,
R.
Timsina
,
E.
Rowe
,
M.
O’Dell
, and
L.
Mainali
,
Biochim. Biophys. Acta, Biomembr.
1863
,
183625
(
2021
).
73.
S.
Chakraborty
,
M.
Doktorova
,
T. R.
Molugu
,
F. A.
Heberle
,
H. L.
Scott
,
B.
Dzikovski
,
M.
Nagao
,
L.
Stingaciu
,
R. F.
Standaert
,
F. N.
Barrera
,
J.
Katsaras
,
G.
Khelashvili
,
M. F.
Brown
, and
R.
Ashkar
,
Proc. Natl. Acad. Sci.
117
,
21896
(
2020
).
74.
H. L.
Nguyen
,
V. H.
Man
,
M. S.
Li
,
P.
Derreumaux
,
J.
Wang
, and
P. H.
Nguyen
,
Phys. Chem. Chem. Phys.
24
,
6225
(
2022
).
75.
Z.
Al-Rekabi
and
S.
Contera
,
Proc. Natl. Acad. Sci.
115
,
2658
(
2018
).
76.
B.
Binotti
,
R.
Jahn
, and
A.
Perez-Lara
,
Arch. Biochem. Biophys.
709
,
108966
(
2021
).
77.
W.
Li
,
W.
Wang
, and
S.
Takada
,
Proc. Natl. Acad. Sci.
111
,
10550
(
2014
).
78.
Y.
Zhang
,
X.
Liu
, and
J.
Chen
,
J. Chem. Inf. Model.
62
,
4523
(
2022
).
79.
C.
Tan
,
A.
Niitsu
, and
Y.
Sugita
,
JACS Au
3
,
834
(
2023
).
80.
S. Y.
Cheng
,
Y.
Cao
,
M.
Rouzbehani
, and
K. H.
Cheng
,
Biophys. Chem.
260
,
106355
(
2020
).

Supplementary Material