Circuit-Based Design of Microfluidic Drop Networks

Microfluidic-drop networks consist of several stable drops—interconnected through microfluidic channels—in which organ models can be cultured long-term. Drop networks feature a versatile configuration and an air–liquid interface (ALI). This ALI provides ample oxygenation, rapid liquid turnover, passive degassing, and liquid-phase stability through capillary pressure. Mathematical modeling, e.g., by using computational fluid dynamics (CFD), is a powerful tool to design drop-based microfluidic devices and to optimize their operation. Although CFD is the most rigorous technique to model flow, it falls short in terms of computational efficiency. Alternatively, the hydraulic–electric analogy is an efficient “first-pass” method to explore the design and operation parameter space of microfluidic-drop networks. However, there are no direct electric analogs to a drop, due to the nonlinear nature of the capillary pressure of the ALI. Here, we present a circuit-based model of hanging- and standing-drop compartments. We show a phase diagram describing the nonlinearity of the capillary pressure of a hanging drop. This diagram explains how to experimentally ensure drop stability. We present a methodology to find flow rates and pressures within drop networks. Finally, we review several applications, where the method, outlined in this paper, was instrumental in optimizing design and operation.


Introduction
Recent advances in microfabrication techniques, such as the 3D printing of new types of resins, have significantly increased the design potential for microfluidic devices, allowing for almost arbitrarily complex interweaving channel networks [1,2]. Optimizing the design can be performed by modeling the behavior of fluids within microfluidic channel networks with intricate numerical methods, such as the finite element method. However, owing to computational limits, 3D models of large microfluidic networks cannot be handled efficiently, particularly when exploring the effects of several design and operation parameters. Therefore, a separate "first-pass" method is required for predicting flow within such networks before establishing more intricate 3D mathematical models.
A powerful method, used for such initial device optimization, is the hydraulic-circuit analogy [3]. This analogy has found various applications in the field of microfluidics [4][5][6]. This approach is based on working with volumetric flow rates Q and fluid pressures p as if they were electric currents and voltages, respectively. While this approach is commonly used for microfluidic channels with completely closed channels, it is rarely used for open fluidic devices with dynamic drop volumes. One of the reasons is that a drop that changes in size cannot be modeled with simple electrical components, such as constant resistors or sources, but require variable components to represent the nonlinear dependency between drop geometry and capillary and hydrostatic pressures.
In this review article, we show the methodology, developed in our group to apply the hydraulic-circuit analogy for the purpose of designing and operating microfluidic hangingand standing-drop networks.
First, we describe the basic useful microfluidic components and all relevant variables to model them. Then, we discuss the practical limitations of drop networks, with a focus on hanging-drop networks. We then lay out a strategy to draw a microfluidic-network scheme that helps in defining the set of equations describing the microfluidic circuit. Solving the set of equations then provides a description of drop volume dynamics and helps to find critical drop behaviors that need to be considered when designing and operating drop networks. Finally, we show how this method can be and has been applied in the context of different and previously published open drop-based microfluidic networks.

Basic Microfluidic-Drop-Network Components
We first present the main elements used for designing and drawing microfluidic networks: hydraulic resistances, fluidic sources, atmospheric pressure, and drop compartments.

Hydraulic Resistance
Hydraulic resistance is directly analogous to resistive elements within electric circuit design. Solving the Navier-Stokes equation for flow through a microfluidic channel results in the velocity vector field v(x, y, z) and the pressure scalar field p(x, y, z). The flow rate Q through the microfluidic channel is found by integrating the velocity vector field over the channel cross-sectional area Ω (Figure 1a). This cross section is typically set at the inlet or outlet of the channel, where the flow should be orthogonal to Ω (n Ω ), which then yields a proportional relationship between Q and the pressure drop between the inlet and outlet of the channel ∆p = p in − p out . scheme that helps in defining the set of equations describing ing the set of equations then provides a description of drop to find critical drop behaviors that need to be considered wh drop networks. Finally, we show how this method can be context of different and previously published open drop-bas

Basic Microfluidic-Drop-Network Components
We first present the main elements used for designing an works: hydraulic resistances, fluidic sources, atmospheric p ments.

Hydraulic Resistance
Hydraulic resistance is directly analogous to resistive el design. Solving the Navier-Stokes equation for flow throug sults in the velocity vector field ( , , ) and the pressure sc rate through the microfluidic channel is found by integra over the channel cross-sectional area Ω (Figure 1a). This cr the inlet or outlet of the channel, where the flow should be then yields a proportional relationship between and the p let and outlet of the channel Δ = − . This proportional relationship is called the hydraulic res is entirely dependent on the channel design and on the vis flowing through the channel. The most general definition of [3] with the perimeter and the area of the cross sect tween the inlet and outlet of the channel. This proportional relationship is called the hydraulic resistance of the channel R and is entirely dependent on the channel design and on the viscosity η of the liquid phase flowing through the channel. The most general definition of this relationship is given in [3] with the perimeter P and the area A of the cross section Ω, and L the length between the inlet and outlet of the channel.
We include useful definitions of hydraulic resistances for circular R O , high-aspect-ratio rectangular R , and rectangular R channel cross sections that remain constant through the channel length ( Figure 1b).
These are the basic geometries that are typically used in the literature to design microfluidic devices, where the width w is defined to be larger than the height h of a channel.

