Compact Empirical Model for Droplet Generation in a Lab-on-Chip Cytometry System

This study describes the construction of a droplet generation speed compact empirical mathematical model for a flow-focusing microfluidic droplet generator. The application case is a portable, low-cost flow cytometry system for microbiological applications, with water droplet sizes of 50-70 micrometer range and droplet generation rates of 500-1500 per second. In this study, we demonstrate that for the design of reliable microfluidic systems, the availability of an empirical model of droplet generation is a mandatory precondition that cannot be achieved by time-consuming simulations based on detailed physical models. When introducing the concept of a compact empirical model, we refer to a mathematical model that considers general theoretical estimates and describes experimental behavioral trends with a minimal set of easily measurable parameters. By interpreting the experimental results for different water- and oil-phase flow rates, we constructed a minimal 3-parameter droplet generation rate model for every fixed water flow rate by relying on submodels of the water droplet diameter and effective ellipticity. As a result, we obtained a compact model with an estimated 5-10% accuracy for the planned typical application modes. The main novelties of this study are the demonstration of the applicability of the linear approximation model for droplet diameter suppression by the oil flow rate, introduction of an effective ellipticity parameter to describe the droplet form transformation from a bullet-like shape to a spherical shape, and introduction of a machine learning correction function that could be used to fine-tune the model during the real-time operation of the system.


