Finsler Geometry Modeling of Phase Separation in Multi-Component Membranes

A Finsler geometric surface model is studied as a coarse-grained model for membranes of three components, such as zwitterionic phospholipid (DOPC), lipid (DPPC) and an organic molecule (cholesterol). To understand the phase separation of liquid-ordered (DPPC rich) Lo and liquid-disordered (DOPC rich) Ld, we introduce a binary variable σ(=±1) into the triangulated surface model. We numerically determine that two circular and stripe domains appear on the surface. The dependence of the morphological change on the area fraction of Lo is consistent with existing experimental results. This provides us with a clear understanding of the origin of the line tension energy, which has been used to understand these morphological changes in three-component membranes. In addition to these two circular and stripe domains, a raft-like domain and budding domain are also observed, and the several corresponding phase diagrams are obtained.


Introduction
Membranes of multiple components, such as 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), dipalmitoylphosphatidylcholine (DPPC) and cholesterol, are receiving widespread attention because of their applications in many fields of science and technology, and numerous studies on the morphological changes have been conducted [1][2][3][4][5][6]. In these membranes, morphological changes are induced by a phase separation. Indeed, the phase separation causes domain formation and domain pattern transition between the liquid-ordered (L o ) and the liquid-disordered (L d ) phases. This domain pattern transition accompanies the morphological changes, such as the two circular domains, the stripe domain, the raft domain and the so-called budding domain [5,6]. The multiplicity of components, as in a glass transition [7], is essential for such a variety of morphologies. To date, these morphologies have been studied on the basis of the line tension energy [3,4] in the context of the Helfrich-Polyakov (HP) model for membranes [8,9]. The line tension energy is defined on the domain boundary and has an important role in the morphological changes [3,4].
However, the origin of the line tension energy is not well understood. In fact, it is unclear what type of internal structure is connected to the line tension energy until now. The problem that should be asked is where the line tension energy originates. Therefore, in this paper, we clarify and discuss the microscopic origin of the line tension energy.
To understand the origin of the line tension energy, we introduce a new degree of freedom σ(= ±1) to represent the L o and L d phases. The Ising model Hamiltonian, which we call aggregation energy, for the variable σ is included in the general HP model Hamiltonian, where the "general" HP model refers to the HP model with a nontrivial surface metric g ab ( δ ab ). Note that the general HP model can be discretized on triangulated surfaces and becomes well defined only when it is treated in the S 1 = √ gd 2 xg ab ∂r ∂x a · ∂r ∂x b , where g is the determinant of the 2 × 2 matrix g ab of the metric function and g ab is its inverse [14]. The symbol n denotes a unit normal vector of the surface. Both S 1 and S 2 are conformally invariant. The conformal invariance is a property in which a scale change g ab (x) → f (x)g ab (x) is not reflected in both S 1 and S 2 for any positive function f . Two metrics g and g are called "conformally equivalent" if a function f (x) exists, such that g ab = f (x)g ab [14]. For the case where g ab (x) is given by the Euclidean metric g ab = δ ab (or the induced metric g ab = ∂ a r · ∂ b r), the surface shape r in 3 is treated from the perspective of statistical mechanics. These are the HP model [8,9] corresponding to polymerized membranes, and the HP model and the Landau-Ginzburg model [15] have been thoroughly investigated [16][17][18][19][20][21][22].