Fluidic Source
A fluidic source either acts as a fluid volume source or drain, as is the case in the flow-through operation of microfluidic chips, or acts as an internal actuator of flow or pressure, as is the case in recirculating operation of microfluidic chips. Each fluidic source uses a force to generate a pressure that drives flow (Table 1). Examples are: (i) surface tension causing a capillary pressure p c that generates a pressure differential between two material phases, e.g., at an air-liquid interface (ALI), or (ii) gravity causing a hydrostatic pressure p g that generates a pressure differential between two elements at different heights h. A summary of common fluidic sources, the forces driving fluid motion in these sources, their definitions, and the physical limitations on the achievable flow rates is provided in Table 1 and schematically illustrated in Figure 2a-c. At their core, all fluidic sources are pressure-driven. "True" pressure-driven sources apply a constant pressure that is a function of their geometry. They are limited by the design properties of the source itself. For example, capillary-driven flow is limited by the minimum achievable fluid surface curvature (Figure 2c), and gravity-driven flow is limited by the maximum height difference generating pressure differentials (Figure 2a). With pressure-driven flow, pressure differentials through the liquid phase of various elements drive fluid motion.
Conversely, with flow-driven flow, liquid motion generates pressure differentials through the liquid phase. Therefore, flow-driven sources apply a variable pressure in order to ensure a constant flow rate. They are limited by the mechanical properties of the chambers generating the flow. For example, peristaltic pumps are limited by the compliance of the peristaltic tubing, leading to unpredictable flow rates for highly resistive microfluidic devices. Syringe pumps are limited by the strength of the syringe itself, which can shatter, or the water tightness of the plunger driving fluid motion in the syringe, which can leak.

Atmospheric Pressure
For the purpose of microfluidic-network design, atmospheric pressure takes the role of the "ground" in the electric circuit design. Atmospheric pressure acts as a reference pressure and ensures that we can accurately compute the absolute pressure in each microfluidic element. Regardless of the driving force of fluid motion, pressure differentials play a crucial role in predicting microfluidic network behavior. However, the absolute value of the pressure within different microfluidic elements can cause various unwanted effects. In closed microfluidic devices, high pressure can cause failure of microfluidic channels, akin to Blaise Pascal's barrel bursting due to hydrostatic pressure [7]. In open microfluidic devices, high pressure can cause a seemingly static system to burst and leak in a highly dynamic manner ( Figure S1).

Drop Compartments
A circular opening in an open microfluidic network is referred to as a "drop compartment". Drop compartments can either be "hanging" (Figure 3a) or "standing" (Figure  3b), depending on the chip concept and the experimental configuration. A drop compartment of radius can be accurately modeled in a microfluidic circuit as a capillary pressure source , linked to the network through a hydrostatic pressure source and a resistance (Figure 3c,d). We approximate the resistance as that of a circular channel of a length equal to the drop height ℎ.
The singular difference between a standing and hanging drop, for the purpose of a hydraulic circuit design, is the sign of the hydrostatic pressure. Flow toward the ALI is therefore driven by Δ = + − for a hanging drop, and Δ = − − for a standing drop (Figure 3c,d).

Atmospheric Pressure
For the purpose of microfluidic-network design, atmospheric pressure takes the role of the "ground" in the electric circuit design. Atmospheric pressure acts as a reference pressure and ensures that we can accurately compute the absolute pressure in each microfluidic element. Regardless of the driving force of fluid motion, pressure differentials play a crucial role in predicting microfluidic network behavior. However, the absolute value of the pressure within different microfluidic elements can cause various unwanted effects. In closed microfluidic devices, high pressure can cause failure of microfluidic channels, akin to Blaise Pascal's barrel bursting due to hydrostatic pressure [7]. In open microfluidic devices, high pressure can cause a seemingly static system to burst and leak in a highly dynamic manner ( Figure S1).

