Shape-Dependent Velocity Based Droplet Routing on MEDA Biochips

Digital microfluidic biochip (DMFB) has attracted attention in the biochemical and medical industries. In particular, a microelectrode dot array (MEDA) biochip, which is composed of a two-dimensional microelectrode array, enables to realize fine-grained manipulation such as dilution, mixing, sensing, and so on in real-time. Unlike existing DMFB biochips, a MEDA architecture allows microelectrodes to control a certain volume of droplet in a fine-grained manner and can vary droplet volume and shape in such a way that it efficiently conducts synthesis and manipulation of droplets. There have been many works in order to improve the efficiency of synthesis of MEDA biochips; however, the synthesis, especially droplet routing, has never considered the shape-dependent velocity of droplets. In this paper, we propose the droplet routing techniques for MEDA biochips with shape-dependent velocity of droplets. The proposed techniques take the advantage of variant velocities of droplets dependent on the shapes and aim to reduce the overall routing time of a droplet from a source to a destination. Simulation results confirm that the proposed techniques can shorten the routing time by 80% compared to the state-of-the-art techniques.


I. INTRODUCTION
Digital microfluidic biochips (DMFBs) have been prevalent as one of Labs-on-a-Chip (LoCs), in a variety of biochemical and pharmacy fields [1], [2] [3]. A DMFB manipulates micro/pico-liter of droplets in biochemical samples on a two-dimensional electrode array based on the electrowetting on dielectric (EWOD) principle [4]. However, existing DMFBs contain several limitations [5], [6]: DMFBs are not capable of controlling the volume and shape of droplets during manipulations and of detecting droplets in real-time, and thus they encounter functionality issues in practical use.
DMFBs based on microelectrode dot array (MEDA) architecture, so-called MEDA biochips, have been developed to The associate editor coordinating the review of this manuscript and approving it for publication was Liang-Bi Chen . overcome the limitations. Unlike the conventional DMFBs, MEDA biochips compose a myriad of microelectrode cells (MCs), where the concept is built on a sea-of-microelectrode array [7], that are fabricated by TSMC 0.35µm CMOS technology [8]. In addition, in contrast to DMFBs, MEDA biochips can dynamically control the size and shape of droplets and detect droplets in real-time during manipulations, resulting in that droplet routing is no longer negligible.
A certain volume of droplets on a MEDA biochip can have a variety of shapes by controlling the MCs. The actuation force on a droplet from a group of MCs works on the droplet to transfer, and the moving velocity is dependent on its volume and shape since the actuation force works to adjacent area of the droplet [9]. However, most of the previous works on MEDA biochips have not assumed the velocity of droplets dependent on the shape, resulting in the inefficient routing in terms of low latency. In addition, it should also be taken into account the reshaping overhead of an intermediate droplet during reshaping [10], [11].
In this paper, we propose droplet routing techniques with shape-dependent velocity. We first formulate the proposed routing techniques that aim to find the best route from a source to destination cell with taking advantage of shape-dependent velocity of a droplet in such a way that the routing time is minimized. Then we propose two routing techniques to efficiently explore feasible solutions. The contributions of this work are presented as follows: 1) We achieve reduction of the routing time by taking advantage of shape-dependent velocity on a MEDA biochip. 2) We propose two efficient routing techniques that solve the proposed routing problem.
The rest of this paper is organized as follows. Section II describes related work. Section III formulates our droplet routing problem that takes into account shape-dependent velocity. Section IV proposes two routing techniques to efficiently solve the problem presented in Section III. Section V describes the simulations and comparison, and Section VI concludes this paper.