Discrete Model
First, in this subsection, let us introduce a new degree of freedom σ, which has only two different values (σ = ±1), on the triangulated lattice (see Figure A1 in Appendix A). We assume that the variable σ i is defined on the triangle ∆ i , and moreover, the values of σ i correspond to two different phases, namely the liquid-ordered (L o ) and the liquid-disordered (L d ) phases, such that: This definition of σ implies that every triangle is labeled by the value of σ, and therefore, σ represents the phase (or domain) to which the triangle ∆ belongs.
We now introduce a discrete Hamiltonian for multi-component membranes. The technical details of the discretization of the continuous Hamiltonian S 1 and S 2 introduced in Subsection 2.1 are described in Appendix A.1, and the discrete expressions for S 1 and S 2 are given in Equation (A8) in Appendix A.1. Using these S 1 and S 2 , we have the total Hamiltonian S, such that: where S (r, σ) denotes that the Hamiltonian depends on the variables r(∈ 3 ) and σ. The three-dimensional vector r denotes the vertex position of the triangulated lattice. The energy λS 0 is called the aggregation energy. When λ → 0, the variable σ becomes random, and this random configuration simply corresponds to the coexistence phase, where L o and L d are not separated. Conversely, when λ becomes sufficiently large, two neighboring σ's have the same σ, and this configuration corresponds the phases where L o and L d are separated. As described above, the second and third terms S 1 and S 2 in S are the discrete Hamiltonians corresponding to the continuous ones introduced in Section 2.1. The coefficient κ of S 2 is the bending rigidity and has units of [1/k B T], where k B and T are the Boltzmann constant and the temperature, respectively. In this paper, we assume that k B T = 1. The symbol n i in S 2 expresses a unit normal vector of the triangle i. The symbols γ ij (σ) and κ ij (σ) denote that γ ij and κ ij depend on the variable σ, and this dependence arises from an interaction between σ and r. The interaction between σ and r is defined by the function ρ (see Equation (A1) in Appendix A.1), such that: where c is a parameter that should be fixed at the beginning of the simulations. From this definition of ρ(∆) and Equation (A9) in Appendix A.1, we have (see Figure 1): correspond to the bonds represented by the duplicated lines.
These expressions represent how the effective surface tension γ ij and bending rigidity κκ ij depend on the position of the bond i j, which is one of the three domain boundary bonds (L o , L o ), (L o , L d ) and (L d , L d ). The symbol (L o , L d ) refers to the bond shared by the two neighboring triangles of the L o and L d phases (see Figure 1). Note that only (L o , L d ) corresponds to the bond on the domain boundary, and the other two correspond to the bonds inside the domains L o and L d . From the expressions of γ ij and κ ij in Equation (5), we understand that the dependence of γ ij and κ ij on the domains and their boundary is automatically determined. Thus, this expression is one of the most interesting outputs of the model in this paper. The values of γ ij and κ ij depend on the parameter c, which is an input parameter.
In Figure 2, γ ij (= κ ij ) for (L o , L o ), (L o , L d ) and (L d , L d ) are plotted as functions of c in the region 1 ≤ c. The expressions of γ ij and κ ij in Equation (5) are symmetric under the exchange c ↔ 1/c, and therefore, we use the value of c rather than 1/c to represent γ ij and κ ij . The curve of γ ij (= κ ij ) against c is almost linear except for the region c 1. The dashed vertical lines in the figure correspond to c = 5 and c = 8.37, which are assumed in the simulations.  The fluid surface model is defined by the sum over all possible triangulations T in the partition function, such that: where the prime in denotes that the center of mass of the surface is fixed at the origin of 3 to protect the surface translation. The dynamical triangulation, denoted by T , is performed using the bond flip technique [23][24][25][26][27][28]. Due to this bond flip, the vertices can freely diffuse over the surface, where two neighboring triangles merge and split into two different ones, and the total number of triangles remains unchanged in this process. Therefore, not only the vertices, but also the triangles diffuse over the surface.
Note that the metric variable, or in other words, the function ρ, is not summed over (or integrated out) in Z; hence, strictly speaking, it is not a dynamical variable. However, the metric variable ρ is effectively considered as dynamical in the sense that ρ changes its value on the surface due to the diffusion of triangles.
Moreover, note that the aggregation energy λS 0 simply corresponds to the line tension energy in [3,4]. Indeed, the energy 1−σ i · σ j at the bond i j has a non-zero positive value only when the bond is on the domain boundary between L o and L d . More precisely, S 0 is proportional to the total number of bonds that form the domain boundary because the mean bond length is constant (or non-zero finite) on the boundary.
We comment on the reason why λS 0 is considered as the line tension energy in more detail. First, the fact that the mean bond length becomes constant is understood from the scale-invariant property of the partition function Z in Equation (6). Indeed, we have S 1 /N = 3/2 [29]. It is easy to see that this relation is satisfied: Z is independent of the scale change r → αr for arbitrary α ∈ , and therefore, we , S(αr) = λS 0 +α 2 S 1 +κS 2 , we have S 1 /N = 1.5 for sufficiently large N. The relation S 1 /N = 3/2 means that γ ij 2 ij is constant. This implies that 2 ij , and hence, ij becomes constant. This constant i j varies depending on whether the bond i j belongs to L o , L d or the boundary between L o and L d because the coefficient γ ij varies depending on these domains and the domain boundary, as in Equation (5). For this reason and because γ ij 2 ij = constant, we understand that the bond length on the boundary between L o and L d becomes well defined (or non-zero finite). Importantly, the mean bond length is expected to be finite, although it fluctuates around the mean value, and the mean value itself varies depending on the domains or the domain boundary. Therefore, λS 0 is considered to be an extension of the line tension energy because S 0 is proportional to the length of the phase boundary if the two phases are clearly separated as the domains L o and L d at least.
The remaining problem to be clarified is how the domain boundary is formed on the triangulated surfaces. During experiments, the area fraction of L o (and L d ) is fixed [3]. Hence, in our model, the total number of triangles N o T for σ = 1 (and N d T for σ = −1) is fixed, where the total number of triangles: is also fixed to be constant (because N T = 2N −4 and the total number N of vertices is fixed). The relation between the area fraction of L o and the fraction of N o T will be described in the next section. Another constraint imposed on the triangles in our model is that the value of σ i of triangle i remains unchanged for all i during the simulations. Therefore, the triangles themselves have to diffuse over the surface to form the L o and L d domains. This triangle diffusion is numerically possible on the dynamically-triangulated surfaces, which are called triangulated fluid surfaces, via the Monte Carlo (MC) technique with dynamical triangulation, as described above [23][24][25][26][27][28].
The function ρ(∆) in Equation (4) characterizes the difference between the phases L o and L d of ∆, and these two different phases are labeled by the variable σ(∆) as in Equation (2). Therefore, the model in this paper is limited to membranes with two-component domains; however, the modeling technique is applicable to membranes with multi-component domains. Here, we comment on how to extend the model to a n-component model. To extend the two-component model, we have to define the value of ρ(∆) for the n-component model, such that ∆ ∈ L i (1 ≤ i ≤ n) (see Equation (4)), where {L 1 , L 2 , · · · , L n } is the set of domains assumed. In this case, the variable σ(∆) should be n components, and therefore, the corresponding energy term λS 0 in Equation (3) should also be extended. The Hamiltonian of the n-states Potts model, for example, can be used for S 0 . Hamiltonians of continuous models, such as the Heisenberg spin model, can also be assumed for S 0 , where the continuous variable σ(∆) should be connected in one-to-one correspondence with the n-component function ρ(∆) (see Equation (4)). In these n-component models, the energy λS 0 is still expected to play the role of line tension energy between two different domains, because λS 0 becomes zero (nonzero) if the phases of two neighboring triangles are identical to (different from) each other. The parameters κ ij and γ ij are given by the same expression in Equation (A9); however, the final expression of these parameters in the n-component model are in general different from those in Equation (5) because of the dependence of the parameters on the definition of ρ i . Indeed, the parameters κ ij and γ ij in the n-component model will be different from those determined by the single parameter c in the two-component model.

