Design of multi-functional microfluidic ladder networks to passively control droplet spacing using genetic algorithms

Computers and Chemical Engineering 60 (2014) 413–425

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Design of multi-functional microfluidic ladder networks to passively control droplet spacing using genetic algorithms Jeevan Maddala, Raghunathan Rengaswamy ∗ Department of Chemical Engineering, Texas Tech University, Lubbock, TX, United States

a r t i c l e

i n f o

Article history: Received 16 January 2013 Received in revised form 8 July 2013 Accepted 15 September 2013 Available online 2 October 2013 Keywords: Droplet-based microfluidics Genetic algorithms Ladder networks Synchronization

a b s t r a c t Accelerated progress in the use of droplet-based microfluidics for high throughput screening and biochemical analysis will require development of devices that are robust to experimental uncertainties and which offer multiple functionalities. Achieving precise functionalities in microfluidic devices is challenging because droplets exhibit complex dynamic behavior in these devices due to hydrodynamic interactions and discontinuities that are a result of discrete decision-making at junctions. For example, even a simple loop device can show transitions from periodic to aperiodic/chaotic behavior based on input conditions. Hence, rational design frameworks that handle this complexity are required to move this field from labs to industrial practice. Two main challenges that need to be confronted in the realization of such a rational design framework are: (i) computational science related to rapid simulation of very large networks; development of predictive models that will form the basis for characterizing droplet motion through interconnected and intricate large-scale networks, and (ii) conceptualization of a design approach that is generic in nature and not very narrowly defined limiting its application potential. In this paper, we develop a GA approach for the design of ladder networks that are used to control the relative droplet distance at the exit. Through several case studies, the potential of the proposed GA approach in designing exquisite ladder structures for multiple functions is demonstrated. A recently proposed network model is used as the basis for all the computational studies reported in this paper. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction The ability to control small volumes of fluids through dropletbased microfluidics has revolutionized conventional biochemical analysis (Ho, Chakrabarty, & Pop, 2011; Su, Chakrabarty, & Fair, 2006). Recent advances in this area have led to the development of devices with basic functionalities such as sorting (Christopher et al., 2009; Cristobal, Benoit, Jaonicot, & Ajdari, 2006; Hatakeyama, Chen, & Ismagilov, 2006), merging (Jin, Kim, Lee, & Yoo, 2010), synchronization (Prakash & Gershenfeld, 2007) and storing (Boukellal, Selimovic, Jia, Cristobal, & Fraden, 2008; Laval, Crombez, & Salmon, 2009; Shim et al., 2007). Using these basic devices it is now possible to design massively parallelized large-scale architectures for lab-on-chip applications. The current state-of-the-art is that these devices are developed using a bottom up approach i.e. by rigorous experimentation to find optimal designs for a given objective. Advanced computational techniques have the potential to accelerate this design discovery process through a top down approach i.e., rigorous simulations to discover optimal designs for the given task, which could then lead to experimental validation. Reliable

∗ Corresponding author. Tel.: +1 806 742 1765. E-mail address: [email protected] (R. Rengaswamy). 0098-1354/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2013.09.009

computational approaches can have a tremendous impact in this area similar to the impact of simulations in the electronics industry. Two main challenges that need to be confronted in the realization of such a computational framework are: (i) computational science related to rapid simulation of very large networks; development of predictive models that will form the basis for characterizing droplet motion through interconnected and intricate large-scale networks, and (ii) conceptualization of a design approach that is generic in nature and not very narrowly defined limiting its application potential. Development of predictive models will have to address the inevitable complex nonlinear dynamics of dropletbased microfluidic systems (Cybulski & Garstecki, 2009; Epstein, 2007; Sessoms, Amon, Courbin, & Panizza, 2010). For example, it has been shown that binary decisions that droplets make at bifurcations can lead to fascinating nonlinear dynamics, including chaos and fractal behavior. Nonetheless, rapid progress is being made in characterizing the dynamic behavior related to the various phenomena that occur in droplet-based microfluidic devices. On the other hand, there is very limited work that can be found on rational design approaches for discovering large-scale microfluidic networks. To address this gap, we propose a genetic algorithm based design framework in this paper. While the ultimate vision is the development of a design discovery process for very general and complex microfluidic networks

414

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Fig. 1. Ladder network with a combination of backward, vertical and forward bypasses.

that incorporate several design elements, as a first step towards this goal, we develop and demonstrate our approach for a novel structure, the ladder networks. We begin with a brief introduction to the ladder networks, the importance of this structure, and the design questions that need to be answered. We then discuss the proposed approach, which consists of two components: (i) a model for ladder networks and (ii) a design algorithm. The network model that is used to simulate the dynamics of the ladder networks and a customized genetic algorithm (GA) approach for design are discussed subsequently. The capability of the proposed GA formulation is demonstrated through three case studies: (i) synchronization of droplets for a particular input spacing (ii) synchronization of droplets for a range of input spacing, and (iii) realization of constant delay between droplets exiting the ladder network. 1.1. Microfluidic ladder networks Several papers have studied the ladder networks in detail. The ladder networks are important in their own right and have been used in microfluidic networks for synchronization. Synchronization is an operation that is used to maintain a zero or relatively small time difference between two events such as droplets exiting two arms of a ladder. That is, droplets entering the two arms of the ladder at different times but leave the ladder at the same time. A microfluidic ladder network (MLN) device contains a pair of parallel channels interconnected with bypasses. Droplets are transported through the parallel channels, whereas the bypass channels are inaccessible to the droplets. Fig. 1 represents an example ladder network, the first bypass slopes backward, and is defined as a backward slant; similarly, the middle and the last bypasses are called as vertical and forward slants, respectively. The bypasses of the MLN are defined based on their structural orientation—in reference to the direction of motion of the droplets. The spacing transformation achieved by a ladder network is defined as its functionality. Fig. 1 shows a pair of droplets entering the ladder network with an inlet spacing of Xin , this spacing is transformed to Xout as the pair of droplets exit the ladder network. The change in droplet spacing is due to the change in the relative velocity between the two droplets. The change in the relative velocity is a result of redistribution of flow across the bypasses. Symmetric ladder networks, i.e., networks with equally spaced vertical bypasses are shown to reduce the droplet spacing (Prakash & Gershenfeld, 2007) whereas asymmetrical ladder networks are shown to expand, contract, flip or synchronize the droplets (Maddala, Wang, Vanapalli, & Rengaswamy, 2012). In general, symmetric ladder networks offer less flexibility in terms of functionality compared to asymmetrical ladder networks. Additionally, asymmetric ladder networks have more parameters to optimize as compared to the vertical ladder networks. While the underlying physical principle of the ladder networks is deceptively simple, yet trying to design a robust ladder network with minimal device foot-print under unavoidable experimental constraints is not trivial. For example, in the pursuit of a ladder

