Simulations of the spreading of a vesicle on a substrate surface mediated by receptor–ligand binding

ARTICLE IN PRESS

Journal of the Mechanics and Physics of Solids 55 (2007) 1166–1181 www.elsevier.com/locate/jmps

Simulations of the spreading of a vesicle on a substrate surface mediated by receptor–ligand binding P. Liua, Y.W. Zhangb,, Q.H. Chenga, C. Lua a Institute of High Performance Computing, Singapore 117528, Singapore Department of Materials Science and Engineering, National University of Singapore, Singapore 119260, Singapore

b

Received 8 May 2006; received in revised form 24 November 2006; accepted 2 December 2006

Abstract A continuum model was introduced for the adhesion of vesicles to substrate surfaces. In the model, the vesicle membrane was assumed to be a closed shell with hyperelasticity. The vesicle cavity is filled with a liquid of fixed volume. The receptors on the membrane are mobile and initially uniformly distributed while the ligands on the substrate surface are fixed and also uniformly distributed. The formation of localized regions of tight binding between receptors and ligands, results in vesicle adhesion to the substrate surface. An adhesive model was introduced to describe the adhesive interaction between the receptors and the ligands. The growth of the adhesion area occurs via recruiting receptors from the non-adhered region through diffusion. Finite-element methods were used to solve the governing equations for the deformation of the vesicle and the receptor diffusion on the membrane surface. Effects of the membrane stiffness, the cohesive parameters and the receptor density on the adhesion kinetics of the vesicle were studied. In addition, the instability of the advancing front of the adhesion was also analyzed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Adhesion; Surface diffusion; Membrane; Vesicle; Finite elements

Corresponding author. Tel.: +65 6516 1542.

E-mail address: [email protected] (Y.W. Zhang). 0022-5096/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2006.12.001

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1167

1. Introduction The adhesion of cells to the extracellular matrix (ECM) is essential to many physiological and pathological processes, such as, neutrophils tracking down pathogens, fibroblasts crawling during wound healing, metastatic cancer cells invading distant parts of the body, formation of thrombosis due to the attachment of platelets on tissues, etc. (Bao and Suresh, 2003; Wolgemuth, 2005; Park et al., 1990). The mechanics and dynamics of cell adhesion and spreading are strongly dependent on receptor–ligand-mediated interactions between cells and the extracellular matrix (Bell, 1978; Bell et al., 1984). The concerted deformations of the cell membrane, membrane-attached proteins, cytoskeleton, ligand type, and the concentration and stiffness of the ECM are believed to affect cell adhesion and spreading (Doebereiner et al., 2004; Hategan et al., 2004; Giannone et al., 2004; Dubin-Thaler et al., 2004; Reinhart-King et al., 2005). Many experimental studies have been performed to investigate the process of cell spreading. However, the observations and conclusions are often inconsistent and detailed adhesion mechanisms and driving forces are still of interest (Dubin-Thaler et al., 2004; Doebereiner et al., 2004; Balaban et al., 2001; Tan et al., 2003; Beningo et al., 2001; Gallant et al., 2005). Recently, Reinhart-King et al. (2005) have used traction force microscopy to measure the forces exerted by bovine aortic endothelial cells (BAECs) during spreading on a flexible substrate. It was observed that cells are capable of exerting significant forces before complete actin polymerization or visible stress fiber formation. Their work seems to suggest that cell adhesion can be reasonably modeled by vesicles. By using vesicles as a biomimetic model for cells, the complexity of cells can be avoided; in vesicles, no cytoskeleton and nucleus effects come into play and also the membrane proteins can be well controlled (Seifert and Lipowsky, 1990; Kloboucek et al., 1999; Boulbitch et al., 2001). Significant theoretical and modeling effort has also been made to understand the energetics, kinetics, and dynamics of cell adhesion and spreading. A thermodynamic framework for cell adhesion mediated by receptor–ligand binding was established in the pioneering work of Bell (1978) and Bell et al. (1984). However, cell adhesion is a nonequilibrium process and can be better understood in terms of a nucleation and growth process. Several kinetic models have been proposed, for example, reaction kinetics and Monte Carlo simulations have been used to describe the binding of ligands to cell surfaces (Zhu, 2000; Zhu et al., 2000; Irvine et al., 2002). A multiscale approach for studying the adhesion of a cell onto a substrate has also been developed (N’Dri et al., 2003). A theoretical model of the motion of the adhesion front based on the receptor diffusion and binding and un-binding reaction was proposed to explain the time dependence of the displacement of the adhesion front on the receptor concentration (Boulbitch et al., 2001). This model was recently extended to include the mechanics of the cell membrane to derive a scaling relationship between the displacement of the adhesion front and the adhesion time (Freund and Lin, 2004; Shenoy and Freund, pffiffi 2005). These authors found that the adhesion front radius increases in proportion to t, where t is the adhesion time, and that the circular front becomes unstable under sinusoidal perturbation for adhesion front radii larger than the nucleation size. This model was subsequently used to study the mechanics of receptor-mediated endocytosis (Gao et al., 2005; Bao and Bao, 2005). Despite the great interest expressed in cellular adhesion and the plethora of experimental and theoretical studies performed thus far, the understanding of specific and non-specific binding in terms of adhesion strength and energy is still rudimentary at best. Van der

