Modeling the phase separation in binary lipid membrane under externally imposed oscillatory shear flow

Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

Contents lists available at ScienceDirect

Colloids and Surfaces B: Biointerfaces journal homepage: www.elsevier.com/locate/colsurfb

Modeling the phase separation in binary lipid membrane under externally imposed oscillatory shear flow Xiao-Bo Chen, Li-Sha Niu ∗ , Hui-Ji Shi Key Laboratory of Failure Mechanics, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, PR China

a r t i c l e

i n f o

Article history: Received 15 February 2008 Received in revised form 9 April 2008 Accepted 9 April 2008 Available online 16 April 2008 Keywords: Modeling Phase separation Binary lipid membrane Oscillatory shear flow

a b s t r a c t By adding external velocity terms, the two-dimensional time-dependent Ginzburg–Landau (TDGL) equations are modified. Based on this, the phase separation in binary lipid membrane under externally imposed oscillatory shear flow is numerically modeled employing the Cell Dynamical System (CDS) approach. Considering shear flows with different frequencies and amplitudes, several aspects of such a phase evolving process are studied. Firstly, visualized results are shown via snapshot figures of the membrane shape. And then, the simulated scattering patterns at typical moments are presented. Furthermore, in order to more quantitatively discuss this phase-separation process, the time growth laws of the characteristic domain sizes in both directions parallel and perpendicular to the flow are investigated for each case. Finally, the peculiar rheological properties of such binary lipid membrane system have been discussed, mainly the normal stress difference and the viscoelastic complex shear moduli. © 2008 Elsevier B.V. All rights reserved.

1. Introduction The hypothesis that cellular membranes are organized in discrete dynamic entities called lipid rafts has been formulated for nearly 20 years [1]. Lipid rafts enriched in cholesterol and sphingolipids could serve many unique biological membrane-associated functions, including dynamical sorting and signaling, transporting vesicle biogenesis, cell surface proteolysis, cell polarization, and viral budding, etc. [2–5]. Actually, it has been demonstrated that for homogenous bilayers formed by purified phospholipids or sphingolipids, a sharp transition between gel and liquid phases will occur at a characteristic melting temperature [6]. In biological membranes containing sphingolipids and phospholipids, the different packing ability of the two lipids probably leads to phase separation [7]. However, the limitation of currently existing research techniques makes it unworkable to perform an investigation into this subject with the veritable biological membrane. On one hand, the composition of the membrane of living eukaryotic cells is extremely complex, consisting of up to 500 different lipid species with numerous proteins and other large molecules. On the other hand, the membrane is structurally and dynamically coupled to the extracellular matrix and the cytoskeleton network [8]. Indeed, for many years such phase behavior has been studied with mixtures of membraneforming lipids, and these early studies are important in developing

∗ Corresponding author. Tel.: +86 10 62772921; fax: +86 10 62781824. E-mail address: [email protected] (L.-S. Niu). 0927-7765/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.colsurfb.2008.04.002

the fluid mosaic model of cellular membranes [9]. A lot of experimental studies of binary or ternary lipid mixtures have been carried out, and the direct observation of raft-like domains in these modeled biological membranes has provided a tangible proof for the coexistence of liquid-ordered and liquid-disordered phases [10–12]. These “artificial rafts” are a reasonable model of raft-containing cell membranes, though being crude undoubtedly [13]. In recent years, more and more theoretical forecasts and numerical simulations on this phase-separation phenomena have come into our sight. Kumar and Rao introduced an off-lattice kinetic Monte Carlo technique for studying the dynamics of multicomponent membranes in vesicle form [14]. Laradji and Kumar recently studied domain coarsening in multi-component bilayer fluid vesicles via dissipative particle dynamics in the absence of spontaneous curvature [15]. Wallace et al. presented the phase separation in a model asymmetric membrane by also using Monte Carlo method [16]. Wada investigated the dynamics of phase separation in two-component fluid membranes confined between parallel and rigid bounding plates utilizing the Cell Dynamical System (CDS) approach [17]. Reigada et al. addressed a simple kinetic model of a two-component deformable and reactive bilayer [18]. As a matter of fact, when lipids cluster into locally planar bilayer structures to form biological membranes, the surrounding must be an aqueous environment which will certainly bring in the effect of shear flow. Moreover, biological cells are mostly living in liquid environment and undergoing shear forces caused by various kinds of flows. Illustratively, it is proved that the membrane of the red blood cell (RBC) is easily deformed by shearing forces [19],

204

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

method [25], the mixing free energy caused by the interaction between the two components can be presented as

  

Gmix =

−˛ ln[cosh()] + ˙

Fig. 1. Schematic representation of a lipid membrane composed of two amphiphiles A and B with an externally imposed oscillatory shear flow (assume that the membrane is placed on the water surface, but the effects from such quiescent support are neglected).

