Optimisation of microfluidic polymerase chain reaction devices

The invention and development of Polymerase Chain Reaction (PCR) technology have revolutionised molecular biology and molecular diagnostics. There is an urgent need to optimise the performance of these devices while reducing the total construction and operation costs. This study proposes a CFD-enabled optimisation methodology for continuous flow (CF) PCR devices with serpentine-channel structure, which enables the optimisation of DNA amplification efficiency and pressure drop to be explored while varying the width (W) and height (H) of the microfluidic (μ) channel. This is achieved by using a surrogate-enabled optimisation approach accounting for the geometrical features of a μCFPCR device by performing a series of simulations using COMSOL Multiphysics 5.4. The values of the objectives are extracted from the CFD solutions, and the response surfaces are created using polyharmonic splines. Genetic algorithms are then used to locate the optimum design parameters. The results indicate that there is the possibility of improving the DNA concentration and the pressure drop in a PCR cycle by ∼2.1 % ([W, H] = [400 μm, 50 μm]) and ∼95.2 % ([W, H] = [400 μm, 80 μm]) respectively, by modifying its geometry.


Introduction
PCR is a molecular technique that enables the fast amplification of a specific DNA segment and has revolutionised biological science and diagnostics [1][2][3]. Following a series of PCR cycles, this technique can create a large number of copies of the target DNA fragment, allowing the detection and identification of gene sequences using appropriate optical techniques [2].
A PCR cycle consists of three stages; denaturation, annealing, and extension. Denaturation takes place at ∼368 K, where the two-stranded DNA molecules denature into pairs of single-stranded ones. Then, the sample enters the annealing stage at ∼328 K, where the primers form primer-template complexes. In the next stage, the temperature is increased to ∼345 K where the polymerase binds to a primer-template complex [4]. After repeating several PCR cycles, DNA is amplified significantly [5]. The efficiency of the reaction depends on the number of PCR cycles performed and the temperature uniformity of the samples [6][7][8].
This molecular technique has provided scientists with the ability to work with "raw" samples (degraded templates, blood, tissue, individual hairs, etc.) [4] and has a variety of clinical applications [1]. Since its invention, several devices that perform PCR have been developed. The first PCR systems were based on devices where the entire chip would be heated and cooled appropriately during the denaturation, annealing, and extension steps of PCR [9]. Many designs developed later utilised microfluidic (μ) technology, considering their significant advantages compared to PCR macroscopic devices (fast heat transfer [10], low * Corresponding author: fotizagl@hotmail.com thermal mass and thermal inertia [11], reduced operation cost due to the small amounts of samples and reagent used [11]).
A great number of PCR applications in microfluidic devices can be found in literature, such as the publications of [12 -17]. Despite the broad use of these devices, μCF-PCR devices present some limitations. The PCR mixture experiences adsorption phenomena at the flow channel interface, which leads to PCR inhibition and carryover contamination, reducing the yield of the reaction. Also, the large channel surface area to sample volume ratio enhances the adsorption of biological/chemical particles. Another drawback is the variation in the dwell times of PCR mixture (PCR mixture moves faster in the channel's centre than it does close to the surface) that increases in the total residence time in the device [18].
Due to these limitations, recent research has also considered the development of droplet-based μPCR devices (DR-PCR). The droplets are characterised by temperature uniformity due to their small size, while at the same time act as separate chemical reactors, providing high reproducibility of reaction conditions. Furthermore, the droplets provide a confined environment, preventing contamination of the samples and any adsorption phenomena at the surfaces of the channel [18]. Detailed descriptions of DR-PCR devices can be found in the publications of [19][20][21].
Considering the complexity and cost associated with the fabrication of some parts of droplet-based devices [18], this work will focus on research related to the optimisation of SPCFPCR devices, to increase their sensitivity, specificity and their ability of multiplexing (amplifying more than one target sequence by using more than one pair of primers [22]) [23].
Several studies have explored the effect of geometrical or operating variables (e.g. thermal conductivity of the substrate materials, channel sizes and spacing, flow rate, residence time in each PCR stage, arrangement of the heaters [24][25]) may have on the PCR efficiency and the thermal cycling. [26] for example studied the effect that residence time can have on DNA amplification, while [26,7] used models of PCR kinetics, heat, and fluid flow transfer to investigate and increase DNA amplification in CFPCR devices.
This paper uses outputs from Computational Fluid Dynamics (CFD) analyses to perform simulation-based optimisation. This approach is commonly adapted in industries, in order to optimise complex flow systems [27]. The number of design variables is considered important; with advanced adjoint and gradient-free surrogate-assisted (such as Gaussian Process Emulators [28] and Moving Least Squares [29]) methods being applied for problems with ~ 1000s and < 100 design variables respectively [30]. Machine learning and surrogate modelling can also be used for achieving temperature control in CFPCR systems [31][32].
The present study is one of the first [31,33] to apply a systematic, CFD-enabled optimisation methodology of the geometrical parameters of μCFPCR devices (such as the one considered recently by [7]) with serpentinechannel structure, that enables the optimisation of DNA amplification efficiency and pressure drop to be explored, through the development of meta-models. Such understanding will enable engineers to maximise DNA amplification and minimise pressure drop, by optimising the geometrical features of a μCFPCR device through a series of simulations. In these simulations, the flow field, heat transfer, and PCR kinetics are included, providing a realistic representation of the performance of one PCR cycle on a channel with a serpentine structure.

Problem Specification
This work is motivated by the publications of [7,31,33], and it simulates the fluid flow, heat transfer, and PCR kinetics that take place in a microfluidic channel where fluid is passing through. The microfluidic channel is characterised by a rectangular cross-section and a serpentine-like structure, while the fluid is considered to have the thermal and fluid physical properties of water (see Figure 1 of [7] and Figure 1). The fluid temperature varies along the microchannel, as it flows through three different temperature zones (368.15, 328.15, and 345.15 K for denaturation, annealing, and extension respectively). The three temperature zones are created using copper wire heaters, and the fluid stays in each temperature zone for a specific residence time, in order to perform successfully one PCR cycle ( Figure 1). This variation in the temperature along the microchannel is responsible for the increase in the DNA concentration by the time the sample exits the microfluidic channel. The substrate material (Kapton, PDMS, and PE) properties are presented in Table 2 of [7], while their heights (HKapton, HPDMS, and HPE) are set equal to 100, 50, and 50 μm respectively [7]. The average inlet (fully developed) velocity is calculated by Equation 1, while the length of each zone is calculated by Equation 7. The width of the microchannel at the extension zone is twice the size of the other temperature regimes, while the residence times in denaturation, annealing, and extension are 3.0s:4.2s:6.2s respectively [7].  [7], B) design case with maximum DNA amplification, C) design case with the minimum pressure drop. The boundary conditions are presented in Figure 1A