Monte Carlo Technique
The canonical Metropolis technique is used [30,31]. The vertex position r is updated, such that r = r+δr. The symbol δr denotes a random three-dimensional vector in a sphere of radius R. The new position r is accepted with probability Min[1, exp(−δS)], where δS = S(new)−S(old). The radius R of the small sphere is fixed such that the acceptance rate for the update of r is approximately equal to 50%.
The triangulation T is updated using the bond flip technique, as described in the previous section [23][24][25][26][27][28]. We use the same technique used in [23][24][25][26][27][28], except for the following constraint. In the bond flip, the two neighboring triangles of the bond change to a new pair of triangles, such that the More precisely, if the two triangles have the same value of σ prior to the bond flip, then the new values of σ for the new triangles are fixed to be the same as the old one. However, if the values of σ are different from each other before the bond flip, then the new values are also fixed randomly to be different. Only through this process is the variable σ updated. Due to this update of σ through the dynamical triangulation, the function ρ changes, and hence, a domain structure of L o (or L d ) is formed on the surface.
We comment on the relation between the fraction φ 0 and the area fraction of L o . As described in the previous section, the mean triangle areas a o and a d in the domains L o and L d are constant because of the scale invariance of Z. Therefore, the area fraction of L o can be written as However, the area fraction of L o is not always reflected in the fraction φ 0 if a o a d .
The initial configuration for the simulations is fixed to be the random phase, where the L o (or L d ) triangles are randomly distributed on the surface under a constant ratio φ 0 . This random state corresponds to the two-phase coexistence configuration.
A single Monte Carlo sweep (MCS) consists of N updates of r and of N updates for the bond flips. The total number of MCS that should be performed depends on the parameters; it ranges from approximately 1 × 10 8 to 8 × 10 8 . The numbers of MCS for almost all simulations are 2 × 10 8 ∼ 3 × 10 8 . The simulations at the phase boundaries are relatively time consuming in general because the domain structure and, hence, the surface shape change very slowly at these boundaries. The total number of vertices N is fixed to N = 5762 in this paper.