I. INTRODUCTION
Bacterial threats have been a noticeable challenge of this century, and a delayed response due to the lack of field-testing options poses risks to human lives and can cause epidemics. Classical microbiological methods are relatively slow, while cytometric methods allow measuring the number and morphology of cells easily, reliably, and quickly. Droplet microfluidics, a new technology developed over the last dozen years, offers breakthrough solutions for creating lowcost, fully portable cytometers for field analysis of bacteria The associate editor coordinating the review of this manuscript and approving it for publication was Kin Fong Lei .
based on very small sample volumes and the possibility of seeking single-cell resolutions.
The present study discusses model-based design of portable cytometer devices based on the concepts of labon-a-chip [1], [2], [3], microfluidics [4], [5], and droplet cytometry [6], [7], [8], [9], [10]. Specifically, we describe the construction of an empirical mathematical model for the calculation of droplet generation rates and dimensions in the water-in-oil flow-focusing-type [11] droplet generation node of a lab-on-a-chip cytometer. This study was partly based on the digital twin model developed by our group [12]. Figure 1 illustrates the topicality trends of the considered research areas based on the publication statistics of the Clarivate Web of Science database [13]. As demonstrated in Figure 1, the concept of lab-on-a-chip has become popular since the beginning of the 21st century, and is presently showing a saturation trend. The overall area of droplet microfluidics became popular slightly later in years [2003][2004] and has demonstrated linear growth until the present time. The most vital concept is droplet cytometry, for which exponential growth with a doubling time of 4-5 years started approximately 12 years ago. An overall comparison of regions over the last decade demonstrated the approximate equality of Western Europe, North America, and the People's Republic of China (PRC) [13].  [13]. Dynamics of the most relevant subfields as Lab-on-Chip, droplet microfluidics and droplet cytometry is compared. Inlet compares contributions from People's Republic of China, North America and European region (incl. Turkey and Israel).
Model-based design has become a mandatory methodology for system design in various applications, including microfluidics [14], [15], [16], [17]. In the present use case, a model of the droplet generation node is required for the prediction of droplet generation parameters, such as generation rates and diameters, both during the cytometer construction and exploitation phases, to improve the control quality via model predictive control.
An important issue in the construction of a mathematical model for any object is the selection of a detailed physical or formal empirical approach [18], [19], [20], [21]. The physical approach can be time consuming for both the computer and developer but can yield reliable results for a wide range of operation conditions, provided that the physical mechanisms and relevant key parameter values are known and modelled correctly. In the case of droplet microfluidics, physical models rely on well-known equations and methods of computational fluid mechanics (CFD), which must be supplemented with less reliable multicomponent fluid flow methods [22], [23]. Thus, in addition to the high computational workload, the major problem of the physical approach is often the presence of non-measurable physical parameters and the hidden influence of numerical factors such as reduced spatial dimensionality and mesh step sizes. Although many authors have illustrated their studies with simulated droplet images, for example, [23], [24], [25], [26], several respected research groups, for example, [24], [27], [28], [29], emphasize the unreliability of numerical physical modelling, particularly if the droplet size, generation rate, and monodispersity characteristics must all be reliably calculated simultaneously. Moreover, considering the three main droplet generation geometry types -co-axial, T-junction and flowfocused [27], [30], the third option, which is also analyzed in the present study, has been estimated most difficult for the point of view of accurate modelling [20], [30], [31]. Our simulation results with the COMSOL Multiphysics R 5.6 Two-Phase Flow Level Set module [32], described below in the simulation section, confirm the unreliability claims regarding the detailed physical modelling approach. To illustrate the relevant difficulties, it is reasonable to point out that to overcome the aforementioned uncertainties and obtain a reliable practical tool for flow-focusing droplet generator design, a largescale experimental study was recently conducted by the group of Boston University [24], [28]. In this study, a generalized flow-focusing structure with an orifice was described using six geometrical parameters: 25 orthogonal structure variants were manufactured using the Taguchi formal scheme, over 30 operation modes for each structure were tested, and a statistical empirical model was obtained to cover a reasonably wide range of droplet diameters and generation rates. In comparison with [24], [28], the present study discusses only one flow-focusing structure without nozzle (orifice) section but, on the other hand, the droplet geometry description includes also the ellipticity factor, droplet generation rates are of 2-3 times higher range and the formula-based analytic formulation having a better transparency and real-time adjustability is used.
In contrast to detailed physical models, the alternative empirical approach is characterized by a formal generalization of experimental data [20], [27], [30]. The empirical approach is usually less labor-intensive and often more accurate, but only for the parameter space covered by the experiments. In practice, the most useful real models are semi-empirical, which means that they combine theoretical principles with available experimental data. In droplet microfluidics, it is reasonable to build all droplet generation models based on the mass conservation principle for the dispersed phase (i.e., droplet fluid-like water), which allows a state connection between three main variables: droplet fluid flow rate, droplet diameter, and droplet generation rate [20], [21], [28], [31]. Some authors, who have investigated the formation of relatively large non-spherical droplets that fill all cross-sections of the generation channel, have added a fourth parameter, the droplet length, for example, [20], [21], [30]. In this study, we introduce an original approach using an effective ellipticity parameter that maintains droplet volume conservation and accounts for experimentally observed droplet shape changes from bullet-like shapes at low continuous-phase flow rates to spherical shapes at high continuous-phase flow rates.
When discussing the droplet diameter empirical models, many authors have used the ratio of dispersed and continuous phase flow rates Q d Q c , for example, [21], [25], [26]. VOLUME 10, 2022 Our droplet image recording results, presented below in the experimental section, do not support the use of this ratio parameter, and demonstrate that for higher water flow values, a proportional oil flow increase is required to achieve a comparable diameter suppression effect.
An important characteristic for the practical applicability of mathematical models is their compactness. We recommend defining compactness based on the following features:1) the minimal number of adjustment parameters, 2) measurability of the adjustment parameters, 3) low computational workload, and 4) transparency of the set of equations [33]. The concept of a compact model is widely used in the field of electronic and semiconductor microchip design [34] for two main reasons: lowering the computational workload and operation with measurable parameters. In droplet microfluidics, the need for compact models has not yet been explicitly recognized and only a few studies have used this term. However these studies have focused solely on estimating the length of droplets based on the ratio of dispersed and continuous phase flowrates [35], [36], [37], [38]. At the same time, nearly 1000 publications (see Figure 1) contain some approximate formulas for the calculation of droplet sizes or generation rates that may be interpreted as compact models for solving some subproblems of droplet microfluidic system design tasks.
An important question in droplet generation model construction is the description of the droplet diameter suppression effect owing to continuous phase (oil) flow. The majority of published results and models, for example, [25], [26], [39], predict a less-than-proportional diameter suppressing effect owing to the increasing continuous phase flow rate Q c . Few studies support either a proportional decrease in diameter, for example, [31], or a stronger than proportional increase [26]. The present experimental study confirmed the applicability of the linear approximation of the dependence of the water droplet diameter on the oil flow rate. Thus, a linear droplet diameter model may be offered that uses only one proportionality factor for a fixed water flow rate and given droplet generation channel width. If completed with two parameters for the description of effective ellipticity changes, a compact 3-parameter model for the calculation of droplet diameters, ellipticities, and generation rate dependencies on the oil flow rate may be constructed.
In recent years, there has been an urgent need to accelerate and simplify the development of microfluidic droplet generators with desired output parameters via automatization and the application of machine learning methods [24], [40], [41]. To realize these goals via empirical statistical modelling by applying artificial neural networks, large-scale experimental testing [24], [28] or sophisticated computer vision methods for additional droplet data collection [40] have been proposed. In this study, a much narrower task scope was considered and only the desired droplet parameters were obtained by adjusting the water and oil flow rates for a fixed microfluidic chip. However, formula-based transparent presentation of mathematical models offers much better possibilities for solving system optimization and real-time model adjustment (i.e., machine learning) tasks. Although the modification of neural-network-based statistical models [24], [40] requires significant effort and time for the collection of additional data and retraining (transfer learning), the empirical model considered here, in the form of mathematical formulas with adjustable coefficients, offers possibilities for the realization of real-time model adjustment and a cytometer system with an extremely simple feedback loop containing an elementary optical sensor.
The remainder of this paper is organized as follows. In Section 2, the microfluidic chip and the measurement setup are described. In Section 3, a short summary of the detailed numerical simulation results and a discussion of the problems that occurred are presented. Section 4 summarizes the experimental results for the different water and oil flow rates. In section 5, the construction principles, formulas, and fitting results of the compact mathematical model are presented. Section 6 discusses the scope of the application of the proposed model. Section 7 presents the main results of the study.

