SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations with Single-site Models

Implicit solvent, coarse-grained models with pairwise interactions can access the largest length and time scales in molecular dynamics simulations, owing to the absence of computationally expensive interactions with a huge number of solvent particles, the smaller number of interaction sites in the model molecules, and the lack of fast sub-molecular degrees of freedom. Without an explicit solvent, the solvent mediated effects are mimicked entirely through the interactions between the model molecules. In this paper, we describe a maximally coarse-grained model for lipids in implicit water. The model is called SiMPLISTIC, which is an acronym for Single-site Model with Pairwise interaction for Lipids in Implicit Solvent with Tuneable Intrinsic Curvature. SiMPLISTIC lipids rapidly self-assemble into realistic non-lamellar and lamellar phases such as inverted micelles and bilayers, the spontaneous curvature of the phase being determined by a single free parameter of the model. Model membrane simulations with lamellar SiMPLISTIC lipids show satisfactory fluid and gel phases with no interdigitation or tilt. Despite being rigid molecules, SiMPLISTIC lipids can generate experimentally relevant values for the bending stiffness of model membrane bilayers with no significant interleaflet coupling. SiMPLISTIC can also simulate mixtures of lipids that differ in their packing parameter or length, the latter leading to the phenomenon of hydrophobic mismatch driven phase separation or domain formation. The model has a large scope due to its speed, conceptual and computational simplicity, and versatility. Applications may range from academic and industrial research in various lipid-based systems, such as lyotropic liquid crystals, biological and biomimetic membranes, vectors for drug and gene delivery, emulsions for cosmetic products, to education, such as interactive molecular dynamics simulations.