Simulation Results
Two types of models, which are denoted as Model 1 and Model 2, are simulated. The Gaussian energy S 1 = ij 2 ij of the canonical surface model is assumed for Model 1. From this assumption, the effective surface tension γ ij for Model 1 is γ ij = 1. Model 2 is the same as the one introduced in Section 2.2. The Gaussian energy S 1 and the parameters γ ij and κ ij for Model 1 and Model 2 are presented in Table 1.
In Model 1, the surface shape is influenced only by κ ij because γ ij = 1. In contrast, in Model 2, the coefficient γ ij influences the surface size because S 1 in Equation (3) has the unit of length squares. Indeed, as described in Section 2, from the scale invariance of Z, we have S 1 /N = 3/2 [29], and therefore, 2 ij deviates from the constant expected from this relation if the constraint γ ij = 1 is not imposed on γ ij . For example, if γ ij is large (small), then 2 ij becomes small (large). Therefore, due to this dependence of 2 ij on γ ij , the size of the triangles in Model 2 depends on the domains. By contrast, there is no dependence of 2 ij on the domains in Model 1, where γ ij = 1 over the entire surface. The input parameters for the simulations are λ, κ, c, and φ 0 , where c is the value of the function ρ in Equation (4) and determines γ ij and κ ij . The parameter φ 0 defined by Equation (8) is identical to the area fraction in Model 1, whereas it differs from the area fraction in Model 2 because the triangle area is not uniform in Model 2. More precisely, the mean triangle area in the L o domain is different from that in the L d domain in Model 2. In Table 2, we show the parameters assumed in the simulations. The values of κ ij corresponding to the input c are listed in Table 3. Table 2. The input parameters λ, κ, c and φ 0 for the simulations. Table 3. The input parameter c automatically defines the values of κ i j and γ i j , where γ i j = 1 for Model 1 and γ ij = κ i j for Model 2. We first show a phase diagram on the λ − φ o plane in Figure 3. The parameters κ and c are fixed to κ = 7 and c = 5, as shown in Table 1, and λ is varied in its relatively small region. The dots (•) are the data points where we perform the simulations to construct the phase diagram. We find that the two phases L o and L d are not separated in the region λ < 0.1; the domain pattern is random; and the surface is almost spherical, as observed in the snapshots. In contrast, in the region λ > 0.1, L o and L d are clearly separated, and the two circular domains and the stripe domain appear. The domain structure depends on the value of φ 0 , and the corresponding surface morphology appears to be almost discontinuously separated on the phase diagram. We observe that the two circular domains change to the stripe domain as the fraction φ 0 , which is identical with the area fraction of L o , increases for constant λ. This result is consistent with the experimental results reported in [3], where the area fraction of L o is changed. The two circular domains and the stripe domain correspond to the L o phase, where κ ij is higher than those of both the L d domain and the boundary, as shown in Table 1 Next, to show the dependence of the surface size on the parameters, we define semi-axis lengths D 1 , D 2 and D 3 of the surface such that D 1 > D 2 > D 3 as in Figure 5. D 1 and D 2 , D 3 correspond to the major and minor axes, respectively. The surface of the stripe domain corresponds to the so-called prolates, where D 1 > D 2 D 3 is expected. It is also expected that D 1 D 2 > D 3 in the so-called oblates, which corresponds to the surface shape of the two circular domains. Therefore, the surfaces with the stripe and two circular domains can be distinguished by the minor axis D 2 .
> > We plot D 2 vs. φ 0 in Figure 4a, where λ = 0.2. As shown, D 2 discontinuously changes against φ 0 at the phase boundary between the two circular and stripe domains. From the plot of D 2 vs. λ in Figure 4b, we also observe that D 2 discontinuously changes at the same phase boundary. The bending energy S 2 /N B in Figure 4c also discontinuously changes at this boundary, and this result indicates that this morphological change is considered as a first-order transition. However, note that the change of the morphology at this phase boundary is relatively smooth. In fact, one circular domain surface, which is not shown as a snapshot in Figure 3, can be observed at the boundary. This implies that the stripe domain surface and one circular domain surface have the same bending energy, or in other words, the bending energy is degenerate. Additionally, note that the phase boundary between the two circular and random domains appears to be continuous. This means that the shape of the two circular domain surfaces continuously changes to the random domain surface. At this phase boundary, the surface shape continuously changes from pancake to sphere.

