The main transition in the Pink membrane model: finite-size scaling and the influence of surface roughness

We consider the main transition in single-component membranes using computer simulations of the Pink model [D. Pink {\it et al.}, Biochemistry {\bf 19}, 349 (1980)]. We first show that the accepted parameters of the Pink model yield a main transition temperature that is systematically below experimental values. This resolves an issue that was first pointed out by Corvera and co-workers [Phys. Rev. E {\bf 47}, 696 (1993)]. In order to yield the correct transition temperature, the strength of the van der Waals coupling in the Pink model must be increased; by using finite-size scaling, a set of optimal values is proposed. We also provide finite-size scaling evidence that the Pink model belongs to the universality class of the two-dimensional Ising model. This finding holds irrespective of the number of conformational states. Finally, we address the main transition in the presence of quenched disorder, which may arise in situations where the membrane is deposited on a rough support. In this case, we observe a stable multi-domain structure of gel and fluid domains, and the absence of a sharp transition in the thermodynamic limit.


I. INTRODUCTION
Lipid membrane bilayers are abundant in nature and to understand their properties is of paramount importance [1][2][3]. One aspect that has received much attention are collective phenomena (phase transitions) taking place in these systems. Among the different phase transitions that can occur [4][5][6][7], the main phase transition is presumably the most important and well studied one [8,9]. This transition, typically driven by the temperature T , is between a "gel" and a "fluid" phase. At low T , the bilayer is in the gel phase (characterized by nematic chain order of the lipid tails), while at high T the bilayer assumes the fluid phase (characterized by the absence of nematic chain order).
Computer simulations have become a well established tool to model the main transition. The challenge in simulations is to strike a balance between the level of detail to include, and the time and length scale one wishes to address [10]. Since collective phenomena involve many molecules and entail large length scales it is clear that, in order to describe the main transition, a significantly coarse-grained particle model is crucial. Strictly speaking, one needs to address the thermodynamic limit (infinite particle number) since only there phase transition properties become properly defined. Indeed, the need for coarse grained modeling of lipid bilayers is well recognized [11][12][13].
An early and highly successful coarse grained approach to study the main transition has been the particle model introduced by David Pink and co-workers [14][15][16]. In this model, the so-called Pink model, only the orientational degrees of freedom of the hydrophobic lipid tails are included, while the positional degrees of freedom of the hydrophilic heads are disregarded. This model, due to its simplicity, allows for the investigation of very large systems, and the nature of the main transition can be probed in great detail. Indeed, key features of the main transition in the Pink model compare well to experiments [17].
However, despite the great success the Pink model has enjoyed, there remain some open questions. In Ref. 18, it was noted that the Pink model at the experimentally determined transition temperature does not undergo any transition. While in systems of finite size there were indications of a transition, these vanished in larger systems. This raises the question as to why no transition could be detected. The aim of this paper is to resolve this issue. As it turns out, to properly model the main transition, a finite-size scaling study is essential. Computer simulations inevitably deal with only a finite number of particles, and their output will depend on the number of particles used, especially near phase transitions. Finite-size scaling provides the framework to systematically extrapolate simulation data to the thermodynamic limit. To date, finite-size scaling studies of the Pink model are scarce, with Ref. 18 being a notable exception. The present paper aims to fill this gap. Our main finding is that, in order to observe the main transition in the Pink model at experimentally relevant temperatures, one of the model parameters needs to be adjusted. This follows quite naturally when one realizes that the universality class of the Pink model is just the one of the two-dimensional (2D) Ising model [15]. As we will show for three lipid species, the "standard" Pink model parameters yield a critical temperature distinctly below the experimental main transition temperature. Consequently, a "re-tuning" of the standard Pink parameters is urgently needed.
As an application, we also address the fate of the main transition in the presence of quenched (immobilized) impurities using the Pink model. The experimental motivation to do so is that this situation may resemble that of a membrane supported on a rough substrate. In binary lipid mixtures, the effect of such impurities on lateral phase separation has recently attracted much attention [19][20][21][22][23][24]. In this paper, we present simulation results for the corresponding scenario in a single component bilayer undergoing the main transition. Within the framework of the Pink model, we find that quenched impurities prevent the main transition from taking place, already at low impurity concentrations. Instead of the formation of macroscopic gel and fluid domains, we now obtain a stable multi-domain structure, which strikingly resembles experimental results. The theoretical justification is that the impurities induce a change in universality toward the 2D random-field Ising class. As is well known, the latter does not support an order-disorder phase transition in the thermodynamic limit [25][26][27][28].

