Biomech Model MechanobiolDOI 10.1007/s1023701002711
ORIGINAL PAPER
Modeling cell entry into a microchannel
Fong Yew Leong
·
Qingsen Li
·
Chwee Teck Lim
·
KengHwee Chiam
Received: 8 April 2010 / Accepted: 8 November 2010© SpringerVerlag 2010
Abstract
Cell entry into a microchannel has potentialapplications in cell sorting and cancer diagnostics. In thispaper, we numerically model breast cancer cell entry into aconstricted microchannel. Our results indicate that the cellvelocity decreases during entry and increases after entry, anobservation in agreement with experiments. We found thatthe cell entry time depend strongly on the cortical stiffnessand is minimum at some critical cortical elasticity. In addition, we found that for the same entry time, a stiff nucleusis displaced toward the cell front, whereas a viscous nucleusis displaced toward the rear. In comparison, the nucleus isless sensitive to the viscosity of the cytoplasm. These observations suggest that speciﬁc intracellular properties can bededucednoninvasivelyduringcellentry,throughtheinspection of the nucleus using suitable illumination techniques,such as ﬂuorescent labeling.
Keywords
Cell entry
·
Microchannel
·
Immersedboundary method
·
Cell sorting
F. Y. Leong
·
K.H. Chiam (
B
)A*STAR Institute of High Performance Computing,1 Fusionopolis Way, #1616 Connexis,Singapore 138632, Singaporeemail: chiamkh@ihpc.astar.edu.sgQ. Li
·
C. T. LimDivision of Bioengineering, Department of MechanicalEngineering, National University of Singapore,Engineering Drive 1, Singapore 117576, SingaporeC. T. Limemail: ctlim@nus.edu.sg
1 Introduction
Cellular deformation in conﬁned environment has attractedmuch attention, due to increasing applications of microﬂuidic devices in medical diagnostics (Suresh 2007). Cellsforced through microﬂuidic constrictions offer the potentialmeansofquantifyingcellmechanicalcharacteristicsinvitro(Doerschuk et al. 1993; Yap and Kamm 2005). For instance,
diseased cells such as cancer cells are known to have differentstiffnessandelasticitycomparedtotheirhealthycounterparts(Gucketal.2005;LeeandLim2007).Suchdifferences
couldbeusedtodistinguishbetweennormalandcancercells(Kamm 2002; Hou et al. 2009).
Generally, continuum cell models can be classiﬁed intotwo broad categories, namely solid models and liquid dropmodels (Lim et al. 2006). For small cellular deformations, a homogeneous viscoelastic solid model is applicable SchmidSchönbein et al. (1981). But for large cellular
deformations, a better model would be the Newtonian liquiddrop model, which describes the cell as a homogeneous liquiddropenclosedbyacortexwithconstant,isotropictension(Yeung and Evans 1989; Evans and Yeung 1989).
To account for the effects of the nucleus on cell deformation (Dong et al. 1991), the compound drop model wasproposed (Hochmuth et al. 1993), which included an encapsulated liquid drop as a model for the nucleus. The compound drop model effectively explains the observed rapidinitial response in micropipette aspiration and fast recoilon recovery (TranSonTay et al. 1998; Kan et al. 1998).
Applications of the compound drop model include cellrecovery under extensional ﬂows (Kan et al. 1999), micro
pipette aspiration and shear ﬂow (Agresar et al. 1998),cell adhesion, and migration (N’Dri et al. 2003) as well
as shear thinning and membrane elasticity (Marella andUdaykumar2004).Morerecentworksincludeshearinduced
1 3
F. Y. Leong et al.
drop breakup and cellsurface adhesion (Khismatullin andTruskey 2005).Other related studies include the work of Yap and Kamm(2005), who used multipleparticletracking microrheology(Tseng and Wirtz 2001) to probe the mechanical properties of neutrophils within a microfabrication channelnoninvasively (McDonald et al. 2000). The experiment wassimulated(Zhouetal.2007)usingadiffuseinterfacemethod(Yue et al. 2004) and was extended to account for effects dueto the nucleus (Zhou et al. 2008a) and viscoelastic cytoplasm (Zhou et al. 2008b). However, it remains unclear howthenucleusandcytoplasmcanindependentlyaffectcellentryin a conﬁned channel.In the present study, we investigate the entry of breastcancerepithelialcellsintoamicrochannelthroughbothsimulation and experiment. Our objective is to relate observedentry behavior to mechanical properties so as to distinguishbetween benign (MCF10A) and malignant cells (MCF7) (Hou et al. 2009). This work could also be relevantto the study of cancer metastasis as malignant cells aresqueezed into microcapillaries or undergo extravasationinto neighboring tissues.
2 Theory and numerical method
We use the immersed boundary method (Peskin 1977) toanalyze the entry of a cell in a model microchannel. Theimmersed boundary method is well known (Peskin 2002),therefore only a brief outline will be provided here. Detailson a more speciﬁc mixed Eulerian–Lagrangian numericalscheme are available elsewhere (Udaykumar et al. 1997;Shyy et al. 2001).2.1 Model geometryFigure 1 shows a schematic sketch of the experimental setupusingamicrochannelwithasquarecrosssection(150
×
10
×
10
µ
m) and a 45
◦
tapered entrance. The length and width of the taper section are
L
1
and
W
1
and those of the channel section are
L
2
and
W
2
, respectively. The computational domainis segregated into three regions, namely, the extracellularﬂuid (f), the cytoplasm (c), and the nucleus (n). The corticaland nuclear diameters of the cell are
d
c
and
d
n
, respectively.Flowinthemicrochannelisdrivenbyadifferentialpressure
p
established between both ends of the setup.A twodimensional model is used to obtain insights atmodest computational costs. Using an axisymmetric testmodel (see Appendix), we found the cell entry behavior tobe similar to one in a twodimensional model and thereforeinsensitive to depth.
Fig. 1
Schematic of a microchannel with a tapered entrance and thecomputational domain
2.2 Governing equationsAssuming that the ﬂuid is Newtonian and incompressible,the momentum and continuity equations within each distinctregion (i
=
f, c, n) are
ρ
i
∂
u
i
∂
t
+
u
i
∇·
u
i
=−∇
p
i
+
µ
i
∇
2
u
i
+
f
j
(1)
∇·
u
i
=
0 (2)where
u
i
is the ﬂuid velocity,
ρ
i
is the ﬂuid density,
µ
i
isthe ﬂuid dynamic viscosity and
p
i
is the ﬂuid pressure.
f
j
isthe force contributed by the interfacial forces
F
j
at adjacentinterface j:
f
j
=
F
j
·
δ(
x
i
−
x
j
(
s
))
ds
(3)where
δ
isthedeltafunctiontodistributetheinterfacialforcespatially,
x
i
isthepositionvectorsoftheﬂuidgridpoints,and
x
j
isthepositionvectorsofthemarkerpointsalonginterface j whose arc length is s. Material properties such as viscosityare distributed spatially using the Heaviside function. Interfacial markers are advected by the ﬂuid using the followingkinematic relation:
∂
x
j
∂
t
=
u
i
·
δ(
x
i
−
x
j
(
s
))
dx
i
(4)The above system of Eqs. (1–4) are discretized on 2D
staggered grid and solved using the projection methodon a ﬁnite difference scheme. The convection terms aresolved explicitly on Adams–Bashforth scheme, and thediffusion terms are solved implicitly on Crank–Nicholsonscheme. Both schemes are secondorder accurate (N’Driet al. 2003). To ensure adequate resolution for tracking of markers, the horizontal and vertical ﬂuid grid spacing arerestricted to

