Phase-field theories for mathematical modeling of biological membranes

G Model

ARTICLE IN PRESS

CPL 4318 1–15

Chemistry and Physics of Lipids xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Chemistry and Physics of Lipids journal homepage: www.elsevier.com/locate/chemphyslip

Phase-field theories for mathematical modeling of biological membranes

1

2

3

Q1

Guillermo R. Lázaro a,∗ , Ignacio Pagonabarraga b , Aurora Hernández-Machado a a

4

b

5

Departament d’Estructura i Constituents de la materia, Universitat de Barcelona, Av. Diagonal 645, E08028 Barcelona, Spain Departament de Fisica Fonamental, Universitat de Barcelona, Av. Diagonal 645, E08028 Barcelona, Spain

6

a r t i c l e

7 16

i n f o

a b s t r a c t

8

Article history: Available online xxx

9 10 11

15

Keywords: Phase-field model Biological membranes Membrane elasticity

17

1. Introduction

12 13 Q3 14

Biological membranes are complex structures whose mechanics are usually described at a mesoscopic level, such as the Helfrich bending theory. In this article, we present the phase-field methods, a useful tool for studying complex membrane problems which can be applied to very different phenomena. We start with an overview of the general theory of elasticity, paying special attention to its derivation from a molecular scale. We then study the particular case of membrane elasticity, explicitly obtaining the Helfrich bending energy. Within the framework of this theory, we derive a phase-field model for biological membranes and explore its physical basis and interpretation in terms of membrane elasticity. We finally explain three examples of applications of these methods to membrane related problems. First, the case of vesicle pearling and tubulation, when lipidic vesicles are exposed to the presence of hydrophobic polymers that anchor to the membrane, inducing a shape instability. Finally, we study the behavior of red blood cells while flowing in narrow microchannels, focusing on the importance of membrane elasticity to the cell flow capabilities. © 2014 Elsevier Ireland Ltd. All rights reserved.

Like many other organelles at the cell scale, membranes are composite structures that exhibit a bewildering complexity. Their basic ingredient is a lipid bilayer composed by hundreds of different lipid species. The bilayer also contains a dense population of transmembrane proteins, which could represent up to 70% of the total mass of the membrane, and in fact these molecules define the functionality of the membrane (Alberts et al., 1994). Many other proteins are anchored to both sides of the bilayer, and membrane composition is balanced by lipid reservoirs which ensure that the physiological properties of the membrane are maintained. Among others, membranes define the cell frontiers, separating the cytosol from the external environment. They also maintain ion gradients which are necessary to produce ATP, and host the proteins that control cell signaling. Focusing specifically on the structural function, membranes present a delicate interplay with the cortex cytoskeleton, a complex mesh formed by filaments of actin preserving cell inner structure and shape, and provides strength and compactness. This picture represents, however, just a rough description of the membrane composition and

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Q2

∗ Corresponding author. E-mail address: [email protected] (G.R. Lázaro).

function, included here to evidence its extreme complexity. The comprehensive understanding of this fascinating system requires of different level of approaches. At the molecular scale, the detailed running of each microstructure can be analyzed thoroughly from the electrochemical interactions between their molecular components. However, the all-encompassing response of the membrane elements invites to a more general description, and physical theories such as the elastic formalism of plates and surfaces offer a formidable tool to characterize biological membranes. At this point, it is convenient to remark that physicists have focused on the specific study of the human erythrocyte (Sackmann, 1995). This cell is unique among the rest of cells of the organism because it lacks any organelle and inner structure, so that all its physical properties are entirely determined by its membrane, representing a much simpler structure than normal cells. The membrane of the erythrocyte is composed by a lipid bilayer with and underlaying spectrin cytoskeleton which anchors to the cytosolic side of the bilayer. This two-dimensional scaffold has a structural function, preventing the cell from vesiculation and large deformations. In order to build a physical theory of the membrane, the complexity of the cell membrane suggests to consider the scales of interest. Our scope is to study phenomena at the cell scale, such as cell morphological deformations or mechanical interactions with the environment. In this context, the atomic description is clearly unaffordable: the difficulty of dealing with such a vast number

http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001 0009-3084/© 2014 Elsevier Ireland Ltd. All rights reserved.

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

2 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

98

99 100 101 102 103 104 105 106 107 108 109 110

111

112 113 114 115 116 117 118 119 120 121 122

of interactions between atoms, even in simple molecules such as lipids, discards any treatment at this scale. Coarse-grained descriptions, which represent each lipid by a number of beads (typically 3–10) that encompass a region of the molecule with similar properties (Marrink et al., 2007; Shillcock and Lipowsky, 2006), offer a path for the study of small sized patches of membranes. The state-of-the-art numerical methods are able to describe the kinetics of typically 106 molecules (Marrink and Tieleman, 2013), involving membrane domains of roughly 100 nm × 100 nm, but still far from the macroscopic cell scale, 10 ␮m × 10 ␮m. It is clear that coarse-grained methods are, still, not appropriate if one pretends to study the overall cell response. For this purpose, it is convenient to invoke mesoscopic theories (Deserno, 2009). By considering the membrane as locally homogeneous and introducing a continuum description, each small part of the membrane is characterized by some certain local properties. These properties must be consistent with the local molecular structure of the membrane, so that a connection between the micro and mesoscales should be derived. This Chapter does not intend to describe in detail the basis of mesoscopic theories in the context of membranes, a subject which has received extensive attention elsewhere (Safran, 1994), but to motivate one of the most relevant methods for studying interface dynamics, the so-called phase field method, as an important tool in the study of membrane elasticity and kinetics. Phase-field models for biological membranes are based on the Helfrich description of membranes (Helfrich, 1973), a modified formulation of the theory of elasticity. We will start presenting this theory and studying its application to lipid membranes, in order to gain some intuition about the physics and elasticity of these structures, and paying special attention to its derivation from the molecular description of the membrane. The main characteristics of membrane elasticity are addressed to be subsequently incorporated to the phase-field model. Then, we will explain the basis of the phase-field methods and how they can be used to model membranes. Finally, we will conclude with some examples of current research on biological membranes that make use of this methodology. 2. The curvature energy The theoretical study of membranes at the cell scale was first performed by Canham (1970), Helfrich (1973), and Evans (1974). They concentrated on the identification of the relevant elastic properties of the erythrocyte membrane by trying to reproduce its distinctive discocyte shape. The main assumption of their approach is that the cell membrane can be described as a two-dimensional sheet, based on its minute thickness compared to the cell length. Helfrich proposed that from the three main type of deformations that a layer can undergo, shear, tilt and bending, only the last does play a relevant role in the membrane elasticity. Accordingly, he proposed a curvature energy to describe the elasticity of cell membranes,  Fb = 2





2

(C − C0 ) dA + G



GdA +



dA +

pdV,

(1)

where C = (c1 + c2 ) and G = c1 c2 are the total an Gaussian curvatures of the membrane surface given by principal curvatures c1 and c2 ,  and G are the bending rigidity and saddle-splay modulus, C0 is the so-called spontaneous curvature that accounts for any asymmetry in the membrane internal structure, whereas  and p generically represent a surface tension and a pressure difference across the membrane. In the Helfrich initial description, these two components are Lagrange multipliers to ensure that cell area and volume, respectively, are conserved. The curvature energy (1) is often expressed in the literature in terms of the mean curvature H = (1/2)(c1 + c2 ), though we follow here the total curvature

notation. It is remarkable that the integral of the Gaussian curvature over a surface is a topologic invariant, and consequently it only plays a role in the membrane elasticity in processes comprising topological transformations. For the case of closed membranes, such as cells, the Gaussian term remains constant and for simplicity it can be ignored. The minimization of (1) for an ellipsoidal shape under the appropriate values of area and volume leads to the biconcave discocyte of the RBC as the equilibrium shape. Ensuing studies investigated the properties and minimal shapes of the Helfrich energy and the theory has been refined to incorporate additional mechanisms such as the area-difference elasticity (Sheetz and Singer, 1974). Subsequently, the elastic contribution of the spectrin cytoskeleton which attaches to the lipid membrane was incorporated, with the aim of explaining the entire phenomenology of the human red blood cell (Evans, 1974; Iglic, 1997). The cytoskeleton provides resistance to in-plane deformations and plays a fundamental role in the cell response under certain deformations (such as morphological changes during crenation (Lázaro et al., 2013) or squeezing during optical tweezers experiments, Li et al., 2005), adding a shear-stretching contribution to the membrane elasticity. In our phase-field model we do not consider the cytoskeleton contribution and it therefore applies to problems in which this network is absent (such as in membranes of cell organelles, e.g. in the Golgi apparatus) or if it plays a subdominant role, as occurs during blood flow Forsyth et al. (2011). In the last years, the Helfrich model has been incorporated to different dynamic theories, offering the possibility of studying new and more complicate phenomena. Many of the results of this theory have proven in good agreement with experiments; nice examples include the theoretical prediction of shapes of the stomatocyte–echinocyte transition (Lim et al., 2002), the study of tubulation when polymers are attached to a lipid vesicle and effectively induce a spontaneous curvature (see Section 4) (Campelo and Hernández-Machado, 2008), or the experiments of stretching of red blood cells with optical tweezers (Li et al., 2005). In this section, we first discuss the microscopic basis of the Helfrich theory in order to show the consistency of the mesoscopic approach. We subsequently present the general theory of elasticity and its main results and applications to biological membranes. 2.1. Microscopic realization

fm =

1 kH 2

AH A0H

2 −1

+

1 kT 2

124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161

162

Any mesoscopic description assumes that one can define domains of small size compared to the length scale of the system size, so that the variables defined at the mesoscale (i.e. at each of these domains) capture the relevant properties of the microscale. With the objective of explaining this fundamental assumption, we present here the simple model proposed by Petrov and Bivas (1984) which, in spite of being highly non-realistic, is useful to naively illustrate the connection between both scales. The model assumes a rough description of the interactions between lipids, and from this simple basis the Helfrich free energy (1) for a bilayer is derived. The model assumes a harmonic approximation of the free energy per molecule,



123



AT A0T

163 164 165 166 167 168 169 170 171 172 173 174

2 −1

(2)

where AH/T are the areas per molecule of the head/tail, respectively, and kH/T are the harmonic constants related to the respective interactions between each group. A0 are the preferred areas, related with the equilibrium intermolecular distance in the relaxed monolayer. Still, the effective constants kH/T should be related with the specific bonds between molecules, but this is difficult to address at this simple level of description. If one defines the neutral surface as the point of the lipid where the forces are balanced, and we call A the

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

175

176 177 178 179 180 181 182 183

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx 184 185 186

187 188 189 190

191

192 193 194 195 196 197