2. Model: total energy functional and governing equations Schematic representation of our model is shown in Fig. 1, we describe such model as an idealized two-component membrane consisting of two incompatible amphiphiles A and B. Assuming that A and B are the local volume fractions of A and B components respectively, we easily get the following relationship A + B = 1 to guarantee the incompressibility. Further, the membrane is defined as a two-dimensional plane surface with a concentration difference order parameter , where  = A − B . So as to satisfy the structural stability of the CDS



dx dy

(1)

where the parameters ˛, ˇ, D depend on the properties of the components and also on the temperature, ˙ is the spatial field of integration. Defining h(x,y) as the off-plane displacement field for the local separation from the flat conformation, the deformable membrane surface can be described by [x,y,h(x,y)]. For surfaces that are nearly flat with only gradual variations of h, this representation is valid actually. The mean curvature H is as following 2

H=

2

1 (1 + h,x )h,yy − 2h,x h,y h,xy + (1 + h,y )h,xx 3/2 2 (1 + h2 + h2 ) ,x

which affects its physiological function of transporting oxygen [20] and determines the hydrodynamic properties of the whole blood [21]. And it has also been demonstrated that endothelial cells (ECs) are subjected to the shear stress resulting from blood flow and are able to convert mechanical stimuli into intracellular signals that affect cellular functions, e.g. proliferation, apoptosis, migration, permeability, and remodeling, as well as gene expression [22]. Actually, Bray et al. have studied an anisotropic version of the Burgers equation, constructed to describe thermal fluctuations of an interface in the presence of a uniform shear flow [23]. Smith et al. discussed the response of a multi-component budded vesicle to an imposed shear flow using dissipative particle dynamics [24]. In the present paper, representatively and for simplicity, we will investigate the phase-separation process of a binary lipid membrane subjected to an externally imposed oscillatory shear flow. This situation would be of great experimental relevance especially for probing the viscoelastic properties of the phase-separating binary mixture. In this work, the kinetics of this process is modeled fundamentally grounding on the two-dimensional timedependent Ginzburg–Landau (TDGL) equations, and an external velocity term on behalf of the flow convection is added respectively. Shear flows with different frequencies and amplitudes are taken into account individually. Such phase-separation process is primarily described via snapshots of the membrane morphology and the simulated scattering patterns. And further onwards, quantitative analysis is performed, concerning the characteristic domain sizes and the rheological properties of the membrane system. The outline of this paper is as follows. In Section 2, the model of our simulation will be specified, including the total energy functional and governing equations. Then in Section 3, the governing equations are numerically solved utilizing CDS technique, details on the method as well as the implementation are given for convenience of the readers. Subsequently in Section 4, simulated results are reported and analytically discussed. Finally in Section 5, some conclusions of this work are presented.

ˇ 2 D 2  + (∇ ) 2 2

(2)

,y

where h,x = ∂h/∂x, h,xy = ∂h/(∂x∂y), etc. According to the present model, we have h,x  1 and h,y  1, so the approximation H =



(1/2)(∇ 2 h/ 1 + |∇ h|2 ) is properly allowed. The rigidity of the membrane leads to an elastic bending energy contribution to the total energy [26,27], that is Gbending =

 2

 = 2

 

2

[2H − C0 ()] dA ˙

 

2

2

2

{[∇ 2 h − C0 ()] + 12 [C0 ()] (∇ h) } dx dy

(3)

˙

where  is the bending rigidity modulus, C0 () is the spontaneous curvature, which reflects the shape asymmetry between the two components, usually simplified to be C0 () = − (where  is called coupling constant). And the area measure dA =  2

1 + (∇ h) dx dy is expanded to second order in h, which has given 2 2 rise to the term (1/2)[C0 ()] (∇ h) . The two-dimensional flat membrane can be assumed to be a small part of a large vesicle exactly. For this reason, surface tension which suppresses membrane deviations from the flat state should also be considered. Here comes another energy contribution [16], i.e. the so-called frame energy

 

Gframe = ˙

 2 (∇ h) dx dy 2

(4)

where  is the surface tension. Thus, with all the energy contributions analyzed above, the total energy functional of the modeled binary lipid membrane reads F = Gmix + Gbending + Gframe

   F = ˙

2  2 (∇ h + ) + 2

− ˛ ln[cosh()] +



()  + 2 4

ˇ 2 D 2  + (∇ ) 2 2

2

 (∇ h)

2

 dx dy

(5)

For simplicity, as has been described in Fig. 1, we consider the case when the binary lipid membrane is subjected to an external oscillatory shear flow with a zero-divergence (i.e. ·V = 0) velocity field, being expressed as V=

dS(t) yex = (yω cos ωt, 0, 0) dt

(6)