II. DESCRIPTION OF DROPLET GENERATION CHIP AND MEASUREMENT SETUP
Microfluidic droplets were generated inside a polydimethylsiloxane (PDMS) chip, as shown in Figure 2. The full thickness of the PDMS chip was 5 mm, and it had a microfluidic channel structure with a depth of 100 µm on one surface (Figure 2a), which was covered by a 1 mm thick glass plate (microscope slide plate). From the three main droplet generation geometries, the T-junction, co-flow, and flow-focusing junction [11], the last geometry variant, where water flows with biological agents, is cut into droplets by a continuous oil flow entering the junction area from the two opposite sides (Figure 2a and 2c). Thus, water droplets were formed in the junction area and in the generation channel with cross-sectional dimensions of 84 µm width and 100 µm height. An overview of the droplet generation unit with the inlet and outlet tubes is shown in Figure 2b. Deionized water was used as the dispersed phase (droplets). For the continuous phase, Sigma-Aldrich 330779 mineral oil [42] with a 2% surfactant [43] was used. The water and oil flow rates were maintained using syringe pumps and software manufactured by SpinSplit [44]. The lighting of the droplet generation junction area was realized from the PDMS side of the chip using a white LED group consisting of two LEDs with cold-color temperatures and two LEDs with warmcolor temperatures. Photorecording was accomplished using a Basler Ace acA640-750um camera in a reduced resolution mode that allowed a frame rate of up to 3300 per second at an exposure time of 100 µs. Thus, as the experiment shows, a droplet per second (dps) generation rate of up to 1600 s −1 can be directly determined from the sequence of the recorded images. Additionally, dps values up to 2300 s −1 can be extrapolated based on the droplet separation distances (see Figure 7).