ARTICLE IN PRESS 1168

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

Waals forces between cells are expected to be attractive and of the long-range type. In addition, the repulsive forces due to the compression of the glycocalyx proteins provide the resistance to cell adhesion. Both van der Waals forces and the compressive forces are characterized as non-specific forces because they are expected to be present between all cells (Bell, 1978). Generally it is believed that cell adhesion is mediated by the specific binding of receptor–ligand pairs and these specific bonds can cause cells to adhere to each other much more tightly than the non-specific electrical forces (Bell, 1978). The membrane–membrane adhesion and separation problem was studied by Evans using a traction separation relation (Evans, 1985). The molecular cross-bridging forces were assumed to be normal tractions continuously distributed between the membranes. Traction separation relations were also used to model the interface decohesion of biomaterials in the past (Tvergaard and Hutchsinson, 1993; Liu et al., 2001). The interaction between the cell and the substrate may also be formulated using a traction separation relation: such an approach should also permit specific and non-specific interactions between the cell surface and the ECM surface to be modeled (Evans, 1985). In this paper, we extend the model developed by Freund and co-workers (Freund and Lin, 2004; Shenoy and Freund, 2005) to study the adhesion of a cell to a substrate surface. The interaction between the surfaces of the cell and the substrate is modeled by a traction separation relation. A fully three-dimensional model of a vesicle is used to mimic a real cell. The membrane of the vesicle is modeled as a thin hyperelastic shell and the cavity inside the vesicle is assumed to be filled by a fluid of fixed volume. Finite-element methods were used to solve the equations governing the mechanics of the cells, the receptor diffusion on the vesicle membrane, and the contact mechanics between the vesicle surface and the substrate surface. Parametric studies of the various system parameters such as membrane stiffness, adhesive parameters and receptor density were also performed to understand the mechanics and kinetics of cell spreading. The instability of the spreading front was also analyzed. 2. Model formulation A vesicle has a relatively simple structure in that its membrane comprises a cellular phospholipid bilayer, in which well-defined hydrophobic proteins can be inserted and retained, and which encloses a fluid of fixed volume and known viscosity. Consequently, a vesicle is often regarded as a biomimetic cellular system (Seifert and Lipowsky, 1990; Kloboucek et al., 1999; Boulbitch et al., 2001). Vesicle adhesion is schematically illustrated in Fig. 1, which shows the adhesion of a vesicle membrane containing receptors to a rigid substrate surface with a uniform and fixed ligand density. It was assumed in this study, that the receptors on the vesicle membrane were mobile. If the attractive forces between the receptors on the vesicle membrane and the ligands on the substrate were assumed to be large enough to overcome the repulsive forces such as steric forces by compressing glycoproteins and solvation forces by squeezing liquid molecules out of the contact gap, the receptor–ligand bonds can form, which cause a reduction of internal energy. The number of such receptor–ligand pairs formed should be proportional to both the receptor and ligand densities when the distances between the receptors and ligands are close enough. This is the case when the receptor density is lower than the ligand density. However, when the receptor density is higher than the ligand density, the excess receptors cannot form receptor–ligand pairs since all the ligands would have already paired up.

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1169

Fig. 1. Schematic of vesicle adhesion mediated by the diffusion of the receptors and the binding of the receptor–ligand pairs. n is the surface unit normal vector and t is the surface unit tangent vector.

Typically, the binder (receptor or ligand) density can be as high as 103–104 binders per square micron meter. Hence, the continuum approach was considered appropriate for biological adhesion. Instead of using the discrete ligand–receptor binding to describe the molecular interaction between the receptors and ligands, the traction and separation relation was used. Here two major parameters were required to describe the adhesion behavior, that is, the adhesion traction, which is the maximum traction required to separate a unit area of the adhesive interface, and the adhesion energy, which is the energy required to separate a unit area of the adhesive interface. In the present work, a traction–separation relation was proposed to describe the adhesive behavior of the membrane receptors and the substrate ligands. It was assumed that the receptors on the vesicle membrane were attracted by the ligands on the substrate surface by the following traction: 8 0; d4dc <   T¼ for cr ocl (1a) d : Hcl cr 1  dc ; dpdc and T¼

8 <



0;



2 d : Hcl 1  dc ;

d4dc dpdc

for

cr Xcl ,

(1b)

where cl is the ligand density on the substrate surface, cr is the receptor density on the vesicle surface, H is a parameter related to the pair interaction between a receptor and a ligand, d is the distance between the receptor and substrate surface, and dc is the cutoff interaction distance beyond which no interaction occurs between a receptor–ligand pair. The maximum adhesive traction is given by T max ¼ Hcl cr . Hence, the maximum adhesive energy is given by E max ¼ T max dc =2. In the present formulation, the force acting on a receptor that is embedded in the vesicle membrane gives rise to two consequences: (1) the force pulls the membrane towards the substrate surface, creating adhesion, and (2) the force attracts the receptors and cause them to diffuse along the membrane wall since receptors are mobile within the bilayer

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1170

