Simulation before fabrication: a case study on the utilization of simulators for the design of droplet microfluidic networks

The functional performance of passively operated droplet microfluidics is sensitive with respect to the dimensions of the channel network, the fabrication precision as well as the applied pressure because the entire network is coupled together. Especially, the local and global hydrodynamic resistance changes caused by droplets make the task to develop a robust microfluidic design challenging as plenty of interdependencies which all affect the intended behavior have to be considered by the designer. After the design, its functionality is usually validated by fabricating a prototype and testing it with physical experiments. In case that the functionality is not implemented as desired, the designer has to go back, revise the design, and repeat the fabrication as well as experiments. This current design process based on multiple iterations of refining and testing the design produces high costs (financially as well as in terms of time). In this work, we show how a significant amount of those costs can be avoided when applying simulation before fabrication. To this end, we demonstrate how simulations on the 1D circuit analysis model can help in the design process by means of a case study. Therefore, we compare the design process with and without using simulation. As a case study, we use a microfluidic network which is capable of trapping and merging droplets with different content on demand. The case study demonstrates how simulation can help to validate the derived design by considering all local and global hydrodynamic resistance changes. Moreover, the simulations even allow further exploration of different designs which have not been considered before due to the high costs.


Introduction
There have been numerous studies reported on droplet micro-uidics over the past decade as summarized by a number of review articles. [1][2][3][4][5] Fundamental studies mainly focus on investigating and elucidating two phase ow and transport phenomena as well as exploring functionalities for droplet manipulation, [6][7][8][9][10][11][12] while application driven studies aim to exploit droplet microuidics to address the challenges associated with the current best practices. [13][14][15][16][17][18] Both active and passive methods have been developed for manipulating droplets such as droplet generation, merging, splitting and trapping. Passive methods rely on the variation of applied pressures, geometries and uid properties to manipulate droplets and thus do not need external components to be integrated with microuidic devices. In general, passive methods are more favorable for large channel networks (i.e. large microuidic networks) where integrating multiple active components becomes challenging.
However, the functional performance of passive methods is highly dependent on the channel dimensions and their fabrication precision as well as the applied pressures because the entire microuidic network is coupled together meaning that any event (i.e. droplet generation, merging, splitting, or exiting of the channel network) occurring in the microuidic network would change the local and global hydrodynamic resistance and thus the pressure drop over different channel sections. Furthermore, design and operational uncertainties such as uctuations in the applied pressures, fabrication defects and dusts or particles entering the microuidic network are inevitable in microuidic studies.
Therefore, the design of the microuidic network needs to be smart enough to be insensitive to these uncertainties. In other words, the designer has to derive the dimensions of the channel network, the applied pressures, etc. so that the resulting speci-cation is as robust as possible. However, in this task, designers rely on their expert knowledge and derive the specication based on simplications as well as assumptions. For example, designers oen simplify or ignore time-dependent effects on the hydrodynamic resistance of microchannels caused by droplets because it is simply impossible to consider all droplet states and positions by hand.
In order to address these challenges, researchers such as Glawdel et al. 19 presented a series of global network design criteria which are supposed to aid designers in determining the correspondingly needed specication. But, while certainly helpful, these criteria only narrow down the microchannel dimensions to a certain range. Using them, still multiple designs are required by varying the channel dimensions i.e. width, length and height within the range set by the network design criteria to gure out the optimal design. Furthermore, recently also an automatic method was presented in ref. 20, which aids the designer in the dimensioning task.
In order to validate the respective designs (which can reach up to dozens depending on the complexity of the microuidic chip), usually physical experiments are conducted. In case that the functionality is not implemented as desired, the designer has to go back, revise the specication, and repeat the prototyping before he/she can conclude whether the revisions eventually lead to the desired result.
This design process apparently comes with major drawbacks: expensive physical experiments are used in order to validate designs. This oen results in multiple iterations for prototyping where each iteration requires the generation of a physical design and the execution of the experimentsresulting in high costs (nancially as well as in terms of time). Moreover, these circumstances frequently prevent the exploration of better and more advanced realizations (i.e. only several values can be picked in order to reduce the total number of prototypes).
In contrast, simulations can signicantly help validate the design before the rst prototype is fabricated. In fact, several simulation approaches have been introduced for simulating microuidic designs in the past. These simulation approaches can be classied in two abstraction levels: Simulations using Computational Fluid Dynamics (CFD): in CFD-simulations, the uid ow is described in the most detailed and accurate way which requires a complex simulation setup (e.g. the generation of a mesh based on the physical design). The high level of physical details causes high computational costs, which makes this kind of simulation most useful for (and also limits it to) small designs or parts of a larger design (e.g. single components). Reviews of corresponding methods are provided in ref. 21-23. Simulations using the one-dimensional (1D) circuit analysis model: this model reduces the microuidic network which is inherently a 3D object to an 1D hydraulic circuit, which has an analogy to an electric circuit. This model can be applied when the ow is laminar, viscous, and incompressible. 24 Furthermore, this model also considers the inuence of droplets on the hydrodynamic resistance of microchannels. 25 But it does not allow designers to predict effects like droplet mixing, formation or splitting. 25,26 The precision of simulations on this level has been systematically evaluated by using physical experiments in e.g. ref. 25-31. Corresponding simulators are especially suited to simulate large microuidic networks, even before starting with physical design.
For the purpose considered here, simulations on the 1D circuit analysis model are most suitable because they can be applied for large microuidic networks and are perfectly applicable for design exploration. This makes them especially suited for addressing the previously discussed challenges: simulations allow users to validate the derived specication by considering all before applied simplications and assumptions. Further, these simulations can be conducted at the beginning of the design process even before a physical design is drawn. If a simulation predicts that the chosen specication is not robust or even invalid, the probability is high that a prototype based on this specication will also show an incorrect functionalitywhich eventually reduces the number of tested prototypes. Furthermore, simulations allow the exploration of different designs.
But the potential of simulations on the 1D circuit analysis model has hardly been exploited by designers yet. In this work, we show the potential simulators can provide during the design of microuidic networks. To this end, we consider the design of a microuidic networkboth, in the "traditional" fashion (i.e. manually, with many prototyping iterations) and in a fashion where a simulation is additionally used. As a case study, we use the microuidic network presented in ref. 32, which consists of pairs of trapping wells. The working principle of the trapping well design is based on the hydrodynamic bypass concept, which has been reported by Takeuchi's group 33 and Vanapalli's group 34 for trapping single droplets. In ref. 32, those concepts have been improved for trapping pairs of droplets. The design does not require any electrodes, magnets, or any other moving parts and is entirely passively controlled by pressure and ow changes. In this study we show that, compared to the original design process which required six fabricated prototypes, one person month of an experienced designer, and nancial costs of USD 1200, using simulation can signicantly help when deriving the specication. More precisely, this study demonstrates that when using simulation, the designer is guided towards the design, which has been nally used in ref. 32. Moreover, the simulations even allow further exploration of new designs and e.g. prediction of their throughputs.
In the following, the applied 1D circuit analysis model is reviewed as well as the details of the considered microuidic network and its implementing application are provided.