approach for a unified atom or a molecular level representation. The various approaches to CG modeling of biomolecular systems and solvated amphiphiles are well-discussed in the literature and the reader is referred to Ref. 22 and 23 for excellent overviews on the topic.
Of the different CG methodologies used for modeling lipids in water, a special class consists of implicit solvent, top-down approaches in which the lipids are modeled as single or multisite molecules interacting via ad hoc potentials that are designed for simplicity and speed of computation. 'Implicit solvent' a means water is never modeled explicitly, the solvent effects being incorporated within effective interactions between lipids. Any vacuum or unoccupied space in the simulation box, therefore, should be considered as water (Figure 1). These particle-based models are top-down in the sense that they are supposed to reproduce the mesoscopic and macroscopic properties such as self-assembly, membrane elasticity etc. without considering the microscopic details explicitly. An accessible review of such implicit solvent coarse-grained (ISCG) models can be found in Ref. 24.
The goal of the present paper is to present a novel ISCG model for self-assembling lipids, where the lipids are modeled as single-site directed ellipsoids interacting via a simple, two-body, anisotropic potential. Being molecular level (i.e. single-site) and implicit solvent, the number of site-site interactions is minimal and the MD timestep can be maximum, because of which this model can simulate lipid systems at the largest scales accessible by any particle-based model for given computational resources. The pairwise interaction potential features a tuneable parameter that determines the spontaneous curvature of the self-assembled phase. For ease of reference, we shall henceforth call this model 'SiMPLISTIC', which is an acronym for 'Single-site Model with Pairwise interaction for Lipids in Implicit Solvent with Tuneable Intrinsic Curvature'.
To put SiMPLISTIC in perspective, let us briefly discuss some of the previous attempts at minimalistic, b ISCG modeling of lipids for offlattice simulations. The pioneering model by Drouffe et al. 25 consisted of directed spheres, which, using a multibody hydrophobic potential, self-assembled into monolayers. To simulate a lipid bilayer with this monolayer, each sphere would be viewed as containing a hydrophobic layer sandwiched between two (top and bottom) hydrophilic layers. This model, therefore, represented a pair of lipids instead of a single a Some authors use the term 'Solvent-free' instead. b By 'minimalistic', we mean models with up to three interaction sites only. For other ISCG models see Ref. 48 lipid. Nogushi and Takasu 26 were the first to give an ISCG model of a single lipid proper. Their model consisted of a rigid, linear trimer of one hydrophilic and two hydrophobic beads. Self-assembly was still due to a multibody hydrophobic potential. Pairwise potentials, however, are far more computationally efficient than multibody ones. This led to the 'Water-free' model by Farago 27 where each lipid was still modeled as a rigid, linear trimer, but now all interactions were pairwise. Brannigan and Brown 28 took a different approach and brought nonspherical components into the picture. Their model consisted of a nematogenic spherocylinder with a hydrophobic site at its one end, all interactions being pairwise. The molecules in the four models mentioned above are essentially rigid linear rotors with only five degrees of freedom. The fluid bilayers produced by them, however, were either too flexible (forming vesicles easily) 25,26 or too stiff (comparable to membranes with cholesterols) 27,28 . To remedy this, Wang and Frenkel 29 introduced flexibility in the bead-based trimer model by using finite extensible nonlinear elastic (FENE) bonds to link the beads. By this time, it had already become clear for the workers in this field that simple Lennard-Jones (LJ) potentials are not suitable for stabilizing a fluid bilayer. To this effect, some had to use densitydependent multibody potentials 25,26,29 and some 28 employed anisotropic, inverse-squared attraction with a longer interaction range. Although Farago 27 was able to stabilize a fluid bilayer with LJ type interactions only, he had to make the pair potentials non-additive and remain limited to only a heavily tuned set of parameters. Farago's strategy required the hydrophobic beads to see the hydrophilic beads as larger than they actually were, thus creating a steric barrier against lipid evaporation from the membrane plane. As expected, this compromised self-assembly and only preassembled bilayers could be studied. Cooke, Kremer and Deserno 30,31 then realized the capability of broad attractive tail potentials to accomplish self-assembly into stable fluid bilayers. Their model 30 was again a flexible (FENE linked) trimer with the hydrophobic beads attracting each other pairwise and the hydrophilic beads providing soft-core, steric interaction only. The range of the pairwise attraction, however, was now a tuneable parameter that affects the thermal stability and the elasticity of the bilayer. For completeness, we must mention the relatively recent model by Noguchi 32 , where each lipid is modeled as a directed sphere with two interaction sites and selfassembly is ensured using a multibody hydrophobic potential. A notable commonality between the above-mentioned ISCG models with pairwise interactions is that none of them incorporates the hydration force 33 , which is a key feature of hydrated amphiphile aggregates. Hydration force accounts for the effective repulsion between two hydrophilic layers in close contact, and decays exponentially with the distance between them. This 'hydration barrier' is the reason, apart from thermal fluctuations such as protrusions and undulations, why membrane fusion is difficult 34 and why a bilayer cannot sit immediately on top of another in multilamellar stacks c (Figure 1a). In 2017, we published a single-site, ISCG model for amphiphiles, with a pairwise, anisotropic potential, that could mimic such a 'hydration barrier'. 35 Amphiphiles were modeled as soft-core directed spheroids that attracted each other, mimicking hydrophobic force, for some relative orientations and repelled each other, mimicking hydration force, for the rest. Because of this, the packing parameter of the molecule could become that of a (truncated) cone even though the steric core was ellipsoidal. The model featured two key parameters: a tuneable range parameter as in the Cooke-Kremer-Deserno model 30 , and a packing parameter tuner that would determine the aggregate morphology. Unlike the rigid anisotropic models of Brannigan and Brown 28 and Noguchi 32 , our single-site 2017 model required no exclusively hydrophobic or hydrophilic interaction site, and yet showed rapid, unassisted self-assembly into fluid bilayers, rods and micelles. As a further test of this model, we used it to build ISCG models of both rigid and flexible bolaamphiphiles 36 , which successfully reproduced many key features of hydrated bolaamphiphile systems in MD simulations.
c Adding water to an anhydrous multilamellar stack introduces a layer of water between consecutive bilayers which leads to 'swelling' 83 . Although this 2017 model 35 fulfilled our requirement of a computationally fast, conceptually and programmatically simple d ISCG model for generic amphiphiles with various packing parameters, it could not generate phases with negative curvatures. Hence, we could not use it to model the non-lamellar lipids. In fact, to the best of our knowledge, none of the above-mentioned minimalistic ISCG models had been shown to produce stable aggregates with negative intrinsic curvatures through spontaneous self-assembly. This is not surprising when we recognize that the formation of inverted morphologies by ISCG models is more involved than merely getting the molecular packing parameters right. To illustrate, an obvious difference between an inverted micelle with its negative curvature and a direct micelle with positive curvature is that the direct micelle does not encapsulate water whereas the inverted micelle does (Figure 1b and c). The hydration barrier due to this captured water should therefore be taken into account appropriately in an implicit solvent model and a simple model that mimics only the packing parameter of the lipids would not suffice e . This is where our new ISCG model SiMPLISTIC fits in. Here, the model particles can efficiently self-assemble into inverted micelles that encapsulate vacuum, which is to be interpreted as water because the model is implicit solvent. SiMPLISTIC can also form bilayers. It can therefore model both lamellar and non-lamellar lipids. In addition, for model membranes, SiMPLISTIC can generate elastic properties consistent with experiments and other simulations, unlike the other rigid models (see discussion above). Other important features implemented in SiMPLISTIC are the hydrophobic mismatch between lipids of unequal lengths for studying demixing and domain formation in multicomponent membranes, mixing rules for studying mixtures of lamellar and non-lamellar lipids, and the ability to form interdigitationfree bilayers regardless of the lipid length. The latter is relevant in view of the fact that phospholipids with symmetric acyl chain lengths do not show interdigitation in bilayers.
SiMPLISTIC, to the best of our knowledge, is the first, truly singlesite model for self-assembling lamellar and non-lamellar lipids using a pairwise potential. SiMPLISTIC's ability to model non-lamellar lipids by simply tuning a parameter, however, came as a by-product. Our original scheme was to design the most minimalistic, single-site lipid model for large-scale model membrane simulations with pairwise interactions only. The basic aim, therefore, was to be able to generate fluid, interdigitation-free bilayers through rapid, spontaneous selfassembly, while also taking into account the interbilayer hydration force. The next section (Sec. II) presents the design of SiMPLISTIC with this original purpose in mind. Sec. III then discusses how the ability to self-assemble into aggregates with different curvatures comes along because of our chosen pair potential. Finally, Sec. IV reports the results of model membrane simulations and extends SiMPLISTIC to include hydrophobic mismatch between lipids of different lengths. We conclude this paper in Sec. V with discussions of SiMPLISTIC's speed, scope and future directions.