membrane. For the receptor diffusion on the membrane surface without considering the adhesion force, the diffusion flux of the receptors on the membrane surface is jc ¼ Drs cr ,

(2)

where D is the receptor diffusivity on the membrane and rs is the surface gradient operator (Weatherburn, 1925). It was assumed that the diffusion flux due to the attractive traction from the ligands on the substrate surface is proportional to the traction tangential to the membrane surface, i.e., jF ¼ MTt3 t, where t is the tangent vector of the membrane surface and t3 is the z-component of t (see Fig. 1), M is the receptor mobility on the membrane surface under the attractive force due to the ligands, and T is the adhesive traction in the negative z-direction. Thus, the total diffusion flux can be written as j ¼ jc  jF .

(3)

Since the number of receptors on the lipid membrane is conserved, the governing equation for receptor density is qcr ¼ Dr2s cr þ rs  jF . qt

(4)

In this study, it was further assumed that once the membrane is in contact with the substrate, the receptors in the contact area become immobile, i.e., receptor diffusivity was set to be zero. The vesicle bilayer membrane is commonly modeled as an incompressible medium (Lim et al., 2004; Dao et al., 2003). The large deformation taking place during vesicle adhesion means that a hyperelastic material model needs to be employed. In our study, the incompressible neo-Hookean hyperelastic material model was chosen. The strain energy potential function is given as (Dao et al., 2003), U¼

 G 2 l1 þ l22 þ l23  3 , 2

(5)

where G is the initial bulk shear modulus and li (i ¼ 1, 2, 3) are the principal stretches. In previous studies (Lim et al., 2004; Dao et al., 2003), the effect of membrane viscosity on cell deformation was found to be negligibly small. Thus, the viscoelasticity of the membrane was ignored in our study. The vesicle cavity is modeled as an incompressible fluid. The key issue here is the modeling of the coupling between the deformation of the membrane and the pressure exerted by the enclosed fluid. To simplify the problem, the fluid in the cavity was assumed to be homogenous in composition and density and also uniform in terms of pressure and temperature. Therefore, the only parameter describing fluid behavior was the reference fluid density rr, which was specified at zero pressure and room temperature. The viscous energy dissipation of the cytoplasm is also known to be small (Evans and Hochmuth, 1976; Dao et al., 2003). The viscosity of the enclosed fluid can therefore be ignored. During adhesion, the membrane will make contact with the substrate surface. For simplicity, it was assumed that the substrate was rigid. It should be noted that during contact, interpenetration of the membrane and the substrate surface is not permitted. Furthermore, the normal traction of the membrane must be balanced by the traction on the substrate surface. Since the membrane is relatively fluidic within the membrane plane, it was assumed that there was no friction between the membrane and the substrate surface.

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1171

3. Numerical procedure A finite-element procedure was designed to solve the governing equations outlined in the preceding section. The procedure involves three steps. Consider the vesicle adhering to a substrate as shown in Fig. 1. Suppose that the shape of the vesicle and receptor and ligand densities at time t are known. Our objective is to calculate the deformation and diffusion that occurs during a subsequent infinitesimal time interval Dt. The first step involves calculating the traction exerted by the ligands on the substrate surface to the membrane at time t using the adhesive model described in the preceding section. In the second step, the deformation of the vesicle under the ligand tractions is computed, while ensuring that the contact condition with the substrate surface is satisfied. This step of the analysis is performed using a commercially available general purpose finite-element package ABAQUS (ABAQUS User’s Manual, 2002). It should be noted that the mechanical equilibrium equations were solved using the configuration at time t as the reference configuration. To predict the mechanical response of the fluid-filled cavity, hydrostatic fluid elements were used. To ensure the stability of the solution, a small damping factor was introduced in the calculations. In the third step, the diffusion equation to determine the change in the receptor density at time t þ Dt is solved. The traction due to the ligands on the substrate surface with the new receptor density can be obtained. This procedure is repeated to compute the change in the shape of the vesicle and the adhesion area as a function of time. Since the mechanical analysis of the deformation of the vesicle (both the membrane and the filled fluid) and its contact with the substrate surface was performed using ABAQUS, the detailed numerical procedure can be found in ABAQUS manuals and will not be repeated here. For the diffusion of the receptors, a finite-element procedure was developed to solve the surface diffusion (Eq. (4)). Let dcr be an admissible variation in the receptor density. Then this equation can be written in a weak form as Z Z  2  qcr dcr dS ¼ Drs cr þ rs  jF dcr dS. (6) qt S S By using the surface divergence theorem (Weatherburn, 1925), one obtains Z Z   qcr dcr dS ¼  Drs dcr  rs cr þ rs dcr  jF dS. S qt S

(7)