Background
In this section, we rst review physical basics of microuidic networks. Aerwards, we describe the considered case study and its application and working principle of the design.

Applied model
The design of microuidic networks is usually conducted on the one-dimensional (1D) circuit analysis model, whose duality to electric circuit is reviewed in ref. 24. In this model, droplets can be considered as described in ref. 25. In the following, we briey review this model.
For microuidic systems at low Reynolds numbers, the relation between two points in a channel can be described by the Hagen-Poiseuille's law 35 with where Q is the volumetric ow rate Q (in [ml min À1 ]), DP is the pressure gradient (in [mbar]), and R is the uidic resistance (in [mbar min min À1 ]). This resistance of channels is constant for low Reynolds numbers and depends on the channel's geometry and the uid viscosity m c (in [mPa s] when given as dynamic viscosity) of the continuous phase. For example, the resistance of a rectangular channel (with length L, width w and height h), where the section ratio h/w is less than 1, is given by where a denotes a dimensionless parameter dened as However, the presence of droplets in channels change the pressures and ow rates in the microuidic network (i.e. the ow state) as they cause additional resistances. When the distance between two adjacent droplets is at least a few channel sections/diameters, their ow perturbations do not interact, 25 which allows the modeling of each droplet by an additional resistance. The overall ow resistance of a microchannel can be calculated by where R is the resistance of the channel, n is the number of droplets inside the microchannel, and R d is the single droplet resistance. The droplet resistance R d is dependent on the capillary number Ca, the viscosity ratio of the dispersed phase to the continuous phase m i m c , the interfacial tension g and also channel dimensions, which is described by 36 where a is a constant and L d is the droplet length.
which has been experimentally veried. Furthermore, this estimation of the droplet resistance is only valid for a conned ow (i.e. where droplet lengths are greater than the channel width). For an unconned ow (i.e. where the diameter of spherical droplets is smaller than the channel width), this equation needs to be changed. In order to determine the ow state in all channels, the mass conservation and the relation described by the Hagen-Poiseuille equation can be employed. 25 In detail, equations similar to the Kirchhoff's laws can be obtained by the following rules: The sum of ow rates into a node is equal to the sum of ow rates out of that node. A node is a point in the microuidic network where the ow splits or merges.
The directed sum of pressure gradients around any closed cycle is zero. The sign of the pressure gradients is dened by the specied direction of the ow rates.
By solving the resulting equation system, the ow state (i.e. DP and Q) in every channel for the current droplet positions can be determined. Note that this ow state becomes invalid as soon as a new droplet is injected, any droplet exits the network, or any droplet enters another channel (as this modies the resistances in the equation system).
The obtained ow rates in the channels allow the determination of the droplet speed by where a s is the slip factor. Under the conditions where the droplet length is between 1.5w and 7.2w, the viscosity ratios is 0.03 or 0.88, and the capillary number between 0.001 and 0.01 without surfactant, Vanapalli et al. 37 found the slip factor to be constant and equal to a s ¼ 1.28. By using the obtained ow state, also pressure gradients can be checked whether or not they exceed the Laplace pressure, which e.g. would cause a droplet to be squeezed through a gap. The Laplace pressure is dened by where g is the interfacial tension (in [N m À1 ]), w gap is the width of the gap, and r d is the droplet radius.