II. RELATED WORK
Much work for DMFBs has been studied since 2000s. The authors in [6] mentioned that routing complexity increased as an increase in the complexity and number of bioassay manipulations mapped on a DMFB, unlike the conventional microfluidic biochips. Cho et al. envisioned that a DMFB design can be done as in VLSI techniques. They provided a three-dimensional representation to schedule manipulations, which describes when and where a manipulation is performed [12]. The goal of this work was to develop a high-performance droplet router with combining conventional high-level synthesis techniques in VLSI. The work in [13] proposed a Tabu search metaheuristic for the synthesis of a DMFB. The literature has been aware of VLSI design techniques, which consist of the allocation, resource binding, scheduling, and placement of the manipulations in the applications. However on DMFBs, there have been several functional limitations: One is that the (1:1) mixing model has basically been employed for mixing and dilution, which may take the long completion time after all [14]. In addition, droplets on the DMFB move along adjacent electrodes and are limited to movement in the x-y axis direction only [15]. Thus, the DMFBs have suffered from the limitations.
A MEDA biochip is smaller than the electrodes used in DMFBs and has MCs with control and detection units on the chip, which is dynamically grouped to efficiently achieve a variety of droplet manipulations [8], [16] [17]. Due to the efficient utilization of MCs, MEDA biochips, for example, can detect a droplet at any position within 10 milliseconds [8]. An error recovery method based on this feature has also existed [18]. Moreover, the MEDA biochips enable diagonal movement by appropriately grouping MCs, and this degree of freedom can be exploited for more efficient manipulations [19], [20]. By taking advantage of the features in MEDA biochips, the manipulations such as mixing and dilution have been improved and the time for each manipulation has been significantly reduced [21]. MEDA biochips have also been successfully used to reduce the size of biochips by using the split operation [22]. In contrast, routing of droplets has been remarkably important than ever, and it needs an efficient technique [11].
There have been studies on MEDA biochips that address droplet routing problems [23], [24] [25], [26] [27], [28]. The authors in [23] assumed that the droplet aspect ratio can be changed only to avoid obstacles on the chip. This work considers the droplet velocity to be dependent on its volume, which controls the thickness of the droplet. In particular, one of challenging issues is that a MEDA biochip handles a larger volume of droplets, the velocity deterioration occurs, resulting in much time for droplet routing [7]. Unfortunately, there is no work to improve the efficiency especially in terms of the routing time by focusing on the velocity dependent on droplet shape.
Some studies have shown that the EWOD force on a droplet varies depending on the droplet volume, shape, and moving direction [7], [29] [30]. In [7], the authors proposed 2D layout design with triangle shaped MCs. In their MEDA biochips, they conducted simulations with considering cross-contamination aware routing for MEDA biochips. A MEDA biochip contains a plenty of microelectrodes, but several cells often become unavailable due to out of order, degradation of EWOD and contamination on the cells. Furthermore, the degradation of EWOD occurs at a certain rate, according to the work in [10]. Under the circumstances, droplet routing has its difficulty with finding an efficient route with satisfying timing requirements for the synthesis. Therefore, this paper attempts to transfer a droplet with a shape-dependent velocity under unavailable cells such as degradation and contamination.

A. PROBLEM DESCRIPTION
In this section, we formulate routing of a single droplet that is transferred from a source cell to destination cell under unavailable some MCs on a MEDA biochip. We are given a reconfigurable W × H MEDA biochip, where the coordinates of each microelectrode are defined as (x, y) with the origin (1, 1) at the lower-left cell. MCs on a MEDA biochip can control a certain volume of a droplet. The droplet is allowed to reshape into a different shape at any time, and the velocity of the droplet is assumed to be depended on its volume and shape. According to the work in [21], the velocity of the droplet depends on the number of actuated MCs adjacent to the droplet.
To introduce the shapes of a droplet, we consider a droplet (A × B) on a MEDA biochip in Figure 1 (a). The horizontal motion requires B of the actuated MCs as shown in   Figure 1 (c). In Figure 1 (d), when the droplet moves diagonally, the number of active MCs is assumed to be (A+B−1). Here, it should be noted that motion in the diagonal direction takes longer distance than either of the horizontal and vertical motion. On MEDA biochips, a droplet can reshape into a different shape.
, the number of actuated MCs is A (Figure 1 (e)). Otherwise, the number of actuated MCs is B (Figure 1 (f)).
Next, we describe the time required for each motion. Basically, the actuation force from MCs can be modeled by profiling the droplet energy as a function of droplet location [21]. For simplicity, we assume that the velocity of a droplet is proportional to the number of actuated MCs adjacent to its droplet but is inversely proportional to the volume, except for diagonal motion. In the real world, such parameters are profiled for each type of droplets such that we should consider chemical and physical properties of droplets. For example, if a droplet of A × B moves horizontally by actuating B MCs, then it takes A time steps. Similarly, the vertical motion needs B time steps and the reshape motion from (A × B) to (B × A) actuates B MCs, taking A time steps. The diagonal motion is specially assumed to be a combination of the horizontal and vertical motion; thus, it takes A + B time steps with actuating A + B − 1 MCs.
A MEDA biochip often contains unavailable MCs due to degradation and contamination. The locations of unavailable MCs are analyzed and known in advance. Given a biochip size, droplet with motions, and the locations of unavailable MCs, our routing tries to find the best route by reshaping under unavailable MCs such that the overall routing time from a source to destination is minimized.