II. THE PINK MODEL
In the Pink model, the lipid bilayer is assumed to consist of two independent monolayers. Each monolayer is represented by a triangular 2D lattice consisting of N sites, and each lattice site contains a single lipid chain. Each lipid molecule is comprised of two independent hydrophobic acyl chains and a hydrophilic polar head. The polar heads are translationally frozen to the lattice, and no particular structure for the polar head groups is assumed. The only degrees of freedom included in the Pink model are the acyl chain conformations. These are not simulated directly (i.e. one does not explicitly model the carbon atoms) but are captured in a coarse-grained fashion whereby the chain conformations are grouped into α = 1, . . . , q discrete states. The original Pink model uses q = 10, but we will consider different values also. These states include the ground state (α = 1), eight lowenergy excitations (α = 2, . . . , q − 1), while all remaining conformations are grouped into a single disordered state (α = q). Each state α is characterized by three coarsegraining parameters, namely an internal energy E α , a cross-sectional area A α , and a degeneracy D α counting the number of chain conformations with energy E α and area A α .

A. coarse graining parameters
To determine the coarse graining parameters, we assume that a single acyl chain consists of i = 1, . . . , M carbon atoms, thereby containing M − 1 carbon-carbon bonds, and that bonds are either in a trans or gauche configuration. The trans configuration yields the lowest energy, while the gauche configuration has a slightly higher energy. The energy difference between the trans and gauche configuration is denoted Γ (Table I). To understand the difference in geometry between trans and gauche bonds consider a chain segment of four consecutive carbon atoms. The positions of the first three atoms define a two-dimensional plane. In the trans configuration, the fourth atom remains in the plane, while in the gauche configuration, it leaves the plane, and it can do so inward or outward. Thus, each gauche bond is two-fold degenerate. In the Pink model, it is assumed that each 2n-th gauche bond takes the chain back to the original plane, and so the gauche degeneracy is given by where n denotes the total number of gauche bonds in the chain, and where the function ceil means "rounding-up" to the nearest integer. It is convenient to mathematically represent the chain conformations on a hexagonal lattice with next-nearest neighbor distance 2a. We emphasize that this lattice is merely an aid to identify the low energy chain conformations which are needed to set the coarse-graining parameters: it should not be confused with the triangular simulation lattice on which the Pink Hamiltonian will eventually be defined. The carbon atoms are placed on the nodes of the hexagonal lattice following certain rules, and nearest-neighbor connections between atoms represent carbon-carbon bonds. The ground state α = 1 corresponds to the chain conformation that is maximally stretched [ Fig. 1(a)]. Note that, in the ground state, the atoms are alternatingly placed on the left and right lattice node, yielding a characteristic "zig-zag" pattern. The ground state by definition contains only trans bonds, its internal energy is set to zero as a reference E 1 = 0, and it is obviously non-degenerate D 1 = 1. The cross-sectional area of the ground state has experimentally been determined as A 1 = 20.4Å 2 [14]. We also introduce the projected length l of the conformation, defined as the difference in z-coordinate between the carbon atom closest to the head group (i = 1) and the one furthest away (i = M ), with the z-direction as indicated in the figure. For the ground state, it follows that l 1 = (M − 1)a.
The eight low-energy excitations (α = 2, . . . , 9) are obtained by systematically incorporating gauche bonds. The effect of such a bond is to disrupt the "zig-zag" pattern of the ground state. That is, one no longer places the atoms alternatingly on left and right nodes, but also allows for "excursions" whereby for two consecutive atoms the same direction is chosen. Each such excursion corresponds to a gauche bond, and has energy cost Γ. The gauche bonds are introduced according to the following rules: (1) The two bonds in the chain closest to the head group must always be in the trans configuration. In Fig. 1, these correspond to the bonds between atoms 1 − 2 and 2 − 3. (2) At most three gauche bonds are allowed, and each time such a bond is included there is an energy cost Γ. (3) The projected chain length l must obey l 1 − l ≤ 3a. (4) The acyl chain cannot fold back onto itself. In the coordinate system of Fig. 1, this means that the z-coordinates of the atoms must obey z i+1 ≥ z i .
Following these rules, we show in Fig. 1(b) the chain conformations (i) and (ii) that form the first excited state α = 2. In (i), a single gauche bond is placed at the very chain end, while in (ii) it is placed at the secondlast position. One immediately sees that both conformations have the same energy E 2 = Γ, and the same projected length l 2 = (M − 2)a. To compute the crosssectional area, one assumes volume conservation for the lipid chains: A α l α = A 1 l 1 . Hence, from the (known) ground state values, the cross-sectional area of the excited state follows. Note that, by placing the gauche bond at the third-last position [ Fig. 1(c)], a shorter projected length is obtained, and so conformation (c) does not belong to the first excited state (even though it has the same energy). The total degeneracy of the first exited state D 2 = 4, which is the total number of conformations, multiplied by the gauche degeneracy of Eq. (1). The coarse-graining parameters of the remaining excited states can be found analogously, and are listed for completeness in Table I. Finally, for the completely disordered state α = q = 10, one assumes E 10 = (0.42M − 3.94) × 10 −13 erg, A 10 = 34Å 2 , and degeneracy D 10 = 6 × 3 M−6 , which have their origins in experimental considerations [16].

