Microﬂuidic Network Simulations Enable On-Demand Prediction of Control Parameters for Operating Lab-on-a-Chip-Devices

: Reliable operation of lab-on-a-chip systems depends on user-friendly, precise, and pre-dictable ﬂuid management tailored to particular sub-tasks of the microﬂuidic process protocol and their required sample ﬂuids. Pressure-driven ﬂow control, where the sample ﬂuids are delivered to the chip from pressurized feed vessels, simpliﬁes the ﬂuid management even for multiple ﬂuids. The achieved ﬂow rates depend on the pressure settings, ﬂuid properties, and pressure-throughput characteristics of the complete microﬂuidic system composed of the chip and the interconnecting tubing. The prediction of the required pressure settings for achieving given ﬂow rates simpliﬁes the control tasks and enables opportunities for automation. In our work, we utilize a fast-running, Kirchhoff-based microﬂuidic network simulation that solves the complete microﬂuidic system for in-line prediction of the required pressure settings within less than 200 ms. The appropriateness of and beneﬁts from this approach are demonstrated as exemplary for creating multi-component laminar co-ﬂow and the creation of droplets with variable composition. Image-based methods were combined with chemometric approaches for the readout and correlation of the created multi-component ﬂow patterns with the predictions obtained from the solver. and agreed to the published of the manuscript.


Introduction
Microfluidics is the backbone of lab-on-a-chip devices [1]. Over the past decade, complex laboratory workflows have been implemented through on-chip technologies [2][3][4][5], ranging from the generation of chemical gradients [6,7], over microfluidic-based nanoparticle fabrication [8], their surface modification [9,10] to point-of-care diagnostics [11] and environmental analysis [12]. Thus, the performance requirements of microfluidic devices are increasing, and reliable control and prediction of flow rate and operating pressure are crucial. Especially pattern-based microfluidics [13] require mechanical designs tailored to their specific application, which require new and optimized approaches for the design [13][14][15] and control of microfluidic networks.
Microfluidic design automation (MDA) tools for modeling laminar flow [16,17] and droplet-based microfluidic systems have already been developed [18][19][20]. This new approach analyzes single-phase microfluidic devices that rely on the pressure-driven operation and multi-phase droplets in a continuous oil phase. A typical pressure-driven microfluidic system is the tree shape [17], with a uniform and identical cross-section over all channels, requiring proportionality between flow resistance and channel length [16]. Thus, the pressure change ∆p, volume flow rate Q, and hydrodynamic resistance R can be related to a voltage change, current, and electrical resistance, respectively [16,17]. These resistive networks predict the flow rate within the microfluidic network (MFN) to obtain

Microfluidic Chip Fabrication and Chip Preparation
The two-layer glass microfluidic chip was prepared by photolithography and wetetching of the glass wafers (thickness 0.7 mm) with a refractive index of 1.47 with hydrofluoric acid with a uniform etching depth of 100 µm mask with a nickel-chromium metallization. Two micro-structured wafers were bonded face-to-face to create the full channels and geometries for the operating unit. The chip dimensions were 16.0 × 12.5 mm. The five fluid interconnection ports arranged in the D-sub scheme were fabricated using ultrasonic drilling. The standard channel width was 140 µm, the width before the flowfocusing unit was 165 µm, and the width of the observation chamber was 200 µm. The full channel height was 100 µm with a planar channel bottom.
For droplet generation, the wetting properties of the microchannel surfaces were adjusted by treatment with PlusOne Repel-Silane ES (a solution of Dimethyldichlorosilane (2% w/v) in Octamethylcyclotetrasiloxane) from GE Healthcare Bio-Sciences AB (Upsala, Sweden). After 10 min of incubation time, the microfluidic chip was rinsed with ethanol to remove residuals. The modification of the glass surfaces increases the hydrophobicity of the channel walls resulting in droplet formation of the aqueous phase components.

Microfluidic Control
The chip was operated by four PX-1 pressure-based flow controllers from Fluigent (Paris, France), ranging from 0 to 1000 mbar. The quality of set pressure and displayed pressure is illustrated in Figure 1a. The sample and buffer solutions were contained in 2 mL Fluiwell-4C vessels with a screw cap held by the Fluigent four Channel Sampling Rack Type 4C.

Microfluidic Control
The chip was operated by four PX-1 pressure-based flow controllers from Fluigent (Paris, France), ranging from 0 to 1000 mbar. The quality of set pressure and displayed pressure is illustrated in Figure 1a. The sample and buffer solutions were contained in 2 mL Fluiwell-4C vessels with a screw cap held by the Fluigent four Channel Sampling Rack Type 4C. The pressure control unit was connected to the microfluidic chip network by PEEK HPLC capillaries with a length of 250 mm ± 5 mm with an outer diameter (OD) of 1/16". The performance of the capillaries without flow sensors is shown in Figure 2a, as well as the resulting inner diameter (ID) of 0.136 mm per channel calculated based on its throughput characteristics in Figure 2b. For control and validation of the inlet pressure settings for each channel, flow sensing units size M with an operating range of 0-80 µL/min ± 2.6 µL/min were used together with a flow sensor board from Fluigent (Paris, France). The quality of the calibrated flow sensing units carrying the buffer solution and dSurf is displayed in Figure 1b. The flow sensors were followed by another set of four PEEK HPLC capillaries with a length of 250 mm ± 5 mm, which were connected to HPLC-Teflon capillaries of length 60 mm, with an The pressure control unit was connected to the microfluidic chip network by PEEK HPLC capillaries with a length of 250 mm ± 5 mm with an outer diameter (OD) of 1/16". The performance of the capillaries without flow sensors is shown in Figure 2a, as well as the resulting inner diameter (ID) of 0.136 mm per channel calculated based on its throughput characteristics in Figure 2b.

Microfluidic Control
The chip was operated by four PX-1 pressure-based flow controllers from Fluigent (Paris, France), ranging from 0 to 1000 mbar. The quality of set pressure and displayed pressure is illustrated in Figure 1a. The sample and buffer solutions were contained in 2 mL Fluiwell-4C vessels with a screw cap held by the Fluigent four Channel Sampling Rack Type 4C. The pressure control unit was connected to the microfluidic chip network by PEEK HPLC capillaries with a length of 250 mm ± 5 mm with an outer diameter (OD) of 1/16". The performance of the capillaries without flow sensors is shown in Figure 2a, as well as the resulting inner diameter (ID) of 0.136 mm per channel calculated based on its throughput characteristics in Figure 2b. For control and validation of the inlet pressure settings for each channel, flow sensing units size M with an operating range of 0-80 µL/min ± 2.6 µL/min were used together with a flow sensor board from Fluigent (Paris, France). The quality of the calibrated flow sensing units carrying the buffer solution and dSurf is displayed in Figure 1b. The flow sensors were followed by another set of four PEEK HPLC capillaries with a length of 250 mm ± 5 mm, which were connected to HPLC-Teflon capillaries of length 60 mm, with an For control and validation of the inlet pressure settings for each channel, flow sensing units size M with an operating range of 0-80 µL/min ± 2.6 µL/min were used together with a flow sensor board from Fluigent (Paris, France). The quality of the calibrated flow sensing units carrying the buffer solution and dSurf is displayed in Figure 1b. The flow sensors were followed by another set of four PEEK HPLC capillaries with a length of 250 mm ± 5 mm, which were connected to HPLC-Teflon capillaries of length 60 mm, with an OD of 1/16" and an ID of 0.5 mm all by JASCO Deutschland GmbH (Pfungstadt, Germany) and connected through ROTILABO ® Silicone Standard sleeves of length 15 mm with an inner diameter of 1 mm from Carl Roth GmbH + Co. KG (Karlsruhe, Germany). The fluid management was controlled with the real-time software tool ALL-in-One (A-i-O) from Fluigent (Paris, France). This software allows to monitor and record the pressure status and flow rates in parallel. For the validation experiments, several measurement modes with different pressure settings were used. The graphical user interface (GUI) allows creating, editing, and loading pre-defined pressure profiles, e.g., calculated by the presented solver, for ease of use during experiments.

