key: cord-0962694-pmfwss65
authors: Zagklavara, Foteini; Jimack, Peter K.; Kapur, Nikil; Querin, Osvaldo M.; Thompson, Harvey M.
title: Multi-objective optimisation of polymerase chain reaction continuous flow systems
date: 2022-03-22
journal: Biomed Microdevices
DOI: 10.1007/s10544-022-00610-6
sha: c94181a477c060075fb5279184bddaeda05dd18c
doc_id: 962694
cord_uid: pmfwss65

A surrogate-enabled multi-objective optimisation methodology for a continuous flow Polymerase Chain Reaction (CFPCR) systems is presented, which enables the effect of the applied PCR protocol and the channel width in the extension zone on four practical objectives of interest, to be explored. High fidelity, conjugate heat transfer (CHT) simulations are combined with Machine Learning to create accurate surrogate models of DNA amplification efficiency, total residence time, total substrate volume and pressure drop throughout the design space for a practical CFPCR device with sigmoid-shape microfluidic channels. A series of single objective optimisations are carried out which demonstrate that DNA concentration, pressure drop, total residence time and total substrate volume within a single unitcell can be improved by up to [Formula: see text] 5.7%, [Formula: see text] 80.5%, [Formula: see text] 17.8% and [Formula: see text] 43.2% respectively, for the practical cases considered. The methodology is then extended to a multi-objective problem, where a scientifically-rigorous procedure is needed to allow designers to strike appropriate compromises between the competing objectives. A series of multi-objective optimisation results are presented in the form of a Pareto surface, which show for example how manufacturing and operating cost reductions from device miniaturisation and reduced power consumption can be achieved with minimal impact on DNA amplification efficiency. DNA amplification has been found to be strongly related to the residence time in the extension zone, but not related to the residence times in denaturation and annealing zones.

The Polymerase Chain Reaction (PCR) has revolutionised biological research and diagnostics since its discovery by Kary Mullis in 1983 (Mullis 1990) . PCR systems perform a thermal cycling procedure to amplify DNA segments, allowing detection and identification of gene sequences using appropriate optical techniques (Does 2013) . They are now used in numerous diagnostic systems, with applications ranging from the rapid detection of infectious diseases (Park et al. 2011 ) to identification of bacteria causing micro-biologically induced corrosion in oil and gas production systems (Zhu et al. 2006; Agrawal and Lal 2009 ). An example of the former is the vital role PCR systems are playing in the public health response to COVID-19 (Abbasi-Oshaghi et al. 2020) . The PCR thermal cycling procedure consists of the three distinct stages of denaturation, annealing and extension. Denaturation takes place at ∼ 95 o C, where the double-stranded DNA denatures into pairs of single-stranded ones. The sample then enters the annealing stage at ∼ 56 o C, where the primers form primer-template complexes. The final stage, extension, generally takes place at ∼ 70 o C and is where the polymerase binds to the primertemplate complexes, catalysing the synthesis of new strands of DNA (Park and Park 2017; Schochetman et al. 1988 ). Small, discrete droplets have been used in conventional PCR devices (DBPCR) as separate chemical reactors. The droplets can provide a highly controlled and contaminantfree reaction environment with much smaller thermal mass than in CFPCR systems (Zhang and Jiang 2016) . Detailed descriptions of DBPCR systems are given in Ma et al. (2019) and Shi et al. (2020) . Despite their advantages, the comparative expense and complexity of DBPCR devices (Zhang and Jiang 2016) have motivated further development and optimisation of single-phase continuous flow PCR (CFPCR) systems, as evidenced by several studies appearing recently (Kaprou et al. 2019; Kulkarni, Goyal et al. 2021; Kulkarni, Salve et al. 2021; Hamad et al. 2021) .

Experimental and numerical studies of CFPCR systems have explored how operating and geometry variables affect the thermal cycling process and, ultimately, the efficiency of the DNA amplification process controlling the PCR yield. These have shown that the most influential parameters include the substrate's thermal conductivity, fluidic channel sizes and spacing, flow rate, while the heating arrangement has also been shown to be very important (Thomas et al. 2014; Chen et al. 2013) . Controlling the residence times in each of the thermal zones is also of key importance since insufficient dwell times can reduce PCR yield significantly; Cao et al. (2011) studied the effect of these factors on DNA amplification efficiency DNAAE), both experimentally and numerically. Combining mathematical models of the kinetics of denaturation, annealing and extension processes with models of the flow and thermal processes has proven to be highly beneficial for understanding and hence improving DNAAEs in CFPCR systems (Wang and Li 2010; Cao et al. 2011; Papadopoulos et al. 2015; Zagklavara et al. 2021) .