III. INTRODUCTORY SIMULATIONS OF UNDERLYING PHYSICS
In general, the prerequisite for the construction of a compact model may be the availability of experimental data or, as an alternative, the availability of a sufficiently reliable detailed physical model with necessary input data. In the case of flow-focused droplet generation junctions, the choice of the detailed physical approach can be complicated by the complex nature of the task, that is, the need for accurate modelling of the balance of competing processes of separation and encapsulation of droplets. Another serious problem is the high computational time required for realistic three-dimensional calculations. Therefore, all affordable two-dimensional calculations, even if the parameters of the physical processes are correctly estimated, can only serve as predictions that need to be confirmed by real experiments. Figure 3 shows the critical competing processes that must be accurately modelled in a flow-focused junction. Figure 3 emphasizes the importance of accurately modelling the surface tension forces, viscosities of both liquids, wall friction effects, channel dimensions, and other factors to obtain a realistic picture of both the liquid flow and droplet formation processes.
To test the possibility of using detailed physical modelling to formulate the basis of the droplet generation model, we performed several numerical simulation series using the COMSOL Multiphysics R 5.6 Two-Phase Flow Level Set Illustration of competing processes of water droplets generation in a flow-focusing cross-junction. The incoming water stream tries to maintain the minimum surface area due to the surface tension forces but is divided into droplets by the ''oil pliers'' acting from both sides. After that the surface tension helps to maintain the size of droplets already formed, provided that the adjacent droplets are at a sufficient distance. At that all flow speeds are decelerated near the walls because of the wall friction effect.
module [32] in the traditional two-dimensional (2D) axisymmetric approximation of geometry [11]. It is important to emphasize that the crucial point for the accuracy of all modelling approaches is the correct handling of water volumes in the task specification. Oil can be considered an auxiliary substance that splits the incoming water stream and suppresses the diameters of the formed water droplets. Because in the 2D-simulation the droplets are cylindrical rather than spherical, the first question in the specification of the 2D-simulation task is to correctly select the effective size of the simulated structure towards the third dimension. Considering the realistic situation of the 3D-experiment (at high oil flow rates), it can be assumed that the droplets are spheres with a volume where D exp is the droplet diameter used in this experiment. In the 2D simulations, the volume of the droplet was defined using the cylinder formula where H eff is the introduced effective size of the structure towards the third dimension (see Figure 4). The diameter and volume of the droplets must be equal to match the water volume counts in the experiments and simulations. Figure 4 shows the methodology for achieving the aforementioned water volume balance conditions when the auxiliary parameter H eff of the 2D-simulation is specified as follows: Specifically, if the actual expected droplet diameter is 60 µm (in the 84 µm channel), then a reasonable measure for the structure depth in 2D-simulations should be 40 µm. Figure 5 summarizes the main results of the COMSOL 2D-simulations with the Two-Phase Flow Level Set module [32] for the droplet generation area described in Figure 2c. For the first adjustment parameter, the effective depth of the structure in the third dimension was specified as 40 µm based on the considerations explained above. For the second essential adjustment parameter, the surface tension coefficient values of σ = 40 ÷ 50mN/m were used to avoid the jetting effect and ensure the acceptable stability of the formed droplets. High surface tension values in a similar range have been recommended for water-tomineral oil interfaces, for example, in [45] and [46]. Next, in presented simulations the ''no-slip'' sub model of high friction walls was used. Other models, such as the Navier slip model with several additional adjustment parameters, did not cause essential changes.
For the main computational parameters, that is, the spatial mesh size, the two standard cases of ''Fine'' with 9536 finite elements and ''Finer'' with 36626 elements were compared. The computational times for the relatively short 20 ms process calculation ranged from 2 h to 14 h on a powerful desktop 16-core Intel i9-computer.
The main results of the COMSOL simulations are shown in Figure 5. The results demonstrate the difficulty of achieving stable droplet diameters and droplet generation rates. When the spatial mesh size was increased to a very high number of final elements, instead of the expected stabilization of the main output parameters, the chaotic behavior of the results demonstrated a remarkable increase, and the definition of certain values of droplet diameters and droplet generation rates became impossible. This emerging instability and uncertainty may be caused by the difficulty of the task, as shown in Figure 3. In summary, detailed physics-based numerical simulations provide supporting explanations for the underlying physical processes. However, the expected results for building a compact model for droplet generation have not been obtained.

