The irreversible thermodynamics of curved lipid membranes

TL;DR

This paper develops a thermodynamic framework for curved lipid membranes, analyzing how elastic bending couples with processes like lipid flow, phase transitions, and protein interactions, providing governing equations and entropy production insights.

Contribution

It introduces a comprehensive irreversible thermodynamics model for curved lipid membranes, linking elastic and dissipative processes with explicit equations and boundary conditions.

Findings

Derived entropy production forms for membrane processes

Identified thermodynamic forces and fluxes in membranes

Provided governing equations for membrane dynamics

Abstract

The theory of irreversible thermodynamics for arbitrarily curved lipid membranes is presented here. The coupling between elastic bending and irreversible processes such as intra-membrane lipid flow, intra-membrane phase transitions, and protein binding and diffusion is studied. The forms of the entropy production for the irreversible processes are obtained, and the corresponding thermodynamic forces and fluxes are identified. Employing the linear irreversible thermodynamic framework, the governing equations of motion along with appropriate boundary conditions are provided.

Tables

Table 1. Table 1: The partial derivatives of the areal mass density ρ 𝜌 \rho , mean curvature H 𝐻 H , and Gaussian curvature K 𝐾 K with respect to the metric tensor a α ​ β subscript 𝑎 𝛼 𝛽 a_{\alpha\beta} and curvature tensor b α ​ β subscript 𝑏 𝛼 𝛽 b_{\alpha\beta} . Each table entry is the partial derivative of the column header with respect to the row header. For example, the third column of the first row indicates ∂ K / ∂ a α ​ β = − K ​ a α ​ β 𝐾 subscript 𝑎 𝛼 𝛽 𝐾 superscript 𝑎 𝛼 𝛽 \partial K/\partial a_{\alpha\beta}=-Ka^{\alpha\beta} .

	$\rho$	$H$	$K$
$a_{\alpha \beta }$	$-\frac{1}{2}\rho a^{\alpha \beta }$	$-\frac{1}{2}{b}^{\alpha \beta }$	$-K{a}^{\alpha \beta }$
$b_{\alpha \beta }$	0	$\frac{1}{2}{a}^{\alpha \beta }$	${\overline{b}}^{\alpha \beta }$

Table 2. Table 2: The chemical species in our example membrane system undergoing protein binding and unbinding reactions, indexed by k 𝑘 k .

k	Species
1	DOPC
2	PI(4,5)P2
3	PI(4,5)P2–epsin-1

Table 3. Table 3: New quantities relevant for reactions within the membrane system.

Symbol	Units	Description
$n_k$	mol / m2	concentration of component $k$.
$n_p^b$	mol / m3	concentration of proteins in the surrounding fluid.
${ℳ}_k$	g / mol	Molar mass of component $k$.
${ℳ}_p$	g / mol	Molar mass of epsin-1 protein.
$ℛ$	mol / m2 s	Rate of reaction per unit area.
$\overrightarrow{ℛ}$, $\overleftarrow{ℛ}$	mol / m2 s	Forward and reverse reaction rates, respectively.
$\overrightarrow{k}$	m3 / mol s	Forward reaction rate constant.
$\overleftarrow{k}$	1 / s	Reverse reaction rate constant.
${\alpha }_k$	–	The coefficient of the $k^{\mathrm{th}}$ species in the chemical reaction.
${\alpha }_p$	–	The coefficient of the epsin-1 protein in the chemical reaction.

Equations

$$\bm{x}=\bm{x}\big{(}\theta^{\alpha},t\big{)}~{}.$$

$${\mathbf{a}}_{\alpha }:={\mathbf{x}}_{,\alpha }.$$

$$\mathbf{n}:=\frac{{\mathbf{a}}_1\times {\mathbf{a}}_2}{|{\mathbf{a}}_1\times {\mathbf{a}}_2|}.$$

$${\mathbf{a}}^{\alpha }\cdot {\mathbf{a}}_{\beta }={\delta }_{\beta }^{\alpha },$$

$$a_{\alpha \beta }:={\mathbf{a}}_{\alpha }\cdot {\mathbf{a}}_{\beta }$$

$$a^{\alpha\beta}:=\bm{a}^{\alpha}\cdot\bm{a}^{\beta}=\big{(}a_{\alpha\beta}\big{)}^{-1}~{}.$$

$${\mathbf{a}}_{\alpha }$$

$${\mathbf{a}}^{\alpha }$$

$$\mathbf{h}=h^{\alpha }{\mathbf{a}}_{\alpha }+h\mathbf{n}=h_{\alpha }{\mathbf{a}}^{\alpha }+h\mathbf{n},$$

$$h^{\alpha }$$

$$h_{\alpha }$$

$${{\sigma }^{\alpha }}_{\beta }$$

$${\sigma }^{\alpha \beta }$$

$$\mathbf{\tau }:={\mathbf{a}}_{\alpha }\frac{\mathrm{d}{\theta }^{\alpha }}{\mathrm{d}\ell }$$

$$\mathbf{\nu }:=\mathbf{\tau }\times \mathbf{n}.$$

$$\mathbf{\nu }={\nu }^{\alpha }{\mathbf{a}}_{\alpha }={\nu }_{\alpha }{\mathbf{a}}^{\alpha }$$

$$\mathbf{\tau }={\tau }^{\alpha }{\mathbf{a}}_{\alpha }={\tau }_{\alpha }{\mathbf{a}}^{\alpha }.$$

$${\mathbf{a}}_{\alpha }={\nu }_{\alpha }\mathbf{\nu }+{\tau }_{\alpha }\mathbf{\tau }.$$

$$\mathbf{i}$$

$$1$$

$$b_{\alpha \beta }:=\mathbf{n}\cdot {\mathbf{x}}_{,\alpha \beta },$$

$$H:=\frac{1}{2}{a}^{\alpha \beta }{b}_{\alpha \beta }$$

$$K:=\frac{1}{2}{\epsilon }^{\alpha \beta }{\epsilon }^{\lambda \mu }{b}_{\alpha \lambda }{b}_{\beta \mu },$$

$${\bar{b}}^{\alpha \beta }:=2H{a}^{\alpha \beta }-b^{\alpha \beta },$$

$${\Gamma }_{\lambda \mu }^{\alpha }:=\frac{1}{2}{a}^{\alpha \delta }\left(a_{\delta \lambda ,\mu }+a_{\delta \mu ,\lambda }-a_{\lambda \mu ,\delta }\right).$$

$$v_{;\alpha }^{\lambda }$$

$$v_{\lambda ;\alpha }$$

$${\mathbf{a}}_{\beta ;\alpha }$$

$${\mathbf{n}}_{,\alpha }$$

$$J:=\frac{\text{d}a}{\text{d}A}.$$

$${\int }_{\mathcal{P}}f\mathrm{d}a={\int }_{{\mathcal{P}}_0}fJ\mathrm{d}A,$$

$$\frac{\text{d}}{\text{d}t}(\cdot ):=(\cdot {)}_{,t}+v^{\alpha }(\cdot {)}_{,\alpha }.$$

$$\mathbf{v}:=\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\dot{\mathbf{x}}=v\mathbf{n}+v^{\alpha }{\mathbf{a}}_{\alpha }.$$

Full text

‡  amaresh.sahu@berkeley.edu

§  sauer@aices.rwth-aachen.de

†  kranthi@berkeley.edu

1. 1 Introduction

2. 2 Kinematics

3. 3 Intra-membrane Flow and Bending

4. 3.1 Balance Laws

5. 3.1.1 Mass Balance

6. 3.1.2 Linear Momentum Balance

7. 3.1.3 Angular Momentum Balance

8. 3.1.4 Mechanical Power Balance

9. 3.2 Thermodynamics

10. 3.2.1 First Law—Energy Balance

11. 3.2.2 Entropy Balance & Second Law

12. 3.2.3 Choice of Thermodynamic Potential

13. 3.3 Constitutive Relations

14. 3.3.1 General Thermodynamic Variables

15. 3.3.2 Helmholtz Free Energy—Change of Variables

16. 3.3.3 Helfrich Energy Density

17. 3.4 Equations of Motion

18. 3.5 Boundary and Initial Conditions

19. 3.6 Coupling to Bulk Fluid

20. 4 Intra-membrane Phase Transitions

21. 4.1 Kinematics

22. 4.2 Balance Laws

23. 4.2.1 Mass Balance

24. 4.2.2 Linear Momentum Balance

25. 4.3 Thermodynamics

26. 4.3.1 First Law—Energy Balance

27. 4.3.2 Entropy Balance & Second Law

28. 4.4 Constitutive Relations

29. 4.4.1 General Thermodynamic Variables

30. 4.5 Application to Lo–Ld Phase Transitions

31. 4.5.1 Helmholtz Free Energy

32. 4.5.2 Equations of Motion

33. 4.5.3 Boundary and Initial Conditions

34. 5 Peripheral Proteins—Chemical Reactions

35. 5.1 Chemistry

36. 5.2 Balance Laws

37. 5.2.1 Mass Balance

38. 5.2.2 Linear Momentum Balance

39. 5.2.3 Angular Momentum Balance

40. 5.3 Thermodynamics

41. 5.3.1 First Law—Energy Balance

42. 5.3.2 Entropy Balance & Second Law

43. 5.3.3 Choice of Thermodynamic Potential

44. 5.4 Constitutive Relations

45. 5.4.1 General Thermodynamic Variables

46. 5.4.2 Helmholtz Free Energy—Change of Variables

47. 5.4.3 Helfrich Energy Density

48. 5.5 Equations of Motion

49. 5.6 Boundary and Initial Conditions

50. 6 Conclusions

51. A Surface parametrizations

52. A.1 Convected Coordinates

53. A.2 Material Time Derivative of In-Plane Basis Vectors

List of important symbols

[TABLE]

1 Introduction

In this paper we develop an irreversible thermodynamic framework for arbitrarily curved lipid membranes to determine their dynamical equations of motion. Using this framework, we find relevant constitutive relations and use them to understand how bending and intra-membrane flows are coupled. We then extend the model to include multiple transmembrane species which diffuse within the membrane, and learn how phase transitions are coupled to bending and flow. Finally, we model the binding and unbinding of surface proteins and their diffusion along the membrane surface.

Biological membranes comprised of lipids and proteins make up the boundary of the cell, as well as the boundaries of internal organelles such as the nucleus, endoplasmic reticulum, and Golgi complex. Lipid membranes and their interactions with proteins play an important role in many cellular processes, including endocytosis [1, 2, 3, 4, 5, 6, 7, 8, 9], exocytosis [10, 11], vesicle formation [12], intra-cellular trafficking [13], membrane fusion [14, 15, 11], and cell-cell signaling and detection [16, 17, 18].

Protein complexes that have a preferred membrane curvature can interact with the membrane surface and induce bending [19], important in processes where coat proteins initiate endocytosis [5, 7, 9, 6, 4] and BAR proteins sense and regulate membrane curvature [3, 20, 21, 2]. In all of these processes, lipid membranes undergo morphological changes in which phospholipids flow to accommodate the shape changes resulting from protein-induced curvature. These phenomena include both the elastic process of bending and irreversible processes such as lipid flow.

Another important phenomena in many biological membrane processes is the diffusion of intra-membrane species such as proteins and lipids to form heterogeneous domains. For example, T cell receptors are known to form specific patterns in the immunological synapse when detecting antigens [17, 18]. In artificial giant unilamellar vesicles, a phase transition between liquid-ordered (Lo) and liquid-disordered (Ld) membrane phases has been well-characterized [22, 23, 24]. Such phase transitions have also been observed on plasma membrane vesicles [25]. Furthermore, morphological shape changes in which either the Lo or the Ld phase domains bulge out to reduce the line tension between the two phases have been observed [22, 26]. These phenomena clearly indicate the coupling between elastic membrane bending and irreversible processes such as diffusion and flow, which must be understood to explain the formation of tubes, buds, and invaginations observed in various biological processes [26, 27, 28, 29].

The final phenomena of interest is the binding and unbinding of proteins to and from the membrane surface, and the diffusion of proteins once they are bound. Protein binding and unbinding reactions are irreversible processes and are ubiquitous across membrane-mediated phenomena [5, 6, 8, 7, 12, 14, 15, 19]. As an example, epsin-1 proteins can bind to specific membrane lipids during the early stages of endocytosis and induce bending [27, 7]. Moreover, antigen detection by T cells can be sensitive to the kinetic rates of T cell receptor binding and unbinding [30, 31, 32]. The kinetic binding of proteins also plays a crucial role in viral membrane fusion [14], where proteins and membranes are known to undergo kinetically restricted conformational changes in the fusion of influenza [33, 34, 35] and HIV [15]. The case of HIV is particularly interesting, as fusion proteins primarily reside at the interface between Lo and Ld regimes and fusion is believed to be favorable because it reduces the total energy due to line tension between the two phases [15].

All of the above phenomena involve elastic bending being fully coupled with the irreversible processes of lipid flow, the diffusion of lipids and proteins, and the surface binding of proteins. Comprehensive membrane models which include these effects are needed to fully understand the complex physical behavior of biological membranes. Our work entails developing a non-equilibrium thermodynamic framework that incorporates these processes.

Previous theoretical developments have modeled a range of lipid bilayer phenomena. The simplest models apply the theory of elastic shells [36] and model membranes with an elastic bending energy given by Canham [37] and Helfrich [38]. Many studies focus on solving for the equations of motion for simple membrane geometries such as the deviations from flat planes [39, 40, 41, 42, 43, 44] and cylindrical or spherical shells [42, 45, 46, 47]. Some of these works also model the coupling between elastic effects and either inclusions [40, 48, 49], surrounding protein structures [44], or fluid flow [42, 43, 50] on simple geometries. More general geometric frameworks based on theories of elastic shells were established to model lipid membranes of arbitrary geometry. Such models used variational methods to determine the constitutive form of the stress components [51, 52, 53, 54, 55, 56]. Models developed from variational methods have been built upon to include protein-induced curvature [57], viscosity [53], and edge effects [58], and are able to describe various membrane processes.

While the formulation of models from variational methods is theoretically sound, the techniques involved are not easily extendable to model the aforementioned coupling between fluid flow and protein-induced spontaneous curvature, phase transitions, or the binding and unbinding of proteins on arbitrary geometries. Recently, membrane models have been developed by using fundamental balance laws and associated constitutive equations [59, 60, 61, 62, 63, 51, 64]. In addition to reproducing the results of variational methods, they have had great success in understanding the specific effects of protein-induced curvature on membrane tension [65, 66] and simulating non-trivial membrane shapes [62]. However, comprehensive models including all of the irreversible phenomena mentioned thus far and their coupling to bending have not been developed for arbitrary geometries.

In this work we develop the general theory of irreversible thermodynamics for lipid membranes, inspired by the classical developments of irreversible thermodynamics by Prigogine [67] and de Groot & Mazur [68]. While these classical works [67, 68] are for systems modeled using Cartesian coordinates, developing this procedure for two-dimensional lipid membranes is difficult because lipid membranes bend elastically out-of-plane and behave as a fluid in-plane. As a consequence, a major complexity arises because the surface on which we apply continuum and thermodynamic balance laws is itself curved and deforming over time, thereby requiring the setting of differential geometry. We address these issues systematically.

The following aspects are new in this work:

A general irreversible thermodynamic framework is developed for arbitrarily curved, evolving lipid membrane surfaces through the fundamental balance laws of mass, linear momentum, angular momentum, energy, and entropy, as well as the second law of thermodynamics. 2. 2.

The contributions to the total entropy production are found and the viscous contribution agrees with earlier variational approaches [53] as well as a balance law formulation using an interfacial flow-based result [60]. 3. 3.

The thermodynamic framework is extended to model membranes with multiple species, the phase transitions between Lo and Ld domains, and their coupling to fluid flow and bending. 4. 4.