Optical Setup and Image Acquisition
The optical setup was a self-built microscope. A self-made LED module with a 3W 3000K LED from CREE Inc. (Durham, USA) with charcoal lighting is used as a transmission light source. The chip holder was integrated into an XY linear stage from Märzhäuser Wetzlar GmbH & Co. KG (Wetzlar, Germany). A microscopy objective from Mitutoyo Corporation (Kawasaki, Japan) with a 10× magnification, numerical aperture (NA) of 0.28, and working distance (WD) of 34 mm was used. In addition, a tube lens with a magnification of 0.5 from Opto GmbH (Gräfelfing, Germany) was used. Thus, the total magnification of the system is 5×. A color camera acA1920-155uc from Basler AG (Ahrensburg, Germany) was used to record the experiments. The validation experiments were recorded at a frame rate of 1 Hz 16-bit for 10 s per pressure profile with an exposure time of 100 µs. An additional long run of 90 s has been performed to demonstrate the stability of the flow pattern. Image processing and data analysis workflow are described in chapter 4.
For a complete measurement, the acquisition of the following reference images was mandatory: IM SAMPLE , IM BRIGHT , IM DARK , IM BLANK , and IM FILLED . For, IM FILLED all channels were uniformly filled with dye-containing buffer to measure the local optical path length d by the absorbance imaging approach [23]. In multi-colored flow, the image has to be acquired for all dye solutions available. The channel mask image IM CHANNEL_MASK was obtained from IM BLANK by thresholding to identify the pixels inside the channel.

Microfluidic Network Solver (MfnSolver)
The presented mfnSolver is a software component for MDA. Based on Kirchhoff's first law, that the sum of currents flowing into a junction or node equals the sum of the outflowing currents, and the second law describing the directed sum of potential differences in a closed loop is zero [24], the resulting system of linear equations of the conductors can be used to calculate the unknown potential differences or currents. In electronics, a conductor is characterized by a resistivity R as shown in Equation (1) (left), where the voltage U denominates the electrical potential difference, driving the transport of electrons through the conductor, also known as Ohm's law. The same can be applied in fluid-dynamics, leading to the hydrodynamic resistance, as shown in Equation (1) (right), where the conductor is a channel, the potential unit is the pressure drop ∆p, which drives the volume flow Q of volume elements through the channel, leading to the Hagen-Poiseuille equation for the hydraulic behavior [17], assuming laminar, viscous, and incompressible flow. The presented mfnSolver overcomes the limitation of uniform cross-section. The initial network representation is a Kirchhoff-graph representation of an incidence matrix with nodes and connecting edges. Afterward, the linearized transport characteristics in terms of an edge-specific resistivity/volume/length (RVL) model are assigned to the individual edges of the network and used to create a conductivity matrix. Our work functionalizes the Kirchhoff-based nodal analysis of MFN simulations for the in-line prediction of control pressures. In particular, the user defines the flow rates along arbitrary edges of the microfluidic network, and the solver calculates the required pressure settings at the IO's to achieve these flow rates. We demonstrate the operation of microfluidic devices at user-defined flow rates utilizing the flow-focusing unit PDG2 as a proof-ofconcept and reference. Finally, the mfnSolver can be utilized to design and optimize microfluidic components, predict network performance, and automate the operation of devices based on pressure-driven fluid management schemes. The mfnSolver solves a network of representative resistive components based on Kirchhoff nodal analysis as schematically illustrated in Figure 3.
Processes 2021, 9, x FOR PEER REVIEW 5 of 23 microfluidic network, and the solver calculates the required pressure settings at the IO's to achieve these flow rates. We demonstrate the operation of microfluidic devices at userdefined flow rates utilizing the flow-focusing unit PDG2 as a proof-of-concept and reference. Finally, the mfnSolver can be utilized to design and optimize microfluidic components, predict network performance, and automate the operation of devices based on pressure-driven fluid management schemes. The mfnSolver solves a network of representative resistive components based on Kirchhoff nodal analysis as schematically illustrated in Figure 3.

MfnSolver Concept and Implementation
The mfnSolver provides the pressure and flow rate at all defined nodes and edges. For a given network, the solution can be obtained manually or semi-automatically using computer algebra systems. However, this approach requires manual problem definition as a system of linear equations. For a single, well-defined configuration, this approach is goal-directed. However, for the iterative investigation and optimization of a large number of network configurations in dialogue with the future user, this approach is inefficient.
Kirchhoff's first law is adapted from the conservation of charge to conserving volume in microfluidic networks. The total sum of all in-and outflowing volumes equals zero at each defined node of a microfluidic network. A node has no volume on its own. Thus, each change of geometric component leading to a change in volumes, such as junction and cross-section changes, is represented by the edge connecting the nodes. Kirchhoff's second energy conservation law is related to the balance of the hydrodynamic potential differences (pressure) around any closed network.
For a simplified demonstration, the exemplary network in Figure 4a is considered. The network consists of the nodes N1, N2, and N3 connected by the edges E3, E4, and E5. E1, E2, and E6 are linked to the microfluidic control (MFC) units via pressure-controlled Input-Output (IO) ports IO1, IO2, and IOref, in Figure 4b. Depending on the pressure settings, they either serve es inlets or outlets. The network consisting of six nodes and six edges is translated into the system of linear equations, resulting in a binary incidence matrix with six columns and six rows. The columns represent the nodes and ports. The rows refer to the edges. The edge directed away from the node is assigned +1 and the arriving edge is assigned −1. The weights of the respective edges correspond to the inverse of the hydrodynamic resistance R.