In deriving Eq. (7), the relations rs cr  n ¼ 0 and t  n ¼ 0 were used. Here n and t are the unit normal and tangent vectors of the membrane, respectively (Fig. 1). To enhance the stability and accuracy, a semi-implicit Euler scheme was used to integrate Eq. (7). In order to calculate the receptor density of the membrane surface, a mesh of triangular finite elements was generated to model the membrane (Zhang et al., 1999). The values of the system parameters were mostly obtained from previous work (Dao et al., 2003; Freund and Lin, 2004; Shenoy and Freund, 2005; Gao et al., 2005, Peeters et al., 2005). For all calculations, the following values of the system parameters were chosen: the receptor diffusivity D ¼ 102 mm2 =s, the mobility M ¼ 5  105 mm=s, and the ligand density cl ¼ 104 =mm2 . The initial configuration is chosen as a spherical vesicle with a uniform distribution of receptors on its surface making a point contact with a flat substrate covered by a uniform and fixed distribution of ligands.

ARTICLE IN PRESS 1172

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

Our simulation results were reported using the following normalization scheme: length was normalized by the diameter of the vesicle, i.e., a0 ¼ D0 ; time by t0 ¼ a20 =D, energy by E 0 ¼ kB T R , where kB is the Boltzmann constant and TR is the room temperature; force by F 0 ¼ E 0 =a0 ; and stress by s0 ¼ F 0 =a20 . 4. Results and discussions The finite-element procedure outlined in the preceding section was developed with an aim to predict the cell adhesion evolution through the adhesive interaction between diffusible receptors on the vesicle surface and the immobile ligands on the substrate surface. Parametrical studies were performed to understand the adhesion kinetics and its dependence on the system parameters. 4.1. Spreading kinetics for a typical case For the typical case, the normalized cutoff distance is dc ¼ 0:1, the initial receptor density is cr0 ¼ cl =8, and the normalized initial maximum traction is T max ¼ 2:416, and the initial in-plane membrane shear modulus, which is defined as the product of the membrane shear modulus and membrane thickness, is K  ¼ 1:81 after normalization. The spreading snapshots are shown in Fig. 2(a)–(d), in which the contour plots of the receptor density distribution are also given. The vesicle starts from the initial state and gradually spreads by increasing its binding area with the substrate surface. The curve of the binding area vs. time during the above spreading process is shown in Fig. 3. Initially the spreading is fast since the conformation of the membrane to the substrate can be easily made by slightly bending the membrane. Once passing this initial stage, a sub-linear relationship between the binding area and time is observed. The spreading finally stabilizes after the sub-linear spreading regime, at a configuration as shown Fig. 2(d). This final state can only be achieved when the adhesive force-driven receptor diffusion flux is balanced by the receptor concentration gradient-driven diffusion flux, that is, no net flux occurs on the vesicle surface. Our calculation indeed showed a negligible diffusion flux on the vesicle surface at the final stage (Fig. 2d). It was noted that the number of receptors is conserved on the vesicle membrane and this eventually leads to the transition from the spreading state to the stable state configuration. 4.2. Parametric studies of system parameters 4.2.1. The effect of initial in-plane membrane shear modulus Different cells may have different membrane molecular structures and thus exhibit different membrane mechanical properties. As a result, the shear modulus and thickness of cell bilayer membranes may be different. The insertion and retaining of membrane proteins and the support from the spectrin network and the cytoskeleton make it difficult to accurately determine the membrane shear modulus and membrane thickness. Here, a combined quantity, that is, the initial in-plane membrane shear modulus, was considered. Parametric studies were performed to investigate its effect on the spreading kinetics. In the calculations, dc ¼ 0:1, cr0 ¼ cl =8 and T max ¼ 1:51. The following four values of the initial in-plane membrane shear modulus were chosen: K  ¼ 0:603, 1.206, 1.810 and 3.018. The curves of the binding area vs. time for the four cases are shown in Fig. 4. It can be seen that

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

a

1173

b 4.0 3.6 3.2 2.8 2.6 2.0 1.6 1.2 0.8 0.4

4.0 3.6 3.2 2.8 2.6 2.0 1.6 1.2 0.8 0.4 t*=0.008

t*=0.044

c

d 4.0 3.6 3.2 2.8 2.6 2.0 1.6 1.2 0.8 0.4

4.0 3.6 3.2 2.8 2.6 2.0 1.6 1.2 0.8 0.4 t*=0.560

t*=0.128

Fig. 2. Snapshots of the vesicle spreading and receptor density contours for a typical case (dc ¼ 0:1, cr0 ¼ cl =8, T max ¼ 2:416, K* ¼ 1.81). (a) t* ¼ 0.008; (b) t* ¼ 0.044; (c) t* ¼ 0.128, and (d) the final state at t* ¼ 0.560.

Contact Area(a20)

0.8

0.6

0.4

0.2

0.0

0

0.3

0.6 Time(t0)

0.9

1.2

Fig. 3. The curve of binding area vs. spreading time for the typical case (dc ¼ 0:1, cr0 ¼ cl =8, T max ¼ 2:416, K* ¼ 1.81). A sub-linear relation is observed.

ARTICLE IN PRESS 1174

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