network for robust synchronization, one might ask what is the optimal number of bypasses that is needed and how does this number depend on droplet hydrodynamic resistance and velocities? What degree of experimental uncertainties can be tolerated while still achieving acceptable synchronizer functionality? Are there alternative and more efficient structural configurations for the synchronizer than that implemented by Prakash and Gershenfeld (2007)? Going beyond the synchronizer, is it possible to design a ladder network for multiple functions such as input delay dependent behavior? Are there any advantages that can be accrued by slanting the arms of the ladder either in the forward or backward directions? To answer these questions, currently neither systematic efforts have been devoted nor rational approaches have been developed. Our proposed approach addresses such questions, specifically in the context of ladder networks. More generally, we view our approach as a first step towards the development of a rational tool for the design of complex microfluidic networks.

2. Proposed approach The ladder design problem comprises both structural and real variables. The structural elements of the problem are related to the topology such as the number of bypasses in the network. The continuous variables are the cross-sections of the bypasses and the positions of the bypasses on the main channels. This combination of structural and continuous variables will be a characteristic of almost all microfluidic design problems. The other important characteristic is the presence of constraints and bounds on the variables. These constraints could arise when minimum bypass separations are specified or if one forbids the bypasses to cross. The objective function has to reflect the goal of the design. The definition of the objective function for general microfluidic designs will be problem-specific. Optimization approaches to process design have been well studied in chemical engineering. Structural optimization problems are solved using combinatorial algorithms and continuous variable optimization approaches are largely gradient-based. Based on the problem characteristics described above, mathematical formulations for microfluidic design in general and the synchronizer problem in particular will be of the Mixed Integer NonLinear Programming (MINLP) variety. Branch-and-bound is an effective approach to solve large-scale MINLP problems (Parker & Rardin, 1988). One difficulty in using a standard MINLP formulation for microfluidic design would be the lack of closed form analytical expression for the objective function for typical design problems. The objective function has to be calculated through numerical simulations using a model such as the one discussed subsequently. This would make repeated gradient evaluations extremely numerically cumbersome. In the paper, we pursue the use of alternate non-gradient approaches for microfluidic design. In particular, we will explore the use of GAs. While other alternatives such as simplex search algorithms exist, GAs have been shown to solve combined structural and continuous variable problems efficiently, for example, in heat and mass exchange network design (Garrard & Fraga, 1998; Ravagnani, Silva, Arroyo, & Constantino, 2005). Another often cited advantage of GAs is that several acceptable designs are identified due to the principle of hyperplane sampling (Whitley, 1994). Further, in general it is reasonably straight-forward to setup a GA search (Goldberg, 1989; Holland, 1975). The general principles that underlie GAs are loosely based on ideas from evolution. In a genetic search, from an initial population of chromosomes, a set of best fit chromosomes (elitist selection) are chosen for reproduction. The parameters in an optimization problem are genes in a chromosome. These genes or parameters can be real numbers, integers or

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

415

even complex data structures. The chosen chromosomes produce newer chromosomes (children) through the use of genetic operators. The children are added to the initial pool. The algorithm is stopped using a termination criterion such as the maximum number of iterations or population stability. A well-defined and efficient genetic search is expected to lead to a final population that is much fitter than the initial random population. From a numerical optimization viewpoint, the fact that the solutions are sampled from different regions combined with the search of promising regions (elitist selection) leads to several acceptable solutions. This is one of the reasons why genetic search is a viable alternative for global optimization problems. A criticism of GAs is that theoretical convergence results are hard to establish. Further, no guarantees on the identification of globally optimal solutions can be made. We will now describe the two components of our design approach: modeling of ladder networks and the GA approach specialized to the ladder network design problem.

proposed GA is presented in Fig. 2. In the first step, MLN properties, such as dimensions, along with fluid properties are specified. The structural design parameters are coded into binary strings. GA requires an initial population to be specified; hence, we generated a random initial population of M networks. Later, as the algorithm progresses, these random structures are modified using genetic operators such as mutation and one point cross-over. The application of the genetic operators on the parent structures lead to an increase in the overall population. When the overall population doubles its initial size i.e., 2M, the selection operator, based on a fitness function, is used to reduce the overall population size back to M—this whole process is defined as one generation. The fitness function is calculated using the network model. Iterations are carried out until the total generations reach a maximum user defined value. The final population consists of several ladder designs that are expected to satisfy the desired design objective. The various components of the algorithm are discussed in detail in Fig. 2.

2.1. Network model

2.2.1. Coding ladder networks The design parameters of the MLN are coded into binary and real values. The coding has to result in a one-to-one mapping. Additionally, it should be possible to span all the possible designs of ladder networks by manipulating the real and binary values of the coding structure. The coding schema is described below:

In our work, the droplet dynamics in the ladder networks are characterized using a network model (Schindler & Ajdari, 2008). Prior studies have shown that the predictions using this model agree reasonably well with the experiments (Maddala et al., 2012; Song, Zhang, & Hu, 2012; Jeevan, Babji, Swastika, Vanapalli, & Rengaswamy, 2012). The basic modeling equations that are used in network approach are described next. Step 1: Resistance of the microfluidic channel is calculated based on its dimensions using the following equation (Bruus, 2008):



 1 192 h 12L R= 3 1− tanh r 5 5 w h w ∞

r

 rw  2h

−1

(1)

where h, w, and L are the height, width and length of the channel, respectively. The bulk phase viscosity is  and r is the series index which takes only odd values. The ladder network has a length of 10.1 mm and the transport channels are separated by 504 ␮m. The total resistance of the channel is calculated as Rt = R + Rd , based on the channel resistances (R) and droplet resistance Rd = 1.55 kg/mm4 S. A node is defined as a point where more than two microfluidic  channels intersect. From mass balance at each node, Qnode = Qi , where Qnode is the flow into or out of the node. Flow in a branch is Q = Rt /P. Step 2: Velocity of the droplet (Vd ) is based on the bulk flow rate in the channel; Vd = Qi —where the slip factor () is set as 1.4. Step 3: In the simulations, droplets move with the velocities calculated in step 2 until they encounter a node. Step 4: As the droplet enters a node, the resistances and velocities are recalculated using steps 1 and 2. The algorithm is repeated until both the droplets exit the system.