B. EXAMPLE
We describe a motivating example that considers a routing problem with a droplet of size 2 in Figure 2. Figure 2 (a) represents the 6 × 6 MEDA biochip. In the problem, we assume that our routing of a droplet starts from a source to destination as shown in the two cells colored in green. As the initial  state, there are two shapes such as horizontal (2 × 1) and vertical (1 × 2) shapes at a source, and this example assumes the former shape. Four cells painted in black in Figure 2 are represented as unavailable MCs due to degradation, contamination, and so on. This problem aims at minimization of the total routing time.
Before presenting our proposed routing, we first derive the overall routing time based on the state-of-the-art (SOTA) technique in [21]. This work allows a droplet to reshape only when it cannot find any route to the destination without reshaping. The droplet moves four cells to the right in Figure 2 (a). Here, the time step required for the transition is 8 (=2 × 4) time steps. Then, the droplet in the shape moves three cells to the upper, as shown in Figure 2 (b). Similarly, the transition takes 3 (= 1 × 3) time steps. In Figure 2 (c), the droplet can hardly reach the goal without reshaping due to the unavailable cells. The droplet is reshaped into the vertical one to avoid unavailable cells. Then, the time in this transition to reshape the droplet is 2 (= 2 × 1). When the droplet moves one cell to the upper as shown in Figure 2 (d), it takes 2 (= 2 × 1) time steps. The overall routing time is 15 (= 8 + 3 + 2 + 2). Recall that the state-of-the-art technique allows the droplet to reshape only for avoiding the unavailable cells.
On the other hand, the proposed routing technique allows the droplet to reshape at any time. Figure 3 (a) shows that the droplet, which has moved one cell to the right, reshapes into the vertical one to take advantage of the force from MCs, resulting in a higher velocity of the shape than the horizontal one. So far, it takes 4 (= 2 × 1 + 2 × 1) time steps. Then, the droplet further moves two cells to the right and reshapes into the horizontal one again in Figure 3 (b), which takes 4 (= 1 × 2 + 2 × 1) time steps. As shown in Figure 3 (c), the horizontal droplet moves two cells to the upper and reshapes into the vertical one, which takes 4 (= 1 × 2 + 2 × 1) time steps. Finally, the vertical droplet reaches the goal by moving one cell, which takes 2 (= 2 × 1) time steps, as shown in Figure 3 (d). Here, the overall routing time is 14 (= 4 + 4 + 4 + 2) time steps, achieving the reduction of the routing time compared to the state-of-the-art. The example demonstrates the potential of reshaping the droplet during routing with taking into account shape-dependent velocity.

C. FORMULATION
This section presents formulation for the droplet routing problem based on integer programming. Although the formulation in this paper is not described as linearized one due to the page limitation, we describe equivalent formulation with logical expressions. We presuppose the notations in the formulation as shown in Table 1.
We express the coordinate and shape of a droplet. Let s denote operation step and (x s , y s ) define the reference point of the droplet at step s. Operation step represents the number of droplet operations. Operation step 1 does not always represent time step 1. Here, the lower-left corner of the droplet is the coordinate called a reference point. The shape of the droplet can be expressed with using (w s , h s ), where w s and h s represent the width and height of the droplet, respectively. Thus, its droplet is located from (x s , y s ) to (x s + w s , y s + h s ). It should be noted a droplet is assumed to occupy the cells formed as a rectangle.
In Formula (1), the volume of a droplet is consistently fixed during its routing. For example, there are either of two states when the volume is 2: (w s × h s ) = (1 × 2) or (2 × 1).
Formula (2) shows the initial location of a droplet. Let X .source and Y .source represent a coordinate the of initial state at step 0.
Next, we express a number of droplet motions. In this paper, we assume that a droplet is allowed to move to several directions or reshape during routing. A droplet can be moved vertical, horizontal, and diagonal directions. In Figure 4 (a), dir s = 1 shows that the droplet horizontally (i.e., to the right or left) moves one cell. As shown in Figure 4 (b), dir s = 2 represents the vertical motion (i.e., the upper of lower). Figure 4 (c) expresses reshaping of the droplet. To move diagonally, dir s = 4 is determined as shown in Figure 4 (d).
Formula (3) shows the horizontal motion when dir s = 1. As shown in Figure 4 (a), the reference point moves in the horizontal direction but maintains in the vertical direction: Similar to Formula (3), Formula (4) shows the vertical motion of dir s = 2: Formula (5) expresses reshaping of the droplet of dir s = 3. When the droplet is reshaped, the shape at step s is different from the one at step s − 1: Formula (6) also shows the case of dir s = 3. Figure 5 shows four ways of reshaping the droplet if it reshapes from (1 × 2) to (2 × 1). As shown in the figure, the reference point can be different, depending on which way to reshape. The expression that allows such possible shapes is given by: Formula (7) shows the diagonal motion of the droplet as shown in Figure 4 Formula (9) expresses the constraints on unavailable cells. Droplets must not enter unavailable MCs.
In Formula (10), sum.dir i represents the number of times each motions.  In this paper, we assume that the time step for each motion depends on the volume of a droplet and number of actuated MCs that force on the droplet. The number of actuated MCs is represented as adjacent cells to the droplet for the direction. active w s ,h s ,motion means the number of active cells adjacent to the droplet. When a droplet of (A × B) moves horizontally, the number of active cells is given by active A,B,hor . Similarly, when the droplet moves vertically, it is given by active A,B,ver , and when the droplet is reshape, it is given by active A,B,res . We can derive the total routing time from a source to a destination cells as shown in Formula (11).
The objective aims at minimization of the total routing time: Minimize : rou.time (12)

