# Development of Coarse-Grained Lipid force Fields Based on a Graph Neural Network

Zhenyu Liao, Ting Si, Tairan Wang, Ji-Jung Kai, Christophe Chipot, Jun Fan

► To cite this version:

Zhenyu Liao, Ting Si, Tairan Wang, Ji-Jung Kai, Christophe Chipot, et al.. Development of Coarse-Grained Lipid force Fields Based on a Graph Neural Network. Journal of Chemical Theory and Computation, 2025, 21 (18), pp.9175-9185. 10.1021/acs.jctc.5c01071. hal-05412283

HAL Id: hal-05412283

https://hal.science/hal-05412283v1

Submitted on 11 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Development of Coarse-Grained Lipid force
Fields Based on a Graph Neural Network

Zhenyu Liao<fnref n="1" />, Ting Si <fnref n="1" />,<fnref n="2" />, Tairan Wang<fnref n="1" />, Ji-Jung Kai <fnref n="3" />,<fnref n="4" />, Christophe Chipot<fnref n="5" />, Jun Fan<fnref n="1" />,<fnref n="3" />,<fnref n="4" /> *

**Corresponding Author**

<footnote n="1">¹ Department of Materials Science and Engineering, City University of Hong Kong, Kowloon 999077, Hong Kong, China;</footnote>

<footnote n="2">² Department of Physics, City University of Hong Kong, Kowloon 999077, Hong Kong, China;</footnote>

<footnote n="3">³ Department of Mechanical Engineering, City University of Hong Kong, Kowloon 999077, Hong Kong, China;</footnote>

<footnote n="4">⁴Centre for Advanced Nuclear Safety and Sustainable Development, City University of Hong Kong, Kowloon 999077, Hong Kong, China;</footnote>

<footnote n="5">⁵ Laboratoire International Associé CNRS and University of Illinois at Urbana-Champaign, UMR no. 7019, Université de Lorraine, BP 70239, F-54506 Vandœuvre-lès-Nancy, France.</footnote>

<footnote>* J. Fan: junfan@cityu.edu.hk</footnote>
**Abstract**

Coarse-grained (CG) lipid models have been employed to effectively simulate important biological events with large spatiotemporal scales on the cell membranes. However, developing fast CG lipid models with atomic-level accuracy remains challenging. A promising strategy is to use graph neural networks (GNN) to learn the CG energy functions from all-atom (AA) simulations, which has proven successful in developing the CG protein force fields. Herein, we leverage the established graph neural network TorchMD-GN, to develop the first GNN-based CG models of lipid bilayers. We trained models based on the AA simulations of the DOPC, DOPS and DOPC/DOPS lipid bilayer. Our CG lipid models accelerate the lipid dynamics by 9.4 times while achieving a similar AA-level accuracy. Temperature-dependent simulations show that the model exhibits some degree of temperature transferability. We also attempted to train the model using bicelles, which enhanced the stability of the lipids and led to better performance in self-assembly and vesicle simulations. This work demonstrates the feasibility of using GNN to build coarse-grained lipid models and provides insights for the improvement of the framework.

**Keywords**

Lipid, Coarse-grained, Force Field, Graph Neural Network, Molecular Dynamics Simulation
Introduction

Cell membranes are essential for all living cells, serving as not only physical barriers but as regulators in substrate transport and signal transduction.¹, ² Understanding the mechanism underlying cellular processes on membranes is of paramount importance. Molecular dynamics (MD) simulations have emerged as a powerful tool for elucidating membrane-related processes at molecular resolution.³⁻⁵ All-atom (AA) simulations provide high-resolution details, but their computational cost becomes prohibitively high, since many biological events on cell membranes take place over large spatiotemporal scales. To bridge this gap, coarse-grained (CG) representations have been developed as an effective strategy in the membrane simulations.⁶

Cell membranes consist mainly of lipids. Therefore, it is essential to develop accurate and computationally efficient CG lipid force fields for CG simulations of the cell membrane. CG force fields are mainly divided into two categories: top-down models and bottom-up models.⁷ The top-down models were parameterized to replicate macroscopic properties. A notable example is the Martini force field, which was fitted by reproducing experimental partition coefficients.⁸, ⁹ Martini is recognized as the most widely used CG lipid model and has been developed to the latest version, Martini 3.¹⁰ Despite its great success, the Martini model still exhibits certain structural and thermodynamic inaccuracies, such as the inability to accurately depict the phase
behavior of lipids and the inaccurate decomposition of entropy and enthalpy in
potential of mean force (PMF).<sup>11-13</sup>

The bottom-up models are fitted by replicating the microscopic properties predicted
by fine-grained simulations, typically at the atomic level.<sup>14</sup> The main task in
developing bottom-up models is to formulate an effective CG Hamiltonian. An ideal
CG Hamiltonian should reproduce the many-body PMF from AA simulations.<sup>15</sup>
However, due to limitations in sampling and the basis set of CG interactions, it is
impractical to achieve an exact match of the many-body PMF. Therefore, the key
challenge is to design effective basis functions to approximate the ideal CG
Hamiltonian. The CG basis functions are typically represented in terms of pair
potentials, such as B-spline basis functions in the force matching method.<sup>15</sup> Several
strategies have been proposed to enhance the accuracy of the CG energies, including
the introduction of virtual particles<sup>16, 17</sup>, incorporating higher-order potential terms to
account for many-body effects<sup>18</sup>, and the inclusion of explicitly expressed
electrostatic terms<sup>19, 20</sup>. Machine learning (ML) techniques provide another viable
strategy, where neural networks are employed to replace conventional pair potentials.
The neural network is trained on datasets obtained from more expensive methods and
serves as the energy function for CG simulations.<sup>21</sup>

Graph neural network exhibits good performance in predicting atomic energies and forces.<sup>22, 23</sup> This network has developed the variational force matching method into a ML framework for the development of CG force fields.<sup>24-29</sup> TorchMD-GN<sup>26</sup> is a
mature GNN architecture that demonstrates outstanding performance in protein
folding prediction. The advantage of the network lies in its multiple interaction
blocks, which can implicitly capture many-body effects by aggregating neighbor
information, thereby preserving the thermodynamics more faithfully.¹⁸ Although
membranes play a critical role in biological systems, there have not been any ML CG
force fields specifically for lipids.

In this work, we trained low-resolution. solvent-free CG models for lipid
membranes with various compositions using TorchMD-GN. We examined two types
of homogeneous membranes, including DOPC and DOPS, as well as a heterogenous
membrane composed of a mixture of DOPC and DOPS. We built a dataset of AA
simulations for these lipid membranes and then mapped the AA lipids to a CG
representation. We employed TorchMD-GN to train GNN-based CG lipid force fields
and evaluated them through several Langevin CG simulations. Our results show that
TorchMD-GN can train fast and accurate CG lipid force fields with some
transferability. We anticipate this work will contribute to the application of ML
methods in CG lipid models.
# Theory and Methodology

## Coarse-graining of Lipids