• Total bypasses (TB): This is a binary coding of the total number of bypasses that are present in the ladder network. If the length of the binary string is p, then the maximum number of bypasses (N) is equal to 2p − 1. As an example, 1011 would be decoded to 11 as the total number of bypasses. • Real value (V): Each bypass is associated with two real values (positions in the top and bottom main channels) resulting in a vector of 2N real values (V), where N is the maximum number of bypasses considered in the algorithm. Although the total number of bypasses in the ladder network can vary from 1 to N, the length of the array is kept constant at 2N for computational flexibility. The real values are between 0 and 1 and represent the positions normalized with respect to the total length. • Choice of top segment (SCt ): This is a binary string with as many entries as the maximum number of allowable bypasses. This group contains binary values which decides the presence or absence of a bypass at a specified location. For example {1, 1, 0, 1} here would imply that there are 4 maximum possible bypasses; however, since there are only three ones the total bypasses are three. Additionally this string also gives information that these bypasses are present on the first, second and fourth nodal positions in the top main channel. • Choice of bottom segment (SCb ): This is the same as the category above but specifies the positions in the bottom channel. The number of ones in the SCb string are exactly equal to the SCt string, this ensures that for every node on the top branch there is one node on the bottom branch.

2.2. Genetic algorithm We will now describe the GA approach for designing ladder networks for any given objective function. GA optimizes the structure of ladder networks by modifying their design parameters, which are the total number of bypasses and their positions defined by where they meet the top and bottom channels. Though the width of the channel and the bypasses can be included in the design parameters, these are fixed values for a particular design exercise. The constraints in the design of ladder networks include: (i) the bypasses of the ladder network are not allowed to intersect and (ii) the minimum spacing between two adjacent nodes is greater than a threshold, which is a fabrication constraint. Both binary and real value coding is utilized. The overall implementation of the

A remark is in order here. It seems like if the coded structure has just the TB and V attributes, the ladder architecture can be decoded using some simple logic such as using the first K (out of maximum N) bypasses based on the number K that the binary string TB represents. One of our initial versions of the GA did use this approach; however, it was very hard to generate diverse solutions using this approach as many of the structures had bypasses confined to the front of the ladder. Further, highly slanted bypass designs are difficult to generate, and even when generated through extreme mutation, there is an inherent limit on the slant due to the required threshold separation between the bypasses. The definition of SCt and SCb binary strings elegantly avoids this difficulty as we will see in the results section. Conversely, once SCt and SCb

416

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Fig. 2. Genetic algorithm implementation flowchart for designing microfluidic ladder networks.

are defined then it seems like TB is not necessary. However, direct mutation of even one digit in the TB string has the possibility of introducing a completely different structure in the population pool and this is useful in generating a diverse solution pool. Fig. 3(a) illustrates the coding scheme for a three bypass ladder network. The TB value of ‘0 1’ means that the ladder network has only one bypass. The ‘V’ group contains six real values (two times maximum three bypasses). The segment choice binary strings, SCt and SCb contain {1, 0, 0} and {0, 0, 1}, respectively. This implies that a single bypass is present between the first node of the top channel and the last node of the bottom channel.

2.2.2. Generating random population Several coded structures are generated randomly for the first generation of the GA. A total of 2N real random numbers between 0 and 1 are generated for the value group (V). The values are generated in such a way that the minimum distance between adjacent nodes is greater than a specified threshold. The maximum number of bypasses is limited to a pre-specified N. Random binary values of length P such that 2P − 1 = N with 1 in at least one position, so that a minimum of one bypass is present in every ladder network are generated for the TB binary string. Similarly, random binary numbers of length N are generated for SCt and SCb while ensuring that the

Fig. 3. Genetic operators implemented on MLNs: (a) coding scheme, (b) mutation and (c) one point crossover.

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

number of ones in the binary SCt string matches with the SCb binary string, and the total number of ones in SCt or SCb binary string is equal to the integer value of the decoded TB string. 2.2.3. Genetic operations Genetic operators are key components of any GA approach. The three types of genetic operations implemented in the proposed algorithm are: (i) mutation, (ii) one point crossover and (iii) selection. 2.2.3.1. Mutation. The mutation operator (changing a bit from 0 to 1 and vice versa) modifies the structure of the ladder network by changing the number and orientation of the bypasses thereby increasing the population diversity. Mutation either increases or decreases the number of bypasses. In the current algorithm, mutation is performed in two stages. In stage 1, mutation is implemented on ‘TB’ string and in stage 2, the mutation operator is implemented on TB string (with probability 0.5), SCt and SCb strings (with probability 0.5). Mutations in the ‘TB’ string are strong mutations that could result in large changes to the number of bypasses. However, stage 2 mutations in SCt or SCb will result in the number of bypasses changing by just one and hence, is a weaker mutation. Stage 1 mutations are implemented for the first 200 generations and after that, to achieve more nuanced variations in ladder structures, stage 2 mutations are implemented. The idea behind this is that at the initial stages of the search, one would like to identify several structures in as many local minima regions as possible and at the later stages, search more carefully within the identified local minima regions. Stage 1 mutations help in the former and stage 2 mutations help in the latter. Fig. 3 explains stage 1 and stage 2 mutations on the ladder networks. In case of stage 1 mutation, the binary string is mutated from ‘0 1’ to ‘1 1’, this results in adding two new bypasses to the MLN. The new bypasses are obtained by randomly selecting the top and bottom positions. This achieves the design shown in Fig. 3(b) (left). In the case of stage 2 mutation, the ladder designs are modified by changing just one bypass. This is achieved by selecting a random location of SCt string and its value is toggled from 0 to 1 or vice versa. Similar manipulation of binary digits is performed on SCb (not necessarily at the same position as SCt ) which ensures that both of the strings, SCt and SCb , have the same number of zeros and ones. Hence, stage 2 mutation adds or removes one bypass, achieving a nuanced operation as shown in Fig. 3(b) (right). 2.2.3.2. One-point crossover. One-point crossover operation combines two ladder networks to create a new pair of ladder networks as shown in Fig. 3(c). In contrast to the mutation operator, this operator modifies all the groups in the coded data structure. Two parent ladder networks are picked randomly from the population and two daughter networks are generated. For two parent ladder networks with the TB strings TB1 and TB2 , the daughter ladder networks’ TB strings, TB and TB are calculated as follows: TB = [TB1 (1 :  − 1) 

TB = [TB2 (1 :  − 1)

TB2 ( : p)]

(2)

TB1 ( : p)]