x
,
y
≤
0
.
05
W
2
. In addition, the initialspacing between the markers is speciﬁed at 0.8 of ﬂuid
1 3
Modeling cell entry into a microchannel
grid spacing. This minimizes the material loss due toinaccuracyoftheimmersedboundarymethodalongtheinterface (Marella and Udaykumar 2004). The error is veriﬁed tobe no more than 1% of the srcinal cell volume.The following boundary conditions are implemented:positive inlet pressure
(
p
in
>
0
)
relative to null outlet pressure
(
p
out
=
0
)
. Noslip condition is imposed on all walls
(
u
wall
=
0
)
. Near the tapered walls, the velocity is approximated as
u
f
−
u
w
=
0 , where
u
f
is the adjacent ﬂuidnode and
u
w
is the adjacent wall node. For other geometrieswhose tapered slope is offdiagonal, the cutcell method isrecommended (Agresar et al. 1998; Udaykumar et al. 1997).
2.3 Interfacial constitutive modelThe compound drop model accounts for two interfaces,namelythecortexandthenuclearmembrane.Thecommonlyassumed model for the cortex is that of constant, isotropictension of a liquid drop (Yeung and Evans 1989; Evans and
Yeung 1989). Here, we include elastic contributions, so thatthe net interfacial force can be expressed as the sum of atension term and an elastic term (Agresar et al. 1998).
F
m
,
c
=
σ
c
·
κ
· ˆ
n
+
η
c
·
λ
·ˆ
t
(5)
F
m
,
n
=
σ
n
·
κ
· ˆ
n
+
η
n
·
λ
·ˆ
t
(6)where
σ
istheinterfacialtensionand
η
istheelasticmodulusfor surface area dilation. The subscripts
c
and
n
refer to thecortex and the nucleus, respectively.
κ
is the interfacial curvature,
λ
is the interfacial strain,
ˆ
n
and
ˆ
t
are the normal andtangential unit vectors, respectively. The vector componentsare evaluated through parameterization of the arc length followed by cubicspline interpolation of the marker points. Asa simpliﬁcation, we assume that both cytoskeletal and membrane properties do not vary in time. Experimental cell entrytimesaretypicallyontheorderof0.1s,atimescalethatmaybe too fast for active processes to remodel the membrane.However,theremaybeotherfactorsthatcausethemembraneproperties to change over the time scale of cell entry. Onesuch scenario could be the slipping of adhesion bonds thatresults in changes to the membrane’s viscoelasticity. However, as we will discuss in Sect. 3.2, the channel surfaces aretreatedwithbovineserumalbumintoavoidadhesion.Therefore, we have neglected transient elastic membrane effectsand assumed that our membrane properties are constant.The cell surface has many folds and wrinkles whichprovide between 80 and 100% excess area (SchmidSchönbein et al. 1980). Based on known stress–strainresponse, nonlinear stiffening only occurs at strains above40%(MarellaandUdaykumar2004).However,suchstiffen
ing was not found to be significant in typical micropipettestudies. We assume that the stress–strain response of thecortex is linear before this limit is reached.In addition, we consider two other known numericalissues. Firstly, the interfacial markers can become eitherclustered or dilated leading to numerical instability (Marellaand Udaykumar 2004). Incidentally, this problem is avoideddue to our inclusion of interfacial elasticity. Secondly, amarker could cross a wall boundary, leading to a violation of nopenetration condition. However, it is assumed that a thinlubricating ﬂuid layer exists between the cell and the wall,and so the cell should not penetrate the wall. To implementthis, virtual compressive springs are introduced between thecell and the wall. This approach is shown to be reasonableby (Dong and Skalak 1992). The equilibrium length of thespring is such that it is between one or two grid spaces toensure that (1) the dimension of the gap is relatively smallcompared to the width of the channel and (2) there is at leastone node to resolve the ﬂuid velocity in that layer. The stiffness issufﬁciently large to ensure that the gap width ismaintained but not so large as to lead to numerical instability. Wefound that, within such limits, the exact magnitude of lengthand stiffness of the springs does not affect the overall rheology significantly. This agrees with the observations madeby (Marella and Udaykumar 2004), who also found that thevirtual spring model does not affect the cell and ﬂuid ﬂowintoamodelmicropipette.Thismethodobviatestheneedtoenforce an artiﬁcial freeslip boundary condition on the wall(Drury and Dembo 2001).2.4 Dimensionless parametersThe characteristic scales for length
l
, pressure
p
, velocity
u
,and time
t
are
L
≡
H
2
,
P
≡
p
·
L L
2
,
U
≡
L
·
P
12
µ
f
,
T
≡
LU
,
(7)where
µ
f
is the viscosity of the suspending ﬂuid. Based onthese characteristic scales, we obtain the following dimensionless parameters,
˜
l
≡
l L
,
˜
p
≡
pP
,
˜
u
≡
uU
,
˜
t
≡
t T
,
(8)
Ca
−
1
c
≡
σ
c
µ
f
·
U
, α
c
≡
d
c
L
, β
c
≡
µ
c
µ
f
, γ
c
≡
η
c
·
L
σ
c
,
(9)
Ca
−
1
n
≡
σ
n
µ
f
·
U
, α
n
≡
d
n
L
, β
n
≡
µ
n
µ
f
, γ
n
≡
η
n
·
L
σ
n
,
(10)where
Ca
−
1
isthereciprocalofthecapillarynumber,
α
istheinterfacialgeometricratio,
β
istheinterfacialviscosityratio,and
γ
is the ratio between interfacial elasticity and tension.Subscripts
c
and
n
refer to the cortex and nuclear membrane,respectively. For the present study, the Reynolds numberis sufﬁciently small
(
Re
≪
1
)
that we do not consider itsinﬂuence on cell entry explicitly.
1 3
F. Y. Leong et al.
3 Results
3.1 Base case modelWe introduce a base case simulation, whose model parameters are
Ca
−
1
c
=
Ca
−
1
n
=
18
,α
c
=
2
,α
n
=
1
.
5
,β
c
=
β
n
=
1, and
γ
c
=
γ
n
=
10
−
4
. The capillary number isbasedonthemeancorticalstiffness(0.01dyn/cm)referencedfrom the range of literature values (0.005–0.035dyn/cm)(Evans and Yeung 1989; Yeung and Evans 1989; Needham
and Hochmuth 1992). The diameters of the cell and nucleusare based on mean experimental values of 20 and 15
µ
m,respectively, with standard deviations of up to 5
µ
m (Fig. 4).Viscosity is initialized as unity for subsequent timescaling(see 3.2). Elasticity tends to be small compared to stiffness,andtheselectedvalueof10
−
4
correspondsto0
.
001dyn
/
cm
2
indimensionalunits.Thisissubsequentlyshowntobeasuitable choice, as we found that an optimal entry time existswithin a range from 10
−
5
to 10
−
3
.The cell is initially stationary, and its center is positionedat a normalized distance of 1.5 from the channel entrance.Figure 2 shows snapshots of cell entry, and Fig. 3 shows the
dependence of cell leading and trailing edge displacement(normalized by
L
) and velocity (normalized by
α
c
·
L
) ontime
˜
t
. The cell entry time
τ
is deﬁned as the time intervalbetween the leading edge and the trailing edge crossing theentrance of the microchannel.After a brief initial transient period (
˜
t
<
0
.
01), the ﬂuidﬂow is nearly steady (Fig. 2a). The cell velocity remainssteady at
˜
t
∼
1
.
5 before decreasing after
˜
t
∼
1
.
5 (Fig. 3b).Midway through the entrance (
˜
t
∼
3
.
9), the cytoplasmicvelocity is observed to uniform (Fig. 2c). Furthermore, thepressure gradient is conﬁned to the cell interior, which suggests a plugged ﬂow behavior (Fig. 2d). Meanwhile, thecell velocity increases until a maximum is reached at
˜
t
∼
6 (Fig. 3b). Beyond the entry region, cell deformation isminimal and is of no further interest.
Fig. 2
Snapshots of cell entry showing velocity (
˜
u
) and pressure (
˜
p
) ﬁelds, using base case parameters:
Ca
−
1
c
=
Ca
−
1
n
=
18,
α
c
=
2,
α
n
=
1
.
5,
β
c
=
β
n
=
1 and
γ
c
=
γ
n
=
10
−
4
1 3
Modeling cell entry into a microchannel
Fig. 3
Cell displacement (
a
) and velocity (
b
) timeplots using base case parameters:
Ca
−
1
c
=
Ca
−
1
n
=
18,
α
c
=
2,
α
n
=
1
.
5,
β
c
=
β
n
=
1, and
γ
c
=
γ
n
=
10
−
4
. Horizontal dotted lines at
˜
l
=
0
.
6 and 2.4 denote the distances of the channel entrance from the initial leading and trailing edges,respectively
Fig. 4
Snapshots of a transfected MCF10A breast epithelial cancercell entering a microchannel. The ﬂuorescentlabeled nucleus is illuminated in (
a
–
c
), whereas the cell membrane outline is distinct underbrightﬁeld (
d
–
f
). The pressure gradient of the microchannel is maintained at constant
p
=
3
.
33Pa
/
µ
m (scale bar=10
µ
m)
3.2 Experimental validationCell entry can be visualized and recorded experimentallyusing a highspeed video camera, as shown in Fig. 4. AMCF10A cell has been pretransfected, so that the nucleusis ﬂuorescent. The nucleus outline is distinct under ﬂuorescence (Fig. 4ac), whereas the cell outline is distinct underbrightﬁeld (Fig. 4df). This allows us to infer the mechanicalproperties of the nucleus and will be further elaborated inSect. 4. To avoid cell adhesion, we use bovine serum albumin (BSA) to prepare the surfaces prior to each experiment.Previously, Zhou et al. (2007) compared their numerical
model with experiments (Yap and Kamm 2005) using an
a priori
speciﬁcation of a viscosity ratio for a real cell (
β
c
=
220), followed by an extrapolation through an assumedpowerlawﬁttedatlowviscosityratios.However,itisunclearhow cell entry depends on the viscosity of the nucleus.Furthermore, cell viscosities are highly variable even for thesame cell type (Yeung and Evans 1989; Evans and Yeung
1989;Hochmuthetal.1993;TranSonTayetal.1998;Need
ham and Hochmuth 1990; Valberg and Albertini 1985; Tsai
et al. 1993).
1 3