Fig. 4. The curves of binding area vs. spreading time at different initial in-plane membrane shear moduli K* (dc ¼ 0:1, cr0 ¼ cl =8, T max ¼ 1:51).

the variation of the in-plane membrane shear modulus affects the spreading velocity of the adhesion front and the final adhesion area. The spreading velocity slows down with an increase in the in-plane shear modulus, and the adhesion area decreases. This result is expected since the increase in the in-plane shear modulus increases the elastic energy required to bend the membrane to conform to the substrate, thus increases the resistance to spreading and reduces the final adhesion area. 4.2.2. The effect of adhesive traction The adhesive traction varies with the ligand density, receptor density and the interaction of the ligand–receptor pair. Here the change of the initial maximum adhesion traction T max ¼ Hcl cr0 is investigated by changing H. In the calculations, dc ¼ 0:1, cr0 ¼ cl =8, and K  ¼ 1:81. The following values of the initial maximum adhesion traction were chosen: T max ¼ 0:302, 0.604, 1.510, and 2.416. The curves of the binding area vs. time for the four cases are shown in Fig. 5. It can be seen that the variation of the initial maximum adhesion traction also affects the spreading velocity of the spreading front and the final adhesion area. Both the spreading velocity and the final adhesion area increase with an increase in the initial maximum traction. Its influence is significant when the dimensionless initial maximum traction is lower than 1.51. However, when the dimensionless initial maximum traction is larger than 1.51, its influence becomes insignificant. This can be explained by the balance between the diffusion flux driven by the concentration gradient and the diffusion flux driven by the adhesive traction. Adhesion spreading is limited by the adhesion-driven diffusion when the adhesive traction is small. Hence, an increase in the adhesive traction can effectively increase the receptor flux to the adhesive front since the reduction in receptor density in the adhesive zone can be quickly compensated by the diffusion flux from outside the adhesive zone driven by the concentration gradient. However, when the

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1175

Fig. 5. The curves of binding area vs. spreading time at different initial maximum adhesion tractions T max (dc ¼ 0:1, cr0 ¼ cl =8, K  ¼ 1:81).

adhesive traction is too large, the diffusion flux from outside the adhesive zone cannot compensate the reduced receptor density in the adhesive zone due to the large adhesioninduced diffusion. As a consequence, a further increase in the adhesion traction will not increase the adhesion velocity due to the limited diffusion driven by the concentration gradient. 4.2.3. The cutoff distance As discussed in the introduction, van der Waals forces between cells are attractive and of the long-range type, and the compression of the glycocalyx layer provides repulsive forces for the intermediate range. These two types of forces are non-specific. On the other hand, specific binding between receptor–ligand pairs causes cells to adhere to substrates much more tightly than the non-specific forces. This type of binding force is strong and yet short ranged. In the present adhesive model, the net effect of these forces is assumed to be attractive and a cutoff interaction distance was introduced for the attraction. Since the exact value of the cutoff interaction distance is unknown, a parametric study was performed. In the calculations, cr0 ¼ cl =8, T max ¼ 1:51 and K* ¼ 1.81. The following three values of the cutoff distance were chosen: dc ¼ 0:06, 0.10 and 0.18. The curves of the contact area vs. spreading time for the three cases are shown in Fig. 6. It is seen that with an increase in the cutoff distance, the spreading rate increases slightly. This is due to the fact that an increase in the cutoff distance results in a linear increase in the adhesion energy. In addition, an increase in the cutoff distance increases the adhesive zone size. As a consequence, more receptors will be attracted by the ligands on the substrate. However, it was found that the final binding area is only weakly dependent on the variation of the cutoff distance when the initial maximum adhesion traction is held fixed. This result suggests that the influence of adhesion energy on vesicle adhesion is much weaker than that of adhesion traction.

ARTICLE IN PRESS 1176

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

Fig. 6. The curves of binding area vs. spreading time at different cutoff separations dc (cr0 ¼ cl =8, T max ¼ 1:51, K  ¼ 1:81).

Fig. 7. The curves of binding area vs. spreading time at different ratios of ligand density to receptor density (dc ¼ 0:1, K  ¼ 1:81, and T max ¼ 2:416). The ligand density has been adjusted to give a constant initial maximum traction T max .

4.2.4. The receptor density Experimental results (Boulbitch et al., 2001) showed that a change in the receptor density can also influence spreading kinetics. The present model was also used to study this influence. In the present calculations, dc ¼ 0:1, K  ¼ 1:81, and T max ¼ 2:416. This means that the initial maximum traction was fixed and the ratio of the ligand density and receptor density was varied. The following four values of the ratio were chosen: cl =cr0 ¼ 1:5, 2.0, 3.0 and 8.0. The curves of the contact area vs. spreading time for the four cases are shown in Fig. 7. For all the cases, a sub-linear relationship between the spreading area and time may

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1177

Fig. 8. Schematic diagram of the configurations of the spreading front. (a) an earlier stage; (b) a later stage. It is noted that with increasing spreading time, the radius of curvature gradually decreases while the spreading angle gradually increases.