The model is expanded to incorporate the coupling between protein binding, diffusion, and flow, and the thermodynamic driving force governing protein binding and unbinding is determined.

Our paper is organized as follows: We develop our model by using the fundamental balance laws of mass, momentum, energy, and entropy to determine the equations of motion governing membrane behavior. We then apply irreversible thermodynamics to determine appropriate constitutive relations, which describe how membrane energetics affect dynamics. Section 2 reviews concepts from differential geometry which are necessary to describe membranes of arbitrary shape and presents general kinematic results. Section 3 models a single-component lipid membrane with viscous in-plane flow, elastic out-of-plane bending, inertia, and protein-induced or lipid-induced spontaneous curvature. In Section 4 we extend the model to include multiple lipid components and determine the equations of motion when phase transitions are possible. In Section 5, we model the binding and unbinding of proteins to and from the membrane surface. Throughout all of these sections, we find membrane phenomena are highly coupled with one another. For each of Sections 3, 4, and 5, we end by giving the expressions for the stresses and moments, writing the equations of motion, and providing possible boundary conditions to solve associated initial-boundary value problems. We conclude in Section 6 by including avenues for future work, both in advancing the theory and in developing computational methods.

2 Kinematics

We begin by reviewing concepts from differential geometry, as presented in [69], which are essential in describing the shape of the membrane and its evolution over time. We model the phospholipid bilayer as a single differentiable manifold about the membrane mid-plane, implicitly making a no-slip assumption between the two sheets of the bilayer. We will follow a similar notation to that presented in [60, 62, 70].

Consider a two dimensional membrane surface $\mathcal{P}$ embedded in Euclidean 3-space $\mathbb{R}^{3}$. The membrane position $\bm{x}$ is a function of the surface parametrization of the patch $\theta^{\alpha}$ and time $t$, and is written as

[TABLE]

Greek indices in equation (1) and from now on span the set $\{1,2\}$. At every location on the membrane surface, the parametrization $\theta^{\alpha}$ defines a natural in-plane basis $\bm{a}_{\alpha}$ given by

[TABLE]

The notation $(\,\cdot\,)_{,\alpha}$ denotes the partial derivative with respect to $\theta^{\alpha}$. At every point $\bm{x}$ on the patch $\mathcal{P}$, the set $\{\bm{a}_{1},\bm{a}_{2}\}$ forms a basis for the plane tangent to the surface at that point. The unit vector $\bm{n}$ is normal to the membrane as well as the tangent plane, and is given by

[TABLE]

The set $\{\bm{a}_{1},\bm{a}_{2},\bm{n}\}$ forms a basis of $\mathbb{R}^{3}$, and is depicted in Figure 1.

At every point $\bm{x}$, we define the dual basis to the tangent plane, $\{\bm{a}^{1},\bm{a}^{2}\}$, such that

[TABLE]

where $\delta^{\alpha}_{\beta}$ is the Kronecker delta given by $\delta^{1}_{1}=\delta^{2}_{2}=1$ and $\delta^{1}_{2}=\delta^{2}_{1}=0$. The covariant basis vectors $\bm{a}_{\alpha}$ and contravariant basis vectors $\bm{a}^{\alpha}$ are related through the metric tensor $a_{\alpha\beta}$ and contravariant metric tensor $a^{\alpha\beta}$, which are defined as

[TABLE]

and

[TABLE]

The metric tensor and contravariant metric tensor describe distances between points on the membrane surface. The covariant and contravariant basis vectors are related by

[TABLE]

where in equations (7)–(8) and from now on, indices repeated in a subscript and superscript are summed over as per the Einstein summation convention.

Any general vector $\bm{h}$ can be decomposed in the $\{\bm{a}_{1},\bm{a}_{2},\bm{n}\}$ and $\{\bm{a}^{1},\bm{a}^{2},\bm{n}\}$ bases as

[TABLE]

where $h^{\alpha}$ and $h_{\alpha}$ are contravariant and covariant components, respectively, and are related by

[TABLE]

In general, the metric tensor $a_{\alpha\beta}$ and contravariant metric tensor $a^{\alpha\beta}$ may be used to raise and lower the indices of vector and tensor components. For a general tensor $\sigma^{\alpha\beta}$, indices are raised and lowered according to

[TABLE]

For a symmetric tensor ${s^{\alpha}}_{\beta}$ with one raised and one lowered index, the order of the indices is not important and the tensor may be written as ${s^{\alpha}}_{\beta}$, ${s_{\beta}}^{\alpha}$, or $s^{\alpha}_{\beta}$, as all forms are equivalent.

When characterizing a membrane patch, it is useful to define a new basis at the membrane patch boundary $\partial\mathcal{P}$. Consider the tangent plane at a point $\bm{x}_{\mathrm{b}}$ on the membrane boundary, with the membrane normal vector $\bm{n}$. We define in-plane orthonormal basis vectors $\bm{\tau}$ and $\bm{\nu}$, where $\bm{\tau}$ is tangent to the boundary while $\bm{\nu}$ is orthogonal to $\bm{\tau}$. If the membrane boundary $\partial\mathcal{P}$ is parameterized by its arc length $\ell$, then the in-plane unit tangent $\bm{\tau}$ and in-plane unit normal $\bm{\nu}$ are defined as

[TABLE]

and

[TABLE]

The basis vectors $\bm{\nu}$ and $\bm{\tau}$ may be expressed in the covariant and contravariant bases as

[TABLE]

and

[TABLE]

Similarly, the basis vectors $\bm{a}_{\alpha}$ may be expressed in terms of the basis vectors $\bm{\nu}$ and $\bm{\tau}$ as

[TABLE]

The orthonormal basis $\{\bm{\nu},\bm{\tau},\bm{n}\}$ at a position $\bm{x}_{\mathrm{b}}$ on the membrane boundary $\partial\mathcal{P}$ is depicted in Figure 1.

The surface identity tensor $\bm{i}$ in the tangent plane and the identity tensor $\bm{1}$ in $\mathbb{R}^{3}$ are given by

[TABLE]

where $\otimes$ denotes the dyadic or outer product between any two vectors. The curvature tensor $b_{\alpha\beta}$ is given by

[TABLE]

and describes the shape of the membrane due to its embedding in $\mathbb{R}^{3}$. Given the contravariant metric and curvature tensors $a^{\alpha\beta}$ and $b_{\alpha\beta}$, the mean curvature $H$ and the Gaussian curvature $K$ can be calculated as

[TABLE]

and

[TABLE]

where the permutation tensor $\varepsilon^{\alpha\beta}$ is given by $\varepsilon^{12}=-\varepsilon^{21}=1/\sqrt{\det(a_{\alpha\beta})}$ and $\varepsilon^{11}=\varepsilon^{22}=0$. The Gaussian curvature may also be written as $K=\det(b_{\alpha\beta})/\det(a_{\alpha\beta})$. The cofactor of curvature $\bar{b}^{\alpha\beta}$ is defined as

[TABLE]

where $b^{\alpha\beta}=a^{\alpha\lambda}a^{\beta\mu}b_{\lambda\mu}$ is the contravariant form of the curvature tensor.

In general, the partial derivative of the covariant or contravariant components of a vector are not guaranteed to be invariant quantities. The covariant derivative, denoted $(\,\cdot\,)_{;\alpha}$, produces an invariant quantity when acting on vector components [69]. To define the covariant derivative, we introduce the Christoffel symbols of the second kind, denoted $\Gamma^{\alpha}_{\lambda\mu}$ and given by

[TABLE]

The covariant derivatives of the contravariant vector components $v^{\lambda}$ and covariant vector components $v_{\lambda}$ are defined as

[TABLE]

where $v^{\lambda}_{;\alpha}$ and $v_{\lambda;\alpha}$ both transform as tensors. To take the covariant derivative of second or higher order tensors, a more complicated formula is required and may be found in [69]. The covariant derivative of the metric tensor as well as the cofactor of curvature is zero. The covariant derivative of a scalar quantity is equal to its partial derivative. It is also useful to note $\bm{v}_{;\alpha}=\bm{v}_{,\alpha}$. The Gauss and Weingarten equations are

[TABLE]

respectively, and provide the covariant derivatives of the basis vectors $\bm{a}_{\alpha}$ and $\bm{n}$.

To model the kinematics of a membrane patch $\mathcal{P}$, we track the patch over time. At a reference time $t_{0}$, we define a reference patch $\mathcal{P}_{0}$. The area of the reference patch $A$ may then be compared to the area of the current patch $a$ at a later time $t$. For an infinitesimal patch $\mathrm{d}a$, the Jacobian determinant $J$ describes the areal dilation or contraction of the membrane and is defined by

[TABLE]

The Jacobian determinant may be used to convect integrals over the current membrane patch $\mathcal{P}$ to integrals over the reference patch $\mathcal{P}_{0}$, as for a scalar function $f$ we can write

[TABLE]

and the same can be written for vector- or tensor-valued functions. The details of the mapping between current and reference membrane configurations, and different coordinate parametrizations, are provided in Appendix A.1 and a detailed description can also be found in [56].

To track how quantities change over time, we define the material derivative $\mathrm{d}/\mathrm{d}t$ according to

[TABLE]

Here $(\,\cdot\,)_{,t}$ denotes the partial derivative with respect to time where $\theta^{\alpha}$ are fixed and $v^{\alpha}$ are the in-plane components of the velocity vector $\bm{v}$, which may be written as

[TABLE]

In equation (33), we have used the shorthand notation $\dot{\bm{x}}$ to express $\mathrm{d}\bm{x}/\mathrm{d}t$, and this notation will be used throughout. The scalar $v$ is the normal component of the membrane velocity given by

[TABLE]

Applying the material derivative (32) to the basis vectors $\bm{a}_{\alpha}$ is nontrivial, and requires convecting quantities to the reference patch $\mathcal{P}_{0}$. The material derivatives of the in-plane covariant basis vectors are calculated in Appendix A.2, as well as in [60], to be

[TABLE]

where the quantities $w_{\alpha}$, $w^{\alpha}$, ${w_{\alpha}}^{\beta}$, and $w_{\alpha\beta}$ are defined for notational simplicity and are given by

[TABLE]

By applying the material derivative (32) to the unit normal $\bm{n}$ (3) and using the relation $\bm{n}\cdot\dot{\bm{n}}=0$, we obtain

[TABLE]

The acceleration $\dot{\bm{v}}$ is the material derivative of the velocity and is calculated as

[TABLE]

The material derivatives of the metric tensor $a_{\alpha\beta}$ and curvature tensor $b_{\alpha\beta}$ are found to be

[TABLE]

where ${w_{\beta}}^{\lambda}$, $w_{\alpha}$, and $w_{\alpha\beta}$ are given by equations (36)–(38). Finally, the time derivative of the Jacobian determinant is found in [60] as

[TABLE]

In three-dimensional Cartesian systems, $\dot{J}/J=\mathrm{div}\,\bm{v}$. Comparing the Cartesian result with the right hand side of equation (44), we see that the two-dimensional analog for $\mathrm{div}\,\bm{v}$ is $v^{\alpha}_{;\alpha}-2vH$.

3 Intra-membrane Flow and Bending

In this section, we develop a comprehensive model of a single-component lipid membrane which behaves like a viscous fluid in-plane and an elastic shell in response to out-of-plane bending. We use the balance law framework for single-component membranes previously proposed in several works [51, 60, 59, 65, 70] and in later sections extend it to model multi-component membranes, phase transitions, and the binding of proteins to the membrane surface.

We begin by determining local forms of the balances of mass, linear momentum, and angular momentum. We then go on to determine the form of the membrane stresses through a systematic thermodynamic treatment. To this end, we develop local forms of the first and second laws of thermodynamics and a local entropy balance. By postulating the dependence of membrane energetics on the appropriate fundamental thermodynamic variables, we determine constitutive equations for the in-plane and out-of-plane stresses. We then use these stresses to provide the equations of motion as well as possible boundary and initial conditions for the membrane, and conclude by briefly discussing how membrane dynamics can be coupled to the surrounding bulk fluid.

3.1 Balance Laws

Our general procedure is to start with a global form of the balance law for an arbitrary membrane patch $\mathcal{P}$, convert each term to an integral over the membrane surface, and invoke the arbitrariness of $\mathcal{P}$ to determine the local form of the balance law. To convert terms in the global balance laws to integrals over the membrane patch, we will need tools to bring total time derivatives inside the integral and convert integrals over the patch boundary to integrals over the membrane surface.

For a scalar-, vector-, or tensor-valued function $f$ defined on the membrane patch $\mathcal{P}$, the Reynolds transport theorem describes how time derivatives commute with integrals over the membrane surface. As described in [60], the Reynolds transport theorem is given by

[TABLE]

Now consider a vector- or tensor-valued function $\bm{f}$, which may be expressed as $\bm{f}=f^{\alpha}\bm{a}_{\alpha}+f\bm{n}$. The surface divergence theorem describes how an integral of $\bm{f}\cdot\bm{\nu}=f^{\alpha}\nu_{\alpha}$ over the membrane boundary $\partial\mathcal{P}$, where $\bm{\nu}$ is the boundary normal in the tangent plane defined in equation (15), may be converted to a surface integral over the membrane patch $\mathcal{P}$. To this end, the surface divergence theorem states

[TABLE]

where $\mathrm{d}s$ is an infinitesimal line element on the membrane boundary.

3.1.1 Mass Balance

Consider a membrane patch $\mathcal{P}$ with a mass per unit area denoted as $\rho(\theta^{\alpha},t)$. The total mass of the membrane patch is conserved, and the global form of the conservation of mass can be written as

[TABLE]

Applying the Reynolds transport theorem (45) to the global mass balance (47) brings the time derivative inside the integral, and we obtain

[TABLE]

Since the membrane patch $\mathcal{P}$ is arbitrary, the local form of the conservation of mass is given by

[TABLE]

As the total mass of the membrane patch is conserved, the mass at any time $t$ is equal to the mass at time $t_{0}$, i.e.,

[TABLE]

where $\rho_{0}=\rho(\theta^{\alpha},t_{0})$ is the areal mass density of the reference patch. Using equation (31) yields

[TABLE]

As the reference patch $\mathcal{P}_{0}$ is arbitrary, the Jacobian determinant $J$ is given by

[TABLE]

in addition to the form provided in equation (30).

Substituting $f=\rho u$ into the Reynolds transport theorem (45), where $u$ is an arbitrary quantity per unit mass, and using equation (49), we obtain

[TABLE]

Equation (53) is a modified Reynolds transport theorem and is useful in simplifying balance laws where quantities are defined per unit mass.

3.1.2 Linear Momentum Balance

It is well-known from Newtonian and continuum mechanics that the rate of change of momentum of a body is equal to the sum of the external forces acting on it. Lipid membranes may be acted on by two types of forces: body forces on the membrane patch $\mathcal{P}$ and tractions on the membrane boundary $\partial\mathcal{P}$. On the membrane patch $\mathcal{P}$, the body force per unit mass is denoted by $\bm{b}(\theta^{\alpha},t)$. At a point $\bm{x}_{\mathrm{b}}$ on the membrane boundary $\partial\mathcal{P}$ with in-plane unit normal $\bm{\nu}$, the boundary traction is the force per unit length acting on the membrane boundary and is denoted by $\bm{T}(\bm{x}_{\mathrm{b}},t;\bm{\nu})$. The global form of the balance of linear momentum for any membrane patch $\mathcal{P}$ is given by

[TABLE]

where the left hand side is the time derivative of the total linear momentum of the membrane patch and the right hand side is the sum of the external forces.

For three-dimensional systems in Cartesian coordinates, one may use Cauchy’s tetrahedron arguments to decompose the boundary tractions and define the Cauchy stress tensor, which specifies the total state of stress at any location [71]. Naghdi [36] performed an analogous procedure on a curvilinear triangle on an arbitrary surface to show boundary tractions may be expressed as a linear combination of the stress vectors $\bm{T}^{\alpha}$ according to