In the above equation, , where the crossover occurs, is a random integer generated from 1 to length of the TB group (P). For a  value of 2 an example of obtaining TB groups is shown in Fig. 3(c). The ‘TB’ group of the parent structures are ‘1 1’ and ‘1 0’, which are interchanged to ‘1 0’ and ‘1 1’. The daughter networks’ real values V , V are generated by weighted mean average of the parent networks’ real values V1 , V2 , and are given as follows: V  = V1 + (1 − )V2 V  = (1 − )V1 + V2

(3)

417

In the above equation,  is a random number between 0 and 1. The choice of segment groups are randomly assigned to match the total number of bypasses as given in the ‘TB’ string. The networks resulting from one point crossover are shown in Fig. 3(c). The new off-springs are tested for the minimum clearance constraint. If the constraint is violated, a new  is generated to obtain different daughter networks. The ladder networks generated through one point cross-over will have the properties of both their parents. Notice that the one point crossover operation increases the overall population by two units at each iteration. 2.2.3.3. Selection. The objective of selection is to terminate nonelitist structures and populate the elite structures for the next generation. Selection operation is performed only after implementing several mutations and one point cross-overs, which results in the overall population increasing to 2M. The objective of selecting elitist and diverse structures is achieved by using two metrics: fitness (F) and diversity (D). Fitness quantifies the structure with respect to the given objective function, whereas diversity quantifies the similarity between a particular network with the rest of the population. (i) Fitness (F): The use of the fitness metric is to quantify the ladder networks based on the given objective function. Networks with highest fitness values are chosen from the population for the next generation. This helps in developing off-springs which have the desired properties. The fitness function for achieving a particular droplet spacing (d) at the exit with minimum number of bypasses is given as follows:

⎛

F1 = ⎝

xu  xl

⎞

|Xout,i − d| + ˛

n−1 ⎠ N

(4)

F = e−ˇF1

In Eq. (4), ˇ and ˛ are normalizing factors and n, N are the total number of bypasses and the maximum number of bypasses, respectively. Eq. (4) represents a fitness function with multiple   objectives. The first term in the exponent

xu xl

|Xout,i − d|

is a summation of error in outlet spacing, for a range of input spacings, with a fixed delay d. The second term ((n − 1)/N) is the normalized value of the total number of bypasses. Therefore, a higher fitness value will have less bypasses and achieve a fixed delay (d). Since the spacing is in microns, ˇ is used as normalization factor and ˛ is determined by trial and error, so that both the terms have similar magnitude. Therefore, reducing the first term results in ladder networks that transform the exit drop spacing to a particular delay ‘d’ for a range of input spacings, whereas the second term in the expression reduces the number of bypasses. (ii) Diversity (D): The objective of diversity is to maintain diverse population at every generation. This helps the GA search in two ways: first, the final population will have several solutions for a given objective function and second, it minimizes the possibility of the GA search converging to one local optimum. The metric to estimate the diversity between two ladder networks is defined as follows: Ti (j) = ||SC t (i) − SC t (j)||∞ + ||SC b (i) − SC b (j)||∞ +||V (i) − V (j)||∞ ∀j = 1 : 2M & j = / i

(5)

D(i) = minimum{Ti } ∀i = 1 : 2M In the above equation, Ti is a vector, with each element (j) quantifying the difference between the ith and jth structures. The diversity is quantified by the differences in the choice of the segment groups, SCt and SCb , along with the differences in

418

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Fig. 4. Average and maximum fitness as a function of generations for achieving synchronization of a fixed inlet droplet spacing, inset contains a plot for first 100 generations.

the actual position of the nodes, V. In the equation, the binary strings, SCt and SCb , and the real values V are considered as column vectors. For the ith structure, the diversity metric, D(i), is defined as the minimum element of the Ti vector using the weakest link in the chain principle. Therefore, a zero value of D(i) implies that an exactly same network exists in the population. Notice that instead of the ∞-norm, other norms such as the 1-norm could also have been used. The next generation population is obtained by choosing networks from three groups. The first 10% are from the elitist group of the previous population based on fitness, i.e., the networks having maximum fitness value irrespective of diversity. Then 80% are chosen based on diversity and fitness, i.e., picking fittest structures that possess a diversity D > 0.03, and the last 10% are randomly assigned. If the population is less than 80% diverse, the number of random structures included in the next generation is increased to maintain the overall population of M for the next generation. 3. Results We now present the results of the proposed GA approach for three prototypical case studies. In the first study, the aim is to obtain ladder network designs that synchronize a pair of widely spaced drops with a minimum number of bypasses. In the second case study, the objective is to design ladder networks that synchronize drops with a range of input spacings ( ± ), where () is the mean of droplet spacing with a variance of . In the final case study, the objective of GA is to uncover ladder networks that transform droplets coming in at a range of input spacings ( ± ) to a constant outlet spacing (d = 2). The three case studies are designed to provide answers to the questions that were discussed at the end of the introduction section in reference to the ladder networks. To recap, the questions were: Are there simple ladder architectures that can achieve synchronization? What degree of experimental uncertainties can be tolerated while still achieving acceptable synchronizer functionality? Is it possible to design a ladder network for multiple functions such as input delay dependent behavior? Interestingly, in the literature, ladder network designs have been generated only for synchronization. It has been shown that the backward slanted bypass is the optimal structure for achieving synchronization of widely spaced drops (Maddala et al., 2012). Ladder designs have not been developed for the two other questions that we seek to answer. In fact, it is not even clear if these two design problems have reasonable solutions. We will show that the proposed GA approach provides insights for these design problems.

3.1. Ladder networks to achieve synchronization for fixed delay input In the first case study, the GA approach is formulated to design ladder networks that achieve synchronization for pairs of droplets entering at a fixed delay of 1400 ␮m. The fitness function that reflects this objective is formulated as follows: F1 = |Xout | +

˛(n − 1) N

(6)

