Curvature-Induced Spatial Ordering of Composition in Lipid Membranes.
The biological lipid bilayer membrane, or in short "biomembrane," is a fundamental building block of the cell. It forms the barrier that separates the interior of the cell from its surroundings but is also responsible for almost all interaction of the cell with its environment, including transport, adhesion, regulation, transduction, and signaling [1-5]. The diverse functionality of the biomembrane is achieved by a seemingly simple structure, two layers that are primarily made from lipid molecules and also some integral proteins, cholesterols, and other functional molecules [6, 7]. This molecular structure of the biomembrane gives rise to the so-called "fluidity" of the membrane ; that is, its constituent molecules are able to move relatively easy within the membrane, which resists bending and stretching but not shear . Consequently, biomembranes have a dynamic structure in the sense that their molecular arrangement (local composition) can change with conditions. For example, depending on temperature (and/or osmotic pressure, acidity, etc.) the biomembrane may possess a uniform mixture of its components or it may segregate into different phases, which are rich in specific constituent and possess different mechanical and chemical properties [10-16].
The fluidity of the biomembrane combined with its spatial heterogeneity brings about a unique coupling between shape (geometry) and composition. For example, lipid phases that possess high bending stiffness highly favor regions with small (magnitude) curvature . Also, the three-dimensional molecular shape of some lipids and proteins results in a nonzero spontaneous curvature that affects the geometry of the biomembrane in their neighborhood. This two-way coupling between shape and composition means that deformations exhibited by biomembranes are strongly influenced by their heterogeneous composition, while the spatial ordering of composition is modulated by the geometry of the membrane [16,18-23].
In the last two decades, much effort has been invested into understanding the consequences of the coupling between shape and composition in biomembranes. Theoretical models have generalized uniform composition models [24-27] to account for multicomponent or multiphase membranes [28-35]. Features such as equilibrium configurations, stability [36-39], interaction with the cytoskeleton [40-42], formation of lipid rafts, anisotropy of the membrane constituents , and even using biomembranes as sensors or actuators [44, 45] have been investigated. The complexity of the problem has forced the usage of sophisticated numerical methods, such as advanced phase field schemes, special nonlinear finite elements, and molecular dynamics simulations [17, 46-50], while analytical derivations have commonly adopted simplifying assumptions, like small deformations, axisymmetry, and so forth. The abovementioned theoretical studies have been motivated by a large body of experimental work, for example, [10-16, 21, 51-53], that demonstrated phenomena such as phase segregation, coexistence of different phases, and formation of domains in vesicles by a variety of methods, for example, fluorescence microscopy, X-ray diffraction, proton microscopy, spin resonance, and NMR imaging.
In a recent work, Parthasarathy et al.  designed an elegant experiment that breaks the two-way coupling between shape and composition and enables direct investigation of the influence of the membrane geometry on the spatial ordering of its composition. To this end, they used a quartz substrate, which was topographically patterned using photolithographic microfabrication techniques. The substrate consisted of continuously alternating high and low curvature contours with one-dimensional periodicity of 2 [micro]m. In order to decouple the main membrane from the underlying substrate a double membrane system was used: first, a supported membrane of uniform composition was deposited on the substrate. Then, the "main" membrane, a giant unilamellar vesicle (GUV), was introduced on top. Parthasarathy et al. showed that beyond a critical temperature, the spatial organization of lipid phases can be directed by gradients of membrane curvature, provided that these gradients are large enough.
In the current paper we analyze this type of experiment by means of a mathematical model combined with numerical simulations. The main goal is to reproduce the experimental observations mentioned above but also to generalize them and motivate new experiments. Accordingly, the structure of the paper is as follows: Section 2 describes the main theoretical considerations, governing equations, and nondimensional quantities that govern the spatial ordering of composition in abiomembrane subjected to imposed geometry. In Section 3, we perform a numerical study aiming at understanding the role of the nondimensional quantities that were identified in Section 2, and in particular their influence on the evolution of the spatial organization of the biomembrane composition. Focus is put on the final (steady state) spatial field of the membrane composition. The main conclusions are discussed in Section 4.
2. Theoretical Considerations
2.1. Governing Equations. Consider a biomembrane composed of two components, for example, two different lipid molecules or two different lipid phases, that lies on a smooth nonflat surface (in their experiment, Parthasarathy et al.  used a double membrane system to decouple the main membrane from the underlying substrate: a "supporting" membrane of uniform composition was deposited on the substrate in order to chemically decouple the "main" membrane from the substrate, and only then the "main" membrane was introduced on top of it) that has a geometry (shape) of a continuously alternating curvature with one-dimensional periodicity; see Figure 1(b). The free energy of the biomembrane takes the form 
[mathematical expression not reproducible], (1)
where integration is performed over the entire surface area of the membrane, A. Above, the first term is the (Helfrich) bending energy [25, 55] which depends on the mean curvature, H, the second term, f, is the specific mixing energy, and the last term describes the energetic penalty for spatial composition gradients. In addition, [phi] : A [right arrow] [0,1] describes the mole fraction of the second component, which we also refer to as local composition or concentration. A few comments are in order: (i) functional (1) does not include a stretching energy term or a Gauss-curvature bending energy term. The reason is that the biomembrane lies freely on a smooth surface; thus its stretching energy vanishes. Also, the Gauss-curvature bending energy vanishes everywhere since the imposed geometry has a 1D periodicity, which results in one of the two principal curvatures being zero. (ii) The two components of the biomembrane differ in their mechanical properties; namely, they have different bending stiffness, [k.sub.H], and spontaneous curvature, [H.sub.0]. Thus, inhomogeneity induces local spontaneous curvature and stiffness that depend on local composition, [phi]. (iii) The specific mixing energy, f, combines the aggregation enthalpy and the entropy of mixing. A simple model that is often adopted for modeling the mixing energy is the so-called "gas lattice" or "regular solution" model, which takes the form [36, 42]
[mathematical expression not reproducible]. (2)
Here, [k.sub.B] is the Boltzmann's constant, [T.sub.0] is the critical temperature defined as the temperature at which the mixing energy changes from single-well to double-well structure, [[rho].sub.0] is the lipid density (number of molecules per unit area), and [bar.T] = T/[T.sub.0] is the nondimensional temperature. Consequently, f is convex (miscible) at temperatures [bar.T] > 1 but nonconvex with double-well structure at temperatures [bar.T] < 1.
It is convenient to define g([phi]) = (1/2)[k.sub.H]([phi])[(H - [H.sub.0]([phi])).sup.2] + f([phi]); that is,
F = [[integral].sub.A](g([phi]) + 1/2[gamma][[absolute value of [nabla][phi]].sup.2])dA, (3)
The chemical potential, [mu], reflects the change in free energy due to a small change in local concentration. In our case, the concentration is a spatial field; thus the chemical potential is a spatial field as well. Accordingly, the chemical potential is defined in terms of the (functional) variation of F with respect to the concentration field ; that is,
[delta]F = [integral][mu][delta][phi]dA [??]
[mu] = [g.sub.,[phi]] - [gamma][[nabla].sup.2][phi], (4)
where [().sub.,[phi]] denotes partial differentiation with respect to [phi].
Equilibrium configurations correspond to local minima of the free energy, subjected to the relevant constraints. In our case, the system is closed so the total number of molecules of each type is preserved:
[[integral].sub.A][phi]dA = const. (5)
In order to calculate equilibrium configurations, one can use the method of Lagrange multipliers and introduce the functional
[F.sup.*] = F - [[lambda].sub.L] [integral][phi]dA. (6)
An equilibrium configuration must satisfy the condition [delta][F.sup.*] = 0; that is,
[mu] = [[lambda].sub.L] = const. (7)
This result implies that, at equilibrium, the chemical potential field must be uniform. Note that the value of this potential is not a priori known (it needs to be calculated) and depends on the specific parameters of the problem at hand, such as the overall (average) concentration and surface geometry.
In what follows, we formulate the equations that govern the evolution of the concentration field [phi]. These equations assume that [phi] evolves such that the free energy decreases at the highest rate. In numerical analysis, this is often termed the "steepest descent" method. Based on this approach, equilibrium configurations can be found by solving these equations (calculating [phi]) for long times. In reality, the details of the time-dependent evolution of the composition field [phi] may not follow exactly this strategy. Nevertheless, our main interest is in the final (equilibrium) configuration rather than in the details of the evolution, and we consider this numerical scheme as a reasonable evolution strategy. The composition, [phi], varies as long as the chemical potential field is nonuniform, or, in other words, gradients exist. Flow is from high potential to low potential; thus the concentration flux, J, takes the form of a generalized Fick's law:
J = -M[nabla][mu], (8)
where M > 0 is the mobility and vectors are identified by bold-face font. Since the overall number of molecules in the biomembrane does not change during the time scale of the experiment, changes in the concentration field are attributed only to flux. Hence, we write the following conservation law:
[[phi].sub.,t] + [nabla] x J = 0. (9)
Plugging (9) into (8) we conclude with the governing equation for the evolution of [phi]:
[[phi].sub.,t] = M[DELTA][mu] = M[DELTA]([g.sub.,[phi]] - [gamma][DELTA][phi]), (10)
[g.sub.,[phi]] = [f.sub.,[phi]] + 1/2[k.sub.H,[phi]][(H - [H.sub.0]).sup.2] + [k.sub.H](H - [H.sub.0])[H.sub.0,[phi]], (11)
[f.sub.,[phi]] = [k.sub.B][T.sub.0][[rho].sub.0][[bar.T]ln([phi]/1 - [phi]) + 2 -4[phi]]. (12)
2.2. Nondimensional Analysis. Next, we rewrite the governing equation in a nondimensional form. Besides the convenient formulation, this procedure enables us to identify of the nondimensional quantities that govern the behavior. To this end, we consider the characteristic scales of energy, length, and time.
The coefficient of the mixing energy [k.sub.B][T.sub.0][[rho].sub.0] sets the typical (specific, per unit area) energy scale. By dividing (11) with the energy scale, [k.sub.B][T.sub.0][[rho].sub.0], we obtain the nondimensional form of [g.sub.,[phi]]:
[mathematical expression not reproducible]. (13)
By introducing the typical bending stiffness, [k.sup.*.sub.H], and a typical mean curvature of the surface, [H.sup.*], the second term in (13) reads
[mathematical expression not reproducible], (14)
where [bar.H] = H/[H.sup.*] and [[bar.k].sub.H] = [k.sub.H]/[k.sup.*.sub.H] denote nondimensional curvature and bending stiffness, respectively, and [H.sup.*] is taken as the maximal curvature of the surface. The bending stiffness of each of the two membrane components (separated phases) is denoted by [k.sub.H][|.sub.[phi]=0] = [k.sub.I] and [k.sub.H][|.sub.[phi]=1] = [k.sub.II]. For specificity, and without loss of generality, we approximate the dependence of [k.sub.H] with [phi] by a linear relation . Thus, [k.sub.H,[phi]] = [k.sub.II] - [k.sub.I] and we choose [k.sup.*.sub.H] = [k.sub.I]. Similarly, we approximate the dependence of [H.sub.0] with [phi] by a linear relation, that is, [H.sub.0,[phi]] = [H.sub.II] - [H.sub.I]. Plugging these relations into (14) we have that
[m.sub.b] = 1/2[k.sub.1][([H.sup.*]).sup.2]/[k.sub.B][T.sub.0][[rho].sub.0]([[bar.k].sub.II] - [[bar.k].sub.I]). (15)
The nondimensional quantity [m.sub.b] represents the ratio between the typical bending energy and the mixing energy and also reflects the contribution of the differences between the bending stiffness of the two components.
Applying a similar procedure to the third term in (13), we find that
[k.sub.H](H - [H.sub.0])[H.sub.0,[phi]]/[k.sub.B][T.sub.0][[rho].sub.0] = [m.sub.h][[bar.k].sub.H]([bar.H] - [[bar.H].sub.0]), (16)
[m.sub.h] = [k.sub.I][([H.sup.*]).sup.2]/[k.sub.B][T.sub.0][[rho].sub.0] x ([[bar.H].sub.II] - [[bar.H].sub.I]). (17)
The nondimensional quantity [m.sub.h] represents the ratio between the bending energy and the mixing energy, similarly to [m.sub.b]; however, it reflects the contribution of the differences in spontaneous curvatures (rather than bending stiffness) between the two phases.
Nondimensional spatial coordinates (location), [[bar.x].sub.i] = [x.sub.i]/[lambda] (i = 1,2), are defined using the length of the 1D period of the surface curvature, [lambda]. The corresponding nondimensional gradient operator is [bar.[nabla]] = [lambda][nabla]. Similarly, by introducing the characteristic time scale, [tau] = [[lambda].sup.2]/M[k.sub.B][T.sub.0][[rho].sub.0], we have that [bar.t] = t/[tau] and [partial derivative]/[partial derivative][bar.t] = [tau][partial derivative]/[partial derivative]t. Using these definitions along with relations (13)-(17), we conclude with the nondimensional form of the governing equation (10):
[partial derivative][phi]/[partial derivative][bar.t] = [bar.[DELTA]]([[bar.g].sub.,[phi]] - [m.sub.1] x [bar.[DELTA]][phi]), (18)
[mathematical expression not reproducible], (19)
[m.sub.1] = 1/[[lambda].sup.2][gamma]/[k.sub.B][T.sub.0][[rho].sub.0]. (20)
Equation (18) and thus the behavior of the biomembrane are governed by three nondimensional quantities; these are [m.sub.b], [m.sub.h], and [m.sub.1] defined in (15), (17), and (20), respectively.
3. Numerical Results
In this section, we present numerical results focusing on the influence of the nondimensional quantities that were identified in the previous section on the response of a biomembrane with imposed geometry.
3.1. Numerical Scheme and Additional Considerations. In our numerical simulations we calculate the evolution in time of the composition field, [phi], in a rectangular portion of the surface. The dimensions of this rectangle are [lambda] x 2[lambda], as illustrated in Figure 1, subjected to periodic boundary conditions. This choice of domain size enables the description of the occurrences while significantly reducing computational effort (compared to simulating the entire biomembrane). The rectangular domain is discretized into [n.sub.1] points in the horizontal ([x.sub.1]) direction and [n.sub.2] points in the [x.sub.2]- direction, with equal spacing in both directions, that is, [DELTA][x.sub.1] = [DELTA][x.sub.2] = 1/([n.sub.1] - 1). Typically, a value of [n.sub.1] = 100 was used to discretize a single length unit (or [lambda] in dimensional length). We adopt a conservative second-order numerical scheme with adaptive time stepping [44, 45, 56, 57]. Figure 2 shows snapshots of the composition field at different times, where a color scale is used to describe the level of local concentration.
All simulations start from random noncorrelated values near [[phi].sub.average], which defines the overall (average) composition. These small (less than 0.01 [[phi].sub.average]) random deviations from the average composition are necessary, in some of the cases, in order to break the "symmetry" and allow for initiation of the evolution. In most simulations we have used [[phi].sub.average] = 0.3 in accordance with the experiment of Parthasarathy et al. . Still, we have also studied the influence of overall composition on the phase behavior.
In the experiments of Parthasarathy et al. , the curvature field was calculated from AFM measurements of the substrate height profile. This imposed geometry (curvature) can be well captured by the simple functional relation
[bar.H]([[bar.x].sub.1],[[bar.x].sub.2]) = [([[bar.x].sub.1] - 0.5).sup.m]/[0.5.sup.m]. (21)
Here, [[bar.x].sub.1] [member of] [0,1] is the nondimensional coordinate, where the geometry of the surface has a period of one. The coefficient m governs the "shape" of the curvature field; higher m values are associated with a larger region, around [[bar.x].sub.1] = 0.5, having flat geometry. A value of m = 2 is typical to the surface profile used in the experiments of Parthasarathy et al. .
The main purpose of the numerical investigation is to study the interplay between the imposed geometry (curvature) and composition of the biomembrane. Hence, we focus our attention on the influence of the nondimensional quantities [m.sub.b], [m.sub.h], and m on the final (steady state) composition field, with dimensional values of [lambda] = 2 [micro]m, [T.sub.0] ~ 300 K, [[rho].sub.0] ~ 5 x [10.sup.4] [[micro][m.sup.-2]], and [gamma] = [10.sup.-19] [J] [36, 54] implying that [m.sub.1] ~ [10.sup.-4]. The quantities [m.sub.b] and [m.sub.h] represent the intensity of coupling between curvature and composition through the difference between the bending stiffness and the spontaneous curvatures of the two phases, respectively. Different values of m, on the other hand, are associated with surfaces having the same maximum curvature but different topography. In the experiments of Parthasarathy et al. , the authors intentionally chose a lipid bilayer with zero spontaneous curvature. Thus, we consider first this case, namely,
[mathematical expression not reproducible]. (22)
3.2. The Effect of [m.sub.b]. Motivated by the experiments of Parthasarathy et al. , we consider first the case of zero spontaneous curvature (22), aiming at reproducing the experimental observations regarding the influence of the surface topography on the phase behavior. In particular, it was suggested that there exists a critical curvature above which the composition morphology is strongly correlated with the surface geometry. Below this critical curvature, the position of domains is rather random and does not seem to register with the geometry of the surface.
In our model, the magnitude of the surface curvature is accounted for by the nondimensional quantity [m.sub.b] (15), where larger values of [m.sub.b] represent higher surface curvature. Note that the mathematical structure of [m.sub.b] indicates that the influence of the surface curvature, through [H.sup.*2], is equivalent to that of_the difference between the stiffness of the two phases, [k.sub.I]([[bar.k].sub.II] - [[bar.k].sub.I]). Thus a higher value of [m.sub.b] reflects either higher curvature of the imposed geometry or a biomembrane with higher stiffness (or a combination of the two). In the experiments of Parthasarathy et al. , the response of biomembranes with similar composition was studied by subjecting them to surfaces with different curvatures. Thus, in these experiments [m.sub.b] was altered by changing the surface feature [H.sup.*].
Figure 3 shows the long-time (steady state) solution of the composition field obtained from numerical simulations for various values of [m.sub.b]. These results demonstrate that [m.sub.b] has a significant influence on the ordering of the two phases as a result of the reciprocity between bending stiffness and surface topography. In particular, energy considerations favor configurations where the stiffer phase is located at regions of lower curvature. From these images, we can learn more on the effect of [m.sub.b] on the biomembrane behavior. As [m.sub.b] decreases, the domains become less organized and show lower correspondence with the surface topography. For example, with m = 2 and [m.sub.b] = [10.sup.-1], the segregated domains of the stiffer (red) phase are oval and perfectly centered at [[bar.x].sub.1] = 0.5, while for [m.sub.b] = [10.sup.-4] the location of the domains is random with no particular preference, and their shape is almost perfectly round due to surface tension which becomes more dominant for small values of [m.sub.b]. These results suggest that the surface geometry affects the phase behavior only when the surface curvature is high enough. In accordance with observations of Parthasarathy et al. , values of [m.sub.b] [equivalent] [10.sup.-4] or smaller correspond to negligible influence of the imposed geometry, while increasingly higher values of [m.sub.b] result in an increasing effect.
Recall that [m.sub.b] also reflects the difference between the bending stiffness of the phases. Hence, we generalize the conclusion of Parthasarathy et al. , which was specific to the lipid phases used in their experiments. By understanding the dual role of [m.sub.b], we conclude that higher stiffness ratios between the two phases decrease the magnitude of the minimal surface features required to couple between the surface topography and the biomembrane composition. This conclusion is somewhat intuitive, but now the mathematical structure of [m.sub.b] describes it quantitatively.
The parameter m governs the width of the flat section of the surface topography. Its influence on the phase behavior is exemplified by comparing the results of Figures 3(a) and 3(b). Due to the energy-related reasoning discussed above, surfaces with smaller values of m constrict the domains of the stiffer phase to a narrower region in the middle of the surface (where the curvature is small). As a result, the circular shape of the domains becomes more oval as m decreases (or [m.sub.b] increases), until the point where constriction is so tight that all domains merge and form a single strip; see for example, Figure 3(a) with m = 2 and [m.sub.b] = 1. Note that such extreme morphology has not been observed in the experiments. The reason is that this type of behavior requires a combination of small m and very large [m.sub.b]. Specifically, the maximum values reached in the experiments of Parthasarathy et al.  were [m.sub.b] [equivalent] 0.01 with m = 2.
Overall Concentration. The stripe morphology can take place for lower values of [m.sub.b] by increasing the overall (average) concentration [[phi].sub.average]. This is exemplified in Figure 4, which shows the effect of overall composition on the ordering of the two phases. For example, with m = 2 and [m.sub.b] = 0.1 the stripe morphology does not appear with overall concentration of [[phi].sub.average] = 0.3 but appears with [[phi].sub.average] = 0.5 or higher. Also, with m = 6 and [[phi].sub.average] = 0.3, the stripe morphology does not appear even for [m.sub.b] = 1; see Figure 3(b) but does appear at [m.sub.b] = 0.1 with [[phi].sub.average] = 0.7; see Figure 4(b). The reasoning for this phenomenon is that the (energetic) avoidance of the stiffer phase from high curvature constricts it to a "stripe" of relatively small curvature in the middle. If the overall concentration of the stiffer phase is low, it forms small round domains that remain separated since they are far enough apart within this stripe. On the other hand, if the overall concentration of the stiffer phase is higher, the stiffer phase almost fills the stripe, which leads to the formation of oval domains, for example, Figure 4(b) with [[phi].sub.average] = 0.5. At even higher concentrations, all domains merge into a single strap in order to minimize line (surface) tension energy, as shown in Figure 4(b) with [[phi].sub.average] = 0.7.
Similarly to Figure 3, the difference between Figures 4(a) and 4(b) stems from the competition between the bending energy and the line tension energy. When the surface has a wider flat section (and the overall concentration is small), the bending energy is less significant within this flat section and the system reduces the line tension energy by forming circular domains that have shorter phase boundaries. However, when the surface has a relatively narrow flat section, the system must restrict the stiffer phase to a narrow strap in order to reduce the bending energy, which becomes significant outside the strap, on the expense of the energy penalty for longer phase boundaries. This leads to the formation of oval domains, and, in cases where the overall concentration is high enough, it leads to the stripe morphology.
Behavior above the Critical Temperature. The simulations presented above, just like the experiments of Parthasarathy et al. , were performed at a temperature lower than the miscibility temperature [T.sub.0]. Hence, the double-well structure of the interaction energy, f([phi]), drives the system towards phase separation and formation of domains. In turn, the topography of the surface (the imposed geometry) affects the spatial ordering and shape of these domains. On the other hand, above the critical temperature, the interaction energy favors miscibility, that is, uniform composition. Such behavior was indeed observed in the experiments as well as in our numerical simulations. Nevertheless, our simulations also show that in cases where the magnitude of the surface features (curvature) is high enough, the imposed geometry can lead to the formation of a spatially nonuniform composition field and to geometry-induced phase separation above the critical temperature; see Figure 5. One must note that this type of nonstandard behavior requires high values of [m.sub.b]. Also, unlike standard phase separation that often exhibits round domains, here the spatial distribution of the composition field is dominated almost completely by the surface features which gives rise only to a stripes-like ordering. Unfortunately, the surface features used in the experiment of Parthasarathy et al.  were much smaller; thus validation of this phenomenon still awaits experimental confirmation.
3.3. The Role of Spontaneous Curvature. A comparison between the definitions of [m.sub.b] (15) and [m.sub.h] (17) suggests that the difference in spontaneous curvatures, [[bar.H].sub.I] and [[bar.H].sub.II], and the difference in bending stiffness, [[bar.k].sub.I] and [[bar.k].sub.II], of the two phases have a similar influence on the biomembrane composition. In particular, (13) indicates that bending energy drives the stiffer phase towards smaller surface curvatures (second term in (13)) and also the phase with higher spontaneous curvature towards locations of higher surface curvature (third term in (13)). Hence, the differences in the bending stiffness and in spontaneous curvature between the two phases make each of the two phases favor (energetically) different locations on the surface. The magnitude of this effect is largely associated with the magnitude of the nondimensional quantities [m.sub.b] and [m.sub.b]. These two quantities have a similar mathematical structure which expresses the relative importance of the bending energy compared to the interaction energy multiplied by the difference in bending stiffness or spontaneous curvatures of the two phases, respectively. Hence, the role of [m.sub.h] and its influence on the biomembrane response when subjected to imposed geometry seems to be qualitatively similar to that of [m.sub.b], which has been studied in the previous section. Nevertheless, we note a fundamental difference between the second and third terms in (13). That is, while [m.sub.b] multiplies [([bar.H] - [[bar.H].sub.0].sup.2], [m.sub.h] multiplies ([bar.H] - [[bar.H].sub.0]). The importance of this difference is twofold: (i) changes in the geometry of the biomembrane, for example, by subjecting it to a surface with higher curvature, affect more significantly the second term compared to the third term in (13). In particular, in the case where [m.sub.b] and [m.sub.h] are comparable, the third term in (13) is more sensitive to the imposed geometry. Thus, roughly speaking, the role of the bending stiffness in ordering the biomembrane composition is more significant than that of the spontaneous curvature. (ii) The sign of the third term in (13) depends on the sign of H - [H.sub.0], while that of the second term does not.
Following the discussion above, we focus our attention below on demonstrating the consequences of the fact that the sign of the third term in (13) depends on the sign of H - [H.sub.0]. To this end, we consider the simple case where [k.sub.II] = [k.sub.I] and [H.sub.II] = -[H.sub.I] are together with the typical dimensional values of [k.sub.I] ~ [10.sup.-19] J, [absolute value of [H.sup.*]] ~ 1[micro][m.sup.-1], [T.sub.0] ~ 300 K, and [[rho].sub.0] ~ 5 x [10.sup.4] [[micro][m.sup.-2]] [36, 54], which imply that [m.sub.b] = 0, [H.sub.0] = (1 - 2[phi])[H.sub.I], and [m.sub.h] ~ -[10.sup.-3][[bar.H].sub.I]. Figure 6 shows the difference in the final (steady state) spatial distribution of the biomembrane composition, [phi], when imposed on convex and concave surfaces and for different values of [[bar.H].sub.I] and of the average concentration. As expected, higher values of [H.sub.I] (which also mean higher difference between the spontaneous curvatures of the two phases) increase the correlation between the observed composition patterns and the surface topography. High values of [H.sub.I] lead to the formation of oval domains and when [H.sub.I] is high enough to the stripe morphology. These features are qualitatively similar to the behavior observed in the previous section where the effect of the difference in bending stiffness was studied. A fundamental difference, however, is in the distribution of the biomembrane composition, [phi], when imposed to convex or concave surfaces. This is illustrated by comparing results in Figure 6 that differ only in the sign of [H.sup.*] (top row compared to the bottom row). A few comments are in order: (i) the different sign of the surface curvature inverts the spatial location favored by each phase. Importantly this is not the case when the two phases differ only in their bending stiffness; (ii) the inversion of preferred locations leads to a different constriction exhibited by the phases; this may result in significant changes in the observed morphology; for example, oval compared to circular domains or stripe morphology compared to oval domains; (iii) the specific example studied here, with [H.sub.II] = -[H.sub.I], is antisymmetric with respect to (sign of) the surface curvature. Thus, cases involving opposite sign of the surface curvature combined with "opposite" average concentrations are expected to yield opposite spatial organization, that is, similar morphology with inverted phase locations. This feature is indeed observed when comparing results in the upper-left and bottom-right "windows" of Figure 6, as well as upper-right and bottom-left "windows" in the same figure.
4. Summary and Conclusions
We formulated a simple mathematical model to study the spatial ordering of composition in a two-component biomembrane that is subjected to prescribed (imposed) geometry. The mathematical model does not account for possible anisotropy of the membrane constituents or for possible interaction between lipids in the two leaflets of the bilayer membrane. In addition, the numerical scheme, which follows the steepest descent method, leads to metastable equilibrium configurations associated with local minima of the free energy. We note, however, that applying the numerical scheme to slightly different initial conditions or to a larger domain of solution did not change the essence (topology) of the solution. Based on this model, we identified key nondimensional quantities that govern the biomembrane response and performed numerical simulations to quantitatively study their influence. Our numerical results show that the geometry-driven ordering of the biomembrane composition is largely governed by the difference between the nondimensional bending stiffness and spontaneous curvatures of each phase, while the magnitude of this phenomenon is proportional to the ratio between the bending energy and the (chemical) interaction energy of the phases. Roughly speaking, energy considerations favor configurations in which the phase that is stiffer and has smaller spontaneous curvature is located at regions having smaller curvature. The numerical simulations reproduced the experimental observation of Parthasarathy et al. , who found that above a critical surface curvature the composition morphology is strongly correlated with the surface geometry, while below this threshold, the position of domains is rather random and does not register with the geometry of the surface. Careful investigation of our model equations enable us to generalize this experimental observation beyond the specific lipid phases used in that experiment.
An important advantage of a mathematical model is that it enables studying the behavior at various settings with minimal effort and resources, before entering the lab. The agreement of our model results with experimental observations strengthens our confidence in the model and numerical scheme and opens the door to examining new and different conditions than those used in the original experiments. For example, we have demonstrated that, if the surface geometry is properly designed, phase separation can occur above the critical temperature. Such curvature-induced phase separation above the critical temperature awaits experimental examination. Also, we propose a systematic procedure to determine which mechanism, the difference in bending stiffness or difference in spontaneous curvatures of the two phases, dominates the coupling between shape and composition. The procedure is based on the observation that the mechanism associated with the difference in bending stiffness depends on the magnitude of the surface curvature but indifferent to the sign (direction) of the curvature. On the contrary, the mechanism related to the spontaneous curvatures strongly depends on both magnitude and sign (direction) of the surface curvature. The consequences of these differences have been demonstrated by a set of simulations.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
 D. E. Ingber, "Cellular mechanotransduction: putting all the pieces together again," The FASEB Journal, vol. 20, no. 7, pp. 811-827, 2006.
 E. Perozo and D. C. Rees, "Structure and mechanism in prokaryotic mechanosensitive channels," Current Opinion in Structural Biology, vol. 13, no. 4, pp. 432-442, 2003.
 D. Stamenovic and N. Wang, "Invited review: engineering approaches to cytoskeletal mechanics," Journal of Applied Physiology, vol. 89, no. 5, pp. 2085-2090, 2000.
 V. Vogel and M. Sheetz, "Local force and geometry sensing regulate cell functions," Nature Reviews Molecular Cell Biology, vol. 7, no. 4, pp. 265-275, 2006.
 C. Zhu, G. Bao, and N. Wang, "Cell mechanics: mechanical response, cell adhesion, and molecular deformation," Annual Review of Biomedical Engineering, vol. 2, no. 1, pp. 189-226, 2000.
 D. M. Engelman, "Membranes are more mosaic than fluid," Nature, vol. 438, no. 7068, pp. 578-580, 2005.
 B. Alberts, Essential Cell Biology: An Introduction to the Molecular Biology of the Cell, vol. 1, Garland, NewYork, NY, USA, 1997.
 S. J. Singer and G. L. Nicolson, "The fluid mosaic model of the structure of cell membranes," Science, vol. 175, no. 4023, pp. 720-731, 1972.
 D. H. Boal, Mechanics of the Cell, Cambridge University Press, New York, NY, USA, 2002.
 T. Baumgart, S. Das, W. W. Webb, and J. T. Jenkins, "Membrane elasticity in giant vesicles with fluid phase coexistence," Biophysical Journal, vol. 89, no. 2, pp. 1067-1080, 2005.
 T. Baumgart, S. T. Hess, and W. W. Webb, "Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension," Nature, vol. 425, no. 6960, pp. 821-824, 2003.
 N. Kahya, D. Scherfeld, K. Bacia, B. Poolman, and P. Schwille, "Probing lipid mobility of raft-exhibiting model membranes by fluorescence correlation spectroscopy," Journal of Biological Chemistry, vol. 278, no. 30, pp. 28109-28115, 2003.
 S. L. Veatch and S. L. Keller, "Separation of liquid phases in giant vesicles of ternary mixtures of phospholipids and cholesterol," Biophysical Journal, vol. 85, no. 5, pp. 3074-3083, 2003.
 S. L. Veatch and S. L. Keller, "Seeing spots: complex phase behavior in simple membranes," Biochimica et Biophysica Acta--Molecular Cell Research, vol. 1746, no. 3, pp. 172-185, 2005.
 M. R. Vist and J. H. Davis, "Phase-equilibria of cholesterol dipalmitoylphosphatidylcholine mixtures--H-2 nuclear magnetic-resonance and differential scanning calorimetry," Biochemistry, vol. 29, no. 2, pp. 451-464, 1990.
 T. Auth and G. Gompper, "Budding and vesiculation induced by conical membrane inclusions," Physical Review E, vol. 80, no. 3, Article ID 031901, 2009.
 L. Ma and W. S. Klug, "Viscous regularization and r-adaptive remeshing for finite element analysis of lipid membrane mechanics," Journal of Computational Physics, vol. 227, no. 11, pp. 5816-5835, 2008.
 U. Seifert, "Curvature-induced lateral phase segregation in two-component vesicles," Physical Review Letters, vol. 70, no. 9, pp. 1335-1338, 1993.
 P. Sens, L. Johannes, and P. Bassereau, "Biophysical approaches to protein-induced membrane deformations in trafficking," Current Opinion in Cell Biology, vol. 20, no. 4, pp. 476-482, 2008.
 V. B. Shenoy and L. B. Freund, "Growth and shape stability of a biological membrane adhesion complex in the diffusion-mediated regime," Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 9, pp. 3213-3218, 2005.
 B. Sorre, A. Callan-Jones, J. B. Manneville et al., "Curvature-driven lipid sorting needs proximity to a demixing point and is aided by proteins," Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 14, pp. 5622-5626, 2009.
 H. Sprong, P. van der Sluijs, and G. van Meer, "How proteins move lipids and lipids move proteins," Nature Reviews Molecular Cell Biology, vol. 2, no. 7, pp. 504-513, 2001.
 H. T. McMahon and J. L. Gallop, "Membrane curvature and mechanisms of dynamic cell membrane remodelling," Nature, vol. 438, no. 7068, pp. 590-596, 2005.
 E. A. Evans, "Bending resistance and chemically-induced moments in membrane bilayers," Biophysical Journal, vol. 14, no. 12, pp. 923-931, 1974.
 W. Helfrich, "Elastic properties of lipid bilayers--theory and possible experiments," Zeitschrift fur Naturforschung C, vol. 28, no. 11, pp. 693-703, 1973.
 U. Seifert, "Configurations of fluid membranes and vesicles," Advances in Physics, vol. 46, no. 1, pp. 13-137, 1997.
 U. Seifert, K. Berndl, and R. Lipowsky, "Shape transformations of vesicles: phase diagram for spontaneous-curvature and bilayer-coupling models," Physical Review A, vol. 44, no. 2, pp. 1182-1202, 1991.
 A. Agrawal and D. J. Steigmann, "Coexistent fluid-phase equilibria in biomembranes with bending elasticity," Journal of Elasticity, vol. 93, no. 1, pp. 63-80, 2008.
 E. L. Elson, E. Fried, J. E. Dolbow, and G. M. Genin, "Phase separation in biological membranes: integration of theory and experiment," in Annual Review of Biophysics, vol. 39, pp. 207-226, Annual Reviews, Palo Alto, Calif, USA, 2010.
 N. Walani, J. Torres, and A. Agrawal, "Anisotropic spontaneous curvatures in lipid membranes," Physical Review E, vol. 89, no. 6, Article ID 062715, 2014.
 J. L. Harden, F. C. MacKintosh, and P. D. Olmsted, "Budding and domain shape transformations in mixed lipid films and bilayer membranes," Physical Review E, vol. 72, no. 1, Article ID 011903, 2005.
 F. Julicher and R. Lipowsky, "Domain-induced budding of vesicles," Physical Review Letters, vol. 70, no. 19, pp. 2964-2967, 1993.
 H. Giang, R. Shlomovitz, and M. Schick, "Microemulsions, modulated phases and macroscopic phase separation: a unified picture of rafts," Essays In Biochemistry, vol. 57, pp. 21-32, 2015.
 I. Fonseca, G. Hayrapetyan, G. Leoni, and B. Zwicknagl, "Domain formation in membranes near the onset of instability," Journal of Nonlinear Science, vol. 26, no. 5, pp. 1191-1225, 2016.
 W. T. Gozdz, N. Bobrovska, and A. Ciach, "Separation of components in lipid membranes induced by shape transformation," The Journal of Chemical Physics, vol. 137, no. 1, Article ID 015101, 2012.
 S. Givli, H. Giang, and K. Bhattacharya, "Stability of multi-component biological membranes," SIAM Journal on Applied Mathematics, vol. 72, no. 2, pp. 489-511, 2012.
 D. Ni, H. J. Shi, Y. J. Yin, and L. S. Niu, "Stability of biphasic vesicles with membrane embedded proteins," Journal of Biomechanics, vol. 40, no. 7, pp. 1512-1517, 2007.
 S. A. Safran, P. A. Pincus, D. Andelman, and F. C. Mackintosh, "Stability and phase-behavior of mixed surfactant vesicles," Physical Review A, vol. 43, no. 2, pp. 1071-1078, 1991.
 N. Walani and A. Agrawal, "Stability of lipid membranes," Mathematics and Mechanics of Solids, 2015.
 N. Gov, A. G. Zilman, and S. Safran, "Cytoskeleton confinement and tension of red blood cell membranes," Physical Review Letters, vol. 90, no. 22, Article ID 228101, 2003.
 Y. Adler and S. Givli, "Closing the loop: lamellipodia dynamics from the perspective of front propagation," Physical Review E, vol. 88, no. 4, Article ID 042708, 2013.
 A. Veksler and N. S. Gov, "Phase transitions of the coupled membrane-cytoskeleton modify cellular shape," Biophysical Journal, vol. 93, no. 11, pp. 3798-3810, 2007.
 N. Bobrovska, W. Gozdz, V. Kralj-Iglic, and A. Iglic, "On the role of anisotropy of membrane components in formation and stabilization of tubular structures in multicomponent membranes," PLoS ONE, vol. 8, no. 9, Article ID e73941, 2013.
 L. Atia and S. Givli, "Biological membranes from the perspective of smart materials--a theoretical study," International Journal of Solids and Structures, vol. 49, no. 18, pp. 2617-2624, 2012.
 L. Atia and S. Givli, "A theoretical study of biological membrane response to temperature gradients at the single-cell level," Journal of the Royal Society Interface, vol. 11, no. 95, 2014.
 I. R. Cooke and M. Deserno, "Coupling between lipid shape and membrane curvature," Biophysical Journal, vol. 91, no. 2, pp. 487-495, 2006.
 Q. Du, C. Liu, and X. Q. Wang, "A phase field approach in the numerical study of the elastic bending energy for vesicle membranes," Journal of Computational Physics, vol. 198, no. 2, pp. 450-468, 2004.
 F. Feng and W. S. Klug, "Finite element modeling of lipid bilayer membranes," Journal of Computational Physics, vol. 220, no. 1, pp. 394-408, 2006.
 H. Agrawal, L. P. Liu, and P. Sharma, "Revisiting the curvature-mediated interactions between proteins in biological membranes," Soft Matter, vol. 12, no. 43, pp. 8907-8918, 2016.
 J. S. Lowengrub, A. Ratz, and A. Voigt, "Phase-field modeling of the dynamics of multicomponent vesicles: spinodal decomposition, coarsening, budding, and fission," Physical Review E, vol. 79, no. 3, article 13, 2009.
 J. C. S. Ho, P. Rangamani, B. Liedberg, and A. N. Parikh, "Mixing water, transducing energy, and shaping membranes: autonomously self-regulating giant vesicles," Langmuir, vol. 32, no. 9, pp. 2151-2163, 2016.
 A. Roux, D. Cuvelier, P. Nassoy, J. Prost, P. Bassereau, and B. Goud, "Role of curvature and phase transition in lipid sorting and fission of membrane tubules," The EMBO Journal, vol. 24, no. 8, pp. 1537-1545, 2005.
 A. V. Samsonov, I. Mihalyov, and F. S. Cohen, "Characterization of cholesterol-sphingomyelin domains and their dynamics in bilayer membranes," Biophysical Journal, vol. 81, no. 3, pp. 1486-1500, 2001.
 R. Parthasarathy, C.-H. Yu, and J. T. Groves, "Curvature-modulated phase separation in lipid bilayer membranes," Langmuir, vol. 22, no. 11, pp. 5095-5099, 2006.
 O.-Y. Zhong-can and W. Helfrich, "Bending energy of vesicle membranes: general expressions for the first, second, and third variation of the shape energy and applications to spheres and cylinders," Physical Review A, vol. 39, no. 10, pp. 5280-5288, 1989.
 S. Givli, "Towards multi-scale modeling of muscle fibers with sarcomere non-uniformities," Journal of Theoretical Biology, vol. 264, no. 3, pp. 882-892, 2010.
 S. Givli and K. Bhattacharya, "A coarse-grained model of the myofibril: overall dynamics and the evolution of sarcomere non-uniformities," Journal of the Mechanics and Physics of Solids, vol. 57, no. 2, pp. 221-243, 2009.
Shimrit Katz and Sefi Givli
Faculty of Mechanical Engineering, Technion-Israel Institute of Technology,
32000 Haifa, Israel
Correspondence should be addressed to Sefi Givli; email@example.com
Received 30 January 2017; Accepted 22 March 2017; Published 4 April 2017
Academic Editor: John Mitchell
Caption: FIGURE 1: Composition field, [phi], is calculated in a rectangular section of dimension [lambda] x 2[lambda], typically discretized into 100 x 200 points, and subjected to periodic boundary conditions. Here, [lambda] is the 1D period of the surface curvature in the experiment. The value of [phi] [member of] [0,1] at each point is shown by different colors (see color bar). The image on the right shows a schematic illustration of the topographically patterned quartz substrate (gray), supported membrane (green), and upper (red) membrane in the experiment of Parthasarathy et al.  (reproduced with permission from ).
Caption: FIGURE 2: Snapshots of typical simulations. (a) m = 2, [m.sub.b] = 0.1, [m.sub.h] = 0, and [[phi].sub.average] = 0.3 at times [bar.t] = (1,10,30, 50,60,80) x [10.sup.4] (from left to right). The simulation was stopped at [bar.t] = 5 x [10.sup.6] with no visible change compared to the last snapshot shown here. (b) Same as (a), but with [[phi].sub.average] = 0.5, which is in the spinodal region of f([phi]). Videos of these simulations are provided in Supplementary Material available online at https://doi.org/10.1155/2017/7275131.
Caption: FIGURE 3: Final (steady state) composition field for different values of [m.sub.b] and [H.sub.0]([phi]) = 0, for (a) m = 2 and (b) m = 6. In both (a) and (b), the leftmost plots show the imposed curvature, while the other four plots illustrate, from left to right, the composition field [phi]([[bar.x].sub.1],[[bar.x].sub.2]) for [m.sub.b] = 1, [10.sup.-1], [10.sup.-3], [10.sup.-4]. [[phi].sub.average] = 0.3.
Caption: FIGURE 4: The effect of overall (average) composition, [[phi].sub.average], on the final (steady state) composition field. (a) m = 2 and (b) m = 6. In both (a) and (b), the leftmost plots show the imposed curvature, while the other four plots illustrate, from left to right, the composition field [phi]([[bar.x].sub.1],[[bar.x].sub.2]) for [[phi].sub.average] = 0.3,0.5, 0.7. [m.sub.b] = 0.1 and [H.sub.0]([phi]) = 0.
Caption: FIGURE 5: Behavior above critical temperature. (a) [m.sub.b] = 0.01, (b) [m.sub.b] = 0.1, (c) [m.sub.b] = 1, (d) [m.sub.b] = 10. m = 2, and [[phi].sub.average] = 0.3.
Caption: FIGURE 6: The effect of spontaneous curvature for different values of [[phi].sub.average] and sign of the surface curvature. Each plot shows the final (steady state) spatial distribution of the composition field In each of the four "windows" the left and right plots, (i) and (ii), correspond to [[bar.H].sub.I] = 0.1 and [[bar.H].sub.I] = 1, respectively. m = 2, [k.sub.I] = [k.sub.II], and [H.sub.II] = -[H.sub.I].
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Author:||Katz, Shimrit; Givli, Sefi|
|Publication:||Computational and Mathematical Methods in Medicine|
|Date:||Jan 1, 2017|
|Previous Article:||Mathematical Modeling and Control of Infectious Diseases.|
|Next Article:||Threshold Dynamics in Stochastic SIRS Epidemic Models with Nonlinear Incidence and Vaccination.|