Considered case study
In this case study, we consider a droplet microuidic system, which can be used to screen drug compounds that inhibit the tau-peptide aggregation. 32 This phenomenon is related to neurodegenerative disorders such as Alzheimer's disease 38 and here protein misfolding and aggregation are considered to play a signicant role. Therefore, the screen process is to gure out the compounds that can inhibit protein aggregation. The use of droplets for this application allows the signicant reduction of the sample consumption volume by a factor of 100 as well as of the reaction time from 2 h to several minutes compared to the traditional 96-well plates.
Due to the droplet's large surface to volume ratio and fast mixing properties, the protein aggregation process reaches plateau within only 30 s. Hence, the reagents such as protein, uorescent dye, and inhibitor, cannot be premixed before droplets get trapped for real-time monitoring, i.e. the reaction needs to be triggered on demand. This real-time monitoring functionality and the therefore required operations are implemented as a passive solution (i.e. no active components are required) in ref. 32. For producing droplets of the two reagents, the design needs two independent droplet generators (i.e. two T-junctions). The key elements of the design are trapping wells (also shortly called traps in the following) as shown in Fig. 1. Each set of trapping wells allows two droplets be trapped, merged, and mixed and, hence, allows the precise control of the reaction time.
More precisely, when a droplet reaches the entrance of a trap and this trap does not already contain a droplet, the droplet has to enter the trap where the droplet should then stay (i.e. should not be pushed into the other trapping well or through any gap downstream). As soon as a droplet of the second reagent is trapped in the adjacent trap, both droplets merge and the reagents mix. On the other hand, when the trap already contains a droplet, following droplets should enter the bypass channel, which is again connected to further trapping wells. In addition to Fig. 1, the working principle of the trapping wells is also nicely illustrated by means of videos available at https:// doi.org/10.1039/c7ra02336g.
In order to implement this desired behavior, two objectives have to be fullled: (1) A droplet has to enter an empty trap: as a droplet always ows along the branch with the highest volumetric ow rate, the ow into the trap has to be larger than the ow into the bypass channel when the trap does not already contain a droplet, i.e. Q trap > Q bypass . On the other hand, when the trap already contains a droplet, succeeding droplets in the stream have to enter the bypass channel. This is guaranteed by drastically decreasing the ow into the trap when the trapped droplet clogs the two narrow gaps, i.e. by ensuring Q trap < Q bypass .
(2) A trapped droplet has to stay in the trap and must not be squeezed through any gap: rst, to prevent the trapped droplet from being squeezed through the gap between the traps, the pressure drop between two connected traps must be less than the Laplace pressure across their intersection. Fig. 1 shows the intersection between both traps and its width w i . Furthermore, when the trap contains a droplet, the droplet's radius is r d # r. This allows the denition of the objective as 32 Second, to prevent a trapped droplet from being pushed out of the trap, a similar objective has to be fullled, 32 i.e.
In this work, we aim to discuss the design of this non-trivial microuidic network and how this can be advanced using simulation. Therefore, we next review the current process of designing microuidic networks.