CG modeling consists of two steps: mapping from atomic structure to CG structure and parameterization of CG interactions. In this work, all coarse-grained lipids were represented by 6-site beads using center-of-mass mapping. The detailed mapping scheme can be found in Figure S1 and Table S2. The phosphate and choline/serine moieties were represented by the headgroup (HG) bead, while the glycerol and ester groups were defined as the middle-group (MG) bead. The first half and the second half of the lipid tail were represented by T1 and T2 beads, respectively. We adopted this strategy to preserve the structural features at the atomic level as much as possible, while ensuring high computational efficiency.¹⁶,³⁰ The 6-site beads model also has a similar resolution to the previously developed ML CG protein models, facilitating their combined use in future CG simulations of membrane proteins.²⁶,³¹

In the parameterization of lipid models, the variational force matching method proves to be effective.¹⁵, ³² It minimizes the mean squared error between the mapped AA forces and the CG forces derived from the CG force field, to find the optimal CG force field:

$$ \chi^2[U_{CG}] = \frac{1}{3N} \sum_{I=1}^{N} \langle |\Xi_I(F(\mathbf{r}_n)) + \nabla_I U_{CG}(\mathbf{R}_N; \theta)|^2 \rangle \quad (1) $$
where $\Xi_f(F(\mathbf{r}_n))$ is the AA forces projected onto the CG beads; $\Xi_f$ is the force mapping operator; $F(\mathbf{r}_n)$ is the instantaneous force collected from all-atom simulations; CG forces $(-\nabla_f U_{CG}(\mathbf{R}_N; \theta))$ are predicted by the CG force field; $\mathbf{R}_N$ is the CG configuration mapped from the AA configuration, $\mathbf{r}_n$.

In the classical bottom-up lipid force fields, spline basis functions are typically used to describe non-bonded interactions, and their coefficients are optimized according to Eq. 1.¹⁶ However, since spline potentials are typically pairwise models, many-body effects are neglected, which can affect the accurate simulation of lipid behavior. Graph neural networks can replace the previous pairwise models, by using a graph structure to learn information from surrounding multilayer beads to predict CG energy and forces, thereby capturing many-body effects.¹⁸, ²⁶ Variational force matching can be transformed into a machine learning approach, where Eq. 1 serves as the loss function, continuously optimizing the network parameters to develop the GNN-based force field.²⁴

## Graph Neural Network Training

Our CG lipid models were built based on the TorchMD-GN²⁶ and trained using TorchMD-Net³³. TorchMD-GN is a graph convolutional network inspired by SchNet³⁴, which demonstrates satisfactory performance in coarse-graining proteins. First, the training data had to be prepared. We collected coordinates and forces obtained from AA simulations of lipids. The coordinates and forces were then
mapped to CG beads. Before feeding the data into the network, we extracted prior energy terms from the training data to avoid exploring physically meaningless regions. The prior energy terms include bond and nonbonded repulsive terms. The bond term was approximated by the harmonic potential:

$$V_{\text{bonded}}(r) = k(r - r_0)^2 + V_0, \qquad (2)$$

where k is the spring constant, r is a pairwise distance, r₀ is an equilibrium distance, and V₀ is a base potential. The nonbonded repulsive term was represented by a function inspired by the repulsive term of Lennard-Jones potentials:

$$V_{\text{repulsive}}(r_{ij}) = 4\bar{U}_{ij}^{-6} + V_0 \qquad (3)$$

where rᵢⱼ is the pairwise distance between bead i and j, U is the repulsive constant and V₀ is a base potential. Prior forces are calculated from the prior energy term. We subtracted the prior forces from the AA forces to obtain the delta force, which is then used for network training.

The input of TorchMD-GN is Cartesian coordinates **R**<sub>N</sub> and the types of CG beads **Z**<sub>N</sub>. First, the input data was passed through an embedding layer to generate a molecular graph. Each node in the graph corresponded to a CG bead and was assigned an embedding feature vector. The edges of the graph are the pairwise distances between beads within the cutoff range. Then, the edge information is used by the interaction block to update the node features. The interaction block typically consists of multiple layers (four layers in this work), which gradually integrate information from distant bead pairs. Lastly, a multi-layer perceptron contracted the node features
into an energy scalar $E_{CG}(\mathbf{R}_N)$. The CG forces $\vec{F}_{CG}(\mathbf{R}_N)$ are derived by calculating
the gradient of the energy scalar with respect to the initial coordinates, which act as
the output of the model. CG forces will also be used to compute the force-matching
loss, thereby continuously optimizing the graph network during training. For
additional details on TorchMD-GN, readers can refer to the original literature²⁶.

All training was carried out using PyTorch Lightning. Most of the hyperparameters are consistent with the original literature, while a hyperparameter search was conducted for several key hyperparameters.²⁶ Since the test loss does not show much difference between successful and failed models, the hyperparameter search was based on the CG simulations produced by the trained models.²⁶, ²⁹, ³⁵ The metrics used for hyperparameters optimization are the stability of the bilayer in the CG simulations, and the ability to reproduce the static properties seen in all-atom simulations. The set of hyperparameters in Table S1 ensure stable CG simulations and closely match the results of AA simulations.

The training input file with the complete set of hyperparameters is available on Zenodo.

The GNN-based force field for each lipid membrane was trained 2 to 4 times with different random seeds. Quick evaluations using 1 ns CG simulations revealed little difference among models with different seeds. One of the models was selected for follow-up validation simulations. The training curves for the selected models are
shown in Figure S2. Each model was trained on 2 Nvidia GeForce RTX 2080 Ti GPUs, taking about 7 minutes per epoch.

## All-atom Simulations

CG models were built for DOPC, DOPS, and a mixture of 1:1 DOPC/DOPS lipids. Corresponding AA simulations were prepared separately for each lipid model as the training sets (Table S2). The initial configuration was generated using CHARMM-GUI Membrane Builder³⁶. The bilayer system consists of 512 lipids arranged in a bilayer within a periodic box of 14 × 14 × 10 nm³. The membrane was solved with TIP3P³⁷ water molecules in a periodic box, and were then neutralized by adding 150mM of sodium and chloride ions. Bicelle is produced using the equilibrated bilayer: the lipid bilayer is placed in a larger water box (18 × 18 × 10 nm³) to make the mirror images invisible. The bilayer will transform into bicelles after hundreds of nanoseconds simulations.

All the simulations were carried out by GROMACS-2020³⁸ along with CHARMM36 force field³⁹, ⁴⁰. Energy minimization was performed on each membrane, with 1000 KJ·mol⁻¹·nm⁻² harmonic restraints on the phosphorus atoms of lipids. The restrained system was heated to 300 K under NPT conditions in 2.5 ns. The restraints were then gradually relaxed, and two equilibrium simulations lasting 2.5 ns each were performed. Afterward, all restraints on the pre-equilibrated membrane were removed, and NPT simulations were continued until the membrane
area reached equilibrium. Each membrane underwent a 300 ns *NPT* simulation. Lastly, a 200 ns *NVT* simulation was run to collect trajectories for CG force field development, storing the coordinates and forces every 0.02 ns. As a result, the database for each model contains 10,000 frames of coordinates and forces. Pressure is kept at 1 bar with a semi-isotropic Parrinello-Rahman barostat⁴¹, while temperature is controlled by a Nose-Hoover thermostat⁴², ⁴³. A smooth cutoff of 10-12 Å and the particle-mesh Ewald⁴⁴ method were applied to treat van der Waals (vdW) and long-range electrostatic interactions, respectively. A time step was set to 2 fs, with the LINCS⁴⁵ algorithm applied to constrain all covalent bonds involving hydrogen.