[TABLE]

The stress vectors $\bm{T}^{\alpha}$ describe the tractions along curves of constant $\theta^{\alpha}$ and are independent of the in-plane boundary unit normal $\bm{\nu}$. Substituting the traction decomposition (55) into the global linear momentum balance (54), applying the surface divergence theorem (46) on the traction term, and applying the Reynolds transport theorem (53) on the left hand side, we obtain

[TABLE]

Since $\mathcal{P}$ is arbitrary, equation (56) yields the local form of the linear momentum balance as

[TABLE]

To recast the traction decomposition (55) into a more familiar form involving the Cauchy stress tensor, we express the stress vectors $\bm{T}^{\alpha}$ in the $\{\bm{a}_{\alpha},\bm{n}\}$ basis without loss of generality as

[TABLE]

where $N^{\alpha\beta}$ and $S^{\alpha}$ are the components of the stress vector $\bm{T}^{\alpha}$ in the $\{\bm{a}_{\alpha},\bm{n}\}$ basis [60, 59]. Substituting the form of the stress vectors $\bm{T}^{\alpha}$ (58) into the traction decomposition (55) allows us to write

[TABLE]

where $\bm{\sigma}$ is the Cauchy stress tensor given by

[TABLE]

Consequently, $N^{\alpha\beta}$ and $S^{\alpha}$ can also be interpreted as the in-plane and out-of-plane components of the stress tensor $\bm{\sigma}$. In specifying $N^{\alpha\beta}$ and $S^{\alpha}$, we will have completely determined the total state of stress at any location on the membrane. The in-plane tension $\hat{\sigma}$ is one-half the trace of the stress tensor (60), and as found in [62] is given by

[TABLE]

The equation for the in-plane tension $\hat{\sigma}$ (61) reinforces the notion that $N^{\alpha\beta}$ describes in-plane stresses and $S^{\alpha}$ describes out-of-plane stresses, as only $N^{\alpha\beta}$ enters equation (61).

When solving for the strong forms of the dynamical equations of motion, we will need to consider the linear momentum balance (57) in component form. In what follows, we decompose the equations of motion in the directions normal and tangential to the surface. To this end, the body force $\rho\bm{b}$ may be expressed as

[TABLE]

where $p$ is the pressure normal to the membrane and $b^{\alpha}$ are the in-plane contravariant components of the body force per unit mass. To express $\bm{T}^{\alpha}_{;\alpha}$ in component form, we apply the Gauss (28) and Weingarten (29) equations to the stress vector decomposition (58) and obtain

[TABLE]

Substituting the body force decomposition (62), divergence of the stress vectors (63), and acceleration (41) into the local form of the linear momentum balance (57), we find the tangential and normal momentum equations are given, respectively, by

[TABLE]

The normal component of the linear momentum balance (65) is usually referred to as the shape equation [72, 45].

Although we do not yet know the form of the stresses $N^{\alpha\beta}$ and $S^{\alpha}$, from equations (64) and (65) we already see coupling between in-plane and out-of-plane membrane behavior. The in-plane stresses $N^{\alpha\beta}$ and the out-of-plane shear $S^{\alpha}$ appear in the in-plane equations (64) and the shape equation (65). In general, we expect in-plane flow to influence out-of-plane bending and vice versa.

The three components of the linear momentum balance (64)–(65) and the mass balance (49) allow us to solve for the four fundamental unknowns: the density $\rho$ and the velocity components $v$ and $v^{\alpha}$. To solve the equations of motion, however, we must first determine the forms of $N^{\alpha\beta}$ and $S^{\alpha}$. We will now systematically determine the form of the in-plane and shear stresses before returning to the equations of motion.

3.1.3 Angular Momentum Balance

In this section, we analyze the balance of angular momentum of the membrane. The rate of change of the total angular momentum of the membrane is equal to the sum of the external torques acting on the membrane patch. In addition to the torques arising from body forces and boundary tractions, the membrane is able to sustain director tractions on its boundary. These director tractions give rise to additional external moments on the membrane boundary which will twist the edges of the membrane patch, as depicted in Figure 2. We will show such moments are necessary to sustain the shear stresses $S^{\alpha}$ introduced in the linear momentum balance. Moreover, in the absence of director tractions the in-plane stresses are shown to be symmetric.

In general, we specify a director field $\bm{d}$ on the membrane patch $\mathcal{P}$ to account for the finite thickness of the lipid membrane [51, 36]. The director $\bm{d}(\theta^{\alpha},t)$ is a unit vector describing the orientation of the phospholipids—when the director $\bm{d}$ does not coincide with the normal to the surface $\bm{n}$, the phospholipids are tilted relative to the normal. At a point $\bm{x}_{\mathrm{b}}$ on the patch boundary $\partial\mathcal{P}$, the director traction $\bm{M}(\bm{x}_{\mathrm{b}},t)$ describes the equal and opposite forces acting on the director $\bm{d}(\bm{x}_{\mathrm{b}},t)$. As the director $\bm{d}$ is dimensionless, the moment per length $\bm{m}$ at the patch boundary is given by the cross product $\bm{d}\times\bm{M}$.

To properly account for the director field, it is necessary to include the director velocity $\dot{\bm{d}}$ in an additional balance law for the director momentum as described by Naghdi and Green [36, 73, 74]. Including the director field would enable us to examine smaller length scale phenomena, for example transmembrane proteins causing phospholipids to tilt and inducing local inhomogeneities in the directors. In this work, however, we choose not to study such phenomena and treat the membrane as a sheet of zero thickness. In doing so, the director is forced to be equal to the normal at every point $\bm{x}$ and is prescribed to be

[TABLE]

Note that equation (66) is equivalent to the Kirchhoff-Love assumption [51]. With this simplification, the moment per unit length of the patch boundary, $\bm{m}$, is given by

[TABLE]

Given equation (67), the global form of the angular momentum balance can be written as

[TABLE]

where $\rho\bm{x}\times\bm{v}$ denotes the angular momentum density at the point $\bm{x}$, and $\rho\bm{x}\times\bm{b}$ and $\bm{x}\times\bm{T}$ denote the torque densities due to body forces and tractions, respectively.

While the director traction $\bm{M}$ may in general have normal and tangential components, the component in the normal direction has no effect on the resulting couple $\bm{m}$ due to equation (67). Thus we restrict $\bm{M}$ to be in the plane of the membrane. Once again using elementary curvilinear triangle arguments described by Naghdi [36], the director traction $\bm{M}$ may be written as

[TABLE]

The couple-stress vectors $\bm{M}^{\alpha}$ in equation (69) must be in the plane of the membrane due to our imposed restriction, and may be written without loss of generality as

[TABLE]

Substituting the couple-stress decomposition (70) into the director traction decomposition (69) allows us to write

[TABLE]

where $\bm{\mu}$ is the couple-stress tensor given by

[TABLE]

Because we require director tractions to not lie in the normal direction, the couple-stress tensor $\bm{\mu}$ in equation (72) does not have any $\bm{a}_{\alpha}\otimes\bm{n}$ component.

Returning to the global form of the angular momentum balance (68), we substitute the director traction decomposition (69) and stress vector decomposition (55) to obtain

[TABLE]

Using the Reynolds transport theorem (53) and the surface divergence theorem (46), equation (73) simplifies to

[TABLE]

Since the membrane patch $\mathcal{P}$ is arbitrary, the local form of the angular momentum balance can be obtained as

[TABLE]

where we have distributed the covariant derivatives and used the Gauss (28) and Weingarten (29) equations.

It is useful to know what constraints the local form of the angular momentum balance (75) imposes in addition to what was known from the linear momentum balance (57). Taking the cross product of $\bm{x}$ with the local linear momentum balance (57) and subtracting it from the local angular momentum balance (75) gives

[TABLE]

Substituting the couple-stress decomposition (70) and traction decomposition (55) into equation (76), we obtain

[TABLE]

Equation (77) indicates the following conditions must be true in order for both the linear momentum balance and the angular momentum balance to be locally satisfied:

[TABLE]

In equation (78), the tensor $\sigma^{\alpha\beta}$ describes the components of in-plane tractions due to stretching and viscous flow only, i.e., $\sigma^{\alpha\beta}$ does not include contributions from moments. This is to say the combination of angular and linear momentum balances impose restrictions between the in-plane stress components $N^{\alpha\beta}$, out-of-plane shear stress components $S^{\alpha}$, and the components of the couple-stress tensor $-M^{\alpha\beta}$. As the boundary moment per length $\bm{m}$ is related to the components of $M^{\alpha\beta}$, equation (79) indicates the relationship between out-of-plane shear stresses and boundary moments. If boundary moments had not been included, there would consequently be no shear stresses at any point on the membrane surface.

Finally, it will be useful to express the boundary moment per length $\bm{m}$ in terms of the in-plane boundary tangent $\bm{\tau}$ and boundary normal $\bm{\nu}$ as

[TABLE]

Using the identity $\bm{a}_{\beta}\times\bm{n}=\tau_{\beta}\bm{\nu}-\nu_{\beta}\bm{\tau}$, which can be derived from the decomposition of the in-plane unit normal $\bm{\nu}$ (16) and in-plane unit tangent $\bm{\tau}$ (17), and substituting the director traction decomposition (69) and couple-stress decomposition (70) into the equation for the moment per length $\bm{m}$ (67), we find the components of the boundary moment per length $\bm{m}$ to be

[TABLE]

At this stage, all previous works using either the balance law formulation [60, 65] or variational methods [53] propose constitutive forms of the in-plane viscous stresses and in-plane velocity gradients to model the irreversible processes of fluid flow. These are then used to determine the equations of motion. In our work, we will naturally find the constitutive form of the in-plane viscous stresses by evaluating the entropy production and proposing relationships between the thermodynamic forces and fluxes in the linear irreversible regime. This framework based on entropy production is naturally extendable to multi-component systems and systems with chemical reactions. In what follows, we proceed to develop such a framework.

3.1.4 Mechanical Power Balance

While a mechanical power balance does not impose any new constraints on the membrane patch, it expresses the relationship between the kinetic energy, internal forces, and external forces, which is useful for the entropy production derivations in subsequent sections. We begin by taking the dot product of the local momentum balance (57) with the velocity $\bm{v}$ and integrating over the membrane patch $\mathcal{P}$ to obtain

[TABLE]

The left hand side of equation (83) is the material derivative of the total kinetic energy, as an application of the Reynolds transport theorem (53) shows

[TABLE]

The first term on the right hand side of equation (83) may be expanded as

[TABLE]

where the second equality is obtained by invoking the surface divergence theorem (46) and the third equality from the boundary traction decomposition (55). By expanding the integrand of the last term in equation (85), we find

[TABLE]

where the first equality is obtained with the relation for $\bm{v}_{,\alpha}$ (35) and $\bm{T}^{\alpha}$ (58), and the third equality by substituting the results of the angular momentum balance (78)–(79). Using the symmetry of $\sigma^{\alpha\beta}$ found in equation (78) and the relation for $\dot{a}_{\alpha\beta}$ (42), the first term in the final equality of equation (86) may be written as $\sigma^{\alpha\beta}w_{\alpha\beta}=\tfrac{1}{2}\sigma^{\alpha\beta}(w_{\alpha\beta}+w_{\beta\alpha})=\tfrac{1}{2}\sigma^{\alpha\beta}\dot{a}_{\alpha\beta}$. Using the product rule on the last term in the final equality of equation (86) gives $M^{\beta\alpha}_{;\beta}w_{\alpha}=(M^{\beta\alpha}w_{\alpha})_{;\beta}-M^{\beta\alpha}w_{\alpha;\beta}$. Using these simplifications, equation (86) may be rewritten as

[TABLE]

The second term on the right hand side of equation (87) is $M^{\mu\alpha}\dot{b}_{\mu\alpha}$, given the relation for $\dot{b}_{\alpha\beta}$ in equation (43). We rewrite the last term in equation (87) as

[TABLE]

With the above simplifications, we find equation (87) reduces to

[TABLE]

Using equation (89), equation (85) can be written as

[TABLE]

where the second equality is obtained by using the surface divergence theorem (46). Substituting equations (90) and (84) into equation (83), we find the total mechanical power balance is given by

[TABLE]

The left hand side of equation (91) contains the material derivative of the kinetic energy (84) and a term describing the internal changes involving the shape and stresses of the membrane, which describe the membrane’s internal power. The terms on the right hand side of the mechanical power balance (91) describe the power due to external forces and moments acting on the membrane.

3.2 Thermodynamics

In this section, we develop the thermodynamic framework necessary to understand the effects of bending and intra-membrane viscous flow on the membrane patch. We develop local forms of the first law of thermodynamics and entropy balance, and we introduce the second law of thermodynamics. We follow the procedure described by de Groot & Mazur [68] to understand the internal entropy production, albeit with one difference. While de Groot & Mazur [68] begin with the local equilibrium assumption and the Gibbs equation, it is technically difficult to write the Gibbs equation for a system which depends on tensorial quantities. In this work, we follow the approach demonstrated in [75] and begin by choosing the appropriate form of the Helmholtz free energy. Following this framework, one can derive an effective Gibbs equation after the analysis is complete.

3.2.1 First Law—Energy Balance

According to the first law of thermodynamics, the total energy of the membrane patch changes due to work being done on the membrane or heat flowing into the membrane. The mechanical power balance (91) describes the rate of work being done on the membrane due to external tractions, moments, and forces. Furthermore, heat may enter or exit the membrane patch in one of two ways: by flowing from the surrounding medium into the membrane along the normal direction $\bm{n}$, or by flowing in the plane of the membrane across the membrane patch boundary. We denote the heat source per unit mass as $r(\theta^{\alpha},t)$, which accounts for the heat flow from the bulk, and the in-plane heat flux as $\bm{J}_{\mathrm{q}}=J_{\mathrm{q}}^{\,\alpha}\,\bm{a}_{\alpha}$. By convention, the heat flux $\bm{J}_{\mathrm{q}}$ is positive when heat flows out of the system across the patch boundary. Defining $e(\theta^{\alpha},t)$ to be the total energy per unit mass of the membrane, the global form of the first law of thermodynamics can be written as

[TABLE]

The total energy per unit mass $e$ consists of the internal energy per unit mass $u$ and the kinetic energy per unit mass $\tfrac{1}{2}\bm{v}\cdot\bm{v}$, and is given by

[TABLE]

Using the Reynolds transport theorem (53) and substituting the expression for the total energy per mass $e(\theta^{\alpha},t)$ (93) into equation (92), we obtain

[TABLE]

Equation (94) shares several terms with the mechanical power balance (91), and by subtracting the two equations, the balance of internal energy can be obtained as

[TABLE]

where the second equality is obtained by using the surface divergence theorem (46). Since the membrane patch $\mathcal{P}$ is arbitrary, the local form of the internal energy balance is given by

[TABLE]

The first two terms on the right hand side of equation (96) describe the heat flow into the system, and the last two terms describe the energy change due to work being done on the system.

3.2.2 Entropy Balance & Second Law

The total entropy of a membrane patch $\mathcal{P}$ may change in three ways: entropy may flow into or out of the patch across the membrane boundary, entropy may be absorbed or emitted from the membrane body as a supply, or entropy may be produced internally within the membrane patch. The local quantities corresponding to such changes are the in-plane entropy flux $\bm{J}_{\mathrm{s}}=J_{\mathrm{s}}^{\,\alpha}\bm{a}_{\alpha}$, the rate of external entropy supply per unit mass $\eta_{\mathrm{e}}(\theta^{\alpha},t)$, and the rate of internal entropy production per unit mass $\eta_{\mathrm{i}}(\theta^{\alpha},t)$, respectively. For the total entropy per unit mass $s(\theta^{\alpha},t)$, the global form of the entropy balance is given by