A. Design Principles
Given the ad hoc nature of any top-down, ISCG model, it is necessary to lay down a set of design principles around which the model would be developed. If certain design principles lead to a successful model, then those principles must have captured some essence of the physics involved at the time and length scales appropriate to the model. The design of SiMPLISTIC is based on the following core principles. d The model potential was essentially a (switched) Gay-Berne potential 40 with a modified well-depth and a broad (cubic spline) tail. e For example, no non-lamellar, inverted phase such as the inverted hexagonal, micellar or cubic phase was formed by the Cooke-Deserno model, even when the lipid profile was made inverted conical by making the headgroup bead smaller than the hydrophobic beads 84 .

Absence of sub-molecular degrees of freedom
We choose a rigid model. This makes the model molecule translate and rotate as a whole with a timescale that is larger than the submolecular degrees of freedom of a flexible molecule would ever allow. A rigid model, therefore, can have the largest MD timestep. Other works in the literature have also found the use of rigid models fruitful 27,28,32,35,37 . Using 31 P-NMR data and MD simulations, Klauda et al. 38 established the rigid body modeling to be a good description of phospholipids in a fluid bilayer. Lipids could be viewed as rigid cylinders with two primary modes of rotationrotation about the long axis and wobble. According to Marsh 39 , rapid long axis rotation kicks in by the main transition. This motivates an axisymmetric model of lipids for our case. Because the lipids are not inversion symmetric, a rigid model of lipids must also have an inherent directedness that may be described by a unit vector û directed from the tail end to the head end (see inset of Figure 2). Such an axisymmetric and directed model has only five physically significant degrees of freedom, viz. the orientation of its unit vector and its centre of mass (c.m) coordinates. The remaining degree of freedom, viz. the rotation about its axis of symmetry is physically inconsequential.

Single-site model with pairwise interaction
A valid single-site model with a simple, anisotropic pair potential is both economic and elegant compared to an equally valid model with multiple interaction sites, if both the models are to be applied at the same length and time scales. The legacy of the Gay-Berne potential (GB) has been particularly influential in this context. GB 40 is a simple, single-site, anisotropic potential that was originally designed to substitute for a rigid tetramer of LJ sites, and has since seen remarkable popularity in modeling mesogenic molecules in liquid crystal research [41][42][43] .

Template for the pair potential
Let us take two model particles i and j with c.m coordinates Meanings of the variables used above are the following. zf : r r at which the potential is minimum and the particles exert zero force on each other; cut : r r at which the potential decays to zero or is cut-off; wd :  the well-depth of the attractive potential; ( ) : sr a spline curve that switches from -1 at zf r to 0 or near zero at cut r ; Steric ( ) : Vr A strictly repulsive potential describing the soft-core. Contextually, it might be noted that the popular LJ and GB potentials also comply with the format given in Eq. (1)- (3). As in GB, any anisotropy in the pair potential is introduced through the dependence of zf r and wd  on the relative orientation of the particles. zfˆˆ( , , ) ij i j r r u u , for example, would correspond to the nonspherical shapes of the molecular cores.
Desirably, the pair potential should be symmetric under swapping of the particles, viz.

Phenomenology
Being top-down, our ISCG modeling approach needs to be guided by some phenomenology of generic lipid bilayers. A preference for certain relative configurations between a pair of nearest neighbor lipids in a hydrated bilayer can be easily ascertained from the schematic in Figure 1a. The lipids labeled 1 and 2 are in the side-side parallel (SSP) configuration, whereas those labeled 1 and 3 are in the tail-tail collinear (TTC) configuration. In a stable bilayer, lipids 1, 2 and 3 are bound to each other. If all the lipids interact only through an implicit solvent pair potential as described in Sec. II.A.3, this implies that the SSP and TTC configurations must have negative pair potential energies or, equivalently, positive well-depths ( wd  ). Figure 1a, on the other hand, belong to adjacent bilayers. They, therefore, should repel each other in their head-head collinear (HHC) configuration so as to mimic the effective hydration force between those bilayers. The potential described in Eq. (3), therefore, needs to be repulsive for the HHC configuration, which is possible if HHC wd 0   . Figure 2 lists all possible side-side and collinear configurations between a pair of model lipids. For the sake of discussion, we shall call these the canonical configurations. Let us now discuss the remaining two canonical configurations that we did not discuss before. If a lipid in a bilayer leaflet is flipped upside down about its centre of mass, it will give rise to the side-side antiparallel (SSA) and the head-tail collinear (HTC) configurations as shown in Figure 3a. Because such a flip is prohibited, as otherwise it would expose the flipped lipid's tail to the surrounding water, the SSA and HTC configurations must be highly unstable. Lipid pairs in SSA or HTC configurations, therefore, must be unbound, which, like the HHC case above, requires SSA contributes to bilayer stability in other ways too. To illustrate, note that the SSA configuration also arises in case of extreme interleaflet interdigitation, viz. when the lipids from one leaflet are pushed into the other leaflet, thus thinning the bilayer down to the width of a monolayer ( Figure 3b). In order to avoid the unstable SSA configuration, therefore, the leaflets start to push each other away. This way, the model monolayers can produce a realistic resistance when forced to interdigitate. On the other hand, when a lipid protrudes out of the bilayer plane or tries to escape from the bilayer, it must pass through a configuration where its tail comes near the head of the neighboring lipids, all the lipids being parallel (Figure 3c). This configuration is very close to HTC and hence, must be unstable. Consequently, the bilayer remains stable against lipid protrusion and evaporation, as desired. So far, we have decided only on the signs of wd  for the canonical configurations listed in Figure 2 by considering their stability or instability. Thermodynamic reasoning enables us to go further. Let fluid T denote the physiological temperature at which the model lipid bilayer must remain stable and fluid. Because cut r is finite, this requires each of the intraleaflet and interleaflet nearest neighbor pairwise binding energies to be comparable in order of magnitude to the thermal energy, viz.