MfnSolver Concept and Implementation
The mfnSolver provides the pressure and flow rate at all defined nodes and edges. For a given network, the solution can be obtained manually or semi-automatically using computer algebra systems. However, this approach requires manual problem definition as a system of linear equations. For a single, well-defined configuration, this approach is goal-directed. However, for the iterative investigation and optimization of a large number of network configurations in dialogue with the future user, this approach is inefficient.
Kirchhoff's first law is adapted from the conservation of charge to conserving volume in microfluidic networks. The total sum of all in-and outflowing volumes equals zero at each defined node of a microfluidic network. A node has no volume on its own. Thus, each change of geometric component leading to a change in volumes, such as junction and cross-section changes, is represented by the edge connecting the nodes. Kirchhoff's second energy conservation law is related to the balance of the hydrodynamic potential differences (pressure) around any closed network.
For a simplified demonstration, the exemplary network in Figure 4a is considered. The network consists of the nodes N1, N2, and N3 connected by the edges E3, E4, and E5. E1, E2, and E6 are linked to the microfluidic control (MFC) units via pressure-controlled Input-Output (IO) ports IO1, IO2, and IO ref , in Figure 4b. Depending on the pressure settings, they either serve es inlets or outlets. The network consisting of six nodes and six edges is translated into the system of linear equations, resulting in a binary incidence matrix with six columns and six rows. The columns represent the nodes and ports. The rows refer to the edges. The edge directed away from the node is assigned +1 and the arriving edge is assigned −1. The weights of the respective edges correspond to the inverse of the hydrodynamic resistance R.  The pressure results from Kirchhoff's modified second low describing the volume flow through edges with where Q E denominates the flow rate vector, MI the incidence matrix, C the conductivity vector and p N is the vector providing the pressure at the nodes N. Thus, MI T ⋅ Q = Q Source follows for the conservation of volume, and MI T ⋅ (MI ⋅ C) ⋅ p N = Q Source as the matrix expression for solving pN. The edges and nodes of the introduced network and their chosen RVL model are implemented into the solver analog to Figure 4b. In the interest of the numerical stability of the solver, a rescaling of the base units is performed to avoid decimal numbers with high exponents. The unit length is defined as 1 mm within the system, the unit mass relates to 1 mg, and the unit time is 1 s. As a result, the parameters and respective units in Table 1 follow the network diagram.
Specific models to represent the channel geometry of the on-chip network and the respective tubing of the MFC unit are implemented. The on-chip hydrodynamic resistance R is calculated based on the Hagen-Poiseuille law where L is the length of the respective edge and dH is the hydraulic diameter with with the area of the cross-section A = WM ⋅ ED + 0.5 ⋅ π ⋅ (ED) 2 and the perimeter P = 2.0 ⋅ (WM + ED) + π ⋅ ED. Consequently, volume is defined as V = A ⋅ L. For the standard tubing, which connects the chip network with MFC, the calculation of the area is adjusted to represent a circular cross-section with A = π 4 ⋅ ID 2 and ID as the inner diameter of the tubing.  The pressure results from Kirchhoff's modified second low describing the volume flow through edges with where Q E denominates the flow rate vector, MI the incidence matrix, C the conductivity vector and p N is the vector providing the pressure at the nodes N. Thus, MI T · Q = Q Source follows for the conservation of volume, and MI T · (MI · C) · p N = Q Source as the matrix expression for solving p N . The edges and nodes of the introduced network and their chosen RVL model are implemented into the solver analog to Figure 4b. In the interest of the numerical stability of the solver, a rescaling of the base units is performed to avoid decimal numbers with high exponents. The unit length is defined as 1 mm within the system, the unit mass relates to 1 mg, and the unit time is 1 s. As a result, the parameters and respective units in Table 1 follow the network diagram. Specific models to represent the channel geometry of the on-chip network and the respective tubing of the MFC unit are implemented. The on-chip hydrodynamic resistance R is calculated based on the Hagen-Poiseuille law where L is the length of the respective edge and d H is the hydraulic diameter with with the area of the cross-section A = WM · ED + 0.5 · π · (ED) 2 and the perimeter P = 2.0 · (WM + ED) + π · ED. Consequently, volume is defined as V = A · L. For the standard tubing, which connects the chip network with MFC, the calculation of the area is adjusted to represent a circular cross-section with A = π 4 · ID 2 and ID as the inner diameter of the tubing.
The mfnSolver is implemented using Python with the numerical library Numpy for the matrix operations and the toolkit Graphviz for visualization and result presentation. The mfnSolver is consistently implemented object-oriented with the class hierarchy shown in Figure 5a. The class mfnObject manages essential attributes and provides them to the derived classes mfnNode (nodes of the network) and mfnEdge (edges of the network). The class mfnDomain manages subgraphs of a fluidic network consisting of nodes and edges, which can be assumed as sub-domains of the complete network, e.g., microfluidic chips with different configurations or higher-level control networks, are derived. The class mfnNetwork forms the core of the solver. It aggregates an arbitrary number of mfnDomain objects into a microfluidic network, creates and solves the Kirchhoff-matrix for the implemented system, and generates the network diagrams in .pdf and output lists in .csv format. Particular microfluidic circuits and subdomains can be implemented as reusable mfnDomain objects common in microelectronics for circuit design and simulation. The mfnSolver is implemented using Python with the numerical library Numpy for the matrix operations and the toolkit Graphviz for visualization and result presentation. The mfnSolver is consistently implemented object-oriented with the class hierarchy shown in Figure 5a. The class mfnObject manages essential attributes and provides them to the derived classes mfnNode (nodes of the network) and mfnEdge (edges of the network). The class mfnDomain manages subgraphs of a fluidic network consisting of nodes and edges, which can be assumed as sub-domains of the complete network, e.g., microfluidic chips with different configurations or higher-level control networks, are derived. The class mfnNetwork forms the core of the solver. It aggregates an arbitrary number of mfnDomain objects into a microfluidic network, creates and solves the Kirchhoff-matrix for the implemented system, and generates the network diagrams in .pdf and output lists in .csv format. Particular microfluidic circuits and subdomains can be implemented as reusable mfnDomain objects common in microelectronics for circuit design and simulation. Thus, the solution shown in Figure 5b considers the microfluidic chip geometry as well as the microfluidic control units. The pressure and flow rate at all defined nodes and edges of the network are provided and can be used as initial operating conditions in pressure-driven experiments. The solution of the solver is obtained in less than 200 ms, even for more complex networks.
In the presented generic mfnSolver, based on the Kirchhoff nodal analysis, each edge on the MFN can be individually configured by choosing the appropriate RVL model. Therefore, the solver provides a set of pre-implemented RVL models representing the pressure-throughput characteristics of standard channel shapes such as a circular tube, capillary slit, rectangular channel, ellipsoidal channel, and half-ellipsoidal channel implementation, and hydraulic diameter-based approximation for channels with an arbitrary cross-sectional shape.