[TABLE]

Applying the Reynolds transport theorem (53) and the surface divergence theorem (46) reduces equation (97) to

[TABLE]

Again, due to the arbitrariness of the membrane patch $\mathcal{P}$, the local form of the entropy balance is given by

[TABLE]

At this point, it is useful to consider the nature of the entropy flux, external entropy supply, and internal entropy production. We define the in-plane entropy flux $\bm{J}_{\mathrm{s}}$ and the external entropy supply per unit mass $\eta_{\mathrm{e}}$ to describe the redistribution of entropy that has already been created. These terms may be positive or negative. We now introduce the second law of thermodynamics by requiring the internal entropy production to be non-negative at every point in the membrane. The second law of thermodynamics is given by

[TABLE]

The internal entropy production (100) is zero only for reversible processes.

3.2.3 Choice of Thermodynamic Potential

The natural thermodynamic potential for the membrane patch is the Helmholtz free energy [59]. The Helmholtz free energy per unit mass, $\psi$, is given by

[TABLE]

where $T(\theta^{\alpha},t)$ is the local temperature of the membrane patch. Taking the material derivative of equation (101), solving for $\dot{s}$, and substituting into the local entropy balance (99), we obtain

[TABLE]

Substituting the local form of the first law of thermodynamics (96) into equation (102) yields the total rate of change of entropy, given by

[TABLE]

Equation (103) will allow us to determine which terms contribute to the internal entropy production, understand fundamental relationships between the stresses, moments, and energetics of the membrane patch, and finally develop constitutive relations between the stresses, moments, and associated kinematic quantities.

3.3 Constitutive Relations

In this section, we choose the fundamental thermodynamic variables for our membrane patch. With this constitutive assumption, we determine the contributions to the entropy flux, external entropy supply, and internal entropy production. We then apply linear irreversible thermodynamics to relate generalized thermodynamic forces to their corresponding fluxes. In doing so, we naturally determine the viscous dissipation due to in-plane fluid flow as well as the dependence of the stresses and moments on the Helmholtz free energy density.

3.3.1 General Thermodynamic Variables

Lipid bilayers have in-plane dissipative flow and out-of-plane elastic bending. The Helmholtz free energy per unit mass $\psi$, as a thermodynamic state function, captures the elastic behavior of lipid membranes. The general thermodynamic variables that the Helmholtz free energy density of a two-dimensional elastic sheet depends on are the metric tensor $a_{\alpha\beta}$, curvature tensor $b_{\alpha\beta}$, and temperature $T$ [36, 59]. The simplest form of the Helmholtz free energy density $\psi$ that captures this behavior is given by

[TABLE]

In this work, we assume the membrane does not thermally expand or chemically swell, so the metric and curvature tensors capture only elastic behavior.

Because the metric and curvature tensors are symmetric, the material derivative of the Helmholtz free energy density $\psi$ is given by

[TABLE]

Substituting equation (105) into the local entropy balance (103), we obtain

[TABLE]

At this stage, we assume the system is locally at equilibrium, and therefore define the entropy as

[TABLE]

where the partial derivative is taken at constant $a_{\alpha\beta}$ and $b_{\alpha\beta}$. Rewriting the heat flux and using equation (107) reduces equation (106) to

[TABLE]

From dimensional arguments, only gradients on the right hand side may contribute to the in-plane entropy flux components $J_{\mathrm{s}}^{\,\alpha}$, which are given by

[TABLE]

The external entropy supply per unit area $\rho\eta_{\mathrm{e}}$ captures entropy being absorbed or emitted across the membrane body. The only term on the right hand side which describes such a change is the heat source $r$. Therefore, the external entropy per unit area $\rho\eta_{\mathrm{e}}$ is given by

[TABLE]

In equations (109) and (110), we obtain the familiar result that heat flow into or out of the system is associated with an entropy change.

As we have determined the terms on the right hand side of equation (108) that contribute to the entropy flux and external entropy, the remaining terms contribute to the internal entropy production. To this end, the rate of internal entropy production per unit area $\rho\eta_{\mathrm{i}}$ is given by

[TABLE]

The terms on the right hand side of equation (111) are a product of a thermodynamic force, which may be imposed on the system, and a thermodynamic flux. Denoting the thermodynamic force as $X_{k}$ and the corresponding flux as $J^{k}$, equation (111) may be generally written as

[TABLE]

In equation (112), the indices $k$ are used as a label, as $X_{k}$ and $J^{k}$ may be scalars, vectors, or tensors. As described by Prigogine [67] and de Groot & Mazur [68], we assume in the linear irreversible regime, i.e., near equilibrium, there is a linear relationship between $X_{k}$ and $J^{k}$ given by

[TABLE]

where $L^{ik}$ are the phenomenological coefficients.

In the internal entropy production (111), there are three thermodynamic forces: the in-plane temperature gradient $T_{,\alpha}$ and the material derivatives of the metric and curvature tensor, $\dot{a}_{\alpha\beta}$ and $\dot{b}_{\alpha\beta}$, respectively. We invoke the Curie principle [76], as done by Prigogine [67] and de Groot & Mazur [68], and propose that the phenomenological coefficients between quantities with different tensorial order must be zero. Therefore, the heat flux $J_{\mathrm{q}}^{\,\alpha}$ is independent of the tensorial forces $\dot{a}_{\alpha\beta}$ and $\dot{b}_{\alpha\beta}$. Similarly, the stresses and moments are independent of the temperature gradients. In the spirit of equation (113), the phenomenological relation for the heat flux is then given by

[TABLE]

where the tensor $\kappa^{\alpha\beta}$ is the thermal conductivity tensor. As the fluid is thermally isotropic in-plane, $\kappa^{\alpha\beta}=\kappa\,a^{\alpha\beta}$, where the constant $\kappa$ is the scalar thermal conductivity. In this case, equation (114) reduces to

[TABLE]

where we use the shorthand $T^{,\alpha}$ to denote $T_{,\beta}\,a^{\alpha\beta}$. We note that in the case of lipid bilayers, there are usually no temperature gradients and equation (115) does not play a major role in describing the relevant irreversible processes.

To obtain the remaining phenomenological coefficients associated with the other irreversible processes in equation (111), we define the thermodynamic fluxes

[TABLE]

for notational convenience. In the linear irreversible regime, the phenomenological relations relating $\pi^{\alpha\beta}$ and $\omega^{\alpha\beta}$ to $\dot{a}_{\alpha\beta}$ and $\dot{b}_{\alpha\beta}$ can be generally written as

[TABLE]

where the fourth-order contravariant tensors $R^{\alpha\beta\gamma\mu}$, $S^{\alpha\beta\gamma\mu}$, $T^{\alpha\beta\gamma\mu}$, and $U^{\alpha\beta\gamma\mu}$ are general fourth-order phenomenological viscous coefficients. The tensors $S^{\alpha\beta\gamma\mu}$ and $T^{\alpha\beta\gamma\mu}$ describe interference between the two irreversible processes driven by $\dot{a}_{\alpha\beta}$ and $\dot{b}_{\alpha\beta}$, and we assume them to be zero for the case of lipid bilayers. The phenomenological relations (118)–(119) then reduce to

[TABLE]

Given the form of the internal entropy production (111), equations (120) and (121) indicate $\pi^{\alpha\beta}$ captures the dissipation due to in-plane flow and $\omega^{\alpha\beta}$ captures the dissipation due to out-of-plane bending. In general, $U^{\alpha\beta\gamma\mu}$ need not be equal to zero and bending can provide another way by which the membrane dissipates energy. However, we assume out-of-plane-bending is not a dissipative process and so $U^{\alpha\beta\gamma\mu}=0$. Consequently, $\omega^{\alpha\beta}=0$, which leads to the constitutive relation for the couple-stress tensor $M^{\alpha\beta}$ being given by

[TABLE]

Because of the in-plane viscous nature of the lipid bilayer, $\pi^{\alpha\beta}$ is nonzero. Lipid membranes are isotropic in-plane, indicating $R^{\alpha\beta\gamma\mu}$ is an isotropic tensor. A fourth-order tensor in general curvilinear coordinates is isotropic when it is invariant to all unimodular transformations of the coordinate system, represented by the tensor ${\Lambda^{\nu}}_{\sigma}$, such that

[TABLE]

For ${\Lambda^{\nu}}_{\sigma}$ to represent a unimodular coordinate transformation, it must satisfy ${\Lambda^{\nu}}_{\sigma}{\Lambda_{\lambda}}^{\sigma}=\delta^{\nu}_{\lambda}$ and ${\Lambda^{\nu}}_{\sigma}{\Lambda_{\nu}}^{\lambda}=\delta^{\lambda}_{\sigma}$, where ${(\Lambda^{-1})^{\nu}}_{\sigma}={\Lambda_{\nu}}^{\sigma}$ [69]. In what follows, we choose forms of ${\Lambda^{\nu}}_{\sigma}$ satisfying these requirements to determine constraints on the form of $R^{\alpha\beta\gamma\mu}$.

First, consider a rotation of the coordinate axes by $\pi/2$ radians about the direction of the normal vector $\bm{n}$. The transformation tensor ${\Lambda^{\nu}}_{\sigma}$ corresponding to this rotation is given by

[TABLE]

Applying equation (124) to the definition of an isotropic tensor (123), we obtain

[TABLE]

reducing the initial 16 variables in $R^{\alpha\beta\gamma\mu}$ to eight.

Next, consider a transformation where we exchange $\bm{a}_{1}$ and $\bm{a}_{2}$. The transformation tensor ${\Lambda^{\nu}}_{\sigma}$ for this operation is given by

[TABLE]

Applying equation (126) to the definition of an isotropic tensor (123) leads to

[TABLE]

reducing the remaining eight variables in $R^{\alpha\beta\gamma\mu}$ to four.

The final transformation we consider is a rotation of the coordinate axes by $\pi/4$ radians about the direction of the normal vector $\bm{n}$. In this case, the transformation tensor ${\Lambda^{\nu}}_{\sigma}$ is given by

[TABLE]

Applying equation (128) to equation (123) yields

[TABLE]

thereby reducing the four remaining degrees of freedom in $R^{\alpha\beta\gamma\mu}$ to three.

Given the independent variables of $R^{\alpha\beta\gamma\mu}$, we now determine its functional form. Due to the linear independence of $R^{1122}$, $R^{1212}$, and $R^{1221}$, we may specify each term arbitrarily. Consider the case where $R^{1122}\neq 0$ and $R^{1212}=R^{1221}=0$. From equations (129) and (125), $R^{1111}=R^{2222}=R^{1122}=R^{2211}$. In Cartesian coordinates, the form of $R^{\alpha\beta\gamma\mu}$ would be expressed with Kronecker delta functions [77]. In curvilinear coordinates, the contravariant form of $\delta^{\alpha}_{\beta}$ is $\delta^{\alpha}_{\mu}a^{\mu\beta}=a^{\alpha\beta}$ and so we may write $R^{\alpha\beta\gamma\mu}=\tfrac{1}{2}\lambda\,a^{\alpha\beta}a^{\gamma\mu}$, where $\lambda$ is a constant and the factor of $1/2$ is included for convenience. The next case to consider is when $R^{1212}=R^{2121}=R^{1111}=R^{2222}\neq 0$ and $R^{1122}=R^{1221}=0$. In a similar fashion to the previous case, with curvilinear coordinates we find $R^{\alpha\beta\gamma\mu}=\zeta_{1}\,a^{\alpha\gamma}a^{\beta\mu}$, where $\zeta_{1}$ is another constant. Finally, if $R^{1221}=R^{2112}\neq 0$ and $R^{1212}=R^{1122}=0$, we obtain $R^{\alpha\beta\gamma\mu}=\zeta_{2}\,a^{\alpha\mu}a^{\beta\gamma}$, where $\zeta_{2}$ is yet another constant. Due to the linear independence of $R^{1122}$, $R^{1212}$, and $R^{1221}$, the general form of $R^{\alpha\beta\gamma\mu}$ is given by the sum

[TABLE]

Equation (130) can be considered as providing the general form of a fourth-order isotropic tensor in curvilinear coordinates.

Substituting the form of $R^{\alpha\beta\gamma\mu}$ in equation (130) into the expression for $\pi^{\alpha\beta}$ (120), we obtain

[TABLE]

where in the third equality we use the symmetry of $\dot{a}_{\gamma\mu}$. Using equation (44) and defining $\zeta=\zeta_{1}+\zeta_{2}$ reduces equation (131) to

[TABLE]

With the form of $\pi^{\alpha\beta}$ in equation (132) and $M^{\alpha\beta}$ in equation (122), the internal entropy production (111) simplifies to

[TABLE]

In the absence of temperature gradients, equation (133) simplifies to

[TABLE]

where use has been made of equation (44). For the inequality in equation (134) to hold, we require $\zeta\geq 0$ and $\lambda\geq 0$. Physically, $\zeta$ describes the internal entropy production from velocity gradients and represents the in-plane shear viscosity, while $\lambda$ describes internal entropy production due to the fluid compressing or expanding and represents the in-plane bulk viscosity. The final form of the total internal entropy production in the linear irreversible regime is given by

[TABLE]

where $\kappa$, $\zeta$, and $\lambda$ are all non-negative.

In this section, we determined how the stresses and moments of the membrane are related to the Helmholtz energy density, and the results can be summarized as

[TABLE]

with $\pi^{\alpha\beta}$ given by equation (132).

At this stage, we derive the Gibbs equation for the single-component membrane system. The Gibbs equation in general relates infinitesimal changes in thermodynamic state functions, and consequently will not account for any dissipation in the system. It is therefore useful to define $(\sigma^{\alpha\beta})^{\text{el}}$ to be the reversible, elastic component of the in-plane stress $\sigma^{\alpha\beta}$ (136) given by

[TABLE]

To derive the Gibbs equation for our membrane system, we start with equation (102) and substitute the material derivative of $\psi$ (105), the local equilibrium assumption (107), the moment tensor $M^{\alpha\beta}$ (137), and the elastic component of the in-plane traction $(\sigma^{\alpha\beta})^{\text{el}}$ (140) to obtain

[TABLE]

Equation (141) relates the rates of change of thermodynamic state functions, and by multiplying both sides of the equation by $\mathrm{d}t$ we find

[TABLE]

Equation (142) is the Gibbs equation for a two-dimensional membrane surface with out-of-plane elastic bending and in-plane elastic compression and stretching.

3.3.2 Helmholtz Free Energy—Change of Variables

We have so far developed general equations of how the membrane stresses depend on a Helmholtz free energy density $\psi$, which in turn depends on the metric tensor $a_{\alpha\beta}$, the curvature tensor $b_{\alpha\beta}$, and the temperature $T$ (104). The energy density $\psi$, in being an absolute scalar field, must be invariant to Galilean transformations. For a fluid film, under Galilean invariance the Helmholtz free energy density may depend on $a_{\alpha\beta}$ and $b_{\alpha\beta}$ only through the density $\rho$, the mean curvature $H$, and the Gaussian curvature $K$, which are functions of the invariants of the metric and curvature tensors [59]. This relationship is written as

[TABLE]

Note equation (143) can also be shown using material symmetry arguments as presented in [78].

When substituting $\bar{\psi}$ into the stress and moment relations (136)–(139), we encounter terms like

[TABLE]

where subscripts including a comma denote a partial derivative—for example, $\bar{\psi}_{,\rho}=\partial\bar{\psi}/\partial\rho$. The variations in $\rho$, $H$, and $K$ due to changes in $a_{\alpha\beta}$ and $b_{\alpha\beta}$ can be easily calculated, and are summarized in Table 1. Substituting the partial derivatives given in Table 1 into equations (144)–(145), we find the stress and moment tensors (136)–(139) expressed in terms of $\bar{\psi}(\rho,H,K,T)$ to be given by