Flow Modelling
Navier Stokes equations are used in our model to describe the motion of fluids [36]. In order to determine the type of flow, an indicative Reynolds number is calculated for a temperature of 345.15 K, for the design case of Qvol = 3·10 −11 m 3 /s, H = 50 µm and W = 100 µm (for the fluid properties of water [7]). According to Equations 1 and 2, the flow in the microfluidic channel is considered to be laminar, since the value of indicative the Reynolds number is found equal to 0.35. (2) where ρ: the fluid density, Uin: the inlet velocity, and µ: the fluid viscosity.

Heat Transfer Modelling
Heat transfer is modelled in steady-state, as presented in Equation 3. The LHS and the first term of the RHS of Equation 3 describe the convective and conductive heat transfer terms respectively. The velocity field (u) and the convective heat transfer term are only non-zero at the fluid domain. The heat generation rate of each j th (j = {1, 2, 3}) heater is described by the second term of the RHS of Equation 3, ℎ , . The different heat generation rates at the heaters are required to accomplish the different target temperatures for the denaturation, annealing, and extension zones. The third and fourth terms of the RHS of Equation 3 describe the heat flux due to thermal radiation for each of the solid substrates (i = {Copper, PDMS, PE, Kapton}) (Equation 4) and the heat losses to the ambient (Equation 5) respectively.
Qrad,i = εiσ(Tamb 4 − T 4 ) (4) Qnat.conv = h(Tamb − T) (5) where Tamb: the ambient temperature, εi: surface emissivity for the i th solid, σ: the Stefan-Boltzmann constant, and h: heat transfer coefficient. A periodic boundary condition is set for the temperature at the inlet and outlet of the unitcell (ensuring that the temperature at the inlet of the n th unitcell is equal to the one at the exit of the (n -1) th cycle). Periodic boundary conditions are also set on the outer sides of the unitcell (zero temperature offset), ensuring that the simulated unitcell is not placed at a corner/side location of the entire PCR device ( Figure  1A). Moreover, three constant temperature boundary conditions are implemented at the interface between the copper wires and the solid substrates (Tden = 368.15 K, Text =345.15 K, and Tann = 328.15 K). The last boundary conditions are set as a simplification to the Joule Heating module, required to describe the function of the copper wire heaters [7], in order to omit the trial and error process required to define the values of the current in each heater. This adjustment in the model is expected to result to greater temperature uniformity and higher values of DNA amplification. Furthermore, natural convection and thermal radiation boundary conditions are implemented through Equations 4 and 5 respectively [7] as presented in Figure 1A.

