A Waste-Efficient Algorithm for Single-Droplet Sample Preparation on Microfluidic Chips
AA Waste-Efficient Algorithm for Single-DropletSample Preparation on Microfluidic Chips
Miguel Coviello Gonzalez
Department of Computer ScienceUniversity of California at Riverside
Marek Chrobak ? Department of Computer ScienceUniversity of California at Riverside
Abstract
We address the problem of designing micro-fluidic chips for sample preparation, which is a crucialstep in many experimental processes in chemical and biological sciences. One of the objectives ofsample preparation is to dilute the sample fluid, called reactant, using another fluid called buffer,to produce desired volumes of fluid with prespecified reactant concentrations. In the model weadopt, these fluids are manipulated in discrete volumes called droplets. The dilution processis represented by a mixing graph whose nodes represent 1-1 micro-mixers and edges representchannels for transporting fluids. In this work we focus on designing such mixing graphs whenthe given sample (also referred to as the target ) consists of a single-droplet, and the objective isto minimize total fluid waste. Our main contribution is an efficient algorithm called
RPRIS thatguarantees a better provable worst-case bound on waste and significantly outperforms state-of-the-art algorithms in experimental comparison.
Discrete Mathematics → Combinatorial Optimization • Theoryof Computation
Keywords and phrases algorithms, graph theory, lab-on-chip, fluid mixing
Digital Object Identifier
Funding
Research supported by NSF grant CCF-1536026.
Microfluidic chips are miniature devices that can manipulate tiny amounts of fluids on asmall chip and can perform, automatically, various laboratory functions such as dispensing,mixing, filtering and detection. They play an increasingly important role in today’s scienceand technology, with applications in environmental or medical monitoring, protein or DNAanalysis, drug discovery, physiological sample analysis, and cancer research.These chips often contain modules whose function is to mix fluids. One applicationwhere fluid mixing plays a crucial role is sample preparation for some biological or chemicalexperiments. When preparing such samples, one of the objectives is to produce desiredvolumes of the fluid of interest, called reactant , diluted to some specified concentrations bymixing it with another fluid called buffer . As an example, an experimental study may requirea sample that consists of 6 µL of reactant with concentration 10%, 9 µL of reactant withconcentration 20%, and 3 µL of reactant with concentration 40%. Such multiple-concentrationsamples are often required in toxicology or pharmaceutical studies, among other applications.There are different models for fluid mixing in the literature and multiple technologiesfor manufacturing fluid-mixing microfluidic chips. (See the survey in [2] or the recent © The copyright is retained by the authors;licensed under Creative Commons License CC-BYLeibniz International Proceedings in InformaticsSchloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany a r X i v : . [ c s . D S ] N ov X:2 Mixing Graphs book [1] for more information on different models and algorithmic issues related to fluidmixing.) In this work we assume the droplet-based model, where the fluids are manipulatedin discrete quantities called droplets . For convenience, we will identify droplets by theirreactant concentrations, which are numbers in the interval [0 ,
1] with finite binary precision.In particular, a droplet of reactant is denoted by 1 and a droplet of buffer by 0. We focus onthe mixing technology that utilizes modules called . A micro-mixer has twoinlets and two outlets. It receives two droplets of fluid, one in each inlet, mixes these dropletsperfectly, and produces two droplets of the mixed fluid, one on each outlet. (Thus, if theinlet droplets have reactant concentrations a and b , then the two outlet droplets each willhave concentration ( a + b ).) Input droplets are injected into the chip via droplet dispensersand output droplets are collected in droplet collectors. All these components are connectedvia micro-channels that transport droplets, forming naturally an acyclic graph that we calla mixing graph , whose source nodes are fluid dispensers, internal nodes (of in-degree andout-degree 2) are micro-mixers, and sink nodes are droplet collectors. Graph G in Figure 1illustrates an example of a mixing graph. w w G G Figure 1
On the left, a mixing graph G that produces droplet set (cid:8) , , , , , (cid:9) frominput set I = { , , , , , } . Numbers on the micro-mixers (internal nodes) represent dropletconcentrations produced by these micro-mixers. If only some of the produced droplets are needed,the remaining droplets are designated as waste. This is illustrated by the mixing graph G thatproduces droplets (cid:8) , , , (cid:9) . Small black circles labeled “w” on micro-mixers represent dropletsof waste. Given some target set of droplets with specified reactant concentrations, the objective isto design a mixing graph that produces these droplets from pure reactant and buffer droplets,while optimizing some objective function. Some target sets can be produced only if we allowthe mixing graph to also produce some superfluous amount of fluid that we refer to as waste ;see graph G in Figure 1. One natural objective function is to minimize the number of wastedroplets (or equivalently, the total number of input droplets). As reactant is typically moreexpensive than buffer, one other common objective is to minimize the reactant usage. Yetanother possibility is to minimize the number of micro-mixers or the depth of the mixinggraph. There is growing literature on developing techniques and algorithms for designingsuch mixing graphs that attempt to optimize some of the above criteria. State-of-the-art . Most of the earlier papers on this topic studied designing mixing graphs forsingle-droplet targets. This line of research was pioneered by Thies et al. [10], who proposedan algorithm called
Min-Mix that constructs a mixing graph for a single-droplet target withthe minimum number of mixing operations. Roy et al. [9] developed an algorithm called . Coviello Gonzalez and M. Chrobak XX:3
DMRW designed to minimize waste. Huang et al. [6] considered minimizing reactant usage, andproposed an algorithm called
REMIA . Another algorithm called
GORMA , for minimizing reactantusage and based on a branch-and-bound technique, was developed by Chiang et al. [3].The algorithms listed above are heuristics, with no formal performance guarantees. Aninteresting attempt to develop an algorithm that minimizes waste, for target sets withmultiple droplets, was reported by Dinh et al. [4]. Their algorithm, that we refer to as
ILP ,is based on a reduction to integer linear programming and, since their integer program couldbe exponential in the precision d of the target set (and thus also in terms of the input size),its worst-case running time is doubly exponential. Further, as this algorithm only considersmixing graphs of depth at most d , it does not always finds an optimal solution (see anexample in [5]). In spite of these deficiencies, for very small values of d it is still likely toproduce good mixing graphs.Additional work regarding the design of mixing graphs for multiple droplets includesHuang et al. ’s algorithm called WARA , which is an extension of Algorithm
REMIA , that focuseson reactant minimization; see [7]. Mitra et al. [8] also proposed an algorithm for multipledroplet concentrations by modeling the problem as an instance of the Asymmetric TSP on ade Bruijn graph.As discussed in [5], the computational complexity of computing mixing graphs withminimum waste is still open, even in the case of single-droplet targets. In fact, it is not evenknown whether the minimum-waste function is computable at all, or whether it is decidableto determine if a given target set can be produced without any waste. To our knowledge,the only known result that addresses theoretical aspects of designing mixing graphs is apolynomial-time algorithm in [5] that determines whether a given collection of droplets withspecified concentrations can be mixed perfectly with a mixing graph.
Our results . Continuing the line of work in [10, 9, 6, 3], we develop a new efficientalgorithm
RPRIS (for
Recursive Precision Reduction with Initial Shift ) for designing mixinggraphs for single-droplet targets, with the objective to minimize waste. Our algorithm wasdesigned to provide improved worst-case waste estimate; specifically to cut it by half formost concentrations. Its main idea is quite natural: recursively, at each step it reduces theprecision of the target droplet by 2, while only adding one waste droplet when adjusting themixing graph during backtracking.While designed with worst-case performance in mind,
RPRIS significantly outperformsalgorithms
Min-Mix , DMRW and
GORMA in our experimental study, producing on average about50% less waste than
Min-Mix , between 21 and 25% less waste than
DMRW (with the percentageincreasing with the precision d of the target droplet), and about 17% less waste than GORMA .(It also produces about 40% less waste than
REMIA .) Additionally, when compared to
ILP , RPRIS produces on average only about 7% additional waste.Unlike earlier work in this area, that was strictly experimental, we introduce a performancemeasure for waste minimization algorithms and show that
RPRIS has better worst-caseperformance than
Min-Mix and
DMRW . This measure is based on two attributes d and γ ofthe target concentration t . As defined earlier, d is the precision of t , and γ is defined as thenumber of equal leading bits in t ’s binary representation, not including the least-significantbit 1. For example, if t = . γ = 4, and if t = . γ = 3. (Both d and γ are functions of t , but we skip the argument t , as it is always understood from context.) Inthe discussion below we provide more intuition and motivations for using these parameters.We show that Algorithm RPRIS produces at most ( d + γ ) + 2 droplets of waste (seeTheorem 1 in Section 5). In comparison, Algorithm Min-Mix from [10] produces exactly d droplets of waste to produce t , independently of the value of t . This means that the waste X:4 Mixing Graphs of RPRIS is about half that of
Min-Mix for almost all concentrations t . (More formally, fora uniformly chosen random t with precision d the probability that the waste is larger than( − (cid:15) ) d vanishes when d grows, for any (cid:15) > DMRW , its average waste isbetter than that of
Min-Mix , but its worst-case bound is still d − O (1) even for small valuesof γ (say, when t ∈ [ , ]), while Algorithm RPRIS ’ waste is at most d/ O (1) in this range.In regard to time performance, for the problem of computing mixing graphs it wouldbe reasonable to express the time complexity of an algorithm as a function of its output,which is the size of the produced graph. This is because the output size is at least as largeas the input size, which is equal to d – the number of bits of t . (In fact, typically it’s muchlarger.) Algorithm RPRIS runs in time that is linear in the size of the computed graph, andthe graphs computed by Algorithm
RPRIS have size O ( d ). Discussion . To understand better our performance measure for waste, observe that theoptimum waste is never smaller than γ + 1. This is because if the binary representation of t starts with γ γ + 1 input droplets 0 and at least onedroplet 1. (The case when the leading bits of t are 1’s is symmetric.) For this reasons, anatural approach is to express the waste in the form γ + f ( d − γ ), for some function f (). InAlgorithm RPRIS we have f ( x ) ≈ x . It is not known whether smaller functions f () can beachieved.Ideally, one would like to develop efficient “approximation” algorithms for waste minimiza-tion, that measure waste performance in terms of the additive or multiplicative approximationerror, with respect to the optimum value. This is not realistic, however, given the currentstate of knowledge, since currently no close and computable bounds for the optimum wasteare known. We use notation prec ( c ) for the precision of concentration c , that is the number of fractionalbits in the binary representation of c . (All concentration values will have finite binaryrepresentation.) In other words, prec ( c ) = d ∈ Z ≥ such that c = a/ d for an odd a ∈ Z .We will deal with sets of droplets, some possibly with equal concentrations. We define a configuration as a multiset of droplet concentrations. Let A be an arbitrary configuration.By | A | = n we denote the number of droplets in A . We will often write a configurationas A = { f : a , f : a , ..., f m : a m } , where each a i represents a different concentration and f i denotes the multiplicity of a i in A . (If f i = 1, then, we will just write “ a i ” instead of“ f i : a i ”.) Naturally, we have P mi =1 f i = n .We defined mixing graphs in the introduction. A mixing graph can be thought of,abstractly, as a linear mapping from the source values (usually 0’s and 1’s) to the sink values.Yet in the paper, for convenience, we will assume that the source concentration vector is partof a mixing graph’s specification, and that all sources, micro-mixers, and sinks are labeled bytheir associated concentration values.We now define an operation of graph coupling. Consider two mixing graphs G and G . Let T be the output configuration (the concentration labels of the sink nodes) of G and I be the input configuration (the concentration labels of the source nodes) for G . Toconstruct the coupling of G and G , denoted G • G , we identify inlet edges of the sinksof G with labels from T ∩ I with outlet edges of the corresponding sources in G . Moreprecisely, repeat the following steps as long as T ∩ I = ∅ : (1) choose any a ∈ T ∩ I , (2)choose any sink node t of G labeled a , and let ( u , t ) be its inlet edge, (3) choose anysource node s of G labeled a , and let ( s , v ) be its outlet edge, (4) remove t and s and . Coviello Gonzalez and M. Chrobak XX:5 their incident edges, and finally, (5) add edge ( u , v ). The remaining sources of G and G become sources of G • G , and the remaining sinks of G and G become sinks of G • G .See Figure 2 for an example. G G G G
14 123838 3814 145858 12 1212 38 121458 381258 12 12 38
Figure 2
Coupling of two mixing graphs G and G . G • G is obtained by identifying inletedges of two sinks of G , one labelled and one , with the outlet edges of the correspondingsources of G . These new edges are shown as dotted arrows. Next, we define converter graphs. An ( i : α, j : β ) -converter is a mixing graph thatproduces a configuration of the form T = { i : α, j : β } ∪ W , where W denotes a set of wastedroplets, and whose input droplets have concentration labels either 0 or 1. As an example,graph G in Figure 1 can be interpreted as a (2 : , )-converter that produces two wastedroplets of concentrations and .If needed, to avoid clutter, sometimes we will use a more compact graphical representationof mixing graphs by aggregating (not necessarily all) nodes with the same concentrationlabels into a single node, and with edges labeled by the number of droplets that flow throughthem. (We will never aggregate two micro-mixer nodes if they both produce a droplet ofwaste.) If the label of an edge is 1, then we will simply omit the label. See Figure 3 for anexample of such a compact representation.
14 12 34 12
14 12 3412 1214 3414 14 34 34
14 12 34121214 3422 22 2 223 3 G G Figure 3 G is a compact representation of G . All nodes in G (except the last intermediatenode with label ) represent an aggregation of at least two nodes from G . X:6 Mixing Graphs
In this section, we describe our algorithm
RPRIS for producing a single-droplet target ofconcentration t with precision d = prec ( t ). We first give the overall strategy and then wegradually explain its implementation. The core idea behind RPRIS is a recursive procedurethat we refer to as
Recursive Precision Reduction , that we outline first. In this procedure, t s denotes the concentration computed at the s th recursive step with d s = prec ( t s ); initially, t = t . Also, by B we denote the set of base concentration values with small precision forwhich we give explicit mixing graphs later in this section. Procedure
RPR( t s ) If t s ∈ B , let G s be the base mixing graph (defined later) for t s , else :(rpr1) Replace t s by another concentration value t s +1 with d s +1 = d s − G s +1 for t s +1 .(rpr3) Convert G s +1 into a mixing graph G s for t s , increasing waste by one droplet. Return G s .The mixing graph produced by this process is G .When we convert G s +1 into G s in part (rpr3), the precision of the target increases by 2,but the waste only increases by 1, which gives us a rough bound of d/ t ∈ [ , ]. To deal with values outside this interval, we map t into t so that t ∈ [ , ],next we apply Recursive Precision Reduction to t , and then we appropriately modify thecomputed mixing graph. This process is called Initial Shift .We next describe these two processes in more detail, starting with Recursive PrecisionReduction, followed by Initial Shift.
Recursive Precision Reduction (RPR) . We start with concentration t that, by applyingInitial Shift (described next), we can assume to be in [ , ]. Step (rpr1): computing t s +1 . We convert t s into a carefully chosen concentration t s +1 forwhich d s +1 = d s −
2. One key idea is to maintain an invariant so that at each recursivestep, this new concentration value t s +1 satisfies t s +1 ∈ [ , ]. To accomplish this, weconsider five intervals S = [ , ], S = [ , ], S = [ , ], S = [ , ], and S = [ , ]. Wechoose an interval S k that contains t s “in the middle”, that is S k = [ l, r ] for k such that t s ∈ [ l + , r − ]. (See Figure 4.) We then compute t s +1 = 4( t s − l ). Note that t s +1 satisfies both t s +1 ∈ [ , ] (that is, our invariant) and d s +1 = d s −
18 14 38 12 58 34 78 S S S S S Figure 4
Graphical representation of intervals S , S , . . . , S . The thick shaded part of eachinterval S k = [ l, r ] marks its “middle section” [ l + , r − ]. Each concentration within interval[ , ] belongs to the middle section of some S k . Step (rpr3): converting G s +1 into G s . Let G s +1 be the mixing graph obtained for t s +1 instep (rpr2). We first modify G s +1 to obtain a graph G s +1 , which is then coupled with an . Coviello Gonzalez and M. Chrobak XX:7 appropriate converter C s +1 to obtain mixing graph G s = G s +1 • C s +1 . Figure 5 illustratesthis process. G s G s+ t s+ waste I s+ … … G s+ t s waste I s+ … … G s+ t s waste I s+ … … C s+ … I s… waste ’ ’’ ’ Figure 5
Conversion from G s +1 to G s . The left image illustrates the computed mixing graph G s +1 with input labels I s +1 (consisting of only 0’s and 1’s) that produces t s +1 along with somewaste. The middle figure illustrates G s +1 , which is obtained from G s +1 by changing concentrationlabels. The last figure illustrates the complete mixing graph G s = G s +1 • C s +1 for t s , shown withina dotted rectangle. Next, we explain how to construct G s +1 . G s +1 consists of the same nodes and edges as G s +1 , only the concentration labels are changed. Specifically, every concentration label c from G s +1 is changed to l + c/ G s +1 . Note that this is simply the inverse of the linearfunction that maps t s to t s +1 . In particular, this will map the 0- and 1-labels of the sourcenodes in G s +1 to the endpoints l and r of the corresponding interval S k .The converter C s +1 used in G s needs to have sink nodes with labels equal to the sourcenodes for G s +1 . That is, if the labeling of the source nodes of G s +1 is I s +1 = { i : l, j : r } ,then C s +1 will be an ( i : l, j : r )-converter. As a general rule, C s +1 should produce atmost one waste droplet, but there will be some exceptional cases where it produces two.(Nonetheless, we will show that at most one such “bad” converter is used during the RPRprocess.) The construction of these converters is somewhat intricate, and is deferred to thenext section. The base case . We now specify the set of base concentration values and their mixinggraphs. Let B = (cid:8) , , , , , , (cid:9) . (Concentrations and are not strictly necessaryfor correctness but are included in the base case to improve the waste bound.) Figure 6illustrates the mixing graphs for concentrations , , , and ; the mixing graphs for theremaining concentrations are symmetric. Initial Shift (IS) . We now describe the IS procedure. At the fundamental level, the idea issimilar to a single step of RPR, although the involved linear mappings and the converter aresignificantly different.We can assume that t < (because for t > the process is symmetric). Thus thebinary representation of t starts with γ ≥ γ − t ∈ [ , ), we coulduse this value as the result of the initial shift, but to improve the waste bound we refinethis choice as follows: If 2 γ − t ∈ ( , ) then let t = 2 γ − t and σ = 1. Otherwise, we have2 γ − t ∈ [ , ], in which case we let t = 2 γ t and σ = 0. In either case, t = 2 γ − σ t ∈ [ , ]and d = d − γ + σ . X:8 Mixing Graphs
14 1214
14 123838
14 1238516516 B B B B
4w ww w w ww
Figure 6
Base mixing graphs B , B , B and B for concentrations , , and , respectively. Let G be the mixing graph obtained by applying the RPR process to t . It remains toshow how to modify G to obtain the mixing graph G for t . This is analogous to the processshown in Figure 5. We first construct a mixing graph G that consists of the same nodes andedges as G , only each concentration label c is replaced by c/ γ − σ . In particular, the labelset of the source nodes in G will have the form I = { i : 0 , j : 1 / γ − σ } . We then construct a( i : 0 , j : 1 / γ − σ )-converter C and couple it with G to obtain G ; that is, G = G • C . This C is easy to construct: The 0’s don’t require any mixing, and to produce the j droplets1 / γ − σ we start with one droplet 1 and repeatedly mix it with 0’s, making sure to generate atmost one waste droplet at each step. More specifically, after z steps we will have j z dropletswith concentration 1 / z , where j z = d j/ γ − σ − z e . In step z , mix these j z droplets with j z j z droplets with concentration 1 / z +1 . We then either have j z +1 = 2 j z , inwhich case there is no waste, or j z +1 = 2 j z −
1, in which case one waste droplet 1 / z +1 isproduced. Overall, C produces at most γ − σ waste droplets. In this section we detail the construction of our converters. Let t s denote the concentrationat the s th recursive step in the RPR process. We can assume that t s ∈ [ , ], because thecase t s ∈ ( , ] is symmetric. Recall that for a t s in this range, in Step (rpr1) we will chosean appropriate interval S k , for some k ∈ { , , } . Let S k = [ l, r ] (that is, l = k · and r = l + ). For each such k and all i, j ≥ i : l, j : r )-converterthat we will denote C ki,j . Our main objective here is to design these converters so that theyproduce as little waste as possible — ideally none. ( i : , j : ) -Converters C i,j We start with the case k = 2, because in this case the construction is relatively simple. Weshow how to construct, for all i, j ≥
1, our ( i : , j : )-converter C i,j that produces at mostone droplet of waste. These converters are constructed via an iterative process. We first giveinitial converters C i,j , for some small values of i and j , by providing specific graphs. Allother converters are obtained from these initial converters by repeatedly coupling them withother mixing graphs that we refer to as extenders .Let J init = { ( i, j ) } i,j ∈{ , } . The initial converters C i,j are defined for the four index pairs( i, j ) ∈ J init . Figure 7 illustrates the initial converters C , , C , and two extenders X , X . . Coviello Gonzalez and M. Chrobak XX:9 Converter C , produces one waste droplet and converter C , does not produce any waste.Converter C , can be obtained from C , by designating one of the outputs as waste.Converter C , is defined as C , = X • C , , and produces one waste droplet of . (Thus C , is simply a disjoint union of C , and X with one output designated as waste.)
14 12 1214 14 1412 1214 12 121214 3412 1212 14 14 C X X C Figure 7
Initial converters and extenders for the case I = (cid:8) i : , j : (cid:9) . The construction of other converters C i,j is based on the following observation: Supposethat we already have constructed some C i,j . Then (i) X • C i,j is a C i,j +2 converter thatproduces the same waste as C i,j , and (ii) provided that j ≥ X • C i,j is a C i +2 ,j − converterthat produces the same waste as C i,j .Let now i, j ≥ i, j ) / ∈ J init be arbitrary. To construct C i,j , using the initialconverters and the above observation, express the integer vector ( i, j ) as ( i, j ) = ( i , j ) + φ (0 ,
2) + ψ (2 , − i , j ∈ J init and integers ψ = d i e − φ = d j + ψ e −
1. Then C i,j is constructed by starting with C i ,j and coupling it φ times with X and then ψ timeswith X . (This order of coupling is not unique but is also not arbitrary, because each extender X requires a droplet of concentration as input.) Since X and X do not produce waste, C i,j will produce at most one waste droplet. ( i : , j : ) -Converters C i,j Next, for each pair i, j ≥ i : , j : )-converter C i,j . These convertersare designed to produce one droplet of waste. ( C , will be an exception, see the discussionbelow). Our approach follows the scheme from Section 4.1: we start with some initialconverters, which then can be repeatedly coupled with appropriate extenders to produce allother converters. Since concentrations and are symmetric (as = 1 − ), we will onlyshow the construction of converters C i,j for i ≥ j ; the remaining converters can be computedusing symmetric mixing graphs.Let J init = { ( i, } i ∈{ , ,..., } ∪ { (2 , } . The initial converters C i,j are defined for all indexpairs ( i, j ) ∈ J init . Figure 8 shows converters C , , C , , . . . , C , and C , . Converter C , can be obtained from C , by designating an output of as waste. Converter C , is almostidentical to X in Figure 9; except that the source labels and are replaced by 0 and1, respectively (the result of mixing is still , so other concentrations in the graph are notaffected). Converters C , and C , are obtained from C , by designating outputs of (cid:8) , (cid:9) and , respectively, as waste. Note that all initial converters except for C , produce at mostone droplet of waste.Now, consider extenders X and X in Figure 9. The construction of other converters C i,j follows the next observation: Assume that we have already constructed some C i,j , with X:10 Mixing Graphs
14 34
12 583838 38 58
14 583838 38 58
14 583838 12 5838 C C C C C w w w w w C Figure 8
Initial converters for the case I = (cid:8) i : , j : (cid:9) . i ≥ j . Then (i) X • C i,j is a C i +1 ,j +1 converter that produces the same waste as C i,j , and(ii) X • C i,j is a C i +8 ,j converter that produces the same waste as C i,j . X X Figure 9 X and X extenders for the case I = (cid:8) i : , j : (cid:9) . Consider now arbitrary i ≥ j ≥ i, j ) / ∈ J init . To construct C i,j , using theinitial converters and the above observation, express the integer vector ( i, j ) as ( i, j ) =( i , j ) + φ (1 ,
1) + ψ (8 , φ, ψ ≥
0, and ( i , j ) ∈ J init − { (1 , } . Then C i,j is constructed by starting with C i ,j and coupling it φ times with X and then ψ times with X (in arbitrary order). Since X and X do not produce waste (and we do not use theinitial converter C , ), C i,j will produce at most one waste droplet.Overall, all converters C i,j , except for C , produce at most one waste droplet. Converter C , produces two droplets of waste; however, as we later show in Section 5, it is not actuallyused in the algorithm. ( i : , j : ) -Converters C i,j In this section, for each pair i, j ≥ i : , j : )-converter C i,j . Most ofthese converters produce at most one droplet of waste, but there will be four exceptionalcoverters with waste two. (See the comments at the end of this section.) The idea of the . Coviello Gonzalez and M. Chrobak XX:11 construction follows the same scheme as in Sections 4.1 and 4.2: we start with some initialconverters and repeatedly couple them with appropriate extenders to obtain other converters.
22 2 2
22 23 2 C C C w C C w w w w C Figure 10
Initial converters for the case I = (cid:8) i : , j : (cid:9) . Let J init = { ( i, j ) } i,j ∈{ , , } ∪ { (4 , , (2 , } . The initial converters C i,j are defined forall index pairs ( i, j ) ∈ J init . Converters C , , C , , C , , C , , C , and C , are shown inFigure 10. Converters C , , C , and C , are obtained from C , by designating outputs of (cid:8) , (cid:9) , and , respectively, as waste. Converter C , is obtained from C , by designatingan output of as waste, and C , is obtained from C , by designating an output of aswaste. Thus, among the initial converters, C , , C , and C , each produces two droplets ofwaste; all other converters have at most one droplet of waste.Next, we provide an observation leading to the construction of other converters C i,j .Consider extenders X and X in Figure 11 and assume that we have already constructedsome C i,j . Then, (i) provided that j ≥ X • C i,j is a C i +3 ,j − converter that producesthe same waste as C i,j , and (ii) provided that i ≥ X • C i,j is a C i − ,j +3 converter thatproduces the same waste as C i,j . We also need the following, less obvious observation: (cid:73) Observation 1. If i, j ≥ i, j ) / ∈ J init ∪{ (6 , } , then ( i, j ) = ( i , j )+ φ ( − , ψ (3 , − φ, ψ ≥
0, and ( i , j ) ∈ J init − { (1 , , (1 , , (3 , } . Proof.
Let i, j ≥ i, j ) / ∈ J init ∪ { (6 , } . We note first that we can represent ( i, j ) as( i, j ) = (˜ ı, ˜ ) + ˜ φ ( − ,
3) + ˜ ψ (3 , − ı, ˜ ) ∈ J init − { (2 , , (4 , } and integers ˜ φ, ˜ ψ ≥
0. If(˜ ı, ˜ ) / ∈ { (1 , , (1 , , (3 , } then we are done. Otherwise, we show how to modify the valuesof parameters ˜ ı , ˜ , ˜ φ and ˜ ψ so that they satisfy the condition in the observation.Case 1: (˜ ı, ˜ ) = (1 , φ, ˜ ψ ≥ i, j ≥
1. Therefore, we can write ( i, j ) as ( i, j ) = (3 ,
3) + ( ˜ φ − − ,
3) +( ˜ ψ − , − ı, ˜ ) = (1 , ψ ≥ i ≥
1. Therefore, we canwrite ( i, j ) as ( i, j ) = (4 ,
2) + ˜ φ ( − ,
3) + ( ˜ ψ − , − ı, ˜ ) = (3 , φ ≥
1, since we could thenwrite ( i, j ) as ( i, j ) = (2 ,
5) + ( ˜ φ − − ,
3) + ˜ ψ (3 , − φ ≥ φ = 0. Then ( i, j ) = (3 ,
2) + ˜ ψ (3 , − ψ ∈ { , } this contradicts that ( i, j ) / ∈ J init ∪ { (6 , } , and for ˜ ψ ≥ j ≥ (cid:74) X:12 Mixing Graphs X X Figure 11 X and X extenders for the case I = (cid:8) i : , j : (cid:9) . Using the observations above, for any pairs i, j ≥ C i,j asfollows. If ( i, j ) = (6 ,
1) we let C , = X • C , (so C , has two droplets of waste). If( i, j ) = (6 , C i,j by starting with C i ,j and repeatedly coupling it with φ copies of X and ψ copies of X , choosing a suitable order of couplings to ensure that eachintermediate converter has at least one output and at least one . (For example, if j = 1then we begin by coupling X first.) As X and X do not produce any waste, these C i,j ’swill each produce at most one droplet of waste.Overall, the converters C i,j we construct have at most one droplet of waste, with theexception of the following four: C , , C , , C , and C , . (It is easy to prove that for theseconverters waste 2 cannot be avoided.) As we show later in Section 5, of these four convertersonly C , is actually used in the RPR process of Algorithm RPRIS , and it is used at mostonce.
In this section we provide the analysis of Algorithm
RPRIS , including the worst-case boundon produced waste, a bound on the size of computed mixing graphs, and the running time.
Bound on waste . We first estimate the number of waste droplets of Algorithm
RPRIS . Let G be the mixing graph constructed by RPRIS for a target concentration t with its correspondingvalues d = prec ( t ) and γ (as defined in Section 1). Below we prove the following theorem. (cid:73) Theorem 1.
The number of waste droplets in G is at most ( d + γ ) + 2 . To prove Theorem 1, we show that the total number of sink nodes in G is at most ( d + γ − σ ) + 3, for corresponding σ ∈ { , } . (This is sufficient, as one sink node is used toproduce t ).Following the algorithm description in Section 3, let G = G • C . From our constructionof C (at the end of Section 3), we get that C contributes at most γ − σ sink nodes to G . (Each waste droplet produced by C represents a sink node in G .) Therefore, to proveTheorem 1 it remains to show that G contains at most ( d − γ + σ ) + 3 sink nodes. This isequivalent to showing that G , computed by process RPR for t (and used to compute G ),contains at most d + 3 sink nodes, where d = prec ( t ) = d − γ + σ . Lemma 2 next provesthis claim. (cid:73) Lemma 2.
The number of sink nodes in G is at most d + 3 . . Coviello Gonzalez and M. Chrobak XX:13 Proof.
Let t b be the concentration used for the base case of the RPR process and d b = prec ( t b ) ≤ d its precision. We prove the lemma in three steps. First, we show that (i) thenumber of sink nodes in the mixing graph computed for t b is at most three. (In particular,this gives us that the lemma holds if t = t b .) Then, we show that (ii) if t = t b then thenumber of converters used in the construction of G is no more than d −
1, and (iii) that atmost one of such converter contains two waste sink nodes. All sink nodes of G are either inits base-case graph or in its converters, so combining claims (i), (ii) and (iii) gives a completeproof for Lemma 2.The proof of (i) is by straightforward inspection. By definition of the base case, t b ∈B = (cid:8) , , , , , , (cid:9) . The mixing graphs for base concentrations are shown in Figure 6.(The graphs for , , and are symmetric to B , B , and B .) All these graphs have atmost 3 sink nodes.Next, we prove part (ii). In each step of the RPR process we reduce the precisionof the target concentration by 2 until we reach the base case, which gives us that thenumber of converters is exactly ( d − d b ). It is thus sufficient to show that d b ≥
2, as thisimmediately implies (ii). Indeed, the assumption that t = t b and the definition of the basecase implies that d ≥
4. (This is because the algorithm maintains the invariant that itstarget concentration is in [ , ] and all concentrations in this interval with precision at most3 are in B .) This, and the precision of the target concentration decreasing by exactly 2 ineach step of the recursion, imply that d b ∈ { , } holds.We now address part (iii). First we observe that converters C k , are not used in theconstruction of G : If we did use C k , in the construction of G then the source labels forthe next recursive step are { , } . Hence, t b = . Now, let t b − be the concentration, and S k = [ l, r ] the interval, used to compute t b . Since t b = , then t b − = ( l + r ). Therefore, bydefinition of S k , t b − ∈ (cid:8) , , , , (cid:9) ⊂ B , so Algorithm RPRIS would actually use a basecase mixing graph for t b − , instead of constructing C k , for t b .So, it is sufficient to consider C ki,j converters that satisfy i + j ≥ i, j ≥
1. Now,from Sections 4.1, 4.2 and 4.3, we observe that the only such converters that contain twowaste sink nodes are C , , C , and C , . Claim 1 below shows that converters C , and C , are not used in the construction of G .Regarding C , , first we note that this converter has exactly six source nodes; see Figure 10,Section 4.3. This implies that C , can not be used more than once in the constructionof G , since the number of source nodes at each recursive step in the RPR process isdecreasing. (Note that there are symmetric converters C , , C , and C , for C , , C , and C , , respectively, where superscript 5 is associated to interval S . Nevertheless, a similarargument holds.) Thus, step (iii) holds. (cid:73) Claim 1.
Converters C , and C , are not used by Algorithm RPRIS in the construction of G for t .We first present the following observations. Consider recursive step s of the RPR process,for which t s is the target concentration. If a converter C i,j is used in this step, then t s ∈ ( , ]must hold; that is t s is in the middle part of interval S (see Figure 4 in Section 3). (Recallthat, by our algorithm’s invariant, t s ∈ [ , ]. Also, note that t s = since otherwise thiswould be a base case and the algorithm would use B from Figure 6 instead.) Further, atthe next step of the RPR process, t s +1 = 4( t s − ) satisfies t s +1 ∈ ( , ].We now prove the claim by contradiction, using the above observations. Assume thateither C , or C , were used in the construction of G . If C , was used in the construction of G , then the concentration labels of the source nodes at the next recursive step are { , } ,and thus, since t s +1 > , there is not enough reactant available to produce t s +1 . X:14 Mixing Graphs
On the other hand, if C , was used in the construction of G , then the concentrationlabels of the source nodes at the next recursive step are { , } . This implies thatthe next step is guaranteed not to be a base case, since all mixing graphs used for basecase concentrations contain at most three source nodes, as illustrated in Figure 6. Now, as t s +1 > , depending on the exact value of t s +1 , the chosen interval for t s +1 must be either S = [ , ], S = [ , ] or S = [ , ]. We now consider these three cases.Case 1: t s +1 ∈ ( , ]. Then the chosen interval is S = [ , ]. The only C i,j converter withsource concentration labels { , } is C , (see in Figure 8 in Section 4.2), whose sinknodes have concentration labels (cid:8) , , (cid:9) . Therefore, the input configuration for the nextrecursive step will be a subset of { , } , which does not have enough reactant to produce4( t s +1 − ) > , thus contradicting the choice of S .Case 2: t s +1 ∈ ( , ]. Then the chosen interval is S = [ , ]. This instance is symmetricto interval S , having source concentration labels { , } , instead of { , } , andtarget concentration t s +1 = (1 − t s +1 ). Thus we proceed accordingly. Since every converterand extender in Section 4.1 adds at least the same number of source nodes with concentrationlabel 0 as source nodes with concentration label 1, then no converter constructed by thealgorithm will have source concentration labels { , } . Hence, we have a contradictionwith the choice of S for t s +1 , and thus also with the choice of S for t s +1 .Case 3: t s +1 ∈ ( , ]. Then the chosen interval is S = [ , ]. The argument here issimple: to produce concentration , at least three reactant droplets are needed, but the inputconfiguration contains only two. Therefore, at the next recursive step, the algorithm will nothave enough reactant droplets to construct a converter C i,j with i, j ≥
1, contradicting thechoice of S for t s +1 .Finally, neither S , S nor S are chosen by our algorithm for t s +1 , contradicting C , being used for the construction of G .This completes the proof of Claim 1 and Lemma 2 (thus also completing the proof ofTheorem 1). (cid:74) Size of mixing graphs and running time . Let G = G • C be the mixing graphcomputed by Algorithm RPRIS for t ; C is constructed by process IS while G is obtainedfrom G (constructed by process RPR) by changing concentration labels appropriately. Weclaim that the running time of Algorithm RPRIS is O ( | G | ), and that the size of G is O ( d ),for d = prec ( t ). We give bounds for G and C individually, then we combine them to obtainthe claimed bounds. (This is sufficient because the size of G , as well as the running time toconstruct it, is asymptotically the same as that for G .)First, following the description of process RPR in Section 3, suppose that at recursivestep s , G s +1 , G s +1 and converter C s +1 = C ki,j are computed. (Note that the algorithm doesnot need to explicitly relabel G s +1 to get G s +1 – we only distinguish G s +1 from G s +1 for thepurpose of presentation.) The size of C ki,j is O ( i + j ) and it takes time O ( i + j ) to assembleit (as the number of required extenders is O ( i + j )). Coupling C s +1 with G s +1 also takestime O ( i + j ), since I s +1 (the input configuration for G s +1 ) has cardinality O ( i + j ) as well.In other words, the running time of each recursive RPR step is proportional to the numberof added nodes. Thus the overall running time to construct G is O ( | G | ).Now, let t be the target concentration for the RPR process, with d = prec ( t ). Then,the size of G is O ( d ). This is because the depth of recursion in the RPR process is O ( d ),and each converter used in this process has size O ( d ) as well. The reason for this bound onthe converter size is that, from a level of recursion to the next, the number of source nodes . Coviello Gonzalez and M. Chrobak XX:15 increases by at most one (with an exception of at most one step, as explained earlier in thissection), and the size of a converter C ki,j used at this level is asymptotically the same as thenumber of source nodes at this level. ( I s and I s +1 in Figure 5 illustrate the idea.)Regarding the bounds for C , we first argue that the running time to construct C is O ( | C | ). This follows from the construction given in Section 3; in step s there are 2 j s dropletsbeing mixed, which requires j s nodes; thus the entire step takes time O ( j s ).We next show that the size of C is O ( d ). Let I be the input configuration for G .From the analysis for G , we get that | I | = O ( d ), so the last step in C contains O ( d )nodes. Therefore, as the depth of C is γ − σ , the size of C is O ( γd ) = O ( d ).Combining the bounds from G and C , we get that the running time of Algorithm RPRIS is O ( | G | ) and the size of G is O ( d ). (The coupling of C with G does not affect the overallrunning time, since it takes O ( d ) time to couple them, as | I | = O ( d ).) In this section we compare the performance of our algorithm with algorithms
Min-Mix , REMIA , DMRW , GORMA and
ILP . We start with brief descriptions of these algorithms, to givethe reader some intuitions behind different approaches for constructing mixing graphs. Let t ∈ (0 ,
1) be the target concentration and d = prec ( t ) its precision. Also, let bin ( t ) be t ’sbinary representation with no trailing zeros. Min-Mix [10]:
This algorithm is very simple. It starts with τ = 0 and mixes it with thebits of bin ( t ) in reverse order, ending with τ = t . It runs in time O ( d ) and produces d droplets of waste. REMIA [6]:
This algorithm is based on two phases. In the first phase, the algorithm computesa mixing graph G whose source nodes have concentration labels that have exactly onebit 1 in their binary representation; each such concentration represents each of the 1 bitsin bin ( t ). Then, in the second phase, a mixing graph G (that minimizes reactant usage),whose sink nodes are basically a superset of the source nodes in G , is computed. Finally, G for t is obtained as G • G . (Although REMIA targets reactant usage, its comparisonto different algorithms in terms of total waste was also reported in [6]. Thus, for the sakeof completeness, we included
REMIA in our study.)
DMRW [9]:
This algorithm is based on binary search. Starting with pivot values l = 0 and r = 1, the algorithm repeatedly “mixes” l and r and resets one of them to their average ( l + r ), maintaining the invariant that t ∈ [ l, r ]. After d steps we end up with l = r = t .Then the algorithm gradually backtracks to determine, for each intermediate pivot value,how many times this value was used in mixing, and based on this information it computesthe required number of droplets. This information is then converted into a mixing graph. GORMA [3]:
This algorithm enumerates the mixing graphs for a given target concentration.An initial mixing graph is constructed in a top-down manner; starting from the targetconcentration t (the root node), the algorithm computes two concentrations x and y (called a preceding pair) such that t = ( x + y ) and both x and y have smaller precisionthan t ; x and y become t ’s children and both x and y are then processed recursively.(Note that a concentration might have many distinct preceding pairs. Each precedingpair is processed.) A droplet sharing process is then applied to every enumerated mixinggraph to decrease reactant usage and waste produced. A branch-and-bound approach isadopted to ease its exponential running time. ILP [4]:
This algorithm constructs a “universal” mixing graph that contains all mixinggraphs of depth d as subgraphs. It then formulates the problem of computing a mixing X:16 Mixing Graphs graph minimizing waste as an integer linear program (a restricted flow problem), andsolves this program. This universal graph has size exponential in d , and thus the overallrunning time is doubly exponential in d .We now present the results of our experiments. Each experiment consisted on generatingall concentration values with precision d , for d ∈ { , , , } , and comparing the outputs ofeach of the algorithms. The results for GORMA and
ILP are shown only for d ∈ { , } , sincefor d ∈ { , } the running time of both GORMA and
ILP is prohibitive.
Concentration W a s t e Concentrations with precision 7
MinMixREMIADMRWGORMARPRISILP
Concentration W a s t e Concentrations with precision 8
MinMixREMIADMRWGORMARPRISILP
Figure 12
The number of waste droplets of algorithms
Min-Mix , REMIA , DMRW , GORMA , ILP , andour algorithm
RPRIS , for all concentrations with precision 7 (top figure) and 8 (bottom figure). Allgraphs are smoothed using
MATLAB ’s smooth function. Figure 12 illustrates the experiments for concentrations of precision 7 and 8. Figure 13illustrates the experiments for concentrations of precision 15 and 20. In both figures, thedata was smoothed using
MATLAB ’s smooth function to reduce clutter and to bring out thedifferences in performance between different algorithms.As can be seen from these graphs, RPRIS significantly outperforms algorithm
Min-Mix , REMIA , DMRW and
GORMA :It produces on average about 50% less waste than
Min-Mix (consistently with our boundof ( d + γ )+4 on waste produced by RPRIS ), and 40% less waste than
REMIA . It also produceson average between 21 and 25% less waste than
DMRW , with this percentage increasing with d .Additionally, for d = 7 , RPRIS produces on average about 17% less waste than
GORMA and . Coviello Gonzalez and M. Chrobak XX:17 only about 7% additional waste than
ILP . Figure 13
The number of waste droplets of algorithms
Min-Mix , DMRW , REMIA , and our algorithm
RPRIS , for all concentrations with precision 15 (top figure) and 20 (bottom figure). All graphs aresmoothed using
MATLAB ’s smooth function. Among all of the target concentration values used in our experiments, there is not asingle case where
RPRIS is worse than either
Min-Mix or REMIA . When compared to
DMRW , RPRIS never produces more waste for precision 7 and 8. For precision 15, the percentage ofconcentrations where
RPRIS produces more waste than
DMRW is below 2%, and for precision20 it is below 3 . GORMA , the percentage of concentrationswhere
RPRIS produces more waste is below 4%.
In this paper we proposed Algorithm
RPRIS for single-droplet targets, and we showed thatit outperforms standard waste minimization algorithms
Min-Mix and
DMRW in experimentalcomparison. We also proved that its worst-case bound on waste is also significantly betterthan for the other two algorithms.Many questions about mixing graphs remain open. We suspect that our bound onwaste can be significantly improved. It is not clear whether waste linear in d is needed forconcentrations not too close to 0 or 1, say in [ , ]. In fact, we are not aware of even a super-constant (in terms of d ) lower bound on waste for concentrations in this range. X:18 Mixing Graphs
For single-droplet targets it is not known whether minimum-waste mixing graphs canbe effectively computed. The most fascinating open question, in our view, is whether it isdecidable to determine if a given multiple-droplet target set can be produced without anywaste. (As mentioned in Section 1, the ILP-based algorithm from [4] does not always producean optimum solution.)Another interesting problem is about designing mixing graphs for producing multipledroplets of the same concentration. Using perfect-mixing graphs from [5], it is not difficultto prove that if the number of droplets exceeds a certain threshold then such target sets canbe produced with at most one waste droplet. However, this threshold value is very large andthe resulting algorithm very complicated. As such target sets are of practical significance, asimple algorithm with good performance would be of interest.It would also be interesting to extend our proposed worst-case performance measure toreactant minimization. It is quite possible that our general approach of recursive precisionreduction could be adapted to this problem.
References Sukanta Bhattacharjee, Bhargab B. Bhattacharya, and Krishnendu Chakrabarty.
Al- gorithms for Sample Preparation with Microfluidic Lab-on-Chip . River Publishers, 2019. Bhargab B. Bhattacharya, Sudip Roy, and Sukanta Bhattacharjee. Algorithmic challengesin digital microfluidic biochips: Protocols, design, and test. In
Proc. International Confer-ence on Applied Algorithms (ICAA’14) , pages 1–16, 2014. Ting-Wei Chiang, Chia-Hung Liu, and Juinn-Dar Huang. Graph-based optimal reactantminimization for sample preparation on digital microfluidic biochips. In , pages 1–4. IEEE, 2013. Trung Anh Dinh, Shinji Yamashita, and Tsung-Yi Ho. A network-flow-based optimalsample preparation algorithm for digital microfluidic biochips. In , pages 225–230. IEEE, 2014. Miguel Coviello Gonzalez and Marek Chrobak. Towards a theory of mixing graphs: a charac-terization of perfect mixability. In
International Conference on Algorithms and Complexity ,pages 187–198. Springer, 2019. Juinn-Dar Huang, Chia-Hung Liu, and Ting-Wei Chiang. Reactant minimization duringsample preparation on digital microfluidic biochips using skewed mixing trees. In
Pro-ceedings of the International Conference on Computer-Aided Design , pages 377–383. ACM,2012. Juinn-Dar Huang, Chia-Hung Liu, and Huei-Shan Lin. Reactant and waste minimizationin multitarget sample preparation on digital microfluidic biochips.
IEEE Transactions onComputer-Aided Design of Integrated Circuits and Systems , 32(10):1484–1494, 2013. Debasis Mitra, Sandip Roy, Krishnendu Chakrabarty, and Bhargab B Bhattacharya. On- chip sample preparation with multiple dilutions using digital microfluidics. In
IEEE Com-puter Society Annual Symposium on VLSI (ISVLSI) , pages 314–319. IEEE, 2012. Sandip Roy, Bhargab B Bhattacharya, and Krishnendu Chakrabarty. Optimization of di-lution and mixing of biochemical samples using digital microfluidic biochips.
IEEE Trans-actions on Computer-Aided Design of Integrated Circuits and Systems , 29(11):1696–1708,2010. William Thies, John Paul Urbanski, Todd Thorsen, and Saman Amarasinghe. Abstractionlayers for scalable microfluidic biocomputing.