## Coarse-grained Simulations

The details of simulations are summarized in Table S2. In the bilayer simulation, the initial membrane conformations are consistent with the reference AA simulations, and the box dimensions correspond to those used in the AA simulations. In the lipid assembly simulation, an equivalent number of lipids as in the training set is randomly placed in an isotropic periodic box, with the side length matching that of the AA membrane plane. The vesicle is modeled using CHARMM-GUI and placed in a periodic box of 24 × 24 × 24 nm³. The bilayer-to-vesicle simulations starts from the lipid bilayer in the training set placed in a large box (24 × 24 × 18 nm³). All the systems were modelled without solvent.
CG simulations were conducted using TorchMD<sup>46</sup>. The force field consists of two components: the force field trained by TorchMD-GN and the prior energy term. Each system underwent Langevin simulations under NVT conditions, at the same temperature 300 K as the AA simulations (unless otherwise noted). The damping constant is set to 0.1 ps<sup>-1</sup>. A time step is typically set to 1 fs or 20 fs.

The CG simulations using the Martini 3 force field followed the simulation protocol outlined in a previous benchmarking study of Martini 3.<sup>11</sup> The Martini 3 beads were mapped into 6-sites models for data analysis, with mapping rules provided in Table S2.

<table><caption>**Table S1.** Summary of MD simulations in this work.</caption><thead><tr><th>Force Field</th><th>System</th><th>Beads</th><th>Replicas</th><th>Time (ns)</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">GN_DOPC</td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>132184</td><td>2</td><td>500</td></tr><tr><td>GN_DOPC<sup>a</sup></td><td>Bilayer</td><td>3072</td><td>10</td><td>2</td></tr><tr><td>GN_DOPC</td><td>Bilayer</td><td>3072</td><td>5</td><td>100</td></tr><tr><td>GN_DOPC</td><td>Self-assembly</td><td>3072</td><td>5</td><td>10</td></tr><tr><td>GN_DOPC</td><td>Vesicle</td><td>5346</td><td>4</td><td>4</td></tr><tr><td>GN_DOPC</td><td>Bilayer-to-vesicle</td><td>3072</td><td>5</td><td>2</td></tr><tr><td>CHARMM36</td><td>Temp. 270 K</td><td>132184</td><td>1</td><td>300</td></tr><tr><td>GN_DOPC</td><td>Temp. 270 K</td><td>3072</td><td>5</td><td>2</td></tr><tr><td>CHARMM36</td><td>Temp. 330 K</td><td>132184</td><td>1</td><td>300</td></tr><tr><td>GN_DOPC</td><td>Temp. 330 K</td><td>3072</td><td>5</td><td>2</td></tr></tbody></table>
<table><tbody><tr><td>Martini 3</td><td>Bilayer</td><td>14356</td><td>2</td><td>400</td></tr><tr><td>Martini 3</td><td>Self-assembly</td><td>11008</td><td>2</td><td>400</td></tr><tr><td colspan="5">GN_DOPS</td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>140182</td><td>1</td><td>500</td></tr><tr><td>GN_DOPS</td><td>Bilayer</td><td>3072</td><td>5</td><td>2</td></tr><tr><td colspan="5">GN_PCPS<sup>b</sup></td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>140622</td><td>1</td><td>500</td></tr><tr><td>GN_PCPS</td><td>Bilayer</td><td>3072</td><td>5</td><td>2</td></tr><tr><td colspan="5">GN_DOPC_BC<sup>c</sup></td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>226182</td><td>1</td><td>500</td></tr><tr><td>GN_DOPC_BC</td><td>Self-assembly</td><td>3072</td><td>5</td><td>20</td></tr><tr><td>GN_DOPC_BC</td><td>Vesicle</td><td>5346</td><td>2</td><td>10</td></tr><tr><td>GN_DOPC_BC</td><td>Bilayer-to-vesicle</td><td>3072</td><td>5</td><td>10</td></tr><tr><td colspan="5">GN_DOPS_BC</td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>232266</td><td>1</td><td>500</td></tr><tr><td>GN_DOPC_BC</td><td>Self-assembly</td><td>3072</td><td>3</td><td>10</td></tr><tr><td colspan="5">GN_PSPC_BC</td></tr><tr><td>CHARMM36</td><td>Training Set</td><td>231512</td><td>1</td><td>500</td></tr><tr><td>GN_DOPC_BC</td><td>Self-assembly</td><td>3072</td><td>3</td><td>10</td></tr></tbody></table>

<sup>a</sup>GN_* refers to the GNN-based CG models corresponding to the lipids. All the models are available in Zenodo.

<sup>b</sup>*_PCPC represents a 1:1 mixture of DOPC/DOPS membrane

<sup>c</sup>*_BC means the model trained on bicelles.
## Analysis Details

All snapshots were visualized and rendered in Visual Molecular Dynamics (VMD)⁴⁷. The radial distribution function is calculated using the measure module of VMD. Z-density was defined as the distribution of the z-coordinates of the corresponding beads. All densities in this work have been area-normalized. The membrane thickness referred to the bin-averaged normal distance between the upper and lower HG beads, with each grid being 10 × 10 Å. Orientational order parameters were obtained from the second-order Legendre polynomial for cos(θ) as follows

$$S(\theta) = \left\langle \frac{1}{2} (3 \cos^2(\theta) - 1) \right\rangle \qquad (4)$$

where θ is the angle between the lipid normal and the vector of the T1-T2 (or HG-MG) bond. The range of order parameters is from -0.5 to 1, where the value of 1 indicates the bond is parallel to the normal direction of the membrane, while the value of -0.5 indicates an orthogonal orientation. The bending modulus was calculated by the real space fluctuation method.⁴⁸, ⁴⁹ The lateral mean square displacement (MSD) was defined by the following equation

$$MSD(t) = \left\langle \left( r_{i,xy}(t+\tau) - r_{i,xy}(\tau) \right)^2 \right\rangle \qquad (5)$$

where rᵢ,xy is the xy coordinates of the MG beads for each lipid, τ is the initial time point and t is the time interval. The diffusion coefficient was defined as a quarter of
the slope of the lateral MSD curve. The order parameters and lateral MSD were both
averaged over all lipids in the trajectories after reaching equilibrium.
## Results

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Flowchart of developing and testing a CG lipid force field based on a graph neutral network, TorchMD-GN.</figcaption></figure>

The workflow for developing and evaluating GNN-based CG lipid forces fields is illustrated in Figure 1. Firstly, AA simulations of lipid bilayers were conducted, and the coordinates and forces were recorded. These coordinates and forces were then mapped to CG representation. The mapped coordinates and their associated CG bead types were fed into a graph neutral network called TorchMD-GN. The input data went through an embedding layer, multiple interaction blocks, and an output layer, then produced the energy and force for each bead. The model was trained iteratively by minimizing the mean squared error loss between the delta forces and the predicted forces from the model, where delta forces referred to the mapped forces obtained after
subtracting the prior energy term. Lastly, extensive CG simulations were performed to systematically evaluate the capabilities and limitations of the ML CG lipids models. More details about the workflow can be found in the Theory and Methodology section.

