Computational Fluid Dynamic Simulation of Aggregation of Deformable Cells in a Shear Flow

wafflecanadianΜηχανική

18 Ιουλ 2012 (πριν από 4 χρόνια και 10 μήνες)

559 εμφανίσεις

Prosenjit Bagchi
Department of Mechanical & Aerospace
Engineering,
Rutgers,The State University of New Jersey,
98 Brett Road,
Piscataway,NJ 08854
Paul C.Johnson
Department of Bioengineering,
University of California,San Diego,
La Jolla,CA 92093
Aleksander S.Popel
Department of Biomedical Engineering,School of
Medicine,
Johns Hopkins University,
720 Rutland Avenue,611 Traylor Bldg.,
Baltimore,MD 21205
Computational Fluid Dynamic
Simulation of Aggregation of
Deformable Cells in a Shear Flow
We present computational fluid dynamic (CFD) simulation of aggregation of two deform-
able cells in a shear flow.This work is motivated by an attempt to develop computational
models of aggregation of red blood cells (RBCs).Aggregation of RBCs is a major deter-
minant of blood viscosity in microcirculation under physiological and pathological con-
ditions.Deformability of the RBCs plays a major role in determining their aggregability.
Deformability depends on the viscosity of the cytoplasmic fluid and on the rigidity of the
cell membrane,in a macroscopic sense.This paper presents a computational study of
RBC aggregation that takes into account the rheology of the cells as well as cell-cell
adhesion kinetics.The simulation technique considered here is two dimensional and
based on the front tracking/immersed boundary method for multiple fluids.Results pre-
sented here are on the dynamic events of aggregate formation between two cells,and its
subsequent motion,rolling,deformation,and breakage.We show that the rheological
properties of the cells have significant effects on the dynamics of the aggregate.A stable
aggregate is formed at higher cytoplasmic viscosity and membrane rigidity.We also show
that the bonds formed between the cells change in a cyclic manner as the aggregate rolls
in a shear flow.The cyclic behavior is related to the rolling orientation of the aggregate.
The frequency and amplitude of oscillation in the number of bonds also depend on the
rheological properties.￿DOI:10.1115/1.2112907￿
1 Introduction
Aggregation of cells plays a key role in many important bio-
logical processes.In the blood of humans and many other mam-
mals,red blood cells ￿RBCs￿ aggregate at low shear rates forming
a linear or fractal-like structure called “rouleaux.” During me-
tastasis,cancer cells often flow through the blood stream as ag-
gregates.Platelets form aggregates to prevent blood clotting.The
present work is primarily motivated by an attempt to develop
computational models and simulations of aggregation of RBCs in
a shear flow.RBCs under normal conditions are highly deform-
able.Under pathological conditions,such as sickle cell anemia
and bacterial infections,deformability of RBCs is greatly reduced.
Clinical studies have shown that under such conditions aggrega-
bility of RBCs is also increased ￿1,2￿.It is however not estab-
lished whether there is a connection between reduced deformabil-
ity and increased aggregability of RBCs.We hypothesize that
deformability of cells,and hence their rheological property,plays
a key role in determining their aggregability.To that end,we
develop a computational fluid dynamic simulation to study aggre-
gation of deformable cells.The simulations considered here are
two dimensional￿2D￿,though the methodology can be extended to
3D.This is a simulation of aggregation of cells that takes into
account cell rheology.The methodology is general and can be
applied to any deformable cell,though the specific rheological
values used here correspond to RBCs.
The mechanism of RBC aggregation in blood flow is as fol-
lows.RBCs are randomly drawn close to each other by the flow of
plasma.If the shearing forces by fluid motion are small,the cells
adhere to each other and form aggregates.Currently there are two
theories that describe the mechanism of aggregation:bridging be-
tween cells by cross-linking molecules ￿3￿,and osmotic force gen-
erated by the depletion of molecules in the intercellular space ￿4￿.
RBC aggregation leads to increased blood viscosity,and hence
elevated resistance to blood flow.Aggregation is common in pa-
tients with peripheral vascular disease.Elevated aggregation is
often associated with a higher risk of cardiovascular disease.RBC
aggregation is elevated after myocardial infarction,ischemic brain
infarcts,in diabetes,and during sepsis.Thus,understanding the
mechanics of RBC aggregation may lead to a better understanding
of cardiovascular diseases and the pathology of blood.Studies on
the effect of RBC aggregation in microcirculation are summarized
in recent reviews ￿1,5￿.
Theoretical approaches to the study of aggregation can be
broadly classified into two categories.The first one is the quasi-
steady approach which predicts the equilibrium shape of the
RBCs forming the aggregate through a balance between the adhe-
sive surface energy and the elastic energy stored in the RBC mem-
brane ￿6,7￿.The motion and breakage of an aggregate in a fluid
flow cannot be modeled in this approach.The second approach is
based on the theory of kinetic modeling of colloidal suspension
￿8–10￿.In the kinetic approach deformation of a cell is not con-
sidered.Further,the details of the flow pattern around each de-
formable cell and aggregate are also neglected which can have a
significant effect on the stability and motion of the aggregate.It
has been shown earlier by in vitro ￿11￿,and recently by in vivo
experiments that aggregation has a significant effect on the veloc-
ity profiles in microvessels ￿12–15￿.
Clearly,a theoretical/computational model of cell aggregation
that would take into account deformability of the cells and their
dynamic motion in a flow is lacking.Aggregation is a multiscale
process.At the nanoscale the molecular bridging between adjacent
cells initiates aggregation.At the microscale,deformability of in-
dividual cells dictates the process.At the macroscale,the flow of
plasma imparts a shear force causing the rolling and breakage of
the aggregates.Here we present a computational fluid dynamic
model that combines the three scales of the problem.Specifically
we study the effect of cell deformation,strength of the adhesion
molecule,and the shear rate of the bulk flow on the motion and
Contributed by the Bioengineering Division of ASME for publication in the J
OUR-
NAL OF
B
IOMECHANICAL
E
NGINEERING
.Manuscript received November 12,2004;revi-
sion received July 27,2005.Review conducted by:C.Ross Ethier.Technical Editor:
Frank Yin.Approved August 15,2005.
1070/Vol.127,DECEMBER 2005 Copyright © 2005 by ASME Transactions of the ASME
breakage of an aggregate formed between two cells placed in a
linear shear flow.Simulations presented here are two dimensional,
but the method can be extended to 3D.
2 Simulation Technique
2.1 Front-tracking Method for Deformable Cells.The
simulation technique considered here is the immersed boundary
method developed by Peskin ￿16￿,and later extended by Univerdi
and Tryggvason ￿17￿ and Tryggvason et al.￿18￿ as the front-
tracking method for multiple fluids with different properties.The
cells are modeled as viscous liquid drops surrounded by elastic
membranes.The viscosity of the liquid inside the cell ￿cytoplasm￿
can differ from that of the liquid outside ￿plasma￿.The main idea
of the front-tracking method is to use a single set of equations for
all fluids,inside and outside.The governing equations for the fluid
flow are solved on a fixed Eulerian grid,and the cell-plasma in-
terface is tracked in a Lagrangian manner by a set of moving grid
used to discretize the interface ￿Fig.1￿.This way the interface
tracking is decoupled from the fluid solver.Thus the method is a
mixed Eulerian-Lagrangian method.A body force term F￿x,t￿ is
introduced in the governing equations that accounts for the inter-
face.It is zero everywhere in the flow except at the interface
F￿x,t￿ =
￿
￿S
f￿x
￿
,t￿￿￿x − x
￿
￿dx
￿
￿1￿
where x is a point in the flow domain,x
￿
is a point on the inter-
face,￿S is the entire interface,and ￿is the Delta function.Here
f=f
s
+f
m
is an interfacial force that arises due to the elastic force f
s
generated in the membrane,and the intercellular bridging force f
m
due to cell-cell aggregation.For incompressible and immiscible
fluids,the governing equations are the continuity and the Navier-
Stokes equations:
￿ u = 0,and ￿
￿
￿u
￿t
+ u  ￿u
￿
= − ￿p + ￿ ￿+ F ￿2￿
Here u￿x,t￿ is the fluid velocity anywhere in the flow,￿ is the
density,p pressure,￿=￿￿￿u+￿￿u￿
T
￿ is the viscous stress,and
￿￿x,t￿ is the viscosity in the entire fluid.For any point within a
blood cell ￿=￿
c
and for any point outside ￿=￿
p
,where ￿
c
is the
viscosity of the liquid interior of the cell,and ￿
p
is the viscosity
of the extracellular fluid.As the cells move and deform,￿￿x,t￿
needs to be updated in the flow field.The details about updating
￿￿x,t￿ can be found in ￿17￿ and are not repeated here.The gov-
erning equations are discretized spatially using a finite difference
scheme,and temporally using a two-step time-split ￿or projection￿
scheme.A typical grid used in the simulation is 128￿128 for the
flow domain and 512 for the cell membrane.Details of the time-
split scheme can be found in ￿19￿.
Once the fluid velocity is obtained by the method described
above,the cells are advected in the Lagrangian manner.The ve-
locity of the cell membrane is obtained by interpolating the veloc-
ity of the fluid as
u￿x
￿
￿ =
￿
S
u￿x￿￿￿x − x
￿
￿dx ￿3￿
where S denotes the entire flow domain.The membrane is then
advected by
dx
￿
dt
= u￿x
￿
￿ ￿4￿
The multidimensional ￿function used above is constructed by
multiplying one-dimensional delta functions,such as ￿￿x−x
￿
￿
=￿￿x−x
￿
￿￿￿y−y
￿
￿,in two dimension.The key to the numerical
implementation of the above formulation of the membrane force
is a smooth representation of the ￿function as ￿16,17￿
D￿x￿ =
1
2￿x
￿
1 + cos
￿x
2￿x
￿
for ￿x￿ ￿2￿x ￿5￿
D￿x￿ = 0 for ￿x￿ ￿2￿x ￿6￿
where ￿x is the grid size.The above representation approaches
the actual delta function as the grid size approaches zero.In the
discrete form,the sharp interface is replaced by a rapid variation
over four grid points.The integrals in Eqs.￿1￿ and ￿3￿ can be
written in the discrete form as
F￿x
j
￿ = ￿
i
D￿x
j
− x
i
￿
￿f￿x
i
￿
￿ ￿7￿
u￿x
i
￿
￿ = ￿
j
D￿x
j
− x
i
￿u￿x
j
￿ ￿8￿
where i and j represent grid points along the interface and in the
interior of the domain,respectively.
2.2 Cell Model.The cells are modeled as viscous liquid drops
surrounded by elastic membranes.The viscosity of the cytoplas-
mic fluid can be different from that of the extracellular fluid.This
difference is taken into account in the front tracking method de-
scribed above.The elastic force f
s
in the membrane can be ob-
tained from a strain energy function.In the present model we
assume that the membrane behaves like a neo-Hookean material.
We note that an RBC membrane is strongly resistant to area dila-
tation.The neo-Hookean model employed here does allow area
dilatation.Membrane models that restrict area dilatation and
hence are more accurate for the RBC have been developed
￿20,21￿.The present methodology does allow incorporation of
these models.The neo-Hookean model is chosen because of its
simplicity.An accurate representation of cell membrane is not
central to this study.For the present purpose,a cell model that
takes deformability into account is sufficient.For a 2D neo-
Hookean membrane,the strain energy function is given by ￿see,
e.g.,￿22￿￿
W= E
s
h￿￿
1
2
+￿
2
2
+￿
1
−2
￿
2
−2
￿ ￿9￿
where E
s
is the shear modulus of elasticity of the membrane,h is
the thickness,and ￿
1
and ￿
2
are the principal stretch ratios.The
tensions T
1
and T
2
in the principal directions are then given by
￿22￿
T
1
=
E
s
h
￿
1
￿
2
￿
￿
1
2