The effect of fluidic channel geometry on PCR performance has been studied widely, with the performance of radial (Schaerli et al. 2009 ) and spiral (Hashimoto et al. 2004 ) geometries having been benchmarked against straight channels (Chiou et al. 2001; Frey et al. 2007 ). The benefits of achieving more uniform flow and thermal conditions have also been explored. The latter was considered experimentally and numerically by Duryodhan et al. (2016) , who showed that employing diverging fluidic channels can create more uniform wall temperatures, while Gui and Ren (2006) showed that flow uniformity can be increased through the use of electro-kinetic flow. A number of studies have focused on the influence of heater arrangement, showing that it is important to control the interference and transition times between the thermal zones, and the thermal 'crosstalk' between adjacent zones, which can require larger gaps between channels and therefore hinder the drive towards device miniaturisation needed to create portable devices for diagnostic testing purposes (Kumar et al. (2013) ; Moschou et al. (2014) ; Papadopoulos et al. (2015) ; Perwez et al. (2019) . Perwez et al. (2019) have recently explored these issues in the context of using a simple, single heater CFPCR chip design. In a similar vein, the lower thermal conductivity of 3D-printable materials has been identified as a major factor limiting its application to CFPCR devices (Park and Park 2017) . Furthermore, when it comes to lab-on-chip devices, the pressure drop requirements can become very important, since they often require sophisticated and expensive microfluidic pumps (Fajrial et al. 2021; Ahn et al. 2004 ) that can be hard to integrate and fabricate (Ahn et al. 2004 ). This paper is motivated by the need to develop an effective multi-objective methodology for CFPCR devices. For example, the development of low-cost and rapid diagnostic devices for use in inaccessible regions requires effective device minaturisation and reduced power consumption, whilst maintaining the required rate of DNA amplification. The aim is to provide a powerful means of striking the appropriate balances between the conflicting performance objectives. A simulation-based optimisation methodology is developed, which uses outputs from Computational Fluid Dynamics (CFD) analyses. This approach is now commonplace in many industries, such as the aerospace and automotive ones, with the continued progress in computing power, numerical schemes and design space exploration methods, making it an increasingly powerful means of optimising complex flow systems (Khatir and Thompson 2019) . The recent review by Haftka et al. (2016) noted that the number of design variables is key. For large problems with O(1000) design variables, employing advanced adjoint methods is vital, whereas for CFPCR systems with < 100 design variables, gradient-free surrogate-assisted methods are effective. Important examples of the latter include Gaussian Process Emulators (Domingo et al. 2020) , and Moving Least Squares, which is effective at minimising the effects of numerical noise (González Niño et al. 2019) . Surrogate modelling using Machine Learning can also be effective for achieving temperature control in CFPCR systems (Lee et al. 2007; Hamad et al. 2021) .

The present study applies optimisation methods on a practical CFPCR flow problem, considered recently by Papadopoulos et al. (2015) and Zagklavara et al. (2021) . The effect of the PCR protocol on the performance of a CFPCR is investigated in detail, examining the importance of the residence time in each temperature zone (denaturation, extension, annealing). Furthermore, the design approach of doubling the channel width in the extension zone to increase the residence time there (see e.g. Papadopoulos et al. (2015) and Zagklavara et al. (2021) ) is also examined, by including it as a design variable. Four objectives of practical interest are studied: the DNA amplification efficiency (DNAAE), the total residence time, the substrate volume and the pressure drop requirements of the unitcell of a microfluidic device. Furthermore, a Pareto front is generated in order to maximise DNAAE, whilst minimising the total residence time and substrate volume of the microfluidic device (Logist et al. 2007; Hashem et al. 2017) . Apart from increasing DNAAE, reducing the total substrate volume and total residence time can lead to significant reductions in cost and processing times of the device. Furthermore, the pressure drop is also minimised in order to facilitate the development of microfluidic pumps for lab-on-chip devices, that are often highly sophisticated and expensive (Fajrial et al. 2021; Ahn et al. 2004) .

The paper is organised as follows. Section 2 describes the PCR problem of interest while Sect. 3 outlines the conjugate heat transfer problem and mathematical and numerical methods employed. Section 4 presents the results of the numerical simulations and optimisation studies, with conclusions drawn in Sect. 5.

Within a single PCR cycle, or unit-cell, the temperature of the flowing fluid through the microchannel changes as it passes through three different temperature zones -typically ∼95, ∼ 55 and ∼ 72 o C in the denaturation, annealing and extension zones respectively (Papadopoulos et al. 2015) . A unit-cell corresponds to one of the N PCR cycles that are placed consecutively in a serpentine arrangement, as presented in Fig. 1 . The temperature changes along the length of the microchannel are designed to increase the DNA concentration significantly by the time the fluid exits the channel. The cases considered here are based on the chip substrate materials, Kapton, PDMS and PE (Table 1) , and the design parameters (Table 2 ) used in Papadopoulos et al. (2015) and Zagklavara et al. (2021) .

The flow is steady and is governed by the incompressible Navier-Stokes equations (Eqs. 1 and 2).

where is the fluid density, the velocity vector, p: the pressure, the viscosity and the external forces applied to the fluid, such as buoyancy force due to gravitational acceleration, Lorentz forces etc (McDonough 2009; Gerbeau and Le Bris 2000) . Flow is laminar since an indicative value of Reynolds number, Re ∼ 0.33 (Eq. 3) can be calculated for Q vol = 3 ⋅ 10 −11 m 3 ∕s , H Fluid = 50 m, W 2 = 400 m, t R,den = (1) where Equations 1 and 2 are solved on the geometry appearing in Fig. 2 , subject to the following boundary conditions: (i) no-slip at the microchannel walls; (ii) fully-developed flow and a value of average inlet velocity, U in , at the inlet of the serpentine channel; (iii) zero (relative) pressure at the exit of the microchannel.

Steady state, conjugate heat transfer is modelled via Eq. 6:

where u=0 everywhere except in the fluid domain. Q heater,j is the heat generation rate of the j th (j = {1, 2, 3}) heater, and is only non-zero at the j th heater-kapton interface. A different heat generation rate is required at each heater to achieve the desired set points of 95, 55 and 72 o C in the denaturation, annealing and extension zones respectively. Q rad,i is the heat flux due to thermal radiation (Eq. 7 (Stefan-Boltzmann law)) of the i th solid substrate (i = {Copper, PDMS, PE, Kapton}), and is only non zero at the outer surfaces of the substrate materials. Q nat.conv is the heat flux due to the heat losses to the ambient, and is given by Eq. 8:

where T amb : the ambient temperature, i : surface emissivity for solid i, : the Stefan -Boltzmann constant and h: heat transfer coefficient.

The boundary conditions are applied on the geometry appearing in Fig. 2 as follows: (i) a periodic boundary condition on temperature at the inlet and outlets of the channel; (ii) the heater temperatures at the copper-solid interface in the denaturation, extension and annealing zones are set to T den = 95 o C, T ext =72 o C and T ann = 55 o C, respectively; (iii) 

Several kinetic models have been developed for the reactions in PCR systems -see for example those of Hunicke-Smith (1998); Athavale et al. (2001) ; Aach and Church (2004) ; ; Papadopoulos et al. (2015) and Chen and Li (2018) . The general equations for the steady state mass conservation of the species are given by Eqs. 9 and 10. The five reactions and the reaction rate constants ( k + A , k − A , k + D , k − D , k E ) considered in this work are presented in Appendix 2, and are described in detail by Papadopoulos et al. (2015) .

where C k is the concentration of the k th species (k = {1,2,..,7} corresponding to S 1 S 2 , S 1 , S 2 , P 1 , P 2 , S 1 P 2 and P 1 S 2 respectively (see Appendix 2)), R k is the reaction rate of the k th species and D k is the diffusion coefficient of the k th species. The reaction rates are presented in Eqs. 24-30 in Appendix 2, while the diffusion coefficients of the species in the set of Eq. 10 are presented in Table 3 (Papadopoulos et al. 2015) .

The implemented boundary conditions are: (i) no flux at the sides of the microchannel, excluding the inlet and outlet; (ii) initial species concentrations are given in Table 4 ; (iii) zero inward species flux at the exit of the microchannel ( ⋅ D k ∇C k ).

The design of the microchannel is based on the design offering the maximum DNA amplification, presented in the publication of Zagklavara et al. (2021) . The coupled series of flow, heat transfer and species transport equations described above are solved subject to the boundary conditions using COMSOL Multiphysics 5.4 (COMSOL 2021), as part of the optimisation study. The material properties and the dimensions of the design parameters of the serpentine channel are presented in Tables 1 and 2 respectively. The properties of the fluid are those of water, while the PCR kinetics are described in detail in Appendix 2. The values of the volumetric flowrate at the inlet ( Q vol ), the ambient temperature ( T amb ), the heat transfer coefficient (h), the gaps between the three temperature regimes ( L 2 , L 4 ) and the heights of Kapton, PDMS and PE ( H Kapton , H PDMS , H PE ) are equal to those used by Papadopoulos et al. (2015) (Table 2) . Natural heat convection occurs at the walls of the channel, as illustrated in Fig. 2 . The ambient temperature and convective heat transfer coefficient are set to T amb = 25 o C and h = 5W∕(m 2 ⋅ K) respectively. As far as the surface-to-ambient radiation is concerned, the surface emissivity of all materials is presented in Table 1 . The channel lengths obtained by Eq. 11 include the 180 o circular arc of R zone = 500 m (when applicable). The effect of the microchannel design variables, W 2 and H Fluid on the PCR amplification efficiency and the pressure drop was studied by Zagklavara et al. (2021) 50 m respectively, in order to further study designs that offer improved DNAAE. The values of the residence times ( t R,den , t R,ext and t R,ann ) and the three channel lengths in the denaturation, extension and annealing zones vary in each simulation, in order to observe the effect that the PCR protocol has on the objectives of interest. More specifically, L 4 , L 8 , and L 6 are calculated by Eq. 11 for the different values of t R,den , t R,ext and t R,ann . Furthermore, W 3 is selected as the fourth variable and is defined according to Eq. 12, where z w3 is a parameter ∈ [0, 1] . The selection of the fourth variable, W 3 , is made in order to study the benefit Table 2   Table 3 Parameter Values (Papadopoulos et al. (2015) Parameter

Values [ m 2 ∕s] Description 

of doubling the width of the microfluidic channel in the extension zone (as originally used by Papadopoulos et al. (2015) ).

The effect of mesh density is considered for the case with W 2 = 200 m and H Fluid = 50 m , with five different mesh densities with 163,517, 321,151, 865,781, 4,035,872 and 6,133,359 elements. The Joule Heating model is implemented to describe the function of the copper wire heaters (Appendix 3), as performed by Papadopoulos et al. (2015) . The effect of mesh density on DNA amplification ( log 2 of the ratio of the average concentration of double stranded DNA at the end of the first cycle to the initial one), pressure drop ( ΔP(Pa)) and power consumption of the heaters ( P h (W)) is given in Table 5 . Table 6 presents the values of the residual errors for the temperature (T), [DNA] and velocity (U) together with the computation times for the five meshes. This shows that the solutions on the mesh with 321,151 elements are effectively mesh independent and all results presented below are obtained on this mesh.· Comparisons with the results obtained here compared to those of Papadopoulos et al. (2015) show that: (i) log 2 of the ratio of the average concentration of double stranded DNA at the end of the first cycle to the initial concentration predicted

is the same value, namely 0.67; (ii) the power requirements of the unit cell for performing 1 PCR cycle is also identical, 0.071 W. Note that the indicative power consumption for the denaturation copper wire heater with a rectangular crosssection is calculated as presented in Eqs. 31-36 in Appendix 3. The final comparison is with the temperature uniformity (T.U.) in each temperature regime while varying the inlet velocity, is shown in Fig. 3 . The agreement is once again generally very good with the results presented in Papadopoulos et al. (2015) . The temperature uniformity values are calculated via Eq. 13.

This optimisation problem focuses on further improving the optimum design of [ 400, 50] presented by Zagklavara et al. (2021) , by studying the effect that the implemented PCR protocol (residence times) and one additional geometrical parameter ( W 3 ) can have on performance objectives of interest in a unitcell. Each unitcell is identical to the next, apart from the species concentrations -hence the periodic boundary conditions mentioned in Sect. 2.2. As a result, improving the performance of 1 unitcell leads to the improvement of all the cycles and hence the entire device.

(13) More specifically, the effect of residence times in the denaturation ( t R,den ), extension ( t R,ext ) and annealing ( t R,ann ) zones together with the channel width in the extension zone ( W 3 ) is investigated on the DNA amplification, pressure drop, total residence time and total substrate volume. The channel lengths are adjusted via Equation 11 to achieve the values of t R,den , t R,ext and t R,ann . A surrogate-enabled approach is adopted and design of experiments is performed. The ranges of the residence times are created by t R,zone | (Papadopoulosetal.,2015) ± 1.5s (Table 7 ). The range of the fourth design variable, W 3 , is set at 400-800 m in order to examine the benefits of increasing (up to twice the width of the microchannel in the other zones, W 2 ) the width in the extension zone in particular, as performed by Papadopoulos et al. (2015) . The material properties and the dimensions of the design parameters of the channel are presented in Tables 1 and 2 respectively, and are based on the design proposed by Papadopoulos et al. (2015) .

The objective functions considered are obtained from the dimensionless measurement of the DNA amplification ( log 2 ( [DNA]

is the average DNA concentration at the end of the channel and [DNA] o the initial DNA concentration), the unitcell pressure drop along the microchannel ( Δ P (Pa)), the total unitcell residence time ( t R,total (s)) and total unitcell substrate volume ( V s,total ( m 3 )). More specifically, COMSOL Multiphysics is used to obtain the values of the four objectives, which are then non-dimensionalised ( obj 1 , obj 2 , obj 3 and obj 4 for −log 2

and V s,total respectively) (scaled to lie between 0-1) in the generated metamodels. obj 1 is defined as the negative of −log 2

[DNA]

[DNA] o in order to switch to four minimisation studies. The Morris Mitchel Latin Hypercube method is used to generate 160 sampling points, using code based on the work of Julie (2012), after modifying it to include the sixteen corner points of the design domain. The 160 sampling points are presented in Appendix 1. The computational model is then evaluated at the 160 sampling points and metamodels for the four objective functions are created using Neural Networks (NN).

As far as the metamodels are concerned, feed-forward NNs (Leijnen and Veen 2020) and Levenberg-Marquardt backpropagation are used for data fitting, based on the matlab function fitnet (MathWorks 2020b). The Mean Squared Error (MSE) performing function is selected together with k-fold evaluation (Manriquez-Sandoval 2021), to test and improve the quality of the NNs. The k-fold method is often used for the evaluation of the performance of classification algorithms, especially for larger datasets (Wong 2015) . Such an example is the work of Abellán-García (2021), that used the k-fold validation method to train an artificial neural network with one hidden layer. The effects of the number of hidden layers, together with the % of testing and training data are investigated for each objective function. NNs with small numbers of hidden layers ([2], [4]) are found to be unable to describe the behaviour of the system adequately, leading to the failure of the optimisation algorithms in obtaining the optimal design solutions. The [4, 4] setup for the hidden layers is selected, since it offers < O(10 −5 ) accuracy in the prediction of the optimum designs. Table 8 presents the designs offering low values of MSE ( < O(10 −5 ) ).

The e05jbc function of the NAG optimisation library (NAG 2020), which is based on the Multi-level Coordinate Search method described in Huyer and Neumaier (1999) , uses the meta-models of the objectives to solve the optimisation problems. Subsection 4.1 presents the metamodels for obj 1 , obj 2 , obj 3 and obj 4 . Subsection 4.2 then describes the optimisation method used to locate the optimum values for the four objectives, which is based on the e05jbc NAG routine (NAG (2020)). Furthermore, the results of a multi-objective optimisation study are also presented in the form of a Pareto front, showing the available compromises between competing objectives.

The fitnet matlab function is used to generate the NNs for −log 2 ( [DNA]

[DNA] o ) , Δp , t R,tot and V S,tot . The values of the objectives are scaled appropriately between 0-1 (see Appendix 6). The 3D response surfaces (for constant values of z w3 ) are developed using the libraries presented by Zhivomirov (2021) , and are presented in Figs. 4, 5, 6 and 7 for log 2 ( [DNA]

[DNA] o ) , Δp , t R,tot and V S,tot respectively. The colorbar is used to present the values of the objectives, while a 3D response surface is printed for a different value of the fourth design variable, z w3 . The sampling data points used to create [DNA] o ) (-) data (colorbar) for (a) z w3 = 0 and (b) z w3 = 1. The optimum solution is presented in a black square in Fig. 4a the response surfaces are provided in Table 14 of Appendix 1. Appendix 5 presents the response surfaces for all four objectives for more values of z w3 . The correlation coefficients between the DNA amplification -total residence time and the DNA amplification -individual residence times are given in Table 9 . 

As part of the single-objective studies, the minima of the metamodels of obj 1 , obj 2 , obj 3 and obj 4 are found at [ t R,den (s), Design 4 (see Table 10 ) offers a 16.42% increase in the value of log 2 Table 10 Details of the designs appearing in Fig. 8 and Table 11 The design variables in each study are presented in bold CW Current Work, Z (Zagklavara et al. (2021) Papadopoulos et al. (2015) . In order to examine the significance of the ∼ 5.7% increase in [DNA] for a single unitcell, ten consecutive PCR cycles are simulated for Designs 2 and 4 (see Table 16 , Appendix 4), using the Joule Heating model for the function of the copper wire heaters. The results are presented in Fig. 8 . According to the data obtained for Design 4, this ∼ 5.7% increase in [DNA] in the first PCR cycle, is expected to increase the concentration of DNA approximately by ∼ 32% in ten cycles (compared to Design 2). Furthermore, by offering a ∼ 51% reduction in the pressure drop requirements, the operating cost of such device is expected to be reduced significantly. However, this design also leads to an increase of 58.7% in the total residence time.

On the other hand, Design 6 (see Table 10 ) leads to a 80.50%, 17.80% and 43.23% decrease in the values of pressure drop, total residence time and total substrate volume respectively. However, this design also comes with a 13.43 % decrease in log 2

[DNA] o (or a 6.6% decrease in the [DNA]) compared to the one presented by Papadopoulos et al. (2015) . Figure 9 presents a comparison between the different unitcell designs, together with their temperature profiles. Figure 10 shows the DNA concentration profiles at the middle plane in the fluid domain for the two designs optimising the four objectives.

The single-objective optimisation results show that conflicts between the objectives results in a complex multi-objective design problem. For the purposes of visualisation, three out of the four objectives ( log 2

, t R,tot , V S,tot ) are selected within a multi-objective optimisation to generate a Pareto front (Fig. 11) . The Pareto front is hence a 3D plot, that is developed using the gamultiobj function (MathWorks 2020a), in order to demonstrate the available compromises between the three objectives. The values of obj 1 , obj 3 and obj 4 are dimensionless and scaled between 0 and 1, to aid visualisation of the multi-objective results. The values of FunctionTolerance and MaxGenerations are adjusted to 1 ⋅ 10 −6 and N DVARS ⋅ 200 , where N DVARS is the number of design variables ( N DVARS = 4). Three of the optimal solutions in the Pareto front plot are validated using the simulation model (Tables 12 and 13 ), deviating less than ∼ 0.15% for all three cases. The Pareto front offers the ability to significantly ameliorate the performance of the device depending on the requirements of the designer/engineer. For example, the design of Point 2 appearing in Tables 12 and 13, illustrates the ability to improve t R,tot and V S,tot by 24.64% and 25.75% respectively when compromising on log 2 for 10 PCR cycles. The details of the four designs are presented in Table 10 . Designs 1 and 2 present the designs by Papadopoulos et al. (2015) and its validation (current work). Designs 3 (Zagklavara et al. (2021) and 4 (current work) present the designs offering maximum DNA amplification 

The development of practical CFPCR devices offers a complex, multi-objective design challenge due to the conflicts between the required DNA amplification and other practical constraints, such as manufacturing and operating costs related to size and power consumption. The latter are particularly important for low-cost devices targeted at lowerincome countries. This paper has developed an effective multi-objective optimisation methodology which allows designers to strike an appropriate balance between the various competing objectives. The methodology uses a series of high fidelity CHT simulations which also account for the kinetics of the DNA amplification to predict the DNA amplification efficiency. As a basis for the chip device, the width and height of the microchannel are constant and (along with parameters of volumetric flowrate, gap lengths and material properties) are consistent with the work of Zagklavara et al. (2021) .

Results indicate that doubling the width of the microchannel in the extension zone, together with the residence time in denaturation and annealing zones does have significant effect on the DNA amplification. The residence time in extension zone however has been found to be strongly related to the Table 12 Validation of three points appearing in the Paretofront plot (Fig. 11 ) Tables 12 and 13 ). The black dots and red triangles represent high fidelity data obtained using COMSOL Multiphysics DNA amplification. From consideration of the Pareto front, several designs are presented, and depending on design priorities, different design solutions can be used to improve the designs of Papadopoulos et al. (2015) and Zagklavara et al. (2021) . The Pareto front includes designs ranging from those with low DNA amplification, low total device volume and operation time to high values of the DNA amplification, high total device volume and operation time or design compromises between the three objectives. The first types of design offers the ability to reduce the total material volume, operation time and pressure drop requirements by up to ∼ 43.2% , ∼ 17.8% and ∼ 80.5% respectively. However, such design modifications can lead up to ∼ 6.6% reduction in the [DNA] in the a unitcell. Single objective optimisation on the DNAAE, shows that it is possible to increase DNA concentration by up to ∼ 5.7% in the first PCR cycle, which simulations show results in an increase of ∼ 32% over ten PCR cycles. At the same time, this design offers a reduction in the total pressure drop ( ∼ 50.5% ) together with a small reduction in the material volume ( ∼ 5.6% ), having however a ∼ 58.7% increase in the total operating time. According to the results obtained, all designs have the potential to minimise pumping requirements for such devices; with reductions in pressure drop allowing for smaller pumps to be used (particularly when building integrated lab-on-chip devices). The smaller size and reduced pumping requirements also minimise power requirements, which is an important consideration when these are used within handheld devices containing their own power-sources. This supports the ongoing efforts to develop field-ready microfluidic systems.

Future research directions include comprehensive experimental validation of the optimisation results, and their extension to a wider range of design variables.

See 

In the denaturation zone, the double-stranded DNA molecules, S 1 S 2 , dissociate into two single strands, S 1 and S 2 (Reaction 14):

where k D is the denaturation constant for melting of D at melting temperature; k −D is the denaturation constant for binding of S at melting temperature; k E is the enzyme inactivation constant. The arrow symbols " ← " and " → " are used to denote net forward and backward reactions. The high

melting temperature causes thermal denaturation of the enzyme responsible for the DNA amplification.

In annealing zone, the single-stranded primer molecules, P 1 and P 2 , bind to S 2 and S 1 respectively, and form the singlestranded template-primer complexes, P 1 S 2 and S 1 P 2 (Reactions 15 and 16):

where k + A is the annealing coefficient of P 1 and P 2 to S 2 and S 1 respectively,; k − A is the dissociation coefficient of P 1 S 2 and S 1 P 2 . In extension zone, the polymerase enzyme binds to P 1 S 2 and S 1 P 2 to form the single-stranded template-primer-enzyme complexes. Then these complexes dissociate into the enzyme and the DNA molecules at the beginning of the subsequent denaturation step (Reactions 17 and 18):

where k E is the addition constant of the enzyme to P 1 S 2 and S 1 P 2 .

The temperature dependence of the various rate constants mentioned earlier, k + D , k − D , k E , k + A and k − A , are given by Eqs. 19-23, as demonstrated also in the work of Papadopoulos et al. (2015) :

where T is the temperature in K, and the k + o , k − o , k + 1 , k − 1 and k 2 constants are presented in Table 15 . The reaction rates are given by Eqs. 24-30.

Appendix 3

Using the Joule Heating model, the power consumption for the heater in the denaturation regime for the design case presented in the work of Papadopoulos et al. (2015) is calculated as follows:

where : the coefficient of thermal expansion of copper (0.0386 K −1 ), 293K : the reference resistivity of copper at 293 K (1.68 ⋅10 −8 Ωm ), (T) den : the resistivity of the ) and length ( L den,heater ) of the copper wires used in the simulations is 10 −4 m, 2 ⋅ 10 −5 m and 0.0586 m respectively, while the value of I den is found to be equal to 0.2345 A by trial and error. The copper wire is bent nine times, resulting to eight straight parts ( L den,straight =0.0071 m) covering the bottom of the heater in a serpentine shape (see Fig. 2 ). The power consumption of the heaters of extension and annealing zones are calculated in the same way, and are equal to 0.02367 and 0.01201 W respectively.

See Tables 16 and 17 Appendix 5 Table 16 Concentrations of all PCR products for 10 cycles for Design 4 (see Table 10 )

Design 4 (see Table 10 )

Start 5.7100 ⋅ 10 −9 0.0000 0.0000 3.0000 ⋅ 10 −4 3.0000 ⋅ 10 −4 0.0000 0.0000

Cycle 1 9.7358 ⋅ 10 −9 5.0622 ⋅ 10 −10 5.0622 ⋅ 10 −10 2.9999 ⋅ 10 −4 2.9999 ⋅ 10 −4 6.7201 ⋅ 10 −10 6.7201 ⋅ 10 −10 Table 17 Concentrations of all PCR products for 10 cycles for Design 2 (see Table 10 )

Design 2 (see Table 10 ) 

All values of the objective functions are scaled between 0−1 for the purposes of multi-objective optimisation. Furthermore, just for the first objective, the values of −log 2

[DNA]

[DNA] o are scaled between 0−1, in order for all optimisation studies to be minimisation problems. Equation 37 is used to scale all values, while Table 18 presents two cases of scaling values of objectives.

Author contributions All authors contributed to the study conception and design. Data collection and analysis were performed by Foteini Zagklavara. The first draft of the manuscript was written by Foteini Zagklavara and Harvey M. Thompson and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Funding The funding is received by EPSRC, as part of the CDT in Fluid Dynamics at the University of Leeds (EPSRC (2021).

The datasets generated during and/or analysed during the current study are included in this published article and its supplementary information files.

Code availability Any code or file developed during the current study are available from the corresponding author on reasonable request.

Ethics approval Not applicable.

Consent for publication All authors gave consent for the publication of the final manuscript.

All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.

(37) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. 

Mathematical models of diffusion-constrained polymerase chain reactions: basis of high-throughput nucleic acid assays and simple self-organizing systems

Diagnosis and treatment of coronavirus disease 2019 (COVID-19): Laboratory, PCR, and chest CT imaging findings

K-fold validation neural network approach for predicting the one-day compressive strength of UHPC

Rapid detection and quantification of bisulfite reductase genes in oil field samples using real-time pcr

Disposable smart lab on a chip for point-of-care clinical diagnostics

Modeling & Simulation Microsystems

Plastic microfluidic chip for continuous-flow polymerase chain reaction: Simulations and experiments

Analytical study of a microfludic DNA amplification chip using water cooling effect

Analysis of PCR kinetics inside a microfluidic DNA amplification system

A closed-cycle capillary polymerase chain reaction machine

Tchobanoglous, MWH's water treatment: principles and design

Simulate real-world designs, devices, and processes with multiphysics software from COMSOL

Polymerase chain reaction

Using ice cores and Gaussian process emulation to recover changes in the Greenland ice sheet during the last interglacial

A simple and novel way of maintaining constant wall temperature in microdevices

EPSRC centre for doctoral training in fluid dynamics at leeds

A frugal microfluidic pump

Autonomous microfluidic multi-channel chip for real-time pcr with integrated liquid handling

A basic remark on some navier-stokes equations with body forces

Computational fluid dynamic enabled design optimisation of miniaturised continuous oscillatory baffled reactors in chemical processing

Numeric simulation of heat transfer and electrokinetic flow in an electroosmosis-based continuous flow PCR chip

Computational fluid dynamics analysis and optimisation of polymerase chain reaction thermal flow systems

Parallel surrogate-assisted global optimization with expensive functions-a survey

Multiobjective optimization of a plug flow reactor using a divide and conquer approach

PCR and cycle sequencing reactions: A new device and engineering model

Global optimization by multilevel coordinate search

Surrogate model optimization toolbox, matlab central file exchange

Ultrafast, low-power, PCB manufacturable, continuous-flow microdevice for DNA amplification

CFD-enabled design optimisation of industrial flows-theory and practice

Miniaturized and IoT enabled continuous-flow based microfluidic PCR device for DNA amplification

Miniaturized thermal monitoring module with CO 2 laser ablated microfluidic device for electrochemically validated DNA amplification

Thermal modeling and design analysis of a continuous flow microfluidic chip

Neuro-genetic optimization of temperature control for a continuous flow polymerase chain reaction microdevice

Multidisciplinary Digital Publishing Institute Proceedings

Control Conference (ECC)

Peanut detection using droplet microfluidic polymerase chain reaction device

K-fold cross validation with & without random shuffle data

Find pareto front of multiple fitness functions using genetic algorithm

Function fitting neural network

Lectures in elementary fluid dynamics: physics, mathematics and applications

The unusual origin of the polymerase chain reaction

All-plastic, low-power, disposable, continuous-flow PCR chip with integrated microheaters for rapid DNA amplification

Advances in microfluidic PCR for point-of-care infectious disease diagnostics

Thermal cycling characteristics of a 3D-printed serpentine microchannel for DNA amplification by polymerase chain reaction

Comparison of continuous-flow and static-chamber PCR devices through a computational study: the potential of flexible polymeric substrates

Numerical investigation of design and operating parameters of thermal gradient continuous-flow PCR microreactor using one heater

Pipe flow: A practical and comprehensive guide

Continuous-flow polymerase chain reaction of single-copy DNA in microfluidic microdroplets

Polymerase chain reaction

Thermal gradient continuous-flow PCR: a guide to design

Optimization of a microfluidic microarray device for the fast discrimination of fungal pathogenic DNA

Multi-physics simulational analysis of a novel PCR micro-device

Kinetic characteristics of continuous flow polymerase chain reaction chip: A numerical investigation

Performance evaluation of classification algorithms by k-fold and leave-one-out cross validation

Optimisation of microfluidic polymerase chain reaction devices. E3S Web Conf

A review on continuous-flow microfluidic pcr in droplets: Advances, challenges and future

5D data visualization with matlab

Rapid detection and quantification of microbes related to microbiologically influenced corrosion using quantitative polymerase chain reaction