## The structural and dynamic properties of the GN_DOPC model

DOPC is a very common phospholipid in cell membranes and has been widely used as a template in the development of the CG force field.<sup>16, 17, 50</sup> Therefore, we first trained a GNN-based CG force field for DOPC lipids, referred as GN_DOPC. To evaluate the GN_DOPC model, we conducted ten replicas of 2 ns CG simulations (with a timestep of 1 fs) and five replicas of 100 ns simulations (with a timestep of 20 fs). We performed a systematic analysis of the simulation results, comparing them with the all-atom data from the training set and the outcomes of Martini 3 CG simulations.

Firstly, the stability of the GN_DOPC model was tested. The snapshot (Figure 2a) and video (Video S1) of the simulation trajectories show that the bilayer structure of remains stable throughout the simulations. We conducted some time-dependent analysis. Figure S3 illustrates that the potential and total energies rapidly stabilize during all the simulations. In addition, we also plotted the time evolution of membrane properties, which shows no significant changes over time (Figure S4a-c). The membrane is also stable in the CG simulations with longer
timestep (Figure S4d-e). While a few lipids drift into the vacuum, the majority of
the lipids maintain the bilayer structure.

Secondly, the structural correlation of the GN_DOPC model between CG and AA
simulations was examined. The distributions of bonds and angles in the CG
simulations are consistent with those in the AA model, except for the T1-T2 bonds
(Figure S5a, b). The inset in Figure S5 shows that due to the bending of the lipid
tails, an additional peak in the distribution of the T1-T2 bonds from the mapped
AA trajectories appears around 4.5 Å. This peak is difficult to capture by the
harmonic potential of the prior bonded term. To address the issue, we attempted
to remove the prior term and trained a network directly from the mapped AA forces.
The results of simulations show that the absence of the prior term leads to significant
distortion of the membrane, with some bonds becoming overly stretched (Figure S6).
We also explored using B-splines potentials to describe the bonded terms. Third-
order B-splines of T1-T2 exhibits a local minimum at the 4.4 Å peak, which is
expected to better capture the non-Gaussian distribution of the T1-T2 bonds (Figure
S5c). These results indicate that it is still necessary to describe the bonded term
with a basis set when training the coarse-grained lipid force field, but a more
flexible functional form, such as B-splines, is recommended to be employed.

The nonbonded structural correlation was assessed by radial distribution
functions (RDFs). The RDFs for beads in both AA and CG (using GN_DOPC)
simulations exhibit a high degree of consistency, which demonstrates that
GN_DOPC can reproduce pair correlations in the AA model. (Figure 2a and S7). In contrast, the RDFs for beads in the CG simulations with Martini 3 are significantly different from those in the AA simulations.

<figure><img src="image_3.png" /></figure>

**Figure 2.** The structural features of the GN_DOPC model. (a) Snapshots of equilibrated DOPC lipid bilayer from CG simulations using the GN_DOPC model. (b) The RDFs between the same type of beads and (c) the z-density distribution of beads from the mapped AA simulations, the CG simulations using the GN_DOPC model, and the Martini 3 simulations.

The membrane structure of the GN_DOPC model was also characterized. We observed that the z-density distributions under CG becomes narrower and higher compared to AA, with the peak position of each bead moving outward from the center of the membrane by 1 Å (Figure 2c). The z-density profile of the T2 beads shows a separation of 1 Å between the upper and lower lipid layers, a feature that is absent in the AA model. Other structural properties are presented in Table 2. The thickness
$(d_{HG-HG})$, tail order parameter $(S_{TT})$ and bending modulus $(\kappa_{Bend})$ in the GN_DOPC model are both greater than in the AA model, suggesting the membrane becomes thicker and more ordered. The head group order parameter $(S_{HM})$ in the GN_DOPC model shows good agreement with the AA model, suggesting that the GN_DOPC model effectively captures the conformation of the head groups. There are also differences in membrane structural properties between the Martini and AA models, such as the peak positions of bead z-density shifting inward by 1 Å, reduced membrane thickness, and lower headgroup ordering. The calculation of the membrane structural properties indicates that there is still room for improvement in the fitting of the bead distribution along the normal direction for GN_DOPC.

<table><caption>**Table 2.** Properties of lipid bilayers, including bilayer thickness $(d_{HG-HG})$, T1-T2 $(S_{TT})$ and HG-MG $(S_{HM})$ order parameter, bending modulus $(\kappa_{Bend})$, and diffusion constant $(D_C)$ from experiments and computational models</caption><thead><tr><th>Status</th><th>d<sub>HG-HG</sub> (Å)</th><th>S<sub>TT</sub></th><th>S<sub>HM</sub></th><th>κ<sub>Bend</sub> (k<sub>B</sub>T)</th><th>D<sub>C</sub> (10<sup>-7</sup> nm/s)</th></tr></thead><tbody><tr><td>Expt.(DOPC)</td><td>44.8<sup>a</sup></td><td>N.A.</td><td>N.A.</td><td>18.3<sup>a</sup></td><td>0.82<sup>b</sup></td></tr><tr><td>AA(DOPC)</td><td>40.2</td><td>0.45</td><td>0.57</td><td>18.8 ± 0.1</td><td>0.86 ± 0.02</td></tr><tr><td>GN_DOPC</td><td>41.2</td><td>0.60</td><td>0.55</td><td>27.8 ± 0.3</td><td>78.53 ± 0.47</td></tr><tr><td>Martini (DOPC)</td><td>38.6</td><td>0.44</td><td>0.26</td><td>15.5 ± 0.1</td><td>5.99 ± 0.05</td></tr><tr><td>AA(DOPS)</td><td>43.5</td><td>0.54</td><td>0.68</td><td>21.7 ± 0.1</td><td>0.21 ± 0.02</td></tr><tr><td>GN_DOPS</td><td>43.0</td><td>0.58</td><td>0.64</td><td>29.6 ± 0.1</td><td>44.64 ± 2.35</td></tr></tbody></table>
<table><tr><td>AA(PC/PS)</td><td>42.3</td><td>0.50</td><td>0.68</td><td>17.6 ± 0.4</td><td>0.29 ± 0.01<sup>c</sup></td></tr><tr><td>CG_PCPS</td><td>42.8</td><td>0.60</td><td>0.65</td><td>29.9 ± 0.3</td><td>42.92 ± 2.75<sup>d</sup></td></tr></table>

<sup>a</sup>Reference<sup>51</sup>

<sup>b</sup>Reference<sup>52</sup>

<sup>c</sup>D<sub>C</sub> (PC) = 0.30 ± 0.01 × 10<sup>-7</sup> nm/s; D<sub>C</sub> (PS) = 0.27 ± 0.01<sup>-7</sup> × 10<sup>-7</sup> nm/s

<sup>d</sup>D<sub>C</sub> (PC) = 42.82 ± 2.01 × 10<sup>-7</sup> nm/s; D<sub>C</sub> (PS) = 43.03 ± 2.05<sup>-7</sup> × 10<sup>-7</sup> nm/s