1
￿￿
1
￿
2
￿
2
￿
and T
2
=
E
s
h
￿
1
￿
2
￿
￿
2
2

1
￿￿
1
￿
2
￿
2
￿
￿10￿
For the 2D simulations considered here,the membrane is a closed
curve.The 2D cell is then equivalent to an actual 3D cell subject
to a stretching in one direction only.That is,T
1
￿0,T
2
=0,where
“1” indicates the “in plane” direction along the membrane,and
“2” indicates the “out of plane” direction normal to Fig.1.The
deformation ￿
2
in the out of plane direction is not zero.But since
T
2
=0,we can express ￿
2
in terms of ￿
1
.Then,for a 2D cell,we
have
Fig.1 Computational domain and grids
Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1071
T =
E
s
h
￿
3/2
￿￿
3
− 1￿ ￿11￿
where T=T
1
and ￿=￿
1
.For a discretized cell,T is the tension
acting along a line segment connecting two adjacent Lagrangian
grid points on the membrane,and ￿is the stretch ratio ￿unde-
formed length by deformed length￿ of the line segment.At any
Lagrangian grid point on the membrane two line segments meet.
The membrane elastic force f
s
at the Lagrangian point is then the
resultant vector of the tensions in the two adjacent segments:f
s
=T
i
e
i
−T
j
e
j
where i and j denotes two adjacent line segments,and
e
i
and e
j
are the unit tangent vectors along them.
2.3 Aggregation Model.Even though the molecular mecha-
nism of aggregation has not been fully resolved ￿bridging versus
depletion￿,here we use a formalism of bond formation that has
been applied to a wide range of cell-cell adhesion phenomena
￿23,24￿.This is a convenient way to simulate cell-cell adhesion
and it does not imply a specific molecular mechanism.The inter-
membrane force f
m
due to molecular crosslinking is governed by
reaction equations with the reaction rates which are functions of
the local distance between the membranes.The reaction term for a
specific segment of a cell membrane is computed by visiting all
other cells and obtaining the minimum distance between the
neighboring cells.If the minimum distance is less than a threshold
length,the segment pair is given an opportunity to form bonds.
Assuming that both cells contribute an equal number of molecules
to bond formation,the reaction equation for bond density n
b
is
given by
￿n
b
￿t
= 2
￿
k
+
￿
n −
n
b
2
￿
2
− k