Lipids 1 and 4 in
, where k denotes the Boltzmann constant. In addition, to ensure stability against thermally induced leaflet interdigitation as well as lipid protrusion and evaporation, and to enforce the interbilayer hydration barrier, we require the three unbound canonical configurations to have sufficient S. Dey, J. Saha SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 4 of 16 potential energy barrier at the physiological temperature. Eq. (4) summarizes our conclusions for later reference. TTC  HHC  SSA  HTC  wd  wd  wd  wd  wd   SSP  TTC  wd  wd  fluid   HHC  SSA  HTC  wd  wd  wd  fluid   ,  0,  ,  ,  0 ,,

B. Model Definition
In view of Sec. II.A.3, to define SiMPLISTIC completely, we only need to specify the steric profile, viz. Steric () Vr and zfˆˆ( , , ) ij i j r r u u , the anisotropic well-depth wdˆˆ( , , ) ij i j  r u u , and the tail () sr .

Steric profile
Mesogenic molecules are frequently modeled as uniaxial, elongated, rigid objects in liquid crystal simulations, and two simplest shapes normally chosen to represent the shape of the molecular cores 43 are-a) the spherocylinder and b) the Gay-Berne spheroid f . The spherocylinder is more efficient for Monte-Carlo (MC) simulations, whereas the GB spheroid is better suited for Molecular Dynamics (MD). This is because, although the contact distance (or distance of closest approach) is simpler to compute for the spheroylinder 46 , the computation of the gradients of the contact distance as required for force and torque determination in MD is far cheaper for the GB spheroid 47 . Because our research group uses MD exclusively, we chose the GB spheroid for the core of SiMPLISTIC. g For a GB spheroid of width 0  and length e  , the contact distance is given as The anisotropy parameter  is related to the aspect ratio 0 :

The tail of the model potential
Deserno and coworkers clearly demonstrated the efficacy of a broad attractive tail potential for generating fluid bilayers through selfassembly 31 . In our 2017 amphiphile model 35 and our subsequent work on bolas 36 , we too showed successful self-assembly, using a simple f A spheroid is an ellipsoid of revolution and hence is uniaxial or axisymmetric. The GB spheroid is a computationally efficient approximation to an actual spheroid. g For a corresponding MC study, a spherocylindrical core should give the same results as what we obtained using the GB spheroid and MD. cubic spline as the broad attractive tail h . In the same vein, () sr for SiMPLISTIC is now defined as Note that the difference cut zf rr  determines the range of pairwise attraction for attractive configurations like SSP and TTC and the range of hydration barrier for HHC. This difference is therefore a parameter of SiMPLISTIC which we shall call range .

Well-depth
In view of the discussion in the last paragraph of Sec. II.A.3, we shall try to form wdˆˆ( , , )  r u r u . Note that these three terms are also sufficient to distinguish between all the canonical configurations in Figure 2. This is illustrated in Table 1 below.  Also, wd  for the SSP and SSA configurations are respectively 0  and 0   , which complies with the sign criterion given in the first row of Eq. (4). Note that, of the three parametersm , a and b , m has no effect on the values of wd  for the different canonical configurations because ˆˆ1 ij  uu for each of them. We, therefore, need to tune a and b in order to make wd  satisfy the remaining conditions in Eq. (4).
It may also be noted that, fixing a and b this way, leaves m as the only free parameter of the model well-depth, which is very much desirable as the number of parameters becomes minimal.
The two unknowns a and b can be solved for if wd  is known for any two of the collinear canonical configurations. To this aim, let us choose wd  for the configurations TTC and HTC simply as 0  and 0   respectively, in compliance with the sign conditions in Eq. (4).
Solving for a and b now, we get 2, 4 3 ab  .  Table 2. Well-depth wd  for the canonical configurations in Figure 2.

C. Model Summary
To summarize, SiMPLISTIC represents lipids in implicit solvent as directed spheroids, which interact using the following potential.  5)], SiMPLISTIC has only two parameters, viz. range and m . The range parameter determines the number of interacting pairs of molecules and hence, quite trivially, affects the overall binding energy of the lipid aggregates which, in turn, affects thermal stability and elasticity 30,31 . m , on the other hand, tunes the degree of alignment between the interacting lipids and therefore, for sufficiently large values, should guarantee the formation of bilayers where the lipids possess high orientational order. In the next section, we shall test this hypothesis and also examine SiMPLISTIC's selfassembly behavior for smaller values of m .