Current design process
In this section, we review and discuss the steps and respective challenges in the current process of designing a microuidic network. Therefore, we use the case study as introduced in the previous section.

Deriving the specication
Taking the desired behavior and the basic structure of the design (cf. Fig. 1), the dimensions of all channels, the applied pressures/ow rates, and the used phases have to be specied next, i.e. the designer has to derive a specication.
Here, some dimensions are given by the application or by fabrication limitations, e.g.: The droplets' size and spacing determine the trapping well radius r and the length L 1 , i.e. r is chosen so that a droplet lls most of the trap and L 1 is chosen so that the following droplets do not contact with the trapped droplet.
The fabrication limits the gap width w gap and gap length L 2 (i.e. the pillars require a minimum size, otherwise they could be peeled off from the silicon wafer).
The properties of the used phases are chosen so that they are suited for forming droplets and for the application. When a droplet reaches the entrance of a trapping well (i.e. point J 1 or J 2 ) and the respective trapping well does not yet contain a droplet, the droplet should flow into the trapping well. Therefore, the flow rate into the trapping well has to be larger than the flow rate into the bypass channel, i.e. Q trap > Q bypass . Furthermore, a trapped droplet should stay in the trapping well and, hence, must not be squeezed into the other trapping well (i.e. P trap 1 À P trap 2 must not exceed the Laplace pressure) and must not be squeezed through the gaps downstream (i.e. P trap 1 À P down must not exceed the Laplace pressure). On the other hand, when the trap already contains a droplet, following droplets should enter the bypass channel. Therefore, a trapped droplet has to decrease the flow rate into the trapping well so that the flow rate into the bypass channel gets larger, i.e. Q trap > Q bypass .
However, besides that, all remaining issues such as the length of the bypass channel or the applied pressures have to be explicitly specied by the designer. Considering that even small droplet-based microuidic networks are composed of several channels which, together with the droplets themselves, yield to plenty of interdependencies, thus this leaves a tedious task. In fact, the authors of ref. 32 spent a signicant portion of their work discussing the respective issues in detail, which is why we refer to this work for a more detailed treatment.
Overall, deriving the specication requires the denition and consideration of plenty of variables which all affect the intended behavior (i.e. whether the objectives are fullled). In this task of deriving the specication, the designer relies on his/ her expert knowledge and oen applies simplications and assumptions. Especially, the time-dependent resistances caused by droplets are oen simplied or completely ignored when deriving the specication as it is impossible to consider all droplet states and positions by hand. For example, in the specication of ref. 32, only for the bypass channel a x number of contained droplets is assumed, which causes an additional resistance.