IV. SIMULTANEOUS ROUTING WITH AN GUIDE DROPLET
Our routing problem presented in the previous section is classified as NP-hard problem, and it can hardly be solved in a practical time [31]. In order to efficiently explore solutions, we propose two routing techniques. As a basic concept of the techniques, we hypothesize that a good route for an unit of droplet is one of efficient routes for routing of a larger size of droplet to result in a shortest path problem; however, there exist multiple best routes for the unit of droplet. Thus, we present simultaneous routing framework utilizing both an unit of virtual droplet (called guide droplet) which is not exist on the biochip and the droplet (called target droplet) which is exist on the biochip. Along with the an guide droplet, the target droplet asks a best route in such a way that the total routing time is minimized. The route should be determined in two steps, such that the route of the guide droplet is determined in advance and then the route of the target droplet is determined. However, the best route of the target droplet does not always match with the best route of the guide droplet, so the method of two steps could not realize. The route of the guide droplet should be determined so that the routing time of the target droplet is minimized. Our technique efficiently finds the solution by prioritizing route determination over reshaping. The details of our techniques are described in the following sections.

A. ROUTE-GUIDED TECHNIQUE
The route-guided (RG) technique indicates that the guide droplet finds good one of routes and the target droplet is guided along the route with changing the shape. Note that this is not a warm start approach but simultaneously solved in the framework. More intuitively, this technique is similar to a route-constrained problem, but our technique does not eliminate the solution space. We are given an example as shown in Figure 6. The guide droplet finds a good route in a practical time as shown in Figure 6 (a), and the target droplet follows the route as shown in Figure 6 (b). The implementation of this technique can be expressed by adding the following formulae. The size of an guide droplet is assumed to be one: (1 × 1). The coordinates of this droplet is represented as (accmp.x s , accmp.y s ). Formula (13) indicates that the guide droplet can only move to adjacent coordinates.
accmp.map (x, y) becomes 1 when the guide droplet passes through the coordinate (x, y), otherwise zero. Formula (14) represents the cells where the guide droplet passes.
In Formula (15), the constraint expresses that the target droplet must be followed along the route where the guide droplet passes. The route-and time-guided (RTG) technique is a further extension of the RG technique presented in the previous section. This technique designates the location at each time step by an guide droplet. Thus, the target droplet must move or reshape with keeping the guidance by the guide droplet. Figure 7 shows the cell through which the guide droplet routes and the operation at each step are determined. A part of the target droplet moves to include the cell through which the guide droplet passed (Figure 3). The RTG technique can be implemented by adding the following formulae to the formulation in SectionIII-C. The RTG technique must guarantee the continuity of the guide droplet.
Formula (16) represents the positional relationship between the sub-droplet and target droplet in each step.