III. SELF-ASSEMBLY, TUNEABILITY, MIXING RULE
In this section we report self-assembly behavior of SiMPLISTIC (as seen in MD simulations for different values of m ), analyze it and extend SiMPLISTIC to model new cases based on that analysis.

Protocol
Self-assembly in an implicit solvent system is usually studied under the NVT ensemble for both MD 26,31,48  In order to study spontaneous, free self-assembly from a disordered, isotropic phase under NVT, choosing an appropriate number density is very important. Implicit solvent systems are far less dense than anhydrous systems such as systems of mesogens, simply because there is an implicit presence of water in the form of unoccupied volume within the simulation box. Too low a number density, however, is not conducive to self-assembly as can be verified easily. One way to get a viable number density for ISCG lipid systems is to start with an approximately valid preassembled structure such as a system of wellseparated bilayers. The disordered, isotropic configuration can then be generated simply by evaporating the preassembled structure at a high temperature under NVT. This way of producing a randomized initialization also eliminates any undesirable core-core overlap between the lipids 28 . As depicted in Figure 5, we generated our disordered, isotropic configuration by evaporating two lattice membranes that were assembled on square grids. To ensure enough space in the simulation box, the two bilayers were kept out of each other's range of interaction, viz. the interbilayer separation was made greater than range . This and selfassembly was studied for that temperature for different values of m . We infer that the lipids were free to arrange themselves during selfassembly because the overall internal pressure or stress, as computed using the virial 54 , remained near zero.

Self-assembly
Self-assembly was studied for 0,0.5,1,...,4 m  at 0 kT   . Figure 6 shows some of the phases obtained. Self-assembly was rapid with structure being discernible by 4000 steps and prominent aggregate formation by 6000 steps. For 03 m , closely packed inverted micelles were obtained, with the micellar curvature decreasing with increasing m . For 3 m  , a preference for bilayer formation was clearly visible, with almost no negatively curved phase remaining in the simulation box at 4 m  . The parameter, m , therefore, determines whether the self-assembled phase will be non-lamellar or lamellar. We also studied self-assembly for some random values of 4 m  and only lamellar bilayers were obtained. To summarize our observations, the maximal curvature for non-lamellar phase is obtained for 0 m  , the curvature of non-lamellar phase decreases with increasing m , and no intrinsic curvature, negative or positive, remains in the system for 4 m  .

B. How m tunes curvature
To explain the aforementioned observations, let us analyze how m affects the energy of lipid splay. We define splay as any configuration between a pair of lipids and ij where each lipid is a mirror image of the other with respect to a mirror plane that bisects and is normal to ij r .
The angle between the lipids is called the splay angle,  (Figure 7). Note that ,        lead to the SSP, TTC and HHC configurations respectively (cf. Figure 2). As seen from the plot in Figure 7, the most S. Dey, J. Saha SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 6 of 16 preferred splay angle max  , viz.  at which wd  is maximum, is negative for each 0 m  . Also, the amount of splay in this configuration, viz. max  , decreases with increasing m . Notice that, because the model lipids are cylindrically symmetric, a preference for non-zero splay implies a conical molecular profile. If the most preferred splay angle is negative, viz. max 0   , the effective molecular shape must be that of an inverted cone with the angle between the cone axis and cone surface given by max 2  . Inverted cones form nonlamellar phases with negative intrinsic curvatures such as the inverted micelles. As max 0   , the effective molecular shape becomes cylindrical which favors the formation of lamellar bilayers with zero spontaneous curvature. From this perspective, it is clear how max  governs the spontaneous curvature of the self-assembled phases. Because max  , in turn, is tuned by m as discussed above (Figure 7), this explains the observations reported in Sec. III.A.2. Note that this also shows how, despite having a fixed spheroidal core, SiMPLISTIC lipids can have their effective shape or packing parameter 55 varied with m and thus exhibit polymorphism on self-assembly.

C. Mixing Rule for Lipids with Different Values of m
A logical extension of the discussion in the last subsection (Sec. III.B) is to ask how to model the interaction between a pair of SiMPLISTIC lipids each of which generates a different spontaneous curvature in single-component self-assembly. In other words, we seek to extend SiMPLISTIC so that mixtures of lipids with different packing parameters, such as lamellar and non-lamellar lipids, can be simulated.
To respectively. To test the validity of this mixing rule, we studied self-assembly for a binary mixture of 800 lamellar and 800 non-lamellar lipids with 4 and 0 m  respectively. Figure 8 depicts a few close-ups of the aggregates formed. Note how the lamellar lipids (shown in red) intermix with the non-lamellar lipids (shown in green) in forming the inverted micelle (Figure 8c). The preference of the lamellar lipids to form bilayers is also discernible from the lamellar clusters of the red lipids just beside the inverted micelle (Figure 8b).

IV. MODEL MEMBRANES, HYDROPHOBIC MISMATCH
This section reports how SiMPLISTIC performs for model membrane simulations. Later in this section, we shall also discuss an extension of SiMPLISTIC to include the effect of hydrophobic mismatch in multicomponent bilayers. Because SiMPLISTIC is required to model only lamellar lipids for membrane simulations, we must take 4 m  in view of the self-assembly results reported in Sec. III.A.2. All simulations, however, are performed with 4 m  only (also see Sec. IV.A.4 for further remarks).