Model 2
In Model 2, not only κ ij , but also γ ij depend on the domain (or the domain boundary) whether it is L o or L d . For this reason, the area of the triangles in the L o domain becomes considerably smaller than that in the L d domain. Therefore, the fraction φ 0 does not reflect the area fraction of L o in this case. In fact, it is easy to see that the area fraction of L o in the snapshots at φ 0 = 0.9 in Figure 6 is much smaller than 90%. Nevertheless, the phase diagram on the λ − φ 0 plane in Figure 6 appears almost the same as that in Figure 3. The only difference between the two phase diagrams is the appearance of one circular domain phase, denoted by "one circ" in Figure 6. This one circular phase is stable, where "stable" means that the surface domain remains unchanged against a small variation of the parameters inside the phase boundary. This is in sharp contrast to the one circular domain surfaces observed at the region close to the boundary between the two circular and stripe domains because these one circular surfaces are very sensitive to the parameter variation and, hence, "unstable". The shape of the one circular surface in the one circular region is almost spherical, such as the one shown in Figure 6, and this result is in contrast to the result in [4], where the one circular phase is separated into two phases: the prolate and oblate phases. One possible reason for why only a spherical surface appears in the one circular domain in Figure 6 is because the L o domain is hardly bent due to the high ratio κ ij (L o , L o )/κ ij (L d , L d ) = 4.24, which is slightly larger than the one 1 < κ ij (L o , L o )/κ ij (L d , L d ) < 3 assumed in [4]. The parameters assumed on this plane are κ = 10 and c = 8.37, which are listed in Table 2. The simulations are also performed on the λ − φ 0 planes for larger κ, such as κ = 15 and κ = 20, and with c = 8.37. The phase diagrams obtained in these simulations are (not shown) relatively close to those shown in Figure 6; however, surfaces with three or four circular domains appear in the lower λ region in the two circular domain phase. The bending energy κS 2 of the three or four domains is lower than that of the two circular domain; moreover, the aggregation energy λS 0 of these multi-circular domains is larger than that of the two circular domain. These are the reasons for the appearance of the three or four domains only in the relatively small λ region in the simulations with relatively large κ.
To observe the variation of the surface size at the phase boundaries, we calculate the size D 2 on the dashed lines in Figure 6 and plot them in Figure 7a,b. We determine that D 2 discontinuously changes against φ 0 and λ at the phase boundaries, similar to that in Model 1 shown in Figure 4a,b. Moreover, the phase boundary is also not as clear because of the same reason as that for Model 1. In fact, the surface shape at the phase boundary between the two circular and stripe domains is not always stable in Model 2, similar to that in Model 1. Figure 7c also shows that S 2 /N B discontinuously changes; however, the gap is very small, and these two phases are hence separated by a weak first-order transition. The phase boundary between the two circular and the random domains is also expected to be continuous in Model 2. The boundaries of one circular to two circular and one circular to stripe are also not as clear, and the boundary of one circular to random is continuous.
Another difference between Model 1 and Model 2, other than the appearance of the stable one circular domain, is the raft-like domain and the budding domain. More precisely, the budding domain can also be seen in Model 1; however, it is more clear in Model 2. The phase diagram of Model 2 on the κ−φ 0 plane is drawn in Figure 8. The parameter λ is fixed to λ = 3, which is relatively large compared to the previous one assumed in the simulations for Figures 3 and 6. Consequently, the energy λS 0 , which is the line tension energy, at the phase boundary between L 0 and L d becomes large in the region where κ is relatively small. This is the reason why the budding domain appears on this κ−φ 0 plane in Figure 8. Note that the budding domain in some of the budding surfaces goes inside the surface, and some of them self-intersect because no self-avoiding interaction is assumed. Figure 8 also shows that the raft-like domain is stable in the relatively large κ region, where the surface hardly deforms. The reason why the raft domains, which are multi-circular domains, appear only at the region of small φ 0 is because the multi-circular domains are more energetically favorable than the stripe domain. Indeed, the effective bending rigidity κκ ij and, hence, κS 2 become very large on the large connected L o domain, such as the stripe domain, where the line tension energy λS 0 is relatively small. Note that S 0 has a non-zero positive value only on the boundary bonds between L o and L d , while S 2 has a non-zero value on all of the bonds. Moreover, note that the boundary length between L o and L d becomes longer (shorter) if the total number of circular domains increases (decreases), whereas the areas of L o and L d remain constant and are independent of the total number of L o domains.
Finally, we show that S 1 /N satisfies the relation S 1 /N = 1.5, which is expected by the scale invariance of Z in Equation (6) [29]. As described in Section 2, the bond length is expected to be well defined in the sense that the mean bond length is constant on the surface, although this constant varies depending on the domains or the domain boundary to which the bond belongs. The data in Figure 9a,b are obtained on the dashed lines in Figure 3, and those in Figure 9c,d are obtained on the lines in Figure 6. These data shown in Figure 9 indicate that the simulations including the energy discretization are successful.    Figure 3 (Figure 6).

