Finite element simulation of multi‐component vesicle morphologies undergoing phase separation

Morphological changes in lipid bilayer vesicles are due to phase transitions and surface deformations occurring in unison. We present a dynamic chemo‐mechanical finite element model to study such changes. This is achieved by coupling two fourth order partial differential equations (PDEs): The Cahn‐Hilliard [2] mass balance equation based on irreversible thermodynamics, and the Kirchhoff‐Love [1] rotation free thin shell equation. The Helmholtz free energy consists of Helfrich energy to model elastic bending, and also includes in‐plane elastic energy with a finite shear modulus for the purpose of regularization. The geometry is discretized by C1‐continuous NURBS shape functions. An implicit second order accurate generalized‐α scheme is used for time integration. Newton‐Raphson iterations are utilized to solve the resulting linearized weak form monolithically. The proposed formulation is illustrated by a numerical example.


Introduction
In this work, we consider multicomponent lipid bilayers and study the phase transitions using a modified Cahn-Hilliard equation. Moreover, lipid bilayers exhibit morphological changes during phase separation. Mechanically, lipid bilayers behave like a fluid in-plane and like a solid out-of plane. We model the mechanical behaviors of lipid bilayers as a Kirchhoff-Love shell based on in-plane visco-elasticity and out-of plane elasticity. Though a shear modulus is incorporated into the formulation specifically for the regularization of the in-plane balance of linear momentum equations, in reality cell membranes do offer considerable amount of shear stiffness when associated with the underlying cortex consisting of cytoskeletal filaments. Finally, the Helfrich model [1] is used for modeling elastic bending, including an important contribution from the preference of mean curvature of the individual phases.

Coupled phase transition and surface deformation
The lipid-bilayer surface S is described by the mapping x = x(ξ α , t), where ξ α , α = 1, 2 are the curvilinear coordinates. We thus consider a 2D surface embedded in 3D space. For a more detailed description of the surface kinematics refer to [1].
If ρ 1 and ρ 2 are the mass densities per unit area of the two phases, the rate of change of phase field parameter, φ := ρ 1 /(ρ 1 + ρ 2 ), is via the mass balance as [2] Where j α = j · a α = a αβ j β is the contravariant component of diffusive flux vector j and ρ = ρ 1 + ρ 2 . The diffusive flux component is a function of chemical potential µ c and can be written as is the degenerate mobility and J is the surface stretch. Eq. (1) is the governing strong form for the phase field and hence can be used to solve for the distribution of φ on the surface at time t.
For the mechanical counterpart, the governing strong form for position x follows from the balance of linear momentum as Where T α is the stress vector, f denotes the body force and v is the velocity. The stress vector is obtained from stress tensor σ through the Cauchy's formula, T α = σ T a α . Using rotation free shell theory, the stress vector can be expressed as a function of stress (σ αβ ) and bending moment (M αβ ) components.

of 3
Section 7: Coupled problems The Helmholtz free energy per reference area (Ψ) is considered to be a function of the surface metric (a αβ ), covariant components of curvature tensor (b αβ ), temperature (T ), phase field parameter (φ), and its surface gradient (∇ s φ). It has three contributions, namely, from membrane deformations, bending and Cahn-Hilliard resulting in the form [3] Here, K, G, c and H 0 represent 2D bulk, shear, bending moduli and spontaneous curvature and are all functions of φ. H is the mean curvature. N , k B , √ λ and ω are constants and denote the number of molecules per reference area, length scale of the phase interface and bulk energy respectively.
Based on the Helmholtz free energy function the chemical potential, stress, and bending moment components follows as Here η represents surface viscosity. The elastic part of the free energy, Ψ el = Ψ memb + Ψ bend , contributes to the chemical potential while Ψ CH contributes to the stress. The surface divergence theorem is applied twice to both Eqs. (1) and (2) to derive the weak form and hence the finite element forces. An adaptive time stepping based on generalized-α scheme [3] is used for time integration. The linearized set of equations are solved monolithically and using Newton-Raphson iterations.

Numerical example
To demonstrate the above formulation, consider an unit sphere with constant internal pressure p int = 0.05 EL −1 0 and initial φ = 0.33 + φ r where φ r ∈ [−0.05, 0.05] is randomly distributed. The geometry is discretized with 2472 cubic NURBS elements and the parameter values used for pure phase states φ = 0 and 1 are listed in Tab. 1. Some of the parameters are expressed in terms of surface Young's modulus E (force per length) and Poisson's ratio ν, which are chosen as E = N ω and ν = 0.3. Fig. 1 shows the initial, an intermediate and the final state of the sphere. The result obtained is similar to axi-symmetric result of [4], a shape also observed in experiments [5]. In order to extend the simulations to viscous fluids with G i = 0, surface ALE based on [6] is required. This will be addressed in future work.