A. Bilayer Characteristics
Model membrane simulations are usually performed under zero lateral tension 30,32 . For such constant tension simulations, we adopted a two-dimensional Nosé-Hoover NPT dynamics 54 in the bilayer plane to regulate the lateral pressure. To prepare the initial state, a preassembled box-spanning bilayer was thermalized at 0 kT   under zero tension. As desired, the surface tension, measured using the diagonal elements of the internal pressure tensor 56 , remained near zero 57 for every production run.

Lateral diffusivity, area per lipid : Fluid and gel phases
The in-plane or lateral diffusivity ( || D ) of lipids in the membrane is computed from the limiting slope of the mean square displacement (MSD) in the bilayer plane, using the Einstein relation for twodimensional diffusion 58,59 . Figure 9 depicts the MSDs plotted for four NVE runs starting with equilibrated tensionless bilayers, with different mean temperatures. The corresponding diffusivities and projected area per lipid ( L A ) are listed in Table 3 A fluid-gel transition is typically associated with significant drop in || D , reduction of L A and increased (hexagonal) bond orientation order 28 . From the MSDs plotted in Figure 9 and a visual inspection of the lipid arrangement in the bilayers depicted there, it is apparent that the membrane is fluid at

Order parameter
The orientational (uniaxial) order parameter for lipids in a bilayer is given as , where n is a unit vector along the average normal of the bilayer. For a tensionless membrane at fluid kT , the order parameter of SiMPLISTIC lipids is 0.967. From this high amount of orientational order, and a visual inspection of the side view of the bilayer as provided in Figure 9, it is safe to infer that the fluid bilayer is not tilted. No tilt was found in the gel phase as well.

Bilayer thickness
The thickness of a bilayer is estimated as follows. Taking the X and Y axes in the bilayer plane, the bilayer thickness lies along the Z direction. The Z coordinate of the head end of a lipid is given by Because our fluid membrane is thicker than this, it implies that there is no interdigitation in the fluid phase, as is also apparent from the side view of the bilayer in Figure 9. No interdigitation was found in the gel phase as well.

Bending modulus
Being designed with a top-down, phenomenological approach, getting the membrane elasticity such as the bending rigidity right is a decisive factor for SiMPLISTIC. All other top-down attempts at ISCG modeling of lipids with rigid models [25][26][27][28] had failed to produce experimentally relevant values for the bending moduli, as was mentioned in Sec. I. There are several methods for determining the bending modulus C  from configurations generated by molecular simulations 60 . Of these, a method based on spectral analysis of the S. Dey, J. Saha SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 7 of 16 height fluctuations of a tensionless bilayer has been used by most previous works on ISCG modeling of lipids. This method is difficult and computationally expensive, which has to do with its requirement of very large system sizes (continuum limit) 31 and a robust sampling of global bilayer undulations to apply the Fourier transform on. In contrast, a relatively recent and computationally simpler method 61 based on real-space fluctuations (RSF) of the splay degrees of freedom is local in nature and, therefore, can be applied to rather small system sizes (as small as ~100 lipids), which makes it more practical for MD simulations. This method has also been successfully validated on a range of lipid systems 62 . For our case of flat bilayer composed of SiMPLISTIC lipids with director û , the RSF method is briefly described as follows.
According to the RSF methodology 61,63 the splay at any point p on the bilayer is given as the local covariant derivative of the vector field  un along any vector ˆ e tangent to the bilayer. Because the bilayer normal n does not vary across a flat bilayer and ˆˆ0 , the splay is given as (11) C kT  can therefore be determined from the variance of a Gaussian with zero mean fitted to the frequency distribution of S  (Figure 10). Employing this method, we computed the bending rigidity of a tensionless SiMPLISTIC bilayer at five different temperatures, as listed in Table 4

Area compressibility modulus and elastic ratio
For constant tension simulations where the membrane area fluctuates about a mean value to regulate the surface tension, the best suited measure for area compressibility modulus A  is given by 56 A denotes the total (projected) area of the membrane and denotes the time average obtained from a production run of the simulation. The elastic ratio is a dimensionless quantity defined as where h denotes the membrane thickness 28,29 . b ranges from 4 12  for bilayers with strong interleaflet coupling and 16 48  for unconnected bilayers 66

