1.1 Polymerase chain reaction

The polymerase chain reaction (PCR) is a highly utilized technique in medical genetics and various fields of biological science. It involves amplifying specific deoxyribonucleic acid (DNA) fragments (amplicons) by subjecting a reaction mixture of nucleic acids to thermal cycling. The mixture consists of oligonucleotide primers complementary to the DNA sequences of interest, deoxynucleotide triphosphates (dNTPs), DNA polymerases, and DNA samples containing the target sequence. The core cycle of the PCR includes 1. denaturation, 2. annealing, and 3. extension, which is repeated multiple times to achieve amplification of the target DNA sequence ^{[14 - 16]} (Fig. 1). The denaturation involves heating the mixture (95) to separate (denaturation) the double-stranded DNA (dsDNA or S1S2) into two separate single-stranded DNA (ssDNA o S1 and S2) strands. During the annealing, the mixture is cooled (55) to bind the oligonucleotide primers (P1 and P2) to the complementary DNA sequences on the ssDNA, forming stable primer-ssDNA complexes (S1P2 and S2P1). Finally, the extension includes raising the temperature to 75, where the DNA polymerases catalyze the synthesis of the complementary DNA strands by extending the primers using deoxynucleotide triphosphates (dNTPs) as building blocks. Natural convection PCR devices exhibit continuous denaturation and renaturation as reactants circulate through distinct reaction zones. Consequently, the concentration dynamics of ssDNA (S1 and S2) and dsDNA (S1S2) in our simulations are anticipated to be similar to the fluorescence intensity behavior observed in quantitative PCR (qPCR) tests.

Figure 1 The PCR reaction involves amplifying specific amplicons by subjecting a reaction of deoxynucleotide triphosphates, DNA polymerases, and DNA samples containing the target to thermal cycles. The sequential cycles of denaturation (heating the mixture to denature the double-stranded DNA or ssDNA), annealing (cooling of the mixture to bind the primers to form stable primer-ssDNA complexes), and extension (synthesis of the complementary DNA strands by extending the primers) are carried out until the amplification of the target DNA sequence is reached. Adapted image from ^[17].  

Quantitative PCR (qPCR) enables the measurement of amplicon accumulation over time by monitoring the fluorescence intensity generated by fluorophores incorporated into the reaction. The fluorescence intensity is directly proportional to the concentration of the amplified DNA. The behavior of the fluorescence intensity during qPCR typically exhibits four distinct phases (Fig. 2), which include an initial stage (basal phase) with minimal fluorescence (baseline signal) before target DNA amplification begins ^[18]. In the second stage (exponential phase), the fluorescence intensity increases exponentially as target DNA amplification progresses rapidly. The linear phase (third stage) involves a steady and linear increase in fluorescence intensity, indicating a constant rate of amplicon generation. The final stage corresponds to the plateau phase, where the fluorescence intensity reaches a maximum, signaling the end point of amplification.

Figure 2 The qPCR technique allows the measurement of the concentration of the amplified DNA during a PCR experiment by measuring the fluorescence emitted during the sequential cycles of the process. The first of the four phases is the basal one, corresponding to the signal before target DNA amplification begins. The exponential signal (second phase) is due to the rapid progress of fluorescence intensity of target DNA amplification. The linear phase (third stage) indicates a constant rate of amplicon generation. The plateau (final phase), with fluorescence maximum intensity, corresponds to the end point of amplification. Adapted image from ^[17].  

This work proposes a methodology framework for simulating and optimizing reaction kinetics in natural convection microfluidic devices using computer simulation and a multi-physics approach. The approach involves coupling heat and mass transfer, fluid flow, and chemistry of the reactants. Visual and regression analyses are performed to evaluate the impact of different operational parameters on reaction speed, aiming to improve microfluidic natural convection systems. The proposed methodology is applied to a PCR microfluidic device that utilizes natural convection for the required reactions.

2.1 Methodology software

COMSOL Multiphysics (ver. 5.3.) was used to simulate a PCR microfluidic device. It is a finite element-based simulation platform capable of simulating fully coupled multi-physics and single-physics phenomena. Nevertheless, other multi-physics simulation software could be used similarly.

2.2 Operational parameters and experiment design

Using pumps is not feasible since we are dealing with natural convection devices. Therefore, a constant temperature gradient is generated to induce a mass flux within the device. This temperature gradient leads to density variations in the fluid, leading to a mass flux driven by the force of gravity. Based on this principle, the operational parameters required to generate a mass flux and promote adequate mixing and reactions among different chemical species include the inner diameters of the devices and fixed temperatures at various locations within the device. By manipulating these parameters, optimal conditions for efficient and accelerated mixing and reactions can be determined. Furthermore, the design and materials of the device can be considered as additional parameters, treated as categorical variables, to evaluate their impact through regression analysis. This type of analysis aids in selecting the appropriate material and design for fabricating the device. It is recommended to establish a complete or fractional factorial design for the simulations to investigate the impact of different parameters on the rate of chemical reactions through a regression analysis. The results of the simulations and the operational parameters will be used to perform visual and regression analyses to evaluate the individual effects of each parameter.