be observed. It is seen that a larger ratio gives rise to a more linear relation while a smaller ratio gives rise to a more sub-linear relation. Steady-state solutions based on an infinite medium (Boulbitch et al., 2001) showed that if the adhesion process is limited by the receptor diffusion, that is, if cl =cr0 is much larger than 1, the spreading area increases linearly with the spreading time; while if the adhesion process is limited by the formation of receptor–ligand pairs, that is, if cl =cr0 is close to 1, the spreading area increases with the square root of the spreading time. In the present simulations, the linear relation is not observed for even a large value of cl =cr0 . This deviation is probably due to the finite size effect. During the vesicle spreading, the radius of curvature of the vesicle membrane at the spreading front decreases gradually as shown in Fig. 8(a)(b). At the same time, the spreading angle as shown in Fig. 8(a)(b) also increases continuously. Thus the steady-state assumption is no longer valid. The gradual decrease in the radius of curvature and the gradual increase in the spreading angle hinder the spreading. As a result, a sub-linear relation between the spreading area and time is observed. The tendency that a decrease in cl =cr0 leads to an increasingly sub-linear relation observed in the present simulations is consistent with that observed experimentally (Boulbitch et al., 2001). 4.3. Instability of the spreading front Experiment results (Boulbitch et al., 2001) and a first-order perturbation analysis (Shenoy and Freund, 2005) showed that the advancing adhesion front of the vesicle is unstable to perturbations in its shape when the radius of the circular spreading area reaches a critical value. During the unstable spreading process, small convex and concave shapes of the spreading front may develop. Experimental results (Boulbitch et al., 2001) also showed that although the spreading front becomes wavy, the wavy amplitude does not grow significantly. The concave parts may become convex and the convex parts may become concave after further growth. To investigate the instability of the advancing adhesion front of the vesicle, a random perturbation was given to the spreading front during spreading. The maximum perturbation amplitude at the spreading front is 2% of the radius of the spreading front. The evolution of the spreading front was then monitored after the perturbation. The present simulations showed that the instability of the vesicle spreading kinetics is influenced by the in-plane membrane shear modulus, the adhesion traction, and the ratio cl =cr0 . It was found that a weaker adhesive traction, or a stiffer membrane, or a smaller ratio of cl =cr0 leads to a more stable spreading front against the small perturbation. For stable cases, the perturbation will decay and the circular spreading front will recover after

ARTICLE IN PRESS 1178

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

the perturbation. If the critical radius for the instability is larger than the radius of the final stable circular spreading area, the spreading may reach a final state, in which the circular contact profile and vesicle shape do not change with time. Examples of such cases are shown in Figs. 2 and 3. To make the advancing adhesion front unstable, the adhesion traction, or ratio cl =cr0 , needs to be increased, or the in-plane membrane shear modulus should be decreased. For example, our calculations show that when cl =cr0 ¼ 8, K  ¼ 0:302, and T max ¼ 1:51, the spreading front is unstable to the perturbation when the radius of the spreading area r X0:416; when cl =cr0 ¼ 8, K  ¼ 0:603, and T max ¼ 2:416, the spreading front is unstable to the perturbation when the radius of the spreading area r X0:332; and when cl =cr0 ¼ 8, K  ¼ 1:81, and T max ¼ 3:02, the spreading front is unstable to the perturbation when the radius of the spreading area r X0:312. Here only the first case above (that is, dc ¼ 0:1, K  ¼ 0:302, cr0 ¼ cl =8 and T max ¼ 1:51) was used to illustrate the details of the unstable spreading kinetics. The other cases share a similar tendency. Fig. 9 shows the snapshots of the evolution of the spreading front before and after the random perturbation. Fig. 9(a) shows the circular contact profile before the perturbation (t* ¼ 0.140, r* ¼ 0.416). Two modulations appear at the spreading front after the perturbation, as shown in Fig. 9(b). More modulations appear with further spreading as shown in Fig. 9(c) and (d). It is seen that the advancing adhesion front of the vesicle is unstable to the perturbation. Small convex and concave modulations of the spreading front were developed during the spreading process. It was observed that although the spreading front becomes modulated,

Fig. 9. Snapshots for the binding area for an unstable case. (a) Before perturbation, the spreading front remains circular at t* ¼ 0.140; (b) after perturbation, two modulations appear at t* ¼ 0.144; (c) more modulations appear at t* ¼ 0.240, and (d) the final state at t* ¼ 0.696.The growth of the modulations appears to be self-limiting, and the positions of the modulations change with spreading time (dc ¼ 0:1, cr0 ¼ cl =8, T max ¼ 1:51, K  ¼ 0:302).

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1179

the modulation amplitude does not grow significantly, indicating that the modulation appears to be self-limiting. In addition, the concave and convex spreading fronts keep moving and changing after further growth. These observations are consistent with the experimental results (Boulbitch et al., 2001). Our simulations also showed that for these unstable cases, a final state can be reached in which the contact area and vesicle shape do not change with time while the spreading front still exhibits small modulations (Fig. 9(d)).