B. Hydrophobic Mismatch
Hydrophobic mismatch is a key player in the physics of multicomponent lipid membranes if the constituent lipids differ in length due to a difference in either chain length or chain saturation [70][71][72] . Such difference in length between neighboring lipids in a leaflet exposes some of the hydrophobic parts of the longer lipid(s) to the surrounding water, which is unfavorable. To avoid this 'hydrophobic mismatch', lipids prefer to be surrounded by lipids with similar length, which may lead to a demixing of dissimilar lipids in the multicomponent system 73 . In other words, hydrophobic mismatch generates a line tension that, if large enough, may generate phase separation through domain formation 72,74,75 .
SiMPLISTIC, as defined so far (Sec. II.C and III.C), can only model systems where all lipids have the same length ( e  in Sec. II.B.1). In order to model multicomponent membranes with hydrophobic mismatch, it therefore needs to handle interactions between lipids with different lengths. We extend SiMPLISTIC in the following by defining just such an interaction and test our model with simulations.
Let us take two types of SiMPLISTIC lipids (A and B) that differ only in their spheroidal lengths, A e  and B e  . The contact distance between two such GB spheroids of unequal lengths had been worked out by Cleaver et al. 76 . Accordingly, the contact distance  ˆ, between a lipid i of type A and another lipid j of type B is given as S. Dey, J. Saha SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 8 of 16 r u r u r u r u u u uu (14) where A  and B  denote the anisotropy parameters of i and j respectively, as defined in Sec. II.B.1. Note that for i and j lipids having identical lengths, AB     , which reduces  in Eq. (14) to GB  in Eq. (5). With the contact distance (  ) obtained as above, we can now choose Eq. (6) as the steric part of the potential between i and j , with R now defined as 00 ( ) / Rr       . Hydrophobic mismatch driven demixing, interpreted as a tendency of lipids of similar length to cluster together within a leaflet, can be modeled through the mechanism of like lipids interacting more strongly than unlike ones. This is most easily achieved by simply reducing the strength of A-B interaction compared to that of A-A and B-B. Care must however be taken to do this for lipids within the same leaflet only. Otherwise, the weaker attraction between unlike lipids from opposing leaflets would lead to direct interleaflet coupling through domain registration. This is unphysical because hydrophobic mismatch, in reality, destabilizes registration rather than favoring it. For further illustration of this last point the reader is referred to Ref. 74 and the Fig. 2 and its corresponding discussion in Ref. 77. In We, therefore, need to define ( , ) AB ee p  now. To this end, note that this p can be interpreted as a dimensionless penalty factor which determines the loss of binding energy due to lipid length mismatch. This loss of binding energy is what gives rise to a line tension  between two domains of unlike lipids within a leaflet, and is directly proportional to it 72 . The line tension between domains of A and B, therefore, is directly proportional to the penalty factor, i.e. p   .
Kuzmin et al. 78 , in Eq. 17 of their paper, have derived a quadratic dependence of the line tension on where  is the amount of hydrophobic mismatch and 0 h is the average thickness of the stretches of leaflet within the two domains. This dependence is exact if the two domains have the same spontaneous curvature, which is a condition that applies well for flat membranes composed of lamellar lipids. In our case,  Figure 11). As desired, no domain registration was observed between the opposing leaflets of the bilayer (Figure 11). SiMPLISTIC therefore, with the inclusion of interaction between lipids of different lengths as described above, succeeds in modeling hydrophobic mismatch.

A. Highlights
This paper described a novel, coarse-grained model called 'SiMPLISTIC' for large scale, particle based simulations of lamellar and non-lamellar lipid systems in implicit water. In absence of explicit water, the hydrophobic effect and the hydration force are mimicked with lipid-lipid attraction and repulsion respectively.

B. Benchmarks
SiMPLISTIC is designed for speed. Here we report some crude benchmarks for SiMPLISTIC and compare it with some other relevant models.
The execution time of an optimized code is not a very meaningful metric for the underlying model's computational efficiency. This is because, apart from the platform, CPU time depends on the S. Dey, J. Saha SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 9 of 16 optimizations performed within the source code as well as during compilation, both of which are subject to variation from one implementation to another. Benchmarks under non-optimal conditions, on the other hand, are more useful because they represent the least performance expected from a model running on any given processor. With this in mind, we benchmarked SiMPLISTIC as follows.
The test code, implemented in Fortran 90 with double precision for the floating-points, consisted of just the SiMPLISTIC force routines and the Nosé-Hoover NVT integrator. To eliminate I/O calls during run, no configuration dumping was done. No parallelization or SIMD vectorization and no domain decomposition or neighbor list were appliedhence, the force routines enumerated every possible pair of lipids serially in the least efficient way. The compilation was done using gfortran with no optimization. The simulated system was the same as reported in Sec. III.A.1, with 1 m  and starting from a selfassembled phase of closely packed inverted micelles at 0 kT   . Although the code ran serially on a single core, no core affinity or task priority was set. Running this test code for several times on a laptop with an Intel i3 processor (clock speed 2.13 GHz) and 32-bit GNU/Linux, we consistently obtained CPU times of 2.05 minutes per 1000 MD steps. We also timed the self-assembly of the same system at 6.06 minutes (6000 MD steps). It goes without saying that with multicore parallelization, SIMD vectorization and GPU accelerations that are routinely performed nowadays, the execution time will be improved manifold. It may also be pointed out that the computation of To compare the computational performance of SiMPLISTIC to that of a standard single-site, rigid model with a similar anisotropic pair potential (viz. Gay-Berne 40 ), and also a multi-site, bead-spring model for lipids in implicit solvent (viz. the Cooke-Kremer-Deserno model 30 ), consider the following. SiMPLISTIC contains less floating point operations (flops) than a switched Gay-Berne (GB) potential and therefore, should be faster than the GB benchmarks on any given platform. SiMPLISTIC is also faster than a single-site amphiphile model previously published by us 35 , again by virtue of requiring fewer flops. Being single-site, this amphiphile model was found to be more than 1.6 times faster than a rigid version of the multi-site Cooke-Kremer-Deserno model, which, in turn, is naturally faster than the original flexible (bead-spring) version 30 . SiMPLISTIC, therefore, should be more efficient than the Cooke-Kremer-Deserno model for large scale phenomena that are not significantly sensitive to lipid chain flexibility.