Thirdly, the dynamic properties of lipids are then evaluated by the lateral mean square displacement (MSD) plot (Figure 3 and S8) and the diffusion constant ($D_C$) (Table 2). We found that the lateral diffusion of lipids in the AA model is relatively slow, with a diffusion constant of 0.86 ± 0.02 × 10<sup>-7</sup> nm/s. The diffusion constant from AA is consistent with the experimental results<sup>52</sup> and another simulation study<sup>53</sup> using the CHARMM36 force field. The lipid diffusion is accelerated in the Martini 3 model, with a diffusion coefficient that is 8 times that of AA. The lipid diffusion in the GN_DOPC model is roughly 100 times faster than in the AA model, due to its highest degree of coarse-graining.

<figure><img src="image_4.png" /></figure>
**Figure 3.** The dynamic properties of the membrane, represented by the lateral mean square displacement (MSD) of lipids obtained from the AA simulations, the CG simulations using the GN_DOPC model, and the CG simulation with the Martini 3 model.

We evaluated the efficiency of the GN_DOPC model and the AA model using the MSD values and wallclock time (Table S3). The wallclock time of the GN_DOPC model is similar to that of the AA model, but the diffusion of lipids is much faster, exhibiting a speedup of 9.4 times.

## Extension to charged and mixed lipids

Developing CG force fields for charged lipids is nontrivial and often requires special treatment of electrostatic interactions.¹⁹ Moreover, since membranes typically consist of various lipid types, it is essential to accurately describe the interactions between these different lipids. Therefore, we used the process of developing the GN_DOPC model to build CG force fields for charged and mixed lipids. DOPS is an anionic lipid with a charged head, which plays a vital role in cellular processes such as viral invasion and protein recognition⁵⁴, ⁵⁵. Thus, we developed a GN_DOPS model for the DOPS lipid bilayer and a GN_PCPS model for the mixed 1:1 DOPC/DOPS membrane. Snapshots and videos reveal that the GN_DOPS and GN_PCPS models can still produce lipid bilayer structures, but their stability is not as good as that of the
GN_DOPC model (Figure 4a, b and Video S2, S3). In particular, these models show defects, such as lipid heads entering the hydrophobic zone, tails protruding into the vacuum, and lipids detaching from the membrane. Time-dependent analysis shows that the membrane is stable overall, without significant changes in properties over time

<figure><img src="image_5.png" /></figure>

**Figure 4.** Extension to charged and mixed lipids. Snapshots of (a) DOPS lipids and (b) 1:1 DOPC/DOPS lipids from the simulations using GN_DOPS and GN_PCPS models, respectively. RDF between selected beads from (c) DOPS and (d) DOPC/DOPS system.

<figure><img src="image_6.png" /></figure>

**Figure 5.** Z-density profiles for beads in (a) GN_DOPS and (b) GN_PCPS model.
We analyzed the structural and dynamic properties of the lipids in the GN_DOPS and GN_PCPS models. The RDF analysis indicates that the GN_DOPS model captures the peak distances but overestimates their magnitudes (Figure 4c and S10a). In contrast, the GN_PCPS model perfectly matches the distribution of beads from AA results (Figure 4d and S10b). The z-density distributions reveal that the results from the GN_DOPS and GN_PCPS models are largely consistent with the AA results, although some overly structured features seen in the GN_DOPC model are also present (Figure 5). The thickness of the DOPS and DOPC/DOPS membranes, as well as the ordering of the head groups, remains consistent between CG and AA (Table 2). However, the ordering of the tails and the bending modulus suggest that the membranes become slightly stiffer. Moreover, similar to the GN_DOPC model, lipids in these two models exhibit faster lateral diffusion.

## Training with bicelles improves the lipid self-assembly behavior

The results above suggest that the GNN-based CG lipid models trained according to our proposed process show satisfactory performance on the training set, efficiently capturing certain structural correlations of AA. It is unclear whether the models are applicable to systems outside the training set. Therefore, we performed self-assembly, vesicle and bilayer-to-vesicle simulations with the GN_DOPC model. The initial configurations for these simulations can be found in
Figure S11, and the system setup is detailed in Theory and Methodology section. The simulation results indicate that GN_DOPC does not fully capture the amphipathic behavior of lipids (Figure 6). In the self-assembly simulation, the lipids failed to entirely form a bilayer structure. The vesicle turns into free lipids in the box. When a bilayer is placed in a larger box to block its view of periodic images, the lipids do not form bicelles or vesicles; instead, they disassemble.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Training CG lipid models with bicelles improve the model's performance in the self-assembly, vesicle and bilayer-to-vesicle simulations. Snapshots at the end frame of each simulation are presented.</figcaption></figure>

The bilayer-to-vesicle simulations offered an insight: the GN_DOPC model has only learned the lipid conformations within the bilayer, and it cannot capture the curved region in the bicelles or vesicles. To enhance the generalization ability of
the CG models, we established a new training set. We placed the lipid bilayer from
the training set into a large water box to make the periodic image invisible. In the
400 ns AA simulations, the bilayer transformed into a disc-like structure known as
bicelles. We trained a new GNN-based CG lipid force field for bicelles
(GN_DOPC_BC) and applied it to the aforementioned simulations (Figure S12).
The GN_DOPC_BC model achieves the self-assembly of lipids, and demonstrated
its stability in vesicles and bicelles (Figure 6 and Video S4-6). In the self-assembly
simulations, the GN_DOPC_BC model can capture the pair correlations under AA
simulations of the lipid bilayer (Figure S13).

We also trained a GN_DOPS_BC model with DOPS bicelles and a GN_PCPS_BC model with 1:1 DOPC/DOPS bicelles, and tested their quality in self-assembly simulations (Figure S12 and S14). Movies of trajectories show that the free lipids quickly aggregate into bilayer under these two models (Videso S7 and S8). These results suggest that training on bicelles, which have richer lipid conformations, seems to be a more reliable method for establishing GNN-based CG lipid force fields. The trained models demonstrate a stronger ability to understand amphipathic behavior of lipids.

The transferability of temperature

Another application of the CG lipid force field is to study lipid phase transition, so
it is essential to explore whether our model can accurately capture lipid phase
behavior.¹³, ⁵⁰, ⁵⁶ Experimental evidence indicates that, DOPC lipids gradually become more ordered until a liquid-to-gel phase transition occurs at 257 K.⁵⁶⁻⁵⁸ Given that the phase transition occurs at a rather low temperature, we carried out several temperature-dependent simulations to probe the temperature transferability of our GNN-based CG lipid force fields within the liquid phase. The GN_DOPC model trained at 300 K was employed to simulate the lipid bilayer at 330 K and 270 K (Figure 7a). The snapshot of simulations shows that in the CG simulation with GN_DOPC, the stable bilayer structure of the membrane is maintained at other temperatures (Figure 7a). The calculation of membrane properties indicates that as the temperature decreases, lipid diffusion slows down, and the membrane becomes thicker and more ordered (Figure 7b-d).

<figure><img src="image_8.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_11.png" /></figure>
**Figure 7.** The transferability performance of the GN_DOPC model in different temperatures. (a) Snapshot derived from CG simulations of the DOPC membrane at 330 K and 270 K. The properties of the membrane, such as (b) diffusion constant, (c) membrane thickness and order parameter, are comprehensively compared between CG with GN_DOPC and AA simulations at various temperatures.
Discussion and Conclusions