Summary and Conclusions
We have studied the phase separation of the three-component membrane with DPPC, DOPC and cholesterol using a Finsler geometry (FG) surface model. The FG model is obtained from the Helfrich-Polyakov (HP) model for membranes by replacing the surface metric with a general one g ab δ ab , which can be called the Finsler metric. In other words, we have extended the HP model to explain the morphological changes of the three-component membranes in the context of FG modeling. This new model includes a new degree of freedom σ, which represents the liquid-ordered (L o ) and liquid-disordered (L d ) domains. The results obtained from Monte Carlo (MC) simulations are consistent with the experimental results that have been reported in the literature. We confirm the phase separation of the L o and L d domains on the surface and that the surface shows a variety of morphologies, such as the two circular domain, the stripe domain, the raft domain and the budding domain.
The line tension energy, which has been used for understanding the morphological changes, simply corresponds to the aggregation energy term λS 0 in our model. Indeed, the value of S 0 is only the total number of bonds on the boundary between L o and L d in our new model. Moreover, the fact that λS 0 is simply the line tension energy implies that the line tension originates from the interaction between the domains because the interaction between the variables σ in S 0 describes the interaction between the domains. This interaction is closely connected to the property of the new model that the surface strength, such as the surface tension and the bending rigidity, is dependent on the bond position on the surface. This property arises from the interaction between σ and r introduced via the Finsler metric.