PCR Kinetics and Species Transport Modelling
The kinetic model presented in the work of [7] is implemented in COMSOL Multiphysics 5.4 ® . The reactions and reaction rate constants are presented in [7]. The Diluted Species Physics module of COMSOL Multiphysics 5.4 ® is used to implement the kinetics of PCR, for denaturation, extension, and annealing. Equation 6 presents the general form of the mass conservation of species in steady-state: (6) where Ck: the concentration of the k th species (k={1,2,..,7}, where 1,2,..,7 describe the doublestranded DNA molecules, two complementary singlestranded DNA molecules, two single-stranded primer molecules, and two single-stranded template primers present in the samples [7]), Rk: the reaction rate of the k th species and Dk: the diffusion coefficient of the k th species. The diffusion coefficients of Equation 6 are presented in Table 3 of [7]. A zero-flux boundary condition is set at the walls of the microfluidic channel while the inlet concentrations of the seven species are given at the input (Table 3 of [7]).

Methodology
The models of the fluid flow, heat transfer and PCR kinetics described in Section 2 are introduced in COMSOL Multiphysics 5.4 ® , using the Laminar Flow, Heat Transfer in Solids and Fluids, and Transport of Diluted Species Physics modules respectively.

Setting up the simulations
The dimensional parameters and the materials of the microchannel presented in this work are similar to the ones presented by [7]. Copper wires are placed at the bottom of the microchannel, to ensure that the three required for PCR temperature regimes are created. The width, length, and height of the microchannel vary in the simulations. The width at the extension zone is twice the width in the annealing and denaturation zones, while the length of the microchannel is adjusted according to Equation 7, to ensure that the same PCR protocol (3s:4.2s:6.2s for denaturation, annealing, and extension respectively) is followed in every simulation despite the changes in the geometrical features of the microchannel (W and H). Therefore, using the same PCR protocol enables studying exclusively the effect of the geometrical parameters of interest on the objectives.  (7) where tR,zone: (zone = {denaturation, annealing, extension}) is the residence time in each temperature regime, defined according to the PCR protocol used in the work of [7], Qvol: the volumetric flowrate, Wzone: the width of the microchannel in each temperature regime, H: the height of the microchannel and Rzone: the curvature ratio for denaturation and annealing. The volumetric flowrate at the inlet, the ambient temperature (Tamb) and the heat transfer coefficient (h) are constant for all simulations and equal to 1.8 μL/min, 5W/(m 2 ·K), and 298.15 K respectively. The gaps between the three temperature regimes and the heights of Kapton, PDMS, and PE (HKapton, HPDMS, HPE) are also constant and equal to 1.670 mm, 1.110 mm, 100 μm, 50 μm, and 50 μm respectively. The surface emissivity (ε) of Kapton, PDMS and PE is presented in Table 2 of [7]. Natural convection occurs at the walls of the unitcell, as illustrated in Figure 2 of [7]. The Joule Heating module is used for the mesh independence study and the validation with the publication of [7] and but not the Design of Experiment (DoE), as discussed below.