[TABLE]

Equations (146)–(149) are identical to those found in [60], obtained without using the formulation of an irreversible thermodynamic framework. When substituting equations (146)–(149) in the equations of motion, it is useful to write the viscous stresses $\pi^{\alpha\beta}$ in a different form which contains the mean and Gaussian curvatures. Substituting equation (42) into equation (132) and using the definition of the cofactor of curvature (24), we obtain

[TABLE]

where $d^{\alpha\beta}$ is the symmetric part of the in-plane velocity gradients defined as

[TABLE]

3.3.3 Helfrich Energy Density

For single-component lipid bilayers, the Helmholtz free energy density contains an energetic cost for bending and an energetic cost for areal compressions and dilations. The energetic cost of bending, called the Helfrich energy [38] and denoted $w_{\textrm{h}}$, is given by

[TABLE]

The constants $k_{\mathrm{b}}$ and $k_{\mathrm{g}}$ are the mean and Gaussian bending moduli, respectively. In equation (152), $C$ is the spontaneous curvature induced by proteins or lipids which the membrane would like to conform to. The compression energy density $w_{\textrm{c}}$ equally penalizes areal compression and dilation, and is given by

[TABLE]

where $k_{\mathrm{c}}$ is the compression modulus. As $k_{\mathrm{c}}$ tends to infinity, the membrane becomes incompressible. The factor of $1/J$ in $w_{\textrm{c}}$ is necessary as areal compressions and dilations are calculated with respect to the reference patch $\mathcal{P}_{0}$ at time $t_{0}$. To see this explicitly, consider the total compression energy $W_{\textrm{c}}$ in the reference frame given by

[TABLE]

The total Helmholtz free energy per unit area is given by

[TABLE]

where $f(T)$ is a function of the temperature such that the entropy $s$ can be calculated using equation (107).

Given the total Helmholtz energy per unit area $\rho\bar{\psi}$ (155) and the forms of the Helfrich (152) and compression (153) energies, we calculate the partial derivatives

[TABLE]

In obtaining equation (158) it is important to realize $J$ is a function of $\rho$, as shown in equation (52). Substituting the partial derivates (156)–(158) into the stress and moment tensors (146) - (149) yields

[TABLE]

with $\pi^{\alpha\beta}$ given by equation (150). Equations (159)–(162) describe the stresses and moments of an elastic, compressible membrane with Helfrich energy density and viscous in-plane flow. The protein- or lipid-induced spontaneous curvature $C$ appears in both the in-plane stresses and the out-of-plane shear. Therefore, proteins or lipids with preferred spontaneous curvature are able to affect both in-plane flow and out-of-plane bending, an important observation noted in other theoretical and computational studies [60, 62, 65].

3.4 Equations of Motion

In this section, we provide the equations of motion for a membrane with elastic out-of-plane bending and intra-membrane viscous flow. The four unknowns in the membrane system are the areal mass density $\rho$ and the three components of the velocity $\bm{v}$. The four equations to be solved for the four unknowns are the local mass balance (49), two in-plane momentum balances (64), and the shape equation (65). We refer to the equations necessary to solve for all the unknowns as the governing equations. Substituting the stress and moment tensors (159)–(162) into equations (64)–(65), we find the governing equations describing the motion of the membrane to be

[TABLE]

The operator $\Delta$ is the surface Laplacian, and is defined by $\Delta(\,\cdot\,)=(\,\cdot\,)_{;\beta\alpha}a^{\alpha\beta}$. In deriving the shape equation (163), we have used the identity $\bar{b}^{\alpha\beta}b_{\alpha\beta}=2K$. We use the form of $\pi^{\alpha\beta}$ in equation (150) to calculate $\pi^{\mu\alpha}_{;\mu}$ and $\pi^{\alpha\beta}b_{\alpha\beta}$, which are found in equations (164) and (163), respectively, and are given by

[TABLE]

As expected, there are multiple modes of coupling among the equations of motion. The spontaneous curvature and its gradients appear in both the in-plane momentum equations and the shape equation. The viscosity provides a medium for coupling, as the viscous stress tensor $\pi^{\alpha\beta}$ appears in all three momentum equations (one out-of-plane (163) and two in-plane (164)). Finally, the presence of curvature and velocity components in the momentum equations couples them with the mass balance.

3.5 Boundary and Initial Conditions

In this section, we provide suitable boundary conditions necessary for solving the governing equations (163)–(165). We begin by deriving the general form of the forces and moments along the membrane edge and then evaluate them for the Helmholtz energy density in equations (155)–(153). We conclude with a discussion of the possible boundary conditions for the membrane system. A major portion of this discussion is developed in [60, 51, 59], but is provided here for completeness.

The total force $\bm{f}$ at any arbitrary boundary on the membrane surface may be decomposed in the $\{\bm{\nu},\bm{\tau},\bm{n}\}$ basis as

[TABLE]

where $f_{\nu}$ and $f_{\tau}$ are the in-plane components of the force and $f_{n}$ is the out-of-plane shear force felt by the membrane in the normal direction. The total force $\bm{f}$ at the membrane boundary is given by

[TABLE]

where $\ell$ is the arc length parametrization of the boundary and $m_{\nu}$ is defined in equation (81) [60, 51, 59, 70]. If the boundaries are piecewise continuous, and the discontinuities are indexed by $i$, the force $\bm{f}_{i}$ on the $i^{\text{th}}$ discontinuity is given by

[TABLE]

where $[m_{\nu}]_{i}$ denotes the change in $m_{\nu}$ as the $i^{\text{th}}$ corner is traversed in the forward direction, as defined by the arc length parametrization $\ell$. To decompose the force $\bm{f}$ in equation (169) in the $\{\bm{\nu},\bm{\tau},\bm{n}\}$ basis, we calculate $\mathrm{d}\bm{n}/\mathrm{d}\ell$. By using the chain rule, the definition of $\bm{\tau}$ (14), the Weingarten equation (29), and the decomposition of $\bm{a}_{\alpha}$ in equation (18), we obtain

[TABLE]

With the result of equation (171), the decomposition of the stress vectors (58), and equation (18) once again, we find

[TABLE]

The moment $M$ which contributes to the elastic behavior of the membrane at the boundary by bending the unit normal in the direction of $\bm{\nu}$ is calculated as

[TABLE]

The three components of the force provided in equations (172)–(174) and the boundary moment $M$ (175) are the general forms of the forces and moments on the membrane boundary. We now determine these quantities for more specific physical situations.

First, consider a Helmholtz energy density of the form $\bar{\psi}(\rho,H,K,T)$ for which the stresses and moments are given by equations (146)–(149). The moment $M$ is obtained by substituting equation (147) into equation (175), which yields

[TABLE]

We next define the normal curvatures in the $\bm{\nu}$ and $\bm{\tau}$ directions, $\kappa_{\nu}$ and $\kappa_{\tau}$, as $\kappa_{\nu}=b^{\alpha\beta}\nu_{\alpha}\nu_{\beta}$ and $\kappa_{\tau}=b^{\alpha\beta}\tau_{\alpha}\tau_{\beta}$, as well as the twist $\xi=b^{\alpha\beta}\nu_{\alpha}\tau_{\beta}$ [60]. Using these definitions along with the identities $H=(\kappa_{\nu}+\kappa_{\tau})/2$ and $K=\kappa_{\nu}\kappa_{\tau}-\xi^{2}$, we substitute the constitutive forms of the stresses (146)–(149) into equations (172)–(174) to obtain

[TABLE]

For the specific case of the Helmholtz energy density provided in equations (155)–(153), we substitute the partial derivatives calculated in equations (156)–(158) into equations (176)–(179) to find

[TABLE]

Now that the boundary force $\bm{f}$ and the boundary moment $M$ have been evaluated for our assumed form of the Helmholtz free energy density (155), we consider the boundary conditions for the governing equations (163)–(165). As the membrane behaves as a fluid in-plane, we may specify either the tangential velocities $v^{\alpha}$ or the in-plane components of the force, $f_{\nu}$ (181) and $f_{\tau}$ (182), at the patch boundary. The shape equation, on the other hand, is an elastic bending equation, and therefore two boundary conditions need to be specified at every point along the boundary. The simplest way to do so is to specify the membrane position and its gradient in the $\bm{\nu}$ direction, or to specify the moment $M$ and the shear force $f_{n}$ at the boundary.

We have now provided the governing equations with possible boundary conditions for both the tangential and shape equations, and close the problem by providing possible initial conditions. To this end, we specify the initial position $\bm{x}$ and the initial velocity $\bm{v}$ everywhere on the membrane surface, and furthermore assume the density $\rho$ is initially constant, the Jacobian $J=1$ everywhere, and the spontaneous curvature $C=0$ for problems initially without proteins. With these example initial conditions, and the suggested boundary conditions, the governing equations (163)–(165) are mathematically well-posed.

3.6 Coupling to Bulk Fluid

We conclude our theoretical developments for single-component lipid membranes by discussing the coupling of the membrane equations of motion with the surrounding fluid. While the bulk fluid provides an additional dissipative mechanism via the bulk viscosity and dominates dissipation for long wavelength undulations [79, 80, 81, 82], the bulk fluid can contribute negligibly in systems with membrane deformations on the order of five microns or less [53, 52, 39]. Consequently, the surrounding fluid can sometimes be excluded when studying small length scale phenomena.

In cases where bulk dissipation is non-negligible, the membrane and surrounding fluid should be modeled together and coupled through interface conditions. The equations of motion for the surrounding bulk fluid are the Navier-Stokes equations. We denote quantities in the bulk with a subscript ‘b’ and label the fluid above and below the membrane with a superscript ‘$+$’ and ‘$-$’, respectively. Thus, $\bm{v}_{\mathrm{b}}^{+}$ and $\bm{v}_{\mathrm{b}}^{-}$ are the bulk fluid velocities on either side of the membrane, and we make the no-slip assumption between the membrane and fluid such that

[TABLE]

at the membrane-fluid interface. Additionally, the stresses in the bulk fluid, denoted by the tensors $\bm{\sigma}_{\mathrm{b}}^{+}$ and $\bm{\sigma}_{\mathrm{b}}^{-}$, enter the membrane equations of motion through the body force $\bm{b}$. A force balance on the membrane yields

[TABLE]

where $\bm{n}$ is the membrane normal. By decomposing the body force in the $\{\bm{a}_{\alpha},\bm{n}\}$ basis as in equation (62), we find the pressure drop $p$ across the membrane to be given by

[TABLE]

The pressure drop and in-plane body forces in equations (187)–(188) enter the equations of motion of the membrane (163)–(164), and the velocities in both the bulk fluid and membrane are solved in a self-consistent manner to satisfy equations (184)–(188). A computational implementation of the aforementioned conditions at a fluid-structure interface is provided in [83] in the context of liquid menisci and elastic membranes. For small deformations of the membrane, one can use the Oseen tensor to couple the bulk fluid and the membrane, as done in [81, 80].

4 Intra-membrane Phase Transitions

Biological membranes may consist of hundreds of different protein and lipid constituents and exhibit complex behavior in which species diffuse to form heterogeneous domains, which then take part in important biological phenomena [17, 18]. In order to better understand such processes, many experiments have been carried out on artificially created giant unilamellar vesicles (GUVs), the composition of which may be precisely controlled. Of particular interest is the phase transition in which a membrane initially in a liquid-disordered (Ld) phase develops, under a suitable change of external conditions, phase coexistence between Ld phases and liquid-ordered (Lo) phases. The nature of these phases is governed by their concentrations—the Ld phase generally consists of a low melting temperature phospholipid such as Dioleoylphosphatidylcholine (DOPC) while the Lo phase generally primarily consists of a high melting temperature phospholipid such as Dipalmitoylphosphatidylcholine (DPPC) as well as cholesterol. A line tension exists at the Lo–Ld phase boundary, and imaging experiments have shown one of the two phases may bulge out to reduce this line tension [22, 26], thus indicating a coupling between membrane bending, diffusion, and flow (Figure 3).

In this section, we extend the irreversible thermodynamic framework of the single-component model to have multiple phospholipid constituents and determine the new forms of the balances of mass, linear momentum, angular momentum, energy, and entropy. We then propose a general constitutive form of the Helmholtz free energy density and determine the constitutive relations for stresses and moments as well as phenomenological relations for the diffusive fluxes. We conclude this general treatment by determining restrictions on the Helmholtz free energy due to invariance postulates. Finally, we study a two-component membrane model in order to better understand the coupling between bending, diffusion, and flow in membranes which exhibit Lo–Ld phase coexistence, as described above. For these cases, we provide the equations of motion and suitable boundary conditions. While we focus on multi-component phospholipid membranes, all the analysis in this section can be applied to the diffusion and segregation of transmembrane proteins as well.

4.1 Kinematics

The membrane is modeled as having $N$ types of phospholipid constituents, which are able to diffuse in the plane of the membrane. We index the phospholipid constituents by $k$, where $k\in\{1,2,\dots,N\}$. In treating multi-component systems, we follow the procedure outlined in [68].

The mass density of species $k$ is denoted by $\rho_{k}(\theta^{\alpha},t)$, and the total mass density $\rho$ is given by

[TABLE]

The barycentric velocity $\bm{v}$, also called the center-of-mass velocity, is defined by the relation

[TABLE]

where $\bm{v}_{k}=v_{k}^{\,\alpha}\bm{a}_{\alpha}+v_{k}\bm{n}$ is velocity of species $k$. While different species may have different in-plane velocity components $v^{\,\alpha}_{k}$, the normal velocity components $v_{k}$ must be the same for all species.

As we track an infinitesimal membrane patch over time, our reference frame moves with the barycentric velocity $\bm{v}$. There exists an in-plane diffusive flux of species $k$ across the patch boundary due to the difference in the species and barycentric velocities. The diffusive flux of species $k$, $\bm{j}_{k}$, is given by

[TABLE]

This flux lies in the plane of the membrane and can be written in component form as

[TABLE]

Summing all the diffusive fluxes in equation (191), we obtain

[TABLE]

using equation (190). We define the mass fraction $c_{k}$ of each species $k$ as

[TABLE]

such that

[TABLE]

As only $N-1$ of the mass fractions are linearly independent, we can completely define the composition of the membrane by specifying the mass density $\rho$ and mass fraction $c_{k}$ for $k=2,3,4,\dots,N$.

4.2 Balance Laws

We follow the same general procedure as in the single-component case to determine local forms of the mass, linear momentum, and angular momentum balances. We employ mixture theory, as presented in de Groot & Mazur [68], to develop a wholistic model of membrane behavior without modeling the individual forces between different phospholipid species. We assume no chemical reactions and defer their study to the last part of this work (Section 5).

4.2.1 Mass Balance

The mass of a single species $k$ in a membrane patch $\mathcal{P}$ may only change due to a diffusive mass flux at the boundary. In this case, the global form of the conservation of mass of species $k$ is given by

[TABLE]

Applying the Reynolds transport theorem (45) and the surface divergence theorem (46) to equation (196), we obtain

[TABLE]

In equation (197) and from now on, $v^{\alpha}$ and $v$ refer to the tangential and normal components of the barycentric velocity $\bm{v}$, respectively, as decomposed in equation (33). Since the membrane patch $\mathcal{P}$ in equation (197) is arbitrary, we obtain the local form of the conservation of mass of species $k$ as

[TABLE]

Summing equation (198) over all $k$ and realizing the right hand side is zero due to equation (193), we recover the local form of the total mass balance to be given by

[TABLE]

Using the local forms of the species mass balance (198) and total mass balance (199) as well as the definition of $c_{k}$ (194), we obtain the balance of the mass fraction of species $k$ as

[TABLE]

4.2.2 Linear Momentum Balance