Drop Compartments
A circular opening in an open microfluidic network is referred to as a "drop compartment". Drop compartments can either be "hanging" (Figure 3a) or "standing" (Figure 3b), depending on the chip concept and the experimental configuration. A drop compartment of radius a can be accurately modeled in a microfluidic circuit as a capillary pressure source p c , linked to the network p n through a hydrostatic pressure source p g and a resistance R c (Figure 3c,d). We approximate the resistance R c as that of a circular channel of a length equal to the drop height h.
The singular difference between a standing and hanging drop, for the purpose of a hydraulic circuit design, is the sign of the hydrostatic pressure. Flow toward the ALI Q c is therefore driven by ∆p = p n + p g − p c for a hanging drop, and ∆p = p n − p g − p c for a standing drop (Figure 3c  . Circuit eq of fluid movement in (c) a hanging and (d) a standing drop. At equilibrium, the equivale ± holds, and no flow is induced. An increase in due to hydrostatic pressure (e.g., the chip or elevating an inlet reservoir) or hydrodynamic pressure (e.g., fluid motion w fluidic network) generates a net flow = Δ ⁄ toward the drop ALI. An increase in dro due to the ALI reduces . A reduction in induces an increase in , which eventually new equilibrium with , stopping flow toward the ALI at a new drop volume.
In practice, since is low relative to the channel hydraulic resistance, the a drop interface acts as a compliant element. Drops react almost instantly to changes in fluid pressure due to hydrostatic and hydrodynamic forces. Howev behavior is particularly difficult to predict when submitted to rapid pressure chan when tilting a chip, elevating an inlet reservoir, or upon step changes in flow rate the higher channel resistance in the microfluidic-drop network. The drop behavi sudden pressure changes needs to be recapitulated with the modeling methodol lined in this section. An exercise to intuit drop dynamics is to follow a step increa drop compartment pressure in a standing drop:  Table 2); Table 2);  Drop geometry changes ( and ℎ as highlighted in Figure 3a,b) induce changes: o A reduction in the ALI radius increases the capillary pressure ; o An increase in the relative height increases the hydrostatic pressure  An increase in and changes the pressure differential Δ = − −  A reduction in pressure differential reduces the flow rate toward the ALI.
This process repeats iteratively until the pressure differential is null and ca scribed by the following differential equation: Looking at the expression for ( , ) and ℎ( , ) in the full-page Table 2, E (6) is difficult to solve analytically, even for a single drop. The highly nonlinear between drop volume and pressure makes it impossible to model the drop us At equilibrium, the equivalence p n = p c ± p g holds, and no flow is induced. An increase in p n due to hydrostatic pressure (e.g., by tilting the chip or elevating an inlet reservoir) or hydrodynamic pressure (e.g., fluid motion within the fluidic network) generates a net flow Q c = ∆p/R c toward the drop ALI. An increase in drop volume due to the ALI reduces r. A reduction in r induces an increase in p c , which eventually reaches a new equilibrium with p n , stopping flow toward the ALI at a new drop volume.
In practice, since R c is low relative to the channel hydraulic resistance, the air-liquid drop interface acts as a compliant element. Drops react almost instantly to gradual changes in fluid pressure due to hydrostatic and hydrodynamic forces. However, drop behavior is particularly difficult to predict when submitted to rapid pressure changes, e.g., when tilting a chip, elevating an inlet reservoir, or upon step changes in flow rates, due to the higher channel resistance in the microfluidic-drop network. The drop behavior upon sudden pressure changes needs to be recapitulated with the modeling methodology outlined in this section. An exercise to intuit drop dynamics is to follow a step increase in the drop compartment pressure p n in a standing drop: • A standing-drop compartment pressure increase drives flow toward the ALI Q c ; Q c > 0 increases the drop volume V → ∆V > 0; ∆V > 0 reduces the ALI radius (equation r(V, a) in Table 2); ∆V > 0 increases the drop height (equation h(V, a) in Table 2); • Drop geometry changes (r and h as highlighted in Figure 3a,b) induce pressure changes: A reduction in the ALI radius increases the capillary pressure p c ; An increase in the relative height increases the hydrostatic pressure p g ; • An increase in p c and p g changes the pressure differential ∆p = p n − p g − p c ; • A reduction in pressure differential reduces the flow rate toward the ALI.
This process repeats iteratively until the pressure differential is null and can be described by the following differential equation: Looking at the expression for r(V, a) and h(V, a) in the full-page Table 2, Equation (6) is difficult to solve analytically, even for a single drop. The highly nonlinear relation between drop volume and pressure makes it impossible to model the drop using only constant electrical components. Instead, hydrostatic and capillary pressures need to be modeled using variable pressure sources, which iteratively change according to differential Equation (6). Therefore, for a more complex problem, e.g., a series of drops in a microfluidic network, we used numerical methods. Table 2. Useful relationships between the various geometric elements of a spherical cap are included herein. Equations for subhemispherical caps are written in blue and equations for suprahemispherical caps are written in green.

Practical Limitation of Drop Systems
The main limitation of drop systems is when the drop becomes too large, causing leakage from the microfluidic network. For the sake of simplification, we chose a hangingdrop setup to explain this phenomenon, referred to herein as "drop crash". In this case, drop crash, explained with Figure 4, is when sustained hydraulic and hydrostatic pressure on the ALI exceed the capillary pressure ( Figure 3a), causing the hanging drop to fall off the chip. However, a standing-drop setup also suffers from this limitation; albeit, works with in keeping drop integrity (Figure 3b). Excessive flow rates, however, can easily generate a that exceeds + . Generally, if a pressure increase causes a flow that increases the drop volume past a