Abbreviations
The following abbreviations are used in this manuscript: To obtain the discrete model from the continuous surface model introduced in Section 2.1, we assume that the surface is triangulated in 3 . The Hamiltonian is defined on the triangulated surfaces, which are composed of three simplexes, such as vertices, bonds and triangles. Thus, all physical quantities, including Hamiltonians and the metric function, are defined on these simplexes labeled by integers. For example, the vertex position r i is defined at the vertex i; the bond length ij is defined on the bond i j; and the elements of g ab are defined on the triangle ∆. Note that the variable r is considered as a mapping from the parameter space M to 3 .
We start with the discrete metric g ab , such that: where ρ is a function on a triangle ∆(⊂ 3 ) (see Figure A1a). More precisely, the elements of g ab are functions on the triangle ∆ M (⊂ M), where M is the aforementioned two-dimensional space M (independent of 3 ). We assume that M is also triangulated by the triangles ∆ M . On this ∆ M , an orthogonal coordinate can be taken for any one of three vertices [14]. For this reason, the metric g ab can be diagonalizable. The inequality ρ > 0 in Equation (A1) is necessary for the positivity of the bond length. This metric depends only on x and is independent of y, and hence, it simply corresponds to the Riemannian metric. Indeed, this metric in Equation (A1) comes from the most general one, such as g ab = E F F G , with the functions of E > 0, G > 0, EG−F 2 > 0. By assuming F = 0, we have , where ρ 2 = G/E and the symbol " " denotes conformally equivalent. Note that this expression of g ab depends on the local coordinates on ∆, and therefore, the expression of g ab implicitly depends on the vertex of ∆ because the coordinate origin is located on one of the vertices of ∆. In the discrete models that have been studied thus far, the Euclidean metric g ab = δ ab (or the induced metric g ab = ∂ a r · ∂ b r) is always assumed as mentioned above, and it has been reported that the model for polymerized membranes undergoes a discontinuous or a continuous transition between the crumpled phase and the smooth phase [32][33][34][35]. Let the vertex r 1 of the central triangle ∆ in Figure A1b be the local coordinate origin in this ∆. By replacing: we have: where n i (i = 1, 2, 3) are the unit normal vectors shown in Figure A1b. The symbol ij (= ji ) is defined by ij = |r j −r i |. Note that the unit normal vector also represents the surface orientation; indeed, n 0 is defined by n 0 = 12 × 13 /| 12 × 13 |, for example. We have three possible coordinate origins in the triangles. For this reason, S 1 and S 2 can be symmetrized by including the terms that are cyclic permutations, such as 1 → 2, 2 → 3, 3 → 1 for ij , n i and ρ i . Summing over all possible terms and multiplying by a factor of 1/3, we obtain: where ρ i (i = 1, 2, 3) are defined on the triangle ∆. The reason for why these three different functions ρ i are included is because the expression for g ab generally depends on the local coordinate as mentioned above. More precisely, ρ i is the element of g ab on ∆ where the coordinate origin is located at vertex i. For arbitrary g ab , we always have the metric of the form in Equation (A1) by the same procedure as described above.
Here, we further simplify the model by assuming that: Thus, we have the expressions: Replacing the sum over triangles ∆ with the sum over bonds ij , we obtain: Note that S 1 and S 2 defined by the sum over triangles ∆ in Equation (A6) are exactly the same as those defined by the sum over bonds ij in Equation (A7), and the difference is only in their expressions. Additionally, note that the suffixes i, j of ij in Equation (A7) denote the bond i j, whereas those of ρ i and ρ j denote the two neighboring triangles i and j of the bond i j. Thus, we finally have: with: where the irrelevant numerical factor 1/3 is replaced by 1/4 in the final expressions for S 1 and S 2 . Note that γ ij (κκ ij ) can be called the effective surface tension (effective bending rigidity) on the bond between vertices i and j. It must be emphasized that the quantities γ ij and κ ij are independent of the bond direction, or in other words, these are symmetric under the exchange of i and j, and for this reason, γ ij and κ ij are considered as the quantities defined on the bond i j. Indeed, from the expression given in Equation (A9), we have: Therefore, the physical quantities γ ij 2 ij in S 1 and κ ij (1−n i · n j ) in S 2 of Equation (A8) are well defined in the sense that these quantities are symmetric under the exchange of i j. The reason why we need this symmetry in the physical quantities γ ij 2 ij and κ ij (1−n i · n j ) is because these quantities correspond to the energies for the expansion and bending of the surface at the bond i j, and these energies are independent of the bond direction, such as the one from i to j or the reverse. Thus, the symmetry property in Equation (A10) allows us to call γ ij and κκ ij the effective surface tension and the effective bending rigidity on the bond i j, respectively. However, as we will show in the next subsection, γ ij and κ ij are not symmetric in general (see Figure A1c). Moreover, note that this problem of whether γ ij and κ ij are symmetric arises only when γ ij and κ ij depend on the functions ρ i and ρ j on the two neighboring triangles i and j. This is in sharp contrast to the case where γ ij and κ ij depend only on the quantity defined on the vertices [10], where γ ij and κ ij are always symmetric. This exchange symmetry/asymmetry reflects the orientation symmetry/asymmetry, which will be discussed in the next subsection.