F = e−ˇF1 The first term in the exponent (Xout ) is the error in synchronization, whereas the second term ((n − 1)/N) penalizes the number of bypasses in the structure. ˇ and ˛ are normalizing factors. As the spacing is in microns, ˇ is chosen to be 104 . Fig. 4 shows the evolution of the ladder networks’ fitness as a function of generations. Starting with an initial population of 100 random structures, the algorithm iteratively updates the population until 10 K generations. Fig. 4 shows the maximum and the average fitness of ladder designs as a function of generations. It is observed that the maximum fitness and the average fitness reach a constant value after the 350th generation. GA was terminated after 10 K generations. The first major change in the evolution of ladder networks occurred at the 10th generation as shown in the inset of Fig. 4. After the 100th generation, the fitness remains constant at 0.992 for a while. The effect of the local mutation started after the 200th generation is observed and at the 280th generation, the fitness has changed to 0.999. It is observed that there are no substantial changes in the elitist structures of the population after the 300th generation. In the previous paragraph, we discussed the evolution of ladder networks’ fitness as a function of generations. We will now describe the structural changes that result in the ladder networks during the fitness evolution discussed above. In the first generation, the structure with 5 bypasses containing a combination of forward and backward slants has the highest fitness. At the 3rd generation, a single bypass structure with slight backward slanted alignment is the most elite structure. A ladder network with two bypasses had the highest fitness in the 4th generation. The two bypass structure was the most fit till the 7th generation. At the 8th generation, this structure is replaced by another single backward slanted bypass structure. After the 8th generation till the last generation, the single backward slanted bypass ladder network remained the most fit with minor modifications in the alignment of the bypass due to local mutation. Finally, the optimal structure with the best alignment

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

419

Fig. 5. Synchronization metric and diversity metric for MLNs in the final population arranged in ascending order of fitness value (index of 100 is most fit and index of 1,as seen in the inset, is least fit).

for this objective emerged in the 9367th generation. Although the single backward slanted bypass was found in 8th generation, the rest of the generations helped in obtaining several other structures along with enhancing the backward slant to the optimal slope. The final ladder structures identified by the GA approach at the end of 10 K generations arranged in ascending order based on their fitness values (F) are shown in Fig. 5. This figure has two y-axes; diversity metric on the right and synchronization error on the left. The index in the figure starts with 10 because the first 10 structures are random and they have very less fitness value. In the plot, the index 10 ladder structure has the lowest F value whereas the index 100 structure has the highest F value. A value of synchronization error (ı) close to zero is the best synchronization and a non-zero value of D implies that the structure is unique in the population. The last 10 structures in the plot have zero D value because they are all chosen based on fitness without considering diversity. Notice that even though the overall fitness value of index 76 structure is less compared to the other higher indices, it has a less ı value (the first term in the objective) compared to the structures with indices of 77–80. As the fitness function has multiple objectives, the best structures need not have the minimum ı value. The index 76 structure achieves better synchronization with 3 bypasses, whereas the higher index structures achieve synchronization with less bypasses. This shows the ability of GA to uncover possible diverse combinations of ladder networks that achieve synchronization for a constant delay. Fig. 6 depicts the diverse ladder networks that achieve timing synchronization. The structure index is shown along with the normalized nodal positions. The 100th structure is the best as

it can achieve synchronization with only one bypass. Maddala et al. (2012) have shown that a single backward slant is the optimum structure to synchronize droplets. Therefore, the proposed GA achieved the best solution for this objective function. The 85th structure achieves synchronization with two bypasses. Even though both these structures (85 and 100) have a backward slanted bypass, the change in the optimal slope of the backward slanted bypass is compensated by the addition of another bypass in the index 85 structure. Of all the two bypass structures that are shown in the figure, the 87th structure synchronizes better than ladder networks 85 and 86. The output drop spacing for the 100th structure with an inlet drop spacing of 1400 ␮m is 0.0 ␮m and for the structures 85, 86, and 87 the outlet spacing of droplets are 1.5, 0.9, and 0.5 ␮m, respectively. Previous studies have shown that ladder network with 10 vertical bypasses reduce the droplet spacing from 160 ␮m to 0.1 ␮m as they exit the network (Prakash & Gershenfeld, 2007; Song et al., 2012). In this work, through the use of multiple objectives, the GA approach revealed several structures that achieve better synchronization (starting from a larger input delay of 1400 ␮m compared to 160 ␮m) when compared to ladder networks identified through laborious experimental trial-and-error procedures. Another interesting outcome of generating multiple solutions is that a detailed understanding on the impact of the number and orientation of bypasses on the functionality of the ladder network is obtained. In the case of the index 100 structure, a ladder network with a single bypass achieved synchronization of the droplets. So, we understand that in the case of single bypass ladder networks, a backward slanted structure is optimum for synchronization.

Fig. 6. Few solutions form the final population produced by genetic algorithm for achieving synchronization of a fixed inlet droplet spacing, the fitness rank is shown for these structures.

420

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Similarly, one might ask as to what are desirable structures for synchronization with two bypasses? Structures 85, 86, 87 answer this question. Structures 85 and 86 have backward slanted bypasses followed by a forward and a backward slanted bypass, respectively. In the case of structure 87, the backward slanted bypass is followed by a vertical bypass. Since higher index structures have higher fitness, we observe that in the case of ladder networks with two bypasses, one desirable structure is a backward slanted network followed by a vertical bypass. Incorporation of the diversity measure in the selection process is of critical importance in discovering such variegated structures. 3.2. Ladder networks to achieve synchronization for a range of input delays

Fig. 7. Normalized average and maximum fitness as function of generations for designing ladder networks that achieve synchronization for a range of input delays.

In the second case study, the genetic algorithm is used to design ladder networks that achieve synchronization for a range of input spacings. Specifically, our objective is to design ladder networks with minimum bypasses that can synchronize pairs of droplets coming at an input spacing of 1050 ± 350 ␮m. The fitness function to achieve the proposed design objective is given below: F1 =

xu  xl

|Xout,i | +

˛(n − 1) N

(7)

F = e−ˇF1 The algorithm is simulated for 15 K generations. It was noticed that the maximum fitness value was a constant value from 600 to 2 K generations and increases substantially at 2 K and 12 K generations. In contrast to the previous case, this problem required several generations for the population fitness to increase. This is due to the complexity of the design requirement. Achieving an almost constant outlet spacing for a range of input spacings using a MLN is challenging because the ladder networks have one-to-one mapping (Maddala et al., 2012) between input delay and output delay. As a result, the fitness function can never become one; it can go close to one when the variability in the outlet spacing is decreased. The fitness value shown in Fig. 7 is normalized to the maximum fitness of the population. Similar to the previous case, Fig. 8 shows the synchronization metric and diversity metric of the MLNs in the final population that are arranged in ascending order based on their overall fitness. The plot has double y-axes, the left axis shows the synchronization metric ı =