4.4. Discussions Recent studies using mouse embryonic fibroblasts in suspension contact with a matrixcoated surface (Dubin-Thaler et al., 2004) showed that in most cells, there is a mixture of cell spreading in the isotropic and anisotropic modes. These cells could spread in an isotropic fashion at some spreading fronts while at the same time at other spreading fronts anisotropic spreading prevails. The cells could also begin by an isotropic spreading manner and then suddenly transition to an anisotropic manner or begin in an anisotropic manner and then suddenly transition into isotropic spreading. The present simulations show that the circular spreading front may become unstable and change into a modulated spreading front. Locally these concave and convex modulations may disappear or reappear during the spreading process. If the circular spreading can be assumed to be isotropic spreading and the modulated spreading, anisotropic spreading, the similarity between the present simulation results and those experimental observations is evident. This suggests that the transition between isotropic spreading and anisotropic spreading may be related to the instability of the spreading front, since a random perturbation should actually exist at the cell-spreading front. The present work has studied the adhesion and spreading of biomimetic cells. However, in order to simulate real cell adhesion, a number of issues still need to be addressed: First, a more accurate traction separation relation is needed to describe the interaction between receptors and ligands. Currently the exact function for the traction separation law is not clear. Van der Waals forces are expected to be attractive and of the long range. In addition, various proteins present in the lipid bilayer of cell membranes may provide resistance to cell adhesion in the intermediate range. Moreover, the adhesion mediated by receptor–ligand bonds is expected to be strong and short-ranged. Thus, a double minimum potential energy profile is expected (Bell, 1978). Hence, the consideration of more complicated function forms for the traction separation law is an important future work. Second, some important cell components, such as cytoskeleton, are needed to be included in the model. The present model can only be used to predict cell adhesion in the early stage. In this stage, the plasma membrane and substrate surface are adhered together by receptor–ligand molecular binding, and the cell adhesion does not so much depend on the reconstruction of cytoskeleton (Reinhart-King et al., 2005). In the later stages of cell adhesion, however, cytoskeletal apparatuses are assembled around the receptors that directly mediate the adhesion, forming focal adhesion points. Third, this continuum model is valid only under specific conditions. For example, at extremely low ligand and receptor concentrations, the discrete nature of binders dominates cell adhesion. Thus contact should be dominated by discrete and random events, leading to the observed anisotropic behavior (Dubin-Thaler et al., 2004). Finally, the cells can sense and respond to external stimuli. How to incorporate such intelligent behavior in cell adhesion is still a challenging

ARTICLE IN PRESS 1180

P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

issue. Much more work is needed to model and simulate cell spreading, migration, interaction and agglomeration. 5. Conclusions In the present paper, a cohesive model was developed to study the biomimetic cell adhesion to a substrate surface. Finite-element methods were used to solve the governing equations of biomimetic cell mechanics, the receptor diffusion on the vesicle membrane and the contact interaction between the vesicle surface and substrate surface. Parametrical studies of various system parameters such as the membrane shear modulus, adhesive parameters and the receptor density have been performed to understand the mechanics and kinetics of cell spreading. It was found that the advancing adhesion front might be stable or unstable when exposed to small perturbations depending on the system parameters. A higher adhesive traction, a softer membrane or a higher cl =cr0 ratio, favors an unstable spreading front. The present simulations suggest that the growth of the modulation of the spreading front is self-limiting for the unstable cases. References ABAQUS User’s Manual, 2002. Version 6.3, ABAQUS Inc., Pawtucket, RI, USA. Balaban, N.Q., Schwarz, U.S., Riveline, D., Goichberg, P., Tzur, G., Sabanay, I., Mahalu, D., Safran, S., Bershadsky, A., Addadi, L., Geiger, B., 2001. Force and focal adhesion assembly: a close relationship studied using elastic micropatterned substrates. Nat. Cell Biol. 3, 466–772. Bao, G., Bao, X.R., 2005. Shedding light on the dynamics of endocytosis and viral budding. Proc. Natl. Acad. Sci. USA 102 (29), 9997–9998. Bao, G., Suresh, S., 2003. Cell and molecular mechanics of biological materials. Nat. Mater. 2, 715–725. Bell, G.I., 1978. Models for the specific adhesion of cells to cells. Science 200 (4342), 618–627. Bell, G.I., Dembo, M., Bongrand, P., 1984. Cell adhesion. Biophys. J. 45, 1051–1064. Beningo, K.A., Dembo, M., Kaverina, I., Small, J.V., Wang, Y.L., 2001. Nascent focal adhesions are responsible for the generation of strong propulsive forces in migrating fibroblasts. J. Cell Biol. 153, 881–888. Boulbitch, A., Guttenberg, Z., Sackmann, E., 2001. Kinetics of membrane adhesion mediated by ligand–receptor interaction studied with a biomimetic system. Biophys. J. 81 (5), 2743–2751. Dao, M., Lim, C.T., Suresh, S., 2003. Mechanics of the human red blood cell deformed by optical tweezers. J. Mech. Phys. Solids 51, 2259–2280. Doebereiner, H.G., Dubin-Thaler, B., Giannone, G., Xenias, H.S., Sheetz, M.P., 2004. Dynamic phase transition in cell spreading. Phys. Rev. Lett. 93 (10), Art. No. 108105. Dubin-Thaler, B.J., Giannone, G., Doebereiner, H.G., Sheetz, M.P., 2004. Nanometer analysis of cell spreading on matrix-coated surfaces reveals two distinct cell states and STEPs. Biophys. J. 86 (3), 1794–1806. Evans, E.A., 1985. Detailed mechanics of membrane–membrane adhesion and separation. Biophys. J. 48, 175–183. Evans, E.A., Hochmuth, R.M., 1976. Membrane viscoelasticity. Biophys. J. 16 (1), 1–11. Freund, L.B., Lin, Y., 2004. The role of binder mobility in spontaneous adhesive contact and implications for cell adhesion. J. Mech. Phys. Solids 52, 2455–2472. Gallant, N.D., Michael, K.E., Garcia, A.J., 2005. Cell adhesion strengthening: contributions of adhesive area, integrin, binding, and focal adhesion assembly. Mol. Biol. Cell 16, 4329–4340. Gao, H., Shi, W., Freund, L.B., 2005. Mechanics of receptor-mediated endocytosis. Proc. Natl. Acad. Sci. USA 102 (27), 9469–9474. Giannone, G., Dubin-Thaler, B.J., Doebereiner, H.G., Kieffer, N., Bresnick, A.R., Sheetz, M.P., 2004. Periodic lamellipodial contractions correlate with rearward actin waves. Cell 116, 431–443. Hategan, A., Sengupta, K., Kahn, S., Sackmann, E., Discher, D.E., 2004. Topological pattern dynamics in passive adhesion of cell membranes. Biophys. J. 87 (5), 3547–3560.