In light of the impressive success of TorchMD-GN in developing CG protein force
fields, we have applied this framework to the development of lipid force fields in this
work. As the first endeavor to develop GNN-based CG lipid force fields, we have
demonstrated the considerable promise of machine learning in the lipid system, while
also recognizing several issues that require improvement.

GNN-based CG lipid force fields demonstrate considerate accuracy in CG
simulations of membranes. The CG simulations almost reproduce the distribution of
bonds and angles, as well as RDFs between beads in the AA simulations. The
accuracy is comparable to that of traditional bottom-up CG models, showing the
feasibly of using GNN as potential functions.¹⁶, ¹⁷, ⁵⁰ However, our models slightly
overestimate the membrane thickness and the ordering of lipid tails. Similar trends
have also been reported in previous bottom-up models.¹⁶, ³⁰ For the GNN-model, we
suggest improving the limitation in two ways. The bonded term should be optimized
by modeling with more flexible potential functions, or being learned directly by GNN.
On the other hand, the environment of the lipid tails is relatively complex, with
neighbors that include particles from both the same layer and the opposite layer.
Therefore, more training data may be required for the GNN to accurately learn the
information of lipid tails.

GNN-based CG lipid force field also ensures computational efficiency. The
wallclock time of CG simulations with our models is similar to that of AA simulations,
since each step requires loading the network and calculating the forces. However, due
to fewer degrees of freedom and a smoother free energy landscape in CG simulations,
each step corresponds to a longer physical time equivalent. When achieving the same
degree of lipid diffusion, our model is significantly faster than the AA force field. The
accelerated lipid diffusion makes our models more suitable for large-scale, multi-
component complex membrane simulations. As TorchMD continues to improve, we
anticipate that our model will become more efficient in the future.⁴⁶

Our training framework is also effective for charged and mixed lipids. The
GN_PCPS model successfully reproduces the pair correlations from AA
simulations distinguish the distributions of DOPS and DOPC lipids, which is
evident in the different peak heights in the HG-related RDFs (Figure 4). These
findings show that TorchMD-GN can easily develop CG lipid models for lipid
membranes of any composition, demonstrating the broad applicability of the
framework.

We also found that training with bicelles improved the self-assembly and membrane
bending behavior of lipids. This finding indicates GNN training require a training set
with a diverse range of lipid conformations. The stability of the GN_DOPC_BC
model in vesicles and bicelles make it promise in large-scale membrane simulations,
such as whole-cell⁵⁹, virion⁶⁰, and nuclear pore complex⁶¹, ⁶² simulations. Before
implementing this in large biological systems, it is essential to model the
interactions of lipids with proteins. The next step is to develop a universal GNN-
based CG force field of protein-lipid. Our model also demonstrates a degree of
temperature transferability, attributed to the fitting of the PMF that incorporates an
entropy term, which captures some temperature dependence. However, developing a
truly temperature-transferable lipid force field is no easy task and requires the
correct partitioning of enthalpic and entropic interactions.⁶³,⁶⁴ The applicability of
this framework for lipids with complex phase behavior, such as DPPC, is worth
exploring.
**Data Availability**

All GNN-based CG lipid force fields, example of the input files for the training and
CG simulations are available at https://doi.org/10.5281/zenodo.15422583

**Supporting Information**

The Supporting Information (SI) is available free of charge at
https://pubs.acs.org/doi/xxx.

AUTHOR INFORMATION

**Corresponding Authors**

**Jun Fan** – Department of Materials Science and Engineering, Department of
Mechanical Engineering, and Centre for Advanced Nuclear Safety and Sustainable
Development, City University of Hong Kong, Kowloon 999077, Hong Kong, China;
orcid.org/0000-0001-8227-9671; Email: junfan@cityu.edu.hk

**Authors**

**Zhenyu Liao** – Department of Materials Science and Engineering, City University of Hong Kong, Kowloon 999077, Hong Kong, China; orcid.org/0000-0001-9500-3420; Email: zhenyliao4-c@my.cityu.edu.hk
**Ting Si** - Department of Materials Science and Engineering and Department of Physics, City University of Hong Kong, Kowloon 999077, Hong Kong, China; orcid.org/0000-0002-4291-9205; Email: tingsi2-c@my.cityu.edu.hk

**Tairan Wang** - Department of Materials Science and Engineering, City University of Hong Kong, Kowloon 999077, Hong Kong, China; orcid.org/0000-0002-0220-6059; Email: trwang3@cityu.edu.hk

**Ji-Jung Kai** - Department of Mechanical Engineering and Centre for Advanced Nuclear Safety and Sustainable Development, City University of Hong Kong, Kowloon 999077, Hong Kong, China; orcid.org/0000-0001-7848-8753; Email: jijkai@cityu.edu.hk

**Christophe Chipot** - Laboratoire International Associé CNRS and University of Illinois at Urbana-Champaign, UMR no. 7019, Université de Lorraine, BP 70239, F-54506 Vandœuvre-lès-Nancy, France; orcid.org/0000-0002-9122-1698; Email: chipot@illinois.edu

Complete contact information is available at:

https://pubs.acs.org/doi/xxx

## Author Contributions

J. F. led the research. Z. L. performed simulations, trained models and wrote the original paper. All authors reviewed the paper.
## Notes

The authors declare no competing financial interests.

## ACKNOWLEDGMENTS

We thank Prof. Lanyuan LU from Nanyang Technological University for insightful discussions. This work was supported by the Research Grants Council of Hong Kong (CityU 11305919 and 11308620) and NSFC/RGC Joint Research Scheme N_CityU104/19. Hong Kong Research Grant Council Collaborative Research Fund: C1002-21G and C1017-22G. This research made use of the computing resources of the X-GPU cluster supported by the Hong Kong Research Grant Council Collaborative Research Fund: C6021-19EF. This project was also supported by CLP Power Hong Kong under grants 9229033 to Prof. Ji-Jung Kai.
(1) Grecco, Hernán E.; Schmick, M.; Bastiaens, Philippe I. H. Signaling from the Living Plasma Membrane. *Cell* **2011**, *144*, 897-909.

(2) Edidin, M. Lipids on the frontier: a century of cell-membrane bilayers. *Nat. Rev. Mol. Cell Biol.* **2003**, *4*, 414-418.

(3) Marrink, S. J.; Corradi, V.; Souza, P. C. T.; Ingólfsson, H. I.; Tieleman, D. P.; Sansom, M. S. P. Computational Modeling of Realistic Cell Membranes. *Chem. Rev.* **2019**, *119*, 6184-6226.

(4) Muller, M. P.; Jiang, T.; Sun, C.; Lihan, M.; Pant, S.; Mahinthichaichan, P.; Trifan, A.; Tajkhorshid, E. Characterization of Lipid-Protein Interactions and Lipid-Mediated Modulation of Membrane Protein Function through Molecular Simulation. *Chem. Rev.* **2019**, *119*, 6086-6161.

(5) Stevens, J. A.; Grünewald, F.; van Tilburg, P. A. M.; König, M.; Gilbert, B. R.; Brier, T. A.; Thornburg, Z. R.; Luthey-Schulten, Z.; Marrink, S. J. Molecular dynamics simulation of an entire cell. *Front. Chem.* **2023**, *11*, 1106495.