PDG2 Chip and Network Model
As an example and reference case, the mfnSolver predicts the settings for pressuredriven flow control for single and multi-phase flow utilizing the flow-focusing unit of the PDG2 chip, in Figure 6. The PDG2 on-chip MFN consists of five ports, of which four are utilized as inlets, and one is used as an outlet to waste. The three inlets P2, P3, and P5 are combined in sequence to create a multi-phase laminar flow when operated at low velocities below the turbulence threshold with Reynolds Numbers Re < π 2 . In fluid dynamics, laminar flow is characterized by a parabular shape of the velocity gradient over the cross-section orthogonal to flow direction. This laminar dispersed phase is focused through the two-arm flow-focusing unit of inlet P1 by a continuous phase. Our PDG2 chip can be utilized for diffusion and reaction experiments at operation based on laminar co-flow. Thus, the solution shown in Figure 5b considers the microfluidic chip geometry as well as the microfluidic control units. The pressure and flow rate at all defined nodes and edges of the network are provided and can be used as initial operating conditions in pressure-driven experiments. The solution of the solver is obtained in less than 200 ms, even for more complex networks.
In the presented generic mfnSolver, based on the Kirchhoff nodal analysis, each edge on the MFN can be individually configured by choosing the appropriate RVL model. Therefore, the solver provides a set of pre-implemented RVL models representing the pressure-throughput characteristics of standard channel shapes such as a circular tube, capillary slit, rectangular channel, ellipsoidal channel, and half-ellipsoidal channel implementation, and hydraulic diameter-based approximation for channels with an arbitrary cross-sectional shape.

PDG2 Chip and Network Model
As an example and reference case, the mfnSolver predicts the settings for pressuredriven flow control for single and multi-phase flow utilizing the flow-focusing unit of the PDG2 chip, in Figure 6. The PDG2 on-chip MFN consists of five ports, of which four are utilized as inlets, and one is used as an outlet to waste. The three inlets P2, P3, and P5 are combined in sequence to create a multi-phase laminar flow when operated at low velocities below the turbulence threshold with Reynolds Numbers Re < π 2 . In fluid dynamics, laminar flow is characterized by a parabular shape of the velocity gradient over the crosssection orthogonal to flow direction. This laminar dispersed phase is focused through the two-arm flow-focusing unit of inlet P1 by a continuous phase. Our PDG2 chip can be utilized for diffusion and reaction experiments at operation based on laminar co-flow. Depending on the viscosity of the liquids and the surface functionalization of the channels, the FFU in Figure 7a can create laminar co-flows of up to four components or droplets with variable internal composition from up to three sample fluids. Depending on the viscosity of the liquids and the surface functionalization of the channels, the FFU in Figure 7a can create laminar co-flows of up to four components or droplets with variable internal composition from up to three sample fluids. The generic mfnSolver is applied to this chip geometry for MFN analysis and calculation of initial experimental pressure parameters for set flow rates to guarantee steady operation. Based on the fabrication parameter mask width WM and etch depth ED, shown in Figure 7b, the hydrodynamic resistance is calculated based on the respective cross-sectional area. The Kirchhoff-graph of the network for laminar co-flow embedded in the MFC, consisting of operation pumps, capillaries, and flow sensing units tubing to connect the control unit with the on-chip PDG2 network, is computed with the mfnSolver. The pumps 1, 2, 3, and 4 of the MFC are connected to the ports P1, P2, P3, and P5, respectively, from left to right. The outlet port P4, connected to the waste, serves as a reference with pressure p = 0 mbar. The solver uses the standard flow rate Q = 0.16 µL/s. While the flow rate of one port is gradually increased according to Table 2, the remaining ports are kept at a standard flow rate. In total, four different sets are modeled to predict the respective pressure outcomes analog to Table 2. The respective mfnSolver output prediction for the feedlines of an operation according to acquisition parameter ID 1 is exemplary shown in the condensed scheme in Figure 8 (top). The detailed output is included in the digital supplement. In between the sets, the port with the varying flow rate is changed. The solver provides the predicted pressure, which is used as the input for pressure-driven operation during the experimental validation. The Kirchhoff-graph allows tracking of the pressure drop between the nodes over the entire network. In addition, all connecting edges provide information on the selected RVL model, including their parameters, the pressure drop, the set flow rate, and the implemented dynamic viscosity η of the liquid the channel carries. The dynamic viscosity is the governing parameter that has to be adjusted to adapt the model for multi-component droplet generation. For the laminar co-flow scenario, the dynamic viscosity is assumed to The generic mfnSolver is applied to this chip geometry for MFN analysis and calculation of initial experimental pressure parameters for set flow rates to guarantee steady operation. Based on the fabrication parameter mask width WM and etch depth ED, shown in Figure 7b, the hydrodynamic resistance is calculated based on the respective crosssectional area. The Kirchhoff-graph of the network for laminar co-flow embedded in the MFC, consisting of operation pumps, capillaries, and flow sensing units tubing to connect the control unit with the on-chip PDG2 network, is computed with the mfnSolver. The pumps 1, 2, 3, and 4 of the MFC are connected to the ports P1, P2, P3, and P5, respectively, from left to right. The outlet port P4, connected to the waste, serves as a reference with pressure p = 0 mbar. The solver uses the standard flow rate Q = 0.16 µL/s. While the flow rate of one port is gradually increased according to Table 2, the remaining ports are kept at a standard flow rate. In total, four different sets are modeled to predict the respective pressure outcomes analog to Table 2. The respective mfnSolver output prediction for the feedlines of an operation according to acquisition parameter ID 1 is exemplary shown in the condensed scheme in Figure 8 (top) and in more detail in Figure S1 in the digital supplement. In between the sets, the port with the varying flow rate is changed. The solver provides the predicted pressure, which is used as the input for pressure-driven operation during the experimental validation. Table 2. Exemplary simulation for laminar co-flow containing five flow and pressure profiles. The model flow rate as the initial solver input is opposed to the predicted pressure for each setting. The full parameter set is provided in Table S1. mixture of the oil phase and the droplets. The viscosity of this two-phase flow has been derived from the experiments by analyzing the pressure drop between the FFU as the droplet generating structure and the outlet for the utilized fixed flow rate of the formed droplet suspension. Here, a combined measured flowrate of 320 µl/s from the inlets P2, P3 and P5 was achieved by a pressure offset of 80 mbar. A shortened example of the mfn-Solver prediction for the feedline pressure settings and the FFU, to create droplet flow, can be seen in Figure 8 (bottom). Figure 8. mfnSolver output with predicted pressure settings for the feedlines to create laminar coflow in the PDG2 chip device (top) (complete mfnSolver output in S1: mfnSolver output Kirchhoffgraph co-flow); mfnSolver output with predicted pressure settings for the feedlines to create multi- Figure 8. mfnSolver output with predicted pressure settings for the feedlines to create laminar co-flow in the PDG2 chip device (top) (complete mfnSolver output in Figure S1: mfnSolver output Kirchhoff-graph co-flow); mfnSolver output with predicted pressure settings for the feedlines to create multi-component droplet flow in the PDG2 chip device (bottom) (also see Figure S2: mfnSolver output Kirchhoff-graph droplet flow).
The Kirchhoff-graph allows tracking of the pressure drop between the nodes over the entire network. In addition, all connecting edges provide information on the selected RVL model, including their parameters, the pressure drop, the set flow rate, and the implemented dynamic viscosity η of the liquid the channel carries. The dynamic viscosity is the governing parameter that has to be adjusted to adapt the model for multi-component droplet generation. For the laminar co-flow scenario, the dynamic viscosity is assumed to match the water in all channels. This assumption also holds during the droplet experiments for the feedlines connected to P2, P3, and P5. P1 is operated with perfluorinated oil with η = ν · ρ = 0.77 cm 2 /s · 1.614 g/cm 3 = 1.24278 g/(cm·s). After the FFU, the fluid is a mixture of the oil phase and the droplets. The viscosity of this two-phase flow has been derived from the experiments by analyzing the pressure drop between the FFU as the droplet generating structure and the outlet for the utilized fixed flow rate of the formed droplet suspension. Here, a combined measured flowrate of 320 µl/s from the inlets P2, P3 and P5 was achieved by a pressure offset of 80 mbar. A shortened example of the mfnSolver prediction for the feedline pressure settings and the FFU, to create droplet flow, can be seen in Figure 8 (bottom).
To validate the presented solver, flow sensors are integrated into the feedlines to monitor the flow rate during pressure-driven operation and compare the solver results with flow measurements. In addition, the on-chip in-line flow at three different positions is measured using an RGB chemometric image-based analysis to validate the solver and compared with the respective analytical predictions of the solver. Both are discussed in the next section.