Prototyping and testing
The task reviewed in the previous section eventually yields a specication which is supposed to realize the intended behavior. However, due to applied simplications and assumptions, the designer cannot be sure whether an implementation based on the specication indeed realizes the desired behavior and whether all objectives are fullled under all settings.
Therefore, as a next step, it is tested whether the specication realizes the desired functionality by using (physical) experiments. To this end, the designer has to fabricate the design based on the derived specication. This rst requires a physical design (also called layout or mask) to be drawn from the specication (e.g. as a vector graphic), which can be used as input for the production process. Aerwards, this is used to fabricate the device using a so-lithography technique, 3Dprinting, or milling. The resulted device then is validated, i.e. experiments are conducted to check whether the device indeed shows the desired behavior. This is the rst time in the design process in which the designer can observe the effects of his/her choices and decisions as well as simplications and assumptions during the derivation of the specication (which provided the basis for the physical realization) on the actual behavior.
In the case that those validations show that the behavior has not been implemented as desired, the specication has to be rened, i.e. the dimensions of some channels, the applied pressure/ow rates, or some of the used phases have to be adjusted. In particular at the beginning of the design process, this of course is likely needed as assumptions might be inappropriate and simplications may have led to an imprecise specication. But then, the entire process of creating a physical design and fabricating another prototypical device has to be repeated in order to test again whether the (now rened) design is correct. This iteration of rening and testing constitutes one of the major drawbacks of today's design process as it requires a signicant amount of time and resources.
In fact, in the case study conducted in ref. 32, a total of six different specications had to be derived, fabricated, and tested until the desired behavior was accomplished. More precisely, prototypes with bypass channel lengths and trapping well gap widths as shown in Table 1 have been fabricated.
Each of these prototypes has experimentally been tested for the trapping efficiency, and the one with ID 2 (i.e. a bypass channel length of 4000 mm and trapping well gap width of 15 mm) eventually showed the desired performance with respect to trapping robustness.
Overall, this resulted in a working time of one month for an experienced designer and nancial costs (including silicon wafer, SU-8 photo resistor, fee charged for clean room, Polydimethylsiloxane (PDMS), and silicone oil, etc.) of approximately USD 200 for a single prototype (i.e. a total of USD 1200 until the desired design eventually worked).

Further missed potential
The complexity of the design process as reviewed above does not only pose a challenge to get a design realizing the desired behavior. Moreover, it also prevents further improvements which, in principle, could be conducted but are too expensive in most cases. Hence, as soon as a prototype shows a correct behavior, the respective design and its specication is usually xed and no more different designs are explored. However, probably different designs would be even more suited, e.g. would be more robust, would be smaller, would increase the throughput or would have a positive effect on the application.
For the case study, no design has been explored to answer the following questions: (1) Question 1: What is the minimal bypass channel length so that droplets still get trapped?
(2) Question 2: What is the maximum pressure over ve sets of trapping wells so that no objective is violated?
(3) Question 3: How many trapping wells can be cascaded and loaded by droplets in a given time, i.e. what is the maximal throughput?
Why these questions are important and how to address them will be discussed later in Section 4.2.2.
Overall, it can clearly be seen that the design process as conducted thus far is certainly not ideal. Deriving a working specication requires several iterations of physical design and Paper prototyping which becomes a time-demanding and costly process. Moreover, as a consequence, oen the rst working design is eventually used even if better and more advanced solutions would, in principle, be possible.

Improved design process using simulation
In this work, we aim for improving the design process reviewed above by using numerical simulations based on the 1D circuit analysis model. These simulations allow the iterations without the need for a physical design nor the testing on an actually fabricated device. How this is accomplished is covered in this section. To this end, we rst review the used simulator as well as its features and application. Aerwards, we revisit the issues discussed in the previous section and describe how simulation can help to addressed them in an efficient way.