The application of the methodology in this work focuses on a PCR device that relies on natural convection. The operational parameters for this PCR device include three specific temperatures: 1) the denaturing temperature (TD), 2) the annealing temperature (TA), and 3) the extension temperature (TE). The values were selected to cover the typical experimental PCR ranks. Additionally, the inner diameter (ID) of the device is considered an operational variable, and the corresponding values for all parameters can be found in Table I. A complete 3⁴-factorial design was used, considering the 81 possible combinations of the parameters listed in Table I, to measure the influence of the parameter on the rate of PCR reactions.

2.3 Device design and meshing

The proposed PCR device design was imported to COMSOL. It was a capillary glass tube shaped like a closed triangular loop with a length of 64 mm [Fig. 3a)]. The CAD models of the device were created using AutoCAD. The diameters are specified in Table I. The device is divided into two domains [Fig. 3b)]. The inner section of the tube was filled with a liquid (in this case, water), where the chemical species will be present and undergo reactions. In contrast, the outer section of the tube is a solid material (in this case, borosilicate glass).

Figure 3 a) Isothermal temperature constraints: denaturation zone (TD), annealing zone (TA), and extension zone (TE). b) Domains of the PCR system: fluid and walls. 

An automatic finer mesh was applied to the imported models into COMSOL. The mesh selection was based on the consideration that larger meshes resulted in warnings or errors within the program. In contrast, smaller meshes did not yield significant changes in the simulations. The properties of the three meshes generated can be found in Table II.

2.4 Mathematical models

2.5 Simulation process

The properties of the different species, including initial concentrations and diffusivity constants, are defined in Table IV.

The simulation process is divided into two primary parts: 1) a steady-state analysis of the termofluidic fields within the device and 2) a time-dependent analysis of the transport and chemical reactions among the different chemical species. Separating simulations into stationary and time-dependent parts simplifies computations and reduces the required computational resources and time. All the mathematical models treated need to be configured in the software along with the boundary and initial conditions. The first study involves calculating the velocity and temperature fields within the device. This was achieved by coupling the Creeping flow, Heat transfer in fluids, and Non-isothermal flow modules. A parametric sweep can be configured to conduct simulations with varying combinations of operational parameters. The average fluid velocity (AvgVel) was computed for each simulation. The second study involves calculating the reaction kinetics of the reactions over time, ranging from 0 to 60 minutes with 1-minute steps. These values could be modified depending on the computing capacities and the accuracy needed. It was achieved by coupling the chemistry and transport of diluted species modules and using the velocity and temperature fields obtained from the first study to calculate the reactions among the chemical species and their corresponding concentrations at each time step. For the PCR device example, the simulations were performed using the complete factorial design with the parameters in Table IV.

2.6 Processing and analysis of the simulation results

The simulations provide valuable information, including the total reaction time (RTIME). However, in the case of the PCR reactions presented as a practical example, we can also calculate additional variables of interest. These variables include the duration of the exponential phase of the reactions (EXPTIME) and the slope of the logarithm of the normalized concentration of S1S2 (SlopeLog-dsDNA) or S1 and S2 (SlopeLog-ssDNA) during the exponential phase. These additional variables are significant as they indicate the suitability of the parameters for one phase, which can extend to subsequent phases. Moreover, the operational parameters, the temperature differences TD-TA, TD-TE, and TE-TA, and the average fluid velocity (AvgVel) are also variables of interest. Once the variables of interest have been collected, a visual analysis of data allows comparison of the different variables to identify those that exhibit a higher degree of explanatory power for a regression analysis. By selecting the most influential variables, we can effectively explain the speed of the reactions. Then, a regression analysis will be conducted to describe the speed of the reactions and quantify the influence of each parameter on the reaction rate.

2.7 Computer specifications

Simulations were performed on a 2019 Acer Predator Helios 300 laptop with 16GB of RAM, an Intel i7-9750H processor, Windows 11 home OS, NVidia GTX 1660Ti graphics card, and x64 architecture.

3.1 Simulation results

The methodology implemented can predict how temperature gradients in natural convection devices affect fluid motion and reaction kinetics in a given problem. We used a PCR to exemplify the application of a natural convection device, with 81 simulations conducted. Table V presents the TD, TA, TE, and ID parameters and AvgVel, Slope_Log_dsDNA, EXPTIME, and RTIME calculations obtained from the 50 simulations that converged. TD, TA, TE, and ID were selected between the lower and higher values considered in an experimental PCR. All stationary analyses converged except for the time-dependent analyses, where sudden negative and infinite concentrations appeared, possibly due to internal calculation errors like rounding decimals, mesh quality, or the magnitude of time steps. The obtained reaction times ranged from 18 to 21 minutes, notably shorter than traditional PCR devices, where reaction times typically range from 30 to 40 minutes ^[8].