n
b
2
￿
￿12￿
where n is the density of the cross-linking molecules on each cell,
and k
+
and k

are the forward and reverse reaction rate coeffi-
cients,respectively.Note that the above equation is slightly dif-
ferent from that for a receptor-ligand interaction,but consistent
with the models of Zhu ￿24￿.The reaction rates are computed as
k
+
= k
+
0
exp
￿

k
ts
￿l − l
o
￿
2
2K
B
T
￿
,if ￿x￿ = l ￿l
t
￿13￿
k

= k

0
exp
￿

￿k
b
− k
ts
￿￿l − l
o
￿
2
2K
B
T
￿
￿14￿
where ￿x￿ is the distance between two facing segments of the cells,
l
t
is the threshold distance below which the bond formation is
initiated,l =￿x￿,and l
o
are the stretched and unstretched bond
length,k
b
is the spring-constant ￿force per stretched length￿,k
ts
is
the transition spring constant,K
B
is the Boltzmann constant,T is
the absolute temperature,and k
+
0
,and k

0
are the ratecoefficients in
equilibrium.The bonds behave like stretched springs,and the
force per bond is given by
f
b
= k
b
￿l − l
0
￿.￿15￿
The aggregation force f
m
per unit length of the cell membrane is
then obtained as
f
m
= f
b
n
b
x
￿x￿
.￿16￿
Note that in this formulation the bond density,rate coefficients,
bond length,and the intermolecular force are not constant and
vary along the contact area of the cells.The values of the param-
eters,such as the spring constants and the reaction rates used in
the present study are given in Table 1.
In the present computations,the threshold distance ￿l
t
￿ is taken
to be about 1.2 times the unstretched length ￿l
0
￿;thus the bonds do
not stretch significantly.We further emphasize that the detailed
molecular mechanism of aggregation is not fully understood,and
the “stretched spring” model assumed here may not be the best
model for RBC aggregation.The theory of polymer bridging is
poorly understood ￿25￿.There is,however,evidence in support of
the “stretched spring” model used here.In experiments,RBC ag-
gregation is mediated by Dextran molecules ￿3,12–15￿.Many ex-
perimental studies,such as those by Chien and Jan ￿3￿ and Brooks
￿26￿,concluded that the Dextran polymers behave like “flexible
linear chains with random coiling.” The more recent works of
Swenson et al.￿25￿ and Barshtein,Tamir,and Yedgar ￿27￿ also
concluded that the polymer bridges are “partly stretched” chains
between the cells.The theoretical work of Skalak and Zhu ￿28￿
also supports the “stretched spring” model for the RBC aggrega-
tion,similar to the one described here.
In the present study we consider two-dimensional simulations.
The flow domain considered here is a Couette flow in a channel
driven by two walls at the top and bottom that move in directions
opposite to each other.In the absence of any cell,the undisturbed
flow is thus a linear shear flow with shear rate denoted by ￿
˙
.The
distance between the top and bottom wall is about 25 ￿m,typical
of a blood vessel dimension in microcirculation.The length of the
channel is also taken as 25 ￿m,and the flow at the left and right
boundaries of the channel is assumed to be periodic.The govern-
ing equations are made dimensionless using a characteristic length
a ￿such as cell radius￿,the applied shear rate ￿
˙
,and the plasma
viscosity ￿
p
.The dimensionless time is then t￿˙
where t is the
physical time.The characteristic scale for velocity is a￿˙
.This
results in the following major dimensionless parameters:
Reynolds number,Re =
￿￿˙
a
2
￿
p
Membrane rigidity,E
*
=
￿
p
￿˙
a
E
s
h
Table 1 Physical and dimensionless parameters used in the simulations.Here a
0
is a characteristic length scale which is often
the radius of an undeformed cell.For an elliptic cell studied here,the semimajor axis of the undeformed cell is taken as the cell
radius.The number of molecules per cell is in the range 50–300.The number of bonds formed in equilibrium is in the range
30–100.The number density n is obtained assuming the cell surface area 100 ￿m
2
.
Physical parameters Dimensionless parameters
Cell radius a 4 ￿m a/a
0
1–0.667
Vessel diameter D 25 ￿m D/a 2￿
￿
˙
25-225 s
−1
Re 0.0003–0.003
E
s
1.2-30￿10
−3
dyn/cm
E
*
0.02–0.5
k
b
10
−6
–10
−8
N/m
K
b
0.125–3.0
￿
p
=1.2 cP ￿
c
=1.2–6 cP ￿
1–5
k
ts
10
−10
–10
−12
N/m
C
1
*
0.3–0.003
k
+
0
10
−7
–10
−9
s
−1
cm
−2
k
+
0
*
0.1–0.001
k

0
10
−11
–10
−13
s
−1
cm
−2
k