B. Pink model Hamiltonian
Having specified the coarse-graining parameters, the Hamiltonian of the Pink model can be written as [30] The first term is the total internal energy of the acyl with the sum over all N sites of the triangular lattice, and s(i) ∈ {1, . . . , q} the conformational state at the i-th lattice site. The second term represents the anisotropic van der Waals interaction between adjacent acyl chains H VDW = −J 0 i,j I s(i) I s(j) , with J 0 the van der Waals coupling constant, and i, j a sum over all 3N nearest-neighboring sites on the triangular lattice. The precise value of J 0 depends on the chain length, and explicit expressions are provided elsewhere [14,15,31]. However, it has been noted that these parameters do not always yield a main transition at the expected temperature [18], and so we will also propose our own values later on. The (dimensionless) variables I α measure nematic chain order, and can be expressed in terms of the cross-sectional areas [16,18] where ω 10 = 0.4 for the disordered state α = 10, and ω α = 1 otherwise. The last term in the Hamiltonian accounts for the interaction between the hydrophilic polar head groups and between them and water and also steric interactions from both head groups and the lipid chains. Although it is possible to consider a more realistic pairwise interaction between the headgroups [30], this interaction can be approximated with a simple pressure term H P = ΠA, where Π is an effective lateral pressure acting on the lipid chains in the bilayer membrane, and A the total cross-sectional area occupied by the lipids chains

III. MONTE CARLO METHODS
To study the phase behavior of the Pink model, we use the Monte Carlo (MC) simulation method. We mostly use triangular lattices of size N = L × L with periodic boundary conditions. The principal MC move consists of randomly picking one of the lattice sites, read-out the conformational state α of that site, and propose a new state β drawn randomly from the set of q possible states. The new configuration is accepted with the Metropolis criterion where D denotes the state degeneracy, ∆H the energy difference between initial and final configuration as given by Eq. (2), k B the Boltzmann constant, and T the temperature. The degeneracy compensates for the fact that some of the states have a much larger entropy, and should therefore appear more often in the ensemble average. By virtue of the MC move, the total projected area A given by Eq. (4) fluctuates during the course of the simulation. In fact, A plays the role of order parameter since it changes abruptly at the main phase transition. Hence, it is instructive to measure the distribution P (A|T, Π), defined as the probability to observe a configuration with projected area A. The distribution depends on the imposed temperature and pressure, as well as on the linear extension L of the triangular simulation lattice. At the main transition, P (A|T, Π) assumes a characteristic bimodal shape, from which a number of important phase properties are obtained (explicit examples are provided in the next section). We note that even with very long simulation runs, distributions P (A|T, Π) of high statistical quality are difficult to obtain, especially in the vicinity of the main transition. The reason is related to free energy barriers that arise from the formation of interfaces [32][33][34]. To overcome this problem, we combine our MC simulations with a biased sampling scheme called successive umbrella sampling [35]; the latter ensures that P (A|T, Π) is sampled accurately over the entire (specified) range in A of interest. A final ingredient to economize simulation time is the use of histogram reweighting [36]. A single simulation run yields P (A|T, Π) at a given temperature T and effective pressure Π; histogram reweighting enables us to extrapolate the measured distribution to different values T ′ , Π ′ . For example, extrapolations in the pressure are performed using P (A|T, Π ′ ) ∝ P (A|T, Π) e −(Π ′ −Π)A/kBT . Extrapolations in the temperature can be performed analogously, but also require storage of the energy histograms; for implementation details see Ref. 37.