Experimental Validation
Within this section, the experimental validation of the introduced mfnSolver is discussed. First, flow sensors were integrated into the interconnecting capillaries. The correlation between predicted target flow rates and the experimental observation for the delivery of the fluids from the pressurized feed vessels to the Inlet-Outlet (IO)-ports of the microfluidic chip device were analyzed.
Secondly, as the experimental design does not simply integrate flow sensors into existing microfluidic devices, we switched to an optical detection method that allows measurements at arbitrary positions of the optically transparent microfluidic device. The flow-focusing chip PDG2, as introduced in the previous section, is utilized for these experiments. We generate laminar co-flows of dyed sample fluids inside the channels and perform a colorimetric analysis along profile lines orthogonal to the channel axis for the measurements. The composition of the dyed sample fluids has been varied according to the parameter profiles in Figure 9 (top) and their respective output composition (middle) and their outcome field of view (FOV) for experimental validation with the exemplary pressure profile of parameter set identifier (ID) 1 (bottom). A selection of the recorded image data for the calculated pressure settings of laminar co-flow and droplet flow is presented in Video S1, S2 and S3.

Validation of Flow Rate Prediction at Feedlines
The flow sensors are integrated into the feedline between the PEEK capillaries to monitor the flow rate during the experiment. Their performance characteristics without chip operation can be found in the Supplementary Information. A standard flow rate of 0.16 µL/s has been chosen for the model. Within each set, the flow rate through one channel is changed in five steps from 0.08 µL/s to 1.28 µL/s, while the remaining ones are set constant to the chosen standard flow rate as input for the model. The respective predicted pressure settings have been implemented as operation pressures for the PDG2 chip. The pressure is gradually increased from 50 mbar to 800 mbar, as predicted in Table 2. Each pressure setting is kept constant for 10 s. A recording of 90 s per pressure setting is performed to demonstrate flow stability of the laminar co-flow. The pressure of a port is increased from 50 mbar to 800 mbar, according to the mfnSolver calculations for laminar co-flow, which illustrates the starting and end condition of the laminar co-flow series for each set. The flow pattern was steady during all experiments, required for image-based analysis and a correlation between the phase and flow fractions.
focusing chip PDG2, as introduced in the previous section, is utilized for these experiments. We generate laminar co-flows of dyed sample fluids inside the channels and perform a colorimetric analysis along profile lines orthogonal to the channel axis for the measurements. The composition of the dyed sample fluids has been varied according to the parameter profiles in Figure 9 (top) and their respective output composition (middle) and their outcome field of view (FOV) for experimental validation with the exemplary pressure profile of parameter set identifier (ID) 1 (bottom). The performance of pressure control and flow sensors is illustrated in the subsection 'Microfluidic control' in the chapter 'Materials and Methods'. The comparison of the pressure-throughput characteristics predicted for laminar co-flow and the experimental parameters of the pressure control and the feedline flow sensors are shown in Figure 10a. We can conclude that the presented model based on the Kirchhoff nodal analysis can predict the overall pressure-throughput of the laminar co-flow scenario. Thus, it can be used to provide pressure-driven flow-control in the feedlines to the microfluidic chip.
For droplet generation, a laminar co-flow pattern of the three sample fluids is created as the dispersed phase. The flow-focusing unit is now functionalized as a droplet generator. The total flow rate of the dispersed phase and the continuous phase is kept constant, while the color pattern of the laminar co-flow is varied, as illustrated in Figure 11. Here, we consider evenly shaped droplets fully embedded and carried by the continuous phase. The comparison of the pressure-throughput characteristics predicted for droplet-flow and the experimental parameters are shown in Figure 10b.
'Microfluidic control' in the chapter 'Materials and Methods'. The comparison of the pressure-throughput characteristics predicted for laminar co-flow and the experimental parameters of the pressure control and the feedline flow sensors are shown in Figure 10a. We can conclude that the presented model based on the Kirchhoff nodal analysis can predict the overall pressure-throughput of the laminar co-flow scenario. Thus, it can be used to provide pressure-driven flow-control in the feedlines to the microfluidic chip. For droplet generation, a laminar co-flow pattern of the three sample fluids is created as the dispersed phase. The flow-focusing unit is now functionalized as a droplet generator. The total flow rate of the dispersed phase and the continuous phase is kept constant, while the color pattern of the laminar co-flow is varied, as illustrated in Figure 11. Here, we consider evenly shaped droplets fully embedded and carried by the continuous phase. The comparison of the pressure-throughput characteristics predicted for droplet-flow and the experimental parameters are shown in Figure 10b.

Validation of Chip Internal Flow
The majority of microfluidic devices are operated at low Reynolds Numbers with Re π as confirmed by W. Wuest [26] for channels with integrated apertures. Therefore, the formation of secondary flow patterns due to inertial forces like vortices, jets, local backflow zones, Dean-Flow patterns, or time-periodic fluctuations (e.g., Karman Vortex Street) can be neglected. The laminar flow regime allows the creation of a laminar co-flow of multiple fluid components, where each fluid component passes through the channel as a separated lamella. Diffusion is responsible for outbalancing the concentrations of the ingredients. This process is comparably slow and leads to a widening of the individual lamella and finally to a complete outbalancing of the concentration of the ingredients over time. An example for the generation of a three-component laminar co-flow in the flowfocusing chip PDG2 with the fluid components Bromophenol blue ('BPB'), Orange G ('OG'), and the water buffer ('W') is demonstrated in Figure 9. Here, the on-chip feedlines come from the top and bottom and feed into the observation channel, flowing from left to right.

Data acquisition and Required Datasets
Image data are recorded by transmission mode microscopy utilizing a white light source with a condenser, a microscopy objective, and an RGB camera. Details are given in the 'Materials and Methods' section. Each pixel of the RGB camera acts as a three-channel photon counting device where the delivered intensity values linearly scale with the number of inclining photons. The camera is applied to measure concentrations. The raw intensities must be converted to a measured quantity, which linearly scales with the concentration c of a light-absorbing dye. According to the Beer-Lambert Law, the absorbance values A as established in UV-VIS spectroscopy and photometry conform with these require-