IV. SUMMARY OF EXPERIMENTAL RESULTS
For the actual microfluidic chip, the expected droplet generation rates were in the range of 500-1500 per second, with droplet diameters of 50-70 micrometer range. Based on these design goals, three test series with constant water flow rates Q w = 4, 8, and 12µL/min and varying oil flow rates from value Q oil = 2Q w to value Q oil = 60 ÷ 88µL/min in steps of 4 µL/min were performed. The selection of droplet images recorded in the beginning section of the 84 µm generation channel is presented in Figure 6. The results in Figure 6 show that at low oil flow rates, the droplets resemble bullets (modelled by the effective ellipsoids below in this study). With an increase in the oil flow rate, the droplets begin to resemble spheres. Simultaneously, the diameter can be suppressed by increasing the oil flow rate. The increasing blurring of droplet fronts and backs at higher droplet formation rates is due to camera exposure settings (100 µs).
An overall summary of the experimental results, including the directly recorded and extrapolated dps values from the droplet separation distances, is shown in Figure 7. The experimental diameters of droplet D were obtained by carefully comparing the droplet lateral sizes with a channel width of 84 µm and smoothing the dependencies with the neighboring points. Thus, the estimated accuracy of the diameters was of the order of ±2 µm. The effective ellipticity numbers, E (approximate ratio of the vertical and horizontal sizes of the droplets in the images), were estimated from the principle of equivalent volumes of the imaginary ellipsoidal and real bullet-like droplets. Additionally, minor smoothing of the experimental diameter and effective ellipticity values was performed to ensure correlation with the real water flow rates.
Owing to the camera frame rate limit, the high droplet rate values over dps > 1600 s −1 were difficult to define from video recordings but were extrapolated on the basis of the observed decrease in droplet distances (see inlet in the upper part of Figure 7). Regarding this extrapolation, it should be mentioned that because of the increasing influence of the friction of the channel walls at higher Q oil rates, the size of the effective high-flow-speed center area of the channel may be smaller at high Q oil values; thus, the extrapolated dps numbers may be underestimated.

V. CONSTRUCTION OF THE COMPACT EMPIRICAL MODEL
In the present microfluidic system design, the main purpose of developing a compact droplet generation model is to obtain a tool for estimating the droplet generation rate dps. The latter depends directly on the droplet volume estimation by the sub-models for the droplet size parameters, such as the diameter D and effective ellipticity E, if the water flow rate Q w is given. The oil flow rate Q oil can be interpreted as an auxiliary factor that suppresses D compact model and can be constructed based on the following approximations: 1) For droplet generation rate dps, recalculation from a single ellipsoidal droplet volume formula V = E (π /6) D 3 can be applied if the diameter and ellipticity are estimated with reasonable accuracy.
2) Initially, it is reasonable to consider all three water flow rate values Q w,i separately. The final result for any Q w value can be interpolated based on three separate results for dps i , D i , and E i .
3) Relying on Figure 7 and seeking the principle of minimal complexity, for diameter dependence on oil flow rate D(Q oil ) the simplest single-parameter linear dependences may be applied; for the zero-oil origin point, the actual channel width value of 84 µm can be used as a common constant. The changes in dps in the 10% range were acceptable for an approximate adjustable model. 4) For E, the decreasing exponent law can be applied with a final level at high oil rate values close to one, which corresponds to the spherical limit form. 5) For machine learning readiness one real-time adjustable correction function C ML (Q w , Q oil ) may be added.
Considering the principles described above, the following set of mathematical equations can be proposed for the compact model (for every water flow rate value Q w,i ,i= 1, 2, 3): where, following the goal of minimizing the number of adjustable parameters, only three fitting parameters, E 8,i , Q E,i , Q D,i were introduced for each of the tested water VOLUME 10, 2022  flow rate values: Q w,1 = 4 µL min Q w,2 = 8 µL min, and Q w,3 = 12 µL min (see Figure 6). In systems (4)-(6), equation (4) is constructed to transform the value of the water flow rate to the number of droplets per second, considering the lateral diameter of droplets D i and the effective ellipticity E i as key parameters for the calculation of a single droplet volume. Equation (5) postulates the simplest linear decrease law for droplet diameters by introducing only one adjustable parameter Q D,i for every water rate, and using an actual channel width of 84 µm as a fixed constant for the low oil flow limit. Equation (6) approximates the exponential decrease in effective ellipticity from the initial high value at Q oil = 8 µL/min to the final unit value using two adjustment parameters: E 8,i and Q E,i .
For the general case of any water flow rate between 4 and 12 µL/min, the simplest reliable piecewise linear approximation may be offered, considering that higher-order approximations such as parabolic approximations may distort the monotony of the dependences. In addition, an advanced feature of machine learning readiness may be included in the real-time empirical adjustment function C ML (Q w , Q oil ) for dps. In the minimal model formulation, the droplet size parameters may be excluded from the real-time adjustment because they droplet size parameters are difficult to measure during real-time operation.
Thus, the piecewise linear interpolation-based generalization of the droplet generation rate calculation for any water rate value can be performed using equations (7) and (8) given below.
The mathematical formulation of the linear approximation with machine learning adjustment for the first interval Q w,1 ≤Q w ≤ Q w,2 can be written as and for the second interval Q w,2 ≤ Q w ≤ Q w, 3 as 3 , The fitting of the three parameters of models (4)-(6) to determine the best agreement between the model and experimental points in Figure 7 was performed by separately minimizing the root-mean-square (RMS) difference between the experiment and simulation for the three water rate values. The weight scalers for the three main output parameters, dps, D, E were 20 s −1 , 2 µm, and 0.1, respectively. In addition, the weights of the low oil rate endpoints for the high water rate curves Q w,2 , Q w,3 were increased to obtain a reasonable balance with the low water Q w,1 curve. The overall results of the fitting are shown in Figure 8. The values obtained for the model coefficients are listed in Table 1.