C. Scope
In view of the features listed in Sec. V.A, SiMPLISTIC should be able to serve as an excellent model for large-scale lipid simulations in academic and industrial research. Apart from research, the model should also be effective in education. For example, SiMPLISTIC's conceptual simplicity, versatility and computational speed make it an ideal model for learning about complex, cooperative phenomena in lipid-based systems through hands-on, interactive simulations.
To inspire future applications, some specific use cases for SiMPLISTIC are suggested below.
-Studying elasticity and nanoscale dynamics in multicomponent membranes as functions of the membrane constitution.
-Studying vesicles consisting of lamellar lipids on the outer surface and non-lamellar lipids on the inner. Hydrophobic mismatch driven demixing may induce budding in multicomponent vesicles. -Interactive molecular dynamics.

Two more phases with negative curvature
Although the free parameter m tunes the spontaneous curvature continuously from inverted micelles to planar bilayers, it is curious that no inverted hexagonal (H II ) or bicontinuous/bilayer cubic (Q II B ) phase was obtained through self-assembly in our simulations. We suspect that a cubic simulation box with periodic boundary conditions, as employed by us, cannot form the H II phase due to the lack of hexagonal symmetry. On the other hand, because it involves curved bilayers, the Q II B phase would probably require a much larger system size. We hope to address these issues in future work(s).

Possible extensions of SiMPLISTIC
SiMPLISTIC, as defined in this paper, is useful for simulating systems consisting of uncharged lipids only. Extending SiMPLISTIC to include interactions with purely hydrophilic and hydrophobic particles, which may be in the form of GB spheroids or spheres, will significantly increase its scope. Such an extension, for example, might be useful in modeling lipid-protein interactions and studying membrane permeability or leakage of cargo from cubosomes and liposomes.
In addition to this, defining electrostatic charges on SiMPLISTIC lipids and other particles will make the model even more effective. Such an extension might enable probing important large scale phenomena such as lipid-DNA complex formation in the presence of counter-ions, fusion of oppositely charged vesicles etc.

An open question
SiMPLISTIC may be seen as a proof of concept, which demonstrates that a single-site, rigid model designed using some simple guidelines (Sec. II.A) can spontaneously self-assemble into fluid phases and generate elastic properties consistent with experiments and other simulations. Moreover, it shows how the intrinsic curvature of the self-assembled phase (or equivalently, the molecular packing parameter 55 ) may be tuned with only a single parameter. Our previous single-site and rigid amphiphile model 35 was also capable of tuning the curvature with a single parameter. However, that model gave phases with positive curvatures, whereas SiMPLISTIC gives phases with negative curvatures. To resolve this dichotomy, it is natural to imagine a single-site, rigid model that can self assemble into phases with all possible curvatures-both positive and negative. Apart from its sheer elegance, such a minimalistic, molecular level model will also be extremely useful, as amphiphiles with any packing parameter may be modeled with it, including mixtures of molecules with different   . How the disordered, isotropic phase was prepared from preassembled lattice bilayers with sufficient interbilayer distance. SiMPLISTIC lipids are depicted as red spheroids with their head-ends marked with blue beads. The sky blue background implies the aqueous environment. This and other such renderings depicted in this work were performed using the molecular graphics software QMGA 79 . Figure 6. Self-assembled phases formed by SiMPLISTIC for different values of m . Lipids are depicted as red spheroids with their head-ends marked blue. The sky blue background implies the water environment. The non-lamellar phases are formed by closely packed inverted micelles (compare with the schematic cross-section in Figure 1). For enhanced clarity, a slice of an inverted micelle is also shown. The micellar radii are greater (and equivalently, the curvature is lesser) for 2 m  compared to 0 m  . For 3 m  bilayers are obtained-a close-up of a box-spanning bilayer is shown [the stray lipids below the lower leaflet of the bilayer are actually part of a separate aggregate in the background (that we have sliced off for clarity) and do not belong to the bilayer in focus].

S. Dey, J. Saha
SiMPLISTIC: A Novel Pairwise Potential for Implicit Solvent Lipid Simulations … arXiv:2010.00561v2 Page 15 of 16 Figure 7. Well-depth wd  of the SiMPLISTIC potential as a function of the lipid splay angle  for different values of m . The peak shifts towards   for increasing m . Hence, increased m leads to a preference for decreased splay. It is easy to see from the spheroid configurations that packing pairs of lipids with   would generate inverted micelles (Figure 1c). Hence, the effective molecular shape (packing parameter) of lipids preferring   is similar to an inverted (truncated) cone.   . For 3   , the number ratio of red to green is 3:1. For  , the number ratio is 1:1.
Although clustering of lipids of the same color in each leaflet is clearly visible for both values of  , a much cleaner phase separation occurs at   due to an increased line tension. Note that there is no domain registration between the leaflets.