lipid section at this point, the head and tail areas can be expanded in terms of the curvature at this neutral point as AH/T = A[1 + ıH/T C

fm

A

200 201 202 203 204 205 206 207 208 209 210 211 212 213 214

2

1 = Km −1 2 A0

(4)

kH /A0H

+ kT /A0T

where Km = is the stretching modulus per molecule. On the contrary, for a positively curved membrane, the head and tail areas are expanded and contracted, respectively, with respect to the neutral one. Thus, introducing (3) in (2), and considering all the molecules contained in the membrane, the energy of the deformation leads to the Helfrich energy (1), with the elastic parameters H

(

H

Gm m

where

T

+

Gm = (A0H − A0T )( C0m =

199

(3)

where ıH/T are the respective distances to the neutral surface. In the case of a flat interface (C = G = 0), when the areas per lipid at the head, tail and neutral points coincide, the energy per molecule reduces to

m = ı2 Km 198

+ ı2H/T G].

T) H

2

−

= ıH ıT Km , T)

H

(

H

+

T T)

, 2

(5)

1 . ıH − ıT

H/T

= kH/T /A0H/T , and the superscript m indicates that it

refers to one of the monolayers. Some properties of the elastic parameters can be inferred from these expressions. First, the bending parameter is always positive, as opposed to the Gaussian modulus which depends on the specific properties of the lipid. Analogously, the sign of the spontaneous curvature depends on the relative position of the neutral surface with respect to the head and chain. Thus far we have just considered a monoloyer. For a bilayer, the elastic parameters depend on whether the monolayers are connected or not. For simplicity, we suppose that there is no interaction and then m

 = 2 ,

G =

2(Gm

− m C0m N ),

c0 = 0.

(6)

where N is the distance between neutral surfaces. In spite of the enormous simplifications considered in this model, it shows how the Helfrich energy can be derived from a simple harmonic description of the lipid–lipid interactions.

3

where  and  are material parameters known as Lamé coefficients. Fel actually represents an excess energy with respect to the undeformed state. The interpretation of the diagonal terms of the deformation tensor u˛ˇ is simple. Considering a stretching deformation x = x (see Fig. 1A), a simple calculation shows that the non-diagonal terms of the strain tensor are zero, leading to an 2 energy of deformation Fel = (/2 + )[ − 1] V , where V is the volume of the underformed object. Conversely, for the shear deformation of Fig. 1B, given by x = x + ˛y and y = y + ˇx, the diagonal terms of the strain tensor vanish. Given that any planar deformation can be decomposed into a pure stretching and a pure shear, we rewrite the elastic energy, Fel =

 

 u˛ˇ −

1 ı u 3 ˛ˇ

2

E 2(1 + )

Fel =

 

u2˛ˇ +

228

u˛ˇ =

218 219 220 221 222 223 224 225 226

229 230

231

1 2



∂u˛ ∂uˇ + ∂xˇ ∂x˛



.

(7)

The energy of an elastic object subjected to a certain deformation, given by the tensor u˛ˇ , reads Fel =

   2



(u )2 + (u˛ˇ u˛ˇ ) dV,

(8)

234 235 236 237 238 239 240 241 242 243

244

245 246 247 248 249



u2 dV, 1 − 2 

(10)

250

251

2.3. The bending mode

252

Cell membranes have a typical thickness of 4 nm, whereas the overall cell length is usually 5–10 ␮m. In order to comprehensively understand the elasticity of membranes, we explore here the elastic properties of plane objects of small thickness compared to their surface. Thus, we represent the membrane as a generic object, such as a plate, subjected to a pure bending deformation. After averaging the elastic parameters over the plate section and adopt the twodimensional description, the Helfrich free energy is recovered. The stress tensor is given by the variations of the elastic energy with respect to the deformation,

˛ˇ (x) =

∂Fel . ∂u˛ˇ (x)

(11)

The specific dependence on the strain tensor (11) reads



1 u ı˛ˇ . 3

plate

F˛

=

∂ ˛ˇ . ∂xˇ

(12)

F˛ext =

∂ ˛ˇ dV = ∂xˇ

(13)



˛ˇ nˇ dA

(14)

Hence, the equation of force balance on the object surface, where the external forces apply, reads

˛ˇ nˇ =

F˛ext ,

254 255 256 257 258 259 260 261 262

263

265

266

The deformation of the object may respond to the action of an external forcing, F˛ext , on the object surface. This force is balanced by the internal stresses of the plate,



253

264



And the force exerted by the object can be derived from

227

217

(9)

233

where E = 9K/(3K + ) and = (3K − 2)/2(3K + ).

2.2. Theory of elasticity Membranes are elastic structures whose mechanical properties are usually described by the general theory of elasticity. The Helfrich free energy (1) with zero spontaneous curvature can be derived from this theory under certain assumptions. Although the strict derivation is cumbersome, a simplified approach is outlined here to facilitate the interpretation of the deformations and elastic modulus of the membrane in terms of its molecular structure. The deformation of any object is described by the deformation  − x , where x are the coordinates in the relaxed vector u˛ = x˛ ˛ ˛  represent the coordinates under deformation. For simstate and x˛ plicity we shall assume an isotropic, homogeneous material. The strain tensor is defined as

216

+ K(u )2 dV,

where  is renamed as shear modulus and K =  + (2/3) is the bulk (or compression) modulus. Within the frame of membranes, often the elasticity is described in terms of two new modulii, the Young modulus E and the Poisson ratio , and the expression for the elastic energy now reads

˛ˇ = K(u ı˛ˇ ) + 2 u˛ˇ − 215



232

(15)

where n is the normal vector to the object surface. Let us consider a flat plate of thickness h. Consider now the bending deformation shown in Fig. 1C, with no in-plane deformation ux = uy = 0 but outof-plane displacement uz = (x, y) in the Monge representation (do Carmo, 1976). We suppose that, given that the plate is very thin, the external forces required to bend it are small compared to the internal tensions across the membrane, so that Fext are neglected in

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

267

268 269 270

271

272 273 274 275 276 277 278 279 280 281

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

4

Fig. 1. Scheme of the main deformations. Bold lines represent the underformed state and dashed lines the deformation. (A) Stretching, when the object increases its area; (B) shear, when the object deforms without a change on its area; and (C) bending, which represents a normal displacement  to the surface.

282 283 284 285 286 287 288 289

(15). Note that this hypothesis is a priori not obvious, but still let us assume its validity under certain conditions (Landau and Lifshitz, 1999). For a strictly flat plate, nz = 1 and therefore xz = yz = zz = 0. If the plate is only slightly bent, these components are strictly nonzero but remain small compared to the rest of components of ˛ˇ , so we can equate them to zero and use the resulting equations to solve the strain tensor components. For this deformation, all the terms of the strain tensor can be calculated explicitly: uxx =

2 −z ∂x ,

312

2

290

uzz =

2 −( /(1 − ))z(∂x 

(16)

2

uxy = uyx = −z ∂xy , uxz = uzx = 0, uyz = uzy = 0. 291 292

By introducing these expressions in (8), the elastic energy of the bent plate is

293

 294

E z 1+

Fel,b =

 +

295

1 2(1 − )

2

2

∂  ∂x∂y

2



2

−

2

2

∂  ∂  + ∂ x 2 ∂ y2

2

2

∂ ∂  ∂ x 2 ∂ y2

dV.

(17)

296 297 298

Integrating across the plate thickness, from −h/2 to h/2, in the z-direction,

299

300

Fel,b

h3 E = 24(1 − 2 )

 

 301

+2(1 − )

2

2

2

∂  ∂  + ∂ x 2 ∂ y2

∂  ∂x∂y

2

2

−

2

2

∂ ∂  ∂ x 2 ∂ y2

 dS.

(18)

302 303 304 305 306 307 308

Some considerations are required to readily identify the dependent terms in (18). In the Monge parametrization, the normal vector to the surface can be obtained from nˆ = ∇ /|∇ |. The total and Gaussian curvatures are defined in terms of the normal vec2 ˆ − ∇ s · nˆ : ∇ s · n], ˆ respectively. tor as C = −∇ · nˆ and G = 1/2[(∇ s · n) Hence, straightforward calculations lead to (do Carmo, 1976)

309

2

310

2

2

2

C = −[(1 + (∂x ) )∂y  + (1 + (∂y ) )∂x ] − 2(∂x )(∂y )(∂xy ) 2

311

2

2

≈ −[∂x  + ∂y ],

  

2

(19)

2

2

2

2

≈ (∂x )(∂y ) − (∂xy ) ,

(20)

313

314 315 316 317 318



C 2 + G G dS,

(21)

where we have introduced the bending rigidity  = Eh3 /12(1 − 2 ) and the saddle-splay modulus G = − (1 − ). This derivation applies generically for objects such as plates or, in the present context, monolayers. However, a bilayer consists of two monolayers which glide each other, so that the bending rigidity of the bilayer, bil , arises as the sum of the rigidities of the monolayers, m . Considering monolayers of thickness h = d/2, the bilayer bending modulus is given by (Deserno, 2014) bil = 1m + 2m =



2 2 2

2 2

for the total (19) and Gaussian (20) curvatures. We have assumed that the out-of-plane displacement  is small relative to the length scale of the plate, so that the gradients are small ∂˛   1. By direct comparison of (18) with expressions (19) and (20), the elastic energy for the bending deformation can be rewritten as Fel,b =

2 + ∂y ),

2

(1 + ((∂x ) ) + ((∂y ) ) )

2

uyy = −z ∂y ,

2

(∂x )(∂y ) − (∂xy )

G=

2

d3 E KA , = 24(1 + ) 48(1 − 2 )

(22)

where we have introduced the area-compression modulus KA = Ed/2(1 + ), which represents the energetic cost of expanding/compressing the area of the plate, and it is related with the volumetric bulk modulus K. Hence, assuming a symmetric plate (C0 = 0) and considering a pure bending deformation, the general elastic energy (8) reduces to the bending contribution (21). The particular elastic modulii used to describe cell mechanics depend on the cell species. In the frame of RBC elasticity, it is usual to characterize the membrane mechanics with the bending and shear modulii. On the contrary, in cells with an inner structure (such as leucocytes or epithelial cells), especially in the framework of atomic force microscope experiments, most studies focus on the Young modulus (Kuznetsova et al., 2007). We hereafter concentrate on the elastic description based on the bending modulus, which is widely extended in the membrane mechanics field. The RBC membrane bending rigidity has been measured by different experimental techniques, for instance flicker spectroscopy and phase contrast microscopy, analyzing the membrane amplitude oscillations due to thermal fluctuations, or disturbing techniques such as micropipette aspiration and by optical tweezers (see Seifert and Lipowsky (1995), and references therein). Typical values fall between 10 and 50kB T with slight deviations depending on the specific technique. The area-compression modulus of flat bilayers has been measured for different lipid species, with typical values of 100–250 mJ/m2 (Evans and Needham, 1987; Rawicz et al., 2000). Taking into account these values for the bending rigidity and area-compression modulus, it can be shown that bilayers

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

319

320 321 322 323 324 325 326 327

328

329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx 356 357 358 359

360

present a high incompressibility. Considering the stretching energy Fst = (1/2)KA (A − A0 )2 /A0 , where A0 is the membrane area in the relaxed state, variations in the membrane area are given by Seifert and Lipowsky (1995) 1 A − A0 = A0 KA



∂Fst ∂A



.

(23)

385

