Tuning-free controller to accurately regulate flow rates in a microfluidic network

We describe a control algorithm that can improve accuracy and stability of flow regulation in a microfluidic network that uses a conventional pressure pump system. The algorithm enables simultaneous and independent control of fluid flows in multiple micro-channels of a microfluidic network, but does not require any model parameters or tuning process. We investigate robustness and optimality of the proposed control algorithm and those are verified by simulations and experiments. In addition, the control algorithm is compared with a conventional PID controller to show that the proposed control algorithm resolves critical problems induced by the PID control. The capability of the control algorithm can be used not only in high-precision flow regulation in the presence of disturbance, but in some useful functions for lab-on-a-chip devices such as regulation of volumetric flow rate, interface position control of two laminar flows, valveless flow switching, droplet generation and particle manipulation. We demonstrate those functions and also suggest further potential biological applications which can be accomplished by the proposed control framework.


Text S.1 8 channel network for a case study
In this section, we introduce a more complex microfluidic network ( Fig. S.1) than the Y-junction network presented in the paper. This case study that uses the complex network will clarify 1) construction of the steady-state gain matrix, 2) selection of controllable channels, 3) implementation of the proposed controller, and 4) stability of the closed-loop system.

Text S.1.1 Construction of the steady-state gain matrix
To construct a steady-state gain, the resistance network ( Fig. S.2) should be used. The network property is as follows:   Using (10), A 1 and A 2 are obtained as follows: p 1 − p n,1 R 1 = p n,1 − p n,4 R 2 + p n,1 − p n,2 R 3 (S.1) Q 2 = Q 4 + Q 5 =⇒ p n,1 − p n,4 R 2 = p n,4 − p 4 R 4 + p n,4 − p n,3 R 5 (S.2) Q 3 = Q 6 + Q 7 =⇒ p n,1 − p n,2 R 3 = p n,2 − p n,3 R 6 + p n,2 − p 2 R 7 (S.3) From (11), A 3 and A 4 also can be obtained as ∆p 1 = p 1 − p n,1 ∆p 2 = p n,1 − p n,4 ∆p 3 = p n,1 − p n,2 ∆p 4 = p n,4 − p 4 ∆p 5 = p n,4 − p n,3 ∆p 6 = p n,2 − p n,3 The conductance matrix is B = diag(1/R 1 , · · · , 1/R 8 ) and can be represented as Using (14) and substituting the nominal resistance values, the steady-state gain can be obtained as (S.8) Text S.1.2 Simultaneously controllable channels and existence of a right inverse of the steady-state gain matrix To prove the existence of a right inverse of G ss , we define controllable channels in a microfluidic network. All networks (or undirected graphs) that have m branches (edges), n nodes (vertices) and l closed loops have a topology relation m = n + l − 1. If a microfluidic network is composed of N c channels, N n nodes, N a input ports and N l closed loops where vertices are nodes and input ports (i.e., n = N n + N a ), then the topology relation yields N c = N n + N a + N l − 1. This can be rewritten as where π is a rank of G ss (i.e., π = rank(G ss )). The rank represents 'the maximum number of simultaneously controllable channels' in a microfluidic network. In the Y-junction network, π = 3−1 = 2, so flow in an arbitrary two of the three channels in the network can be simultaneously controlled, and flow in the other channel is automatically determined by mass conservation (e.g., Q 1 = Q 2 + Q 3 in the Y-junction network). The number of controllable channels is two, so three input sources (one inlet and two outlets (Figs. 2, 3, 6), or two inlets and one outlet (Figs. 4, 5)) are required. Here, an outlet that is not selected as a controllable channel can be open to the atmosphere (Fig. 4, 5). Therefore, inlet and outlet should be selected appropriately depending on the purpose of the application. In matrix G ss , the j-th row concerns the j-th output flow, and the rows of G ss must be reduced to < π because the number of controllable channels cannot exceed π. Here, if π < N a , then the reduced steady-state gain G ss has a full-row rank and therefore always has a right inverse matrix. If π = N a , then G ss is a square and invertible matrix.
In the case study, the maximum number of simultaneously controllable channels is Therefore, three channels of the network can be simultaneously controllable and the matrix G ss always has a right inverse matrix (∵ π = 3 < 4 = N a ).
If we would like to simultaneously control two of the three controllable channels, then π = 2 = N a − 1 and N a = 3, (S.11) so three pressure sources are required.