(6) Ingólfsson, H. I.; Lopez, C. A.; Uusitalo, J. J.; de Jong, D. H.; Gopal, S. M.; Periole, X.; Marrink, S. J. The power of coarse graining in biomolecular simulations. *WIREs Comput. Mol. Sci.* **2014**, *4*, 225-248.

(7) Noid, W. G. Perspective: Coarse-grained models for biomolecular systems. *J. Chem. Phys.* **2013**, *139*, 090901.

(8) Marrink, S. J.; Monticelli, L.; Melo, M. N.; Alessandri, R.; Tieleman, D. P.; Souza, P. C. T. Two decades of Martini: Better beads, broader scope. *WIREs Comput. Mol. Sci.* **2023**, *13*, e1620.

(9) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. *J. Phys. Chem. B* **2007**, *111*, 7812-7824.

(10) Souza, P. C. T.; Alessandri, R.; Barnoud, J.; Thallmair, S.; Faustino, I.; Grünewald, F.; Patmanidis, I.; Abdizadeh, H.; Bruininks, B. M. H.; Wassenaar, T. A.; et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. *Nat. Methods* **2021**, *18*, 382-388.

(11) Loose, T. D.; Sahrmann, P. G.; Qu, T. S.; Voth, G. A. Changing Your Martini Can Still Give You a Hangover. *J. Chem. Theory Comput.* **2024**, *20*, 9190-9208.

(12) Carpenter, T. S.; López, C. A.; Neale, C.; Montour, C.; Ingólfsson, H. I.; Di Natale, F.; Lightstone, F. C.; Gnanakaran, S. Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined Martini Coarse-Grained Force Field. *J. Chem. Theory Comput.* **2018**, *14*, 6050-6062.

(13) Sharma, P.; Desikan, R.; Ayappa, K. G. Evaluating Coarse-Grained MARTINI Force-Fields for Capturing the Ripple Phase of Lipid Membranes. *J. Phys. Chem. B* **2021**, *125*, 6587-6599.
(14) Jin, J.; Pak, A. J.; Durumeric, A. E. P.; Loose, T. D.; Voth, G. A. Bottom-up Coarse-Graining: Principles and Perspectives. *J. Chem. Theory Comput.* **2022**, *18*, 5759-5791.

(15) Noid, W. G.; Chu, J.-W.; Ayton, G. S.; Krishna, V.; Izvekov, S.; Voth, G. A.; Das, A.; Andersen, H. C. The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models. *J. Chem. Phys.* **2008**, *128*, 244114.

(16) Pak, A. J.; Dannenhoffer-Lafage, T.; Madsen, J. J.; Voth, G. A. Systematic Coarse-Grained Lipid Force Fields with Semiexplicit Solvation via Virtual Sites. *J. Chem. Theory Comput.* **2019**, *15*, 2087-2100.

(17) Sahrmann, P. G.; Loose, T. D.; Durumeric, A. E. P.; Voth, G. A. Utilizing Machine Learning to Greatly Expand the Range and Accuracy of Bottom-Up Coarse-Grained Models through Virtual Particles. *J. Chem. Theory Comput.* **2023**, *19*, 4402-4413.

(18) Wang, J.; Charron, N.; Husic, B.; Olsson, S.; Noé, F.; Clementi, C. Multi-body effects in a coarse-grained protein force field. *J. Chem. Phys.* **2021**, *154*, 164113.

(19) Srivastava, A.; Voth, G. A. A Hybrid Approach for Highly Coarse-grained Lipid Bilayer Models. *J Chem Theory Comput* **2013**, *9*, 750-765.

(20) Srivastava, A.; Voth, G. A. Solvent-Free, Highly Coarse-Grained Models for Charged Lipid Systems. *J. Chem. Theory Comput.* **2014**, *10*, 4730-4744.

(21) Sahrmann, P. G.; Voth, G. A. On the emergence of machine-learning methods in bottom-up coarse-graining. *Curr. Opin. Struct. Biol.* **2025**, *90*, 102972.

(22) Schütt, K. T.; Arbabzadah, F.; Chmiela, S.; Müller, K. R.; Tkatchenko, A. Quantum-chemical insights from deep tensor neural networks. *Nat. Commun.* **2017**, *8*, 13890.

(23) Ruza, J.; Wang, W.; Schwalbe-Koda, D.; Axelrod, S.; Harris, W. H.; Gómez-Bombarelli, R. Temperature-transferable coarse-graining of ionic liquids with dual graph convolutional neural networks. *The Journal of Chemical Physics* **2020**, *153*, 164501.

(24) Wang, J.; Olsson, S.; Wehmeyer, C.; Pérez, A.; Charron, N. E.; de Fabritiis, G.; Noé, F.; Clementi, C. Machine Learning of Coarse-Grained Molecular Dynamics Force Fields. *ACS Cent. Sci.* **2019**, *5*, 755-767.

(25) Köhler, J.; Chen, Y.; Krämer, A.; Clementi, C.; Noé, F. Flow-Matching: Efficient Coarse-Graining of Molecular Dynamics without Forces. *J. Chem. Theory Comput.* **2023**, *19*, 942-952.

(26) Majewski, M.; Pérez, A.; Thölke, P.; Doerr, S.; Charron, N. E.; Giorgino, T.; Husic, B. E.; Clementi, C.; Noé, F.; De Fabritiis, G. Machine learning coarse-grained potentials of protein thermodynamics. *Nat. Commun.* **2023**, *14*, 5739.

(27) Wellawatte, G. P.; Hocky, G. M.; White, A. D. Neural potentials of proteins extrapolate beyond training data. *J. Chem. Phys.* **2023**, *159*, 085103.

(28) Charron, N. E.; Musil, F.; Guljas, A.; Chen, Y.; Bonneau, K.; Pasos-Trejo, A. S.; Venturin, J.; Gusew, D.; Zaporozhets, I.; Krämer, A. Navigating protein
landscapes with a machine-learned transferable coarse-grained model. *arXiv preprint arXiv:2310.18278* **2023**.

(29) Durumeric, A. E. P.; Charron, N. E.; Templeton, C.; Musil, F.; Bonneau, K.; Pasos-Trejo, A. S.; Chen, Y.; Kelkar, A.; Noé, F.; Clementi, C. Machine learned coarse-grained protein force-fields: Are we there yet? *Curr. Opin. Struct. Biol.* **2023**, *79*, 102533.

(30) Sahrmann, P. G.; Voth, G. A. Understanding the coarse-grained free energy landscape of phospholipids and their phase separation. *Biophys. J.* **2025**, *124*, 620-636.

(31) Ugarte La Torre, D.; Takada, S. Coarse-grained implicit solvent lipid force field with a compatible resolution to the Cα protein representation. *The Journal of Chemical Physics* **2020**, *153*, 205101.

(32) Izvekov, S.; Voth, G. A. A Multiscale Coarse-Graining Method for Biomolecular Systems. *J. Phys. Chem. B* **2005**, *109*, 2469-2473.