where x denotes the flow direction and y denotes the shear gradient direction, ex is the unit vector of the x direction, dS(t)/dt represents the shear rate with strain S(t) =  sin ωt, ω is the angular frequency, and  is the strain amplitude.

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

Here dimensionless parameters are introduced to simplify all the notations mentioned above, just to avoid the complex work of analyzing dimensions: f =

F kB T

t = , t0

,

y˜ =

y , a0

˜ = D

D , kB T

˜ = 

a20 kB T

= ,

˝ = ωt0 ,

˜ = a0 , 

,

h h˜ = , a0

V˜ =

˛ ˜ =

 , ˜ = kB T

˛a20 kB T

˜ = ˇ

,

kB T

,

V t0 = ( y˜ ˝ cos ˝ , 0, 0). a0

f = ˜ ˙



˜ ˜ 2+ (∇ 2 h˜ + ) 2

−˛ ˜ ln[cosh()] +

˜ 2 ˜ ( ˜ )  + 2 4

˜ ˜ ˇ D 2 2 + (∇ ) 2 2



˜ 2 (∇ h)

 d˜x d˜y

(7)

with  and 2 representing the gradient operator and Laplacian operator in the dimensionless unit plane. To describe the time evolution of the order parameter  and ˜ via adding external velocity that of the off-plane displacement h, terms to the two-dimensional time-dependent Ginzburg–Landau equations, the convection–diffusion equations will be formed as following [28,29] ∂ +∇ · (V˜ ) = M ∇ 2 ∂ ∂h˜ + ∇ · (h˜ V˜ ) = −N ∂





ıf ı

ıf ıh˜



⇒

⇒

∂ = −V˜ · ∇  + M ∇ 2 ∂

∂h˜ = −V˜ · ∇ h˜ − N ∂



ıf ıh˜





ıf ı



(8)

(9)