Text S.1.3 Selection of controllable channels
We can arbitrarily select controllable channels by selectively eliminating rows of G ss . If we want to control Q 2 , Q 3 and Q 5 , then the 1-st, 4-th, 6-th, 7-th and 8-th rows of G ss should be removed as where the flow rates inside channel 2, 3, 5 can be independently controlled while the others will be determined by their dependency on the controlled flows. Any three channels in the network can be simultaneously controlled by this method.
By eliminating rows that will not be controlled, G ss becomes a full row rank matrix that always has a right inverse matrix. For this case, the right inverse of G ss is 6.0 5.5 0.5 6.0 −11.5 5.5 −6.0 5.5 −11.5 −6.0 0.5 5.5     . PressureRegulator(u c ); 15: end while; measured flow rate is allocated to y m . At line 6, control input is computed as presented in (4). From line 7 to line 12, the control input is evaluated by the saturation function in (6). At line 14, the computed control input is sent to the pressure regulator to insert control input pressure to the microfluidic chip. where the transfer matrix function is constructed by synthesizing eight T-models as presented in the paper. The transfer matrix function of the proposed controller is By constructing a feedback loop, which consists of the controller C(s) in negative feedback with the plant G(s), the transfer matrix function of the closed-loop system is obtained. If poles of the closed-loop system are all negative real, then the closed-loop system is internally stable. Moreover, if the system is internally stable, then the origin of state is asymptotically stable [ZDG + ]. Here, ∆t scales the feedback gain in the controller. Thus, ∆t should be bounded in a gain margin of the system. Near the origin, as ∆t decreases, the poles go to the right half plane and the closed-loop becomes unstable (Fig. S.4a). In our case study, the closed-loop system becomes unstable when ∆t < 0.01s in the discrete implementation. The closed-loop stability of the discrete system can be directly evaluated. The discrete state-space model of a microfluidic network system is where A, B and C are constant matrices with appropriate dimensions; the model obtained by minimal realization of G(s). v and w are unknown constant vectors that can represent constant model uncertainties or slowly-varying process disturbances [DR12]. By substituting the control input into (S.16) and subtracting the previous time-step equation from (S.16) and (S.17), the model can be rewritten as Therefore, closed-loop dynamics to which the proposed controller is applied can be represented by the augmented state model as (S.20) If all eigenvalues ofÃ are < 1, then the closed-loop system is asymptotically stable [Che95] (Fig.  S.4b).
In conclusion, when the proposed robust controller is applied, the output flow rate converges to the desired output flow rate.

Text S.2 Tuning of transient responses using a single parameter
The transient response (not steady-state response) can be tuned by using a single scalable parameter α to alter the original gain matrix G † ss as follows: The closed-loop transient response increases as α increases, and slows as α decreases. This tuning parameter can be used if fast response is required. Also, as we mentioned in Text S.1.6, the closedloop response can be unstable in the fast discrete implementation. For this case, α can be used to adjust unstable poles to be stable.

Text S.3 Measurement of image intensity
To measure concentrations of a dye solution, we used image intensity computed using a selfdeveloped image-processing algorithm. The works are in three steps: 1) convert the color image to gray-scale; 2) define the centerlines of each channel that will be evaluated ( Fig. S.5); 3) compute the intensity of each line. The mean intensity was used to describe the concentration of dye in the liquid solution.

Text S.4 Analysis method for size distribution of droplets
The droplet size distribution was obtained using a self-developed image-processing algorithm. Droplets generated in the Y-junction network have the same vertical length as the channel width (w), but their horizontal length (L) depended on both flow rates. The image processing algorithm to compute the droplet size (horizontal length) proceeds as follows: 1) Binarizing ('im2bw' function in MATLAB)  Theoretical analysis of droplet size with respect to flow rates in a T-junction which is similar to the Y-junction was published previously [GFSW06]. That study identified a scaling law for droplet length (L) with respect to flow rate: L ∼ 1 + aQ oil /Q water (S.21) In our experiments (Table S.1), L increased as oil flow rate increased (Fig. S.6a, b) and decreased as water flow rate increased (Fig. S.6c). Table. S.1 summarizes all droplet experiments.