Practical Limitation of Drop Systems
The main limitation of drop systems is when the drop becomes too large, causing leakage from the microfluidic network. For the sake of simplification, we chose a hangingdrop setup to explain this phenomenon, referred to herein as "drop crash". In this case, drop crash, explained with Figure 4, is when sustained hydraulic p n and hydrostatic p g pressure on the ALI exceed the capillary pressure p c (Figure 3a), causing the hanging drop to fall off the chip. However, a standing-drop setup also suffers from this limitation; albeit, p g works with p c in keeping drop integrity (Figure 3b). Excessive flow rates, however, can easily generate a p n that exceeds p c + p g .   Generally, if a pressure increase causes a flow that increases the drop volume past a hemisphere, then the ALI radius starts to increase. Therefore, instead of the capillary pressure increasing to balance internal hydraulic pressure, it is reduced and diverges from equilibrium. This nonlinear divergence drives flow into the drop at an increasing rate. The divergence can then cause a drop to grow past the limit, given by the critical capillary length l c = γ/ρg [8], which is characteristic of a fluid of a given density ρ subjected to a gravitational acceleration g. We call this phenomenon a "crash" of the system, as gravitational forces overtake capillary forces on the drop ALI. Figure 4 explains the dynamics of drop crashing. Figure 4a shows, on the y-axis, the pressure, normalized by the maximum capillary pressure of a hanging drop, versus, on the x-axis, the volume within the hanging drop, normalized by the volume of a hemispherical drop. The normalization and equations giving the plotted curves are presented in the Supplementary Material section "Hanging drop crash". The black line, showing the normalized drop capillary pressure p c , is drawn with Equation (S2). The dashed black line, showing normalized internal pressure p n + p g , is drawn with Equation (S4).
An exercise to interpret Figure 4a is to follow the p c starting at a V = 0.02, i.e., an empty drop (Figure 4d). In this case, the V increases stably until p c = p n + p g . Significantly increasing ∆h from Figure 4d leads to the state plotted in Figure 4e. In this case, V increases slowly until V = 1 (minimal difference between p c and p n + p g ); then, V rapidly increases until the drop crashes (divergent difference between p c and p n + p g ). Reducing ∆h from Figure 4d leads to the state plotted in Figure 4f. In this case, V decreases slowly until p c = p n + p g .
An unstable volume increase (Figure 4e) can be recuperated by reducing the drop height to a technical minimum of ∆h = 0. However, if this is performed too late, as would be the case with a large extrahemispherical drop, it leads to the state plotted in Figure 4g. In this case, a volume decrease is impossible, and the drop is irrecuperable. A pseudostable drop condition exists, but any perturbation will cause a drop crash.
In practice, Figure 5 can be consulted, where the state of a drop is entirely defined by experimental pressure and volume. Figure 5a shows a linear representation of the normalized pressure. Acting experimentally on a drop network changes the pressure p n resulting in varied internal pressures (Figure 5a dashed lines). We plot p n on the y-axis of Figure 5b. The observable dimension of a drop is its volume V on the x-axis of Figure 5. The response of a drop is seen as a color axis with the pressure difference ∆p = p n + p g − p c . The drop dynamics is either a volume increase (red) or decrease (blue) until it reaches the stable black curve. This results in the phase diagram for a drop compartment of a given diameter (2a = 3 mm). The critical drop volume above which gravitational forces cause a drop crash is also shown with a dashed red line.
The drop dynamics explained in this section are crucial for developing an intuition for manipulating drop systems. We include other phase diagrams for several drop compartment apertures in Figure S2. . The drop dynamics is either a volume increase (red) or decrease (blue) until it reaches the stable black curve. This results in the phase diagram for a drop compartment of a given diameter (2a = 3 mm). The critical drop volume above which gravitational forces cause a drop crash is also shown with a dashed red line.
The drop dynamics explained in this section are crucial for developing an intuition for manipulating drop systems. We include other phase diagrams for several drop compartment apertures in Figure S2.

Methodology for Modeling Drop Networks
We will describe the method that we elaborated to model drop platforms with the basic elements of microfluidic-drop networks that we defined and detail their limitations. We first define the known and unknown variables of the model that we then use to define the microfluidic network scheme. This microfluidic network scheme allows us to define the set of equations describing the microfluidic circuit. Solving the set of equations with a matrix approach allows us to gain valuable insights into microfluidic-drop-network design and operation. The described methodology is generally applicable to any microfluidic-drop network, but we will demonstrate it on a specific example, referred to as the "scalable microphysiological system".

Variable Definition
Known and unknown variables must be properly defined and enumerated to define the set of equations that describe fluid flow within the microfluidic network and to numerically solve the unknown variables. All variables are summarized in Table 3. Table 3. Variable definition for each microfluidic element. Subscripts in the nomenclature are descriptive of the element, where c, r, v, and g denote capillary, reservoir, valve, and gravity. Subscripts are also enumerative, where i denotes the i'th node, ij denotes an element between the i'th and j'th nodes, and in and out denote the inlet and outlet. Known variables are either design-or operation-related. They include hydraulic resistances, prescribed flow rates, capillary pressures, and hydrostatic pressures.

Microfluidic
Hydraulic resistance is defined by the channel design. Microfluidic channels typically have a constant resistance given by Equations (2), (3), or (4), depending on their cross section and channel length. However, reservoirs [9], drops [10][11][12][13], and other elements that can fill or empty over time will lead to varying channel lengths and, therefore, are modeled as continuously variable hydraulic resistances. The filling and emptying dynamics of such reservoirs is discussed in Supplementary Materials section "Filling or emptying reservoirs". Additionally, valving mechanisms [14] will induce a step decrease or increase in the resistance of a channel, either allowing or restricting flow along a specific path.
Flow rates can be prescribed during the operation of a microfluidic device. Mechanical fluidic sources (Table 1) are the only sources able to supply a steady flow rate. The flow rate can be constant when a continuous flow through the microfluidic network is defined, e.g., for continuous sampling [10,15], substance dosing [16], or medium replenishment [11]. The prescribed flow rate can also vary, e.g., for applying varying shear stresses [17] or periodically increasing liquid turnover [11].
Capillary pressures are an operational variable, given by the geometry of a curved ALI, which are determined by liquid cohesion and interfacial adhesion forces. Capillary pressure varies throughout an experiment to comply with liquid pressure within a drop compartment. However, if a steady state of flow is reached, e.g., in flow-through applications [11], then the capillary pressure can be constant.
Hydrostatic pressure is an operational variable that is given by the relative vertical position of a microfluidic element. Hydrostatic pressure can be constant throughout an experiment, where, if the chip is kept static, the level of the chip returns a hydrostatic pressure map through the device, which induces a change in volume from drop to drop. Hydrostatic pressure can also vary, as is the case with gravity-driven flow, e.g., upon tilting the chip, upon unequally filling the chip's inlet and outlet reservoirs, or upon connecting an inlet to an elevated reservoir [9].