1 T

T Xout,i i Xin,i , and the right axis shows the diversity metric.

Ideally, ı should be zero but because the ladder devices transform the input spacing bijectively (Maddala et al., 2012), the ı value is always non-zero. Similar to previous case, the first 10 networks

are randomly designed at every generation, whereas the last 10 networks are the best fit networks that are generated over several populations. The structure with index 83 has a ı value less than the best fit structure because it could achieve better synchronization by increasing the number of bypasses. Fig. 9 shows two structures from the pool of the final population after 15 K iterations. The 100th index structure has 10 bypasses with a combination of forward and backward slanted bypasses. The exit spacing of drops is −77 ␮m and 79 ␮m for an inlet spacing 700 ␮m and 1400 ␮m, respectively. The negative value of the exit spacing implies that the lagging droplet at the entrance has overtaken the leading droplet. This ladder structure reduces the variation of inlet drop spacing of 700–156 ␮m at the exit. This shows that cleverly designed ladder networks with only 10 bypasses could synchronize droplets entering at variable inlet spacings. The network not only synchronizes droplets but also reduces the droplet spacing variance by 78%. Previous studies (Prakash & Gershenfeld, 2007; Song et al., 2012) show that MLNs with 10 vertical bypasses achieve synchronization of drops coming at 190 ± 60 ␮m. The exit spacing of these droplets is reduced by 93 ± 5%. In the present ladder design, the droplets entering an initial spacing of 1050 ± 350 ␮m is reduced to 96 ± 2%. Further, if these bypasses are made vertical, the exit spacing is reduced to only 52 ± 20%. This shows the importance of the orientation of bypasses to achieve synchronization for a range of input spacings. In the case of the structure 83, the exit spacing of droplets is −71 ␮m and 81 ␮m for an inlet spacing of 700 ␮m and 1400 ␮m, respectively. It is observed that due to increase in three bypasses from the 100th index structure, the variance has reduced, i.e., from 156 ␮m to 152 ␮m. Further, the GA, approach identifies 80 different

Fig. 8. Synchronization metric and diversity metric for the final population arranged in ascending order of fitness.

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

421

Fig. 9. Two structures from the final population that have the best fitness and best synchronization metrics.

structures that could achieve the task of synchronizing droplets coming at a range of input spacings. Other interesting insights can be gleaned from the structures that are identified by the GA approach. Though the objective is to synchronize droplets, the droplet spacing is not reduced by all the bypasses in an optimal structure. For example, the first 4 bypasses of the ladder networks expand the droplet spacing while the remaining bypasses reduce the spacing to zero as shown in Fig. 10. Fig. 10 shows the relative distance between the droplets as they travel through the nodes of the ladder network. As flows are calculated by the network model only when the droplets are at the nodes, the position of the droplets are captured at those instances. There are 21 nodal positions plotted on the x-axis, because of 10 bypasses and one entrance value. It is observed that in contrast to synchronizer designs studied in literature (Prakash & Gershenfeld, 2007; Song et al., 2012) the manner in which asymmetric networks achieve synchronization is quite different. The reason is that the asymmetric bypasses change flows at different location on the transport channels compared to vertical bypass ladder networks. The objective of GA is to not only reduce the mean but also the variance; therefore, to reduce the variance in the inlet, the first few bypasses in the ladder networks act as an amplifier and the latter bypasses act as a filter as shown in Fig. 10. Thus the first few bypasses amplify the signal to a particular value where the signal can be completely filtered out by the rest of the bypasses. This is an outcome of complex flow variations that occur in these bypasses as they vary due to droplet motion in the ladder networks. This approach to synchronization of a range of input delays was never

conceived by us at the start of these computational experiments but is a result of the analysis of the ladder networks uncovered by the proposed GA approach. 3.3. Ladder network to achieve a fixed delay for a range of input spacing In the third case study, the objective is to design ladder networks which produce a fixed droplet spacing at the outlet for a range of input spacings. The optimal ladder networks should increase the distance between pairs at the outlet. The fitness function to achieve this objective is given below: F1 = F=

xu 

|Xout,i − d| +

xl e−ˇF1

˛(n − 1) N

(8)

In the above equation, the value of d is 2200 ␮m. The simulation is terminated after 25 K generations. Similar to the previous case studies, the maximum fitness increases at discrete intervals as shown in Fig. 11. Interestingly, in this case, the average fitness changed before the maximum fitness. Fig. 12 shows the diversity and the delay metric, ı, where ı =

1 T

T Xout,i −d i Xin,i ; similar to the previous

cases there are structures in between 10th and 90th indices, which have a better delay metric compared to the 100th structure. Fig. 13 shows the ladder designs of indices 100 and 76. The design with index 100 expands the droplets to a specific delay d,

Fig. 10. Transformation of inlet spacing as a pair of droplets move in the ladder network.

422

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Fig. 11. Normalized average and maximum fitness as a function of generations for achieving a fixed delay at the outlet for a range of input spacings.

with minimum number of bypasses. The best structure does not have the least delay metric (ı), the 76th index ladder design has a delay metric that is less. To understand the phenomena of expansion, consider structure with index 100. A pair of droplets enter the ladder network with 700 ␮m inter-droplet spacing. The bottom droplet enters the ladder network and slows down whereas the top droplet comes with the initial speed expanding the relative spacing between the droplets. By the time the top droplet comes to the first node, the bottom droplet is caught in several bypasses such that the velocity of the bottom droplet is always lower than the top droplet. Therefore, by the time the top droplet exits the ladder the spacing is increased to 2000 ␮m. That is an increase of three