V. SIMULATIONS A. SIMULATION SETUP
We have conducted simulations to demonstrate the effectiveness of our proposed techniques. We were given a biochip that includes its size, unavailable MCs, a source and destination for routing. In addition, we assumed that the size of a target droplet is two given in the simulations. The goal of the routing problem is to minimize the overall routing time from a source to destination. We compared the routing time with the following two techniques: • SOTA: It allows to reshape only when a droplet is faced to blockage by unavailable MCs, which is presented in [21].
• Proposal: Our proposed routing technique that a droplet can reshape at any time as mentioned in Section III.
• SOTA-based RG: The SOTA technique based on the route-guided technique proposed in Section IV-A.
• SOTA-based RTG: The SOTA technique based on the route-and time-guided technique.
• Proposal-based RG: The proposed technique based on the route-guided technique.
• Proposal-based RTG: The proposed technique based on the route-and time-guided technique. There have appeared very little work considering a shape-dependent velocity model on droplet routing, and the work in [21] is a state-of-the-art technique. Although we could easily compare to the techniques which do not take advantage of shape-dependent velocity, it should be intuitively obvious that the results of routing time will take longer than those considering shape-dependent velocity. The techniques, which have a velocity model with or without taking into account shape dependency, also have smaller solution space than our proposed technique. Therefore, we did not employ others than the SOTA in [21] for comparison in this paper.
We show the simulation scenario as follows: We were given MEDA biochip assumed to be width = 10 and height = 15. Unavailable MCs accounted for all the cells were assumed to be 0%, 5%, 10%, and 20%, distributed randomly on MCs. For each failure rate of unavailable MCs except for 0%, we were given randomly 10 different arrangements of unavailable MCs. The failure rate is determined as one of the problem settings. Unavailable cells should be assumed for MEDA biochips because they can occur due to contamination and electrode failure. The target route was from (1, 1) to (10,15) as the source to destination.
The simulations have been conducted on Ryzen Threadripper 3970X (3.7 GHz, 32 cores and 64 threads in total) with 256 GB memory. IBM ILOG CPLEX Optimization Studio 20.1.0 is used to solve both the state-of-the-art and the proposed texhniques. The runtime is limited by up to 10 hours in CPU time. If an optimal solution is failed to be found within the runtime, the solution at the time limit is employed for comparison.

B. RESULTS
We show the simulation results for each failure rate of unavailable MCs as shown in Figures 8 (a), 8 (b) and 8 (c). The horizontal axis shows the problem instances that are randomly generated and the vertical axis shows the total routing time from the source to destination. Note that the results shown as Baseline on the left-most in each graph represent the failure rate is 0%, which means that all the MCs are available.
The overall results show that the proposed techniques achieves shorter routing time on average than the SOTAbased techniques. There was no significant difference in the computational time spent to determine the routing time. All results required nearly 10 hours in CPU time. As shown in Figure 8 (a), the proposed techniques can find good routes by taking advantage of reshaping during routing. In the 5th instance, the proposal-based RG obtains the longer routing time than the SOTA-based one because the solution space of our techniques is relatively larger than the SOTA-based ones. If the solution space is very large, it takes a great deal of time to determine the solution. Therefore, processing to the upper limit of the computation time worse the solution. Figure 8 (b) shows the similar tendency towards Figure 8 (a). The results present that the proposed technique shorten the routing time 15% on average. In Figure 8 (c), there are some cases that the proposed techniques obtain worse solutions compared to the SOTA-based ones. As increasing failure rate of MCs, the solution space shrinks and the best solutions of the proposed techniques become closee to the SOTA technique. Furthermore, the increased number of unavailable MCs enables the SOTA technique to have chances to reshape a droplet. In such a case, the effectiveness of the proposed techniques can hardly be observed if the unavailable cells are increased; however, the cases with 20% of failure rate are extremely special since such a biochip is generally discarded. Therefore, our proposed techniques can shorten the routing time in the practical failure rates.
We evaluate the RG-based and RTG-based methods. At 0%, 5%, and 10% (Figure 8 (a),8 (b)), the RTG-based method succeeded in reducing routing time by 10% to 20% compared to the SOTA method. On the other hand, at 20% (Figure 8 (c)), the RG-based method succeeded in reducing routing time by about 10% compared to the SOTA method. In all simulations, either RG-or RTG-based methods performed equally well or better than SOTA and Proposal-only methods. These results indicate that the RG-and RTG-based methods are sufficiently useful.

VI. CONCLUSION
In this paper, we have proposed droplet routing techniques with taking advantage of shape-dependent velocity to minimize the total routing time on a MEDA biochip. In addition, we have also presented two efficient exploring techniques called RG and RTG. Our simulation results have shown that our techniques significantly shorten the overall routing time 15% on on average, compared to the state-of-the-arts. In future, we will extend this work so that multiple droplets can be controlled with a combination of mixing, dilution, and so on at the same time. The evaluation of its scalability of this work is also left as our future work.