For a membrane with multiple components, each species $k$ is acted on by a body force $\bm{b}_{k}(\theta^{\alpha},t)$. The entire membrane patch is subjected to a traction $\bm{T}(\bm{x}_{\mathrm{b}},t;\bm{\nu})$ at the boundary, and the total linear momentum of the membrane patch is the sum of the linear momentum of each species as shown in equation (190). In this case, the global form of the linear momentum balance is given by

[TABLE]

Defining the mass-weighted body force $\bm{b}$ as

[TABLE]

and using the barycentric velocity $\bm{v}$ in equation (190), the global linear momentum balance (201) can be reduced to

[TABLE]

Equation (203) has the same form as the single-component global linear momentum balance (54), where now the body force and velocity are mass-averaged quantities. The boundary traction $\bm{T}$ may once again be decomposed into a linear combination of stress vectors as in equation (55), and by using the same set of procedures as in the single-component section we find the local form of the linear momentum balance is given by

[TABLE]

As equation (204) is identical to the single-component result (57), equations (58)–(60) continue to be valid descriptions of the membrane stresses.

We have so far seen the local mass balance and local linear momentum balance are identical for single-component and multi-component membranes. The director traction $\bm{M}$ may also be decomposed in the same manner as in the single-component case (69), and thus the results of the single-component angular momentum balance—namely the symmetry of $\sigma^{\alpha\beta}$ (78) and the form of $S^{\alpha}$ (79)—are appropriate for multi-component membranes as well. Furthermore, the mechanical power balance is unchanged from the single-component case (91) since it only depends on the linear momentum balance.

4.3 Thermodynamics

In this section, we follow the same procedure as in the single-component section to develop the local form of the first law of thermodynamics. The local form of the entropy balance and the second law of thermodynamics are unchanged. We end with the expression for the internal entropy production.

4.3.1 First Law—Energy Balance

The global form of the first law of thermodynamics for the multi-component membrane system is given by

[TABLE]

Equation (205) is identical to the one-component result (92), except for the force term which is now a sum over the individual species. The difference in body force contribution between the global form of the first law of thermodynamics (205) and the mechanical power balance (91) is

[TABLE]

where we have used the definition of the mass averaged body force $\bm{b}$ (202) and the definition of the in-plane diffusive flux $\bm{j}_{k}$ (191). With equation (206), the mechanical power balance (91), and assuming the decomposition of total energy $e$ into internal energy $u$ and kinetic energy given in equation (93), the global form of the first law of thermodynamics (205) simplifies to

[TABLE]

In writing the covariant component of the in-plane body force acting on species $k$, $(b_{k})_{\alpha}$, in equation (207), we have explicitly included parenthesis to avoid confusion with the curvature tensor. Comparing equation (207) with the single-component result (95) shows the only difference is the $\bm{b}_{k}\cdot\bm{j}_{k}$ term calculated in equation (206). Since the membrane patch $\mathcal{P}$ is arbitrary, the local form of the first law of thermodynamics is given by

[TABLE]

Comparing equation (208) to equation (96), one additional term has appeared—namely the sum of $(b_{k})_{\alpha}\,j_{k}^{\,\alpha}$.

4.3.2 Entropy Balance & Second Law

We define the in-plane entropy flux $\bm{J}_{\mathrm{s}}$, the external entropy supply $\rho\eta_{\mathrm{e}}$, and the internal entropy production $\rho\eta_{\mathrm{i}}$ as in the single-component model. The global entropy balance is unchanged from the single-component case (97), and consequently the local form of the entropy balance is given by equation (99). The second law of thermodynamics is unchanged as well, and is given by equation (100).

At this stage, we invoke the Helmholtz free energy density $\psi$ for the multi-component membrane system, which is of the form given in equation (101). Following an identical procedure to the single-component case described in Section 3.2.3, we find the local form of the entropy balance can be expressed as

[TABLE]

4.4 Constitutive Relations

In this section, we extend the framework developed in the single-component section. When we apply linear irreversible thermodynamics to relate thermodynamic forces to their corresponding fluxes, we find diffusive species fluxes are driven by gradients in chemical potential and non-uniform body forces across different species. We learn that the form of the stresses and moments are identical to the single-component case.

4.4.1 General Thermodynamic Variables

With $N$ types of phospholipids composing the lipid membrane, there are an additional $N-1$ degrees of freedom relative to the single-component case, which are the $N-1$ mass fractions $\{c_{k}\}_{k=2,\dots,N}:=\{c_{2},c_{3},c_{4},\dots,c_{N}\}$. For this case, the Helmholtz free energy density per unit mass for the membrane system is formally written as

[TABLE]

We define the chemical potential of species $k$, $\mu_{k}$, as

[TABLE]

where the partial derivative holds $a_{\alpha\beta}$, $b_{\alpha\beta}$, and $c_{j\neq k}$ constant, with $j$ in $\{2,3,4,\dots,N\}$. The chemical potential $\mu_{1}$ is not defined in equation (211) as $c_{1}$ is not an independent variable due to the relation (195). For the membrane system, the chemical potential describes the change in Helmholtz free energy when species $k$ is exchanged with species one, such that all other mass fractions remain constant.

Taking the material derivative of equation (210) and multiplying by $\rho$, we obtain

[TABLE]

where in obtaining the second equality we have substituted equations (200) and (211) and used the local equilibrium assumption given by

[TABLE]

We substitute equation (212) into the entropy balance (209) to find

[TABLE]

We rewrite the gradient terms in equation (214) to obtain

[TABLE]

From equation (215), we see only the first term contributes to the in-plane entropy flux $J_{\mathrm{s}}^{\,\alpha}$, given by

[TABLE]

Again, the heat source per unit mass is the only term which contributes to the external entropy supply $\rho\eta_{\mathrm{e}}$, given by

[TABLE]

The terms on the right hand side of equation (215) which contribute to neither the entropy flux nor the external entropy contribute to the rate of internal entropy production per unit area $\rho\eta_{\mathrm{i}}$, yielding

[TABLE]

Using equation (193), we simplify the second term on the right hand side of equation (218) as

[TABLE]

Substituting equation (219) into the rate of internal entropy production per unit area (218), we obtain

[TABLE]

Each term on the right hand side of equation (220) is the product of a thermodynamic force and a corresponding flux, as in equation (112). In the linear irreversible regime, one can now assume linear phenomenological relations between the thermodynamic forces and fluxes as in equation (113). We invoke the Curie principle [76, 68] and assume the fluxes $J_{\mathrm{q}}^{\,\alpha}$ and $j_{k}^{\,\alpha}$ are independent from the tensorial thermodynamic forces $\dot{a}_{\alpha\beta}$ and $\dot{b}_{\alpha\beta}$, and the tensorial fluxes are independent from the vectorial thermodynamic forces. Therefore, in the linear irreversible regime the most general relations between the vectorial fluxes and forces are given by

[TABLE]

where the phenomenological coefficients $F_{k}^{\,\alpha\beta}$ and $G^{\alpha\beta}$ describe interference between the heat flow and in-plane species diffusion. For the lipid bilayers under consideration, we assume these to be zero. In the case where there is no cross-coupling between the heat and diffusive fluxes, the in-plane fluxes simplify to

[TABLE]

We assume, as before, the thermal conductivity tensor $\kappa^{\alpha\beta}$ is isotropic and obtain

[TABLE]

As noted in the single-component case, we assume lipid bilayers do not sustain in-plane temperature gradients and equation (225) is not of relevance in our modeling of irreversible processes. Assuming further the diffusion tensors $D_{\ell}^{\,\alpha\beta}$ only affect $j_{k}^{\,\alpha}$ when $k=\ell$ and $D_{k}^{\,\alpha\beta}$ is isotropic, we obtain $D_{\ell}^{\,\alpha\beta}=D_{k}\,a^{\alpha\beta}$ when $k=\ell$ and $D_{\ell}^{\,\alpha\beta}=0$ otherwise. With these assumptions, the linear phenomenological relation for the in-plane species flux $\bm{j}_{k}$ simplifies to

[TABLE]

Equation (226) indicates gradients in chemical potentials and differences in the in-plane forces $\bm{b}_{k}$ between individual species drive the diffusive fluxes.

Finally, the tensorial terms in the rate of internal entropy production $\rho\eta_{\mathrm{i}}$ are identical to those of the single-component case (111). Therefore, following similar arguments to the single-component case, we find the form of the stresses and moments in equations (136)–(139) remain unchanged for multi-component systems.

For lipid bilayers, it can again be argued that Galilean invariance [59] or material symmetry arguments [78] lead to a change in the fundamental variables of the Helmholtz free energy density, as given by

[TABLE]

Given equation (227), we note there is again no change to the form of the stresses and moments given by equations (146)–(150), in terms of their depedendence on $\rho$, $H$, and $K$.

We have now determined the governing equations for a general multi-component membrane patch, for which the unknowns are the density $\rho$, the three components of the velocity $\bm{v}$, and the $N-1$ mass fractions $c_{k}$ for $k\in\{2,3,\dots,N\}$. The corresponding equations are given by the total mass balance (199), the three linear momentum balances (64)–(65) in which we know the equations for the stresses, and the $N-1$ linearly independent mass fraction balances for $c_{k}$ (200), where $k\in\{2,3,\dots,N\}$ and $j_{k}^{\,\alpha}$ is given by equation (226).

As the form of the tractions $\bm{T}$ and director tractions $\bm{M}$ on the patch boundary are unchanged, the boundary conditions outlined for the single-component system in Section 3.5 are again appropriate for the mechanical equations of motion of multi-component systems. As the total mass balance (199) and mass fraction balances (200) contain only first order derivatives in time, we need to provide their initial values to solve for the dynamics of the membrane patch. The required boundary conditions on $c_{k}$ will depend on the form of the Helmholtz free energy, as will be seen in the following example.

4.5 Application to Lo–Ld Phase Transitions

Now that we have found the governing equations for a general multi-component membrane, we return to the specific case of Lo–Ld phase coexistence observed in GUVs [22, 26]. As mentioned earlier, a membrane initially in an Ld phase can form interfaces between Lo and Ld regions when, for example, the temperature is quenched (Figure 3). Experimentally, such a phase transition has been observed in three-component systems [22, 23, 24].

We consider a scenario in which the three components of a membrane patch are DOPC, DPPC, and cholesterol as in [22]. The Ld phase consists of DOPC, which is a low melting temperature lipid, while the Lo phase consists of DPPC, which is a high melting temperature lipid, as well as cholesterol. Because the mass fraction of DOPC is an order parameter in that it can be used to specify which phase a region is in, we define the dimensionless order parameter $\phi$ as

[TABLE]

In doing so, we choose to model our three-component system with only two components, in which one component is DOPC while the other component is the lumped DPPC and cholesterol.

4.5.1 Helmholtz Free Energy

With two membrane components, the Helmholtz free energy density contains four terms. The first of these terms captures the energetic cost of bending, through a modified form of the single-component Helfrich energy density $w_{\text{h}}$ (152) in which the bending rigidity $k_{\mathrm{b}}$ and spontaneous curvature $C$, if it exists, depend on the local mass fraction $\phi$. Such modifications reflect experimental observations that the bending modulus of Lo and Ld phases differ [26, 84], and also that local species concentrations may affect the spontaneous curvature of the membrane [19]. The Helfrich energy density for our two-component membrane system is

[TABLE]

The Helmholtz free energy density also contains a modified form of the single-component cost of areal dilation and compression, $w_{\text{c}}$ (153), in which the compression modulus $k_{\mathrm{c}}$ is again modified to depend on the local mass fraction $\phi$. Different phases may be relatively more or less compressible, and the multi-component compression energy $w_{\text{c}}$ is written as

[TABLE]

Next, for our Helmholtz free energy to give rise to phase separation, we include a double-well potential $w_{\text{dw}}$ by considering a mean-field model of molecular interactions [85, 86, 87], given by

[TABLE]

The double-well potential contains a competition between the energetic term $\chi\phi(1-\phi)$ and the entropic term $k_{\textrm{B}}\,T\,[\phi\ln\phi+(1-\phi)\ln(1-\phi)]$, where $k_{\textrm{B}}$ is Boltzmann’s constant and $\chi$, which has units of energy, describes the mean-field interaction energies of the two components. In equation (231), we assume the temperature is low enough such that the energetic term dominates and the potential has two minima—for our model, this assumption is valid for $\chi>2\,k_{\textrm{B}}\,T$. The parameter $k_{\phi}$ tunes the relative energetic penalty of concentrations which are not energetically favorable. The factor of $1/J$ is required for $w_{\text{dw}}$ to have no density dependence when convected to the reference patch $\mathcal{P}_{0}$, in the same way as the compression energy $w_{\text{c}}$.

The final term in our two-component Helmholtz free energy accounts for the line tension at the interface between two phases. We consider a simplified view in which $\phi=1$ in the Ld phase and $\phi=0$ in the Lo phase, such that the gradient of $\phi$ is a delta function which is nonzero at every point along the phase boundary. This simple example motivates the gradient term in the Helmholtz free energy, $w_{\text{g}}$, written as

[TABLE]

In equation (232), we have chosen the simplest form of $w_{\text{g}}$ which penalizes phase boundaries and is invariant to coordinate transformations. The constant $\gamma$ describes the energetic penalty for phase boundaries, and equation (232) is the curvilinear form of the energetic contribution in standard Cartesian coordinates [85].

Combining all of the energy contributions, the total Helmholtz free energy density $\psi$ for our two-component membrane system is

[TABLE]

To determine the stresses, we note that $w_{\text{g}}=\tfrac{\gamma}{2}\phi_{,\mu}\phi^{,\mu}=\tfrac{\gamma}{2}\phi_{,\mu}\phi_{,\lambda}a^{\lambda\mu}$ and thus $\psi$ may not be completely expressed as a function of only $\rho$, $H$, $K$, $T$, and $\phi$—consequently, in equation (233) we have written $\psi$ instead of $\bar{\psi}$. We therefore use the equations for the stresses in terms of the more general variables $a_{\alpha\beta}$, $b_{\alpha\beta}$, $T$, and $\phi$ (136)–(139). Using the partial derivatives found in Table 1, we find the membrane stresses and couple-stresses are given by

[TABLE]

In equations (234)–(237), we can determine which components of the total free energy $\rho\psi$ (233) contribute to the in-plane and out-of-plane stresses through the coefficients $k_{\mathrm{b}}(\phi)$, $k_{\mathrm{g}}$, $k_{\mathrm{c}}(\phi)$, $\gamma$, and $k_{\phi}$. For example, the couple-stresses $M^{\alpha\beta}$ (235) and shear stress $S^{\alpha}$ (237) are only affected by the Helfrich bending energy $w_{\text{h}}$ (229), while the in-plane stresses $N^{\alpha\beta}$ (236) contain contributions from the bending energy $w_{\text{h}}$, the compression energy $w_{\text{c}}$ (230), and the gradient energy $w_{\text{g}}$ (232). The gradient energy contributions are Korteweg-like, and describe in-plane momentum transfer due to line tensions and concentration gradients in the system [88, 89, 90, 91].

The double-well potential $w_{\text{dw}}$ (231) does not enter the stresses, but affects the equations of motion through the chemical potential $\mu_{\phi}=\partial\psi/\partial\phi$, calculated as

[TABLE]

The chemical potential $\mu_{\phi}$ thus contains contributions from the bending energy, compression energy, and gradient energy.

4.5.2 Equations of Motion

For our two-component model, there are five unknowns: the total mass density $\rho$, the mass fraction $\phi$, and the three components of the velocity $\bm{v}$. The five corresponding governing equations are the total mass balance, the mass fraction balance of $\phi$, and the three components of the linear momentum balance. The total mass balance (199) is given by

[TABLE]