Thus, substituting the stretching energy by the bending energy (21), the variations in area induced by bending deformations are found to be disregardable, /KA R2 ∼ 3 ×10−8 , for the given values  = 2.0 × 10−19 J, KA = 10−1 J/m2 and cell radius R = 8 ×10−6 m. This huge energetic penalty imposed to the expansion of area implies that any membrane deformation driven by bending is required to effectively maintain a constant area. For practical purposes, it is usual to remove the elastic contribution of the area-compression and strictly impose the constraint of constant area, for instance by introducing a Lagrange multiplier, as proposed by Canham and Helfrich and shown in (1). The molecular basis of this phenomenon is found in the strong hydrophobicity of the lipid tails. An expansion of the membrane area implies that the area per molecule increases from its relaxed state and the lipid tails are therefore exposed to the water molecules. To avoid this situation, the attractive interaction between molecules enforces the bilayer incompressibility. Due to the fluidic nature of the bilayer, the shear modulus is, by definition, zero. Finally, it must be considered that cells contain specific regulatory systems to maintain their enclosed volume constant, based on ion gradients which control the osmotic pressure (Helfrich, 1973). Accordingly, another Lagrange multiplier, accounting for this volume conservation, is incorporated to the elastic membrane energy. The addition of the bilayer incompressibility and constant enclosed volume to the bending contributions of (21) allow to recover the complete Helfrich free energy (1).

386

3. Phase-field models

361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384

387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416

The Helfrich theory provides a basis for the elastic characterization of cell membranes. Still, a geometric description is needed to deal with the membrane morphology. Minimal shapes of model vesicles can be found by different analytic parameterizations (Deuling and Helfrich, 1977; Seifert et al., 1991), but the study of more complex problems, including any kind of dynamic phenomena, usually requires of a more flexible framework. Given that the representation of the membrane as a two-dimensional layer is reasonably accurate, the simplest and most direct formulation consists in defining a mesh of points which represents the membrane neutral surface, and from here extract the local mean curvature or deformation tensor necessary to compute the elastic energy. This must be combined with a minimization procedure (such as a Monte Carlo free minimization, Lim et al., 2002) to obtain equilibrium shapes, or introducing the elastic mesh in a dynamic theory if one is interested in membrane kinetics. Most important examples include the immersed boundary methods (Peskin, 2002; Kaoui et al., 2012), integral boundary methods (Pozrikidis, 1992, 1995) or multiparticle collision dynamics (Malevanets and Kapral, 1999; McWhirter et al., 2008). Methods in this direction have been successfully applied to the study of many membrane related topics (Li et al., 2005; Peng et al., 2013). However, these methods require a explicit tracking of the membrane position and the calculation of the deformation variables, i.e. the curvature. This can potentially cause problems related to the need of a good resolution in the mesh grid. For instance, the curvature calculation in a discrete tessellation is an old problem in discrete differential geometry (Gatzke and Grimm, 2006). Even the most accurate techniques fail to achieve a good definition, and the problem must be generally simplified to avoid large numerical calculations. A different approach, based on

5

an Eulerian rather than a Lagrangian description, are the phasefield models. The membrane is identified from an auxiliary scalar field defined in the entire space, and the method details the dynamics of the field, rather than specifically dealing with the membrane position. This formulation also avoids the problematic of defining the boundary conditions at the membrane surface. The application of phase-field methods to amphiphilic systems was extensively investigated in the past (Gompper and Schick, 1994), although it has not been until recently that these models have been used in the study of cell morphology and dynamic response. One of the main advantages of the phase-field model is that the evolution and shape of the membrane does not need to be tracked, as in the explicit methods, but it spontaneously evolves with the phase field dynamics. Phase-field methods also invite to a deep analytical exploration as they have a robust physical basis. 3.1. Ginzburg–Landau free energy

419 420 421 422 423 424 425 426 427 428 429 430 431

433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449



=

418

432

Phase field models have been applied in the study of a wide different phenomena of phase transitions and interfacial problems, such as superconductivity and solidification (Steinbach, 2009). The origins of phase-field models are found in the mean-field approach to phase transitions. The Landau theory consists in a free energy which is expanded in powers of a scalar field, called order parameter , which receives different interpretations depending on the system, as we will see below, but, for instance, in a ferromagnetic system it is readily identified as the magnetization. The symmetries of the system specify the value of the coefficients of the expansion. However, this model does not account for the presence of interfaces which could have an energetic cost associated with the interaction between the components of each phase. Ginzburg and Landau generalized this expression in their studies about superconductors and they incorporate to the free energy of the system, F[ ], powers of the gradients of the order parameter,

F[ ] =

417

L( , ∇ , ∇ 2 )dV

 

450

2

f ( ) + g( )(∇ ) + c(∇ 2 )

2

 dV.

(24)

where L is the energy density, f( ) is the bulk potential so that in equilibrium f ( eq ) = 0 and eq are the stable bulk phases, and g( ) and c represent the energetic cost of having an interface (we will go in deep on this below). Typically, f( ) is chosen to form a double-well potential, so that two phases are present. The interface is characterized by a smooth profile characterized by a width , and it is usual to fix the interface position at the isosurface = 0. The method does not necessarily need that operates on the same scale of the width of the real interface, which in some systems can be the atomic scale, but it is sufficiently to require that (i) is much smaller than any other length of the system and (ii) the interface incorporates the relevant information from the microscale trough the effective constants of the model. The construction of a Ginzburg–Landau free energy often responds to a purely phenomenological basis and it generally attends to the symmetries of the system. The order parameter is related to a characteristic physical property depending on the specific system; in the particular case of a membrane, it can be associated with the concentration of lipids, (x) = 1 ∓ lip (x)/0 , where 0 is the lipid density of maximum packing and the ∓ refers to the inner or outer regions of the membrane. Since lipids are indissoluble in water, lip = 0 in the fluid bulk that surrounds the membrane, and therefore = ±1 in these outer and inner regions, whereas it achieves the value lip = 0 in the membrane, giving rise to = 0. It is noteworthy that is defined in a spatial scale sufficiently large

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

451

452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

6 477 478 479

480

to have a statistically good average of the molecular densities, but yet small enough to capture the mesoscopic spatial variations of the density.



3.2. Thermodynamics of phase-field models

ıF = −

485

In thermodynamics, the chemical potential of a species in a mixture is defined as the free energy change for deviations in the concentration from the reference value. Analogously, within the phase-field framework the chemical potential is obtained from the free energy (24),

486

[ ] =

481 482 483 484

Introducing expressions (28) in (29), and after several straightforward integrations of those terms containing second and third gradients of ıx˛ , one finds

ıF ∂L ∂L ∂L . = − ∇ˇ + ∇2 ı ∂ ∂(∇ ˇ ) ∂(∇ 2 )

(25)

523 524 525



∂L ∂L ∂L ∇˛ + ∇˛∇ˇ + ∇ ˛ ∇ 2 ıx˛ dV ∂ ∂(∇ ˇ ) ∂(∇ 2 )

526

∂L ∂L ∂L ∇ ˛ ıx˛ dV − ∇ˇ + ∇2 ∂ ∂(∇ ˇ ) ∂(∇ 2 )

527

 −

 −

522



∂L ∂L ∂L ∇˛ + ∇˛∇ˇ − ∇ˇ ∇˛ 2 ∂∇ ˇ ∂∇ ∂∇ 2

× ∇ ˇ ıx˛ dV.

528

(30)

529 530

487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505

In equilibrium, the chemical potential is constant and in addition to the boundary conditions in the bulk phase (x → ± ∞) = ±1, this equation can be solved obtaining the equilibrium profile of the interface. Typically, the relaxed profile is a smooth tanh-like function with interface width . As we will see later on, the chemical potential is a relevant variable as it governs the diffusion of the order parameter (Section 3.6). In the presence of a non-uniform phase field, there is a thermodynamic local force at each point (e.g. this force can be naively understood in the case in which is a concentration, and then gradients of concentration force diffusive fluxes toward a uniform concentration). Given that in the bulk  = const ., the force density is localized at the interface. Although this elastic force can be derived by different procedures, here we follow up the elastic description of Section 2.2. The stress tensor is an important magnitude to characterize the deformation of membranes, though it is not obvious how to compute this tensor from a phase-field free energy of the form (24). The work required to deform an object by a small displacement ıx˛ is given by (Landau and Lifshitz, 1999)

 506

ıF = −

 ∇ ˇ ˛ˇ ıx˛ dV =

˛ˇ ∇ ˇ ıx˛ dV.

(26)

510

In the phase-field framework, we need to specify the variations in the order parameter due to the deformation ıx. Assuming that these small variations only correspond to convective fluxes (Brannick et al., 2014),

511

∂t + ∇ · ( v) = 0,

507 508 509

512 513

(27)

and writing v˛ = ıx˛ /ıt, we obtain the variation of the order parameter and, after differentiation, its derivatives ı = − ∇ ˛ ıx˛ − ∇ ˛ ıx˛ . ı∇ ˇ = −∇ ˇ ∇ ˛ ıx˛ − ∇ ˇ ∇ ˛ ıx˛ − ∇ ˇ ∇ ˛ ıx˛

514

−∇ ˛ ∇ ˇ ıx˛ .

(28)

ı∇ 2 = −∇ 2 ∇ ˛ ıx˛ − 2∇ ˇ ∇ ˇ ∇ ˛ ıx˛ − ∇ ˛ ∇ 2 ıx˛

516 517 518

From these expressions we express the change induced in the order parameter when the interface is perturbed with the general deformation ıx˛ . The work necessary to induce such a deformation of an interface with free energy (24) reads

519

 

 520 521

ıF =

ıLdV =



˛ˇ =

L− −

ıL ı



∂L ∂L ∂L ı + ı∇ ˇ + ı∇ 2 dV. ∂ ∂∇ ˇ ∂∇ 2 (29)

531 532 533 534 535 536

 ı˛ˇ −

∂L ∂L ∇˛ + ∇ˇ ∇˛ ∂(∇ ˇ ) ∂(∇ 2 )

∂L ∇ ˛ ∇ ˇ . ∂(∇ 2 )

537

(31)

538 539

From (31), it can be shown that the divergence of the stress tensor reduces to

∇ ˇ ˛ˇ = − ∇ ˛ ,

(32)

which provides an expression for the local force density of the interface in terms of the chemical potential. The expressions obtained so far are valid for any free energy of the form (24), so that the elastic properties of the interface are solely determined by this free energy. The dependence of the force density on the chemical potential can be understood from a thermodynamic perspective, taking into account that the stress represents the internal reaction to an external pressure (strictly speaking, the pressure and stress tensor are related by P˛ˇ = − ˛ˇ ). Hence, Eq. (32) may be related with the Gibbs–Duhem equation, which reads VdP =



Ni di

(33)

540 541 542 543 544 545 546 547 548 549 550 551 552 553

554

i

where Ni is the amount of matter of the species i and taking into account that ∼ N/V, leads to dP = d. The force density arises as the free energy change per unit volume, ı, due to the transport of matter concentration for a change in the chemical potential ı. 3.3. Elastic properties of phase-field interfaces