A. the "standard" Pink model revisited
We first consider the main transition in a membrane consisting of DPPC lipids to settle a controversy when this system is being simulated using the Pink model. The acyl chains in DPPC consist of M = 16 carbon atoms, and the experimentally obtained main transition temperature T DPPC = 314.0 K [31]. However, simulations based on the Pink model could not detect a transition at this temperature [18]. The latter simulations used the "standard" Pink parameters as listed in Table I, van der Waals coupling constant J 0 = 0.710 × 10 −13 erg, and pressure Π = 30 dyn/cm. Hence the question arises as to why no transition could be detected. To answer this question we perform additional DPPC simulations using the Pink model, with the same parameters as in Ref. 18, but over a wider range in temperature and pressure. The picture that emerges is the following: At high temperature the distribution P (A|T, Π) is always singlepeaked (corresponding to one phase) for all value of the lateral pressure Π. At low temperature, P (A|T, Π) is doubled-peaked for a special value of the lateral pressure, Π = Π COEX , corresponding to two-phase coexistence [ Fig. 2(a)]. Here, the left peak reflects the gelphase, the right peak the f luid-phase. The numerical criterion to locate Π COEX is to vary Π until the fluctuation A 2 − A 2 reaches a maximum [38], with the ther-

mal averages computed as
At the temperature T = T c where the transition from a single to doubled-peaked distribution occurs, the system becomes critical. To locate the critical temperature a finite-size scaling analysis is performed, whereby we plot the Binder cumulant U 1 = ∆ 2 / |∆| 2 , ∆ ≡ A − A , versus temperature T for different system sizes L. In the thermodynamic limit while in systems of finite size, curves for different L intersect at T = T c [39,40]. In Fig. 2(b), we show the result for DPPC obtained using the "standard" Pink model parameters: the data scale as expected, and from the intersection the critical temperature T c can be accurately "read-off".
The corresponding estimates of T c as well as the coexistence pressures Π COEX for three lipid species are collected in Table II. For all lipid species considered, the computed critical temperature T c is distinctly below the experimental melting temperature T m . In other words: if one simulates the Pink model at the experimental melting temperature T m , one is always inside the one-phase region, where P (A|T, Π) is single-peaked! This, apparently, is the reason why no phase transition could be seen in previous studies [18]. One possibility to get the proper value for the transition temperature, i.e. such that T c coincides with T m , is to re-tune the value of J 0 . This has been done for the three lipid species by systematically changing the coupling constant J 0 using histogram reweighting and finite-size scaling. Our proposed values J * 0 and corresponding pressures Π * COEX for the three lipid species are summarized in Table II. For completeness, we still confirm the universality class of the critical point, which for the Pink model is expected to be the one of the 2D Ising model [15]. To this end, we consider the susceptibility χ = ∆ 2 − |∆| 2 /(k B T L 2 ) [41], which diverges at the critical point χ ∝ |t| −γ , t = T /T c − 1, with critical exponent γ. In systems of finite size, the divergence is rounded, but γ can still be obtained using the standard finite-size scaling procedure of plotting χ L −γ/ν versus t L 1/ν [42], where ν is the correlation length critical exponent. Provided suitable values γ, ν, T c are used, data for different L collapse. The result for DPPC is shown in Fig. 3, where the "standard" parameters of the Pink model were used. Indeed, by using the 2D Ising values {γ = 7/4, ν = 1}, and T c = 291.7 K of Table II, an excellent data collapse is observed (similar good collapses are obtained for DMPC and DSPC also). The order parameter critical exponent has also been measured, and the 2D Ising value β = 1/8 was confirmed (scaling plot not shown). Therefore, even though the Pink model is a 10-state model, its critical behavior remains in the universality class of the 2D Ising model. This further motivates the idea of reducing the q = 10 states in the Pink model to an effectively two-state description as is frequently done [15,17,[43][44][45].

B. modified Pink model with fewer states
We now consider the effect of lowering the number of states in the Pink model. For this purpose, an appropriate number of intermediate states was removed, based on the maximum number of gauche bonds. In the "standard" 10-state Pink model at most three gauche bonds are allowed. We now consider the case where at most two gauche bonds are permitted, by removing states α = 8, 9 from Table I, yielding an 8-state model (to keep the total number of states constant the degeneracy of the removed states was added to the disordered state, but we emphasize that this correction is small). We apply our previous finite-size scaling analysis to the resulting 8-state model for DPPC, using the "standard" value J 0 = 0.710 × 10 −13 erg. As expected, the critical point remains in the universality class of the 2D Ising model, but it is "shifted" to T c = 309.4 K and Π COEX = 26.0 dyn/cm. Similarly, by allowing at most one gauche bond, a 5-state model is obtained. In this case, the DPPC critical point is located at T c = 351.5 K and Π COEX = 87.7 dyn/cm. Therefore, lowering the number of states in the Pink model while leaving the other parameters untouched, one finds that both the critical temperature and pressure increase. This trend is consistent with the Pink model simulations of Ref. 46 for DPPC performed at the experimental melting temperature T m but with a lower number of states. In these simulations, hysteresis loops indicating a first-order transition are clearly visible around T m for q < 6. Indeed, as our scaling analysis shows, by lowering the number of states q, T c eventually exceeds T m , resulting in a genuine phase transition at T m .
To conclude: lowering the number of states q does not affect the universality class of the Pink model, which remains 2D Ising (provided q ≥ 2, of course). Hence, the topology of the phase diagram remains the same, merely the critical point gets shifted. Depending on the parameters used, the main transition in the Pink model is either first-order (T < T c ), or it is 2D Ising critical (T = T c ). We do not claim that the main transition as observed in experiments necessarily conforms to this scenario (we return to this point in Section V).

C. Pink model with quenched disorder
As a final illustration, we consider the main transition in a solid-supported membrane, which has received considerable attention in experiments [6,[47][48][49][50][51]. A striking feature observed in one of these studies is the formation of coexisting gel and fluid domains that do not coalesce with time, but instead form a multi-domain structure that is stable over hours [50]. To understand the stability of this structure is not trivial, due to the large amount of line interface it contains. Here we attempt to reproduce such a multi-domain structure within the framework of the Pink model. Our hypothesis is that the solid support onto which the membrane is deposited has a certain roughness. Since surface roughness is random and time independent, it constitutes a form of quenched disorder. We assume that this gives rises to regions on the surface where certain lipid tail conformations are preferred over others. We capture this effect in the Pink model by randomly labeling a fraction p of the lattice sites as "pinning sites". At the pinning sites, the corresponding lipid chain is fixed into the ground state conformation. This extension is trivially incorporated into our MC simulations: we simply do not apply the MC move to pinning sites. We specialize to DSPC, using the "standard" Pink parameters of Tables I and II. In Fig. 4(a), we show ln P (A|T, Π) at T = 291 K and Π COEX = −18.9 dyn/cm in the absence of quenched disorder (p = 0). At this temperature, which is well below T c , the main transition is strongly first-order. Consequently, there is a significant line tension σ between gel and fluid domains; the latter is related to the free energy barrier [34,52] ∆F ≡ k B T ln (P max /P min ) = 2σL , indicated by the vertical double-arrow in Fig. 4(a). Here, P min is the value of P (A|T, Π) at the minimum "between the peaks", while P max denotes the average peak value. The physical motivation for Eq. (7) is that, for crosssectional areas "between the peaks", the bilayer reveals a coexistence between two slab domains where the total interface length equals 2L [ Fig. 4(b)]. For DSPC, and assuming the lattice constant to be 1 nm, we obtain σ ∼ 1.1 pN, which is compatible with experimental values [53]. Next, we consider a DSPC bilayer with a fraction p = 0.03 of the lattice sites marked as "pinning sites". In Fig. 5(a), we show distributions P (A|T, Π) for T = 291 K obtained for 4 different random positions (samples) of pinning sites. Even though the temperature is the same as in Fig. 4(a), a unique double-peaked distribution can no longer be identified. In contrast, P (A|T, Π) is strongly sample dependent, and a multitude of rather exotic shapes is revealed. This behavior is characteristic of systems that belong to the universality class of the FIG. 6. The (disorder-averaged) Binder cumulant U1 as a function of temperature T for DSPC with a fraction p = 0.03 of pinning sites. In contrast to Fig. 2(b), an intersection of the curves for different L can no longer be identified. Instead, as the system size L increases, U1 → π/2, indicative of the one phase region. For each system size, the disorder average comprised 200 samples of pinning sites.
2D random-field Ising model (2D-RFIM) [21]. Hence, by introducing the pinning sites, we have changed the universality class of the Pink model from ordinary 2D Ising toward 2D-RFIM (the pinning sites essentially correspond to a field of infinite strength acting at random locations).
There are two features of the 2D-RFIM universality class that are remarkably consistent with experimental results for the main transition in supported membranes. First of all, 2D-RFIM universality implies the absence of macroscopic coexistence between gel and fluid domains [25][26][27][28]. Indeed, inspection of simulation snapshots [ Fig. 5(b)] reveals an equilibrium multi-domain structure, that is highly anisotropic, strongly resembling experimental AFM images [50]. A second (related) feature is that the 2D-RFIM has no true phase transition in the thermodynamic limit. In finite systems, there may be signs of a transition (or even several transitions; note that some of the distributions in Fig. 5(a) are triplepeaked), but they will be "smeared" over a wide temperature range, and do not persist in the thermodynamic limit. Precisely this behavior has also been reported in experiments [48,50]. Simulation evidence that the pinning sites prevent a sharp transition in the thermodynamic limit follows from the (disorder-averaged) Binder cumulant denotes an average over different samples of pinning sites. As shown in Fig. 6, U 1 → π/2 as L increases, consistent with only a single phase.

V. CONCLUSION
In this paper, the main phase transition in singlecomponent phospholipid membranes was investigated using the Pink model. Our simulations of the pure membrane (i.e. without quenched disorder) confirm the formation of macroscopic gel and fluid domains below a critical temperature T c , and at the coexistence pressure Π COEX . We also demonstrated that, using the accepted values of the Pink model parameters, T c falls below experimentally measured transition temperatures. This explains why no phase transition was detected at the experimental transition temperature in the simulations of Ref. 18. To resolve this issue, we propose that the strength of the van der Waals coupling in the Pink model be increased. By using the values proposed in this work, T c of the Pink model in the thermodynamic limit coincides with the experimental main transition temperature. In addition, finite-size scaling was applied to confirm the universality class of the critical point in the Pink model, which was shown to be 2D Ising. This result holds irrespective of the number of conformational states q (as long as q ≥ 2, of course). Hence, to capture the generic features of membrane phase behavior, a highly-detailed model is not always needed (which is consistent with the findings of Ref. 24).
We have also used the Pink model to describe the main transition in the presence of quenched disorder, which may arise in case the membrane is deposited on a rough support. Assuming that this induces regions in the membrane where certain tail conformations become preferred, the universality class changes toward that of the 2D random-field Ising model. In the presence of quenched disorder, the Pink model reveals a stable multi-domain structure, and the absence of a sharp transition; these findings are indeed consistent with some of the experimental observations. Although the Pink model (both the pure version and the one containing quenched disorder) seems well suited to describe the main transition, we wish to end with a warning. By using the Pink model, one inevitably casts the main transition into the Ising universality class. This may not be entirely appropriate, as the main transition is essentially a melting transition leading to the forma-tion of nematic chain order. A liquid-crystal model may therefore be more suitable, such as the Maier-Saupe approach followed in Ref. 54. If, indeed, the main transition occurs close to a critical point [43,55] it may well be necessary to replace the discrete set of states of the Pink model by a continuous one [56]. In that case, one enters the regime of Heisenberg-type models (which, provided certain conditions are met, do support 2D phase transitions [57][58][59]). The investigation of the main transition in terms of such a continuous model could be an interesting topic for future work.