To date, there are four main approaches for the construction of compositionally asymmetric bilayers for MD simulations:

• Ensure equal numbers of lipids in the two leaflets (EqN).
• Match the surface areas (or lipid packing densities) of the two leaflets to those from cognate symmetric bilayers (SA).
• Eliminate differential stress, i.e., ensure zero leaflet tension (0-DS).
• Emulate biological asymmetry (EmBioAsym).

Most of these methods require performing at least one simulation to arrive at the desired asymmetric membrane model. Once the relative leaflet compositions and abundances are determined, the bilayer can be simulated with any conventional software to investigate its properties in detail. Below we briefly describe each of these approaches and discuss a few studies in which they have been applied.

Since in symmetric bilayer simulations, the two leaflets generally have equal numbers of lipids (Figure 2A), one approach to building asymmetric bilayers is to ensure that same condition (Figure 2B). This can be accomplished by building a symmetric bilayer and replacing individual lipids with different ones, or by specifying the same total numbers of lipids in the two leaflets and letting software generate the asymmetric lipid compositions. Here, all lipids are treated identically, thus disregarding any structural differences or packing preferences between species. This approach has been used to study the physical properties of both simpler lipid mixtures [27,30] and asymmetric plasma membrane models of increasing complexity [75,76].

An alternative to the EqN method is choosing the relative numbers of lipids in the two leaflets so that at the initial stage of bilayer construction, their average areas per lipid (or surface areas) match those of cognate symmetric bilayers (Figure 2C). This can be accomplished by first simulating two symmetric bilayers with the lipid compositions of the asymmetric membrane leaflets, followed by either stitching their leaflets together or using the obtained equilibrium areas per lipid to calculate the respective numbers of lipids and building the asymmetric bilayer from scratch. The latter approach is more general and allows for the construction of asymmetric bilayers of different sizes. This method ensures that the relative packing densities of the two leaflets remain fixed even as the two individual leaflet areas may simultaneously increase or decrease in comparison to their symmetric counterparts. This protocol has been applied to the analysis of the membrane potential [24,25,26], interleaflet coupling [29,32,34], the permeability of plasma membrane models [77], and cholesterol interleaflet distribution [39,42].

To avoid simulating two symmetric bilayers, one can alternatively use reported areas per lipid (APL) for the individual lipids and, if the leaflets contain more than one component, assume ideal mixing and calculate the respective mole-fraction-weighted packing densities. For instance, CHARMM-GUI [54] uses individual lipid APLs to estimate the number of lipids in each leaflet when only lipid composition and lateral box dimensions are provided. This approach was recently used in a study comparing different construction methods for asymmetric bilayers [8].

Recent experiments have found that preferred lipid packing densities in symmetric bilayers can be altered in asymmetric bilayers, presumably due to the effects of the interleaflet coupling [78]. A corollary is that lipid areas in asymmetric bilayers cannot be estimated a priori from their values in symmetric bilayers. It follows that the EqN and SA methods described above may produce bilayers with differentially stressed leaflets [67,73]. To ensure zero leaflet tension (Figure 2C), one can follow the protocol outlined in [67], which is based on the following principle. An asymmetric bilayer is first constructed (for example, using the EqN or SA methods) and simulated until the APL is converged over the last ~200 ns. Then, the lateral pressure profile is calculated from the converged portion and used to obtain the leaflet tension. If the tension is non-zero, the number of lipids is adjusted by either removing lipids from the leaflet with negative tension or adding lipids to the leaflet with positive tension. The exact number of lipids to add or remove can be chosen on a trial-and-error basis or estimated from the relationship between the bilayer area compressibility modulus and leaflet tension (Equation 3 in [67]). A new asymmetric bilayer is then built from scratch with the updated leaflet lipid numbers and the same steps are repeated (i.e., simulation, calculation of leaflet tension, adjustment, and rebuilding if necessary) until zero leaflet tension is reached. This approach has been used to examine the effects of asymmetry on the mechanical properties of anionic asymmetric bilayers [74] and the ability of gramicidin to scramble lipids [46], as well as the effects of differential stress on the interaction of small molecules with asymmetric membranes [73].

The asymmetry in cell plasma membranes (PM) is maintained via the activity of flippase and floppase enzymes which, when active, move lipids between leaflets against their concentration gradients [79]. Since the lipid specificity of these enzymes is arguably restricted to certain lipid types [80], one hypothesis is that cells regulate their PM lipid organization by restricting the asymmetry of some lipids while letting others equilibrate between leaflets according to their chemical potential. One notable example is cholesterol, a major component of mammalian cell plasma membranes that can rapidly flip between leaflets [6,81,82,83]. Interestingly, while cholesterol is capable of alleviating stresses in the membrane [84], its strong preference for interaction with saturated lipids may dominate over elastic and entropic forces and drive its distribution in a way that increases the differential stress [10,36,85]. This illustrates both the natural tendency of a bilayer constituent to equilibrate its distribution based on its chemical potential and the fact that realizing this tendency may produce rather than eliminate stresses in the membrane. In that respect, two methods have emerged to examine the equilibration of lipid redistribution in a simulated bilayer in the presence of imposed asymmetry [8,9].

The first approach involves simulations in the NPT ensemble and utilizes P2₁ boundary conditions (Figure 2E) [8]. These PBCs allow lipids to sample both leaflet environments by freely exchanging between leaflets during the simulation [63]. To mimic the activity of flippases and simulate asymmetric membranes, the method involves constraining some lipids to stay in one leaflet while allowing others to equilibrate their leaflet concentrations via an interleaflet redistribution [8]. Thus, it is possible to start with a bilayer constructed with one of the methods described above, then restrict some lipids to their respective leaflets and simulate the system with P2₁ PBCs to examine the preferred lipid distribution of the unconstrained membrane components in the presence of the imposed asymmetry. Consequently, since the relative numbers of lipids in the two leaflets are not constrained, they can dynamically change during the simulation. As noted by the authors, while opening transient pores in the membrane can also accelerate the exchange of lipids between leaflets, the advantage of P2₁ PBCs is that the chemical equilibrium reached by the freely diffusing lipids is a property of the asymmetric leaflets in the absence of any mechanical perturbations such as those imposed by a pore [8].

The second approach starts with a compositionally symmetric bilayer and replaces some of the lipids with new ones to generate the initial asymmetry (if the bilayer contains the same lipid numbers across leaflets, this is equivalent to the EqN method). It then proceeds with a simulation not in the NPT ensemble (as discussed above), but instead in a semi-grand canonical ensemble that allows dynamic changes in lipid identity (specifically their saturation or headgroup type) during the simulation (Figure 2F) [9]. This approach emulates the action of lipid-translocating enzymes by imposing a chemical potential difference between some molecules in one leaflet while letting the leaflet lipid compositions adjust in accordance with the chemical potential of all lipid species subject to the imposed constraints. In these simulations, the lipid number asymmetry is fixed, but the changes in the degree of lipid saturation or type across species dynamically alter the relative packing densities of the leaflets by virtue of their changing lipid compositions. This method is well suited for investigating how some asymmetries might naturally arise from others and has thus helped explain certain experimental observations of the leaflet lipid compositions in erythrocyte membranes [9].