Validation of Chip Internal Flow
The majority of microfluidic devices are operated at low Reynolds Numbers with Re < π 2 as confirmed by W. Wuest [26] for channels with integrated apertures. Therefore, the formation of secondary flow patterns due to inertial forces like vortices, jets, local backflow zones, Dean-Flow patterns, or time-periodic fluctuations (e.g., Karman Vortex Street) can be neglected. The laminar flow regime allows the creation of a laminar co-flow of multiple fluid components, where each fluid component passes through the channel as a separated lamella. Diffusion is responsible for outbalancing the concentrations of the ingredients. This process is comparably slow and leads to a widening of the individual lamella and finally to a complete outbalancing of the concentration of the ingredients over time. An example for the generation of a three-component laminar co-flow in the flow-focusing chip PDG2 with the fluid components Bromophenol blue ('BPB'), Orange G ('OG'), and the water buffer ('W') is demonstrated in Figure 9. Here, the on-chip feedlines come from the top and bottom and feed into the observation channel, flowing from left to right.

Data acquisition and Required Datasets
Image data are recorded by transmission mode microscopy utilizing a white light source with a condenser, a microscopy objective, and an RGB camera. Details are given in the 'Materials and Methods' section. Each pixel of the RGB camera acts as a three-channel photon counting device where the delivered intensity values linearly scale with the number of inclining photons. The camera is applied to measure concentrations. The raw intensities must be converted to a measured quantity, which linearly scales with the concentration c of a light-absorbing dye. According to the Beer-Lambert Law, the absorbance values A as established in UV-VIS spectroscopy and photometry conform with these requirements. Additionally, the transmitted optical path length d must be taken into account For the conversion of the received image data to absorbance data, the following equation can be derived from Equation (5) where IM SAMPLE denominates the recorded image of the chip during an experimental run, IM DARK is the recorded image with the light source switched off as the dark reference, and IM BLANK is the recorded image of the chip for all channels filled with the buffer solution.
Care must be taken in the implementation to avoid exceptions due to division by zero and logarithm of zero values. Our implementation manages these cases by substituting values below the limit of detection (LOD) in the numerator and denominator of Equation (6) by the LOD value, calculated as two times pixel noise standard deviation as measured in the acquired IM BRIGHT , which is captured from the pure light source without the chip.

Single-Pixel Chemometric Analysis
For the quantitative chemometric analysis, the following quantities and definitions are utilized.
The area-related amount of substance N describes the amount of a pure component detected at a given area element. The physical meaning results from the Beer-Lambert Law in Equation (5) with where c denominates a volumetric concentration and d the effective optical path length (NA > 0). The parameter can be interpreted as the amount of substance of a component, given in mol/m 2 , detected by a single camera pixel. Per definition, the amount of a pure component detected for d = d max , at the nominal channel height H, equals 1. The phase fraction α i denominates the fraction of a pure component i contained in a mixture of n components and is given by: Assuming a volumetric mixture of n components, the phase fraction α i can also be calculated from the volume fractions V i of the mixture. The same applies to the volume flow Q i of pure components in a laminar co-flow scenario.
The measured RGB absorbance A PIXEL:RGB of a single pixel reflects the absorbance of a mixture of the fluid components where each component i is characterized by unique spectral absorbance properties A i:RGB . For analysis, we want to know the phase fraction α i for each of the contributing components i. Therefore, we have the spectral properties of the pure components A i.:RGB available as a reference. The observed single-pixel RGB absorbance spectrum A OBS:RGB of a mixture of n pure components where each component contributes with an amount N i is described by the following system of equations: A OBS:RGB and A i:RGB are retrieved from the measurements. Thus, Equation (10) can be solved to obtain Ni by the numerical solver for chemometric and spectral unmixing. The solver is based on a least-square method, which has been recently proposed and discussed by N. Coca [27].
RGB values are retrieved from single-pixel readout and converted to absorbance values A according to Equation (6). The area-related amount of substance N i for the components is obtained from these RGB absorbances. For each pixel measurement, the phase fractions α i for Bromophenol blue ('BPB'), Orange G ('OG'), and the water buffer ('W') fluid components are calculated from the respective amount N i divided by the total amount N T , exemplary demonstrated for 'BPB': The sum of phase fractions of the pure components accumulates to 1.0 for the respective component as the mixture analysis yields α i for all the contributing components 'OG', 'BPB', and 'W':

Accumulated Phase Fractions in a Channel
In the previous section, we have demonstrated the calculation of the phase fractions α i for a single pixel measurement. For laminar co-flow experiments, the phase fractions α are measured for each acquired sample image at different channel positions. Therefore, the data are analyzed along measurement boxes, spanning the channel width, as illustrated in Figure 12. Box positions and dimensions can be specified for an arbitrary number of boxes.

Accumulated Phase Fractions in a Channel
In the previous section, we have demonstrated the calculation of the phase fractions αi for a single pixel measurement. For laminar co-flow experiments, the phase fractions α are measured for each acquired sample image at different channel positions. Therefore, the data are analyzed along measurement boxes, spanning the channel width, as illustrated in Figure 12. Box positions and dimensions can be specified for an arbitrary number of boxes. The phase fraction is calculated by accumulation over all pixels of the box inside the channel. Due to the particular shape of our channels shown in Figure 12, the local channel height h(y) changes depending on the position. Thus, the optical path length d becomes a positional dependent variable. While the channel height h(y) can be calculated precisely from the known channel shape, the optical path length d is additionally influenced by the NA of the illumination and the refraction of the illumination beam at the curved channel walls. The prediction of these effects is challenging. Therefore, the effective optical path length d(y) is derived from the absorbance of the reference IMFILLED, according to Equations (5) and (6), for a user-definable color-channel. To compensate for the variation of the channel height, d(y) is normalized towards maximum channel height dmax. Thus, the phase fraction αi accumulated over the channel results in The phase fraction is calculated by accumulation over all pixels of the box inside the channel. Due to the particular shape of our channels shown in Figure 12, the local channel height h(y) changes depending on the position. Thus, the optical path length d becomes a positional dependent variable. While the channel height h(y) can be calculated precisely from the known channel shape, the optical path length d is additionally influenced by the NA of the illumination and the refraction of the illumination beam at the curved channel walls. The prediction of these effects is challenging. Therefore, the effective optical path length d(y) is derived from the absorbance of the reference IM FILLED , according to Equations (5) and (6), for a user-definable color-channel. To compensate for the variation of the channel height, d(y) is normalized towards maximum channel height d max . Thus, the phase fraction α i accumulated over the channel results in where y 1 and y 2 denominate the wall positions of the channel. The accumulated phase fractions α i provide the phase fractions of all fluid components, passing through the channel at the boxed regions. For comparison of the in-line image data and the mfnSolver model prediction, the phase fraction of each line over the channel is calculated analog to Equation (13), while the mfnSolver output takes the ratio of the sample flowrate over the total in-line flow rate of the respective channel. Five measurement boxes are distributed over the main channel length, as shown in Figure 12. The correlation of measured phase fractions with solver-predicted flow fractions is given in Figure 13a. A slope of 1.0 denominates the perfect fit. The figure shows good agreement between the measured phase fractions in the channel sections and the initially predicted flow settings from the mfnSolver, with a slope of 1.03 (R 2 = 0.97). Details on the composition of the three-component co-flow for the individual subsets of the measurement and its relation to the mfnSolver predicted compositions are given in Figure 13b as a phase diagram. The predicted flow fractions show good agreement with the measured ones, and also their point of intersection almost coincides. -sf