Unknown Variables
Unknown variables need to be defined and computed to recapitulate microfluidic network dynamics. They include flow rates through each microfluidic element and pressures at each node between microfluidic elements.
Flow rates need to be defined between each node of the microfluidic circuit. The direction of the flow and the corresponding flow rate must be set, where a negative flow rate indicates a reversed flow. The best strategy is to define flow rate variables according to each hydraulic resistance, mirroring the resistance nomenclature.
Pressures need to be defined at each node and are used to compute flow rates through the device. Pressure decreases linearly in the direction of the flow across each microfluidic channel with a constant cross section.

Microfluidic-Network Scheme
Knowing the common microfluidic elements, defined in Table 3, allows for drawing the electric equivalent circuit of a microfluidic-channel network. We show an example of a circuit design for a published standing-drop chip (Figure 6a) [9]. The standing-drop chip consists of 10 standing-drop compartments (Figure 6b), arranged in series and flanked by reservoirs on both sides (Figure 6a). The standing drops have a defined geometry (Figure 6c) allowing the definition of flow resistances through the chip. A single-compartment circuit is shown (Figure 6d) and repeated N = 10 times. The full circuit is schematically shown (Figure 6e), where each open surface is connected to atmospheric pressure. Figure 6e is used to define the set of circuit equations. The strategy to draw the electric equivalent circuit is elaborated in the Supplementary Material section named "Circuit design strategy".

Circuit Equations
As soon as a circuit for a microfluidic network is fully defined, the corresponding set of equations can be written down in a matrix form. Equations take two forms: (i) equations describing the flow through resistive elements, linking pressure drops to flow rates, and (ii) equations describing mass conservation at nodes.

Resistive Pressure Drop
The relationship between hydraulic resistance, flow rate, and pressure must be satisfied in every part of the microfluidic network that features flow. This results in a number of equations that is equal to the number of resistive elements in the device. Equations of the resistive pressure drop are of the form p i − p j = Q ij R ij , where the pressure drop from node i to node j drives a flow rate Q ij through the resistive element R ij .

Fluidic Nodal Rule
The nodal rule, also known as Kirchhoff's current law, ensures conservation of charge within a circuit. Equivalently, the fluidic nodal rule ensures the conservation of volume within a microfluidic network. This rule looks at each node and equates all influxes in that node to its outfluxes. This results in a number of equations that is equal to the number of nodes in the device. Equations of the fluidic nodal rule are of the form Q hi = Q ci + Q ij , where the flow rate Q hi flowing from node h to node i equals the sum of the flow rate Q ci flowing from node i to drop i and the flow rate Q ij going from node i to node j (Figure 6d, blue arrows).

Circuit Matrix
The resistive pressure drop and fluidic nodal rule are used to define a set of equations.
The number of defined equations should be equal to the number of unknown variables. If there are more equations, then the system is overdetermined, which happens when equations are linearly dependent. Linearly dependent equations should be identified and combined. The set of equations are then written in a matrix notation Ax = y. The known vector y contains the known pressures or flow rates. The unknown vector x contains all unknown flow rates and all unknown pressures. The characteristic matrix A typically contains hydraulic resistances R i in most of the diagonal elements and −1 or 1 in the relevant off-diagonal elements. Here, we show the resulting matrix equation Ax = y for the example of a gravity-driven tilting chip, shown in Figure 6, with N = 3.

Circuit Equations
As soon as a circuit for a microfluidic network is fully defined, the corresponding set of equations can be written down in a matrix form. Equations take two forms: (i) equations describing the flow through resistive elements, linking pressure drops to flow rates, and (ii) equations describing mass conservation at nodes. The set of equations describing a pressure-driven microfluidic-drop network typically yields a symmetric sparse matrix. Rigorously defining the circuit, variables, and equations, which describe the flow within a microfluidic network, gives significant freedom in redesigning the channels, reconfiguring the network, and changing the device operation conditions. For example, the chip from Figure 6, described by Equation (7), can be operated with a continuous flow rate instead of being gravity-driven [16]. In this case, the known variable p g−in from Equation (7) becomes an unknown inlet pressure p in , and the unknown variable Q in from Equation (7) becomes a known inlet flow rate Q in . The resulting sparse matrix is no longer symmetric: Reconfiguring the matrix and solving it with the new inlet resistances R in that included the tubing for the fluidic source evidenced that driving the chip with a continuous flow made it much more sensitive to leaks. This finding led to a redesign of the chip and its setup to be able to drive it with a continuous flow for the purpose of applying pharmacokinetic drug concentration exposures [16].

Insights from Solving Circuit Dynamics
The circuit matrix, previously defined in Equation (7), results in a set of differential equations relating the rate of change in the volume of the i'th drop Q ci = dV i /dt to the microfluidic-network design and operation parameters. Numerically solving this set of differential equations resulted in various useful values:

•
The flow rate through each element Q ij was computed and could be used to estimate shear stresses within a channel [12,13]. In turn, by modifying the channel design of a microfluidic chip, shear stresses could be tuned to ensure that physiological shear stress levels are achieved in adherent cell cultures; • The flow rates at nodes and their ratios could be used to estimate molecule mixing within a microfluidic network [18]; • Reservoir volumes V ri were computed and could be used to optimize the tilting scheme to ensure that flow rates remained relatively constant during the operation of a microfluidic chip [12,13]; • Inlet pressure p in for flow-rate-driven chips could be evaluated and minimized to avoid excessive pressure on fluidic sources; • Drop volumes V i were computed as a function of time and could be used to optimize the operation of self-controlled hanging-drop setups [10]; • Drop volumes V i were computed as a function of time and could be minimized to obviate catastrophic drop failures [10][11][12].

Applications
The method outlined in this review has been applied to several open microfluidic devices designed and published within our laboratory. Previously, the modeling methodology was only described superficially, as it was out of scope of the respective papers. Here, we will give detailed information and elaborate on the insights gained by the modeling of three of these devices. The first was a standing-drop microfluidic network, developed and redesigned by Lohasz et al. [9,16,18], herein referred to as "Scalable microphysiological system". The second was a hanging-drop microfluidic network, developed by Boos et al. [12,13], herein referred to as "Placenta-on-a-chip". The third was a hanging-drop microfluidic network, first developed by Misun et al. [15] and redesigned by Wu Jin et al. [10], herein referred to as "Hanging-drop-based islet perifusion system".

Scalable Microphysiological System
The scalable microphysiological system was the subject of the model presented in Figure 6, which was used to predict flow rates through the device that then were validated experimentally (Figure 7). This system was operated by tilting it back and forth, inducing continuous gravity-driven bidirectional flow. We cultured microtissues within the compartments that represented different organ models. The constant perfusion enabled continuous interaction and interorgan communication. With the validated model, various operational conditions could be simulated to predict device operation and optimize tilting schemes for gravity-driven flow at room (20 °C) and incubator (37 °C) temperatures. Temperature changes induced change in the medium's viscosity and density, which impacted flow. First, the maximum achievable flow rate at 0 min, when the prescribed volume was first added to the device, was plotted as a function of the tilting angle ( Figure 7d). Then, by computing the time before the inlet reservoir was fully emptied, the maximum tilting interval was plotted as a function of the tilting angle (Figure 7d). Figure 7g, like Figure 7d, was also generated to predict The time dependence of the flow rate could easily be determined with the model, as it was a direct output Q out from solving Equation (7). Experimentally, this flow rate was more tedious to validate. Setting the chip at a given tilting angle α of 10 • , 20 • , and 30 • , a prescribed volume of 150 µL was deposited in the inlet. The liquid was subsequently completely sampled at the outlet at a given time point. The process of liquid prescription and sampling was repeated for every time point indicated by the measured values in Figure 7a-c.
With the validated model, various operational conditions could be simulated to predict device operation and optimize tilting schemes for gravity-driven flow at room (20 • C) and incubator (37 • C) temperatures. Temperature changes induced change in the medium's viscosity and density, which impacted flow. First, the maximum achievable flow rate Q in at 0 min, when the prescribed volume was first added to the device, was plotted as a function of the tilting angle α (Figure 7d). Then, by computing the time before the inlet reservoir was fully emptied, the maximum tilting interval was plotted as a function of the tilting angle α (Figure 7d). Figure 7g, like Figure 7d, was also generated to predict the effect of various initial inlet reservoir liquid volumes V rin (100, 150, and 200 µL) on the tilting intervals. As the volumes increased, the different liquid levels increased the initial hydrostatic pressure. These tilting intervals were used to optimize chip operation to ensure continuous flow during microtissue culture and to prevent the complete draining of the elevated (source) reservoir. The optimized tilting cycle can be seen in Figure 7h and was used to generate the published results [9].

Placenta-on-a-Chip
The placenta-on-a-chip was the second iteration of the microfluidic multitissue platform for advanced embryotoxicity testing, invented by Boos et al. [12], that recapitulated the maternal-placental-embryonic axis [13]. The chip, shown in Figure 8a, featured two liquid phases, a "maternal" and "embryonic" side, separated by a semipermeable membrane on which placental cells grew. Placental cells were grown on the membrane in the closed maternal compartment, which was accessible by two hanging drops on either side. Embryoid bodies were cultured in hanging drops underneath the placental barrier on the embryonic side. Gravity-driven flow was induced through the network by continuously tilting the device by ±5 • .
The maternal side featured a large culture chamber (Figure 8a) in which placental cells were seeded to form a confluent layer as a model of the placental barrier. We modeled the shear stresses on the placental barrier along the black cut line shown in Figure 8b. With this shear stress model, we showed that the placental barrier on the maternal-embryonic circular boundary, delineated by the dashed red lines in Figure 8b, was exposed to a uniform shear stress. We optimized the microfluidic design of the maternal side in order to minimize the shear stress by using the modeling techniques elaborated in the previous section to predict flow as a function of a tilting angle of ±5 • (Figure 8c). The optimization consisted of including a resistive serpentine channel (Figure 8b) and adding additional reservoir drops to increase the volume throughput at each tilting step. We ensured that the intervals of the maximum flow rates within the maternal side were sustained as long as possible within each tilting cycle (Figure 8c) to promote cell polarization on the barrier. We achieved a maximal flow rate during chip operation of 4.6 µL min −1 , which we used to compute the shear stress in Figure 8b. timization consisted of including a resistive serpentine channel (Figure 8b) and adding additional reservoir drops to increase the volume throughput at each tilting step. We ensured that the intervals of the maximum flow rates within the maternal side were sustained as long as possible within each tilting cycle (Figure 8c) to promote cell polarization on the barrier. We achieved a maximal flow rate during chip operation of 4.6 μL min −1 , which we used to compute the shear stress in Figure 8b.