VI. DISCUSSION AND APPLICATION AREA
The definition of the application area is an important issue in empirical models. The general principle is that the reliability of the results can only be assumed in the range of parameter values covered by the experimental results. The application area of the proposed compact empirical model is illustrated in Figure 9. In addition, we compared our experimental results with existing droplet length estimation models and found that existing models are also suitable for our microfluidic setup after some modifications to microfluidic chipdependent parameters [35], [36], [37], [38]. It is important to emphasize that previous studies have solely focused on estimating the length of droplets.
As shown in Figure 8, the droplet rate calculation accuracy of the proposed compact model remained in the range of 20% when considering the all-parameter area. However, it is   important to emphasize that the trends of changes due to water and oil flow rate changes were modelled correctly. Additionally, for the central region of the planned operation around dps = 500-1500 s −1 the accuracy is much better and is already in the 5-10% range. Moreover, this number can be improved by machine-learning adjustments during real operations if the droplet generation rates are measured using optical measurements.
The reason for the moderate accuracy of the proposed simple 3-parameter model is the simplicity of modelling the droplet diameter using the simplest 1-parameter linear dependence. Since the droplet formation rate depends on the droplet diameter according to the cubic law, small differences in diameters of about 3% were increased to 10% when the formation rates were taken into account. It is possible to introduce a sub-model of a more precise diameter; however, the accompanying increase in the number of model parameters may require additional complex measurements.
From a technical viewpoint, it seems more reasonable to use a simpler model with the possibility of real-time adjustment.

VII. CONCLUSION
Compact models are a well-established approach in electronics and microelectronics but are not yet sufficiently appreciated in the relatively young field of droplet microfluidics. Although for design of any technical system, like droplet cytometry portable apparatus in the present use case, the availability of compact models for all subsystems is a highly VOLUME 10, 2022 desirable precondition for the successful design of the system as well for later exploitation of the system. In addition, as demonstrated in the present study, the alternative approach of detailed physical modelling may not yield usable results in the case of the droplet generation task of microfluidics, where the competing balances of different physical mechanisms must be accurately modelled.
The original new results presented in the present study may be summarized as follows: 1) For the first time, a compact empirical model of droplet generation speed has been developed.
2) The applicability of the linear approximation of the dependence of the droplet diameter on the oil flow rate for the actual flow-focusing microfluidic water droplet generator task (droplet sizes in the 50-70 micrometer range and generation rates in the 500-1500 per second range) was demonstrated.
3) The concept of effective ellipticity was introduced to describe a unified model for the change of droplet geometry from a bullet-like to a spherical shape. 4) The methodology for the construction of a minimized compact 3-parameter droplet generation rate model with 5-10% accuracy for the calculation of the oil flow rate dependence at fixed water flow rates for the desired operation region was described. 5) A machine learning extension to the basic model for further adjustment using real-time measurement results was proposed. 6) The droplet-volume-based equivalence condition to make 2D-simulations comparable to the real 3D experimental geometry is discussed.