times the inlet spacing. In case of 1400 ␮m, the spacing increases to 2400 ␮m which is an increase of 1000 ␮m. The objective was to increase drops coming in between 700–1400 ␮m and 2200 ␮m. Therefore a variation of 700 ␮m at the inlet is reduced to a variation of 400 ␮m at the exit, that is a 42% reduction in variation as shown in Fig. 14. Every bypass in the ladder network is important to achieve this particular transformation but the extent to which they affect the functionality varies. For example, in case of ladder network 100, removing the last bypass changes the functionality by 5% whereas removing the first bypass changes the functionality by 20%. In case of ladder network with index 76 at the lower limit of the inlet spacing, the exit spacing is 2100 ␮m, which is 100 ␮m higher than the 100 structure. At 1400 ␮m, the exit spacing is 2300 ␮m, that is a decrease of 100 ␮m from the 100th structure. In the case of structure 76, the exit variation is just 200 ␮m for a 700 ␮m variation in the inlet, that is a reduction of 71% in variation. Interestingly, in this case, similar to the previous example, the first bypass reduces the inlet spacing and the other bypasses increase the spacing. Further, due to the diversity factor there are 80 different MLNs that perform reasonably well for the same objective. 3.4. Sensitivity analysis The objective function for synchronizing drops for fixed delay (cf. Section 3.1) is used to evaluate the influence of some of the genetic parameters and operators on the performance of the genetic algorithm. The sensitivity of this algorithm to variations

Fig. 12. Synchronization metric and diversity metric for MLNs arranged in ascending order of fitness value.

Fig. 13. Two structures from the final population with best fitness (100) and best synchronization metric (76).

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

423

Fig. 14. Transformation of inlet spacing as a pair of droplets move in the ladder network.

in the population size, diversity and the multi-objective weighting factor are presented. Additionally, to understand the effect of novel encoding on the final solutions, the solutions obtained by genetic algorithm are compared with that of a random search algorithm. As GAs are based on stochastic rules, we would like to understand how the performance of the algorithm changes for different trials with the same objective function. Therefore, to analyze the performance of the proposed algorithm, 10 trials are carried out with the same objective function. After 104 iterations the results are compared and we observed that in all these attempts, the proposed GA approach obtained a single bypass ladder network as the optimal solution with a difference of 0.001% in the slope of the bypass, which led to an insignificant change in the fitness function. 3.4.1. Population size To study the effect of population size on the performance of the GA, simulations are performed varying the population size in the range of 10–100 networks. Each simulation is terminated after 1000 generations. Fig. 15 shows the number of iterations required to reach the maximum fitness as a function of the population size. It is observed that there is no defined pattern on the performance of the GA with respect to population size. A population size of 100 could reach a fitness value of 1 in 360 iterations, whereas when the algorithm is simulated with other population sizes the algorithm

Fig. 15. Influence of the population size on the performance of genetic algorithm.

never reached a fitness value of 1. Hence a population size of 100 is used in the simulations. 3.4.2. Effect of genetic operations In the proposed algorithm, either a mutation or a cross-over operation is performed on the population at every iteration with an equal probability. The effect of mutation and cross-over on the performance of the GA are analyzed by varying the percentages for these operators at every iteration. The mutation operator percentage is varied from 10% to 100%, which changes the cross-over percentage from 90% to 0%. Fig. 16 shows the iterations required for the GA to reach a steady fitness value. It is observed that for equal probability of mutation and cross-over the genetic algorithm attains the best fitness using a minimum number of iterations. Additionally it is observed that for 100% mutation the algorithm has the least fitness value. 3.4.3. Effect of encoding The proposed encoding scheme efficiently handles constraints on the network design as well as potentially span all the possible designs of the ladder networks. The solutions obtained by GA are compared with a random search using the proposed encoding scheme. This explores if the best designs obtained by GA are a result of the encoding scheme or due to a combination of the encoding scheme with genetic operators. As the GA is simulated for 104 generations, a random search is performed by generating an equal number of networks used in GA (≈106 ). The ratio of the

Fig. 16. Influence of genetic operations on the performance of genetic algorithm.

424

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425

Table 1 Influence of random population on the performance of genetic algorithm. % elite networks

% random networks

Generation number

Fitness

0 20 40 80 100

100 80 60 20 0

500 181 471 86 500

0.5067 0.9948 0.9996 0.9997 0.7173

Table 2 Influence of ˛ on the performance of genetic algorithm. % no. of branches

% weighing factor (˛)

Iteration number

Fitness

6 8 1 1 1

0 10−6 10−5 10−4 10−3

52 283 257 86 500

0.9990 0.9953 0.997 0.9998 0.9988

synchronization error between the solution obtained by GA to that of random search is 0.00075. This implies that the synchronization error in the case of GA is substantially less compared to that of the random search. Therefore, optimal solutions are identified not only because of novel encoding but also because of the genetic operations. 3.4.4. Effect of random structures in the population At every generation, the reproduction operator picks 90 structures from a pool of 200 structures and combines this with 10 random structures to increase the diversity in the population. Several combinations of random structures to elite structures are simulated to study the effect of random structures on the overall performance of GA. The GA is terminated after 500 generations and the final fitness and the generation number are presented in the Table 1. It is observed that a population with only random structures or only elite structures has the least performance in terms of fitness. We also observed that a combination of elite and random structures with a bigger percentage of elite structures seem to give better performance. Hence we choose 90% elite and 10% random structures at every generation. 3.4.5. Effect of multi-objective weighing parameter The normalizing factor ˛ in the objective function (Eq. (4)) is varied to study its effect on the overall performance of the GA. The parameter ˛ normalizes the error in bypasses to that of the synchronization error. This value is varied in the range of 0 to 10−3 as shown in Table 2. Similar to the previous case, the algorithm is terminated after 500 iterations and the fitness values along with the number of bypasses present in the optimal design are reported in Table 2. It is observed that when ˛ value is zero the optimal solution has 6 bypasses but as the value of ˛ is increased the optimal ladder design proposed by the GA has only one bypass. We chose an ˛ value of 10−4 as this resulted in the best solution in a minimum number of iterations. 4. Conclusions A GA approach was formulated for designing microfluidic ladder networks. Using this approach, three different ladder network design problems were studied. The first design problem was one of synchronization of exit timings of the droplets. The second problem addressed the reduction of variance along with synchronization when the input droplets arrive at a range of input delays instead of a fixed delay as in the first design problem. The third design problem was to have the droplets leave at a fixed delay. It was shown