Hanging-Drop-Based Islet Perifusion System
A platform to perform glucose-stimulated insulin secretion (GSIS) assays with single pancreatic islets was presented by Misun et al. [15]. The platform enabled studying the dynamics of insulin release by islet microtissues at high temporal resolution. Perifusion of single-islet microtissues produced biologically more relevant results than the typical pooling and analysis of several microtissues. The platform enabled the detection of complex pulsatile insulin secretion at single-islet resolution, while pooling of microtissues only would yield average values. The platform was adapted by Wu Jin et al. [10] to parallelize insulin sampling, and allowed for the simultaneous investigation of four pancreatic islets in separate hanging drops (Figure 9). of single-islet microtissues produced biologically more relevant results than the typical pooling and analysis of several microtissues. The platform enabled the detection of complex pulsatile insulin secretion at single-islet resolution, while pooling of microtissues only would yield average values. The platform was adapted by Wu Jin et al. [10] to parallelize insulin sampling, and allowed for the simultaneous investigation of four pancreatic islets in separate hanging drops (Figure 9).  Like in the original device [15], sampling was carried out by applying a continuous flow through hanging drops hosting the individual islet microtissues and sampling at the outlets (green region in Figure 10A). An additional drop adjustment fluidic branch was added with a control drop to set the drop height (blue region in Figure 10A). This branch allowed automatically regulating and stabilizing drop heights throughout the network without the need of a microscopy-based feedback (in principle, setting an outlet with a needle-type valve allows keeping the network drop height constant). The self-regulation of the drop height, depicted in Figure 10B, was demonstrated experimentally [10]. The control outlet aspirated liquid, as long as the control drop ALI was lower than the needle outlet. Once the drop receded above the needle outlet, the control outlet aspirated air. This procedure repeated throughout the experiment, controlling the drop height.
allowed automatically regulating and stabilizing drop heights throughout the network without the need of a microscopy-based feedback (in principle, setting an outlet with a needle-type valve allows keeping the network drop height constant). The self-regulation of the drop height, depicted in Figure 10B, was demonstrated experimentally [10]. The control outlet aspirated liquid, as long as the control drop ALI was lower than the needle outlet. Once the drop receded above the needle outlet, the control outlet aspirated air. This procedure repeated throughout the experiment, controlling the drop height. There was a continuous flow of liquid from the inlet toward the islet drops due to active sampling from the sampling outlets at the rate of out sampling . An irregular flow was observed from the inlet toward the control drop due to the alternating withdrawal of liquid and air from the needle-type valve at a rate of out control . (B) Self-regulation of drop heights between interconnected hanging drops. Liquid was constantly added into the system at a rate of and withdrawn through the sampling outlet at a rate of out S and through a needle-type valve at a rate of out C . (i) Two hanging drops in equilibrium with identical Laplace radius and pressure. The inflow was evenly distributed to the two drops. As long as the tip of the needle-type valve was immersed in the drop, liquid removal through the valve occurred. Consequently, (ii) the size and the Laplace pressure of the two drops became different with drop radius and pressure being inversely correlated. Due to the pressure difference, there was an increased flow toward the control drop (left) with the lower pressure. As soon as the tip of the needle-type valve was exposed to air, only air was withdrawn. The control drop remained stable at the height defined by the valve needle length. (iii) In the next step, the two drops reached equal Laplace radius and pressure through equilibration through the liquid phase. The system was now stable, and the drop height in control and tissue drop was maintained constant. Reproduced from Wu Jin et al. [10] (CC BY 4.0).
However, a phenomenon known as capillary-wetting hysteresis [19,20] caused certain irregularities in the drop height control. This phenomenon, described in detail in Figure S3 in the Supplementary Materials section "contact-angle hysteresis", induced an inhomogeneous wetting of the needle depending on whether the drop size was increasing or decreasing. The wetting hysteresis of an advancing or receding air-liquid-solid tripoint induced more complex wetting dynamics. The wetting dynamics were modeled by including a hysteresis in the drop height (dashed red line in Figure 11) and successfully recapitulated the dynamics observed experimentally. This result showed that, although The inflow (Q in ) was split into two flows, one toward the islet drops (right) and the other toward the control drop (left). There was a continuous flow of liquid from the inlet toward the islet drops due to active sampling from the sampling outlets at the rate of Q out sampling . An irregular flow was observed from the inlet toward the control drop due to the alternating withdrawal of liquid and air from the needle-type valve at a rate of Q out control . (B) Self-regulation of drop heights between interconnected hanging drops. Liquid was constantly added into the system at a rate of Q in and withdrawn through the sampling outlet at a rate of Q out S and through a needle-type valve at a rate of Q out C . (i) Two hanging drops in equilibrium with identical Laplace radius and pressure. The inflow was evenly distributed to the two drops. As long as the tip of the needle-type valve was immersed in the drop, liquid removal through the valve occurred. Consequently, (ii) the size and the Laplace pressure of the two drops became different with drop radius and pressure being inversely correlated. Due to the pressure difference, there was an increased flow toward the control drop (left) with the lower pressure. As soon as the tip of the needle-type valve was exposed to air, only air was withdrawn. The control drop remained stable at the height defined by the valve needle length. (iii) In the next step, the two drops reached equal Laplace radius and pressure through equilibration through the liquid phase. The system was now stable, and the drop height in control and tissue drop was maintained constant. Reproduced from Wu Jin et al. [10] (CC BY 4.0).
However, a phenomenon known as capillary-wetting hysteresis [19,20] caused certain irregularities in the drop height control. This phenomenon, described in detail in Figure  S3 in the Supplementary Materials section "contact-angle hysteresis", induced an inhomogeneous wetting of the needle depending on whether the drop size was increasing or decreasing. The wetting hysteresis of an advancing or receding air-liquid-solid tripoint induced more complex wetting dynamics. The wetting dynamics were modeled by including a hysteresis in the drop height (dashed red line in Figure 11) and successfully recapitulated the dynamics observed experimentally. This result showed that, although the drop height of the control drop changed by ±25 µm, as the control drop withdrew liquid, the height of the hanging drops hosting islets changed substantially less by approximately ±2.5 µm.
the drop height of the control drop changed by ±25 μm, as the control drop withdrew liquid, the height of the hanging drops hosting islets changed substantially less by approximately ±2.5 μm.  Figure 10B. The blue oscillatory pattern indicates the constant readjustment of the control drop height due to the pulsed withdrawal of excess medium through the needle-type outlet. The inflow rate was 80 μL/min, outflow rate was 15 μL/min for each drop, and the control needle outflow rate was 45 μL/min. A spontaneous wetting instance was simulated at 30 s and 60 s, where the aperture of drop 1 increased by 5% (10 μm) each time.

