Abstract
This paper details the infinite dimensional dynamics of a prototype microfluidic thermal process that is used for genetic analysis purposes. Highly effective infinite dimensional dynamics, in addition to collocated sensor and actuator architecture, require the development of a precise control framework to meet the very tight performance requirements of this system, which are not fully attainable through conventional lumped modeling and controller design approaches. The general partial differential equations describing the dynamics of the system are separated into steady-state and transient parts which are derived for a carefully chosen three-dimensional axisymmetric model. These equations are solved analytically, and the results are verified using an experimentally verified precise finite element method (FEM) model. The final combined result is a framework for designing a precise tracking controller applicable to the selected lab-on-a-chip device.
1. Introduction
Biomedical microdevices usually require the very tight performance specification to be satisfied so that they become properly functional. Carefully designed controllers can guarantee certain performance specifications such as robustness and stability. Successful controller design process requires a precise mathematical model. A perfect mathematical model is the one that captures the highest amount of the dynamics of the system. Likewise, numerical tools can provide precise insights into a system’s dynamics. Finite element method (FEM) is a common numerical tool that can provide precise simulation of the steady-state and transient behaviour of a microdevice. However, FEM models are not suitable for controller design because they are not in a compact mathematical form. In this work, we present a mathematical model that captures infinite dimensional dynamics of the heat distribution inside a prototype system. The selected system is a lab-on-a-chip (LOC) device that is designed and built to perform the Polymerase Chain Reaction (PCR) and capillary electrophoresis (CE) processes. The three-dimensional complete FEM model of the system is used as the reference to derive the mathematical model. The FEM model is experimentally verified with the actual setup.
Lab-on-a-chip (LOC) devices are finding growing attention in various applications ranging from chemical analysis and water contamination sensors to genetic analysis [1]. A LOC device can be thought of as a miniaturized instrument capable of performing single or multiple laboratory processes that integrate various functionalities such as sample pretreatment, sample transportation, genetic amplification, mixing, reaction, separation, and detection. LOC devices offer advantages beyond conventional biochemical and clinical laboratory settings. Dramatic reduction of sample volume and reduced consumption of costly reagents are two distinct benefits of miniaturized laboratory processes. Most significantly, miniaturization improves analysis characteristics, such as shorter analysis time and higher separation performance, and can even enable innovative applications that are not otherwise attainable.
Genetic and molecular analyses are among the major applications of LOC devices and have become feasible by the adoption of the Polymerase Chain Reaction (PCR) technique. Using PCR, specific regions of even a single DNA molecule can be replicated up to about 107 to 108 times. This process is also called selective amplification of DNA. The amplified DNA segments can be later analyzed and compared to a standard pattern to find out if they correspond to some disease condition. PCR has been traditionally a lengthy molecular biology operation, but the LOC concept, proposed in the early 1990s [2, 3], has pushed this technology to attain ever shorter reaction times. DNA amplification using PCR involves repeatedly cycling the temperature of the reaction chamber in three predefined temperatures. Optimum reaction temperature and exposure times of the sample to each temperature are key factors for successful DNA amplification [4]. The PCR-LOC research community emphasizes the importance of ±1°C temperature precision and rapid transition time for high amplification efficiency. At the same time, overshoots and undershoots, which are likely to happen at higher speeds, have a detrimental impact on the amplification yield. For example, overshoots at denaturation temperature will reduce the lifetime of replicating enzyme (Taq polymerase). Undershoots at annealing temperature may cause the amplification of unspecified products. Nevertheless, rapid thermal cycling of reaction chamber in DNA amplification not only reduces the amplification time, but it is also important to ensure the stability of biological reagents during amplification as the Taq polymerase has a half lifetime that decreases rapidly at temperatures close to 100°C [5]. It is important to note that about half of the time in a standard PCR process is spent in transitions which are undesirable. As a result, researchers have focused their attention to reduce transition time between temperature steps in PCR process. To achieve this goal, both structural design of the device and thermal control of the process need to be optimized. Consequently, various platforms for PCR-LOC devices have been reported in the literature [6] which can be roughly divided into two major categories: continuous-flow PCR and static PCR [7]. Furthermore, different heating techniques have been chosen, where thermoelectric modules (TEMs) and thin-film resistive heaters [8–11] are the most commonly used techniques. Other techniques, such as infrared radiation [12] and microwaves [13, 14], have also been reported. These different structural designs are discussed in detail in the excellent review papers by Zhang et al. [7], Roper and coworkers [6], Verpoorte [15], and DeMello [16]. The latest achievements in PCR-LOC devices design can be found in the recently published survey paper by Zhang and Ozdemir [1].
Most of the published work on PCR-LOC devices has been devoted to the structural design aspects, and less attention has been paid to optimize the thermal control of the process. As a result, there is very little work in the literature on the dynamics and control of these microsystems. Although neural network based controllers [17] and hybrid controllers [18] are used for thermal control of PCR-LOC devices, PI/PID controllers are still the most commonly preferred controllers to deal with the thermal control problem of PCR-LOCs [5, 19, 20]. What makes the thermal control of most of the PCR-LOC devices a challenging task is the fact that the temperature inside the PCR-LOC chips not only depends on time but is also spatially distributed. Thus lumped modeling and controller design cannot exactly address the problem, fundamentally because in this approach the distributed parameter nature of the system is neglected. We observed this critical problem to some extent in the formerly designed PCR-LOC chip [21] where switching controllers were employed [8]. The PCR-LOC system studied in this paper is the next generation of the one described in our previous work [22]. The new design has undergone significant improvement in structural design with respect to its predecessor.
While lumped modeling can precisely describe the dynamics of lumped systems, it usually cannot provide a fair description of the dynamics of the distributed parameter systems. Those systems that can be described by partial differential equations are often referred to as distributed parameter systems or PDE systems which feature distributed dynamics and are categorized under a general set of the systems known as complex dynamical systems [23]. Distributed parameter systems are subsets of a more general set called infinite dimensional systems because they can be mathematically formulated as differential equations on an abstract linear vector space of infinite dimension [24]. In contrast, lumped systems can be formulated mathematically as ordinary differential equations on a finite dimensional linear vector space. Designing controllers for distributed parameter systems is called infinite dimensional controller design or PDE controller design [24–26].
This paper characterizes the heat distribution dynamics of a microfluidic thermal process. The work begins with a description of the process and experimental setup. Then results of FEM simulations on the system were analyzed to choose a suitable structure for the model. Next the PDE equations describing the model in steady-state and transient parts were derived. These equations were solved analytically from which the dynamic temperature distribution in the whole domain can be expressed by an infinite dimensional model. Next, we optimize the calculated model to minimize the error due to the simplification of the actual multidomain system to an abstract single-domain system. FEM simulation and laboratory experiments are used to verify the model. The result of this work is a framework for boundary control design for thermal control of the aforementioned PCR-LOC device.
2. Materials and Methods
2.1. Chip Configuration
Genetic analysis usually consists of three steps: sample preparation, amplification, and detection. In the first step the DNA is extracted from a raw sample such as blood. Then, a selected fragment of this DNA is amplified through the PCR process to obtain sufficient DNA copies to perform detection. Next, the amplified region of DNA is labelled with a fluorescent dye, and finally it is detected through various analysis techniques [27].
The integrated genetic analysis platform that has been developed in the Applied Miniaturization Laboratory (AML) at the University of Alberta is capable of performing complete genetic analysis on a single microfluidic chip [28]. This system is depicted in Figure 1. The microchip is composed of two layers of Borofloat glass with a polydimethylsiloxane (PDMS) layer in the middle. Reaction chamber and fluid channels are carved on one of the glass layers, and a Pt heater/sensor is patterned on the bottom glass plate using standard microfabrication techniques. A close-up view of the microfluidic chip, depicted in Figure 2, clearly shows the configuration of the platinum heater-sensor ring. A schematic representation of the components within the chip is shown in Figure 2(a). Perspective view of the glass/PDMS/glass chip with horizontal laser beam is shown in Figure 2(b) and the chip atop a heatsink and light detector is depicted in Figure 2(c). To improve the cooling rate during thermal cycling, the microchip is mounted on a copper heat sink. The heat sink has a circular opening aligned to the center of the PCR chamber and heater/sensor, ring axis. The radius of this opening is optimized to increase the cooling rate without having negative effect on the heating rate [21]. We are mostly interested in heat transfer dynamics around the reaction chamber and heater-sensor ring. The chamber, heater/sensor and heat-sink opening are circular and share the same axis. This suggests that we should look for an axisymmetric cylindrical structure for our model. The half cross-section view of the chip along a vertical plane parallel to the smaller side and passing through the center of the chip is shown in Figure 3. The irregular and multilayer multimaterial structure of the PCR-LOC system is not suitable to form a set of solvable analytic equations; additional simplifications are needed. To find a valid simplified model describing the dynamics of the chip where the heater and chamber are located, we first have to analyze the temperature profile around the heater and chamber area. Figure 4 depicts the temperature profile along the chip thickness for some selected distances from the center of the chip, ranging from 2.5 mm to 5.0 mm. Referring to Figure 4, the temperature profile along the thickness of the chip at the upper glass layer (1.354 mm to 2.454 mm from the bottom of the chip) shows a semilinear trend. By choosing the radius of 3.5 mm from the center of the chip, the slope is negligible, as the temperature varies for less than 0.8°C, which is even less than the acceptable error for the chamber temperature. Therefore, with reasonable approximation, we can consider a constant temperature wall as the boundary, where this boundary is positioned at a distance of 3.5 mm from center of the chip. These calculations, in addition to several other FEM simulations, have been performed using COMSOL Multiphysics 3.5a on a comprehensive FEM model of the chip that describes its layered geometry, material properties, and boundary conditions with a high degree of accuracy. More FEM simulations disclosed some other important facts about this chip which can be summarized as follows.(1)The amount of heat going to the upper portion of the chip through the PDMS and top glass layer has an approximately constant ratio to the total heat being dissipated by the heater.(2)If we consider a disk-shape bounded volume at the center of the chip, the amount of heat lost by natural convection in air through the upper surface of the disk is negligible compared to the amount of heat lost by conduction through the wall.(3)If the PDMS layer was replaced by a glass layer, only a small percentage of the dynamic temperature distribution would be affected, due to the alteration of thermal properties (thermal conductivity and specific heat capacity). This small effect on the dynamics can be compensated later by using a correction factor.(4)The vertical temperature profile in the top glass layer at a distance of 3.5 mm from the center is approximately constant.
(a)
(b)
(c)
2.2. Model Structure
Using these observations we can conclude that, by picking a single layer glass disk, we can model the temperature distribution at a closed disk in the upper side of the chip with a good level of approximation relative to the temperature accuracy of ±1°C demanded by PCR-LOC systems. This model structure is shown in Figure 5 where the heater is defined by the circular strip at the bottom surface with width and distance to the vertical axis. It is worth noting that we only modelled the area of interest in the chip, and we have put aside the rest of the chip. Furthermore, we defined a comparison aperture that can be characterized by a fixed size rectangular area in cross-section view of the actual PCR-LOC microchip and our model. The comparison aperture contains the reaction chamber, and it is discussed in detail in the simulation results section. We chose this scenario to arrive at a model structure which is simple enough to find its analytical solution meanwhile it is accurate enough to characterize the actual thermal process. A simple estimate of the thermal capacity and radial thermal resistance of the PDMS layer shows that they are small compared to those of the glass layers. It is expected that treating the PDMS layer as if it has the same thermal properties of glass will make the systems dynamic behavior 9% too fast (because the radial resistance is predicted to be 10% too low and the thermal capacity is estimated to be about 1% too high). It is possible to correct for this by scaling the thermal parameters of the glass. The assumptions made will result in underestimating the vertical gradient, but that gradient cannot be controlled by control mechanism and must instead be dealt with at the design and simulation level (i.e., dynamic simulations of the present design ensured that the PCR chamber radius and height were sufficiently small that the vertical gradient could not lead to temperature variations greater than ±1°C). The FEM model by itself was verified by extensive experiments, wherein thermochromic liquid crystal (TLC) sensitive to different temperatures was deposited in the chamber to read the temperature within its volume in steady-state conditions.
3. Heat Distribution Model
3.1. Problem Formulation
The partial differential equation describing temperature distribution in the chosen domain at any time is given by the heat conduction equation as follows [30]: where is the temperature of every point of the domain at any given time. is thermal conductivity, is density, and is heat capacity. Since the chosen model features axial symmetry, the temperature in the cylinder does not depend on . This means that . Therefore can be dropped and the Laplace operator can be written in a 2D cross-section domain using and as follows:
To avoid any confusion between and the notation used for density, we introduce here notation as follows: which is thermal diffusivity of the material in . Finally the partial differential equation giving the temperature at the half part cross-section of our cylindrical domain with and as its horizontal and vertical axes can be written as
The definition of our domain implies that and are limited to and . The initial and boundary conditions are given by where is the initial distribution of temperature. is the temperature of the outside wall of the cylinder. The boundary condition (7) comes from the axisymmetrical cylindrical structure of the model. The boundary condition (8) is based on the assumption that the upper surface of the chip is isolated. describes the geometric distribution of the control action signal at the bottom surface, which is defined by the heater location. A general solution of (4) can be found by adding the solution of homogeneous and nonhomogeneous equations describing the steady-state and a transient temperature distribution, respectively:
Substituting (10) in (4) we get
Without the loss of generality, we can assume as we defined to describe the steady-state temperature. By applying (10) to the initial condition (5) and boundary conditions (6), (7), (8), and (9), we arrive at two PDEs:
The first one is an initial value problem with homogeneous boundaries, while the second one is a nonhomogeneous boundary value problem.
3.2. Steady-State Part
To solve (13), we employ the method of separation of variables [31]: where and are functions of and , respectively. The constant part is used to normalize the boundary conditions and to achieve the homogeneous boundary condition. Using (14) in (13) and dividing we arrive to
The left side contains functions of alone, while the right side contains functions of alone. Since this equality must hold for all and in the given interval, the common value of the two sides must be a constant; say varying neither with nor with . The constant here was chosen as a negative value to avoid a trivial solution [31]. Now we have two ODEs for the two factor functions:
The boundary conditions of (13) can also be stated in the product form (14) which results in the following:
The solutions of (16) are given by [31] where is the Bessel function of the first kind of order and is the Bessel function of the second kind of order . From the boundedness condition at stated by (18), we must have because is unbounded at zero. Thus the solution (21) becomes
From the boundary condition (19) and the fact that we see that which implies , so the solution (22) becomes
The boundary condition (17) can be employed to determine : which can be satisfied only if so that where is the th positive root of . Thus ’s can be expressed as
So a solution satisfying the boundary conditions is where and . The general solution of (13) can be obtained by employing the superposition principle [31]. By replacing with and summing up all eigenfunctions, it follows that
From and the last boundary condition (20) we can write
The right side is the Bessel series expansion of with coefficients equal to [31], so we get
3.3. Transient Part
We employ the separation of variables again to solve (12) where , , and are functions of , , and , respectively. Following the same procedure as the steady-state part and defining , , and as eigenvalues for , , and , respectively, we arrive at the following ODEs: where . By applying (35) into (12), the boundary conditions become
The solutions of (36), (37), and (38) are given by
From the boundedness condition at expressed by (40), we must have . Thus the solution of (43) becomes
From the boundary condition (41), we have which implies , so the solution (44) becomes
From boundary condition (39) we can determine . Consider so we have where is the th positive root of . A solution of (36) satisfying the boundary conditions is where . From the boundary condition (42) we can determine . We have which can be satisfied only if or Using (53) in (48) and the fact that the solution (44) becomes
Now from (53) and (50) it follows that
So a solution satisfying all boundary conditions is given by where , , and . Replacing with and summing up and we obtain the solution by the superposition principle
From the initial condition we have
This can be written as where is defined by (34). To simplify the equation, without the loss of accuracy, we consider that which is reasonable because we start the experiment when the system is in a uniform temperature distribution. We can write (60) in this form
Clearly because of orthogonality of the Bessel functions, these coefficients have to be equal; that is,
Written in standard form, we have
Recalling the Fourier cosine series expansion [31], the left side of (63) is the Fourier cosine series expansion of the right side with respect to . Thus substituting and replacing and with and , respectively, we have
So for are given by
Equation (65) can be written as and solving the integral we have
Since takes only positive integers, we can simplify (67) and arrive to
Thus we have and is defined by (68) as where . Hence, it follows that the complete solution to the problem is given by (10), (31), and (58), where their coefficients are defined by (34), (69), and (70).
3.4. Stability Analysis
Stability analysis of the calculated model can be easily performed by studying the transient part of the general solution. A quick visit of transient equation (58) reveals that the infinite poles or eigenvalues of the system are given by the coefficients of the exponential term in infinite series as follows: where and . Because thermal diffusivity, , takes only positive real values for different materials, all the poles are negative real poles, and therefore the system is stable. Calculating the slowest poles of the system using (71) and the parameters given in Table 1, we arrive to the following results:
4. Simulation Results
The parameters with the values listed in Table 1 are used in simulation. We assume that C and that an input power of one Watt is being applied to the heater. The input to the system is given by the heat flux transferred from the heater to the system. This heat flux is considered in definition of which by itself determines the coefficients. Definition of is given based on the chosen model configuration by the following equation: where , , is heat flux from the heater to chip, and is the thermal conductivity of the chip. For discontinuities we have . Now for we can write
For heat flux , we have where is the amount of the power transferred from the heater to the chip and is the area of the heater which can be calculated as follows:
Recalling the recurrence relations for the Bessel function [31], we have
For , (77) can be written as follows:
Using a change of variable, , for the integral part of we can write
Taking into account (75) and (80) in (74), we have
Now we have all the required equations to calculate the temperature distribution using our analytic solution.
The next important step before start running simulations is how to compensate the error caused due to the simplification of the actual system into a single-domain model with Borofloat glass as the only material involved. Substituting the PDMS layer with a glass layer could result in a significant amount of error and mismatch between model and the actual system. To overcome this issue, we constructed an optimization problem to optimize our model parameters arriving to a best fit between single domain analytic model and multidomain actual system. Chip radius, chip thickness, and input power in the mathematical model are selected as the optimization degrees of freedom. We focus our attention on the two regions that have the highest importance. The first region is a rectangular area which starts from the left bottom side of the chamber and continues up to the top of the constant temperature wall and is all made from glass in the analytical model. We called this region comparison aperture. The second region which is a subset of comparison aperture and has higher importance is reaction chamber. To get quantitative measure of the error in both the chamber and the selected aperture, we define the following parameters in both regions: absolute average error, maximum absolute error, and uniformity of the error. The definitions of the first two measures are obvious from their naming, and the last measure is defined as follows: where is error at selected aperture indicated by and is error at chamber region indicated by . Next we define our weighted-sum objective function including norm of the error at chamber location, uniformity measure at chamber location, and uniformity measure at selected aperture as follows: where , , and are weight coefficients and is area of chamber region. By selecting proper weights, optimization results in a maximum absolute error of 0.14°C in chamber temperature, which is less than 0.03 of the error when compared to the results of a nonoptimized model ( mm, mm, and W), satisfy the acceptable error criteria for reliable PCR operation. Error distribution in whole aperture area is depicted in Figure 8. It is clear that, although the maximum absolute error reaches 3.6°C in the area close to the heater, the error in the chamber area is kept lower than 0.14°C.
The resulted steady-state temperature distributions based on optimized analytical solution and the FEM simulation are shown in Figures 7 and 6, respectively. Comparison aperture is indicated by dash lines. The solutions of the two results show an excellent degree of correlation, verifying accuracy of the steady-state solution.
Dynamic simulation of temperature evolution inside the reaction chamber is the next important result that we would like to study about our analytic model. Precise FEM simulations disclosed that lower center point of the chamber and upper edge of the chamber have the highest and the lowest temperatures inside the chamber, respectively (SW and NE edges of the chamber in Figure 5). The results follow the fact that upper center point of the chamber (NW edge of the chamber in Figure 5) has a temperature close to mean temperature of the chamber. We selected upper center point of the chamber to study the dynamic evolution of temperature in our model and in actual system. An input power which consisted of three predefined values that are required to take the chamber temperature to denaturation, annealing, and extension temperatures is applied to both models. Both systems were powered up after being settled in ambient temperature and kept in each power level for 60 seconds. The experiment ran over a whole PCR cycle in a 240-second interval in open-loop condition. The results are depicted in Figure 9. The optimization was performed for the lowest temperature in PCR cycle. However the error at the highest temperature is still under 1°C. Model dynamics nicely follows the FEM results, and as we expected, the transient response is about 1% faster than that of the actual system, due to the difference between the thermal properties of PDMS and Borofloat glass.
5. Conclusion
This paper presents the derivation and the solution of the heat distribution equation of a microfluidic chip designed and manufactured at the Applied Miniaturization Laboratory at the University of Alberta, Canada. The partial differential equation featuring three-dimensional domain was used to describe the heat distribution of the chip during a DNA amplification process. The axisymmetric architecture of the chip helped reduce the domain to two-dimensional cylindrical form. The PDE equation describing the heat distribution of the chip was divided into steady-state and transient equations and then solved explicitly. Precise FEM simulations in addition to laboratory measurements were used to verify the developed model. The modeling approach is considered to be appropriate for similar thermal processes with collocated sensors and actuators. Future work is being carried out toward developing closed-loop tracking control based on boundary control techniques.
A dynamic simulation of the system was made using the developed model to confirm that the scaling of the glass thermal properties could be adjusted to match the temperature predicted by the FEM simulations at all points in the PCR chamber with less than 1°C error at all times during a PCR temperature cycle.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors would like to thank S. Poshtiban, G. V. Kaigala, and J. Jiang for many helpful technical discussions.