Martini Coarse-Grained Nitrogen Gas Model for Lipid Nanobubble Simulations
Lipid nanobubbles have shown a great potential to be used for ultrasound molecular imaging and as biocompatible drug and gene delivery carriers, which integrate the advantages of both the biocompatibility of lipids and potent physicochemical properties of nanobubbles. Molecular dynamics (MD) simulation provides a powerful tool to investigate fundamental scientific problems related to lipid nanobubbles. With coarse-grained (CG) models, the system can be simulated with longer time scale and larger length scale. However, there are very few CG gas models for lipid nanobubble simulations. Hence, in this work, we developed a simple CG nitrogen gas model for the Martini force field by adjusting the Lennard-Jones interactions of 𝑁 with itself, water, and lipids. Our model can reproduce the density of pure 𝑁 , the density of 𝑁 within nanobubbles, as well as the partitioning thermodynamics of 𝑁 in DPPC bilayers. Further lipid nanobubble selfassembly simulations validate the reliability of our CG 𝑁 parameters. Besides, using threecomponent lipid nanobubbles, our simulations indicate that the nanobubble monolayers could also have clear liquid-liquid phase separation, the degree of which was proportional to the radius (𝑟) of the lipid nanobubble, and reached the maximum when 𝑟 → ∞ (i.e. a planar monolayer at the gaswater interface). By comparing the planar lipid monolayer and lipid bilayer systems, we furthermore find that the latter shows much less obvious phase separation, pointing to a thus far unexplored inter-leaflet coupling effect. In short, our CG MD simulations of lipid nanobubbles, monolayers and bilayers, in the presence of explicit gas molecules, provides useful insights into the role of membrane curvature and interleaflet coupling in the phase separation of multicomponent lipid membranes, and opens the way to more realistic simulations of gas/water or gas/lipid interfaces in general.
Introduction<!>Model and Methods<!>Results and Discussion<!>G1-G1 G1-P4 G1-Q0 G1-Qa G1-Na G1-C1 G1-C4 G1-SP1 G1-SC1 G1-SC3<!>Lipid Nanobubble and Monolayer Simulations Indicate the Critical Role of Membrane<!>Lipid Monolayer and Bilayer Simulations Demonstrate Inter-leaflet Couplings Modulate the Membrane Phase Separation.<!>Conclusions
<p>Nanobubbles have many unique physicochemical properties and are reported to have better stability than microbubbles [1][2] , and show great potential in a series of biomedical applications including ultrasound molecular imaging [3][4] , drug/gene delivery [5][6] , water treatment [7][8] , sonoimmunotherapy [9] , and so on. Amphiphilic molecules such as lipids can self-assemble along the air-water interface of nanobubbles. Similar to roles of pulmonary surfactant [10] , these molecules can greatly reduce the surface tension of the interface and increase the stability of nanobubbles. As the most abundant component of cell membranes, lipids play vital roles in maintaining its proper structure and functions. Hence, lipid nanobubbles may have excellent biocompatibility. On the other hand, lipid nanobubbles can fuse with cell membranes under a certain intensity of ultrasound [11] . Thus, cell membranes can incorporate lipids from lipid nanobubbles [11] and encapsulated gas molecules can re-distribute into the hydrophobic region of the cell membrane [12] . The former will change the local membrane composition, while the latter can decouple the two membrane leaflets and modify membrane structural properties. These two aspects may both dramatically affect the dynamics and functions of membrane proteins, which is so far largely unexplored. In other words, there is plenty of room in revealing the effects of lipid nanobubbles on structure and dynamics of cell membranes on the molecular level, which may greatly expand current biomedical applications of lipid nanobubbles.</p><p>Molecular dynamics (MD) simulations provide a powerful tool to investigate the interactions between biomolecules at atomic/near-atomic resolution [13] . Many computational efforts have been made to study the stability and dynamics of lipid nanobubbles as well as their interactions with model cell membranes [11,[14][15][16] . Especially, coarse-grained (CG) models (e.g. Martini [17][18][19] ) allow MD simulations with much larger length scale and longer time scale than all-atom models, and thus can better capture the complexity of realistic lipid nanobubbles and cell membranes. However, in the current state-of-art CG MD simulations, the gas phase is modeled as a vacuum in the isothermal-isochoric (NVT) ensemble [11,[20][21][22] . This will greatly hinder precise studies on the interactions between lipid nanobubbles and model cell membranes. On the one hand, in the NVT ensemble, the lateral dimension is fixed, which makes the description of large membrane deformations unreasonable. On the other hand, gas molecules can permeate into the hydrophobic region of a lipid bilayer [12] , therefore the absence of explicit gas molecules could leave possible potent physical effects of gas molecules on membrane lipids and proteins unclear. Hence, it will be critical to develop CG models of gas molecules for more precise description of the molecular-level interactions between lipid nanobubbles and biological systems (Fig. 1).</p><p>In this work, based on previous available results of 𝑁 gas from atomistic MD simulations [12,23] and experiments [24][25] , we developed a Martini-compatible CG model of 𝑁 gas (Fig. 1), which can well reproduce the lipid nanobubble self-assembly process and, at the same time, overcome possible artifacts of lipid monolayer (e.g. lipid nanobubble [11] , pulmonary surfactant [22] , tear film [26] ) simulations due to the vacuum treatment of gas molecules. With this model, we performed CG MD simulations of lipid nanobubbles with three different lipid components, 1,2-dipalmitoyl-snglycero-3-phosphocholine (DPPC), 1,2-dilinoleoyl-sn-glycero-3-phosphocholine (DUPC), and cholesterol, and found that the expected liquid-liquid phase separation could appear in the nanobubble monolayer. Besides, by comparing with planar lipid monolayer and lipid bilayer systems, our results indicated that membrane curvature and interleaflet coupling both have a strong impact on the phase separation.</p><!><p>In this work, the parametrization of 𝑁 gas within the CG Martini force field is mainly based on reproducing 𝑁 gas density (1.091 g/L) [24] at body temperature (~310K), as well as previous atomistic simulations of 𝑁 nanobubbles [23] and 𝑁 -lipid membrane interaction systems [12] . The obtained CG 𝑁 gas model was further validated in lipid nanobubble self-assembly simulations, which needs to be consistent with the experiments [4,27] . Then, CG MD simulations of lipid nanobubble, lipid bi-monolayer and lipid bilayer systems with the optimized 𝑁 gas parameters were performed to reveal effects of membrane curvature and inter-leaflet coupling on the membrane phase separation.</p><p>Martini Force Field. As a popular CG model, the Martini force field (version 2.1) [17][18] was used in the current work. In this model, generally 4 heavy atoms are mapped into 1 interaction site, including four main types: polar (P), nonpolar (N), apolar (C), and charged (Q). For fragments of molecules that cannot be adequately mapped with a 4 to 1 resolution (e.g., ring-like molecules), smaller, S-type particles are used. Four different subtypes (d=donor, a=acceptor, da=both, 0=none) are introduced to bead types of N and Q to mimic hydrogen bonding capacities and allow fine representation of the chemical nature. For bead types of P and C, five different subtypes (from 1, low polarity to 5, high polarity) are used to describe the degree of polarity. Ten levels of nonbonded interactions using different parameters ( ε , σ ) for Lennard-Jones (LJ) 12-6 potential energy function ( 𝑈 4𝜀</p><p>) are employed to describe the differential nonbonded interactions between different beads (Table S1). In addition, a Coulomb potential is used to describe electrostatic interactions between charged beads. For Martini DPPC/DUPC lipids [28] , the head-group consists of two charged bead types (Q0, +1e and Qa, -1e), the glycerol ester backbone is represented by two nonpolar bead types (Na), and each tail contains four apolar bead types (C1 for the saturated parts, C4 for the unsaturated parts). Five different bead types (SP1, SC3, SC1, SC1, C1) make up the cholesterol (Chol) molecules [29] . A polar bead type (P4) represents one CG water bead, and corresponds to four water molecules. In this work, one new CG bead type (G1) is introduced to represent one 𝑁 gas molecule. Systematic parameterization of the nonbonded interactions between G1 and the other beads (lipids and water) was performed to reproduce the key physiochemical properties mentioned above. Molecular Dynamics Simulations. The CG MD simulations of all systems were performed using the GROMACS program v2016.5 [30] and Martini force field [17][18] , while the visualization of system snapshots was done using VMD [31] . For all simulations, periodic boundary conditions were applied in three dimensions. The v-rescale thermostat [32] with a relaxation time τ = 1 ps was used to maintain a constant temperature of 310 K and a constant pressure of 1 bar was kept by Parrinello-Rahman pressure coupling [33] (coupling constant of 5 ps and compressibility is 5×10 -5 bar -1 ) in the NPT ensemble. For systems of pure 𝑁 gas, 𝑁 -water and lipid nanobubbles, an isotropic pressure coupling method was used, whereas semi-isotropic pressure coupling was applied to planar lipid bi-monolayer and lipid bilayer systems. A standard 1.2 nm cutoff was applied for van der Waals interactions, where the LJ potential was shifted to zero smoothly from 0.9 to 1.2 nm to reduce the cutoff noise. For the Coulumbic potential, a 1.2 nm cutoff was used for short-range electrostatic interactions while shifting to zero from 0 to 1.2 nm smoothly. The neighbor list for nonbonded interactions was updated every 10 steps with a cut-off of 1.2 nm. A 20 fs timestep was used to integrate the equations of motion.</p><p>System Setup. In order to parameterize the self-interaction parameters between 𝑁 molecules, 1728 G1 (𝑁 ) beads were evenly placed in a box of 6.0×6.0×6.0 nm 3 for all simulations, where 9×25 combinations of LJ interaction parameters 𝜎 ∈[0.28, 0.44] (nm) (∆𝜎=0.02 nm) and 𝜀 ∈[0.2, 2.6] (kJ/mol) (∆𝜀=0.1 kJ/mol) were investigated. Each simulation was run for 50 ns. For the crossinteraction parameters between 𝑁 and water molecules, 264 independent simulations of 2520 (box size: 6.9×6.9×6.9 nm 3 ), 5040 (box size: 8.5×8.5×8.5 nm 3 ), 7560 (box size: 9.7×9.7×9.7 nm 3 ) or 10080 (box size: 10.7×10.7×10.7 nm 3 ) molecules and different scaling factors (Table S2) and 𝑁 /water ratios (Fig. S1) using initial LJ parameters from SAFT-γ CG force field [34] were performed, which indicated that a proper scaling factor or different 𝜀 will be necessary to reproduce the gas bubble formation within Martini CG force field. Thus, we further performed systematic benchmark simulations with 13 different G1-G1 and 10 different G1 (𝑁 )-P4 (water) parameter sets (Table S3) to reproduce the minimal critical 𝑁 concentration required for gas bubble formation [23] . For each LJ parameter set, a 50 ns simulation of 3375 molecules with gradually increased 𝑁 /water ratio (box size: 7.0×7.0×7.0 nm 3 , Table S4) was run. In order to optimize LJ parameters between 𝑁 and lipids, a series of benchmark simulations with different parameters (𝜀, Table S5) were performed. Each system was simulated for 300 ns and consists of 128 DPPC, 2576 water and varying amount of 𝑁 molecules, with the initial box of 6.4×6.4×6.4 nm 3 . Four cases of 20, 200, 300 or 400 𝑁 molecules were considered to directly compare with the atomistic simulations by Zhang et al. [23] To further validate the LJ parameters obtained from the above benchmark simulations, the lipid nanobubble self-assembly process was studied. Different 𝑁 bubble size and density as well as number of DPPC molecules were considered (Table S6). For each simulation system, the 𝑁 bubble was placed at the center of the box, and DPPC together with water molecules were evenly distributed in the remaining space of the box. The initial simulation box is about 12.0×12.0×12.0 nm 3 , and each simulation was run for 3 µs to capture the whole self-assembly process. To investigate phase separation of lipid nanobubbles, systems with a lipid ratio of 5:3:2 (DPPC:DUPC:Chol) and different sizes were studied (Table S7). In the initial system, the lipid nanobubble was placed at the center of the box, which was 1 nm away from the edge of the box. The initial lipid distribution of the lipid nanobubble is approximately uniform. Each simulation was run for 3 µs to quantify the degree of phase separation. As control, we also simulated planar lipid monolayer and bilayer systems with comparable lipid composition to system N8 (Table S7-S9), which would provide insights into the roles of membrane curvature and interleaflet coupling in the phase separation of lipid membranes.</p><p>Trajectory Analysis. To quantify the degree of phase separation, we performed a contact analysis. To do so, we first obtained the total number of lateral contacts, 𝑁 , among unsaturated lipids DUPC in the phase separated systems; a contact was defined based on a distance cutoff of 0.6 nm between any two CG beads of the specified lipid type. Then we normalized 𝑁 by 𝑁 , the total number of lateral contacts in a pure bilayer of DUPC obtained from the last 200 ns trajectory of a 600 ns run. Hence, 𝑁 represents the maximum number of contacts between unsaturated lipids in a fluid bilayer at the same temperature. The normalized number of lipid contacts (𝑁 𝑁 /𝑁 ) was used as a proxy for quantifying relative domain size. Cholesterol preference was determined based on the number of contacts (cutoff 0.6 nm) of cholesterol with saturated (𝑁 ) and unsaturated lipids (𝑁 ) as 𝜒 𝑁 /𝑛 𝑁 /𝑛 𝑁 /𝑛 , 𝜒 𝑁 /𝑛 𝑁 /𝑛 𝑁 /𝑛 where 𝜒 and 𝜒 are the fraction of cholesterol in contact with saturated and unsaturated lipids, and 𝑛 and 𝑛 are the total number of CG beads of saturated and unsaturated lipids, respectively.</p><!><p>Parametrization of the Nitrogen Gas Model in the Martini Force Field. In order to determine the Martini-compatible parameters for 𝑁 gas molecules, we firstly calibrated self-interaction parameters between 𝑁 CG beads (G1) by reproducing the experimental 𝑁 gas density (1.091 g/L) at body temperature [24] . According to Cao's work on the parameterization of 𝑁 gas molecules within SAFT-γ CG force field [34] , the LJ interaction parameters between 𝑁 gas molecules are 𝜎=0.36 nm and 𝜀=0.7 kJ/mol. In order to evaluate the transferability from SAFT-γ 𝑁 parameters to Martini, the LJ interaction parameters 𝜎 ∈[0.28, 0.44] (nm) (∆𝜎=0.02 nm, 9 systems) and 𝜀 ∈ [0.2, 2.6] (kJ/mol) (∆𝜀 =0.1 kJ/mol, 25 systems) for nonbonded interactions between G1 beads are systematically benchmarked for the pure 𝑁 gas systems containing 1728 𝑁 gas molecules. Each simulation system was run for 50 ns at T=310 K, which added up to a total simulation time of 9×25×50 ns=11.25 µs. The density of 𝑁 gas molecules was evaluated over the last 20 ns trajectories as shown in Fig. 2a. Generally, stronger LJ interactions will induce larger 𝑁 gas density. The gas density with 𝜎 =0.36 nm and 𝜀 =0.7 kJ/mol, i.e. corresponding to the SAFT-γ parameters, is 1.087 g/L, and slightly underestimates the exact experimental 𝑁 gas density (Fig. 2a). Considering the value of energy parameter 𝜀 is especially critical to the nanobubble formation process, while the gas density appears less sensitive to the value of the distance parameter 𝜎, we finally chose 𝜎=0.36 nm and narrowed the range of 𝜀 ∈[0.7, 1.9] kJ/mol, which will be further determined in the following nanobubble formation benchmark simulations. The triangle and five-pointed star symbols represents the final choices of 𝑁 -𝑁 nonbonded interaction parameters for Cao's work [34] and our work. (b) Optimized parameters (𝜀 is 1.6 and 1.4 kJ/mol for G1-G1 and G1-P4 LJ interactions respectively) enable better reproduction of the critical concentration for gas nucleation and gas density within the nanobubble in Zhang et al.'s work [23] . (c) Optimized parameters (𝜀 is 1.4, 1.4, 1.7 and 2.0 kJ/mol for G1-Q0, G1-Qa, G1-Na and G1-C1 LJ interactions, respectively) reproduce the preferred localization of 𝑁 inside a DPPC lipid bilayer [12] . 𝑁 and DPPC headgroups are shown as red and blue spheres respectively, DPPC tails as green bonds, and water molecules as gray dots.</p><p>In order to obtain the nonbonded interaction parameters between coarse-grained 𝑁 beads (G1) and water beads (P4), the Lorentz-Berthelot (LB) combination rule (𝜎 𝜎 𝜎 /2, 𝜀 𝑇 1 𝑘 𝜀 𝑇 𝜀 𝑇 ; 𝑘 is the scaling factor) was used. With the standard LB combination rule (𝑘 =0, 𝜎 =0.42 nm and 𝜀 =1.9 kJ/mol for the nonbonded interactions between G1 and P4 beads), it is very difficult to reproduce the nanobubble formation process in water from the CG MD simulations (Fig. S1). However, this can be corrected by introducing a proper scaling factor 𝑘 for 𝜀 (Tables S2, Fig. S1). Hence, LJ parameters 𝜎=0.36 nm, 𝜀 ∈[0.7, 1.9] kJ/mol for G1-G1 and 𝜎=0.42 nm, 𝜀 ∈[1.1, 2.0] kJ/mol for G1-P4 (Table S3) were carefully calibrated to obtain the optimized energy parameters which can reproduce the critical 𝑁 gas concentration (22 kg/m 3 ) necessary for the nanobubble formation [23] . For each LJ parameter set (13×10=130 in total), 50 ns benchmark simulations (11 different initial gas concentrations from 10 kg/m 3 to 30 kg/m 3 , ∆=2 kg/m 3 , Table S4) were performed of a system composed of 𝑁 CG beads in CG water. On average, eight simulation systems of different initial gas concentrations were performed for each parameter set. In other words, the total simulation time for this part is 13×10×8×50 ns=52 µs. As shown in Table S3, LJ parameters of 𝜎=0.36 nm, 𝜀=1.6 kJ/mol for G1-G1 and 𝜎=0.42 nm, 𝜀=1.4 kJ/mol for G1-P4 enable the reproduction of the critical 𝑁 gas concentration for the nanobubble formation observed in all-atom MD simulations [23] . Using this parameter set, we could also obtain a comparable 𝑁 gas density within the nanobubble as reported in Zhang et al. [23] (Fig. 2b). Note that the density of N2 inside a nanobubble is much higher than the bulk gas density, in line with experimental findings [25]. It is worth mentioning that LJ parameters 𝜎=0.36 nm, 𝜀=1.6 kJ/mol for G1-G1 can also achieve a value of 𝑁 gas density (1.090 kg/m 3 , Fig. 2a) in close agreement to the experimental 𝑁 gas density at body temperature (1.091 kg/m 3 ). Hence, the interaction parameters for 𝑁 -𝑁 (G1-G1) and 𝑁 -water (G1-P4) are fixed.</p><p>In order to expand the usage of the Martini 𝑁 CG model to lipid systems, we further determined the LJ interaction parameters between the 𝑁 bead (G1) and the lipid beads. The benchmark simulations are mainly based on the united-atom MD simulations of 𝑁 -DPPC lipid bilayer systems [12] . In this work, Li et al. [12] quantitatively described the partitioning thermodynamics of 𝑁 gas molecules into the hydrophobic regions of a DPPC lipid bilayer. For our benchmark simulation systems, 𝑁 gas molecules (G1 bead) were initially evenly distributed in the bulk water (P4 bead type) region of a DPPC (comprising Q0, Qa, Na, and C1 bead types) bilayer system. According to the Martini interaction parameters (Table S1), Q0 and Qa share almost the same 𝜀 as P4 interacting with other beads. Hence, the 𝜀 values for G1-Q0 and G1-Qa nonbonded interactions are also fixed as 1.4 kJ/mol. Based on the LB combination rule, the corresponding 𝜎 values for G1-Q0 and G1-Qa are 0.49 nm, while for G1-C1 and G1-Na 0.42nm (the effective size of 𝜎 0.47 nm is assumed, except for interaction between charged and most apolar beads which share 𝜎 0.62 nm in the Martini model). According to the nonbonded interaction rules for beads with different polarities in Martini (Table S1), 𝜀 ∈[1.6, 2.0] kJ/mol for G1-C1 and 𝜀 ∈[1.4, 2.0] kJ/mol for G1-Na LJ interactions (Table S5) were finally evaluated in order to reproduce the distribution of 𝑁 gas molecules in the DPPC bilayer systems as indicated by the reference simulations [12] . The energy parameter 𝜀 of the G1-C1 nonbonded interaction should always be larger than that of the G1-Na nonbonded interaction, to reflect the more polar nature of the Na compared to the C1 bead type. Therefore, a total of 42 different sets of 𝜀 values were tested for G1-C1 and G1-Na interactions (Table S5). For each parameter set, four lipid bilayer systems with 𝑁 gas molecules of different numbers dispersed in bulk water were performed. The simulation time of each benchmark system was 300 ns, with the last 100 ns used for the analysis of density profiles. Hence, the total simulation time for this part is 42×4×300 ns=50.4 µs. We finally chose 𝜀 =2.0 kJ/mol for G1-C1 and 𝜀 =1.7 kJ/mol for G1-Na, which properly reproduced the preferred localization of 𝑁 gas molecules (Fig. 2c, S2) [12] . Besides, we further obtained the nonbonded interactions parameters between 𝑁 gas molecules (G1 bead) and DUPC, cholesterol (C4, SP1, SC1 and SC3 beads) based on the above benchmark simulations and parameterization rules for nonbonded interactions in Martini model (Table 1).   Validation with Lipid Nanobubble Self-Assembly Simulations. As we know, the lung alveoli provides the key air-water interface for the gas exchange process. The synthetic lipid lamellar bodies and surfactant proteins can be transported to the interface and re-assemble into the lipid monolayer structure [10,22,35] . Hence, from the perspective of bionics, the generation of gas bubbles in lipid solution can be one simple experimental method to synthesize lipid micro/nano-bubbles [4,[36][37] . In other words, the feasibility of the above-obtained 𝑁 gas CG model can be further evaluated by lipid nanobubble self-assembly simulations. In order for this, we set up the initial simulation systems as follows: (1) an 𝑁 gas nanobubble, with varying radius and number of molecules, was placed at the center of the simulation box (12×12×12 nm 3 ); (2) 180 or 256 DPPC molecules were evenly dispersed around; (3) the whole system was then solvated by water molecules. In total, 12 self-assembly simulation systems (Table S6) were set up. During the preequilibrium stage of each simulation, the position of 𝑁 gas molecules was constrained, while DPPC and water molecules were fully relaxed. Then, the simulation went through a subsequent 3 µs production run. For amphiphilic molecules such as DPPC lipids, it is easy to form micelles in aqueous solution [38] , which may convert into spherical vesicles in the absence of gas bubbles [39][40] .</p><!><p>However, in the presence of a gas nanobubble, dispersed lipids or small lipid clusters could reassemble along the gas-water interface with hydrophobic lipid tail toward the gas core, and thus form a lipid nanobubble. Fig. 3 shows a typical assembly process for the lipid nanobubble formation in our simulations. First, within a few ns, small lipid clusters and dispersed DPPC molecules close to the 𝑁 nanobubble quickly spread out along the nanobubble surface forming a lipid nanobubble. Lipids far from the 𝑁 nanobubble initially self-assembled into a micelle capturing a few free 𝑁 molecules. Eventually, after about one or two hundreds of ns, the lipid micelle fused with the lipid nanobubble. Finally, the self-assembled DPPC nanobubble relaxed to a stable spherical structure in bulk water for the remaining time of the 3 µs CG MD simulations. It is worth mentioning that insufficient or excessive DPPC molecules can lead to semi-coated lipid nanobubbles or tubular shaped ones, respectively (Table S6), which is consistent with the principles for self-assembly of lipids [38] . In other words, we are capable of properly simulating the self-assembly process of the lipid nanobubble formation, which further validates the application feasibility of our 𝑁 gas CG model in lipid-related biological systems.  S7). All lipid tails are shown as green bonds, head-groups of DPPC, DUPC and Chol as blue, yellow and white spheres correspondingly, and N 2 as red spheres. Water molecules are not shown for clarity. Scale bar: 5 nm. (b) Normalized lateral contacts of unsaturated lipids (DUPC) for system N1-N8. (c) Cholesterol preference to different lipids for system N1-N8. Error bars are standard deviations derived from the statistics over four 500 ns blocks of the last 2 µs.</p><!><p>Curvature in Membrane Phase Separation. The generation of membrane curvature, which can be mediated by lipids, proteins or an external stimulus, is necessary in a series of important biological processes for the proper functions of cells such as protein sorting, membrane fusion, organelle shaping and enzyme activation. [41] The presence of local membrane curvature could induce lipid sorting or redistribution [42][43][44][45] , which would change the local membrane properties to fulfill certain functions. On the other hand, lipid domains, which are driven by liquid-liquid phase separation of membrane lipids and proteins, also play critical roles in membrane-related biological processes. However, whether and how membrane curvature modulates the size and stability of lipid domains is still unclear. As discussed above, the spherical lipid nanobubble, which has only one lipid leaflet, may serve as an ideal model system for this purpose. As is known, for the phase separated lipid bilayer system, membrane domains have both intra-leaflet [46][47][48][49] and inter-leaflet [50][51][52][53][54][55][56] dynamics. Usually, these two kinds of membrane domain dynamics are closely related to each other. In other words, our lipid nanobubbles allow the investigation of effects of membrane curvature on intraleaflet membrane domain dynamics exclusively. Hence, we set up eight three-component lipid nanobubble systems (system N1-N8) containing DPPC, DUPC and CHOL molecules with the initial inner nanobubble radius ranging from 4 nm to 11 nm. The system details can be found in Table S7. As shown in Fig. 4a, clear membrane phase separation appeared in our lipid nanobubbles, with formation of liquid-disordered (𝐿 ) domains enriched in DUPC and liquid-ordered (𝐿 ) domains enriched in DPPC and cholesterol, in line with experimental data for this three-component system as well as with previous Martini-based CG MD simulations of ternary bilayer systems [48,57] . Larger lipid nanobubbles (smaller membrane curvature) showed more obvious membrane phase separation, which was further validated by quantifying the normalized lateral contacts of unsaturated lipids (Fig. 4b) as well as cholesterol preferences (Fig. 4c). We rationalize the curvature effect by acknowledging the larger bending rigidity of 𝐿 domains [58] , which makes it energetically unfavorable to form such domains on highly curved surfaces.  S8). All lipid tails are shown as green bonds, head-groups of DPPC, DUPC and Chol as blue, yellow and white spheres correspondingly, N2 as red spheres, water molecules as gray dots. Scale bar: 5 nm. (c) Normalized lateral contacts of unsaturated lipids for system N8 and M8. (d) Cholesterol preference to DPPC, DUPC and their differences for system N8 and M8. Error bars are standard deviations derived from the statistics over four 500 ns blocks of the last 2 µs.</p><p>As the radius of the lipid nanobubble increases, MD simulations need much more computational resources. Hence, the largest lipid nanobubble we investigated in this work was limited to an inner radius of 10 nm. In order to overcome this shortcoming, we set up planar lipid bi-monolayer systems (Fig. 5a-b) to simulate the case with the radius → ∞. System M8 (Table S8) has the same number of lipids and area per lipid in each monolayer as that of system N8. As shown in Fig. 5cd, the results clearly indicated that system M8 had much stronger membrane phase separation than that of system N8, which further validated the point that the degree of membrane phase separation decreased with the degree of membrane curvature (reciprocal of the radius of curvature). Note that 𝑁 gas molecules were introduced into our lipid bi-monolayer setup to directly model the gas phase, which was thus far modeled by a vacuum slab in previous state-of-art CG MD simulations of monolayer systems [22,[59][60][61][62] . When the vacuum is used, the z-compressibility of the simulation box has to be zero. With the real gas molecules in the simulations, this setup is no longer necessary. Besides, the explicit interactions between gas molecules and lipids may facilitate the better description of the dynamics of interfacial lung surfactant lipids during the exhalation (monolayer compression) and inhalation (monolayer expansion) processes.  S9). Lipid tails are colored with green bonds, while lipid headgroups are depicted as follows: DPPC (blue), DUPC (yellow), and CHOL (white). Water molecules are not shown for clarity. Scale bar: 5nm. (b) Normalized lateral contacts of unsaturated lipids for system B8 and M8. (c) Cholesterol preference to DPPC, DUPC and their differences for system B8 and M8. Error bars are standard deviations derived from the statistics over four 500 ns blocks of the last 2 µs.</p><!><p>As is discussed above, in the absence of inter-leaflet couplings, obvious membrane phase separation appears in both lipid nanobubble and the monolayer systems. Previous studies have indicated that inter-leaflet couplings, which can be regulated by many physicochemical factors [56] , modulate membrane domain registration/anti-registration dynamics [50-52]   . Whether inter-leaflet couplings can affect the intra-leaflet membrane domain dynamics is still not clear. Here, lipid bi-monolayer and bilayer systems with the same lipid number and comparable area per lipid will allow us to resolve this question directly. As shown in Fig. 6, compared to the lipid bilayer system (system B8, Table S9), the lipid monolayer system (system M8, Table S8) shows the stronger membrane phase separation. In other words, the inter-leaflet couplings appear to negatively impact the intra-leaflet membrane domain separation dynamics in our simulation systems. A possible explanation of this observation is the entropic cost associated with the (anti)alignment of the domains across the two leaflets in the bilayer systems, or enthalpic cost in case of non-alignment.</p><!><p>In this work, based on the density of 𝑁 gas molecules [24] as well as previous atomistic MD simulations of 𝑁 gas nanobubbles [23] and 𝑁 -lipid bilayer systems [12] , we obtained a Martinicompatible coarse-grained 𝑁 gas model through a series of benchmark simulations, which defined the detailed LJ nonbonded interaction parameters between 𝑁 and other molecules (Table 1). The subsequent lipid nanobubble self-assembly simulations further validated the feasibility of our coarse-grained 𝑁 gas model. Hence, we applied this model to study the lipid dynamics of three-component lipid nanobubbles. The results indicate that the mixture of saturated and unsaturated lipids as well as cholesterol can undergo clear 𝐿 -𝐿 phase separation in lipid nanobubbles. The larger the lipid nanobubble size, the more obvious the membrane phase separation will be. The lipid monolayer simulations further confirmed the role of membrane curvature in the phase separation of a single membrane leaflet. It also showed the feasibility of our gas model for pulmonary surfactant simulations, which may overcome the possible artefacts of using vacuum for the gas phase in previous coarse-grained MD simulations [22,[59][60][61][62] . Besides, by comparing lipid monolayer and bilayer simulations, our results indicated that inter-leaflet coupling could hinder the intra-leaflet membrane domain dynamics. In short, our work shows that the developed coarse-grained 𝑁 gas model can be well applied to simulate lipid nanobubbles with complex lipid compositions and in the presence of an explicit gas-water interface, which are essential for studying the interaction molecular mechanism between lipid nanobubbles and biological systems. However, further efforts are needed to extend the parameterization of 𝑁 toward its interaction with other bead types of the Martini force field, and to develop models of other common gas molecules for lipid nanobubbles. In light of the recently released new version of Martini (v3.0) [63] , it would be worthwhile exploring whether one of the new particle types could be used to represent 𝑁 , with the added benefit that all cross-interactions with other bead types are already validated. Based on the current model, we and others may start exploring the interactions between lipid nanobubbles (which may include proteins and encapsulated drug molecules) and biological systems, which are essential for promoting the wider biomedical applications of lipid nanobubbles. Table S3. Critical concentration (g/L) for nitrogen nanobubble formation obtained from simulations of nitrogen-water mixtures with different LJ interaction parameter combinations.</p>