ARTICLE IN PRESS P. Liu et al. / J. Mech. Phys. Solids 55 (2007) 1166–1181

1181

Irvine, D.J., Hue, K.A., Mayes, A.M., Griffith, L.G., 2002. Simulations of cell–surface integrin binding to nanoscale-clustered adhesion ligands. Biophys. J. 82 (1), 120–132. Kloboucek, A., Behrisch, A., Faix, J., Sackmann, E., 1999. Adhesion-induced receptor segregation and adhesion plague formation: a model membrane study. Biophys. J. 77, 2311–2328. Lim, C.T., Dao, M., Suresh, S., Sow, C.H., Chew, K.T., 2004. Large deformation of living cells using laser traps. Acta Mater. 52, 1837–1845. Liu, P., Cheng, L., Zhang, Y.W., 2001. Measuring interface parameters and toughness—a computational study. Acta Mater. 49 (5), 817–825. N’Dri, N.A., Shyy, W., Tran-Son-Tay, R., 2003. Computational modeling of cell adhesion and movement using a continuum-kinetics approach. Biophys. J. 85 (4), 2273–2286. Park, K., Mao, F.W., Park, H., 1990. Morphological characterization of surface-induced platelet activation. Biomaterials 11, 24–31. Peeters, E.A.G., Oomens, C.W.J., Bouten, C.V.C., Bader, D.L., Baaijens, F.P.T., 2005. Mechanical and failure properties of single attachd cells under compression. J. Biomech. 38, 1685–1693. Reinhart-King, C.A., Dembo, M., Hammer, D.A., 2005. The dynamics and mechanics of endothelial cell spreading. Biophys. J. 89 (1), 676–689. Seifert, U., Lipowsky, R., 1990. Adhesion of vesicles. Phys. Rev. A 42 (8), 4768–4771. Shenoy, V.B., Freund, L.B., 2005. Growth and shape stability of a biological membrane adhesion complex in the diffusion mediated regime. Proc. Natl. Acad. Sci. USA 102 (9), 3213–3218. Tan, J.L., Tien, J., Pirone, D.M., Gray, D.S., Bhadriraju, K., Chen, C.S., 2003. Cells lying on a bed of microneedles: an approach to isolate mechanical force. Proc. Natl. Acad. Sci. USA 100, 1484–1489. Tvergaard, V., Hutchsinson, J.W., 1993. The influence of plasticity on mixed-mode interface toughness. J. Mech. Phys. Solids 41, 1119–1135. Weatherburn, C.E., 1925. On differential invariants in geometry of surfaces, with some applications to mathematical physics. Q. J. Pure Appl. Math. L (3), 230–269. Wolgemuth, C.W., 2005. Lamellipodial contractions during crawling and spreading. Biophys. J. 89 (3), 1643–1649. Zhang, Y.W., Bower, A.F., Xia, L., Shih, C.F., 1999. Three dimensional finite element analysis of the evolution of the voids and thin films by strain and electromigration induced surface diffusion. J. Mech. Phys. Solids 47, 173–199. Zhu, C., 2000. Kinetics and mechanics of cell adhesion. J. Biomech. 33, 23–33. Zhu, C., Bao, G., Wang, N., 2000. Annu. Rev. Biomed. Eng. 2, 189–226.