Appendix A.2. Finsler Geometry Model
In this subsection, we show that the discrete surface models constructed above are well defined only in the context of Finsler geometry modeling [10]. For this purpose, we should first remind ourselves of the fact that the symmetry properties in Equation (A10) can be observed in the model only under the condition of Equation (A5). This symmetry is not present in the model of Equation (A4). To show the breakdown of the symmetry of Equation (A10) in the model of Equation (A4) in more detail, we replace the sum over triangles ∆ of S 1 and S 2 in Equation (A4) with the sum over bonds ij before the condition of Equation (A5) is imposed. In this new expression of S 1 , which is expressed by the sum over bonds ij , the Gaussian bond potential for the bond 12, which is shared by the triangles ∆ + and ∆ − as in Figure A2a, for example, is given by: In this expression, the former half S + 1 ( 12 ) = (1/3) ρ + 1 +1/ρ + 2 2 12 is the contribution from ∆ + , and the latter half S − 1 ( 12 ) = (1/3) ρ − 2 +1/ρ − 1 2 12 is the contribution from ∆ − . However, it is clear that γ 12 , and hence, S 1 ( 12 ) in Equation (A11) are not symmetric under the change of surface orientation. In fact, we have:S for the opposite orientation (see Figure A2b). In Equation (A12), we write the coefficient of 2 12 by γ 21 because it is obtained from γ 12 in Equation (A11) by exchanging the suffixes 1 and 2. It is also easy to show that γ 12 γ 21 , and hence, S 1 ( 12 ) are not always identical toS 1 ( 12 ) in general.
Thus, we find that the asymmetry γ 12 γ 21 means that S 1 ( 12 ) is not invariant under the orientation exchange. From this, we can see that the Gaussian bond potential energy (and also the bending energy) of the bond 12 of one surface configuration differs from that of the opposite orientation configuration. However, we have no reason for the difference in S 1 for two surfaces with different orientations. Thus, the model defined by Equation (A4), which is orientation asymmetric, is ill defined in the context of conventional surface modeling. Moreover, we have to remark that the model defined by Equation (A8), which is orientation symmetric, is also ill defined. The reason for this ill definedness is that the bond length squares calculated with ρ + in ∆ + is not always identical to the one calculated with ρ − in ∆ − in Figure A2a, where ρ ± 1 = ρ ± 2 in the model of Equation (A8). Indeed, the metric on ∆ + is given by g ab = 1/ρ + 0 0 ρ + , where the coordinate origin is at the vertex 1. Then, we have the bond length squares (1/ρ + ) 2 12 for the bond 12 with respect to the metric g ab , and changing the vertex origin from 1 to 2, we also have ρ + 2 12 . Thus, summing over these two expressions without the coefficient 1/2, we have (1/ρ + +ρ + ) 2 12 for the bond length squares. Through the same procedure, we have (1/ρ − +ρ − ) 2 12 from ∆ − . These two square lengths of the bond 12 must be the same. However, we have: because ρ + ρ − in general. The edge length of triangles should uniquely be given as the basic requirement even in the discrete models. Therefore, in a model construction on the triangulated lattices, we always obtain an ill-defined discrete model if we start with an arbitrary Riemannian metric in which the elements are defined on the triangles. Note that the bond "length" used here is the length with respect to g ab on ∆ ± and is different from the Euclidean bond length ij (also note that g ab is simply a Riemannian metric at this stage). However, these ill-defined models in Equations (A4) and (A8) become well defined in the context of Finsler geometry [10][11][12]. In this context, the bond length calculated with g ab on the triangle ∆ + can be considered as the direction-dependent length from 1 to 2, and the one calculated with g ab on ∆ − can be considered as the length from 2 to 1 ( Figure A2a). Therefore, the inequality in Equation (A13) is satisfied in general. Moreover, the aforementioned quantities S + 1 ( 12 ) and S − 1 ( 12 ) in Equations (A11) and (A12) are meaningful because these quantities are also considered as direction dependent in the Finsler geometry context ( Figure A1c).