(33) Pelaez, R. P.; Simeon, G.; Galvelis, R.; Mirarchi, A.; Eastman, P.; Doerr, S.; Thölke, P.; Markland, T. E.; De Fabritiis, G. TorchMD-Net 2.0: Fast Neural Network Potentials for Molecular Simulations. *J. Chem. Theory Comput.* **2024**, *20*, 4076-4087.

(34) Schütt, K.; Kindermans, P.-J.; Sauceda Felix, H. E.; Chmiela, S.; Tkatchenko, A.; Müller, K.-R. Schnet: A continuous-filter convolutional neural network for modeling quantum interactions. *Advances in neural information processing systems* **2017**, *30*.

(35) Fu, X.; Wu, Z.; Wang, W.; Xie, T.; Keten, S.; Gomez-Bombarelli, R.; Jaakkola, T. *Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations*; 2022. DOI: 10.48550/arXiv.2210.07237.

(36) Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; Dávila-Contreras, E. M.; Qi, Y.; Lee, J.; Monje-Galvan, V.; Venable, R. M.; et al. CHARMM-GUI Membrane Builder toward realistic biological membrane simulations. *J. Comput. Chem.* **2014**, *35*, 1997-2004.

(37) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. *J. Chem. Phys.* **1983**, *79*, 926-935.

(38) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **2015**, *1-2*, 19-25.

(39) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O'Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. *J. Phys. Chem. B* **2010**, *114*, 7830-7843.

(40) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Jr. Optimization of the Additive CHARMM All-Atom Protein Force Field
Targeting Improved Sampling of the Backbone φ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. *J. Chem. Theory Comput.* **2012**, *8*, 3257-3273.

(41) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. *J. Appl. Phys.* **1981**, *52*, 7182-7190.

(42) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. *Mol. Phys.* **1984**, *52*, 255-268.

(43) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Physical Review A* **1985**, *31*, 1695-1697.

(44) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. *J. Chem. Phys.* **1993**, *98*, 10089-10092.

(45) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. *J. Comput. Chem.* **1997**, *18*, 1463-1472.

(46) Doerr, S.; Majewski, M.; Pérez, A.; Krämer, A.; Clementi, C.; Noe, F.; Giorgino, T.; De Fabritiis, G. TorchMD: A Deep Learning Framework for Molecular Simulations. *J. Chem. Theory Comput.* **2021**, *17*, 2355-2363.

(47) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. *J. Mol. Graphics* **1996**, *14*, 33-38.

(48) Khelashvili, G.; Kollmitzer, B.; Heftberger, P.; Pabst, G.; Harries, D. Calculating the Bending Modulus for Multicomponent Lipid Membranes in Different Thermodynamic Phases. *J. Chem. Theory Comput.* **2013**, *9*, 3866-3871.

(49) Doktorova, M.; Harries, D.; Khelashvili, G. Determination of bending rigidity and tilt modulus of lipid membranes from real-space fluctuation analysis of molecular dynamics simulations. *Physical Chemistry Chemical Physics* **2017**, *19*, 16806-16818, 10.1039/C7CP01921A.

(50) Sahrmann, P. G.; Voth, G. A. Enhancing the Assembly Properties of Bottom-Up Coarse-Grained Phospholipids. *J. Chem. Theory Comput.* **2024**, *20*, 10235-10246.

(51) Pan, J.; Tristram-Nagle, S.; Kučerka, N.; Nagle, J. F. Temperature Dependence of Structure, Bending Rigidity, and Bilayer Interactions of Dioleoylphosphatidylcholine Bilayers. *Biophys. J.* **2008**, *94*, 117-124.

(52) Lindblom, G.; Orädd, G. Lipid lateral diffusion and membrane heterogeneity. *Biochimica et Biophysica Acta (BBA) - Biomembranes* **2009**, *1788*, 234-244.

(53) Venable, R. M.; Ingólfsson, H. I.; Lerner, M. G.; Perrin, B. S., Jr.; Camley, B. A.; Marrink, S. J.; Brown, F. L. H.; Pastor, R. W. Lipid and Peptide Diffusion in Bilayers: The Saffman–Delbrück Model and Periodic Boundary Conditions. *J. Phys. Chem. B* **2017**, *121*, 3443-3457.

(54) Yang, S.-T.; Kiessling, V.; Tamm, L. K. Line tension at lipid phase boundaries as driving force for HIV fusion peptide-mediated fusion. *Nat. Commun.* **2016**, *7*, 11401.

(55) Shahidullah, K.; London, E. Effect of Lipid Composition on the Topography of Membrane-Associated Hydrophobic Helices: Stabilization of Transmembrane Topography by Anionic Lipids. *J. Mol. Biol.* **2008**, *379*, 704-718.
(56) Pan, J.; Tristram-Nagle, S.; Kucerka, N.; Nagle, J. F. Temperature dependence of structure, bending rigidity, and bilayer interactions of dioleoylphosphatidylcholine bilayers. *Biophys J* **2008**, *94*, 117-124.

(57) Garvey, C. J.; Bryant, S. J.; Elbourne, A.; Hunt, T.; Kent, B.; Kreuzer, M.; Strobl, M.; Steitz, R.; Bryant, G. Phase separation in a ternary DPPC/DOPC/POPC system with reducing hydration. *J. Colloid Interface Sci.* **2023**, *638*, 719-732.

(58) Leonenko, Z. V.; Finot, E.; Ma, H.; Dahms, T. E.; Cramb, D. T. Investigation of temperature-induced phase transitions in DOPC and DPPC phospholipid bilayers using temperature-controlled scanning force microscopy. *Biophys J* **2004**, *86*, 3783-3793.

(59) Stevens, J. A.; Grünewald, F.; van Tilburg, P. A. M.; König, M.; Gilbert, B. R.; Brier, T. A.; Thornburg, Z. R.; Luthey-Schulten, Z.; Marrink, S. J. Molecular dynamics simulation of an entire cell. *Front. Chem.* **2023**, *Volume 11* - 2023, Perspective.

(60) Yu, A.; Pak, A. J.; He, P.; Monje-Galvan, V.; Casalino, L.; Gaieb, Z.; Dommer, A. C.; Amaro, R. E.; Voth, G. A. A multiscale coarse-grained model of the SARS-CoV-2 virion. *Biophys. J.* **2021**, *120*, 1097-1104.

(61) Hudait, A.; Voth, G. A. HIV-1 capsid shape, orientation, and entropic elasticity regulate translocation into the nuclear pore complex. *Proceedings of the National Academy of Sciences* **2024**, *121*, e2313737121.

(62) Yu, M.; Heidari, M.; Mikhaleva, S.; Tan, P. S.; Mingu, S.; Ruan, H.; Reinkemeier, C. D.; Obarska-Kosinska, A.; Siggel, M.; Beck, M.; et al. Visualizing the disordered nuclear transport machinery in situ. *Nature* **2023**, *617*, 162-169.

(63) Jin, J.; Pak, A. J.; Voth, G. A. Understanding Missing Entropy in Coarse-Grained Systems: Addressing Issues of Representability and Transferability. *The Journal of Physical Chemistry Letters* **2019**, *10*, 4549-4557.

(64) Lu, L.; Voth, G. A. The multiscale coarse-graining method. VII. Free energy decomposition of coarse-grained effective potentials. *The Journal of Chemical Physics* **2011**, *134*, 224107.