Algorithm for a microfluidic assembly line

5.5 MB·

by  Tobias M. Schneider

Algorithm for a microfluidic assembly line
Tobias M. Schneider, Shreyas Mandre, and Michael P. Brenner
School of Engineering and Applied Sciences, Harvard University, 29 Oxford Street, Cambridge, MA 02138
Microfluidic technology has revolutionized the control of flows at small scales giving rise to new
possibilities for assembling complex structures on the microscale. We analyze different possible al-
gorithmsforassemblingarbitrarystructures,anddemonstratethatasequentialassemblyalgorithm
can manufacture arbitrary 3D structures from identical constituents. We illustrate the algorithm
byshowingthatamodifiedHele-Shawcellwith7controlledflowratescanbedesignedtoconstruct
the entire English alphabet from particles that irreversibly stick to each other.
1 Developing novel methods for assembling complex (a) inlets with
1 structures from small particles has been the focus of controlled flowrates
0 much recent investigation [1, 2]. Traditional approaches
2
have mostly revolved around designing selective interac-

n tions between constituent elements of the assembly. The
a structureisthenassembledintheabsenceofanyexternal
J
controlifthermalfluctuationscandrivethesystemtoits
9 free particles
energetic ground state [1, 3], although non-trivial energy
1
landscapes often render this approach challenging.

] Hereweconsideradifferentpossibilityforassemblyon (b) (c)
n
small scales, in which microfluidic flow control is used to
y
steerand assemble smallparticlesintostructuresof high
d
- complexity. The basic idea follows from the observation
u
that if we could construct an arbitrary time dependent
l
f flow field (cid:126)v((cid:126)x,t), then particles in the flow could be ad-
.
s vectedalongarbitrarypathsandmovedtoarbitraryloca-
c
tionsatafixedtime. Thisapparentlyallowstoconstruct
i
s any complex structure, with the individual components
y
binding irreversibly upon contact. 30 1
h
p Of course, the flow field(cid:126)v((cid:126)x,t) cannot be arbitrary; it 20 0.5
[ must conserve mass and momentum [4], 1− 10 1−
v1 ∇·(cid:126)v =0, −∇p+µ∇2(cid:126)v+ρ(cid:126)b=0. (1) /fπaHδτk-100 /fπaHδτk-0.05
1 -20
9 where p is the pressure, µ the liquid viscosity and(cid:126)b is a -30 -1
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
7 volumetricforce. Thequestionofwhatstructurescanbe t/τ t/τ
3 builtthushingesonwhatflowfieldscanbeproducedus- FIG. 1: (Color online) (a) Schematic: Five particles are ad-
. vected by the flow field in a circular cell. The flow field,
1 ingcurrenttechnology[5],andwhatarethelimitsforthe
visualizedbyitsstreamlines,issetupbyelevenflowratesim-
0 possible structures that can form with such flow fields.
posed on the boundary (arrows indicating the strength and
1 Volumetric forces can be produced using e.g. magnetic direction). (b/c)Twotrajectoriestransportingparticlesfrom
1
fields [6] or optical tweezers [7]; alternatively, the flow afixedinitialtotheirdesiredfinalpositionwiththerequired
:
v can be generated by inlets specifying (cid:126)v at the boundary flowrates. The linear trajectory (b) requires flowrates almost
Xi of a cell. 2ordersofmagnitudehigherthanthetrajectoryin(c). This
optimized trajectory minimizes dissipation. Flowrates are
We analyze the case where pressure inlets around the
ar boundary force the flow (Fig. 1), and (cid:126)b = 0. Fluid giveninunitsoftheraterequiredtotransportasingleparticle
in one time unit across the cell.
mechanical constraints prohibit simultaneous control of
manyparticleswiththisdevice;however,wedemonstrate
that a sequential assembly algorithm allows the assem-
bly of arbitrary structures. The device can be designed scribed velocity (cid:126)v((cid:126)x,t), (cid:126)x ∈ ∂Ω. The linearity of (1)
to manufacture the entire English alphabet using 7 con- implies the velocity of the suspended particles is linear
trolled flow rates in two dimensions. in the boundary forcing [4], i.e.
Linear Response to Boundary Forcing. Consider N
particles suspended in the flow domain Ω with their in-
N (cid:90)
stantaneous positions being (cid:126)xj, j = 1,2,...,N. Let the (cid:88)Rjk(cid:126)x˙k = Kj((cid:126)x1,x(cid:126)2,...x(cid:126)N;ξ(cid:126))(cid:126)v(ξ(cid:126),t) dS+F(cid:126)j,(2)
flow be forced on the boundary of the flow cell by a pre- k=1 ∂Ω
2
where F(cid:126) is the force acting on the jth particle, possibly in the far field.
j
duetonon-hydrodynamicinterparticleinteractions. The WeconsiderthedevicedepictedinFig. 1(a),acircular
responsecoefficientsK andR dependonthegeometry domain of radius a with flowrates f prescribed at the
j jk k
of the flow cell [8] and can be computed numerically. If, boundary at M discrete inlets with positions R(cid:126) . The
k
the boundary forcing occurs through M discrete inlets velocity field at any position (cid:126)x is then given by
of area ∆S, located at ξ(cid:126) with prescribed velocities (cid:126)v ,
k = 1,2,...,M, then Ekqn. (2) can be written in thke (cid:126)u((cid:126)x)=− 1 (cid:88)M (cid:126)x−R(cid:126)k f ≡B((cid:126)x)·f. (4)
discrete form πH ((cid:126)x−R(cid:126) )2 k
k=1 k
R(x) x˙ =M(x)·f +F, (3) Thus, the matrix M in Eqn. (3) may be constructed
by combining the N position dependent matrices B(x(cid:126))
i
with x = [(cid:126)x1,(cid:126)x2,...,(cid:126)xN], flowrates f = corresponding to each particle.
[(cid:126)v ,(cid:126)v ,...,(cid:126)v ]∆S and M =K ((cid:126)x ,x(cid:126) ,...x(cid:126) ;ξ(cid:126) ). We first test whether this flow cell will allow simulta-
1 2 M jk j 1 2 N k
neous control of all N particles by boundary inlets. The
Eqn. (3) is an instantaneous linear relation between
flow rates scale inversely with the duration of the assem-
the imposed flowrates and the particle velocities. This
bly so that the scale for the flux depends on the chosen
implies that prescribed particle trajectories x(t), may be
timescale. Wescaleflowratebythefluxrequiredtomove
realizedbyimposingsuitablef obtainedbyinverting(3).
a single particle across the cell F = πaHδ/τ, where δ is
Thefeasibilityofthismethodthenhingescruciallyupon
the width of the inlet and the duration τ of the pro-
invertibility of M, which in general needs not be square.
cess. Fig. 1(b) shows that if we require the particles to
Condition for Controlling Particle Trajectories. A
betransportedinstraightlinesatconstantvelocityfrom
necessary condition for inverting M is that the number
their initial position to their final position, the required
of independent controls exceed the number of degrees of
flowrates reach up to 30 times this value. This happens
freedom. With N particles in 3 dimensions, there must
because when two or more particles are moving towards
be at least M = 3N +1 flow inlets; 3N inlets control
each other a distance (cid:15) apart, a strain rate field of (cid:15)˙/(cid:15) is
particle degrees of freedom, with the additional inlet en-
required. If the approach velocity is constant, the strain
forcingvolumeconservation. Similarlyintwodimensions
diverges as (cid:15) → 0, implying diverging flowrates. Thus,
at least 2N +1 flow inlets are required.
whenever particles are brought together at constant ve-
Ingeneral, thesealgebraicconditionsarenotsufficient
locity, large fluxes are required.
for the practicality of this assembly method. We will see
OptimizedTrajectories. Wecantrytoreducetheflow
belowthatinpractice,theflowratesrequiredtoindepen-
ratesbychoosingtheparticletrajectoriesandspeedscon-
dentlysteerlargenumberofparticlescanbetoolargeto
necting the initial and the final states to minimize this
bepractical,owingtothepoorconditioningofthematrix
effect. Intuitively, we can decrease the speed of the par-
M.
ticles as they approach each other to minimize the re-
Hele-Shaw Flow as a Specific Example. To explore
quired flowrates. To find out if this is sufficient we com-
this in more detail, we consider a specific example. We
putetheoptimaltrajectoriesconnectingtheinitialtothe
consider the fluid mechanics of typical microfluidic de-
desired particle configuration by finding the trajectories
vices [9] where the vertical gap thickness H is much
that minimize dissipation. The dissipation rate is given
smallerthanthelateralextent. InHele-Shawapproxima-
by
tion [8], the velocity profile across the gap is a parabola,
with the gap-averaged velocity (cid:126)u being proportional to H2 (cid:90)
w = |∇p|2dA. (5)
the gradient of the pressure field p satisfying Laplace’s
12µ
A
equation. A particle at position (cid:126)x responds to the flow
k
by moving at a speed proportional to the local fluid ve- For the discretized forcing w can be expressed as a
locity, (cid:126)x˙ = β(cid:126)u((cid:126)x ), with β depending on the particle quadraticformintheflowratesw =f†·D·f withmetric
k k
sizeandshape[8]. Thislineardependenceoftheparticle D beinganalyticallygivenintermsoftheinletpositions.
velocity on the gap-averaged fluid velocity is a general We minimize the dissipation and thereby the flowrates
consequence of the linearity of Stokes flow. The quanti- under the constraints that the dynamics move the parti-
tativevalueforβcanbecalculatedintwodifferentlimits: cles according to their equation of motion Eqn. (3) from
Iftheparticleismuchsmallerthanthegap,itisadvected their chosen initial state to their final state. This re-
by the local velocity, so β depends on the location of the quires that we find the trajectories that minimize the
particle relative to the walls. Alternatively, if the par- Lagrangian
toifclsecaalpepurnodxiemrlaytinelgytshpeaHnselteh-Sehgaawpawpipdrtohx,itmhaetsioenpabrraetaiokns L=(cid:90) 1dt(cid:110)f†·D·f −λ†·(cid:2)x˙ −M ·f(cid:3)−γ e†·f(cid:111) (6)
down near the particle. The factor β can now be calcu- 0
latedbysolvingtheStokesequationsclosetotheparticle, where λ and γ are Lagrange multipliers to enforce the
matchingthesolutionwiththeparabolicHele-Shawflow constraints and e† = (1,...,1). λ enforces the equation
3
of motion, whereas γ requires that the fluxes f satisfy particle positions in the aggregate’s frame of reference
volume conservation. and R is the two-dimensional rotation. How many de-
TominimizeLweconsiderasmallperturbationofthe grees of freedom are required for the assembly? Control-
flowrates in the direction of steepest descent, i.e. f → ling the position of the aggregate requires two degrees
f +∆f with of freedom. Rotating it requires a stagnation point flow
superimposedonthelocalmeanflowneartheaggregate.
(cid:18)δL(cid:19)† (cid:110) (cid:111) Thisrequiresindependentcontroloftwomoredegreesof
∆f =−(cid:15) =−(cid:15) 2D·f +M†·λ−γe . (7)
δf freedom, the strength and orientation of the stagnation
point flow. Finally, the free particle location demands
To fix the unknown fields γ and λ we require the vari-
twoadditionaldegreesoffreedom. Includingvolumecon-
ation of L with respect to x,λ and γ to vanish for the
servation we thus need 7 inlets to absolutely control all
updated flowrates and trajectory x → x + ∆x. Eval-
degrees of freedom.
uating the stationarity condition at linear order in the
Sequentially adding particles allows to assemble struc-
increments allows to eliminate γ from Eqn. (7) and se-
tures of arbitrary shape and thus spell a word (Fig. 2)
lects a unique solution of the adjoint equations
withoutanyfeedbackcontrol. Technically,thetrajectory
∂M for a particle to be added is constructed by defining the
λ˙ =−f mnλ (8)
k n ∂x m desiredpositionwhichistakenfrom[10]andthedirection
k
of approach in the frame of reference of the aggregate.
whicharenumericallysolvedtogetherwiththeevolution
Spline interpolation between the initial position of the
of x using MATLAB. [13]
newparticleanditsfinalpositionwhentheaggregatehas
The minimization starts from the linear trajectory in
reachedthedesiredconfiguration,thenyieldstheparticle
Fig.1(b)alongwhichwechoosevelocityvariationsdamp-
path. The aggregate can also be moved arbitrarily; we
ing large peak flowrates. Panel (c) shows the same 5
choose to keep its orientation fixed and move the center
particles following an optimized trajectory for which the
of mass to the device center in each step. The velocity
flowrates are reduced by more than an order of magni-
alongthistrajectoryischosensuchthatitvanishesatthe
tude,tothedesiredrange. Thereductionisduetochoos-
ing (cid:15)˙/(cid:15) to approach a finite number as (cid:15)→0. Thus, op-
timizing trajectories can lead to simultaneous control of
a modest number of particles using boundary flowrates.
Nonetheless, this approach does not scale to larger
number of particles; e.g, we could not find trajectories
allowing simultaneous control of 13 particles required to
spell the letter ‘B’ at moderate fluxes. This reflects a
implementation independent physical limit of directing
flow fields with boundary fluxes, resulting from the fact
that flow modes forced at the boundary decay in ampli-
tude as one moves away from the wall. The character-
istic length scale governing the decay of the boundary
modesisthedistancebetweentheinjectionpointswhich
scales inversely with the number of inlets. But as the
particlenumberincreases,weneedmoreinletstocontrol
the particle motion. This limits the number of particles
that can be controlled simultaneously. In our numerical
experiments we could not find trajectories transporting
more than 6 individual particles at moderate flowrates.
Sequential Assembly. To overcome this physical limi-
tation, an algorithm is required that decouples the num-
ber of controlled degrees of freedom from the num-
ber of particles in the desired structure. This can be
achieved with a sequential approach, where one parti-
cle after the other is attached to an aggregate which
moves as a rigid body subject to force- and torque bal- FIG. 2: (Color online) Sequentially particles are added
to form aggregates of arbitrary shape. The three columns
ance. Construction of the corresponding M matrix re-
presentsnapshotsalongtheassemblyofthreeexamplestruc-
quires explicitly accounting for the translation and rota-
tures presented in the last row. The required time-variation
tion of the cluster consisting of Z particles at positions of the seven controlled flowrates are computed a priori (eg.
(cid:16) (cid:17)
(cid:126)x =(cid:126)x +R(ϕ) ξ(cid:126) −ξ(cid:126) whereξ(cid:126) aretheprescribed Fig. 3). See movie M1 (online) for an animation.
α cm α cm α
4
1
0 plement feedback control [12], though this significantly
-1
1 complicates the process especially for three dimensional
0
-1 structures. Another intriguing possibility is to simply
10 embrace the existence of stochasticity and construct the
-1
most stable trajectories that maximize the probability
1
0
-1 of formation of the desired structures in the presence of
10 noise. Finding algorithmic methods for carrying out this
-1
optimizationareimportantdirectionsforfutureresearch.
1
0
-1
1
0 We acknowledge support by the National Science
-1
0 1 2 3 4 5 6 7 8 9 10 11 12 13 Foundation; the Kavli Institute for Bionano Science and
t
FIG. 3: (Color online) Flowrates fk/πaHδτ−1 for the seven Technology at Harvard; and the German Science Foun-
inletsasafunctionoftimet/τ. Theseallowtospelltheletter dation through grant Schn 1167/1 (TMS).
‘B’. At each integer time a new particles enters the cell.
initialandfinalconfigurationsothatflowratesreachzero
aftereachassemblystep. Withthischoiceoftrajectories,
[1] G. Whitesides and B. Grzybowski, Science 295, 2418
the flowrates are computed inverting Eqn. (3). No ex- (2002).
tensive optimization was required to limit the fluxes (see [2] H. - Y. J. Yeh, and J. S. Smith, IEEE Photon.
Fig.3)toreasonablevaluesofthesameorderasthosein Technol. Lett. 6, 706 (1994); G. Whitesides and
Fig. 1(c). Optimizing trajectories might however still be M. Boncheva, Proc. Natl. Acad. Sci. USA 99, 4769
(2002); R. R. A. Syms, et al., J. Microelectromech. Sys.
useful either to further reduce the required flowrates or
12,387(2003);C.Dobson,Seminarsincellanddevelop-
to minimize internal forces within the growing aggregate
mentalbiology15,3(2004);M.BonchevaandG.White-
so that the forming structure can sustain the stresses in
sides,MRSBull30,736(2005);G.Meng,etal.,Science
the flow. 327, 560 (2010).
We envision a typical mode of operation of the mi- [3] M. Boncheva, et al., Proc. Natl. Acad. Sci. USA 102,
crofluidic device similar to that of a macroscopic robotic 3924 (2005); M. Fialkowski, A. Bitner, and B. Grzy-
assembly line: Once the constituents of an assembly bowski, Nat Mater 4, 93 (2005); P. W. K. Rothemund,
Nature 440, 297 (2006); S. Douglas, et al., Nature 459,
and the detailed assembly sequence are decided on, the
414 (2009).
flowrates on the inlets of the device and the precise in-
[4] J. Happel and H. Brenner, Low Reynolds number hydro-
stants at which each constituent is introduced are cal-
dynamics: with special applications to particulate media
culated. These pre–computed flowrates are then set up (Kluwer, 1991).
in a time–periodic fashion, and a train of the assembly [5] S.QuakeandA.Scherer,Science290,1536(2000);B.D.
constituents are introduced at the inlets of the device to Gates, et al., Annu. Rev. Mater. Res. 34, 339 (2004);
accomplish a train of assembled products. T.SquiresandS.Quake,Rev.Mod.Phys.7,977(2005);
G.Whitesides,Nature442,368(2006);J.-H.Jang,etal.,
In summary, we have shown that microfluidic assem-
Angew. Chem. Int. Ed. 46, 9027 (2007); D. Dendukuri
bly can be an efficient strategy if structures are built
andP.S.Doyle,Adv.Mater.21,1(2009);D.K.Hwang,
sequentially. Usingthisapproach,weshowthatthetem-
D.Dendukuri,andP.S.Doyle,LabChip8,1640(2008);
poral control of only 7 flowrates allows to build arbi- K. Sung, et al., J. Am. Chem. Soc 130, 1335 (2008).
trarily shaped particle aggregates in 2 dimensions. In [6] M.Jullien,J.Paret,andP.Tabeling,Phys.Rev.Lett.82,
3 dimensions, the same argument implies that 11 dif- 2872 (1999); N. Ouellette, P. O’ Malley, and J. Gollub,
ferent flowrates are required. Different chamber geome- Phys. Rev. Lett. 101, 174504 (2008).
[7] E. Dufresne and D. Grier, Rev. Scient. Inst. 69, 1974
triesandhydrodynamicinteractionscanbeincorporated
(1998); D. Grier, Nature 424, 810 (2003).
into this framework, which rests solely on the linearity
[8] C. Pozrikidis, J. Fluid Mech. 261, 199 (1994).
of Stokes flow. We anticipate that similar algorithms
[9] Y. Xia and G. M. Whitesides, Annu. Rev. Mater. Sci.
can be constructed using other forcing mechanisms such 28, 153 (1998); M. Hashimoto, et al., Soft Matter 4,
as electrokinetics [11], where non-hydrodynamic electri- 1403 (2008).
cal forcing contributions need to be included into the [10] D. M. Eigler and E. K. Schweizer, Nature 344, 524
transfer matrices. A challenge for practical implementa- (1990).
[11] A. E. Cohen, Phys. Rev. Lett. 94, 118102 (2005).
tion is to quantify the sensitivity of particle trajectories
[12] A. E. Cohen and W. E. Moerner, Proc Natl Acad Sci
to noise in the imposed flow rates. This is particularly
USA 102, 4362 (2006).
relevant when scaling down the assembly to submicron
[13] Since the trajectory increment preserves constraints to
scale, where the Pecl´et number corresponding to Brow- linear order only, we explicitly enforce in each iteration
nian motion is small and hence potentially disruptive. that particles reach their prescribed final positions ex-
One option for dealing with errors and noise is to im- actly. This corrects for neglected higher order effects.