−∇ ˛ ∇ 2 ıx˛ − 2∇ ˇ ∇ ˛ ∇ ˇ ıx˛ − ∇ ˛ ∇ 2 ıx˛ . 515

In the first term of the right hand side of (30) the gradient ∇ ˛ L can be recognized. Analogously, the second term contains the expression of the functional derivative of L, given by (25). Hence, after integrating by parts the first term, and comparing with (26), the stress tensor is identified,

In the Helfrich theory, the physical meaning of each elastic modulii is clear, as each modulus is explicitly associated to the elastic deformation that penalizes (i.e. in (24), the surface tension represents the energetic cost of having an interface of surface dA, and the bending modulus corresponds to the cost of an interfacial bending given by the curvature C). In the phase-field representation of the interface free energy, however, the information about the geometrical properties of the interface deformation is implicitly contained in the gradients of the order parameter, and the elastic properties of the interface cannot be directly identified. It is therefore necessary to establish a connection between the phase-field coefficients,

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

555 556 557 558 559

560

561 562 563 564 565 566 567 568 569 570 571

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592

593

594 595 596

g and c, of the phase-field free energy (24), and the interface elastic parameters expressed in the Helfrich energy (1). The comparison between both free energies for simple geometries, such as spheres or cylinders, is a useful method for obtaining a mapping between the Helfrich and phase-field representations. In this Section, we derive the expressions of the elastic coefficients of (1) from the parameters of the free energy (24), as first proposed by Gompper and Zschocke (1991), and Gompper and Zsckocke (1992). Thereby, the method allows to specify the elastic properties of the phasefield free energy and conciliate this expression with the classic elastic description. Let us suppose a generic interface with free energy (1), with contributions of surface tension and bending, but disregard the pressure difference term assuming that the interfacial surface is open. Note that this free energy can be locally understood as an expansion in terms of the radius of curvature 1/R, since C ∼ 1/R and G ∼ 1/R2 . Hence, the surface tension is associated to the zeroth order, 1/R0 ; the spontaneous curvature corresponds to the first order of the expansion, 1/R1 ; and the bending and saddle splay modulii are associated with the second order, 1/R2 . The elastic energy per area (1) for a sphere of radius R is given by (s) Fe

A

=



+

  2 2

c0

598 599 600 601 602

603

604 605

606

2 2 + G − . c0 + R R2

607 608 609 610

Fe   =  − c0 + . A R 2R2

612

(36)

where g = ∂g/∂ and f = ∂f/∂ . Multipliying by ∂z 0 , and integrating



3

f ( 0 ) − g( 0 )( 0 ) + 2c ∂ 0 ∂ 0 −

1 2 (∂ 0 ) 2



= const.

(37)

Note that, formally, (36) is the stationary condition of the Euler–Lagrange equation and then (37) is its first integral. For the specific geometry of spheres and cylinders, the free energy (24) reads F

A

[ ]

=

1





∞

drr d f ( R ) + g( R )(∂r R )

Rd+1

0 2

+c ∂r R +

d ∂r R r

2

2

,

(38)

619

where d = 1 and 2 for cylinders and spheres, respectively, and R is the phase-field profile for these two configurations. If the radius R is large, the phase-field profile R can be expressed in terms of the relaxed planar solution, 0 (z), by identifying the position of the interface at z = r − R. This approximation, known as locally flat interface, allows to expand the profile in powers of 1/R,

620

R (r − R) = 0 (r − R) +

614 615 616 617 618

621 622

1 (r − R) + ··· R

dr 1 + 0



2

+c ∂r 0 +



R

2

2g( 0 )(∂r 0 ) + 3c(∂r 0 )

d ∂r 0 r+R

∞

×

2

−

2 625

1 ı(d − 2)2c R2

2

dr(∂r 0 ) .

626

(40)

627

0 628

The expansion of this expression in powers of 1/R leads to F(c,s) [ ] = A

629 630



∞

2

2

2

dr[2g( 0 )(∂r 0 ) + 4c(∂r 0 ) ]

631

0



∞

+

dr 0



∞

+ 0

2 2r 2 2 [2g( 0 )(∂r 0 ) + 4c(∂r 0 ) ] R

632

2 r2 2 2 [2g( 0 )(∂r 0 ) + 4c(∂r 0 ) ] R2



∞

dr 0

2c 2 (∂r 0 ) . R2

633

(41)

Comparing this expression with the free energy of cylinders (35) and spheres (34), and identifying the coefficients of the corresponding terms of the expansion, one obtains



634

=



(39)

By introducing this expression in (37), and using this to remove f( R ) from (38), to leading order the free energy for spheres and membranes now reads,

637 638

s(z)dz, +∞

−c0 =



636

+∞

−∞

ıF[ ] 2 2 4 = f  − 2g( 0 )∂z 0 − g  (∂z 0 ) + 2c ∂z 0 = 0, ı



613

(35)

The inclusion of the cylinder is important to identify the bending modulus, which appears coupled to the saddle-splay modulus in the spherical geometry. The equilibrium condition for the free energy (24) determines the relaxed planar interface profile 0 (z) in the normal direction z,

(c,s)

=

  r d

635

For a For simplicity, the first term is rewritten  =  cylinder, the Gaussian curvature vanishes and then the free energy reduces to

611

[ ]

A



∞

(34) + c02 /2.

2

F

623 624



+ ı(d − 1)

(c)

597

(c,s)

7

zs(z)dz,

−∞ +∞

=

(42)

639

2

2c(∂z 0 ) dz, −∞



+∞

z 2 s(z)dz,

G + 2 = −∞

where we have recovered the planar profile notation as a function of the normal coordinate z and we have introduced the function s(z) = 2g( 0 )(∂z 0 )

2

2 2 + 4c(∂z 0 ) .

(43)

The physical meaning of this expression is found within the elasticity theory framework. Helfrich (1981) established that the surface tension, spontaneous curvature and Gaussian bending modulus can be obtained as the zero-th, first and second moments of the lateral stress profile across the membrane, respectively. Thereby, the functional form of (42) suggests to interpret s(z) as the stress profile, though it is remarkable that in the expression for the second moment, G and  are coupled, a feature not captured in the Helfrich theory. The particular microstructure and chemical composition of the interface dictate the interactions and internal tensions of the interface, ultimately determining its elastic behavior. For instance, in the case of two immiscible fluids, such as the water–oil coexistence, the unique contribution to the free energy corresponds to the surface tension, associated solely to the cost of having a surface in which both species of molecules interact. The benefit of homotypic interactions when a molecule is surrounded by others of the same species is lost in the interface, where the van de Waals interactions between oil and water molecules are weaker. Accordingly, the system tries to minimize this contact surface, but the bending of the surface does not have any energetic cost. Other systems, such as monolayers formed by amphiphilics, present a more complicated interface internal physics, reflected in an extremely rich phenomenology. In the presence of water and oil, amphiphiles

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

640 641

642

643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

8 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693

form several different structures ranging from micelles to lamellar (plane) arrays of amphiphiles separating regions of each fluid (Gompper and Schick, 1994). The principles driving this kind of organization respond to the polar nature of the amphiphile, with a hydrophilic head, which prefers to interact with water, and the hydrophobic tail, which prefers the oil. The system evolves to minimize the contact between the hydrophobic tail and water. In these systems, the lateral tensions between the amphiphiles induce a more complex elastic behavior, and all the terms in (1) contribute to the interface elasticity. Still, in the low curvature regime (when the interface thickness  is much smaller than the typical radius of curvature, /R  1; this regime applies to most systems, since  is of the order of the lipid length), the dynamics of the interface is generally dominated by the surface tension term, since it corresponds to the leading term of (1). Lipid bilayers present an even more complex internal structure than monolayers, and are known to induce a strong reduction in the surface tension of the interface, sometimes of up to 5 orders of magnitude (Gompper and Schick, 1994), as a result of the internal balance of the lateral stresses between lipids. Accordingly, the dynamics of the membrane is driven by bending (the subsequent non-vanishing term in the curvature expansion (1), assuming membrane symmetry). Since the scope of this article is the modeling of biological lipid membranes, hereafter we focus on tensionless interfaces. It is noteworthy that the elimination of the surface tension contribution in phase-field models is highly non-trivial, as it requires of a subtle balance of the lateral tensions of the interface.

where the subscript b indicates that this model corresponds to a bending free energy. The free energy can be rewritten as Fb [ ] =

∗



2

( [ ]) dV,

2

(45)

where we have introduced the functional = − + 3 − 2 ∇ 2 . From this expression the chemical potential as defined in (25) reads ∗

2

b = ıFb [ ]/ı =  [(3 − 1)

− ∇ 2

2

].

(46)

The relaxed profile is obtained by solving the equilibrium condition, b = 0. Although this equation does not have a unique solution, the trivial one = 0 represents the minimal energy solution, given that the bending energy is always positive. Other potential solutions may arise as metastable solutions. This equilibrium profile of the model √ (44), = 0, can be analytically integrated, leading to 0 = tanh(z/ 2 ), connecting the bulk phase -1 with the phase +1. If this expression is introduced in (42) with the energy coefficients√from (44), the resulting elastic parameters are  = c0 = 0,  = (2 3 /3 2)∗ and G = 0, consistent with a symmetric membrane. Hence, the model describes the cell as two fluid domains separated by an interface with the elastic properties characteristic of membranes. From the bending free energy (44) we can study the interface properties that determine its elastic response. The stress tensor is computed from (31), obtaining b

˛ˇ

=−

 ∗  2

726 727

728

729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747

{[− + 6 − 5 − 2 (3 − )∇ 2

4

6

2

2

3

2

748

2

− (3 2 − 1) 2 (∇ ) − 6 2 2 (∇ ) + 2 2 ∇ 2 ∇ 2

749

2

694

− 2 (∇ 2 ) ]ı˛ˇ + 2 2 (3 2 − 1)∂˛ ∂ˇ

3.4. Cell membrane model

695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722

The theory presented thus far offers the necessary ingredients to build a model for cell membranes. The simplest approach to describe the cell is to follow the Ginzburg–Landau spirit, considering two domains, the inner fluid (e.g. cytosol) and the outer aqueus environment (e.g. plasma), and associate each medium to one of the equilibrium phases of the order parameter. The coefficients associated to the order parameter gradients determine the elastic properties of the interface, and hence they must be chosen to capture the two main characteristics of membranes: resistance to bend and vanishing surface tension. The first condition is achieved by c > 0, given that a second derivative is necessary to introduce a bending contribution. The second condition requires a negative constant value of g, or alternatively an inhomogeneous function g = g( ). However, the choice of the free energy coefficients to obtain a vanishing surface tension is delicate. The coefficients determine the equilibrium profile of the order parameter, and in general its solution requires of numerical integrations. The profile enters in the elastic parameters calculation, (42). Thus, the elastic properties depend on the coefficients via the equilibrium profile which in turn can be very sensitive to the values of the coefficients. A priori, there is not a unique solution for g and c that produces the prescribed elastic properties of the membrane. We consider the particular case fb ( ) = 2 − 2 4 + 6 , gb ( ) = 2 2 (3 2 − 1), and cb = 4 , which is of particular interest because its equilibrium equation has analytical solution and therefore it allows a more fine control of the interfacial behavior. The interface free energy reads (Du et al., 2004; Campelo and Hernández-Machado, 2006)

723

724 725

∗ Fb [ ] = 2

 2

2

( − 2 + + (3 − 1) (∇ ) + (∇ ) )dV, 2

4

6

2

2

4

2

(44)

750

− 2 ∂˛ ∂ ∂ˇ ∂ − 2 [∂˛ ∂ˇ ∇ + ∂ˇ ∂˛ ∇ ] 4

4

2

+ 2 4 ∂˛ ∂ˇ ∂ ∂ + 4 4 ∂˛ ∂ˇ ∇ 2 }.