where M and N are the phenomenological lateral mobility parameters. Substituting these variables and expressions into Eqs. (8) and (9), we obtain the governing equations for this phase˜ ∇ h) ˜ 2 in Eq. (10) and separation process. (The items (1/2M ˜ ( 2˜ 2 2 ˜ ˜ 1/2N ˜  [ ∇ h + (∇ ) · (∇ h)] in Eq. (11)) behaving small quantity have been ignored.) ∂ ∂ ˜ ˜ 2 + ˇ) = − ˝ cos ˝ · y˜ + M ∇ 2 [(˜  ∂ ∂˜x ˜ ˜ ∇ 2 h] ˜ ∇ 2  + ˜  −˛ ˜ tanh() − D

(10)

∂h˜ ∂h˜ ˜ ∇2 −  ˜ ˜ ∇ 2 h] = − ˝ cos ˝ · y˜ − N[˜ ∇ 2 ∇ 2 h˜ + ˜  ∂ ∂˜x

(11)

Noting that in the governing equations obtained in Section 2, the temporal field and the spatial field are coupling each other, ˜ Moreover, it is and so are the two focused parameters  and h. obvious that the governing equations are highly nonlinear. Hence, it appears really impossible to procure any analytical solution, and numerical-solving methods become the only alternative. Firstly, we discretize the temporal field into small segments with an increment step  , and the integration form of the equations in the temporal field would be



− ˝ cos ˝ · y˜



(13)

To handle the highly nonlinear differential expressions in the two-dimensional space, we would like to employ Cell Dynamical System approach proposed by Oono and Puri [25]. Indeed, CDS has been certified to be an efficient numerical integrating method to solve complicated nonlinear differential equations. In the twodimensional CDS model, the system is discretized into a L × L square lattice (with L = 128 in the current simulation) of the cell size a0 . Therefore, the order parameter field (n, ) and the off-plane dis˜ placement field h(n, ) are assigned to each lattice point n = (nx ,ny ) respectively, with nx ,ny = 1, 2, 3, . . ., L. Here in our model, periodic boundary condition has been imposed in the direction parallel to the flow [25], that is (nx + P(L − 1), ny , ) = (nx , ny , )

(14)

˜ x , ny , ) ˜ x + P(L − 1), ny , ) = h(n h(n

(15)

whereas in the direction perpendicular to the flow, it should be dissimilar [28] (nx + S( )Q (L − 1), ny + Q (L − 1), ) = (nx , ny , )

(16)

˜ x , ny , ) ˜ x + S( )Q (L − 1), ny + Q (L − 1), ) = h(n h(n

(17)

where P and Q are arbitrary integers, S( ) is the instantaneous strain at the moment . Being the quintessence of the CDS approach, the gradient operator and the Laplacian operator in such two-dimensional square lattice are discretely approximated by [25,30]

∇ (nx , ny ) =

∇2

=

1 ( 2

(nx + 1, ny ) −

−

(nx , ny − 1))

(18)

−

(19)

1 i

6

+

1

ii

12

i

(nx − 1, ny ),

(nx , ny + 1)

ii

˜ i represents nearest-neighbor where is standing for either  or h, lattice points and ii represents next-nearest-neighbors. With these fundamental principles, the discrete forms corresponding to Eqs. (10) and (11) finally become (nx , ny , +  ) = (nx , ny , ) +  ·



− ˝ cos ˝ ·

1 y˜ [(nx + 1, ny , ) 2



− (nx − 1, ny , )] + M ∇ 2 I (nx , ny , )

3. Implementation of numerical solution

( +  ) = ( ) +  ·

∂h˜ − N[˜ ∇ 2 ∇ 2 h˜ ∂˜x

− ˝ cos ˝ · y˜

˜ ∇2 −  ˜ ˜ ∇ 2 h] + ˜ 

where kB is the Boltzmann constant, T is the temperature of the system, a0 and t0 are the units of length and time scales, respectively. Then Eq. (5) turns to be

  

 ˜ +  ) = h( ) ˜ h( +  ·

x x˜ = , a0

ˇa20

205

∂ ˜ ˜ 2 + ˇ) + M ∇ 2 [(˜  ∂˜x



˜ ˜ ∇ 2 h] ˜ ∇ 2  + ˜  −˛ ˜ tanh() − D

(12)

(20)

˜ x , ny , +  ) h(n ˜ x , ny , ) +  · {− ˝ cos ˝ · 1 y˜ [h(n ˜ x + 1, ny , ) = h(n 2 ˜ x − 1, ny , )] + N ∇ 2 I ˜ (nx , ny , )} − h(n h

(21)

where I (n, ) and Ih˜ (n, ) are the discrete thermodynamic forces given by ˜ ˜ 2 + ˇ)(n, I (n, ) = (˜  ) − ˛ ˜ tanh[(n, )] ˜ ˜ ∇ 2 h(n, ˜ ∇ 2 (n, ) + ˜  −D ) ˜ ˜ ˜ h(n, ˜ Ih˜ (n, ) = ˜ ∇ 2 h(n, ) + ˜ (n, ) −  )

(22) (23)

˜ = 1, D ˜ = Standard parameters have been fixed to be ˛ ˜ = 1.3, ˇ ˜ = 0.2,  ˜ = 0.01, and M = N = 0.05 in the present 0.5, ˜ = 10/3, 

206

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

simulations [17,25], and the dimensionless time increment is set to be  = 0.5. For the current set of parameters, the interface width  ˜ ˛ D/ ˜ − 1 = 1.29. Such val between two domains should be = ues of the parameters are just coinciding with the experimental values for real biological membranes. Exactly speaking, the typical value of the bending rigidity modulus  is about the order of 10kB T [31], the surface tension may approach  ∼ kB T/(nm)2 , and the coupling constant is  ∼ (100 nm)−1 [32]. When experimentally imaging the domains in binary lipid membranes, the suitable scanning scale approximates 10 ␮m [33], and the diffusion coefficient of lipids is reported as M ∼ 104 (nm)2 /s [34]. With such quantities, the previously mentioned space unit could be determined to be a0 ∼ 100 nm, and the typical time unit can be given as t0 ∼ 0.01 s. When we change the shear frequency as ˝ = 0.0002, 0.0005, 0.002, 0.005, 0.008, 0.01, the shear amplitude is fixed at = 0.2. In contrast, when we focus on the variation of the shear amplitude = 0.05, 0.2, 0.4, 0.4, 0.6, 0.8, the shear frequency is fixed at ˝ = 0.002. The initial distribution of  is specified by random numbers uniformly distributed in the range [−0.01, 0.01] around the average value 0 = 0, corresponding to a completely disordered uniform state. And the initial membrane shape is supposed to be ˜ perfectly flat which indicates h(n, 0) = 0. By especially monitor˜ < 0.05 will be ing variations of the height field, the condition ∇ h successfully satisfied.

4. Results and discussions 4.1. Dynamical evolution of phase morphologies Due to the fact that the average strain would be zero in each time period of shearing, differently from the case with steady shear,

the effect of oscillatory shear on the morphology of domains will not be so easy to understand intuitively. In this subsection, we present the visualized results of our simulation via snapshot figures of the membrane shape. As is shown in Fig. 2, the four snapshots in (A) enable us to follow the whole time history of the phaseseparation process without external shear flow. In this case, mainly by force of the thermodynamic potential, the system is quenched from a disordered phase state into a two-phase regime (the saddlebacking represents A-rich phase, the valley denotes B-rich phase). With time going on, domains become more and more coarsened and an isotropic bicontinuous structure emerges in the late stage of the phase-separation process. In Fig. 2(B), external oscillatory shear flow is applied, with a fixed shear amplitude = 0.2, different shear frequencies are compared. Membrane shapes at = 100000 for each frequency indicate that the larger ˝ is, the more obvious domains elongate in the direction perpendicular to the flow. If the frequency is high enough, a distinct anisotropic morphology will be produced. In the present model, the instantaneous strain is behaving as a sinusoidal function of time . When the frequency is high enough, the convection caused by the flow will change very fast between positive and negative cyclically. Therefore, unstable fluctuations could not be convected away in the flow direction and the domains will grow mainly along the velocity gradient direction. Since the gradient of the strain is in the y direction, it displays that in the region near y = 0 the domain alignment is very weak, leading to an approximately axisymmetrical structure. Interestingly, this is totally different from the case when a strong steady shear flow is imposed on such a phase-separating system. Actually in that case, string-like domains may completely elongate along the direction of the flow throughout the whole region [35–37]. Then in Fig. 2(C), with a fixed shear frequency ˝ = 0.002, we consider shear amplitudes of varied values. From the three images upside,

˜ (A) Time evolution of the phase morphology for = 0, ˝ = 0, i.e. without external shear flow. (B) The Fig. 2. Snapshot figures of the phase morphologies (performed by h). membrane shape at = 100,000 for different shear frequencies with = 0.2. (C) The membrane shape at = 100,000 for different shear amplitudes with ˝ = 0.002.

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

we can see that when is not too large, at the same moment = 100,000, larger shear amplitude will mean more coarsened net regime. When is set to be appropriately large, domain morphologies oriented in the direction perpendicular to the flow will also be

207

observed, just like the case of high frequency. In the two figures underside, the shear amplitude is so large that in this late stage string-like domains in most regions exhibit a quite straight and slim feature. Yet, in the region near y = 0, because of the stacking induced

Fig. 3. Configurations of structure factor C(q, ) in the q space at different moments for typical cases. (A) Without external shear flow. (B) Different shear frequencies with

= 0.2. (C) Different shear amplitudes with ˝ = 0.002.

208

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

by the oscillatory shear flow, the domains are in much larger scale. 4.2. Simulated scattering patterns It is indeed necessary to more quantitatively discuss about this kinetic ordering process. We would like to stress that for the reason that the field h˜ displays a same evolution law as the field  does [17], from this section on, our discussions will only focus on the field . Here firstly the discrete Fourier transform of (n,t) would be take into account, which is defined as q ( ) =

n

(n, ) exp(iq · n)

(24)

with the wave-number q = 2 n/L and n ∈ {0, 1, . . ., L − 1}2 . Based on this treatment, the structure factor C(q, ), which is the main observable parameter for the description of the phase-separation kinetics, will be determined by C(q, ) = q ( )−q ( ), where the average is over the ensemble of the system. Being the directly obtained result of small-angle light scattering experiment, the structure factor is in fact the “power spectrum” of the field (n, ) in reciprocal space. Scattering patterns performed by C(q, ) is approved to reveal statistical trends and correlations, despite of the spatially stochastic two-phase constitution at specified instant and location. As is seen in Fig. 3, when there is no external shear flow imposed, the scattering pattern keeps its isotropic feature as a circular ring all through the whole phase-separation process, with the corresponding radius becoming smaller along with the increasing time,

Fig. 4. Plots of time evolution of the characteristic domain sizes for each case in double logarithmic scale. (A) Without external shear flow. (B) Different shear frequencies with = 0.2. (C) Different shear amplitudes with ˝ = 0.002.

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

209

Fig. 4. (Continued ).

which indicates that the maximum of the scattering intensity shifts to smaller wave-numbers. Next, we turn to the comparison of shear flows with different shear frequencies and a fixed amplitude

= 0.2. It appears that when ˝ is not large enough, saying not larger than 0.002, the scattering ring could be preserved but transforms to be an ellipse, with the symmetry axes rotating. Then, as ˝ is set large, in the early evolution stage, elliptical ring with its symmetry axes rotated can be clearly observed, whereas such rings will be destructed later. For ˝ = 0.01, the scattering pattern finally turns into two isolated speckles. For large frequencies, we note that the patterns at adequately late stage are locating extremely near the qx axis, which actually corresponds to the domain elongation in the direction perpendicular to the flow. In fact, this interesting relationship between the morphology and the scattering pattern is in consistent with the experimental result of Hashimoto et al. [35], where the two isolated scattering speckles are along the direction perpendicular to the elongated string phase induced by a steady shear flow. In Fig. 3(C), in respect of the fundamental influencing mechanism of such external oscillatory shear flow, the phenomena of being shaped elliptical, rotating symmetry axes, ring destruction as well as the eventual double-speckles would be observed undoubtedly. For each case of large enough , since the shear frequency is set the same, the angle between the symmetry axis and the qx axis at a typical moment will be of no difference correspondingly. Further, though has not been shown out, it could also be found

that for each individual evolution process, the maximum of C(q, ) keeps in step with the increasing . This is just in accordance with the coarsening of the morphology. 4.3. Time growth law of characteristic domain sizes The anisotropic phenomenon reflected by the domain elongation implies the systematic existence of two characteristic scales in the size distribution of domains, namely in the two orientations parallel and perpendicular to the imposed flow. Here in this subsection, from the knowledge of the structure factor C(q, ), we are allowed to compute the characteristic domain sizes in both directions with the following formulas Rx ( ) =



where q2x ( ) =

2 q2x ( )

,

Ry ( ) =

 2 q C(q, ) dq x , C(q, ) dq



2

(25)

q2y ( )

q2y ( ) =

 2 q C(q, ) dq y C(q, ) dq

(26)

In the absence of external shear flow, as is seen in Fig. 4(A), it is quite clear that the Rx ( ) and Ry ( ) curves are superposed, indicating the isotropic phase evolution throughout the whole process in accordance with the morphology snapshots shown in

210

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

Fig. 2(A). Actually in this case, thermodynamics-induced diffusion is the main physical mechanism responsible for phase-separation process. Apparently, a critical time point c ≈ 3000 classifies the phase-separation event into a short early stage with a very gentle growing speed and a later stage with an exponential growth law. This crossover behavior can be explained by a deterministic surface diffusion mechanism [38,39]. In the later stage, there is a remarkable time evolution law Ri ( ) ∼ 0.24 (i represents x,y). It is clear that the power law exponent is smaller than 1/3 for diffusive growth predicted by the Lifshitz–Slyozov theory, which has been reported in many previous studies [25,39]. Firstly, in this model, phase coarsening is a result of the strong coupling between the diffusion-and-coalescence growth mechanism and the elastic nature of the membrane. The curvature can considerably slow down the coarsening segregation process [18]. Secondly, a basic assumption of our model is that the two membrane components favor different intrinsic membrane curvatures, which is controlled ˜ (see ˜ It will mean to enhance the parameter ˇ by the parameter . Eq. (12)) which actually acts to prevent the system from phase separation. When the shear flow is applied, as is mentioned in the foregoing analyses, here the growth of characteristic domain sizes also becomes more and more anisotropic. When the external shear effect is enhanced, meaning that either the frequency or the amplitude is increased, a much earlier starting moment of the shear effect could be predicted. Judge from Fig. 4(B), it is evident that when the frequency is too small (say the case ˝ = 0.0002), in the previously mentioned later stage, the contrasting relationship between Rx ( ) and Ry ( ) curves behaves quite different from cases with larger frequency. The reason is that in each semi-period of the oscillatory shear flow, the shear effect is very similar to a steady shear flow,

which may generally reinforces the domain coarsening in the flow direction and inhibits that in the perpendicular direction [36]. Then, it is clearly observed that when ˝ is large enough, Rx ( ) curve will stay upon Ry ( ) curve, and the larger ˝ is, the further the two curves depart away from each other. As is depicted in each image, the curve slopes actually represent exponents of the growth laws Ri ( ) ∼ ˛i . It should be noted that in these cases with external oscillatory shear, such growth law just expresses the evolution tendency, whereas the actual evolution curve is indeed oscillatory. Significantly, as ˝ becomes larger, the phenomenon of decreasing ˛x and increasing ˛y in each stage implies that the applied external oscillatory shear flow suppresses the domain coarsening in the flow direction and reinforces that in the perpendicular direction. In Fig. 4(C), with the same ˝ = 0.002, there is no doubt that when is set larger, the influence of the external shear flow on the domain growth will be more evident, just as the plots show. Interestingly, in the considered range of , concerning the later stage, ˛y stands at 0.24 for all cases while ˛x falls for cases with sufficiently large . Another important point is that the oscillation amplitude of Ry ( ) curve is keeping its step with the increasing , yet Rx ( ) curve does not perform a so conspicuous oscillation for all cases of . 4.4. Rheological behavior of such membrane system The domain microstructures which depend on their thermal and deformation history conditions can lead to diverse rheological properties of the membrane system, which will be of fundamental importance for many biological applications. Related to the presence of the flow, we would like to turn to the study of the behavior of the stress response induced by the interface in such

Fig. 5. Time dependence of the normal stress difference N1 . (A) Different shear frequencies with = 0.2. (B) Different shear amplitudes with ˝ = 0.002.

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

211

Fig. 6. Storage modulus G and loss modulus G : (A) shear frequency dependence for storage and loss moduli at a fixed shear amplitude = 0.2. (B) Shear amplitude dependence for storage and loss moduli at a fixed shear frequency ˝ = 0.002.

a phase-separating system. From the rheological viewpoint, the excess stress tensor is introduced as [40,41]

G =

 e =− xx

1

|q|≤qc

 e yy



|q|≤qc

q2x C(q, ) dq,

(2 )d 1

=−

e xy =−

(2 )

G =

q2 C(q, ) dq, d y

1

|q|≤qc

qx qy C(q, ) dq

(2 )d

(27)

With this, here the normal stress difference is given by

 N1 =

where

e xx

e − yy

= |q|≤qc

1 (2 )d

(q2y − q2x )C(q, ) dq

(28)

where qc is a phenomenological cutoff, d is the spatial dimension and equals to 2 in the present model. We show in Fig. 5 the time evolution of the normal stress difference N1 . As can be observed in (A), with no exception for all cases of ˝, in the early stage N1 exhibits a decline from zero till a minimum. And after that, a gentle rise is quite clear before N1 begins to perform a steady oscillation, keeping in step with the external oscillatory shear. It is worth to point out that for ˝ = 0.0002 and 0.0005, in the later stage the curves are staying above zero while for large ˝ it will be totally under zero. The same reason which has been explained in the previous subsection will be responsible for this result, namely that oscillatory shear with a too small ˝ may produce an analogous effect as a steady shear. Then in (B) (for image clarity, = 0.8 is not plotted), it is certain that the substratal evolution discipline is the same as that in (A). However, for large enough shear amplitudes, taking = 0.6 for example, we obviously find that before N1 gets into a steady oscillation, there is a pseudo wave-crest within each period of the external oscillatory shear (see the interval between the two vertical dashed lines). This should be attributed to the interaction and competition between the diffusion mechanism and the flow convection in such stage. So far, the most useful way of characterizing such complex fluids is to investigate rheological properties by means of the linear viscoelastic experiment. When the storage modulus and loss modulus are directly measured, all the other linear viscoelastic properties can be calculated from them. It is customary to introduce the complex shear modulus G* , being expressed as [42,43] G∗ = G + iG

(29)

1

T 1

T



T e xy ( ) sin(˝ ) d ,

0



T e xy ( ) cos(˝ ) d

(30)

0

In the above expressions, G is the elastic contribution to the complex modulus, which is a measure of stored energy without phase difference between the shear stress and the shear strain, and represents the elastic component of the material. G is the loss or viscous modulus, which is a measure of the energy per unit volume dissipated by the material as heat, and represents the viscosity of the materials. According to Fig. 6, first of all it is noted that for all cases G is much larger than G , which indicates that after such a phaseseparation process the system behaves to be elastic more than viscous. It is of interest that for G and G , the dependences of shear frequency and shear amplitude are all exhibiting a linear relationship, three of which are up-climbing except the slightly descending G ∼ ˝ line. 5. Conclusions In this study, solving the modified two-dimensional TDGL equations by means of CDS approach, we have modeled the phase separation in binary lipid membrane under an external oscillatory shear flow. It is revealed that in comparison to the case without any external shear effect, the presence of such an oscillatory shear flow will significantly alter the phase morphology and the domain growth behaviors, although the average deformation is zero for each period of the oscillatory shear. We have studied the effects of the applied flow by changing the frequency at a fixed shear amplitude = 0.2, and also by changing the amplitude at a fixed shear frequency ˝ = 0.002. For the case of very small amplitude and frequency, within the time rang of our simulation, it can be easily imagined that the domain structure will be still isotropic in appearance, just as the case without external shear flow. When either the amplitude or the frequency is increased, microphase stripes will emerge, orienting along the direction perpendicular to the external shear flow. It is totally different from the case of steady shear flow where the string-like phase elongates along the flow direction. Corresponding to the phase morphology, the scattering pattern defined by structure factor displays a very intuitionistic and interesting timedependent evolution. More exactly, such an external oscillatory shear flow has mainly induced two impacts on the scattering pat-

212

X.-B. Chen et al. / Colloids and Surfaces B: Biointerfaces 65 (2008) 203–212

terns. On one hand, the shape alters from an original circular ring to an elliptical ring, which then destructs or even degenerates into two isolated speckles in strong shear regime. And the other impact is the rotation of the symmetry axes. To more precisely analyze the role of the external oscillatory shear flow in the phase separation, the time growth laws of the characteristic domain sizes in both directions are concerned. It is shown that such a shear flow (not too weak) would mainly reinforce the domain coarsening in the direction perpendicular to the flow and inhibit that in the parallel direction. Furthermore, rheological properties interrelated to the external shear are taken into account, including the time evolution of the normal stress difference for each case, and the shear frequency dependence as well as the shear amplitude dependence of both storage modulus and loss modulus. As far as we are aware of, scarcely any experiments reporting these features of such binary lipid membrane system with externally imposed oscillatory shear flow have been carried out. Consequently, it would be interesting and necessary to devise an experimental setup for testing and supporting the computer simulated results. It is expected that the present simulation might provide beneficial suggestions for studying and controlling the domain structure in binary lipid membrane and even in veritable biological membrane, by suitably imposing an external oscillatory shear flow. Acknowledgements

[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

We greatly appreciate the financial support from the Chinese National Natural Science Foundation No. 10672086 and Special Funds for the Major State Basic Research Projects 2004CB619304. References [1] M.B. Sankaram, T.E. Thompson, Biochemistry 29 (1990) 10670. [2] K. Simons, E. Ikonen, Nature 387 (1997) 569.

[35] [36] [37] [38] [39] [40] [41] [42] [43]

D.A. Brown, E. London, Annu. Rev. Cell Dev. Biol. 14 (1998) 111. J.R. Sevinsky, L.V.M. Rao, W. Ruf, J. Cell Biol. 133 (1996) 293. E. Stang, J. Kartenbeck, R.G. Parton, Mol. Biol. Cell 8 (1997) 47. F. Gisou van der Goot, H. Thomas, Semin. Immunol. 13 (2001) 89. D.A. Brown, E. London, J. Biol. Chem. 275 (2000) 17221. S. Mayor, M. Rao, Traffic 5 (2004) 231. S.J. Singer, G.L. Nicolson, Science 175 (1972) 720. D.A. Brown, E. London, J. Memb. Biol. 164 (1998) 103. D.A. Brown, Proc. Natl. Acad. Sci. U.S.A. 98 (2001) 10517. A.V. Samsonov, I. Mihalyov, F.S. Cohen, Biophys. J. 81 (2001) 1486. C. Dietrich, Z.N. Volovyk, M. Levi, N.L. Thompson, K. Jacobson, Proc. Natl. Acad. Sci. U.S.A. 98 (2001) 10642. P.B.S. Kumar, M. Rao, Phys. Rev. Lett. 80 (1998) 2489. M. Laradji, P.B.S. Kumar, Phys. Rev. Lett. 93 (2004) 198105. E.J. Wallace, N.M. Hooper, P.D. Olmsted, Biophys. J. 88 (2005) 4072. H. Wada, J. Phys. Soc. Jpn. 72 (12) (2004) 3142. R. Reigada, J. Buceta, K. Lindenberg, Phys. Rev. E 71 (2005) 051906. T. Nakajima, K. Kon, N. Maeda, K. Tsunekawa, T. Shiga, Am. J. Physiol. Heart Circ. Physiol. 259 (1990) H1071. C.H. Wang, A.S. Popel, Math. Biosci. 116 (1993) 89. S. Chien, S. Usami, R.J. Dellenback, M.I. Gregersen, Am. J. Physiol. 219 (1970) 136. J.L. Yi-Shuan, H.H. Jason, S. Chien, J. Biomech. 38 (2005) 1949. A.J. Bray, A. Cavagna, R.D.M. Travasso, Phys. Rev. E 64 (2001) 012102. K.A. Smith, W.E. Uspal, J. Chem. Phys. 126 (2007) 075102. Y. Oono, S. Puri, Phys. Rev. A 38 (1988) 434. W.T. Gozdz, G. Gompper, Europhys. Lett. 55 (4) (2001) 587. U. Seifert, Adv. Phys. 46 (1997) 13. O. Takao, E. Yoshihisa, Macromolecules 26 (1993) 4928. S. Ramaswamy, Phys. Rev. Lett. 69 (1) (1992) 112. S. Puri, Y. Oono, Phys. Rev. A 38 (1988) 1542. W. Helfrich, Z. Naturforsch. C 28 (1973) 693. J.L. Harden, F.C. MacKintosh, P.D. Olmsted, Phys. Rev. E 72 (2005) 011903. V.D. Gordon, P.A. Beales, Z. Zhao, C. Blake, F.C. MacKintosh, P.D. Olmsted, M.E. Cates, S.U. Egelhaaf, W.C.K. Poon, J. Phys.: Condens. Matter 18 (2006) L415. N. Kahya, D. Scherfeld, K. Bacia, B. Poolman, P. Schwille, J. Biol. Chem. 278 (30) (2003) 28109. T. Hashimoto, K. Matsuzaka, E. Moses, A. Onuki, Phys. Rev. Lett. 74 (1995) 126. E.N.M. Cirillo, G. Gonnella, G.P. Saracco, Phys. Rev. E 72 (2005) 026139. Z.Y. Shou, A. Chakrabarti, Phys. Rev. E 61 (2000) R2200. H. Furukawa, Adv. Phys. 34 (6) (1985) 703. S. Puri, Y. Oono, J. Phys. A: Math. Gen. 21 (1988) L755. A. Onuki, Phys. Rev. A 35 (1987) 5149. F. Corberi, G. Gonnella, A. Lamura, Phys. Rev. E 61 (2000) 6621. ´ Phys. Rev. E 53 (1996) 2588. H. Komatsugawa, S. Nose, J.C. van der Werff, C.G. de Kruif, C. Blom, J. Mellema, Phys. Rev. A 39 (1989) 795.