where we have expanded the material derivative to show how the in-plane velocities enter the mass balance. The mass fraction balance of $\phi$ is given by equation (200), where we substitute the form of the diffusive flux found in equation (226) to obtain

[TABLE]

where the chemical potential $\mu_{\phi}$ is given by equation (238). The chemical potential contains information from all four components of the total Helmholtz free energy density, and is coupled to the concentration and curvatures. Substituting the stresses and moments (234)–(237) into the tangential equations (64) and the shape equation (65), we obtain

[TABLE]

and

[TABLE]

From the governing equations (239)–(242), we can clearly see the highly nonlinear coupling between species diffusion, bending, and intra-membrane viscous flow. We note a similar model was proposed by Agrawal and Steigmann [61], in which viscous stresses were neglected. As the presence of multiple components introduces diffusive fluxes, evident from equation (240), and therefore velocity gradients, our model indicates viscous stresses cannot be omitted in describing multi-component membranes.

4.5.3 Boundary and Initial Conditions

In this section, we specify boundary conditions on the governing equations (239)–(242), as appropriate, for them to be mathematically well-posed.

The mass fraction balance for $\phi$ (240) has a fourth-order spatial derivative, as it contains the surface Laplacian of the chemical potential, which itself contains the surface Laplacian of the mass fraction $\phi$ (238). This is similar in structure to the shape equation, and for a well-posed problem we must specify two boundary conditions in $\phi$ at every point on the boundary. For example, we may specify $\phi$ and its gradient in the $\bm{\nu}$ direction, or the chemical potential $\mu_{\phi}$ and the diffusive flux in the $\bm{\nu}$ direction, $j_{\phi}^{\,\alpha}\,\nu_{\alpha}$, where the diffusive flux $j_{\phi}^{\,\alpha}$ is given by equation (226). Such boundary conditions are studied for systems described by Cartesian coordinates in [92].

The momentum equations are of the same general form as in the single-component case. For the in-plane momentum equations (242), we may specify either the velocities $v^{\alpha}$ or the in-plane force components $f_{\nu}$ and $f_{\tau}$ of the total force $\bm{f}$ at the boundary, given in equations (168)–(169). For the shape equation (241), we may specify the position and its gradient or the shear force $f_{n}$ and the moment $M$ (175) at the boundary. We find the functional form of the forces and moments by substituting the stresses (234)–(237) into the general relations (172)–(175) and obtain

[TABLE]

We close the problem by providing example initial conditions. As in the single-component case, we specify the initial membrane position $\bm{x}$, velocity $\bm{v}$, and assume initially the density $\rho$ is constant and the Jacobian $J=1$ everywhere. One may also specify an initial distribution of the mass fraction $\phi$ everywhere, which determines the spontaneous curvature $C(\phi)$. With these initial conditions and the aforementioned boundary conditions, our problem is well-posed.

5 Peripheral Proteins—Chemical Reactions

The binding and unbinding of peripheral proteins plays a crucial role in many cellular processes involving lipid membranes. In vitro experiments have shown dramatic shape changes can occur as a result of protein binding and unbinding events involving the membrane and the surrounding fluid [14, 15, 13, 19, 93, 94]. Once proteins bind to the membrane, they are able to diffuse in-plane as well. The complex interplay between protein binding, in-plane diffusion, in-plane lipid flow, and membrane bending is currently not well-understood. In this section, we extend the multi-component model to include binding and unbinding chemical reactions of proteins onto membranes. We determine the new form of the balance laws and equations of motion, and learn the thermodynamic driving force governing the binding and unbinding reactions. We motivate our theoretical development with the following example.

5.1 Chemistry

Consider the membrane to consist of two types of phospholipids: DOPC and Phosphatidylinositol-4,5-biphosphate (PI(4,5)P2). In the fluid surrounding the membrane, there is a reservoir of epsin-1 proteins which bind only to PI(4,5)P2 lipids [27, 7]. The concentration of PI(4,5)P2 is usually around 2% by mass under physiological conditions. Because we model the membrane as a continuum, at every point in space there may be multiple binding and unbinding events (Figure 4). In this entire theoretical development, electrostatic effects are neglected. While we study this scenario in detail, another motivating example is the binding of BAR proteins onto membranes, which has been extensively studied in [20].

Our system includes the membrane patch $\mathcal{P}$ and any epsin-1 proteins bound to it. In the foregoing example, there exist three chemical species in our system. They are labeled by the index $k$, and are provided in Table 2. The index $k=3$ refers to an epsin-1 protein bound to a PI(4,5)P2 lipid. The binding and unbinding reaction between epsin-1 and PI(4,5)P2 lipids can be written as

[TABLE]

The overall reaction rate $\mathcal{R}$ for equation (247) has units of a molar flux and is given by

[TABLE]

where $\overset{\rightarrow}{\mathcal{R}}$ and $\overset{\leftarrow}{\mathcal{R}}$ are the forward and reverse reaction rates, respectively. Defining $\mathscr{M}_{\mathrm{p}}$ as the molar mass of epsin-1 proteins, $\mathscr{M}_{\mathrm{p}}\mathcal{R}$ is the mass flux from the reservoir of proteins either above or below the membrane, herein called the bulk phase, to the membrane surface.

The concentrations $n_{k}$ on the membrane patch, with units of moles per area, are defined as

[TABLE]

where $\rho_{k}$ is the mass density of species $k$ and $\mathscr{M}_{k}$ is the molar mass of species $k$. The concentration of espin-1 proteins in the bulk phase is denoted $n_{\mathrm{p}}^{\mathrm{b}}$ and has units of moles per volume. The forward and reverse chemical reactions in equation (247) are assumed to be elementary steps, for which the forward and reverse reaction rates may be written as

[TABLE]

where $\overset{\rightarrow}{k}$ and $\overset{\leftarrow}{k}$ are the forward and reverse rate constants, respectively. At this stage, we introduce the stoichiometric coefficients $\alpha_{k}$, which are the coefficients of the $k^{\mathrm{th}}$ species in the chemical reaction (247) and from inspection are given by $\alpha_{1}=0$, $\alpha_{2}=-1$, and $\alpha_{3}=1$. It is notationally convenient to define the stoichiometric coefficient of the epsin-1 protein in the chemical reaction (247) as $\alpha_{\mathrm{p}}=-1$. A summary of the new variables introduced in this section is provided in Table 3.

5.2 Balance Laws

We continue to use our continuum mechanical framework, including mixture theory, to describe the membrane patch. All of the kinematic results derived in Section 4.1 continue to hold. The binding and unbinding of proteins, however, affects the overall mass balance of the membrane patch and consequently the global forms of the linear and angular momentum balances as well.

5.2.1 Mass Balance

For a membrane patch $\mathcal{P}$, the mass of species $k$ can change due to a diffusive mass flux $\bm{j}_{k}$ at the patch boundary or the binding and unbinding of proteins along the membrane surface, which is a mass flux from the bulk phase to the membrane surface. In this case, the global form of the conservation of mass of species $k$ is given by

[TABLE]

The form of the diffusive flux of species $k$ is given in equation (192) such that the sum of diffusive fluxes satisfies equation (193). Comparing equation (252) to its multi-component analog (196), we see the binding and unbinding of proteins affects the total mass of each species through the reaction rate $\mathcal{R}$. Again, we use the Reynolds transport theorem (45) and the surface divergence theorem (46) to obtain

[TABLE]

As the membrane patch $\mathcal{P}$ is arbitrary, the local form of the mass balance of species $k$ can be obtained as

[TABLE]

While we have accounted for bound proteins through the index $k=3$, which refers to epsin-1 proteins bound to PI(4,5)P2 lipids, it is possible to determine the continuity equation for proteins alone by multiplying equation (254), for $k=3$, by $\mathscr{M}_{\mathrm{p}}/\mathscr{M}_{3}$.

Summing equation (254) over all $k$, we find the local form of the total mass balance is given by

[TABLE]

indicating binding and unbinding reactions change the mass of the membrane patch only due to the addition or removal of proteins. When proteins are binding, $\mathcal{R}>0$ and the mass of the membrane patch increases, while when proteins are unbinding $\mathcal{R}<0$ and the mass of the patch decreases. Because the mass of an infinitesimal patch changes over time due to protein binding, the relation $J=\rho_{0}/\rho$ from equation (52) is no longer true. Consequently, $a_{\alpha\beta}$ determines the Jacobian determinant $J$, but does not determine the areal mass density $\rho$.

To calculate the appropriate form of the modified Reynolds transport theorem for the case of binding and unbinding proteins, we apply the Reynolds transport theorem (45) where $f=\rho u$ and substitute the local mass balance (255) to obtain

[TABLE]

As before, in equation (256) $u$ can be a scalar-, vector-, or tensor-valued function.

5.2.2 Linear Momentum Balance

As our system now consists of the membrane patch and any proteins bound to it, the binding and unbinding of proteins will affect the total linear momentum of the membrane patch. The mass flux of proteins from the bulk to the membrane, $\mathscr{M}_{\mathrm{p}}\mathcal{R}$, carries momentum to and from the membrane surface. Assuming a no-slip condition between the membrane surface and the adjacent fluid, the proteins travel at an average velocity $\bm{v}$ just before binding and just after unbinding. Correspondingly, the momentum flux from the bulk to the membrane is $\bm{v}\mathscr{M}_{\mathrm{p}}\mathcal{R}$, and the global form of the linear momentum balance is given by

[TABLE]

Using the Reynolds transport theorem (256) and the definition of the mass-weighted body force provided in equation (202) with $N=3$, equation (257) can be reduced to

[TABLE]

Equation (258) is identical to the global form of the linear momentum balance of a multi-component system, found in Section 4.2.2 to be given by equation (54). Since the membrane patch $\mathcal{P}$ is arbitrary, the local form of the linear momentum balance may be written as

[TABLE]

5.2.3 Angular Momentum Balance

Similar to the linear momentum analysis, the additional change in the angular momentum due to the chemical reactions is $\bm{x}\times\bm{v}\mathscr{M}_{\mathrm{p}}\mathcal{R}$. Therefore, the global form of the angular momentum balance is given by

[TABLE]

Applying the Reynolds transport theorem (256) to the left hand side, we obtain

[TABLE]

which is identical to the global form of the multi-component angular momentum balance found in Section 4.2.2 to be given by equation (68). The local form of the angular momentum balance for our membrane system is therefore identical in form to the single-component case (75). As the local linear momentum balance (259) is also unchanged, the symmetry of $\sigma^{\alpha\beta}$ (78), form of $S^{\alpha}$ (79), and mechanical power balance (91) are valid in the presence of chemical reactions as well.

5.3 Thermodynamics

In this section, we develop local forms of the first and second laws of thermodynamics as well as the local entropy balance. Thermodynamic balance laws must now take into account energy and entropy changes due to chemical reactions.

5.3.1 First Law—Energy Balance

In addition to the mechanisms described in the multi-component case, the total energy of the membrane patch can change due to the flux of proteins from the bulk phase to the membrane surface. Assuming the changes in energy due to the binding and unbinding of proteins can be decomposed into kinetic and internal energy components and there is no slip between the membrane surface and the surrounding fluid, these changes are given by $(\tfrac{1}{2}\bm{v}\cdot\bm{v}+u)\mathscr{M}_{\mathrm{p}}\mathcal{R}$. Correspondingly, the global form of the energy balance can be written as

[TABLE]

In equation (262), only the last integral term is new compared to the multi-component form (205). Applying the Reynolds transport theorem (256) and substituting the definition of the total energy (93), we obtain

[TABLE]

The right hand side of equation (263) is equal to the right hand side of the multi-component global energy balance (205). Following similar arguments to those of the multi-component case described in Section 4.3.1, the local form of the first law of thermodynamics is found to be

[TABLE]

Equation (264) is identical to the local energy balance of the multi-component membrane (208), for $N=3$.

5.3.2 Entropy Balance & Second Law

As in the momentum balances and the energy balance, in the global entropy balance we introduce a new term to account for the entropy changes due to the chemical reaction. To this end, the entropy is modified by $s\mathscr{M}_{\mathrm{p}}\mathcal{R}$, and accordingly the global form of the entropy balance is given by

[TABLE]

Applying the Reynolds transport theorem (256) and the surface divergence theorem (46) to equation (265), we obtain

[TABLE]

As the membrane patch $\mathcal{P}$ is arbitrary, the local form of the entropy balance is

[TABLE]

The second law of thermodynamics is still given by equation (100).

5.3.3 Choice of Thermodynamic Potential

As in the previously considered single- and multi-component membranes, we express the entropy balance in terms of the Helmholtz free energy density $\psi$. Using equations (101) and (267) we obtain the total rate of change of entropy to be

[TABLE]

Substituting the local form of the first law of thermodynamics (264) into equation (268), we obtain

[TABLE]

Equation (269) is identical to the multi-component entropy production, given by equation (209), for $N=3$.

5.4 Constitutive Relations

In this section, we extend the framework developed in the multi-component analysis to determine the internal entropy production equation. We again apply linear irreversible thermodynamics to relate thermodynamic forces and fluxes, and determine the stresses, moments, and diffusive species fluxes. We also determine the thermodynamic driving force governing the binding and unbinding of proteins in the linear irreversible regime, and then extend our analysis to apply to systems arbitrarily far from chemical equilibrium.

5.4.1 General Thermodynamic Variables

As mentioned earlier, due to the binding and unbinding of proteins, the metric tensor $a_{\alpha\beta}$ determines the Jacobian determinant $J$ but not the areal mass density $\rho$. Consequently, the Helmholtz free energy per unit mass $\psi$ may depend on all the species densities $\rho_{1}$, $\rho_{2}$, and $\rho_{3}$ in addition to $a_{\alpha\beta}$, and may be written as

[TABLE]

Taking the material derivative of $\psi$ in equation (270) and multiplying by $\rho$, we obtain

[TABLE]

where in obtaining the second equality, we used the local equilibrium assumption

[TABLE]

defined the chemical potential as

[TABLE]

and substituted the expression for $\dot{\rho}_{k}$ given by equation (254). The chemical potential $\mu_{k}$ (273) has units of energy per unit mass, and is now defined for all three species, as $\rho_{1}$, $\rho_{2}$, and $\rho_{3}$ are linearly independent.

Substituting equation (271) into the entropy balance (269) and rearranging terms yields

[TABLE]

As outlined by Prigogine [67], we introduce the chemical affinity $\mathscr{A}$ for the chemical reaction in equation (247), which was first introduced by de Donder [95] as

[TABLE]

where $\mu_{\mathrm{p}}^{\mathrm{b}}$ is the chemical potential of proteins in the bulk, again defined to have units of energy per unit mass. As described by de Groot & Mazur [68], because the diffusive fluxes $\bm{j}_{k}$ sum to zero (193), for any scalar quantity $f_{k}$ we can write

[TABLE]

Equation (276) also holds for vector or tensor quantities. Substituting equation (275) into equation (274), reorganizing terms using equation (276), recognizing $(v^{\alpha}_{;\alpha}-2vH)$ contains $\dot{a}_{\alpha\beta}$ due to equation (44), and rearranging terms gives

[TABLE]

By inspection of equation (277), the in-plane entropy flux $\bm{J}_{\mathrm{s}}$ is

[TABLE]

The external entropy supply now contains a new contribution from the bulk phase, namely $\alpha_{\mathrm{p}}\mu_{\mathrm{p}}^{\mathrm{b}}\,\mathscr{M}_{\mathrm{p}}\mathcal{R}/T$, in addition to $\rho r/T$. The rate of change of mass of the membrane is given by $\mathscr{M}_{\mathrm{p}}\mathcal{R}$ and is equivalent to the diffusive flux of proteins from the bulk in the direction normal to the membrane. Therefore, $\alpha_{\mathrm{p}}\mu_{\mathrm{p}}^{\mathrm{b}}\,\mathscr{M}_{\mathrm{p}}\mathcal{R}/T$ is the entropy flux at the boundary of the bulk phase [68]. Accordingly, the external entropy supply $\rho\eta_{\mathrm{e}}$ is given by