2

751

(47)

752 753

In Fig. 3C, the lateral stress profile sb (z) for the membrane free energy (44) is shown, in addition to the stress profile st (z) for the classic surface tension interface which applies, for instance, for a water–oil phase separation, given by the coefficients c = 0 (not bending penalization) and g( ) = g0 > 0. In this last case, the model includes an energetic cost for the existence of the interface. This area penalty is the surface tension. The stress through the interface are represented by a unique positive term, implying that the pressure is negative and thus the interfacial molecules are compressed, trying to minimize the surface area. In the case of the membrane, the lateral stress of the phase-field model does not capture a realistic profile, such as the one shown in Fig. 3B, but it condensates the information in two terms of repulsion and attraction. One could interpret the forces at the middle of the interface as the entropic repulsion between tails, whereas the two lateral attractive regions may correspond to the cohesive stresses acting at the hydrophilic–hydrophobic transition region. In general terms, the profile recovers the behavior proposed in the simple microscopic model studied in Section 2.1. In this approach, the balance between attractive and repulsive interactions means that the interface is not energetically penalized for having a certain surface (i.e., it has zero surface tension) but instead it is penalized for having a non zero curvature, as the first non-vanishing term in the expanded free-energy corresponds to 1/R2 . The parameter that determines this energetic cost is identified as the bending rigidity of the membrane. The membrane model should include both the bending contribution and the compression of the membrane. The areacompression effect is modeled by directly imposing a constant surface area, as first proposed by Helfrich. This is done by adding

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

Fig. 3. (A) Scheme of a lipid bilayer. The multiple interactions between the different chemical groups of the bilayer (including repulsion between polar groups, hydrophobic attraction, cohesive, repulsive effects between tails, etc.) leads to a complicate stress profile in real cell membranes. However, they are characterized by a balance between internal tensions, so that the surface tension vanishes in this system. (B) Example of lateral stress profile sexp (z) for a lipid bilayer, based on the results of Hu et al. (2013) obtained from MARTINI simulations of DMPC bilayers. Q4 (C) Lateral stress profile of the membrane phase field model, sb (z) (blue line), and for a tension interface, st (z) (red line). The phase-field model for membranes does not reproduce the exact lateral stress profile of a realistic membrane but effectively concentrates the interactions in two contributions, a central term of repulsion and two symmetric attractions, and their balance recovers the tensionless nature of the membrane. A physical interpretation of this simple profile is that the central term corresponds to the entropic repulsion between the lipid tails, whereas the lateral attractions would represent the attractions between head and tails, following the spirit of the model of Petrov and Bivas (1984). In comparison, the lateral stress profile of a tension interface, involves a single term at the frontier between the two phases, penalizing the presence of the interface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

784 785

786

787 788

789

790 791

792

793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817

a Lagrange multiplier in the bending free-energy, so that the complete membrane free energy is 3 Fmem [ ] = Fb [ ] + √ 2 2



2

A (∇ ) dV.

 x  √

2

(49)

and the volume integral converges to a surface integral over the interface, given by the coordinates x , 3 √ 2 2



2

(∇ 0 ) dV −→ → 0

obtained, as shown in Fig. 2 B. If r < 0 the interface is stable and the system will decay to the initial planar interface. For r > 0 the system is unstable. Note that in these calculations an explicit time dependence is required, though it has not been presented yet. The time evolution is discussed in Section 3.6 and we do not make here further comments on this, but let us concentrate on the stability of the system for each equilibrium phase. The phases ±1 are both stable for the entire range of wavelengths of the perturbation,  = 2/q. The phase 0 = 0 is, however, stable for a short range  < . Thus, eq = 0 is not macroscopically stable since it can only be present as a microemulsion, and thus being irrelevant in our membrane description. 3.5. Generalized membrane model

= ı(x),



Fig. 2. (A) Bulk potential of the curvature free energy (44), with three phases of equal energy. (B) Linear dispersion of the decay of a phase when subjected to a weak sinusoidal perturbation. The ±1 phases are stable, whereas the 0 phase is only stable for small domains of length < .



ı(x − x )dV =

818 819 820 821 822 823 824 825 826 827 828 829

(48)

where  A represents the Lagrange multiplier for the area conservation. If the interface is small, then (∇ 0 )2 behaves as a ı-function, 3 lim √ sech4

→0 4 2

9

dS.

(50)

Although the last term in (48) has the expression of a surface tension, note that  A is shape dependent and it varies with the deformation. The expressions of the chemical potential mem and stress tensor mem for the complete membrane model must be obtained. Finally, it should be noted that a term 6 is present in (44). The system actually presents three equilibrium phases, eq = 0, ± 1, evoking the case of amphiphilic systems, where often a third phase is introduced accounting for the presence of the lipid-rich domain (Gompper and Schick, 1994). Although the three phase description is actually more realistic (since it considers explicitly the presence of lipids), the two phase model described above is more consistent with the Helfrich theory, in the sense that it treats the membrane as a sheet characterized by its local curvature, disregarding the microstructure of the membrane. By means of a linear stability analysis, it can be shown that the phase eq = 0 is not macroscopically stable. Fig. 2 A shows the form of the bulk potential fb , with three equilibrium phases. The stability analysis reveals that the phase eq = 0 presents a different behavior with respect to eq = ±1. The exact derivation requires of a detailed description, but we briefly outline the main steps here. Let us suppose a planar interface separating two of the stable phases. A sinusoidal small perturbation is introduced, 0 = eq + eiqx , and one assumes that it will decay as = 0 ert , an ansatz generally valid in the linear regime q  1. Introducing these expressions in the model, the growth rate r(q) is

830

The model presented in this section corresponds to a symmetric homogeneous membrane. However, the phase-field methodology allows to extend the model to cover more complex membranes. Two cases are considered here: asymmetric membranes with a non zero spontaneous curvature that affects to the membrane balance, and multicomponent membranes in which the elastic properties vary along the membrane surface due to an inhomogeneous composition. 3.5.1. Asymmetric membranes Cells often present asymmetric membrane composition, and benefit from the control of this property in a number of ways, as will be discussed for the case of pearling and tubulation in Section 4. The phase-field free energy (44) can be modified in order to account for a membrane asymmetry, effectively captured in the Helfrich model by a positive spontaneous curvature, C0 . The resultant free-energy reads (Campelo and Hernández-Machado, 2007b) FSC =

∗ 2

 



831 832 833 834 835 836 837 838

839 840 841 842 843 844 845 846

2

(− + 3 − ∇ 2 ) − c0 (1 − 2 )

dV.

(51)

The term (1 − 2 ) represents a ı-function centered at the interface position, so that the interface is forced to accommodate its √ surface to the spontaneous curvature C0 = c0 / 2. In principle, this spontaneous curvature can be nonuniform and present local domains in the membrane with variable asymmetry. 3.5.2. Multicomponent vesicles The different lipid species present in the membrane often pack forming monospecific aggregates, with important implications in the cell functioning as they are related with membrane trafficking and signaling (Simons and Vaz, 2004). The experiments carried out by Baumgart et al. (2003) with lipidic vesicles showed that the bilayer lipid composition can be controlled to form these

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

847

848 849 850 851 852

853 854 855 856 857 858 859

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

10

866

aggregates, and modulate the vesicle shape. Depending on the specific lipids of the domain, the vesicle presents inhomogeneous membrane properties. The morphology assumed by the vesicle responds to the balance of this particular membrane composition. The problem was theoretically approached by Wang and Du (2008), who developed a model in which both the bending rigidity and the spontaneous curvature depend on the local domain,

867

Fmc =

860 861 862 863 864 865



˝

i (C − C0,i )2 dA + 2



i dl,

3.6. Membrane dynamics and hydrodynamic coupling

∂˝

883

884

Fmc

869 870 871 872 873 874 875 876 877 878 879 880 881 882



=

2 ()  (− + 3 − ∇ 2 ) − c0 ()(1 − 2 ) dV 2

+

l( , )dV,

(53)

888

and added the usual constraints of constant total area and volume for the order parameter . The functional l( , ) accounts for the line element of the frontiers of the lipid domains,

889

l( , ) = L|∇ | |∇ | ,

886 887

890 891 892 893 894

895 896 897 898 899 900 901 902 903 904 905 906

907

908 909 910 911 912 913 914

2

2

(54)

with normalization constant L. Minimization of (53) leads to the minimal shapes of the vesicle. The model nicely reproduces the shapes obtained in Baumgart et al. (2003), and highlights the subtle control of membrane shape mediated by the formation of local aggregates with suitable elastic properties. 3.5.3. Alternative membrane models The model presented here exploits the balance of stresses along the smooth interface in order to obtain a tensionless membrane, whose only elastic contribution is then the bending. The goal of this approach is that the information of the interface geometry (e.g. the local curvature) is intrinsically contained in the gradients of the order parameter, so that we only need to specify the dynamics of the order parameter and forget any other treatment of the surface membrane. A different approach has been proposed by Biben and Misbah (2003) and Biben et al. (2005). The basis of this model can be conceived as intermediate between phase-field and explicit methods. They proposed a bending free energy of the form  F= 2



dV (Cˆ − C0 )

2 |∇ |

2

,

916 917 918 919 920 921 922

(55)