Mesh independence study
A mesh independence study is performed in order to ensure that the results of the simulations do not only converge but are also independent of the mesh resolution.
All five simulations of the different meshes converge successfully, and the values of the  [7,33]), the pressure drop (∆P (Pa)) and the power consumption of the heaters (Ph (W)) are recorded (the Joule Heating [7] module is used for describing the function of the copper wire heaters).
Results show that mesh independence is achieved with the selection of the ∼321,000 element-mesh, which is then used to generate the DoE [33].

Validation with [7]:
The current model is validated against the work of [7] while using the Joule Heating [7] module for the function of the copper wire heaters: • The value of [DNA] is found equal to 0.67, equal to the one presented in the work of [7].
• The power requirements of the three heaters is found to be 0.071 W for one PCR cycle, which is equal to the one presented by [7].

Optimisation Methodology
After validating the model, an optimisation problem is formed in order to optimise the performance of the microfluidic device by modifying the geometry of the microchannel. More specifically, the width (W) and the height (H) are selected as the two design variables of the optimisation problem, varying within the ranges of 100 -400 μm and 50 -80 μm respectively. These two design variables are modified (while at the same time adjusting the channel length in each temperature regime according to Equation 7), to maximise the DNA amplification (obj1) and minimise the pressure drop, ΔΡ (obj2). In order to form an optimisation problem where both objectives are minimised, obj1 is (re)defined as -2 [DNA] [DNA] . Also, in order to visualise the behaviour of the two objectives over the entire design domain, dimensionless and scaled (0-1) response surfaces are created using a polyharmonic spline. More specifically, the model described earlier (does not include the Joule Heating module) is implemented in COMSOL Multiphysics 5.4 ® for 80 DoE points, each of which corresponds to a different design case scenario. After generating the response surfaces for obj1 and obj2, genetic algorithm is used to locate the minimum values of obj1 and obj2.

Design of Experiments
The Morris Mitchel Latin Hypercube method is used to generate 80 DoE points in MATLAB ® . The code developed is based on the work of [35], and is modified to include the corner points of the design domain.

Optimisation Algorithm
The genetic algorithm, ga (MATLAB ® function [39]), is used to obtain the optimum values of obj1 and obj2. The default parameters of the ga are implemented [39].

Optimisation
The designs generating the optimum solutions for obj1 and obj2 are tested using the Joule Heating module (Tables 2 and 3 respectively). One can observe a ∼2.1 % increase in the [ ] in one PCR cycle for W = 400 μm and H = 50 μm ( Figure 1B). As far as the pressure drop is concerned, the design case of W = 400 μm and H = 80 μm ( Figure 1C) results in a ∼95.2% decrease in the pressure drop.

Conclusions
The results of this work indicate that when compared to the work of [7], there is the possibility of increasing the [ ] by ∼2.1 % in the first PCR cycle, for the design case of W = 400 μm and H = 50 μm. The design case of W = 400 μm and H = 80 μm results in a 95.2% decrease in the value of the pressure drop, compared to the one presented in the original design of [7] (W = 200 μm and H = 50μm). From a commercial perspective, varying the geometrical features of this device is expected to reduce the overall production and operational cost.
For both objective functions, the optima reside on the design space boundaries. Future work is expected to focus on the more practically relevant and interesting multi-objective problems that arise from the study of these two objectives since no single optimal design exists that minimises both obj1 and obj2. Furthermore, extra work is expected to take place on the optimisation of the gaps between the heaters, minimising their length and/or the residence time, while at the same time maintaining the improved DNA amplification efficiency and reduced pressure drop. Extra work can be performed in examining designs with greater design domains, considering that both optimal solutions fall on the boundaries of the design domain.