Conclusions
This review describes a method used to model microfluidic standing-and hangingdrop networks. The technique can be used with any open microfluidic network to analyze subtle ALI dynamics and optimize chip operation with different flow sources. We show that, with basic microfluidic components and rigorous variable and circuit definitions, intricate, open microfluidic circuits can be modeled successfully. The nonlinear dynamics of an ALI can be represented by using a behavioral model, which describes the complex dependency between drop geometry, capillary and hydrostatic pressures, and flow rate with a differential equation.
The same strategy can be applied to various microfluidic chips, for example, to predict flow behavior, to optimize chip operation, to compute shear stresses on cells inside chambers and channels, to design microfluidic channels, or to model hysteresis behavior on wetted needles.
The versatility of this method has been demonstrated in several articles and conference papers [10,11,13]. Ultimately, circuit-based open-microfluidic-network design is a "first-pass" modeling technique that predicts flow. With an optimized circuit design, more intricate 3D mathematical models can be used to characterize fluid dynamics and species transport within final designs.

Supplementary Materials:
The following supporting information can be downloaded at: www.mdpi.com/xxx/s1. Figure S1: Modeling and experimental results of the operation of the scalable microfluidic chip elaborated further in the manuscript text; Figure S2 Figure 10B. The blue oscillatory pattern indicates the constant readjustment of the control drop height due to the pulsed withdrawal of excess medium through the needle-type outlet. The inflow rate was 80 µL/min, outflow rate was 15 µL/min for each drop, and the control needle outflow rate was 45 µL/min. A spontaneous wetting instance was simulated at 30 s and 60 s, where the aperture of drop 1 increased by 5% (10 µm) each time.

Conclusions
This review describes a method used to model microfluidic standing-and hangingdrop networks. The technique can be used with any open microfluidic network to analyze subtle ALI dynamics and optimize chip operation with different flow sources. We show that, with basic microfluidic components and rigorous variable and circuit definitions, intricate, open microfluidic circuits can be modeled successfully. The nonlinear dynamics of an ALI can be represented by using a behavioral model, which describes the complex dependency between drop geometry, capillary and hydrostatic pressures, and flow rate with a differential equation.
The same strategy can be applied to various microfluidic chips, for example, to predict flow behavior, to optimize chip operation, to compute shear stresses on cells inside chambers and channels, to design microfluidic channels, or to model hysteresis behavior on wetted needles.
The versatility of this method has been demonstrated in several articles and conference papers [10,11,13]. Ultimately, circuit-based open-microfluidic-network design is a "firstpass" modeling technique that predicts flow. With an optimized circuit design, more intricate 3D mathematical models can be used to characterize fluid dynamics and species transport within final designs.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/mi13071124/s1. Figure S1: Modeling and experimental results of the operation of the scalable microfluidic chip elaborated further in the manuscript text; Figure S2: Phase diagram of hanging drops of various aperture diameters; Supplementary section on "Circuit design strategy"; Supplementary section on "Contact-angle hysteresis" with Figure S3: Advancing and receding contact angle measurements of a drop; Supplementary section on "Filling or emptying reservoirs" with Figure S4: Geometric properties of the standing-drop compartments described by Figure 6, Figure S5: Filling and emptying dynamics of a standing-drop port for a hydrophobic advancing capillary front or for a hydrophilic receding capillary front, Table S1: Calculating volumes at each filling step (advancing capillary front), Table S2: Calculating volumes at each empty step