that the GA approach identified remarkable and diverse structures for each one of these design problems. In order to evaluate the performance of the algorithm, different variations of the proposed genetic operators are examined along with the sensitivity analysis. While the objective functions considered in this work are based on the exit spacing, the proposed approach can be used for any generic ladder structures that involve more complex spatial and temporal control of droplets. Further, in this paper we considered only pairs of droplets; however, this work could be easily extended to trains of droplets. As microfluidic ladder designs are widely observed in nature (Chau, Rolfe, & Cooper-White, 2011; Hosoya, Okado, Sugiura, & Kohno, 1991; Roth-Nebelsick, Uhl, Mosbrugger, & Kerp, 2001; Skalak, Özkaya, & Skalak, 1989), the proposed approach could be used to understand natural systems better. Our future work will focus on generalizing the proposed approach to any complex microfluidic network design. References Boukellal, H., Selimovic, S., Jia, Y., Cristobal, G., & Fraden, S. (2008). Simple, robust storage of drops and fluids in a microfluidic device. Lab on a Chip, 9, 331–338. http://dx.doi.org/10.1039/b808579j Bruus, H. (Ed.). (2008). Theoretical microfluidics. New York: Oxford University Press, Inc. Chau, L. T., Rolfe, B. E., & Cooper-White, J. J. (2011). A microdevice for the creation of patent, three-dimensional endothelial cell-based microcirculatory networks. Biomicrofluidics, 5, 034115(1)–034115(14). Christopher, G. F., Bergstein, J., End, N. B., Poon, M., Nguyen, C., & Anna, S. L. (2009). Coalscence and splitting of confined droplets at microfluidic junctions. Lab on a Chip, 9, 1102–1109. http://dx.doi.org/10.1039/b813062k Cristobal, G., Benoit, J. P., Jaonicot, M., & Ajdari, A. (2006). Microfluidic bypass for efficient passive regulation of droplet traffic at a junction. Applied Physics Letters, 89(034104), 1–3. http://dx.doi.org/10.1063/1.2221929 Cybulski, O., & Garstecki, P. (2009). Dynamic memory in a microfluidic system of droplets traveling through a simple network of microchannels. Lab on a Chip, 10, 484–493. http://dx.doi.org/10.1039/b912988j Epstein, I. R. (2007). Can droplets and bubbles think? Science, 315(5813), 775–776. Garrard, A., & Fraga, E. (1998). Mass exchange network synthesis using genetic algorithms. Computers & Chemical Engineering, 22(12), 1837–1850. Goldberg, D. (1989). Genetic algorithms in search, optimization, and machine learning. Reading, Massachusetts: Addison-Wesley. Hatakeyama, T., Chen, D. L., & Ismagilov, R. F. (2006). Microgram-scale testing of reaction conditions in solution using nanoliter plugs in microfluidics with detection by MALDI-MS. Journal of the American Chemical Society, 128, 2518–2519. http://dx.doi.org/10.1021/ja057720w Ho, T. Y., Chakrabarty, K., & Pop, P. (2011). Digital microfluidic biochips: Recent research and emerging challenges. In Hardware/Software Codesign and System Synthesis (CODES + ISSS), 2011 Proceedings of the 9th International Conference on (1st ed., pp. 335–343). Holland, J. (1975). Adaptation in natural and artificial systems. Cambridge Massachusetts: The MIT Press. Hosoya, Y., Okado, N., Sugiura, Y., & Kohno, K. (1991). Coincidence of “ladderlike patterns” in distributions of monoaminergic terminals and sympathetic preganglionic neurons in the rat spinal cord. Experimental Brain Research, 86, 224–228. Jeevan, M., Babji, S., Swastika, B, Vanapalli, S. A., & Rengaswamy, R. (2012). Design of a model-based feedback controller for active sorting and synchronization of droplets in a microfluidic loop. AIChE Journal, 58(7), 2120–2130. http://dx.doi.org/10.1002/aic.12740 Jin, B. J., Kim, Y. W., Lee, Y., & Yoo, J. Y. (2010). Droplet merging in a straight microchannel using droplet size or viscosity difference. Journal of Micromechanics and Microengineering, 20(035003), 1–3. http://dx.doi.org/ 10.1088/0960-1317/20/3/035003 Laval, P., Crombez, A., & Salmon, J. B. (2009). Microfluidic droplet method for nucleation kinetics measurements. Langmuir, 25(3), 1836–1841. http://dx.doi. org/10.1021/la802695r Maddala, J., Wang, W., Vanapalli, S., & Rengaswamy, R. (2012). Traffic of pairs of drops in microfluidic ladder networks with fore-aft structural asymmetry. Microfluidics and Nanofluidics, 1–8. http://dx.doi.org/10.1007/s10404-012-1054-z Parker, R. G., & Rardin, R. L. (1988). Discrete optimization. Reading, Massachusetts: Academic Press. Prakash, M., & Gershenfeld, N. (2007). Microfluidic bubble logic. Science, 315, 832–835. http://dx.doi.org/10.1126/science.1136907 Ravagnani, M., Silva, A., Arroyo, P., & Constantino, A. (2005). Heat exchanger network synthesis and optimisation using genetic algorithm. Applied Thermal Engineering, 25(7), 1003–1017. Roth-Nebelsick, A., Uhl, D., Mosbrugger, V., & Kerp, H. (2001). Evolution and function of leaf venation architecture: A review. Annals of Botany, 87, 553–556.

J. Maddala, R. Rengaswamy / Computers and Chemical Engineering 60 (2014) 413–425 Schindler, M., & Ajdari, A. (2008). Droplet traffic microfluidic networks: A simple model for understanding and designing. Physical Review Letters, 100(044501), 1–4. http://dx.doi.org/10.1103/PhysRevLett.100.044501 Sessoms, D. A., Amon, A., Courbin, L., & Panizza, P. (2010). Complex dynamics of droplet traffic in a bifurcating microfluidic channel: Periodicity, multistability, and selection rules. Physical Review Letters, 105, 154501. Shim, J., Cristobal, G., Link, D. R., Thorsen, T., Jia, Y., Piattelli, K., et al. (2007). Control and measurement of the phase behavior of aqueous solutions using microfluidics. Journal of the American Chemical Society, 129, 8825–8835. http://dx.doi.org/10.1021/ja071820f

425

Skalak, R., Özkaya, N., & Skalak, T. C. (1989). Biofluid mechanics. Annual Review of Fluid Mechanics, 21, 167–204. Song, K., Zhang, L., & Hu, G. (2012). Modeling of droplet traffic in interconnected microfluidic ladder devices. Electrophoresis, 33(3), 411–418. Su, F., Chakrabarty, K., & Fair, R. (2006). Microfluidics-based biochips: Technology issues, implementation platforms, and design-automation challenges. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 25(2), 211–223. http://dx.doi.org/10.1109/TCAD.2005.855956 Whitley, D. (1994). A genetic algorithm tutorial. Statistics and Computing, 4(2), 65–85.