0
*
10
−3
–10
−5
1072/Vol.127,DECEMBER 2005 Transactions of the ASME
Viscosity ratio,￿ =
￿
c
￿
p
Spring constant,K
b
=
k
b
￿
p
￿
˙
a
Values of the dimensionless parameters used in this study are
given in Table 1.The “membrane rigidity” E
*
is the inverse of the
shear modulus E
s
.For a normal RBC,E
s
=6￿10
−3
dyn/cm,and
￿=5.However,under pathological conditions,such as sickle cell
disease and malaria,significant deviation from these values may
occur.Further,the parameters governing the kinetics of the RBC
aggregation under pathological conditions are not clearly known.
Clinical studies do show a tendency of high aggregation under
such conditions.Thus considering a range of the parameters as
used in our simulations ￿Table 1￿ is of biological interest.
3 Results and Discussion
3.1 Simulation of Single Cell.We first discuss some results
on the simulation of deformation of a single cell subject to a linear
shear flow.The geometry of the flow domain is the same as that
described in Fig.1.A single,isolated cell is first placed at the
center of the channel.The initial cell shape is assumed to be
elliptic with an aspect ratio of 0.375.As the flow starts,the cell
deforms,and the cell surface and the cytoplasmic fluid rotate in a
tank-treading fashion.The trajectory of a fixed point on the cell
surface and streamlines obtained from the velocity field are shown
in Figs.2￿a￿ and 2￿b￿,respectively.In Fig.2￿c￿,we plot the steady
shape of the cell at various shear rates ￿￿
˙
=25,75,and 225 s
−1
￿.
The cell elongates and aligns itself with the direction of the flow
as the shear rate is increased.The effect of varying shear modulus
is described in Fig.2￿d￿ which also shows that the cell becomes
more elongated with decreasing rigidity.The elongation of the cell
is accompanied by a slight increase in the surface area which is
found to be at most 2%.In Fig.2￿e￿,the result for ￿=5 is shown.
A rotation of the cell like a rigid body at higher cytoplasmic
viscosity is observed.The direction of rotation is consistent with
the direction of the applied shear.In the figure only one instanta-
neous orientation is shown after the cell starts rotating.There is
almost no visible change in the cell shape.Though the cell mem-
brane undergoes tank-treading motion,its velocity is much less
than the earlier cases of ￿=1 shown in Figs.2￿b￿–2￿d￿.These
results agree qualitatively with those of previous researchers,e.g.,
Pozrikidis ￿29￿ and Ramanujan and Pozrikidis ￿30￿.These studies
considered 3D cell deformation at zero inertia,whereas the
present study is 2D and inertia is not zero,though small.
It should be noted that the smoothing of the cell membrane
using the discrete delta function in the immersed boundary
method introduces a numerical viscosity.While such a numerical
viscosity is extremely bad for flow computations at high Reynolds
numbers,it is not detrimental for the low Reynolds number cases
considered here Re￿0.001￿.The projection method that is used
here to solve the Navier-Stokes equations satisfies the incompress-
ibility condition up to the precision of the computer ￿10
−14
￿.Thus
mass conservation is exactly satisfied on a local basis.We also
keep track of the cell volume during the simulations.The change
in the cell volume is less than ±0.1% for a grid resolution of
128￿128.The cell does not show any net increase/decrease in
volume over a long time.Thus the small error is probably due to
the numerical error in computing the cell volume,rather than an
error in the front-tracking method.Since the simulations presented
here are 2D,there is no direct way of validating the results on cell
deformation.However,we have been able to extend the method to
3D,and thus been able to compare the deformation of a 3D
spherical capsule obtained from our immersed boundary simula-
tion with the boundary integral calculation of Pozrikidis ￿29￿.The
comparison is shown in Fig.2￿f￿ for a low Reynolds number of
0.001.Note that Eggleton and Popel ￿31￿ had also used a similar
immersed boundary method for large deformation of capsules.A
detailed comparison of the immersed boundary results,boundary
integral results,and asymptotic results in the limit of small defor-
mation,is given in their paper.
3.2 Simulation of Aggregation.Simulations are performed
for aggregation of two elliptic cells in Couette flow at shear rates
typical of physiological conditions.Simulations with biconcave
shape of an RBC are also performed and described later.Table 2
lists the range of shear rate,spring constant,membrane rigidity,
and viscosity ratio used in the simulations.The simulation strat-
egy is as follows.First we let the cells aggregate under a no-flow
condition as shown in Fig.3￿a￿.The static,no-flow condition is
chosen as the starting condition since aggregation is believed to be
the highest under such a condition.To initiate aggregation,we
first place the cells in close proximity such that the minimum
distance between the adjacent cells is below the threshold distance
l
t
.The number of bonds and the contact area between two cells
increase until they reach steady values.The cells reach final
steady shapes,and an aggregate is formed.Once the equilibrium
is reached under no-flow condition,the shear flow is imposed.
Under the action of the shear flow,the aggregates roll,deform,
and depending on the spring constant,shear rate and other param-
eters,they either remain intact or break apart.In the following
section we describe such dynamic events as revealed by our
simulations.
Figures 3￿b￿–3￿e￿ show the instantaneous orientations of the
rolling aggregate at various ￿,and K
b
.Figure 3￿b￿ corresponds to
￿
˙
=75 s
−1
,K
b
=2,￿=1.The aggregate rolls almost like a rigid
body.Consequently,the tank-treading motion of an individual cell
in the aggregate is not present any more.The trace of a point on
the membrane surface is shown which clearly indicates no tank
treading.The individual cells undergo some degree of deforma-
tion by the fluid.The contact area and the total number of bonds
remain nearly steady.Consequently,the aggregate does not break
and remains stable.Interestingly the contact area does not remain
flat;it changes shape during each turn of the aggregate.In Fig.4,
the instantaneous flow field around the rolling aggregate is pre-
sented.We also note that the rolling of the aggregate induces
fluctuations in the surrounding flow as evident from the color
contours in the figure.Such disturbed flow field can have an effect
on the cells,such as leukocytes and platelets adhering to the wall,
and can either dislodge them or stabilize the cell-to-wall adhesion.
The folded shape of the cell interface in these figures is some-
what peculiar.There is little experimental work reporting the
shape of a rotating aggregate.Most previous works have reported
aggregate shapes under static conditions ￿3,28,32￿.However,
Skalak et al.￿6￿ reported the shape ￿Fig.8 in their paper￿ of an
aggregate obtained via electron microscopy which shows the sig-
moidal shape of the interface and thus is very similar to the ones
found here.It is thus possible to see such shapes in 3D.The
sequence of four snapshots given in Fig.3￿b￿ suggests that the
interface is rotating like a flexible but nonextensible element ￿or,
thread￿.This is expected since the membrane can easily deform
but not extend.Simulations of deformable and nonextensible ele-
ments in shear flow have shown similar nature of folding ￿see,
e.g.,￿32￿￿.Our simulations also show strong curvatures near the
edge of the interface,particularly at high K
b
.Skalak and Chien
￿33￿ showed that as the attractive force increases,the curvature
near the edge also increases.Their work considered a two-cell
rouleaux similar to the present work,and the semicircular shape
of individual cell obtained in their work is similar to our result
corresponding to the no-flow condition ￿Fig.3￿a￿￿.It is shown
later ￿Sec.3.3￿ that inclusion of bending resistance in the model
would prevent the folding of the interface and the membrane.
Results at a lower spring constant K
b
=0.875 are shown in Fig.
3￿c￿.The contact area is now less than before as many bonds are
broken by the fluid force.However,the bond breaking rate is
slow,and the aggregate does not break instantly.We have per-
formed this simulation up to t￿˙
=110￿￿1.5 s￿ at which point
Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1073
about 35% of bonds still exist.At K
b
=0.175,shown in Fig.3￿d￿,
the aggregate deforms significantly from its static shape.The as-
sembly very quickly breaks,and the cells are completely detached
from each other.The effect of increasing ￿ is shown in Fig.3￿e￿.
Increased cytoplasmic viscosity reduces the deformability of the
cells as seen before in Fig.2￿e￿ for a single cell.The aggregate
shape is now nearly circular,and individual cells almost retain
their semicircular shape that they obtain under no-flow conditions.
Increased viscosity makes the bonds very stable,and hence no
disaggregation is observed.
Fig.2 Simulation of single cell deformation.„a… An elliptic cell in a shear flow:The trajectory of a point on the cell
surface is indicated to show tank treading.„b… Streamlines constructed fromthe velocity field in and around the cell.
The cytoplasmic fluid rotates with the cell surface.Also shown are the effects of varying shear rate „c…,membrane
rigidity „d…,and cell-to-plasma viscosity ratio „e….Initial cell shape is shown by the solid line.In „c… Œ ￿˙
=25 s
−1
,¯￿˙
=75 s
−1
,--- ￿
˙
=225 s
−1
.In „d… Œ E
*
=0.02,¯E
*
=0.1,--- E
*
=0.5.In „e… ￿=5.In „f…,the Taylor deformation parameter for
an initially spherical capsule deforming in a shear flow is shown.The solid line corresponds to the present front-
tracking method extended to 3D,and the symbols represent the boundary integral results †29‡.
1074/Vol.127,DECEMBER 2005 Transactions of the ASME
In the absence of bending resistance,as in the above figures,the
cell shape is not unique at zero transmural pressure,i.e.,the pres-
sure difference across the RBC membrane,￿inside pressure - out-
side pressure￿ under hydrostatic conditions.A detailed discussion
on the effect of transmural pressure on the resting shape of a cell
was given by Pozrikidis ￿34,35￿.Thus in the simulations presented
here,an initially elliptic cell that undergoes deformation in a shear
flow may not return to the same initial shape if the flow is
stopped.Note that in our simulations,the transmural pressure is
not a parameter that we specify.The pressure inside and outside
the cell evolves naturally.Furthermore,the shear flow considered
in our case is not a pressure driven flow.Thus in absence of a cell,
the pressure is constant everywhere in the flow field and can be
taken to be zero.Under the static no-flow condition,the transmu-
ral pressure is zero.As the flow starts,the cell deforms,and an
internal pressure higher than the external pressure is established.
Thus the transmural pressure is positive for a deformed cell in a
flow.The cell shape ￿and the transmural pressure￿ under the influ-
ence of the external shear is then unique for a given shear rate and
membrane properties.When two cells form an aggregate,the
transmural pressure is again found to be positive in both static and
rolling conditions.If the aggregate breaks and the cells detach
completely from each other ￿with the shear flow still acting￿ the
cells will return to the deformed shape that they assumed when
placed in a shear flow without aggregation.If the aggregate does
not break at all,but the shear flow is switched off,the shape of the
aggregate returns to the initial steady shape that was obtained
under static conditions.A negative transmural pressure is not ob-
served in our simulations.The conclusions remain the same even
when bending resistance is included.
In Fig.5,we study the effect of varying K
b
,￿,￿
˙
and E
*
on the
evolution of total number of bonds between the cells.The increase
in the number of bonds under no-flow condition is shown in Fig.
5￿a￿.The equilibrium is reached at around t￿˙
￿95.At this point
the shear flow is started,and the number of bonds starts changing
as the aggregate begins to roll.The number of bonds decreases
very rapidly for K
b
=0.175.In this case 95% of the bonds are
broken in 0.25 s ￿t￿˙
￿19￿.At K
b
=0.875 and 1,the decrease is
quite slow.Finally,at K
b
=2 and 3 the aggregate is stable,and the
number does not decrease at all.The effect of viscosity ratio is
shown in Fig.5￿b￿,where results for ￿=1 and 5 are plotted.
Higher cytoplasmic viscosity increases the stability of the aggre-
Fig.3 Simulation of aggregate formation between two cells.
„a… Evolution of cell shape under the no-flow condition.Three
time instants are t￿˙
=0,50,and 97.„b… ￿˙
=75 s
−1
,K
b
=2,￿=1.
The four instants are approximately at equal time interval dur-
ing a rotation of the aggregate.„c… Effect of reduced bond
strength at K
b
=0.875,￿=1,￿
˙
=75 s
−1
.„d… K
b
=0.175,￿=1,￿
˙
=75 s
−1
.„e… Effect of cell viscosity:K
b
=0.875,￿=5,E
*
=0.096,￿˙
=75 s
−1
.
Fig.4 Fluid velocity vectors and color contours showing the
fluid velocity component in the x direction during aggregate
rolling
Table 2 Simulation runs for two cell aggregate
Simulation runs ￿˙
￿s
−1
￿ K
b
￿ E
*
1 75 0.175 1 0.096
2 75 0.875 1 0.096
3 75 1.0 1 0.096
4 75 2.0 1 0.096
5 75 3.0 1 0.096
6 75 0.175 5 0.096
7 75 0.875 5 0.096
8 75 0.875 1 0.02
9 75 0.875 1 0.5
10 25 2.625 1 0.096
11 225 0.292 1 0.288
Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1075
Fig.5 Evolution of the number of bonds with time.„A… Effect of spring constant
K
b
.Curves „a… to „e… represent K
b
=0.175,0.875,1.0,2.0,and 3.0.„B… Effect of
viscosity ratio ￿.The four curves are for „a… K
b
=0.175,￿=1,„b… K
b
=0.875,￿=1,
„c… K
b
=0.175,￿=5,„d… K
b
=0.875,￿=5.„C… Effect of varying shear rate „a…,„b…,and
„c… in the plot… and varying membrane rigidity „b…,„d…,and „e….„a… ￿
˙
=225 s
−1
,„b…
￿˙
=75 s
−1
,and „c… ￿˙
=25 s
−1
.„b… E
*
=0.096,„d… E
*
=0.5,and „e… E
*
=0.02.
1076/Vol.127,DECEMBER 2005 Transactions of the ASME
gate.For K
b
=0.175 at which a rapid aggregate breakage was ob-
served earlier at ￿=1,now results in a much slower bond break-
ing at ￿=5.
The effects of shear rate and membrane rigidity are considered
in Fig.5￿c￿.Curves ￿a￿,￿b￿,and ￿c￿ are the results for ￿
˙
=225,75,
and 25 s
−1
,respectively.At the highest shear rate the aggregate is
quickly broken.At the intermediate shear rate,bond breaking
takes place over a longer time.At the lowest rate,the aggregate is
completely stable.The effect of membrane rigidity is elucidated
by the curves ￿b￿,￿d￿,and ￿e￿ in the same figure.The curves
represent E
*
=0.096,0.5,and 0.02 with 0.096 corresponding to
the normal RBCs.The bond breaking decreases for both cases
when E
*
increases or decreases from this value.Thus it appears
that E
*
=0.096 is near an optimum value at which the aggregate
breaks at the highest rate.We also note that the aggregate breaking
does not happen instantly.For different combinations of param-
eters,the breakage can take place over a few seconds.
Figures 5￿a￿–5￿c￿ also suggest that the number of bonds oscil-
lates as the aggregate rolls in a shear flow.The shape and orien-
tation of the aggregate at five different time instants are shown in
Fig.6 together with the number of bonds.The oscillations are
related to the rolling orientation of the aggregate.A rapid bond
breakage is initiated when the contact area is along the flow di-
rection.On the contrary,the number of bonds increases when the
contact area is normal to the flow direction.In the first case,the
bonds break primarily by the shearing force of the fluid acting
parallel to the contact area.In the latter orientation,the bond
breakage is by pulling the cells in opposite directions by the fluid
forces normal to the contact area.The results are in agreement
with the experimental observation by Chien and Jan ￿3￿ and Chien
et al.￿36￿ that the shearing force is more effective to break an
aggregate than the normal force.Figures 5￿a￿–5￿c￿ also suggest
that the frequency and amplitude of the oscillations in bond num-
bers are not the same for all cases.For easily deformable cells the
amplitude of oscillation is higher and the frequency is lower.With
decreasing deformability,￿e.g.,at higher ￿ and K
b
￿ the amplitude
of oscillation is lower but the frequency is higher.
It should be noted that the total number of bonds depends both
on the change in bond density and change in contact area.If the
aggregate is stable,the contact area is found to remain nearly the
same,but the bond density varies with time giving rise to the
observed oscillatory behavior of the total number of bonds.If,on
the other hand,the aggregate is breaking,the contact area and the
bond density both are decreasing.At any given time,the bond
density is found to be nearly constant over the contact area,except
near the two edges where it smoothly goes to zero.Thus it seems
that instead of using Eq.￿12￿,one can assume a constant bond
density and contact area to get the similar result with an appropri-
ate choice of the spring constant.Such an approach would be
correct for quasi-steady simulations.However,the bond density
and contact area evolve in time,and they are not known a priori.
Thus Eq.￿12￿ is a convenient way to model the bond kinetics.
Figures 5￿a￿–5￿c￿ show that with different choices of the spring
constant one can either get a stable or breaking aggregate for a
constant shear rate.Alternatively one can fix the spring constant to
a specific value and get either a breaking or a stable aggregate
under varying shear rates.
The instantaneous orientation of an aggregate defined by an
angle ￿that a line joining the center of mass of each cell forms
with the direction of the shear flow ￿x axis￿ is shown in Fig.7￿a￿.
The location of the center of mass of each cell can be obtained
from the simulation data.As the aggregate rolls,￿varies periodi-
cally from 0 to 2￿.When ￿=0 and ￿,the contact area is normal
to the flow direction,and when ￿=￿/2 and 3￿/2,it is parallel.
The results show that as the spring constant or cell viscosity in-
creases the period of oscillation decreases.This can be understood
from the theory of rotation of an ellipsoidal particle in a shear
flow ￿37￿.In the limit of a 2D circular particle,the angular orien-
tation changes linearly with time,and the period of rotation is
minimum.A departure from the circular shape means nonmono-
tonic behavior in angular orientation and angular velocity,and an
increase in period of rotation.At ￿=5 at which the aggregate
shape is nearly circular,the result is close to Jeffery’s solution.
However,at reduced ￿ or K
b
,￿is nonmonotonic,and the period
of rotation is higher.The angular velocity ￿
˙
shown in Fig.7￿b￿ is
also periodic with the period half that of ￿,as expected from the
Jeffery’ solution.The angular velocity is low when the contact
Fig.6 Explanation of oscillations in the total number of bonds as the aggregate rolls in
shear flow.The shape and orientation of the aggregate are shown at a few time instants.
Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1077
area is normal to the flow direction,and high when it is parallel.
The simulations presented here are 2D.Quantitative differences
may arise if 3D simulations are performed.However,the major
conclusion of the paper,that aggregability increases with reduced
deformability,is expected to hold even in 3D.
3.3 Effect of Viscoelasticity and Bending Resistance.Vis-
coelastic nature of the membrane of a red blood cell has been
reported in the literature ￿e.g.,￿38￿￿ and measured to be
about 0.001 dyn s/cm.The results described in the previous sec-
tion do not include the effect of membrane viscoelasticity.In order
to include this effect,Eq.￿11￿ is modified by adding a term
￿￿
m
/￿￿￿d￿/dt￿,where ￿
m
is the membrane viscoelasticity.Several
dimensionless parameters then can be formed,such as ￿
m
/￿
p
a
and ￿
m
￿˙
/E
s
,where a is the characteristic length of the cell.The
parameter ￿
m
/￿
p
a is about 100 for a normal RBC implying a
strong effect of viscoelasticity.The second parameter,also called
Deborah number ￿39,40￿,is in the range 0.01 and 1 for shear rates
between 1 s
−1
and 100 s
−1
.The Deborah number compares the
membrane viscous and elastic forces,and could be more relevant
for the present problem.Since RBC aggregation is dominant at a
low shear rate,a corresponding low value of the Deborah number
indicates that the membrane viscous force may be small compared
to the elastic force justifying the exclusion of the viscoelastic
effect in the results presented so far.Nevertheless,we have in-
cluded the viscoelasticity in our model and repeated several simu-
lations for elliptic cells at ￿
m
/￿
p
a=20.Sample results are shown
in Figs.8 and 11 for one representative case.The results suggest
that the rotation rate of the rolling aggregate and the bond break-
ing rate are reduced when viscoelasticity is included.However,
the shape of the cells,the oscillating nature of the bond breaking,
and the overall decrease in the number of bonds remain nearly the
same compared to the case without viscoelasticity.The decrease
in rotation rate and the bond breaking rate is probably due to the
memory effect of the membrane viscosity.Adetailed investigation
of this effect is beyond the scope of this paper.It should be noted
that there has been very little theoretical work done on the defor-
mation of viscoelastic membranes in shear flow.
The RBC membrane also has a bending modulus E
B
of about
1.8￿10
−12
dyn cm.The bending resistance prevents formation of
high curvature during deformation.The RBC also has a biconcave
shape under resting condition.The shape of the cell is expected to
have an effect on aggregation.The bending resistance is also re-
sponsible for the biconcave shape of the cell.To include the bend-
ing resistance in our simulation,we follow the approach by
Pozrikidis ￿41￿.Equation ￿11￿ is modified as T=￿t+qn,where ￿
is the membrane tension obtained by differentiating the strain-
energy function with respect to the extension ratio,q=dm/dl is
the transverse shear tension,m=E
B
￿￿−￿
r
￿ is the bending moment,
￿is the local curvature,￿
r
is the reference curvature in resting
configuration,l is the arc length along the membrane,and t and n
are the unit tangent and normal vectors.The relevant dimension-
less parameter is E
B
/￿a
2
E
s
￿;for a normal RBC and a￿4 ￿m,this
is about 10
−3
.First we have validated our simulations by repeating
some results given in Pozrikidis ￿41￿ which are shown in Fig.9.
This figure shows the stages of evolution of an initially elliptic
cell and a biconcave cell toward the reference state of a circle.The
2D immersed boundary simulations show qualitatively similar be-
havior of the 3D boundary integral simulations by Pozrikidis ￿41￿.
We have repeated several simulations of cell aggregation pre-
sented in the last sections with the bending moment included.Two
representative cases are shown in Figs.10 and 11.Both the static
Fig.7 „a… Angular orientation and „b… velocity of an aggregate
at ￿˙
=75 s
−1
.The dash line is Jeffery’s †37‡ solution.-￿- K
b
=0.875,￿=5;-￿- K
b
=0.875,￿=1;-￿- K
b
=3,￿=1.
Fig.8 Effect of membrane viscoelasticity on the instanta-
neous shape of the rolling aggregate.„a… Viscoelasticity in-
cluded.Four different time instants are shown:t￿˙
=42,48,53,
and 58.„b… Viscoelasticity not included.Shapes at t￿˙
=42 and
48 are shown.
Fig.9 Relaxation of an initially elliptic „left… cell and a biconcave „right… cell toward their resting circular
shape under the action of bending resistance at E
B
/„a
2
E
s
…=0.003,￿=1
1078/Vol.127,DECEMBER 2005 Transactions of the ASME
and rolling cases for elliptic and discoid shapes are considered.
For the elliptic cells,the results under the no-flow condition show
that the inclusion of bending moment removes any high curvature
near the edge of the contact area.However the contact area re-
mains nearly the same in the two cases.The shape of the rolling
aggregate is similar to the one observed without bending moment,
except that the folding of the contact line is reduced.The variation
in the number of bonds ￿Fig.11￿ also shows similar trends in
terms of bond breaking and oscillations for the two cases.Similar
results were found for all cases that were repeated.Thus we con-
clude that for an aggregate rolling in a shear flow,bending resis-
tance may not have a significant effect.Of course,under static
condition,as seen in Fig.10,bending resistance has an effect on
cell shape.We have also performed one simulation using bicon-
cave shape ￿Fig.10￿.The overall trends of aggregation,aggregate
rolling,and bond breakage are similar to those for elliptic cells.
As expected,the bond breaking rate for the biconcave cells is
slower due to higher contact area.
4 Conclusion
In this paper we present 2D numerical simulations on the ag-
gregation of deformable cells in a Couette flow.The simulation
technique couples the bulk hydrodynamics with the rheology of
the cells and the kinetics of intercellular bond formation.The
focus of the paper is to describe ￿1￿ the effect of cell rheological
Fig.10 Effect of bending resistance and biconcave shape on the shape evolution of an aggre-
gate.Both elliptic and biconcave cells are considered.Figures at left and right show the formation
of the aggregate under static and rolling conditions,respectively.In „a…,the solid line is for bend-
ing moment included,and the dash line is without bending moment.
Fig.11 Effect of viscoelasticity and bending resistance on the
number of bonds.Bold line is for elliptic cells without bending
resistance and viscoelasticity.…,viscoelastic membrane;—,
elliptic cell with bending resistance;–––––,biconcave discoid
cell with bending resistance.
Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1079
properties on the stability of an aggregate and ￿2￿ the dynamics of
the rolling motion of the aggregate in the shear flow.
The major conclusion from this work is that the deformability
of the cells as determined by cell-to-plasma viscosity ratio and the
shear modulus of the cell membrane,has a significant effect on
the stability and motion of the aggregate.Aggregates made of
deformable cells are easily breakable by a shear flow,while those
made by less deformable cells are not.These observations are
consistent with many clinical studies which show that increased
aggregation of red blood cells is often accompanied by reduced
cell deformability.Such is the case in sickle cell disease,malaria,
and hypertension.The present study,however,is the first theoret-
ical study to establish a direct link between deformability and
aggregability.
We also observe that in a shear flow,the number of bonds
oscillates with time.The oscillations are related to the rolling
orientation of the aggregate.The bond breakage is maximum
when the contact area is parallel to the flow direction,and mini-
mum when the contact area is normal.This result implies that the
shearing force is more efficient to break the bonds than the normal
pulling force.The frequency and amplitude of the oscillations are
dictated by the rheology of the cells and the bond strength.A
deformable aggregate has higher amplitude and longer period of
rotation than a stable aggregate.This is evident when angular
orientation and velocity are taken into consideration which are
found to match qualitatively with the analytical prediction of Jef-
fery ￿37￿.
Acknowledgment
This work was supported by a grant from the National Heart,
Lung,and Blood Institute HL 52684.We also thank the referees
for their comments which have significantly improved the quality
of the paper.The computations were performed on IBM p690 at
NCSA,UIUC and HP V2500 at CACR,Caltech.
References
￿1￿ Rampling,M.W.,Meiselman,H.J.,Neu,B.,Baskurt,O.K.,2004,“Influence
of Cell Specific Factors on Red Blood Cell Aggregation,” Biorheology 41,pp
91–112.
￿2￿ Ami,B.,Barshtein,G.,Zeltser,D.,Goldberg,Y.,Shapira,I.,Roth,A.,Keren,
G.,Miller,H.,Yedgar,S.,2001,“Parameters of Red Blood Cell Aggregation
as Correlates of the Inflammatory state,” Am.J.Physiol.,280,pp.H1982–
H1988.
￿3￿ Chien,S.and Jan,K.-M.,1973,“Ultrastructural Basis of the Mechanism of
Rouleaux Formation,” Microvasc.Res.5,pp.155–166.
￿4￿ Neu B.,and Meiselman H.J.2002,“Depletion-Mediated Red Blood Cell
Aggregation in Polymer Solutions,” Biophys.J.,83,pp.2482–2490.
￿5￿ Popel,A.S.and Johnson,P.C.,2005,“Microcirculation and Hemorheology,”
Annu.Rev.Fluid Mech.37,pp.43–69.
￿6￿ Skalak,R.,Zarda,P.R.,Jan,K.-M.,and Chien,S.,1981,“Mechanics of
Rouleau Formation,” Biophys.J.,35,pp 771–781.
￿7￿ Evans,E.A.,1985,“Detailed Mechanics of Membrane-Membrane Adhesion
and Separation.” Biophys.J.48,pp 175–183.
￿8￿ Murata,T.,and Secomb,T.W.,1988,“Effects of Shear Rate on Rouleau
Formation in Simple Shear Flow,” Biorheology,25,pp.113–122.
￿9￿ Chen,J.,and Huang,Z.,1996,“Analytical Model for Effects of Shear Rate on
Rouleau Size and Blood Viscosity,” Biophys.Chem.96,pp.273–279.
￿10￿ Lim,B.,Bascom,P.A.J.,and Cobbold,R.S.C.,1997,“Simulation of Red
Blood Cell Aggregation in Shear Flow,” Biorheology,34,pp.423–441.
￿11￿ Goldsmith,H.L.,and Karino,T.,1980,“Physical and Mathematical Models of
Blood Flow:Experimental Studies,￿ Erythrocyte Mechanics and Blood Flow,
Cokelet,G.R.,Meiselman,H.J.,and Brooks,D.E.,eds.,John Willey & Sons.
￿12￿ Bishop,J.J.,Nance,P.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001a,
“Effect of Erythrocyte Aggregation on Velocity Profiles in Venules,” Am.J.
Physiol.280,pp.H222–H236.
￿13￿ Bishop,J.J.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001b,“Effects
of Erythrocyte Aggregation and Venous Network Geometry on Red Blood Cell
Axial Migration,” Am.J.Physiol.281,pp.H939–H950.
￿14￿ Bishop,J.J.,Nance,P.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001c,
“Erythrocyte Margination and Sedimentation in Skeletal Muscle Venules,”
Am.J.Physiol.281,pp.H951–H958.
￿15￿ Bishop,J.J.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2002,“Effect of
Aggregation and Shear Rate on the Dispersion of Red Blood Cells Flowing in
Venules,” Am.J.Physiol.283,pp.H1985–H1996.
￿16￿ Peskin,C.S.1977,“Numerical Analysis of Blood Flow in the Heart,” J.
Comput.Phys.,25,pp.220–233.
￿17￿ Unverdi,S.O.and Tryggvason,G.1992,“A Front-Tracking Method for Vis-
cous Incompressible Multi-Fluid Flows,” J.Comput.Phys.100,pp.25–37.
￿18￿ Tryggvason,G.,Bunner,B.,Esmaeeli,A.,Juric,Al-Rawahi,N.,Tauber,W.,
Han,J.,Nas,S.,and Jan.,Y.-J.,2001,“A Front Tracking Method for the
Computations of Multiphase Flow,” J.Comput.Phys.169,pp 708–759.
￿19￿ Bagchi,P.,and Balachandar,S.,2002,“Steady Planar Straining Flow Past a
Rigid Sphere at Moderate Reynolds Number,” J.Fluid Mech.,466,pp.365–
407.
￿20￿ Skalak,R.,Tozeren,A.,Zarda,P.R.,and Chien,S.,1973,“Strain Energy
Function of Red Blood Cell Membranes,” Biophys.J.,13,pp.245–264.
￿21￿ Evans,E.A.,and Skalak,R.1980,Mechanics and Thermodynamics of
Biomembranes,CRC Press,Boca Raton,FL.
￿22￿ Barthes-Biesel,D.,Diaz,A.,and Dhenin,E.,2002,“Effect of Constitutive
Laws for Two-dimensional Membranes on Flow-induced Capsule Deforma-
tion,” J.Fluid Mech.,460,pp.211–222.
￿23￿ N’Dri,N.A.,Shyy,W.,and Tran-Son-Tay,R.,2003,“Computational Model-
ing of Cell Adhesion and Movement Using a Continuum-Kinetics Approach,”
Biophys.J.,85,pp.2273–2286.
￿24￿ Zhu,C.,1991,“A Thermodynamic and Biomechanical Theory of Cell Adhe-
sion.Part 1:General Formulism,” J.Theor.Biol.150,pp 27–50.
￿25￿ Swenson,J.,Smalley,M.V.,and Hatharasinghe,H.L.M.,1998,“Mechanism
and Strength of Polymer Bridging Flocculation,” Phys.Rev.Lett.,81,pp
5840–5843.
￿26￿ Brooks,D.E.,1973,“The Effect of Neutral Polymers on the Electrokinetic
Potential of Cells and other Charged Particles.III:Experimental Studies on the
Dextran/Erythrocyte System,” J.Colloid Interface Sci.,43,pp.701–713.
￿27￿ Barshtein,G.,Tamir,H.,and Yedgar,S.,1998,“Red Blood Cell Rouleaux
Formation in Dextran Solution:Dependence on Polymer Conformation,” Eur.
Biophys.J.,27,pp.177–181.
￿28￿ Skalak,R.,and Zhu,C.,1990,“Rheological Aspects of Red Blood Cell Ag-
gregation,” Biorheology,27,pp 309–325.
￿29￿ Pozrikidis,C.,1995,“Finite Deformation of Liquid Capsules Enclosed by
Elastic Membranes in Simple Shear Flow,” J.Fluid Mech.,297,pp.123–152.
￿30￿ Ramanujan,S.,and Pozrikidis,C.,1998,“Deformation of Liquid Capsules
Enclosed by Elastic Membranes in Simple Shear Flow:Large Deformations
and the Effect of Fluid Viscosities,” J.Fluid Mech.361,pp.117–143.
￿31￿ Eggleton,C.D.,and Popel,A.S.,1998,“Large Deformation of Red Blood
Cell Ghosts in a Simple Shear Flow,” Phys.Fluids 10,pp.1834–1845.
￿32￿ Derganc,J.,Bozic,B.,Svetina,S.,and Zeks,B.,2003,“Equilibrium Shapes of
Erythrocytes in Rouleau Formation.” Biophys.J.,84,pp.1486–1492.;Torn-
berg,A.-K.,and Shelley,M.,2004,“Simulating the Dynamics and Interactions
of Flexible Fibers in Stokes Flows,” J.Comput.Phys.,196,pp.8–40.
￿33￿ Skalak,R.,and Chien,S.,1983,“Theoretical Models of Rouleau Formation
and Disaggregation,” Ann.N.Y.Acad.Sci.,416,pp.138–148.
￿34￿ Pozrikidis,C.,2003a,“Membrane Theory for Capsules and Cells,￿ in Model-
ing and Simulation of Capsules and Biological Cells,edited by C.Pozrikidis,
Chapman and Hall/CRC Press Mathematical Biology and Medicine Series.
Boca Raton,Florida.
￿35￿ Pozrikidis,C.,2003b,“Deformed Shapes of Axisymmetric Capsules Enclosed
by Elastic Membranes,” J.Eng.Math.,45,pp.169–182.
￿36￿ Chien,S.,Sung,L.A.,Sinchon,S.,Lee,M.M.L.,and Skalak,R.,1984,“En-
ergy Balance in Red Cell Interactions.” Ann.N.Y.Acad.Sci.416,pp.138–
148.
￿37￿ Jeffery,G.B.1922,“The Motion of Ellipsoidal Particles Immersed in a Vis-
cous Fluid,” Proc.R.Soc.London,Ser.A,102,161–170.
￿38￿ Evans,E.A.,and Hochmuth,R.M.,1976,“Membrane Viscoelasticity,” Bio-
phys.J.,16,pp 1–11.
￿39￿ Barthes-Biesel,D.,1991,“Role of Interfacial Properties on the Motion and
Deformation of Capsules in Shear Flow,” Physica A,172,pp.103–124.
￿40￿ Barthes-Biesel,D.,and Sgaier,H.,1985,“Role of Membrane Viscosity in the
Orientation and Deformation of a Spherical Capsule Suspended in Shear
Flow,” J.Fluid Mech.,160,119–136.
￿41￿ Pozrikidis,C.,2001,“Effect of Membrane Bending Stiffness on the Deforma-
tion of Capsules in Simple Shear Flow,” J.Fluid Mech.,440,pp.269–291.
1080/Vol.127,DECEMBER 2005 Transactions of the ASME