Validation of Predicted Flow Rates Based on Droplet Composition
In this section, the total flow rate of the dispersed phase and the continuous phase is kept constant, as shown in Figure 14 (top). Only the phase fractions of the dispersed phase are varied, as illustrated in Figure 14 (middle). The FOV with the respective measurement boxes and the region for droplet detection and measurement are shown in Figure 14 (bottom). Details on the composition of the three-component co-flow for the individual subsets of the measurement and its relation to the mfnSolver predicted compositions are given in Figure 13b as a phase diagram. The predicted flow fractions show good agreement with the measured ones, and also their point of intersection almost coincides.

Validation of Predicted Flow Rates Based on Droplet Composition
In this section, the total flow rate of the dispersed phase and the continuous phase is kept constant, as shown in Figure 14 (top). Only the phase fractions of the dispersed phase are varied, as illustrated in Figure 14 (middle). The FOV with the respective measurement boxes and the region for droplet detection and measurement are shown in Figure 14 (bottom). Boxes 1, 2, and 3 mentioned in Figure 14 (bottom) show a laminar co-flow, analyzed analogously to the previously described section. The correlation of the co-flow phase fractions compared to the predicted flow fractions show good agreement in Figure 15a, which fits the description in the previous section. However, it can be seen that the non-symmetric color-combination from the inlets leads to more outliers and higher deviations between the measured and calculated phase fractions from the ideal slope.
The droplet-flow phase fractions are calculated by accumulating all pixels inside a droplet using the scheme described for a box measurement. The droplet regions are detected by binarization of pre-processed input images with a given threshold, followed by filling the holes and eroding the binary mask to the inner volume of the droplets. The color components inside the drops are qualitatively analyzed. While the overall linearity can still be seen from the experimental data in Figure 15b, systematic deviations are observed, further discussed in the following section. Boxes 1, 2, and 3 mentioned in Figure 14 (bottom) show a laminar co-flow, analyzed analogously to the previously described section. The correlation of the co-flow phase fractions compared to the predicted flow fractions show good agreement in Figure 15a, which fits the description in the previous section. However, it can be seen that the nonsymmetric color-combination from the inlets leads to more outliers and higher deviations between the measured and calculated phase fractions from the ideal slope.
The droplet-flow phase fractions are calculated by accumulating all pixels inside a droplet using the scheme described for a box measurement. The droplet regions are detected by binarization of pre-processed input images with a given threshold, followed by filling the holes and eroding the binary mask to the inner volume of the droplets. The color components inside the drops are qualitatively analyzed. While the overall linearity can still be seen from the experimental data in Figure 15b

Precision of the Chemometric Analysis in Optofluidics
We confirm a straight correlation of the measured phase fractions with the solver predicted flow fractions in the measured data. At the same time, we recognize 'color'specific systematic deviations for the different colored fluids for the laminar co-flow patterns. For the measurements of the droplet contents, significantly higher deviations are observed. In general, image data acquisition is performed in transmission mode. Obviously, this raises the question of whether distortions of the optical beams during sample fluid and channel transmission may be assigned to these deviations. The following effects can be considered: 1. Optical refraction at discrete refractive index transitions at the curved channel sidewalls; 2. Optical refraction at dynamic refractive index changes at the boundary of the two fluid samples, with different refractive indices. Here, one has to note that adding a dye to a buffer changes its refractive index. 3. Chromatic aberration (chromatic distortion) due to refraction, which applies to items 1 and 2.
These effects will influence the correctness of the measured transmitted optical path length of the fluid related to local channel heights and droplet heights. Furthermore, they may cause a shift in the spectral properties detected by the camera pixels at the affected positions, which leads to flawed absorbance spectra as the input for the chemometric analysis. The impact of these effects on the accuracy of the measured values can be assessed using intermediate data from the data processing chain available for each experiment.
A dynamic refractive index change at the boundary between 'W' (n1) and 'BPB' (n2) generates a curved dynamic refractive index transition. The resulting refraction is related to chromatic aberration and a shift in the detected spectral signatures for this region, which can be observed in Figure 16. These deviations lead to misleading results in the chemometric analysis. In this case, we observe false-positive amounts for the magnitudes of 'OG' and phase-contrast ('PC') with opposing signs, as indicated with the arrows in Figure 16.

Precision of the Chemometric Analysis in Optofluidics
We confirm a straight correlation of the measured phase fractions with the solver predicted flow fractions in the measured data. At the same time, we recognize 'color'-specific systematic deviations for the different colored fluids for the laminar co-flow patterns. For the measurements of the droplet contents, significantly higher deviations are observed. In general, image data acquisition is performed in transmission mode. Obviously, this raises the question of whether distortions of the optical beams during sample fluid and channel transmission may be assigned to these deviations. The following effects can be considered:

1.
Optical refraction at discrete refractive index transitions at the curved channel sidewalls; 2.
Optical refraction at dynamic refractive index changes at the boundary of the two fluid samples, with different refractive indices. Here, one has to note that adding a dye to a buffer changes its refractive index.
These effects will influence the correctness of the measured transmitted optical path length of the fluid related to local channel heights and droplet heights. Furthermore, they may cause a shift in the spectral properties detected by the camera pixels at the affected positions, which leads to flawed absorbance spectra as the input for the chemometric analysis. The impact of these effects on the accuracy of the measured values can be assessed using intermediate data from the data processing chain available for each experiment.
A dynamic refractive index change at the boundary between 'W' (n 1 ) and 'BPB' (n 2 ) generates a curved dynamic refractive index transition. The resulting refraction is related to chromatic aberration and a shift in the detected spectral signatures for this region, which can be observed in Figure 16. These deviations lead to misleading results in the chemometric analysis. In this case, we observe false-positive amounts for the magnitudes of 'OG' and phase-contrast ('PC') with opposing signs, as indicated with the arrows in Figure 16.  The curved shape of the channel sidewalls results in refraction and chromatic aberration, leading to false amounts of the components at near-wall positions, as shown in Figure 17. The valid area of analysis is the center of the channel. In addition, we observe an alternating sign of the deviation, which might be related to an alignment error of the illumination beam to the optical axis of the imaging system. The curved shape of the channel sidewalls results in refraction and chromatic aberration, leading to false amounts of the components at near-wall positions, as shown in Figure 17. The valid area of analysis is the center of the channel. In addition, we observe an alternating sign of the deviation, which might be related to an alignment error of the illumination beam to the optical axis of the imaging system.
The same combination of refraction and chromatic aberration is observed in proximity to the curved interface of the investigated droplets, as illustrated in Figure 18. The refractive index change between the continuous oil phase and the droplet components enhances these optical effects. The curved shape of the channel sidewalls results in refraction and chromatic aberration, leading to false amounts of the components at near-wall positions, as shown in Figure 17. The valid area of analysis is the center of the channel. In addition, we observe an alternating sign of the deviation, which might be related to an alignment error of the illumination beam to the optical axis of the imaging system. The curved shape of the channel sidewalls results in refraction and chromatic aberration, leading to false amounts of the components at near-wall positions, as shown in Figure 17. The valid area of analysis is the center of the channel. In addition, we observe an alternating sign of the deviation, which might be related to an alignment error of the illumination beam to the optical axis of the imaging system. The curved shape of the channel sidewalls results in refraction and chromatic aberration, leading to false amounts of the components at near-wall positions, as shown in Figure 17. The valid area of analysis is the center of the channel. In addition, we observe an alternating sign of the deviation, which might be related to an alignment error of the illumination beam to the optical axis of the imaging system.
The same combination of refraction and chromatic aberration is observed in proximity to the curved interface of the investigated droplets, as illustrated in Figure 18. The refractive index change between the continuous oil phase and the droplet components enhances these optical effects. The same combination of refraction and chromatic aberration is observed in proximity to the curved interface of the investigated droplets, as illustrated in Figure 18. The refractive index change between the continuous oil phase and the droplet components enhances these optical effects.
All these effects are contributing to the correctness of the measured values. They could be outbalanced by using refractive index-matched special fluids, as proposed by D. Malsch [23]. However, these special fluids are incompatible with any biological assay or chemical reaction. Therefore, the decision was made towards fluid compositions, typical for a broad range of biological and molecular biology assays, such as Phosphate Buffer at pH 8.0 and perfluorinated oils for droplet experiments. The results of this section demonstrate how optofluidic effects are related to the quality and precision of the on-chip chemometric measurements and are intended to provide valuable hints for implementing these measurement approaches in real lab-on-a-chip applications. All these effects are contributing to the correctness of the measured values. They could be outbalanced by using refractive index-matched special fluids, as proposed by D. Malsch [23]. However, these special fluids are incompatible with any biological assay or chemical reaction. Therefore, the decision was made towards fluid compositions, typical for a broad range of biological and molecular biology assays, such as Phosphate Buffer at pH 8.0 and perfluorinated oils for droplet experiments. The results of this section demonstrate how optofluidic effects are related to the quality and precision of the on-chip chemometric measurements and are intended to provide valuable hints for implementing these measurement approaches in real lab-on-a-chip applications.

Discussion and Outlook
While much work has already been done towards microfluidic design automation [13,[16][17][18] and control, to our knowledge, a general and easily implementable solution is still unavailable to the community. The developed mfnSolver can be understood as one core component of the MDA architecture, which provides a reliable system simulation for the fast and iterative optimization of microfluidic networks. Thus, it is applicable for the design of lab-on-a-chip systems and the in-line prediction of the correct settings for their operation.
For implementing the microfluidic network graph of a particular chip design, it might be difficult to exactly access and measure the geometrical parameters of commercially available microfluidic devices. However, the mfnSolver is comparably robust to length variations, as long as the majority of the pressure drop can be assigned to the interconnecting tubes. For example, a variation of four mm at different positions of the onchip network only leads to an average difference across pressure predictions of 2.337 mbar with a sample standard deviation of 0.364 mbar. On the one hand, the solver calculates the hydrodynamic resistance as displayed in Equation (3). Thus, the length is only factored in once in the nominator of the equation.
On the other hand, if the network is designed so that the main pressure drop occurs along the feedlines, the network behavior is dominated by the throughput characteristics of the interconnecting system. In this case, the Kirchhoff-graphs show that only less than 5% of the pressure drop across the entire system drops over the microfluidic chip itself. The most significant pressure drop is contributed by the capillaries connecting the microfluidic control with the chip. The respective manufacturer usually provides measures of the inner diameter of the capillaries. It is also the crucial parameter in calculating the hydraulic diameter, which appears with a power of four in Equation (3). Exemplarily, a variation of 20 µm across the four capillaries can lead to a deviation of over 15% of the pressure output provided by the solver. Therefore, it is necessary to calibrate the inner diameter of the capillaries before usage.

Discussion and Outlook
While much work has already been done towards microfluidic design automation [13,[16][17][18] and control, to our knowledge, a general and easily implementable solution is still unavailable to the community. The developed mfnSolver can be understood as one core component of the MDA architecture, which provides a reliable system simulation for the fast and iterative optimization of microfluidic networks. Thus, it is applicable for the design of lab-on-a-chip systems and the in-line prediction of the correct settings for their operation.
For implementing the microfluidic network graph of a particular chip design, it might be difficult to exactly access and measure the geometrical parameters of commercially available microfluidic devices. However, the mfnSolver is comparably robust to length variations, as long as the majority of the pressure drop can be assigned to the interconnecting tubes. For example, a variation of four mm at different positions of the on-chip network only leads to an average difference across pressure predictions of 2.337 mbar with a sample standard deviation of 0.364 mbar. On the one hand, the solver calculates the hydrodynamic resistance as displayed in Equation (3). Thus, the length is only factored in once in the nominator of the equation.
On the other hand, if the network is designed so that the main pressure drop occurs along the feedlines, the network behavior is dominated by the throughput characteristics of the interconnecting system. In this case, the Kirchhoff-graphs show that only less than 5% of the pressure drop across the entire system drops over the microfluidic chip itself. The most significant pressure drop is contributed by the capillaries connecting the microfluidic control with the chip. The respective manufacturer usually provides measures of the inner diameter of the capillaries. It is also the crucial parameter in calculating the hydraulic diameter, which appears with a power of four in Equation (3). Exemplarily, a variation of 20 µm across the four capillaries can lead to a deviation of over 15% of the pressure output provided by the solver. Therefore, it is necessary to calibrate the inner diameter of the capillaries before usage.
While we only present the mfnSolver for applications utilizing pressure-driven flow control, the solver is not limited to this application range. The mfnSolver can be further applied to iterative MFN design optimizations to tailor on-chip channel characteristics to meet desired control parameters. Furthermore, several networks can be implemented into the solver and stitched together to analyze and simulate parallel operation. Due to the high computational speed, the mfnSolver can be easily integrated into closed-loop control applications.

Final Conclusions
In this work, we confirmed the outstanding potential of the microfluidic network solver (mfnSolver) as a reliable tool for fast and precise prediction of pressure settings for steady-state operation conditions and getting insights into on-chip pressure drop and flow characteristics. In contrast to common CFD analyzers, the mfnSolver can analyze a complete microfluidic network within a few milliseconds. The mfnSolver can be implemented and tailored based on the chip fabrication parameters, etch depth, mask width, and on-chip lengths between the defined nodes and available microfluidic control units. It precisely predicts the operation parameters for laminar co-flow and droplet operation within the established microfluidic framework.
Furthermore, the solver has been experimentally validated to create and analyze laminar co-flow and droplets with systematically variable compositions. In combination with chemometric analysis, absorbance imaging has been found as a precise and contactless optical method for quantitative monitoring of transport processes in microfluidic devices. The impact of optofluidic effects on the precision of the chemometric analysis has been discussed.