Figure 4 illustrates the velocity and temperature profiles within the PCR device for simulation SIM78, with TD = 98, TA = 65, TE = 74, and ID = 797.5 μm. The velocity profile [Fig. 4a)] shows a parabolic-like behavior, satisfying the no-slip boundary condition. The temperature profile [Fig. 4b)] indicates the fluid flow direction, which extends from the denaturation zone (TD zone) to the extension zone (TE zone) and ultimately to the annealing zone (TA zone). These characteristics hold for the remaining stem simulation as well. However, the current direction of fluid flow is a problem as it slows down the reaction rate. Ideally, the fluid should flow from the denaturation zone (TD zone) to the annealing zone (TA zone) and then to the extension zone (TE zone). This directional issue is due to the geometry of the system and the positioning of the temperature constraints.

Figure 4 a) Velocity and b) temperature profile from the simulation SIM78 (ID = 797.5 μm, TD = 98, TA = 65, and TE = 72). 

Figures 5 show the reaction kinetics in two simulations, SIM3 [Figs. 5a)], ID = 797.5 μm, TD = 92, TA = 55, and TE = 70) and SIM78 [Fig. 5b), ID = 797.5 μm, TD = 98, TA = 65 , and TE = 72 ], presenting continuous processes of denaturation and renaturation. This characteristic holds for the remaining simulations as well. For SIM3, the concentration of dsDNA (S1S2) and ssDNA (S1 and S2) were similar to the fluorescence intensity observed in qPCR tests, whereas SIM78 showed denaturation and renaturation rates exhibiting variations over a time interval with a decrease in dsDNA (S1S2) concentration at certain time intervals. It can be attributed to a fast amplification rate that exceeded the denaturation rate until the reactants were sufficiently depleted to reduce the amplification rate, allowing the denaturation process to enter a more pronounced plateau phase.

Figure 5 Kinetics of the PCR reactions of the chemical species of selected simulations. a) SIM3 with parameters ID = 797.5 μm, TD = 92, TA = 55, and TE = 70 and b) SIM78 with parameters ID = 797.5 μm, TD = 98, TA = 65, and TE = 72. 

3.2 Target variable selection

Minimizing the total reaction time (RTIME) is needed to increase specific reaction rates. Other significant variables in PCR reactions include the duration of the exponential phase (EXPTIME) and the slope of the logarithm of the normalized concentration of dsDNA (S1S2) during the exponential phase (SlopeLog-dsDNA).

These variables provide insights into the suitability of parameters for one phase, which can extend to subsequent phases (Fig. 6). Figure 6a) displays the positive correlation between EXPTIME and RTIME, while Fig. 6b) exhibits the negative correlation between EXPTIME and SlopeLog-dsDNA. Then, there is also a negative correlation between RTIME and SlopeLog-dsDNA. Consequently, it is necessary to maximize Slope_Log_dsDNA or minimize RTIME or EXPTIME. A regression analysis was performed after selecting one of these variables as a response to the system. A histogram of the three mentioned variables was plotted to make the proper selection [Fig. 6c) for RTIME, 6d) for EXPTIME, and Fig. 6e) Slope_Log_dsDNA]. EXPTIME was chosen because it is closer to a constant distribution than the other two variables, which is useful when a simple regression is applied.

Figure 6 Correlation between response variables of the PCR system: a) EXPTIME vs RTIME and b) EXPTIME vs Slope_Log_dsDNA. Frequency distribution of the response variables: c) RTIME, d) EXPTIME, and e) Slope_Log_dsDNA.  

4.1 Regression analysis

The (TD-TA), (TE-TA), and ID data were standardized for the regression analysis by subtracting the mean and dividing by the standard deviation to visualize the regression residuals. The resulting model corresponds to Eq. (15)

where the temperature difference (TD-TA) is the most significant parameter, followed by the inner diameter ID and the temperature difference (TE-TA), respectively. To minimize EXPTIME, TA temperature should be decreased, thereby simultaneously increasing (TD-TA) and (TE-TA) temperature differences. Also, we have the flexibility to manipulate the TD and TE temperatures for better control over the process.

Figure 10a) shows a scatterplot of the residuals with no discernible pattern indicating any regression issues and, for instance, a reliable regression model. Figure 10b) shows a histogram of the residuals, which shows a closer normal shape distribution, considered acceptable in the analysis context.

Figure 10 Residuals distribution of the EXPTIME vs (TD-TA), (TE-TA), and ID regression. a) Scatterplot of the regression residuals and b) frequency distribution of the residuals. 

Nevertheless, the natural convection PCR has certain limitations, especially regarding the TD temperature. Going beyond 98 can result in the degradation of chemical components in real-life applications. However, we can still optimize the reaction rates by adjusting the parameters ID, TE, and TA. Another approach is reducing the heating area size (TD zone) to minimize the time available for potential chemical degradation, potentially allowing it to surpass the 98 limit. Despite the challenges of the fluid flow direction and the rates of renaturation and reverse annealing, the regression analysis allows us to quantify the influence of several parameters in the system. Then, through this analysis, we can get insights into mitigating these undesirable effects and enhancing the overall system performance. Moreover, the developed model helps to optimize the PCR natural convection system. However, further enhancements can be achieved by incorporating multiple designs or materials as parameters (categorical variables) to quantify the benefits of different designs or materials and make informed decisions based on their effectiveness.