[TABLE]

where we have substituted $\alpha_{\mathrm{p}}=-1$.

Next, the terms on the right hand side of equation (277) which have not contributed to the entropy flux or the external entropy contribute to the internal entropy production. To this end, the rate of internal entropy production per unit area $\rho\eta_{\mathrm{i}}$ is given by

[TABLE]

Equation (280) is the sum of terms which are a product of a thermodynamic force and a corresponding flux, where the forces and fluxes are either scalars, vectors, or tensors. We invoke the Curie principle [76, 68] and assume the phenomenological coefficients connecting terms of different tensorial order are zero. As the vectorial and tensorial terms in equation (280) are similar to the multi-component internal entropy production (218), we follow the procedure described in Section 3.3.1 to find the relations between the vectorial and tensorial thermodynamic forces and fluxes as

[TABLE]

with $\pi^{\alpha\beta}$ given by equation (150). Equation (282) holds for $k\in\{2,3\}$, and $j_{1}^{\,\alpha}$ is calculated as $j_{1}^{\,\alpha}=-j_{2}^{\,\alpha}-j_{3}^{\,\alpha}$. It can be seen from equation (283) that the in-plane stresses now contain additional contributions from the adsorption of proteins onto the membrane.

There is only a single scalar thermodynamic term in equation (280), namely that involving the chemical affinity $\mathscr{A}$, and the phenomenological relationship between the chemical affinity and the reaction rate is obtained as

[TABLE]

where $L$ is a positive constant. Equation (287) is analagous to the result for chemical reactions in three-dimensional Cartesian systems, as determined by Prigogine [67]. While linear irreversible thermodynamics may provide us with phenomenological relations with a wide range of validity, for chemical reactions nonlinear effects become significant soon after departure from equilibrium. Established theory allows us to operate outside the linear irreversible regime [67, 87]. To this end, the chemical potential of an individual species $k$ is defined as

[TABLE]

Here $\mu_{k}^{\circ}(T)$ is a chemical potential under a set of standard conditions, which is a function of only the temperature, and $R$ is the ideal gas constant. The activity of component $k$, denoted as $a_{k}$, is simply the concentration $n_{k}$ for ideal systems. The equilibrium constant $K_{\mathrm{eq}}(T)$ for the reaction is given by

[TABLE]

and is related to the forward and backward reaction rate constants, $\overset{\rightarrow}{k}$ and $\overset{\leftarrow}{k}$, through the equation

[TABLE]

We now have several relations which relate the reaction rates to rate constants (248)–(251), chemical affinity (275), and thermodynamic quantities between the membrane and bulk (288)–(290). Combining these equations, we obtain

[TABLE]

In developing equation (291), none of the arguments required the system be close to equilibrium. The result is therefore generally valid and describes how the reaction rate is related to a thermodynamic force. In the limit

[TABLE]

we Taylor expand the reaction rate (291) about equilibrium to find

[TABLE]

where $\overset{\rightarrow}{\mathcal{R}}_{\mathrm{e}}$ is the forward reaction rate at equilibrium ($\overset{\rightarrow}{\mathcal{R}}_{\mathrm{e}}=\overset{\leftarrow}{\mathcal{R}}_{\mathrm{e}}$). Comparing equation (293) with equation (287), we find the phenomenological constant $L$ is given by $L=\overset{\rightarrow}{\mathcal{R}}_{\mathrm{e}}/R$ in cases where equation (292) is valid.

5.4.2 Helmholtz Free Energy—Change of Variables

As in the single- and multi-component cases, the requirement for the Helmholtz free energy density $\psi$ to be an absolute scalar field places restrictions on its functional form. Under Galilean invariance the Helmholtz free energy density may depend on $a_{\alpha\beta}$ and $b_{\alpha\beta}$ only through $J$, $H$, and $K$, where $\partial J/\partial a_{\alpha\beta}=\tfrac{1}{2}Ja^{\alpha\beta}$ and $\partial J/\partial b_{\alpha\beta}=0$. In this case, the Helmholtz free energy density may be written as

[TABLE]

Substituting equation (294) into the stresses found in equations (283)–(286) and using an analogous procedure to the derivation of equations (146)–(149), we obtain

[TABLE]

5.4.3 Helfrich Energy Density

For the membrane system under consideration, we assume the bound proteins do not exhibit a phase transition. Consequently, the Helmholtz free energy contains neither the gradient contribution (232) nor the double-well contribution (231) present in the multi-component model. We further assume the spontaneous curvature in the membrane, $C$, is only due to the presence of bound proteins, and may be written as $C=C(\rho_{3})$. Finally, we assume there is an energetic penalty for protein-bound and unbound PI(4,5)P2 on the membrane patch.

Given the above assumptions about membrane energetics in the case of protein binding reactions, the total Helmholtz free energy consists of three terms. The first is the now familiar Helfrich bending energy $w_{\text{h}}$ given by

[TABLE]

where we assume the bending modulus $k_{\mathrm{b}}$ is independent of concentration. The cost of areal dilation and compression, $w_{\text{c}}$, is again written as

[TABLE]

where we also assume the compression modulus $k_{\mathrm{c}}$ is independent of concentration. Finally, we model the energetic penalty of the minority species PI(4,5)P2 and PI(4,5)P2–epsin-1 as

[TABLE]

where $k_{2}$ and $k_{3}$ are constants. Equation (301) may be considered as a simple model for understanding the energetic penalty of minority species, while a general derivation from microscopic considerations is left to future work. The total energy density is given by

[TABLE]

where as in the single-component case (155) $f(T)$ is a function of the temperature such that the local equilibrium assumption (272) holds.

Given the Helmholtz free energy density in equation (302), the chemical potentials $\mu_{k}$ defined in equation (273) can be obtained as

[TABLE]

Substituting the Helmholtz free energy density (302) and chemical potentials (303)–(305) into equations (295)–(298), we find the stresses and moments of the membrane to be

[TABLE]

5.5 Equations of Motion

In this section, we provide the equations of motion for a membrane undergoing protein binding and unbinding reactions. For the specific membrane patch under consideration, there are seven unknowns: the total mass density $\rho$, the species densities $\rho_{2}$ and $\rho_{3}$, the reaction rate $\mathcal{R}$, and the three components of the velocity $\bm{v}$. We choose the total mass density $\rho$, rather than the DOPC mass density $\rho_{1}$, as a fundamental variable because we are interested in the overall membrane behavior rather than that of the chemically unreactive DOPC. The reaction rate $\mathcal{R}$ may be expressed as

[TABLE]

where we have substituted equation (250) for $\overset{\rightarrow}{\mathcal{R}}$ into equation (291). In equation (310), we assume the concentration of bulk proteins $n_{\mathrm{p}}^{\mathrm{b}}$ is known, and $n_{2}$ is given by equation (249) for $k=2$.

The conservation of total mass is given by equation (255), where we expand the material derivative to obtain

[TABLE]

The mass balance of species two and three are given by equation (254), where we can substitute the diffusive flux obtained in equation (282). Since $\mu_{2}-\mu_{1}=2\,k_{2}\,\rho_{2}$ from equations (303)–(304), the mass balance for species two is given by

[TABLE]

Calculating the chemical potential difference $\mu_{3}-\mu_{1}$ from equations (303) and (305), we obtain the mass balance for species 3 as

[TABLE]

Substituting the stresses and moments (306)–(309) into the tangential equations (64) and the shape equation (65), we obtain

[TABLE]

Equations (310)–(315) solve for the seven unknowns of the membrane patch. We note that in developing these equations, we have assumed a phase transition does not exist—and we enforced this assumption by our choice of $\psi$ in equation (302). If we were to consider a different system in which proteins exhibit phase separation, we would use the Helmholtz energy density provided in the multi-component case (233) with $\phi$ given by $\phi=\rho_{3}/\rho$. As we have already calculated the stresses and moments for such an energy density, we could easily determine the equations of motion and the consequences of a system with such a phase separation.

5.6 Boundary and Initial Conditions

We conclude the problem of protein binding and unbinding reactions by providing possible boundary conditions. As the momentum equations (314)–(315) are of the same structure as in the multi-component example (241)–(242), the momentum balance boundary conditions from the multi-component analysis of Section 4.5.3 are appropriate for a membrane patch with chemical reactions. The tangential equations require either the in-plane velocities $v^{\alpha}$ or the in-plane components of the force, $f_{\nu}$ and $f_{\tau}$, to be specified at the boundary. The shape equation becomes mathematically well-posed if we specify the position and its gradient in the $\bm{\nu}$ direction, or the moment $M$ and the shear force $f_{n}$, at the boundary. We determine the force and moment at the boundary by substituting the stresses (306)–(309) into equations (172)–(175) to obtain

[TABLE]

The species balances for PI(4,5)P2 and PI(4,5)P2–epsin-1 are given by equations (312) and (313), respectively. Both species balances contain the Laplacian of the corresponding mass density. We therefore specify either $\rho_{k}$ or the species flux at the boundary in the $\bm{\nu}$ direction, $j_{k}^{\alpha}\,\nu_{\alpha}$, where $k\in\{2,3\}$.

Finally, we note the equation for the reaction rate $\mathcal{R}$ (310) does not contain any spatial derivatives, and we only need to know the initial concentrations to determine the value of $\mathcal{R}$ at any later time. Given these initial concentrations, the initial membrane configuration and velocity, and boundary conditions, our problem is mathematically well-posed.

6 Conclusions

In this paper, we developed an irreversible thermodynamic framework for arbitrarily curved lipid membranes. We began by modeling a compressible single-component lipid membrane, including out-of-plane elastic bending and in-plane viscous fluid flow. Using the balances of mass, linear and angular momentum, energy, and entropy, we determined the entropy production for single-component membranes in terms of the viscous stresses and in-plane velocity gradients through the framework of linear irreversible thermodynamics. This framework provided a natural way to develop constitutive laws including the viscous stresses, and the resulting equations of motion are identical to the results of earlier studies which proposed constitutive forms of the viscous stresses [53, 60]. We then extended the model to include multiple components that could diffuse in the plane of the membrane. We modeled phase transitions between Lo and Ld domains, and learned how phase transitions are coupled to fluid flow, diffusion, and bending. Finally, we extended our multi-component model further to include the binding and unbinding of peripheral proteins in a biologically relevant example. We found how bending, flow, diffusion, and binding are coupled, and determined the thermodynamic driving force governing protein binding.

The current theory could be expanded in several ways to develop a more complete description of lipid membranes. First, we could model lipid membranes as two monolayer leaflets rather than a single sheet as in [46, 39, 52]. By accounting for individual monolayer leaflet behavior, we would understand how inter-monolayer friction and asymmetries between the monolayers affect the equations of motion. Furthermore, we could model the bulk viscosity of the fluid surrounding the membrane and understand how an additional dissipative mode affects membrane dynamics as in [42, 96]. Finally, we could model phase transitions involving proteins in which proteins, once bound, can separate into high-density and low-density phases. While this is a simple theoretical extension from our current Lo–Ld phase transition model, phase transitions due to proteins are important in many biological phenomena. Computational microscopic studies have recently demonstrated the importance of protein-lipid and protein-protein interactions in understanding membrane bending [97, 98, 99, 100], and have also shown a force of assembly between proteins which favor the disordered phase, but which have been placed in the ordered phase [101]. To this end, it may be of interest to understand how protein interactions in the presence of Lo and Ld phases can affect bending, flow, and the collective reorganization of the membrane.

In addition to further theoretical advances, the theoretical framework presented in this paper can be used to develop numerical methods such as finite element methods to simulate membrane behavior in different biological processes such as endocytosis, intra-cellular trafficking, and cell-cell signaling. An initial effort along these lines is presented in [62].

Finally, we note the theoretical framework presented here may be generally applied to systems beyond lipid membranes. For example, we could apply this procedure to understand the dynamics of the cell wall, which plays a central role in providing structure as well as interacting with the cell membrane. As the cell wall has elastic in-plane behavior in addition to elastic out-of-plane bending, the constitutive assumptions will differ from those of lipid membranes but can be easily included into the framework developed here.

Acknowledgements

We thank Eva Schmid, David Limmer, and Clay Radke for useful discussions and Joel Tchoufag for carefully reading the manuscript. K.K.M. acknowledges the support of the University of California, Berkeley and National Institutes of Health Grant No. R01-GM110066. R.A.S. acknowledges the support of the German Research Foundation (DFG) through Grant No. GSC 111. A.S. is supported by the Computational Science Graduate Fellowship from the U.S. Department of Energy.

Appendix A Surface parametrizations

A.1 Convected Coordinates

To better understand the connection between kinematic quantities on the current and reference membrane configurations, as done in [60, 56], we begin by introducing the convected coordinates $\xi^{\alpha}$ which parameterize the reference configuration—defined to be the membrane patch at some fixed time $t_{0}$. The convected coordinates $\xi^{\alpha}$ are defined as

[TABLE]

At later times the membrane patch will in general occupy a different configuration, yet any such configuration can be mapped back to the reference patch. It is therefore appropriate to talk about how a point with constant $\xi^{\alpha}$ moves in time, as these coordinates are convected along with material points—hence our choice of the name “convected coordinates.” Since a point with constant convected coordinates $\xi^{\alpha}$ always refers to the same material point, $\xi^{\alpha}$ can also be called a Lagrangian or referential parametrization.

To understand how a point with constant fixed surface coordinates $\theta^{\alpha}$ moves in time, we consider the relationship between the two parametrizations $\theta^{\alpha}$ and $\xi^{\alpha}$. The fixed surface coordinates $\theta^{\alpha}$ change over time as the membrane deforms, and can be formally written as

[TABLE]

The above equation allows us to relate velocities in terms of the two parametrizations. To begin, the membrane position $\bm{x}$ is expressed in terms of the surface coordinates as $\bm{x}(\theta^{\alpha},t)$ and in terms of the convected coordinates as $\hat{\bm{x}}(\xi^{\alpha},t)$, with

[TABLE]

where the ‘hat’ denotes a quantity expressed in terms of the convected coordinates $\xi^{\alpha}$. The velocity $\bm{v}$ is formally defined as the rate of change of position of a material point, written as

[TABLE]

By expressing $\hat{\bm{x}}(\xi^{\beta},t)$ in terms of the fixed surface coordinates as $\bm{x}(\theta^{\alpha},t)$ and applying the chain rule, we obtain

[TABLE]

We define the in-plane contravariant velocity components $v^{\alpha}$ as the rate of change of the fixed surface coordinates for a given $\xi^{\beta}$, given by

[TABLE]

and moreover require the fixed surface coordinates $\theta^{\alpha}$ to be chosen such that

[TABLE]

so points of constant $\theta^{\alpha}$ move only in the direction normal to the membrane. The velocity can accordingly be decomposed as

[TABLE]

as done in equation (33), on substitution of equations (A.6)–(A.7) into equation (A.5).

To further explain the relationship between the fixed surface and convected parametrizations, we consider a scalar function $f$ which can be written as $f(\theta^{\alpha},t)$ or $\hat{f}(\xi^{\beta},t)$. The material derivative of $f$ is formally defined as

[TABLE]

which on substitution of $f(\theta^{\alpha},t)$ and application of the chain rule leads to

[TABLE]

in agreement with equation (32).

A.2 Material Time Derivative of In-Plane Basis Vectors

In what follows, we derive explicitly the expression for the material time derivative of $\bm{a}_{\alpha}$ using the convected coordinates $\xi^{\alpha}$. We first express $\dot{\bm{a}}_{\alpha}$ using the formal definition of the material derivative in terms of the convected coordinates, which yields

[TABLE]

which is the statement of equation (35).