You are an expert at summarizing long articles. Proceed to summarize the following text: 
this work is motivated by the formation of lipid rafts in lipid membranes . lipid bilayer membranes are of utmost importance for the survival of cells . they separate the interior of cells from the extracellular environment and compartmentalize subcellular organelles so suitable micro - environment can be maintained in the enclosed domains for various vital biochemical and biophysical reactions . they are material basis for morphological changes such as budding , tubulation , fission and fussion that occur during cell division , biological reproduction , and intracellular membrane trafficking . they also provide a physical platform to store and transduce energy as electrochemical gradients , to segregate or disperse particular membrane proteins , and to act as messengers in signal transduction and molecular recognition processes @xcite . while most of these functionalities depends on the fluidity of the lipids and thereby the free diffusion of lipids and proteins in the bilayer , accumulated evidences show that lipids and proteins on bilayer membranes segregates into discrete domains of distinct composition and various sizes @xcite . the domain boundaries can appear as the barriers of free lateral diffusion of lipids and proteins , as the measured lateral diffusion coefficients of lipids and proteins _ in vivo _ are less than the measured coefficients in artificial pure bilayer by more than one order of magnitude @xcite . inside the domains and on the domain boundaries particular proteins may aggregate , to cause various membrane curvature as the consequence of the modification of local membrane composition @xcite or to complete specific signal transduction @xcite . some of these domains are transient , with a duration ranging from seconds to minutes , some can persist for the entire life of the cell , and the domains themselves can diffuse on the membrane surface as well @xcite . the composition , location , size , configuration , duration of these domains and the dynamics of these characteristics are of functional and structural significance to the associated biological processes . efforts integrating direct microscopic measurements , biophysical modeling , and computational simulations have been invested to elucidate the underlying physics of the dynamics of these lipid domains and predict their biological consequences @xcite . before introducing our approach based on the geodesic curvature energy of the lipid domain boundaries we first review four most representative theoretical studies on the dynamics of lipid domains . lipid domains may appear as a result of lipid phase separation caused by distinct spontaneous curvatures . when bilayer membranes have multiple lipid species of distinct spontaneous curvatures , individual lipid species may be localized to regions where the local mean curvatures best approximate the corresponding spontaneous curvatures of the residing lipid species @xcite . wang and du formalized this reasoning by summing up the classical canham - evans - helfrich energy @xcite for each individual lipid species and the line tension energy to generate a multi - component lipid membrane model @xcite . by representing the membrane bending energy using the phase field formulation , they have obtained rich patterns of membrane morphology and the generation of lipid membrane domains of different mean curvatures , where lipid species of the approximate spontaneous curvatures are concentrated . this model was also extended to simulate the open membrane thanks to the line tension energy , and the closing of membrane pores was simulated corresponding to the vanishing linear tension energy . these permanent domains have sizes that are determined by local mean curvatures of the membrane necks or bumps . these sizes in general do not match the measured sizes of mobile lipid rafts @xcite . the classical phase separation model based on the ginzburg - landau ( gl ) free energy could also be directly applied on a membrane surface to generate surface phase separation , and the results can be related to the lipid domains . a surface cahn - hilliard equation can be derived for the gradient flow of the gl free energy , and the numerical simulations will produce large separated domains as a result of the coarsening dynamics @xcite . in order to generate small domains at spatial and temporal scales comparable to experimental results , camley and brown couples the gl free energy for quasi two - dimensional binary lipids mixtures to the random hydrodynamics and thermal fluctuations @xcite . the random in - plane velocity field of the membrane is given by saffman - delbruck hydrodynamic model @xcite . this velocity field is added to the cahn - hilliard equation for the gradient flow of the gl free energy to produce an advection - diffusion equation , above which a gaussian white noise is added , modeling the thermal fluctuation as a random source to the order parameter . complete phase separation shall occur as the result of a sequence of coarsening dynamics when the gl free energy is minimized , while the domain boundaries flicker as a result of random hydrodynamic and thermal perturbations , with a flickering magnitude depending on the competing between the random perturbation and the persisting linear tension . under high line tension small domains will merge to form large separated domain , but small domains under a critical size could remain separated for a long time during the course of coarsening if the line tension is not large enough to suppress the random perturbation , giving rise to lasting microdomains . various dynamical scaling rates were summarized to related the microdomain size and the time when the domain size is far way from the saffmann - delbruck length @xmath0 determined by the relative viscosity of the lipid membrane with respect to the surrounding fluid field . this approach has been recently extended to model multicomponent membranes with embedded proteins @xcite . it is also possible to simulate lipid microdomains in the deterministic setting . arguing that lipid rafts are microdomains of lipids compactly organized around embedded protein receptors , witkowski , backofen and voigt proposed to supplement the classical ginzburg - landau free energy with gaussian potentials localized at specified positions where the membrane proteins are supposed to be embedded @xcite . by specifying the center and modulation of these external potentials they were able to produce lipid microdomains of arbitrary size at arbitrary position . coarsening dynamics were reproduced by solving the cahn - hilliard equation for the gradient flow of the total energy , and a scaling law was deduced for the growth of the microdomains . in contrast to the above approaches , one is not able to drive the lateral diffusion coefficient of the microdomains as their positions are specified in the construction of the free energy . in general this model lacks a biophysical interpretation of the external potential and the related parameters that could justify the striking generation of lipid microdomains in the absence of line tension . in addition to the above continuum approaches , particle - based discrete methods have also been developed to simulate the lipid microdomains . molecular dynamics simulations , fully atomic or coarse grained , as reviewed in @xcite , have been able to generate lipid membrane domains that could interpret some experiments on complex model membranes . two discrete methods that gives particular valuable insight into the structure and dynamics of lipid microdomains are dynamical triangulation monte carlo ( dtmc ) @xcite and dispersive particle dynamics ( dpd ) @xcite . dtmc neglects the solvent hydrodynamics and approximate a bilayer membrane as a randomly triangulated sheet . each vertex is described by a three - dimensional position vector , and all vertices are connected by flexible tethers , which flip during the course of dynamical triangulation to simulate phase segregation . dpd adopts a coarse - grained representation of amphiphilic lipids as connected head ( h ) and tail ( c@xmath1 ) beads , with an variable number of tail beads . the geometric and molecular mechanics representations of lipids in dpd differ from the coarse - grained molecular dynamics simulations ( those based on the martini force field @xcite , for example ) in that all beads are soft with interaction defined by effective forces that reproduce the hydrodynamic behavior of fluid bilayer membrane rather than the classical intermolecular interactions . dpd allows asymmetric lipid composition in the two leaflets , where different lateral sizes of the lipid domains can be simulated . while dpd can simulate the membrane properties to length and time scales that are unattainable by md and coarse - grained md simulations , there still exist gaps for these discrete methods to model large systems containing millions of amphiphiles over biologically relevant time scales , for the handling of which the continuum models have intrinsic advantages .    here we work to develop a lipid micro - organization model by recognizing the intrinsic geodesic curvatures of the boundaries between lipids of different species ( saturated and unsaturated , ordered and disordered lipids , too ) . experimental observations and molecular dynamics simulations show that lipid microdomains have relatively stable size ranging from below 10 nm to larger than 200 nm , depending on the lipid composition and temperature @xcite . modulo the perturbations to the boundaries of these microdomains caused by random forces due to the thermal fluctuations , hydrodynamic interaction etc , the size of an individual microdomain can be characterized by the curvature of boundary , which is the geodesic curvature since the microdomain is a patch on the three - dimensional surface . consequently , one can identify an _ intrinsic geodesic curvature _ for a binary mixture of lipid species , similar to the identification of the spontaneous curvature of an individual lipid species . a geodesic curvature energy can be defined for a given binary mixture of lipid species on a membrane surface , in a way similar to the definition of canham - evans - helfrich curvature energy . it is expected that the minimization of this geodesic curvature energy will generate the optimal interfaces of the binary mixture of lipids . the minimization is conducted on the 2-manifold that represents the membrane surface , and these optimal interfaces define the boundaries of the lipid microdomains . since the optimal interfaces are unknown _ a priori _ , and there will be topological changes as the initial interfaces evolve during the course of curvature energy minimization , we adopt an eulerian formulation of the geodesic curvature energy by relating the geodesic curvature to a phase field function on the membrane surface . we shall show that the phase field representation of the geodesic curvature energy can be regarded as a generalization of the phase field modeling of the elastic bending curvature energy of bilayer membrane deformation @xcite . the equivalence between the geodesic curvature energy and the phase field representation can be justified when the phase field function is given by a properly scaled hyperbolic tangent function of the signed geodesic distance to the boundaries of lipid microdomains . the minimizer of the geodesic curvature energy in the eulerian representation is obtained by evolving the gradient flow of the phase field function . the forth - order evolutionary partial differential equation for the gradient flow are numerically solved by using a @xmath2 interior penalty discontinuous galerkin method on the triangulated membrane surface . computational simulations using the geodesic curvature energy model on different surfaces demonstrate that our model can produce microdomains with specified intrinsic geodesic curvatures . the rest of the article is organized as follows . in section [ sect : model ] a brief review of the structure of lipid bilayer membrane is followed by the introduction of geodesic curvature of the interfaces between saturated and unsaturated phases of lipids . geodesic curvature energy is then defined in lagrangian and eulerian formulations , respectively . technical details of the justification of some critical properties of geodesic distance and geodesic curvature are given in appendix [ sect : append ] . the gradient flow of the geodesic curvature energy is derived using the energetical variational principle in section [ sect : numeric ] . we split the linear and nonlinear components of the governing equation for the gradient flow so that the linear components can be solved using an implicit time marching method and the convergence at each time step is obtained by efficient iterations over the nonlinear components . the proposed mathematical model and the related numerical methods are applied in section [ sect : simulations ] on a set of static surfaces with varying curvatures to examine their ability in simulating the generation of microdomains of specified geodesic curvatures . we summarize this work and outline its perspectives in the final section [ sect : summary ] . we consider a bilayer membrane with binary mixture of lipids . these can be two species of lipids , or lipids in saturated and unsaturated states , or lipids in ordered or disordered states . these two species of lipids have their own regular molecular geometry at the equilibrium state for a given temperature , modulo random thermal and hydrodynamic fluctuations . by regular we mean the uniform lateral spacing of lipid molecules , the ordered match of hydrophobic tails of two leaflets , among other features . this regular molecular geometry dictates the intermolecular interactions among lipids of same species . at the interface between domains of distinct lipid species , the regular molecular geometry of either species has to be relaxed in a way such that the intermolecular interactions in the transitional region near the interface will fit the different molecular geometry of the other species . this relaxation shall create curved interface between two species in a manner similar to the creation of the surface tension at a fluid surface . [ fig : surface_pattern_model ] illustrates the mismatch of lipid molecular structure at the interface between saturated ( s ) and unsaturated ( u ) lipids . because the lipid domains are to be modelled as patches on membrane surface , the domain boundaies will be curves on a two - dimensional surface embedded in @xmath3 . notice that if the interface is a geodesic , then it is a locally straight line that does not curve to either domain it separates . how far an interface is from being a geodesic is characterized by the geodesic curvature , which by definition is the curvature of the interface projected onto the tangent plane of the surface . the geodesic curvature will be an intrinsic property of _ a binary lipid mixture _ , similar to the spontaneous mean curvature being an intrinsic property of _ an individual lipid species_. we are then motivated to define the curvature energy of the lipid domain boundary to be @xmath4 where @xmath5 is the one - dimensional domain boundary contour on two - dimensional surface embedded in @xmath3 , @xmath6 is the geodesic curvature of the interface , @xmath7 is the spontaneous geodesic curvature of the lipid mixture to be separated , and @xmath8 is the geodesic curvature energy coefficient . this form of curvature energy is identical to the classical canham - evans - helfrich energy in which the integration in over the two - dimensional membrane surface and the mean curvature is adopted . the spontaneous geodesic curvature quantifies the geometry and molecular mismatch at the lipid interface . in the transitional region near the interface two species of lipids coexist in a hybrid state , and a free energy for the hybrid packing of the lipids ( denoted by the subscripts 1 and 2 below ) was proposed in recent theoretical studies @xcite : @xmath9 where @xmath10 is the length of the lipid chains in the transitional region and @xmath11 is the chain length in the equilibrium bulk . coefficients @xmath12 are the free energetic costs of mismatch between two species and their respective hybrids at the interface , and @xmath13 is the energetic cost of mismatch hybrid chains of different species . we choose @xmath14 for the purpose of simplification , but the generalization to distinct @xmath8 s can be made easily . the domain curvature can be related to the lipid geometrical properties by the following relation @xcite : @xmath15 where @xmath16 is the volume of the lipid chains , @xmath17 is the length that characterizes the molecular spacing of the lipid head groups , and @xmath18 is the head group spacing in the hybrid region near the interface . here the subtraction sign is chosen if the species is included in the microdomain , otherwise the addition sign shall be chosen . the chain length in the equilibrium bulk state , @xmath11 , can be computed by dividing the molecular volume using the head group area in the bulk state @xmath19 eqs.([eqn : vi]-[eqn : li0 ] ) allow us to represent the bending energy @xmath20 in eq.([eqn : bending_energy ] ) as a quadratic function of the geodesic curvature @xmath6 . the minimizer @xmath7 can be analytically calculated , and can be approximated to the linear order of @xmath21 and @xmath22 by @xmath23,\ ] ] where @xmath24 and @xmath25 . this approximate minimizer @xmath7 will appear in a more complicated form if @xmath26 , and calculations will show that @xmath7 is not sensitive to the value of @xmath27 in its biologically relevant range @xcite . the minimizer @xmath7 is identified as the _ spontaneous geodesic curvature _ of the corresponding binary lipid mixture .      for a given binary lipid mixture in bilayer membrane , one can work to minimize the energy in eq.([eqn : geocurv_eng ] ) to get the optimal boundary contour of surface lipid microdomains . this could be done , for example , by computing the shape derivative following the general analytical approach outlined in @xcite and using this shape derivative as the boundary force to drive the motion of the initial boundary contour . practical application of this approach , analytically or numerically , might be hindered by the topological changes of the boundary contour , during which there are moments that the surface is not smooth and its principal curvatures are not well defined . one can resort to phase field or level set methods with which the interface can be traced by evolving a function in a higher dimensional space . here we adopt the phase field approach , for which we define a phase field function @xmath28 associated with the boundary contour @xmath5 on the membrane surface @xmath29 as @xmath30 here @xmath31 is the signed geodesic distance from @xmath32 to the contour @xmath5 , @xmath33 in the interested domain enclosed by @xmath5 , @xmath34 on @xmath5 , and @xmath35 outside . correspondingly , @xmath36 and @xmath37 inside and outside , respectively , with the width of the transitional layer near @xmath5 determined by small @xmath38 . this function @xmath31 can be regarded as the generalization of the signed distance function used in euclidean space @xmath39 to the surface ; some of the critical critical properties of the latter can be shown true for @xmath31 . in appendix [ sect : append ] we provide details of relating interface normal and geodesic curvature to this signed geodesic distance . with this definition of @xmath28 we can compute the geodesic curvature @xmath6 of its arbitrary level set , including @xmath5 where @xmath40 . to facilitate the computation we first define @xmath41 which gives @xmath42 , ~ q''(x ) = -\frac{1 } { \epsilon } \tanh \left ( \frac{x}{\sqrt{2}\epsilon } \right ) \operatorname{sech}^2 \left ( \frac{x}{\sqrt{2}\epsilon } \right ) .\ ] ] then @xmath43 where @xmath44 and @xmath45 are surface gradient , surface divergence , and laplace - beltrami operators on @xmath29 , respectively . the proof of the last equality ( [ eqn : lap_phi ] ) can be found in appendix [ sect : append ] as well , where we also proved that the tangent normal vector field for the level sets of @xmath31 is @xmath46 and the geodesic curvature @xmath47 . we then compute @xmath48 it follows that the geodesic curvature is @xmath49 replacing @xmath50 with @xmath28 and noticing that @xmath51 ( to be proved in appendix [ sect : append ] ) , we finally get the phase field representation of the geodesic curvature : @xmath52 correspondingly , we define a phase field representation of the geodesic curvature energy on the entire surface @xmath29 , as follows : @xmath53 where @xmath54 . for the phase field function @xmath28 defined in eq.([eqn : phi_def ] ) one can show that this phase field curvature energy formulation is equivalent to eq.([eqn : geocurv_eng ] ) . this nature is analogous to the equivalence between the canham - helfrich - evans curvature energy and the phase field representation of the membrane elastic energy @xcite . the proof of the equivalence proceeds rather similarly , and is omitted here . our geodesic curvature energy model of surface microdomain formation can be regarded as a generalization of the mean - curvature driven membrane bending model from three - dimensional euclidean space to two - dimensional surface . and the other way around , the mean - curvature driven membrane bending model can be viewed as as a special case of the curvature energy driven model on manifolds to euclidean spaces . we expect our phase field based geodesic curvature model can be applied to other physical systems that occur on manifolds where the formation and the topological change of surface patterns are driven by appropriate curvatures . these may include the change of spacetime topology curved by matter and energy in general relativity @xcite . we follow the gradient flow of the curvature energy ( [ eqn : g_phase_field ] ) to find its minimizer @xmath28 whose zero level set shall give the the optimal lipid microdomain boundary . to derive the governing equation for the gradient flow we compute the first variation of @xmath55 with respect to @xmath28 : @xmath56 , \label{eqn : var_g}\ ] ] where @xmath57 the full expansion of this variation reads @xmath58 during the energy minimization the amounts of each species of lipids shall be conserved ; for this reason we consider the following constraint @xmath59 whose first variation with respect to @xmath28 is @xmath60 we then introduce the following governing equation for the gradient flow of @xmath28 @xmath61 where @xmath62 is the pseudo - time , @xmath13 is the diffusion coeffcient , and @xmath63 is a lagrangian multiplier used to ensure the conservation of @xmath28 . the representation of @xmath63 can be derived by integrating eq . ( [ eqn : gradient_flow_surface_pattern ] ) and noting that @xmath64 , thus @xmath65 and consequently , @xmath66 which yields a forth - order nonlinear surface diffusion equation . @xmath67 alternatively , one could derive a cahn - hilliard equation for the surface phase field function @xmath28 as @xmath68 which guarantees the conservation of @xmath28 and thus does not need a lagrangian multiplier . however , it involves a sixth order surface derivative and thus is more complicated when the equation is to be solved numerically on a discretized surface @xmath29 .    to begin the time discretization of eq.([eqn : gradient_flow_surface_pattern_2 ] ) we first average the function @xmath69 over the current and the next steps @xmath70 and @xmath71 , upon which the crank - nicolson scheme is applied : @xmath72 where the averaged function is defined by @xmath73 with @xmath74    to achieve convergence of this nonlinear implicit scheme at each time step , we define an inner iteration for computing @xmath75 . the solution during these inner iterations are denoted by @xmath76 , which is expected to converges to @xmath75 as @xmath77 . we then replace @xmath75 in eq.([eqn : scheme_1 ] ) with @xmath78 and @xmath79 to get @xmath80 where the new average function is defined by @xmath81 with @xmath82    to efficiently solve the nonlinear eq.([eqn : scheme_2 ] ) we split the average functions into linear and nonlinear components , with the linear component being the function of @xmath79 only and the nonlinear component being the function of @xmath83 and @xmath78 , as follows : @xmath84 the nonlinear implicit scheme is now given by @xmath85 which , after all terms of @xmath79 collected to the left - hand side , turns out to be @xmath86 \delta t + \phi_n \nonumber \\ = & \frac{1}{2 } \big ( \lambda(\psi_{m } ) + \lambda(\phi_n ) \big ) \delta t - \frac{\delta \gamma k \epsilon}{2 } \delta_s^2 \phi_n + \nonumber \\   & \delta t \gamma \delta_s \left ( \frac{k}{4 \epsilon } ( \psi_m^2 + \phi_n^2 -2)(\psi_m + \phi_n + 2 h_c \epsilon ) \right ) + \nonumber \\   & \delta t \gamma \delta_s \left ( \frac{k}{2\epsilon^2 } ( \psi_m^2 + \psi_m \phi_n + \phi_n^2 + h_c \epsilon(\psi_m + \phi_n ) \right ) \cdot \nonumber \\   & \left ( f(\psi_m ) + f(\phi_n ) \right )   + \phi_n.\end{aligned}\ ] ] this last equation is iterated over the inner index @xmath87 until the convergence is reached as signified by @xmath88 for some desired tolerance . @xmath78 is initialized as @xmath83 for inner iterations at each time step . the convergent @xmath89 will be passed to @xmath75 to advance the computation to the next time step . we adopt a @xmath90 interior penalty discontinuous galerkin ( dg ) method to solve the forth - order linear equation ( [ eqn : scheme_4 ] ) for @xmath79 . this method was initially developed for solve forth - order equations in two - dimensional euclidean space @xcite . when it is applied to solve equations on triangulated two - dimensional surface , an additional term arises on all edges , representing the mismatch of fluxes due to the different outer normal directions of two neighbouring triangles at the shared edge . the order of convergence of the surface dg method with respect to the mesh refinement matches that observed on the original two - dimensional numerical experiments . a large time step is allowed in the computational simulated presented below owning to the implicit time discretization . the inner iterations in general converge within five steps . the solutions of a surface cahn - hilliard equation in @xcite is based on the reduction to two second - order equations , one for the order parameter and the other for the chemical potential . the equations are then solved using standard surface finite element method along with an explicit time discretization , the latter suffers from a strict discrete energy stability condition and only very small time increments are allowed . an even larger time step might be chosen by applying a recently developed exponential time differencing method @xcite . with this method one needs to deliberately estimate the nonlinear component of @xmath91 at every step so its slowly varying part can be combined with the linear component of @xmath91 , possibly further improving the discrete energy stability of the scheme . in this section we shall apply the geodesic curvature energy variational model to simulate the lipid microdomain formation on three surfaces with different intrinsic geodesic curvatures . the dynamics of the microdomain formation in each set of simulations will be compared to those generated by the surface cahn - hilliard equation @xcite , the latter will produce large surface domain as a result of minimizing the arc length of the domain boundaries . we first consider a unit sphere on which a trigular mesh with @xmath92 approximately uniformly distributed nodes . we choose @xmath93 , a time increment @xmath94 , and a random field of @xmath28 . the results are compared side by side with those of the classical surface cahn - hilliard equation in fig . [ fig:1atm_3963 ] . using a x - means clustering method @xcite we are able to identify a number of microdomains whose radii are then calculated . the radius associated with each microdomain is approximately @xmath95 . this means the curvature is approximately @xmath96 , close to the specified spontaneous geodesic curvature . nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ]   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ] +   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ]   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ] +   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ]   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ] +   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ]   nodes show the formation of local microdomains predicted by the geodesic curvature energy ( left column ) and domain separation predicted by the allen - cahn equation associated with the classical ginzburg - landau energy ( right column ) from the same initial random field ( top row ) . sampling time from top to bottom : @xmath97 and @xmath98.,title="fig:",width=143 ]    the second set of simulations is conducted on a more complicated surface as shown in fig . [ fig:3atm ] . we choose this surface to be the molecular surface of three particles of unit radius respectively centered at @xmath99 and @xmath100 . such molecular surfaces are obtained by tracing the centers of spherical probe of the radius of the water molecule ( @xmath101 , we are using dimensionless length is this work , though ) as the probe is rolling over the surfaces of the molecules @xcite . this surface is quasi - uniformly meshed with @xmath102 nodes and we set @xmath103 and @xmath94 . starting with a random initial field we finally identified six microdomains using the x - mean clustering method at the equilibrium state , whose radii are estimated . as seen in fig . [ fig : domain_radii](left ) , the radii of the microdomains approximate the given spontaneous geodesic curvatures .     and @xmath98.,title="fig:",width=143 ]   and @xmath98.,title="fig:",width=143 ] +   and @xmath98.,title="fig:",width=143 ]   and @xmath98.,title="fig:",width=143 ] +   and @xmath98.,title="fig:",width=143 ]   and @xmath98.,title="fig:",width=143 ] +   and @xmath98.,title="fig:",width=143 ]   and @xmath98.,title="fig:",width=143 ]        in the last set of simulations we choose the molecular surface of six particles of unit radius respectively centered at @xmath104 and @xmath105 . the quai - uniform surface mesh has @xmath106 nodes and we choose @xmath107 and @xmath94 for the simulation . the dynamics of the microdomain formation is presented in fig . [ fig:6atm ] along with that of the domain separation described by the surface cahn - hilliard equation . one can see from fig . [ fig : domain_radii](right ) that the largest raft radius obtained by the simulation is about @xmath108 which means the curvature of that raft is about @xmath109 , a value close to given spontaneous geodesic curvature .     and @xmath110.,title="fig:",width=143 ]   and @xmath110.,title="fig:",width=143 ] +   and @xmath110.,title="fig:",width=143 ]   and @xmath110.,title="fig:",width=143 ] +   and @xmath110.,title="fig:",width=143 ]   and @xmath110.,title="fig:",width=143 ] +   and @xmath110.,title="fig:",width=143 ]   and @xmath110.,title="fig:",width=143 ]    the radii of the microdomains generated in our simulations are not exactly the given spontaneous geodesic curvature . rather they are distributed around the given curvature . apart from the numerical error in simulation and in x - mean clustering and radii estimate , this non - uniform distribution of domain radii is mostly related to the total quantity of the lipid phases in the initial random field . the initial quantity may not be exact to cover an integer number of microdomains with the given radius . however , the overall distribution of radii around the given radius of curvature demonstrates that our geodesic curvature model is capable of predicting the formation of microdomains that are caused by the geometrical and molecular mechanical mismatch of lipid mixtures . the predicted microdomains can be compared to the observed lipid rafts , and the boundaries of these microdomains can be identified to provide locations where specific proteins can aggregate . coupling of our model of geodesic curvature driven microdomains formation to the localization of proteins will provide a very useful quantitative technique for studying the crucial roles of these proteins in high - fidelity signal transmission in cells @xcite . in this work we investigated formation of microdomains on lipid membrane surfaces . we examine the geometrical and mechanical mismatch of lipids of different species on domain boundaries and introduce the geodesic curvature of these boundaries in the energy characterization of the equilibrium state . this gives rise to the geodesic curvature energy model of the surface microdomain formation . a phase field representation of the model is derived to ease the computational simulations of the model . we develop a @xmath90 interior penalty surface finite element method along with an explicit / implicit time splitting scheme for solving the forth order nonlinear partial differential equation on surfaces . numerical simulations on different surfaces demonstrate that our model and numerical techniques are able to produce the circular raft - like microdomains we are interested in . the generated microdomains have the radii corresponding to the given spontaneous geodesic curvature for the binary lipid mixture . there are accumulated evidences that lipid rafts aid in vesicle budding and the deformation of the membrane . currently many models are being produced to model vesicle budding and membrane deformation . this opens up a whole research area of connecting predicted microdomains using the developed geodesic curvature model to other biophysical processes . this work has been partially supported by the simons foundation and national institutes of health through the grant r01gm117593 as part of the joint dms / nigms initiative to support research at the interface of the biological and mathematical sciences . this section is devoted to the justification of a number of surface calculus identities used in the phase field formulation in this article . although the counterparts of these identities in euclidean spaces are well known , a rigorous derivation of these surface identities , to our knowledge , does not exist . we believe the justification of these identities is useful to scientists working on the various problems involving surface patterns and surface energies . since the phase field function @xmath111 is an algebraic function of the signed geodesic function @xmath31 , it is sufficient to work with @xmath31 in the derivation below .    let @xmath29 be an oriented smooth surface without boundary in three - dimensional euclidean space @xmath3 and let @xmath112 be a curve in @xmath29 , parameterized by the arc length @xmath113 . the darboux frame of the curve @xmath112 is defined by @xmath114 where @xmath115 are respectively the unit tangent vector , the unit normal vector , and the unit tangent normal vector . let @xmath31 be the signed geodesic function from the curve @xmath112 , following the convention that @xmath33 in the surface region of interests enclosed by @xmath112 , zero on @xmath112 , and positive outside . one can define the similar darboux frame on an arbitrary level set of @xmath31 on @xmath29 , giving rise to three vector fields @xmath116 and @xmath117 on @xmath29 . we shall have    the tangent normal vector field @xmath117 on @xmath29 is the normalized surface gradient of the signed geodesic distance function @xmath31 : @xmath118    let @xmath119 be an arbitrary point on @xmath29 and let @xmath120 be the level set of @xmath31 passing through @xmath119 for a constant @xmath121 , parameterized by the arc length @xmath113 . associated with @xmath112 at @xmath119 are the three vectors , @xmath122 defined by the darboux frame . let @xmath123 be the two - dimensional orthogonal coordinate of the tangent plane of the surface @xmath29 originated at @xmath119 , with @xmath124 aligned with @xmath125 . furthermore , we consider a local monge parameterization of @xmath29 around @xmath119 , with the coordinate variables @xmath126 . the linear map @xmath127 from @xmath128 to @xmath123 is such that @xmath129 the signed geodesic distance function @xmath31 can be given locally by either @xmath130 or @xmath131 . the specific level set @xmath132 can be locally given by @xmath133 in the @xmath128 plane . it follows that @xmath134 using the mapping @xmath127 we will have at point @xmath119 @xmath135 suggesting that @xmath136 is orthogonal to the tangent vector @xmath125 of @xmath137 at @xmath119 . since @xmath136 is on the tangent plane it is orthogonal to surface with the normal vector @xmath138 , and hence it must be in the direction of the tangent normal vector @xmath139 . finally we check @xmath140 . since @xmath132 is a level set and @xmath124 is aligned with the tangent direction of this level set it follows that @xmath141 . now that @xmath142 is in the direction of the tangent norm , and given that @xmath132 is a level set of the signed geodesic distance function , @xmath142 is indeed in the direction of the geodesic passing through @xmath119 , thus @xmath143 . this suggests that @xmath144 and the theorem follows . after the establishment of the tangent normal vectors field @xmath117 on @xmath29 , we can compute the geodesic curvature of any level set of the signed geodesic distance function @xmath31 by using the surface divergence of @xmath117 . let @xmath112 be a level set @xmath145 , parameterized by the arc length @xmath113 . one can associate a frenet - serret frame to @xmath112 , defined by @xmath146 where @xmath147 are respectively the unit tangent vector , the frenet normal vector , and the frenet binormal vector . the frenet normal vector is also called the principal normal vector . the curvature @xmath8 of the curve @xmath112 is @xmath148 the plane spanned by the orthogonal vectors @xmath149 is called osculating plane . let @xmath150 be the angle between the osculating plane and the tangent plane . the geodesic curvature is defined to be the projection of the curvature to the tangent plane , i.e. , @xmath151 it can be shown that the geodesic curvature is indeed the ordinary curvature of the curve obtained by projecting @xmath112 on to the tangent plane ( p.261 , @xcite ) . this last fact will be used to prove the following theorem :    the geodesic curvature @xmath152 of the space curve @xmath112 given by the level set @xmath153 of the signed geodesic distance function equals to the surface divergence of the tangent normal vector field @xmath117 : @xmath154    let @xmath123 be the orthogonal coordinate of a local monge parameterization of the surface @xmath29 near the point @xmath155 such that @xmath156 are respectively in the directions of the tangent and tangent normal vectors @xmath157 at @xmath155 , and @xmath155 itself is mapped to the origin of the @xmath123 plane . let the height of the surface be locally @xmath158 . under this parameterization the space curve @xmath112 will be mapped to a curve on @xmath123 plane , denoted by @xmath159 . the space curve @xmath112 can be parameterized by @xmath124 only by setting @xmath160 accordingly , the tangent vector @xmath161 can be given as a function of @xmath124 by @xmath162 since this tangent vector is identical to the element vector @xmath163 at @xmath164 by the choice of the monge parameterization , it follows that at @xmath155 @xmath165    we compute the geodesic curvature of the surface curve @xmath112 at @xmath155 by evaluating the curvature of the plane curve @xmath166 at @xmath164 : @xmath167 where the negative sign is assigned so the orientation of the plane curve gives outer - pointing normal vector that is consistent with the orientation of @xmath123 coordinate . it remains to show that @xmath168 . we represent the tangent normal vector field @xmath169 using the local monge parameterization . for any point @xmath170 on the surface curve @xmath112 we have @xmath171 for some @xmath124 and then @xmath172 where @xmath173 to compute surface divergence we notice that @xmath174 where @xmath175 are the components of the vector @xmath117 respectively in @xmath176 directions . on one hand we shall have @xmath177 at @xmath155 where @xmath178 for that @xmath179 achieves its maximum value as a component of a unit normal vector owning to @xmath180 at @xmath155 . on the other hand , @xmath181 - \\    & \frac{f_{\xi}}{f^2 } \big [ 2(h_{\xi \xi } + h_{\xi \eta } f'(\xi ) ) h_{\xi } + 2 ( h_{\xi \eta } + h_{\eta \eta } f'(\xi ) ) h_{\eta } \big ]   \big [ ( h_{\eta}h_{\xi } + h^2_{\eta } f'(\xi ) ) - f'(\xi ) \big ] \end{aligned}\ ] ] it follows again from @xmath180 that @xmath182 and @xmath183 as requested .    the last fact to justify is the identity eq.([eqn : lap_phi ] ) . recalling that for a scalar @xmath184 and a vector field @xmath185 on the surface @xmath29 with outer normal vector @xmath186 it holds true that @xmath187 \\    & =   a ~ n \cdot \curl ( n \times v ) + n \cdot \big [ \nabla a \times ( n \times v ) \big ] \\    & =   a \nabla_s v + n \cdot \big [ ( \nabla a \cdot v ) n - ( \nabla a \cdot n ) v \big ] , \end{aligned}\ ] ] we shall have @xmath188 99 g.  van meer , d.  r. voelker , g.  w. feigenson , membrane lipids : where they are and how they behave . , nat . rev . cell biol . 9 ( 2008 ) 112124 . [ ]        s.  l. veatch , s.  l. keller , separation of liquid phases in giant vesicles of ternary mixtures of phospholipids and cholesterol , biophys . j. 85  ( 5 ) ( 2003 ) 3074  3083 . http://dx.doi.org/http://dx.doi.org/10.1016/s0006-3495(03)74726-2 [ ] . a.  kusumi , c.  nakada , k.  ritchie , k.  murase , k.  suzuki , h.  murakoshi , r.  s. kasai , j.  kondo , t.  fujiwara , paradigm shift of the plasma membrane concept from the two - dimensional continuum fluid to the partitioned fluid : high - speed single - molecule tracking of membrane molecules , annu . 34 ( 2005 ) 351378 . http://dx.doi.org/10.1146/annurev.biophys.34.040204.144637 [ ] . f.  s. vieira , g.  correa , m.  einicker - lamas , r.  coutinho - silva , host - cell lipid rafts : a safe door for micro - organisms ? , biology of the cell 102  ( 7 ) ( 2010 ) 391407 . http://dx.doi.org/10.1042/bc20090138 [ ] . t.  wang , s.  d. cady , m.  hong , nmr determination of protein partitioning into membrane domains with different curvatures and application to the influenza m2 peptide , biophys . j. 102  ( 4 ) ( 2012 ) 787  794 . http://dx.doi.org/http://dx.doi.org/10.1016/j.bpj.2012.01.010 [ ] . w.  b. huttner , j.  zimmerberg , implications of lipid microdomains for membrane curvature , budding and fission : commentary , current opinion in cell biology 13 ( 2001 ) 478484 . http://dx.doi.org/10.1016/s0955-0674(00)00239-8 [ ] . l.  d. renner , d.  b. weibel , cardiolipin microdomains localize to negatively curved regions of escherichia coli membranes , p. natl . usa 108  ( 15 ) ( 2011 ) 62646269 . ryu , d.  yoo , n.  j. wittenberg , l.  r. jordan , s .- d . lee , a.  n. parikh , s .- h . oh , lipid membrane deformation accompanied by disk - to - ring shape transition of cholesterol - rich domains , j. am . 137  ( 27 ) ( 2015 ) 86928695 . http://dx.doi.org/10.1021/jacs.5b04559 [ ] . a.  hinderliter , p.  f.  f. almeida , c.  e. creutz , r.  l. biltonen , domain formation in a fluid mixed lipid bilayer modulated through binding of the c2 protein motif , biochemistry 40  ( 13 ) ( 2001 ) 41814191 , pmid : 11300799 . http://dx.doi.org/10.1021/bi0024299 [ ] . t.  yeung , g.  e. gilbert , j.  shi , j.  silvius , a.  kapus , s.  grinstein , membrane phosphatidylserine regulates surface charge and protein localization , science 319  ( 5860 ) ( 2008 ) 210213 . [ ] . d.  v. nicolau , k.  burrage , r.  g. parton , j.  f. hancock , identifying optimal lipid raft characteristics required to promote nanoscale protein - protein interactions on the plasma membrane , mol . cell biol . 26 ( 2006 ) 313323 . http://dx.doi.org/10.1128/mcb.26.1.313-323.200 [ ] . x.  wang , q.  du , modelling and simulations of multi - component lipid membranes and open membranes via diffuse interface approaches , j. math . ( 2008 ) 347371 . http://dx.doi.org/10.1007/s00285-007-0118-2 [ ] . c.  yuan , j.  furlong , p.  burgos , l.  j. johnston , the size of lipid rafts : an atomic force microscopy study of ganglioside gm1 domains in sphingomyelin / dopc / cholesterol membranes , biophys . j. 82 ( 2002 ) 25262535 . http://dx.doi.org/10.1016/s0006-3495(02)75596-3 [ ] . q.  du , l.  ju , l.  tian , finite element approximation of the cahn - hilliard equation on surfaces , comput . methods appl . 200  ( 29?32 ) ( 2011 ) 2458  2470 . http://dx.doi.org/http://dx.doi.org/10.1016/j.cma.2011.04.018 [ ] . b.  a. camley , f.  l.  h. brown , dynamic simulations of multicomponent lipid membranes over long length and time scales , phys . ( 2010 ) 148102 . http://dx.doi.org/10.1103/physrevlett.105.148102 [ ] . t.  witkowski , r.  backofen , a.  voigt , the influence of membrane bound proteins on phase separation and coarsening in cell membranes , phys . ( 2012 ) 1450914515 . http://dx.doi.org/10.1039/c2cp41274h [ ] . j.  c. shillcock , r.  lipowsky , equilibrium structure and lateral stress distribution of amphiphilic bilayers from dissipative particle dynamics simulations , j. chem . phys . 117  ( 10 ) ( 2002 ) 50485061 . http://dx.doi.org/http://dx.doi.org/10.1063/1.1498463 [ ] . h.  li , a.  a. gorfe , aggregation of lipid - anchored full - length h - ras in lipid bilayers : simulations with the martini force field , plos one 8  ( 7 ) ( 2013 ) e71018 . q.  du , c.  liu , r.  ryham , x.  wang , modeling the spontaneous curvature effects in static cell membrane deformations by a phase field formulation , comm . anal . 4 ( 2005 ) http://dx.doi.org/10.3934/cpaa.2005.4.537 [ ] . q.  du , c.  liu , x.  wang , simulating the deformation of vesicle membranes under elastic bending energy in three dimensions , j. comput . 212  ( 2 ) ( 2006 ) 757  777 . http://dx.doi.org/doi : 10.1016/j.jcp.2005.07.020 [ ] . r.  brewster , s.  a. safran , line active hybrid lipids determine domain size in phase separation of saturated and unsaturated lipids , biophys . j. 98  ( 6 ) ( 2010 ) l21l23 . http://dx.doi.org/http://dx.doi.org/10.1016/j.bpj.2009.11.027 [ ] . s.  c. brenner , s.  gu , t.  gudi , l.  yeng sung , a quadratic c0 interior penalty method for linear fourth order boundary value problems with boundary conditions of the cahn - hilliard type , siam j. numer . 50  ( 4 ) ( 2012 ) 20882110 . http://arxiv.org/abs/http://dx.doi.org/10.1137/110847469 [ ] , http://dx.doi.org/10.1137/110847469 [ ] . x.  wang , l.  ju , q.  du , efficient and stable exponential time differencing runge - kutta methods for phase field elastic bending energy models , j. comput . ( 2016 ) 21  39 . http://dx.doi.org/http://dx.doi.org/10.1016/j.jcp.2016.04.004 [ ] . d.  pelleg , a.  w. moore , x - means : extending k - means with efficient estimation of the number of clusters , in : proceedings of the seventeenth international conference on machine learning , icml 00 , morgan kaufmann publishers inc . , san francisco , ca , usa , 2000 , pp . 727734 . h.  j. hocker , k .- j . cho , c .- y . k. chen , n.  rambahal , s.  r. sagineedu , k.  shaari , j.  stanslas , j.  f. hancock , a.  a. gorfe , andrographolide derivatives inhibit guanine nucleotide exchange and abrogate oncogenic ras function , p. natl . usa 110  ( 25 ) ( 2013 ) 1020110206 . http://arxiv.org/abs/http://www.pnas.org/content/110/25/10201.full.pdf [ ] , http://dx.doi.org/10.1073/pnas.1300016110 [ ] .
lipid bilayer membranes are not uniform and clusters of lipids in a more ordered state exist within the generally disorder lipid milieu of the membrane . these clusters of ordered lipids microdomains are now referred to as lipid rafts . 
 recent reports attribute the formation of these microdomains to the geometrical and molecular mechanical mismatch of lipids of different species on the boundary . here 
 we introduce the geodesic curvature to characterize the geometry of the domain boundary , and develop a geodesic curvature energy model to describe the formation of these microdomains as a result of energy minimization . 
 our model accepts the intrinsic geodesic curvature of any binary lipid mixture as an input , and will produce microdomains of the given geodesic curvature as demonstrated by three sets of numerical simulations . 
 our results are in contrast to the surface phase separation predicted by the classical surface cahn - hilliard equation , which tends to generate large domains as a result of the minimizing line tension . 
 our model provides a direct and quantified description of the structure inhomogeneity of lipid bilayer membrane , and can be coupled to the investigations of biological processes on membranes for which such inhomogeneity plays essential roles . 
 lipid bilayer membrane ; lipid domain ; geodesic curvature ; surface phase separation ; phase field method ; variational principle ; numerical simulations 53b10 , 65d18 , 92c15