where Cˆ is the total curvature of each isosurface of constant . Thus, Cˆ is defined in the entire volume rather than in the membrane surface, albeit the presence of the delta-like function | ∇ | reduces the volume integral to a surface integral in the limit of thin interface, ˆ − xm ) → C at the membrane position so that in this limit C(x)ı(x xm . The mean curvature is numerically computed from the norˆ where nˆ = ∇ /|∇ |. Hence, mal vector to the isosurface, Cˆ = −∇ · n,

923

The theory of Ginzburg–Landau provides a basis for the energetic characterization of membranes. However, the time dependence of the interface is obviously a critical ingredient in the modeling of the membrane phenomenology. In non-equilibrium dynamics it is usual to assume that each small volume element is locally in thermodynamic equilibrium so that the whole system evolves toward a global equilibrium. This mesoscopic hypothesis is generically known as local equilibrium and constitutes a fundamental assumption for the subsequent theory. In the case of membranes, the separation in time scales between the rearrangement and diffusion of lipids, ∼10−9 s, and typical cell deformation times and mechanic response, ∼10−3 s, ensures that the local equilibrium approach applies correctly. In the framework of the Cahn–Hilliard theory (Cahn and Hilliard, 1959), the dynamic evolution of the order parameters is dictated by a diffusive equation,

∂ =∇· ∂t



885

915

(52)

where i refers to the lipid species present in the domain ˝i and then the membrane rigidity i and the spontaneous curvature C0,i depend on the specific domain ˝i . The whole membrane surface obviously arises from the composition of all these domains, ˝ = ˝1 ∪ ˝2 ∪ · · · The last term in (52) represents a line tension term accounting for the border between the different lipid domains. Wang and Du (2008) considered the problem of two coexistent species and introduced a phase field formulation based on two coupled order parameters, and . The main field represented the vesicle surface, whilst the auxiliary field  formed a perpendicular surface to the vesicle, so that the regions where both surfaces superpose define the domain of one of the species. Furthermore, the parameter coefficients of the free energy of depend on the field, () and c0 (), with two constant values at  = +1 and  = −1 and a smooth transition in the interface. Thus, they propose a free energy of the form

868

although in this scheme the interface is not tracked either, the mean curvature must be specifically computed, implying a number of stability numerical problems that require of a very fine tuning of the model. For instance, ∇ → 0 in the bulk, given rise to numerical divergences of the normal vector far from the interface. This model of specific curvature calculation reproduces the desired elastic behavior of the membrane, but it deviates from the original spirit of phase-field models.



M∇

ıFmem ı



.

(56)

It is not difficult to show that the total amount of order parameter is constant in the system. If the interface is locally in equilibrium, this implies that the total amount of each equilibrium phase is also constant. This property is specially convenient in the context of cell modeling, given that the conserved evolution ensures that the amount of fluid both inside and outside the cell frontiers remains constant, and the enclosed volume restraint of the Helfrich free energy (1) is directly fulfilled. Alternative models which consider a non-conserved dynamics, such as an Allen–Cahn dynamic equation (Du et al., 2004, 2005), explicitly introduce the volume conservation by adding the correspondent Lagrange multiplier. The Cahn–Hilliard equation (56) dictates the dynamics of the interface but, in many systems, hydrodynamic effects of the aqueous environment are also crucial in the membrane evolution. A paradigmatic example is the study of red blood cells and lipid vesicles while flowing along capillaries forced by an external flow. To model the interaction of the membrane with the surrounding fluid, it is usual to incorporate the Navier–Stokes equation to describe the dynamics of the fluid, and both equations are coupled describing the interaction between the membrane and the fluid. The complete model is

∂ + v · ∇ = M ∇ 2 mem . ∂t





924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939

940

941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961

(57)

962

(58)

963



∂v + (v · ∇ v) = −∇ P + fmem + ∇ 2 v + fext . ∂t

where fext is an external forcing driving the flow, such as a pressure difference (P)/L between the edges of a channel of length L, or a gravity force g. The advection term v · ∇ in (57) describes how the fluid accelerates the membrane, and the response exerted by the membrane is given by the force density fmem . The elastic force density fmem is obtained from the divergence of the stress tensor (31), f = ∇ · mem = − ∇ mem .

(59)

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

964 965 966 967 968 969 970 971

G Model

ARTICLE IN PRESS

CPL 4318 1–15

G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

979

The complexity of these equations and the geometries that are usually studied in membrane problems avoid an analytic integration of the model, and the four numerical methods are required. The Navier–Stokes equation can be integrated by different methods which must incorporate the coupling with the order parameter. Among others, new hybrid formulations have been developed for immerse boundary methods (Du and Li, 2011; Shao et al., 2013) and Lattice–Boltzman methods (Kendon et al., 2001).

980

4. Applications

972 973 974 975 976 977 978

992

In this section, we overview two different membrane related problems that have been studied by means of phase-field methods. The first case, the membrane instabilities triggered by hydrophobic insertion of polymers, concerns the local shape deformation induced by conformational changes in the membrane microstructure. The second case studies the morphological response of red blood cells to an external force, the flow induced in a highly confined channel by a pressure gradient, a problem of relevance in blood microcirculation. In spite of the different nature of the two problems, the flexibility and robustness of phase-field methods offers a good approach to both, by just incorporating slight modifications to the model.

993

4.1. Polymer induced membrane instabilities

981 982 983 984 985 986 987 988 989 990 991

994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014

1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032

Cell membranes do not only assume structural functions but also contribute to a wide variety of fundamental processes of the cell, including transport of nutrients and energy (Cai et al., 2007; Matteis and Luini, 2008). Molecules are enclosed in a lipid vesicle which can be transported from one region of the cell to other, allowing exchange of matter between the Golgi apparatus, the endoplasmic reticulum and the outer environment of the cell (i.e. exocitosis). A different mechanism of transport is tubulation, in which the formation of long tethers can connect two distant organelles, serving as conduits of direct communication. The spontaneous deformation of the membrane necessary to develop such specific shapes, both lipid buds and tubes, is a fascinating phenomena. Among other mechanisms, cells are able to induce an instability that triggers the formation of these structures by controlling the concentration of certain polymers that anchor to the membrane, generating an imbalance in the inner structure of the membrane that forces it to deform into new equilibrium shapes, relaxing the membrane tensions introduced by the anchored polymers. In this context, phase field models have proven to be powerful and elegant approaches for the understanding and characterization of these dynamic instabilities. 4.1.1. Pearling The pearling instability observed when monocomponent lipid vesicles are exposed to the presence of certain amphiphilic molecules represents an interesting model system to understand the formation of lipid vesicles in cells. Tsafrir et al. (2001) realized a series of experiments in which they introduced amphiphilic polymers in the surrounding region of a tubular vesicle, and after the molecule diffusion they observed the formation of arrays of pearls of different size. The molecular basis of this phenomena is found in the anchoring of the amphiphilic molecule that introduces its hydrophobic head into the lipid bilayer of the vesicle. The disruption induced by the polymer introduces an intrinsic asymmetry to the bilayer, and afterwards the equilibrium shape of the vesicle is no longer stable, giving rise to the development of pearls. Campelo and Hernández-Machado (2007a) made use of a phase field methodology to study the dynamic pearl formation. The model of Campelo and Hernández-Machado (2007a) assumes that spontaneous curvature depends on the polymer

11

concentration , given that a higher concentration implies a higher density of anchored polymers along the vesicle surface which force larger bending deformations. They therefore considered the free energy (51). The dependence between spontaneous curvature and polymer concentration is assumed to be linear, c0 (x, t) = (0) (1) c0 + c0 (x, t). The model assumes that the polymer diffusion is very fast, show that the membrane interacts with an homogeneous polymer concentration. Additionally, based on the good qualitative and quantitative agreement between numerical results and experiments, it can be inferred that the assumption of an spatially homogeneous concentration of polymers represents actually a good description. The free energy (51) is introduced in a conserved dynamic formalism (56), allowing an study of the temporal evolution of the deformation. The results presented in Campelo and Hernández-Machado (2007a) show that initially a first pearl is formed at the cap of the tube, connected to the main part by a thin neck, although the pearl does not pinch-off. Further addition of polymer triggers the gradual formation of pearls, all with a radius determined by the relevant spontaneous curvature, Rpearl = 1/c0 , so that the bending energy of the mebrane vanishes. However, beyond a critical value of spontaneous curvature, c0 = 2A/3V, where A and V are the surface area and volume of the tube, respectively, the membrane cannot support such a high asymmetry and the regular configuration, in which all the pearls have the same size, desestabilizes into a chain of increasingly smaller size (being the largest pearl the one situated at the tube cap), as shown in 6A. The size gradient increases for higher polymer concentrations, as it imply an important relaxation of the bending energy with respect to more homogeneous configurations. The results highlight the subtle control of the vesicle shape that can be achieved by the tuning of polymer concentration (Fig. 4). Q5

4.1.2. Tubulation The formation of long tubes in highly oblate vesicles was experimentally studied by Tsafrir et al. (2003). They found that the hydrophobic insertion of polymers in the lipid membrane triggers a sequence of complicate shape deformations. Initially, the addition of the polymer solution in the surrounding of a fluctuating vesicle 5a induces the formation of a bud and the suppression of the surface fluctuations. Then, the buds grow deforming into cylindrical tubes. After the polymer source is depleted, tubes stabilize in a bud-like shape, and eventually the vesicle absorbed the buds recovering a compact shape. The spontaneous curvature phase field model (51) can be applied to face this problem (Campelo and Hernández-Machado, 2008). In this case, the vesicle is subjected to the presence of an inhomogeneous gradient of polymer concentration, (x). An important assumption of the model is that the vesicle acts as a lipid reservoir, so that the tube has no lipid mass limitation during its growth. The phase-field results show that the tube is formed and grows for large concentration gradients, when the polymers are localized at the tube cap, but the homogenization of the polymer distribution leads to the membrane deformation into buds, structures with a circular body and thin neck connecting with the main vesicle body. The results show a qualitative nice agreement with experimental results, see Fig. 5B, capturing the growth of tubes when sharp polymer concentrations are present and then stabilizing in the mesostable configuration of small buds for more homogeneous concentrations after the polymer diffusion. The phase-field modeling demonstrates that the complicate processes of tubulation and budding can be entirely understood as the shape response of the vesicle to an inhomogeneous distribution of spontaneous curvature given by the anchoring of polymers.

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064

1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095

G Model CPL 4318 1–15 12

ARTICLE IN PRESS G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

Fig. 4. Red blood cell modeling with phase-field methods. In this case, the order parameter is defined in a two-dimensional domain an axisymmetry is applied in the horizontal axis. The initial configuration was an ellipse with the relevant reduced volume of these cells, vred = 0.6, and we impose rotational symmetry around the minor axis. The inner region of the ellipse was initialized at = +1 and the outer one corresponds to = −1. The interfaces relaxes in a fast time scale, acquiring its equilibrium, tanh-like smooth profile; then, the surface tension of the interface vanishes and the minimization of the shape is driven by curvature. In equilibrium, the shape obtained is the characteristic biconcave profile, the original result of Canham and Helfrich.

1096

1097 1098 1099 1100 1101 1102 1103

4.2. Red blood cells dynamics in microchannels Erythrocytes or red blood cells (RBCs, hereafter) are responsible for the transport of oxygen along the circulatory system, and their functionality critically depends on their deformation capability as they cross through very thin capillaries (when oxygen is delivered) of smaller section than the RBC size. RBCs are known to exhibit a fascinating behavior when flowing in confined channels, as a result of their complex elastic response to the external flow perturbation

(Abkarian et al., 2008). The RBC anucleated nature implies that its elastic properties are controlled by the plasma membrane. Several diseases, such as malaria and sickle cell anemia (Gallagher, 2005) modify the membrane elasticity of the cell affecting to the cell functionality and circulatory properties of the individual. The behavior of RBCs in confined channels is therefore an important subject for the understanding of blood circulation and it has applications in blood handling and pathology diagnosis (Toner and Irimia, 2005). For these reasons, the study of RBC flow has raised interest in recent years and different theoretical and numerical studies have focused on the subject (Noguchi and Gompper, 2005; Vlahovska et al., 2009; Kaoui et al., 2009). It has been proven that the dynamic behavior of RBCs during blood flow is not dominated by the presence of the cytoskeleton (Forsyth et al., 2011), and then the Helfrich theory represents an accurate description of the elastic behavior of the RBC membrane at these conditions. Therefore, phase-field methods provide an interesting pathway for the study of blood flow (Lázaro et al., in press). The phase-field model (44) reproduces the elastic behavior of a cell membrane. The equilibrium shape of a RBC (the so-called discocyte) can be obtained by minimizing the free-energy (44) for an ellipsoid with the appropriate surface area and volume, as both quantities are fixed in RBCs due to the membrane incompressibility and constant enclosed volume, respectively. Although the geometric characteristics of human RBCs present high interindividual variation, typical values are A = 140 m2 and V = 90 m3 , leading to a reduced volume red = V/(4R3 /3) ∼ 0.58, where R = (A/4)1/2 . The ellipsoid minimization can be performed by integration of equation (56), providing the equilibrium shape of the RBC, the socalled biconcave discocyte. The RBC must then be introduced in a flow driven by a pressure gradient, included in our model in the Navier–Stokes equation (58). The key parameters of the model are the bending rigidity of the cell membrane, , and the external forcing (P)/L that induces a flow velocity v. RBCs present a typical value of  = 2.0 × 10−19 J, and the relevant flow velocities observed in microchannels are of v∼0.1 cm/s. The competition between the viscous effects of the flow and the elastic ones of the cell membrane are given by the capillary number

C =  / =

a3 v , b

(60)

where   = a3 / and  = b/v are the relaxational time of the cell and the inverse of the shear rate, respectively. a is the RBC diameter

Fig. 5. Curvature driven instabilities: comparison between experiments and numerical results from the phase-field model. (A) Pearling instability. The upper plot shows the time evolution of the pearling process: the first pearl is formed at the cap and the instability extends subsequentially along the tube. In the bottom, an example of an inhomogeneous distribution of pearls, a favorable configuration for high polymer concentrations in comparison with the homogeneous pattern. (B) Tubulation. Comparison between the tube growth characteristic of the beginning of the process, when the polymer distribution is highly inhomogeneous and concentrated in the tube budding, when polymer concentration is more homogeneous and the structures deform into buds. The polymer concentration was fixed to a Gaussian profile cap, and the √ (x) = 0 /( 2) exp(−|x − xappl |2 /2 2 ), which increases from the vesicle surface to the point where the polymer solution is applied, xappl . The increase of the standard deviation  implies a more homogeneous polymer distribution, which can be related with the diffusion of the polymers in the bulk. Reprinted figures with permission from (A) F. Campelo and A. Hernández-Machado, Phys. Rev. Lett. 99 (2007) 088101, and (B) F. Campelo and A. Hernández-Machado, Phys. Rev. Lett., 100 (2008) 158103. Copyright 2014 by the American Physical Society.

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141

1142

1143 1144

G Model CPL 4318 1–15

ARTICLE IN PRESS G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

Fig. 6. (A) Comparison between phase-field results and experiments of RBCs flowing in highly confined channels: slippers (left) and parachutes (right). Experimental figure from Ref. (Tomaioulo et al., 2009). (B) RBC morphologies for a cell initially placed at the channel axis, as a function of the capillary number C = a3 v/b. The crosses represent the cell center of mass. RBCs initially retain a centered position, with a shape similar to the equilibrium discocyte. For intermediate flow velocities, RBCs are localized at an off center position, where they assume a slipper shape. For high flow velocities, RBCs adopt a symmetric parachute shape centered in the channel. (C) Bending energy of the RBC morphological sequence. The plot shows the three main morphological regimes. Centered discocytes are characterized by a gradually increasing bending energy, until the energy peaks before RBCs adopt the slipper shape. This deformation allows a relaxation of the bending energy, as shown in the inset which represents a zoom up of the curve. Parachutes are strongly penalized shapes. (D) Ratio between the bending and incompressibility terms of the elastic membrane energy (48). At low flow velocities the bending is dominant and completely determines the morphological response. On the contrary, at high flow velocities, when the cell is subjected to severe stretchings, the incompressibility contribution is dominant limiting the RBC deformation.

1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172

and b represents the channel width. For a fixed bending rigidity, the capillary number is proportional to the flow velocity. Depending on the relative RBC rigidity with respect to the incoming flow, RBCs flowing in microchannels assume different shapes. The complete morphological sequence is shown in Fig. 6B. At low flow velocities, RBCs flow unperturbed maintaining a discocytic shape, and showing a slight deformation that enables them to accommodate their shape to the profile of the external flow. For RBCs located at the axis this implies that they increasingly bend with the flow velocity. At some point, RBCs flowing close to the axis are forced off center and migrate laterally to an asymmetric position in the channel, where they acquire the so-called slipper morphology. RBCs assume an asymmetric shape that allows them to orient its axis with the local flow profile. As the flow velocity increases, however, RBCs gradually acquire an horizontal inclination, with their axis parallel to the channel walls. Further increasing the flow velocity leads to instability of the slipper shape, and the cell recovers a symmetric shape returning to the channel axis and developing a parachute-like morphology. The presence of slippers and parachutes is restricted to very thin channels, as they are critically influenced by the hydrodynamic interactions with the walls. In thicker channels, where velocity gradients are lower, at the RBC scale the cell interacts with a more planar flow profile. Accordingly, they are not forced to bend, maintaining a discocytic shape even at high flow velocities. The sequence of RBC morphologies can be interpreted from the bending energy of each shape, as shown in Fig. 6C. Initially, RBCs are subjected to slight but gradual deformations, and their bending

13

energy increases linearly with the flow velocity. At some point, the lateral position of the RBC allows the cell to deform into a slipper, a more symmetric shape than the previous curved discocytes, and the cell deformation energy decreases. The RBC off center position is characterized by a lower bending energy along a considerable range of flow velocities, until eventually the cell shape deforms into a parachute. The parachute morphology implies a high deformation from the equilibrium discocyte, and consequently it is strongly energetically penalized. The relative contributions to the elastic energy of the bending and incompressibility terms of (48) is shown in Fig. 6D. At low flow velocities, the weak flow forces induce negligible stretching to the membrane and hence the incompressibility effect vanishes. As the flow increases and cell deformations are larger, the incompressibility becomes dominant. The critical value separating both regimes corresponds to the off-center slipper. The results reproduce the observed shapes of RBCs (see Fig. 6A) and highlight the relevance of the cell elasticity and its interplay with the flow, and especially the effect of the channel walls. Accordingly, rigid RBCs, such as those infected by malaria, could show a different morphological behavior thus affecting to the flow properties of blood, such as a higher viscosity. This problem constitutes a good example of how small changes in the structure of the cell membrane can impair the functionality of a much larger, macroscopic system, as the results of the phase-field model prove.

5. Conclusions The complexity of cell membranes composition forces to refuse a detailed description of the membrane microstructure and adopt a mesoscopic approach, in which the internal physics of the membrane is captured by the effective variables at the mesoscale. The mechanical response of the cell membrane is typically described by the theory of elasticity. The small thickness of the membrane relative to the cell length scale allows the treatment of the membrane as a two-dimensional elastic layer. In particular, Helfrich proposed that membrane elasticity is dominated by bending rather than the in-plane contribution of shear. He also incorporated to the model a spontaneous curvature, which accounts for asymmetries in the internal structure of the membrane, reflecting that membranes sometimes show a preferential curvature. Finally, the strong bonds formed between lipids lead to an incompressibility condition in the membrane surface area. Based on these main assumptions, the Helfrich theory has proven to be extraordinarily successful, with an outstanding number of cell membrane problems studied by means of this model. Based on the Helfrich theoretical framework, we have presented the phase-field formalism as a powerful tool for approaching complex phenomena related with dynamics and morphology of biological membranes. Phase field models have been applied to the study of different interface problems, but only recently for membrane modeling. The phase-field theory makes use of an order parameter which has two stable phases, and the interface connects both phase domains by a smooth profile. Thereby, one needs to solve the dynamics of the field, avoiding the complex treatment of the moving boundary condition of the interface. Phase-field methods allow a fine tuning of the interface properties, as demonstrated in this article, a basic property when dealing with membranes given the specific elastic characteristics of these structures. We have shown the link of the phase-field model with the elastic theory of membranes, providing a flexible methodology for the study of membranes with different properties. The membrane model presented here effectively describes the internal stresses of the membrane, which in turn determine the elastic response of the cell. The phase-field can be coupled to a velocity field describing the hydrodynamics of the surrounding fluid, since in

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196

1197

1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235

G Model CPL 4318 1–15 14

ARTICLE IN PRESS G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

many situations the hydrodynamic effect plays a critical role in the membrane dynamics. The robust physical basis of the complete 1238 phase-field Navier–Stokes model allows to deal with very differ1239 ent membrane phenomena, as shown in the Applications section, 1240 where we have presented two examples of membrane modeliza1241 tion. First, we have explained vesicle pearling and tubulation. These 1242 phenomena are triggered by the insertion of amphiphile polymers, 1243 which anchors to the membrane and induce a preferential curva1244 ture. In the phase-field approach, the interface profile is modified 1245 introducing an asymmetry, forcing the vesicle to deform adopt1246 ing the pearl (or tube) configurations. The second case presented 1247 concerns the mechanics and morphological response of RBCs when 1248 flowing through highly confined channels. The interplay between 1249 membrane elasticity and external flow leads to a complex sequence 1250 of deformed cells. The phase-field model allows to understand the 1251 Q6 relevance of membrane properties (e.g. rigidity) in the fluid behav1252 ior of blood. 1236 1237

1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 Q7 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313

References Abkarian, M., Faivre, M., Horton, R., Smistrup, K., Best-Popescu, C.A., Stone, H.A., 2008. Cellular-scale hydrodynamics. Biomed. Mater. 3 (3), 034011. Alberts, B., Bray, D., Raff, M., Roberts, K., Watson, J.D., 1994. Molecular Biology of the Cell. Garland Science, New York. Baumgart, T., Hess, S., Webb, H., 2003. Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension. Nature 425, 821824. Biben, T., Kassner, K., Misbah, C., 2005. Phase-field approach to three-dimensional vesicle dynamics. Phys. Rev. E 72, 041921. Biben, T., Misbah, C., 2003. Tumbling of vesicles under shear flow within an advectedfield approach. Phys. Rev. E 67, 031908. Brannick, J., Liu, C., Qian, T., Sun, H., 2014. Diffuse interface methods for multiple phase materials: an energetic variational approach., pp. 1, arXiv:1402.5375. Cahn, J.E., Hilliard, J.E., 1959. Free energy of a nonuniform system: II. Nucleation in a two-component incompressible fluid. J. Chem. Phys. 31, 688–699. Cai, H., Reinisch, K., Ferro-Novick, S., 2007. Coats, tethers, Rabs, and SNARES work together to mediate the intracellular destination of a transport vesicle. Dev. Cell. 12, 671–682. Campelo, F., Hernández-Machado, A., 2006. Dynamic model and stationary shapes of fluid vesicles. Eur. Phys. J. E 20 (1), 37–45. Campelo, F., Hernández-Machado, A., 2007a. Model for curvature-driven pearling instability in membranes. Phys. Rev. Lett. 99, 088101. Campelo, F., Hernández-Machado, A., 2007b. Shapes instabilities in vesicles: a phasefield model. Eur. Phys. J. Special Topics 143, 101–108. Campelo, F., Hernández-Machado, A., 2008. Polymer-induced tubulation in lipid vesicles. Phys. Rev. Lett. 100, 158103. Canham, P.B., 1970. The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J. Theor. Biol 26, 6181. Deserno, M., 2009. Mesoscopic membrane physics: concepts, simulations, and selected applications. Macromol. Rapid. Commun. 30, 752–771. Deserno, M., 2014. Chemistry and physics of lipids: special issue: membrane mechanics: from the molecular to the cellular scale. In: Jesús, P.-G., Francisco, M. (Eds.), Ch. Fluid Lipid Membranes: From Differential Geometry to Curvature Stresses. Elsevier, pp. 716–755. Deuling, H.J., Helfrich, W., 1977. Theoretical explanation for myelin shapes of red blood cells. Blood Cells 3 (3), 713720. do Carmo, M.P., 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, NJ, USA. Du, Q., Li, M., 2011. On the stochastic immersed boundary method with an implicit interface formulation. Discrete Contin. 15 (2), 373–389. Du, Q., Liu, C., RYham, R., Wang, X., 2005. A phase-field formulation of the Willmore problem. Nonlinearity 18, 1249–1267. Du, Q., Liu, C., Wang, X., 2004. A phase field approach in the numerical study of the elastic bending energy for vesicle membranes. J. Comp. Phys. 198, 450–468. Evans, E., Needham, D., 1987. Physical properties of surfactant bilayer membranes: thermal transitions, elasticity, rigidity, cohesion, and colloidal interactions. J. Phys. Chem. 91, 42194228. Evans, E.A., 1974. Bending resistance and chemically induced moments in membrane bilayers. Biophys. J. 14 (12), 923. Forsyth, A.M., Wan, J., Ristenpart, W.D., Stone, H.A., 2011. The dynamic behavior of chemically stiffened red blood cells in microchannel flow. Microvas. Res. 80, 37–43. Gallagher, P.G., 2005. Red blood cell membrane disorders. Hematol. Am. Soc. Hematol. Educ. Program 13, 8. Gatzke, T.D., Grimm, C.M., 2006. Estimating curvature on triangular meshes. Int. J. Shape Model. 12, 129. Gompper, G., Schick, M., 1994. Self Assembling Amphiphilic Systems. Academic Press Limited, London. Gompper, G., Zschocke, S., 1991. Elastic properties of interfaces in a Ginzburg–Landau theory of swollen micelles, droplet crystals and lamellar phases. Europhys. Lett. 16, 731.

Gompper, G., Zsckocke, S., 1992. Ginzburg–Landau theory of oil–water–surfactant mixtures. Phys. Rev. A 46 (8), 4836. Helfrich, W., 1973. Elastic properties of lipid membranes: theory and possible experiments. Z. Naturforsch. C 28, 693703. Helfrich, W., 1981. Liquids at Interfaces, Les Houches Session XXXV. In: Balian, R., Kleman, M., Poirier, J. (Eds.), Ch. Amphiphilic mesophases made of defects. North Holland Publishing Company, Amsterdam, pp. 716–755. Hu, M., de Jong, D.H., Marrink, S.J., Deserno, M., 2013. Gaussian curvature elasticity determined from global shape transformations and local stress distributions: a comparative study using the martini model. Faraday Discus. 161, 365–382. Iglic, A., 1997. A possible mechanism determining the stability of speculated red blood cells. J. Biomech. 35-40, 923. Kaoui, B., Biros, G., Misbah, C., 2009. Why do red blood cells have asymmetric shapes even in a symmetric flow? Phys. Rev. Lett. 103, 188101. Kaoui, B., Krueger, T., Harting, J., 2012. How does confinement affect the dynamics of viscous vesicles and red blood cells? Soft Matter 8, 9246–9252. Kendon, V.M., Cates, M.E., Pagonabrraga, I., Desplat, J.-C., Blandon, P., 2001. Inertial effects in three-dimensional spinodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study. J. Fluid Mech. 440, 147–203. Kuznetsova, T.G., Starodubtseva, M.N., Yegorenkov, N.I., Chizhik, S.A., Zhdanov, R.I., 2007. Atomic force microscopy probing of cell elasticity. Micron 38, 824–833. Landau, L.D., Lifshitz, E.M., 1999. Theory of Elasticity. Butterworth-Heinmann, Oxford. Lázaro, G.R., Hernández-Machado, A., Pagonabarraga, I., 2014. Rheology of Red Blood Q8 Cells Under Flow in Highly Confined Microchannels (in press). Lázaro, G.R., Melzak, K.M., Toca-Herrera, J.L., Pagonabarraga, I., Hernández-Machado, A., 2013. Elastic energies and morphologies of the frst stages of the discoechinocyte transition. Soft Matter 9, 6430–6441. Li, J., Dao, M., Lim, C.T., Suresh, S., 2005. Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Biophys. J. 88, 3707–3719. Lim, G.H.W., Wortis, M., Mukhopadhyay, R., 2002. Stomatocyte–discocyte–echinocyte sequence of the human red blood cell: evidence for the bilayer couple hypothesis from membrane dynamics. Proc. Natl. Acad. Sci. U. S. A. 99 (26), 16766–16769. Malevanets, A., Kapral, R., 1999. Mesoscopic model for solvent dynamics. J. Chem. Phys. 110, 8605–8613. Marrink, S.J., Risselada, H.J., Yefimov, S., Tieleman, D.P., de Vries, A.H., 2007. The MARTINI force field: grained model for biomolecular simulations. J. Phys. Chem. B 111, 7812–7824. Marrink, S.J., Tieleman, D.P., 2013. Perspective on the Martini model. Chem. Soc. Rev. 42, 6801–6822. Matteis, M.A.D., Luini, A., 2008. Exiting the Golgi complex. Nat. Rev. Mol. Cell. Biol. 9, 273–284. McWhirter, J.L., Noguchi, H., Gompper, G., 2008. Flow-induced clustering and alignment of vesicles and red blood cells in microcapillaries. Proc. Natl. Acad. Sci. U. S. A. 106 (15), 6039–6043. Noguchi, H., Gompper, G., 2005. Shape transitions of fluid vesicles and red blood cells in capillary flows. Proc. Natl. Acad. Sci. U. S. A. 102 (40), 14159–14164. Peng, Z., Li, X., Pivkin, I.V., Dao, M., Karniakidis, G.E., Suresh, S., 2013. Lipid bilayers and cytoskeletal interactions in a red blood cell. Proc. Natl. Acad. Sci. U. S. A. 110, 13356–13361. Peskin, C.S., 2002. The immersed boundary method. Acta Numer. 11, 1–39. Petrov, A.G., Bivas, J., 1984. Elastic and flexoelectric aspects of out-of-plane fluctuations in biological and model membranes. Prog. Surf. Sci. 18, 389. Pozrikidis, C., 1992. Boundary Integral & Singularity Methods for Linearised Viscous Flow. Cambridge University Press, New York. Pozrikidis, C., 1995. Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow. Soft Matter 297, 123–152. Rawicz, W., Olbrich, K.C., McIntosh, T., Needham, D., Evans, E., 2000. Effect of chain length and unsaturation on elasticity of lipid bilayers. Biophys. J. 79 (1), 328–339. Sackmann, E., 1995. Handbook of biological physics. In: Lipowsky, R., Sackmann, E. (Eds.), Ch. Biological Membranes Architecture and Function, vol. 1. Elsevier. Safran, S.A., 1994. Statistical Thermodynamics of Surfaces, Interfaces, and Membranes. Addison-Wesley, Reading. Seifert, U., Berndel, K., Lipowsky, R., 1991. Shape transformations of vesicles: phase diagram for spontaneous-curvature and bilayer-couple models. Phys. Rev. A 44 (2), 11821202. Seifert, U., Lipowsky, R., 1995. Structure and Dynamics of Membranes, Handbook of Biological Sciences, vol. 1, Ch. The Morphology of Vesicles. Elsevier, Amsterdam. Shao, J.Y., Shu, C., Chew, Y.T., 2013. Development of an immersed boundary phase field-lattice Boltzmann method for Neumann boundary condition to study contact line dynamics. J. Comp. Phys. 234, 8–32. Sheetz, M.P., Singer, S.J., 1974. Biological membranes as bilayer couples. A molecular mechanism of drug–erythrocyte interactions. Proc. Natl. Acad. Sci. U. S. A. 71, 4457–4461. Shillcock, J.C., Lipowsky, R., 2006. The computational route from bilayer membranes to vesicle fusion. J. Phys.: Condens. Matter 18, 1191–1219. Simons, K., Vaz, W., 2004. Model systems, lipid rafts, and cell membranes. Annu. Rev. Biophys. Biomol. Struct. 33, 269–295. Steinbach, I., 2009. Phase-field models in materials science. Mod. Simul. Mater. Sci. Eng. 17, 073001. Tomaioulo, G., Simone, M., Martinelli, V., Rotoli, B., Guido, S., 2009. Red blood cell deformation in microconfined flow. Soft Matter 5, 3736–3740. Toner, M., Irimia, D., 2005. Blood-on-a-chip. Annu. Rev. Biomed. Eng. 7, 77–103.

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397

G Model CPL 4318 1–15

ARTICLE IN PRESS G.R. Lázaro et al. / Chemistry and Physics of Lipids xxx (2014) xxx–xxx

1398 1399 1400 1401 1402

Tsafrir, I., Caspi, Y., Guedeau-Boudeville, M.A., Arzi, T., Stavans, J., 2003. Budding and tubulation in highly oblate vesicles by anchored amphiphilic molecules. Phys. Rev. Lett. 91, 138102. Tsafrir, I., Sagi, D., Arzi, T., Guedeau-Boudeville, M.A., Frette, V., Kandel, D., Stavans, J., 2001. Pearling instabilities of membrane tubes with anchored polymers. Phys. Rev. Lett. 86, 1138–1141.

15

Vlahovska, P.M., Podgorski, T., Misbah, C., 2009. Vesicles and red blood cells in flow: from individual dynamics to rheology. C. R. Phys. 10, 775–789. Wang, X., Du, Q., 2008. Modeling and simulations of multi-component lipid membranes and open membranes via diffuse interface approaches. J. Math. Biol. 56 (3), 347–371.

Please cite this article in press as: Lázaro, G.R., et al., Phase-field theories for mathematical modeling of biological membranes. Chem. Phys. Lipids (2014), http://dx.doi.org/10.1016/j.chemphyslip.2014.08.001

1403 1404 1405 1406 1407