Used simulator
In this section, we describe the used simulator, which is applicable for simulating droplet microuidics operated in a laminar, viscous, and incompressible ow regime. The descriptions are provided in two parts: rst, we review its general working principle, which is supposed to provide a background on the internals of the simulator but can easily be omitted by readers who are only interested in a user-perspective of the simulator. Aerwards, we summarize the interface and core features which result from this working principle and, hence, are available to designers. 4.1.1 Internal working principle. In this work, we employ the simulator presented in ref. 39, which has been implemented in Java (which makes the tool platform-independent) and can be downloaded at http://iic.jku.at/eda/research/ microuidics_simulation.
The simulator automatically derives, applies, and solves microuidic equations based on the model described in Section 2.1. Therefore, the microuidic network is represented as a directed graph. The edges in the graph represent channels and pumps and their directions represent the counting direction of the ow rates. The nodes represent points in the microuidic network where edges are connected. This directed graph is used to automatically derive an equation system by applying the Kirchhoff's laws as described in Section 2.1. More precisely, the algorithm generates an equation for each node (i.e. in each node the sum of all incoming ow rates is equal to the sum of all outgoing ow rates) as well as for each cycle (i.e. the sum of all pressure gradients around a closed cycle is equal zero). For determining the cycles in the graph, a variant of Paton's algorithm 40 is used. In this equation system, also the effects of droplets are considered. However, the ow state becomes invalid as soon as a new droplet enters a channel/exits the network. Consequently, the entire equation system has to be frequently re-evaluatedwhich would be impossible by hand. The simulator conducts this automatically and solves the equation system by applying a lower-upper (LU) decomposition. 41 Overall, the simulation algorithm performs the following steps in a loop: (1) Compute ow state: the simulator computes the ow state (i.e. the pressures and ow rates in all channels) by considering all droplets and their positions. Therefore, the simulator automatically derives an equation system which considers the resistances caused by all droplets and solves it. The resulted ow rates allow the determination of the droplet speeds, i.e. by eqn (7).
(2) Check objectives: based on the obtained ow state, the objectives (e.g. a pressure must not exceed the Laplace pressure) are checked and in case of any violation, the simulator informs the designer.
(3) Compute next event-time: the next event time is computed by the minimum time until a new droplet is injected or any droplet enters the next channel/exits the network (here the droplets' speeds are used).
(4) Update system state: nally, the simulator updates the system state to this event time (i.e. the droplet positions and their resistances). For the droplet causing the event, the ow state decides which channel this droplet enters next (i.e. the channel with the highest ow rate).
Summarizing, the simulator re-calculates the ow state when the old one becomes invalid. This allows the simulator to consider all interdependencies caused by droplets. Furthermore, these event-based calculations make the algorithm efficient. More precisely, computing the ow state can be done in polynomial time with respect to the number of variables in the equation system and checking the objectives, computing the next event-time, as well as updating the system states is computationally inexpensive. This algorithmic efficiency allows to simulate large microuidic networks in negligible computation times.
4.1.2 User interface and features. The designer can use the simulator as a black box and, as already stated, he/she does not need to understand the algorithm in detail. Further, he/she does not need to provide a physical design yet as all simulations are conducted based on the information which is available to the designer anyway when he/she has completed the speci-cation. More precisely, for setting up a simulation, the designer only has to provide the specication which describes the channel dimensions, the structure of the network (i.e. how are the channels connected), the applied pressures or volumetric ow rates from the inlets, the properties of uids (i.e. the viscosities, densities, interfacial-tension), and the objectives which have to be fullled. Aer providing the specication, the simulator allows the analysis of the instantaneous ow state of all channels in the micro-uidic network, the trajectories of droplets (i.e. the droplet switching at bifurcations or whether they enter a trapping well) and their spacing and patterns, the droplet speeds and the residence time (i.e. the time for passing a channel), as well as the objectives whether they are fullled. The simulator's features can be observed in the video at http:// iic.jku.at/eda/research/microuidics_simulation, which shows the graphical output when simulating the chip of ref. 32 (cf. also considered in our case study). Additionally, Fig. 2 provides a comparison of the simulator and the physical experiment. On the le side of this video and gure, the graphical output of the simulation is shown which shows the droplet positions. On the right side of this video and gure, the simulation is compared to a physical experiment. Therefore, we used images of a physical experiment captured with a frequency of 50 fps. When comparing the simulation with the experiment, we can observe equal droplet trajectories as well as similar droplet speeds and residence times. Furthermore, the simulator constantly checks the objectives introduced in Section 2.2. The obtained results prove the precision with respect to a physical implementation.

Application of the simulator
The simulator and its features as reviewed above can now be utilized in order to address major shortcomings of the current design process for microuidic networks as reviewed in Section 3. Recall that the designer derives the specication based on manual calculations, simplications, as well as assumptions and can, thus far, only validate the specication by fabricating and testing the resulted prototype. Utilizing the simulator, many of these tests can now be conducted much earlier in the design process and without the need for either a physical design or a fabricated chip. This additionally allows for a much more elaborated design exploration as variances of the specication can be validated much easier, with signicantly less manpower, and basically no cost. In the following, this is demonstrated by (1) revisiting the design process reviewed in Section 3.1 as well as Section 3.2 with respect to the features provided by the simulator and (2) illustrating the further potential which can be gained by the improved possibilities with respect to design exploration (and, thus far, has been missed as discussed in Section 3.3).
4.2.1 Utilizing the simulator in the design process. In ref. 32, the designer came up with six different specications to be tested using physical experiments. Instead of drawing six physical designs, fabricating respective prototypes, and conducting the experiments, in the following we evaluate these six specications by utilizing the simulator.
The simulation of the specication with bypass length L bypass ¼ 3000 mm and gap width w gap ¼ 15 mm predicts that the ow into an empty trap is equal to Q trap ¼ 2.054 ml min À1 and the ow into the bypass channel is equal to Q bypass ¼ 1.871 ml min À1 . Due to the fact that a droplet always ows along the path with the highest ow rate, a droplet enters the empty trap under perfect conditions. However, the ow rate ratio does not allow a robust decision of the droplet path. For example, even a small particle blocking the ow into the trapping well would reverse the ratio and, therefore, would violate the Objective 1 (Q trap > Q bypass , cf. Section 2.2).
Using simulations, we found that this ow rate ratio between the trapping well and the bypass channel increases when the bypass channel length is increased and/or the gap width is increased (which is the case for all other specications). However, an increase of the bypass channel length also causes an increase of the time that a droplet requires to pass this bypass. For example, the simulation showed that, when increasing the bypass channel length from L bypass ¼ 4000 mm to L bypass ¼ 5000 mm (using a gap width of w gap ¼ 15 mm), the time that a droplet requires to pass a single bypass increases from 0.25 s to 0.32 s (which is an increase of 28%). This effect adds up when multiple trapping wells are cascaded. When additionally a certain upper limit for the loading-time of droplets (i.e. especially relevant for cells) has to be fullled, less trapping wells can be cascaded using a longer bypass channel, which decreases the throughput.
On the other hand, the time for droplets to be trapped can be reduced by increasing their speed, which can be achieved by applying higher pressures at the inlets. Using the simulator, we tested different pressures for all six specications and the simulator predicts that the pressure is limited by the Objective 2 (cf. Section 2.2). More precisely, too high pressures at the inlets cause too high pressure gradients between the trapping wells and the narrow gaps so that the condition described in eqn (10) is violated. This causes the trapped droplet to be pushed through the two narrow gaps. The simulation shows that droplet speeds can be higher for the prototypes which have a gap width of only 15 mm. This is because a smaller gap width increases the Laplace pressure. † These simulation results now allow the designer to evaluate the robustness and performance of the different specications. Table 2 summarizes the obtained insights. Based on these results, the designer can evaluate the specications with respect to their robustness. The results show a clear preference for the specication with ID 2. Exactly this specication is the one which was eventually realized in ref. 32 (cf. Section 3.2). The appealing features of this simulator are clear: instead of six fabricated prototypes, 1 person month of manual labor, and a total of USD 1200 of further costs, utilizing the simulator would allow for obtaining the same result by only fabricating one single design, spending only 1/6 of the previously spent time of manual labor plus few hours for simulating, and a total of USD 200. Apparently, this is a signicant improvement compared to the design process applied in ref. 32. 4.2.2 Utilizing the simulator for design exploration. As discussed in Section 3.3, the complexity of the design process is not only a burden to get a design realizing the desired behavior, but also prevents the designer from exploring even better solutions. Utilizing the simulator, this burden is signicantly reduced. In fact, as described next, all questions raised in Section 3.3 can now efficiently be addressed.
Question 1: What is the minimal bypass channel length?
Here, we explore the design with respect to the limits of the bypass channel length. Recall, a short bypass channel length decreases the time for a droplet to be trapped and hence affects the throughput of the design. Furthermore, a short bypass channel length also minimizes the area of the physical design, which is an important criterion of the design due to the limited space of a typical microuidic chip.
The bypass channel lengths determine the ow rate into the bypass Q bypass and also into the trapping well Q trap . Therefore, the bypass channel length is limited by the Objective 1 (cf. Section 2.2), i.e. the ow rate into an empty trapping well has to be larger than the ow rate into the bypass channel.
In order to explore the limits, we conducted simulations where we stepwise reduce the bypass channel length in the specications and stop as soon as the rst droplet does not enter an empty trapping well (which violates the objective Q trap > Q bypass ).
For a gap width of w gap ¼ 15 mm and for a gap width of w gap ¼ 25 mm, the simulation predicts a minimal bypass channel length of L bypass ¼ 2800 mm and L bypass ¼ 900 mm, respectively. Both values are only theoretical limits and respective designs would be sensitive to any imperfections as e.g. dust or imprecisions caused by the fabrication process. Therefore, a prototype will never be pushed to these limits.
Instead these limits allow the designer to estimate the robustness of the design. For example, this design exploration would have prevented a prototype with L bypass ¼ 3000 mm and w gap ¼ 15 mm (i.e. the specication with ID 1 from Table 2) to be considered in the rst place.
Question 2: What is the maximum pressure over ve sets of trapping wells so that no objective is violated?
By increasing the pressure at the inlets, also the droplet speed is increased and, accordingly, a higher droplet speed decreases the time required for a droplet to travel from its injection until it gets trapped. This might be crucial since, the droplet loading time is limited for some biological experiments and, therefore, the droplet speed is an important factor for the throughput of a design.
In this design exploration, we explore for the considered design with ve trapping wells (as proposed in ref. 32) the maximal possible pressures across these trapping wells. Again, we increment the pressures until any objective is violated. The simulation results reveal that the Objective 2 limits the maximal pressures. More precisely, for too high pressures, the simulation predicts that the pressure across the trapping well and the narrow gaps exceeds the Laplace pressure, which would cause a droplet to be squeezed out of the trap. Table 3 shows the obtained results for the six specications proposed in Section 3.2. These results show that a smaller gap width allows for a higher pressure drop. This can be explained because the Laplace pressures is higher for smaller gaps. Interestingly, also the shorter the bypass channel, the higher the possible pressure. These results again conrm that a gap width of only 15 mm is more robust for higher pressures.
Question 3: How much trapping wells can be cascaded and loaded by droplets in a given time?
The design proposed in ref. 32 cascades ve trapping wells. However, it would be possible to cascade more trapping wells in order to increase the throughput. This number of cascaded trapping wells determines the time until all droplets are trapped, i.e. the loading time. The maximal allowed loading time  depends on the bio-assays and can particularly be relevant to cells. Furthermore, the pressure drop over the trapping wells must not exceed operating settings, i.e. exceeding pressures will bow the PDMS channels. Typically, the pressure applied to PDMS microuidic devices is limited to 5 bar.
In the current design process, it would be costly to explore designs with different numbers of trapping wells and measure the required loading time. Therefore, it is currently unexplored how many trapping wells can be cascaded so that all droplets can be trapped within a certain maximal loading time. Utilizing the simulator, now also this question can be easily explored.
Therefore, we create specications containing between 15 and 45 pairs of cascaded trapping wells. For each of these designs, three different pressures (100 mbar, 200 mbar, and 300 mbar) over the trapping wells are applied. These pressures and especially the maximal pressure have been selected so that no objective is violated. Then, the simulator is used to measure the time which is required for all droplets to be trapped, i.e. the overall loading time. Fig. 3 summarizes the obtained results for the three different pressures. We can see that the loading time increases with the number of cascaded trapping wells and is generally lower when higher pressures are applied. This increase can be explained because (1) the distance to the last set of trapping well increases and (2) the droplet speed decreases as the overall resistance caused by the trapping well increases. Overall, also here the simulator allows for obtaining results that were not available beforehere, to implement a bio-assay with a maximal loading time and a maximal throughput.

Conclusion
In this study, we demonstrated how simulation can help in the design process of droplet microuidics. Therefore, we conducted a case study comparing the design processes of the microuidic design proposed in ref. 32 with and without simulation. When no simulation was used, the designer validated the specication by using physical experiments, which resulted in six prototypes, took one person month and produced nancial costs of USD 1200. If the designer is not experienced, the costs could even be much higher. Instead, when a simulation was used (the used simulator is available at http://iic.jku.at/eda/research/microuidics_simulation), the designer was able to validate the specication before any prototype and even before any physical design was made. These simulations allow the prediction of the robustness of the respective specications and, in fact, show a clear preference for the specication which was eventually realized in ref. 32. Hence, the simulations allow for selecting the most robust design without the need to explicitly fabricate and test them. Furthermore, the simulations even allowed to explore further designs, which were too costly to consider in the current design process without simulation.
Overall, the use of simulation on the 1D circuit analysis is perfectly suited to address the wide variety of microuidic networks and their respective design challenges. Simulations allow for a quick validation and exploration of arbitrary microuidic designs (the setup of the simulation is hardly any work compared to physical experiments as well as the computation time sums up to at most a few seconds). Furthermore, they signicantly help to increase the robustness of the design as well as to accelerate the design process and, by this, reduce the overall costs.

Conflicts of interest
There are no conicts to declare.