Numerical framework for simulating bio-species transport in microfluidic channels with application to antibody biosensors

Diagnosis is a fundamental stage in health care and medical treatment. Microfluidic biosensors and lab-on-a-chip devices are amongst the few practical tools for achieving this goal. A new computational code, specifically for designing microfluidic-integrated biosensors is developed, the details of which is presented in this work. This new approach is developed using control-volume based finite-element (CVFEM) method and solves bio-recognition chemical reactions and full Navier–Stokes equations. The results of the proposed platform are validated against the experimental data for a microfluidic based biosensor, where excellent agreement is achieved. The properties of the biosensor, sample, buffer fluid and even the microfluidic channel can easily be modified in this platform. This feature provides the scientific community with the ability to design a specific biosensor for requested point-of-care applications.• A new approach is developed using control-volume based finite-element (CVFEM) method for investigating flow inside a microfluidic-integrated biosensor. It is also used to study the influence of surface functionalization on binding cycle.• The proposed model solves bio-recognition chemical reactions as well as full Navier–Stokes and energy equations. Experimental-based or personalized equations of the chemical reactions and flow behaviour are adoptable to this code.• The developed model is Fortran-based and has the potential to be used in both industry and academia for biosensing technology.


Introduction
Biosensors are affordable real-time monitoring devices that make dynamic monitoring of specific targeted molecules (environmental, agricultural or body micro samples) possible. According to the world health organization, in the past two decades, ischemic heart diseases, stroke, chronic obstructive, pulmonary disease, Alzheimer, lower respiratory infections are top five global causes of death [21] , while they can be cured if they are detected at early stages. These examples and a lot of deadly contagious diseases, cancers and environmental causes can be detected and analysed by biosensors which would make them one of the most interesting devices in the 21st century [8] . Microfluidic-integrated biosensors are functionalized surfaces installed in a channel with the scale of micron. A schematic view of the detection process in these biosensors is illustrated in Fig. 1 . Targeted molecules, which can be coronavirus or any kind of virus, DNA or protein that is needed to be detected or measured, is first mixed with the buffer fluid, which helps with a better detection of the targeted molecules. Then the buffer fluid flows into the microfluidic channel, and as it passes the functionalized surface, chemical reactions take place at this stage. The transducer converts these chemical reactions into measurable signals, which show the presence of the targeted molecule in the sample and its amount. The binding cycle of this detection process is presented in Fig. 1 , which shows that in the first stage, a small portion of biological samples is inserted into the microfluidic channel (e.g. body sweat which contains ions, metabolites and hormones). In the next stage the samples react with the antibodies while passing by the surface of the biosensor. As saturation takes place in the third stage, a wash out flow goes into the channel (stage 4) and completes the binding cycle. Signals of these binding cycles get processed and reveal the results of the detection procedure.
For the purpose of designing biosensors and improving their efficiencies, various experimental works and limited number of numerical analyses have been conducted [4 , 16,17 ]. In recent years, computational fluid dynamics has effectively been applied to various biomedical-related projects involving design, validation and proof-of-concept [2 , 6 , 7 , 11-13 , 18] . Previous numerical analyses are mainly based on finite element methods (FEM) or finite volume methods (FVM). In the FVM approach, diffusion coefficient is reduced by 10 −10 in order to reduce numerical instability [3] . On the other hand, in the FEM approach, the analytical velocity profile is used instead of solving the whole Navier-Stokes equation. Both cases lead to an almost non-realistic detection time. In the previous works, different software packages and codes were utilized to model the flow inside the microfluidic channel. The discretization method is a key component in these numerical solutions. Table 1 presents three different discretization methods; finite element method (FEM), control-volume method (CVM) and control-volume based finite-element method (CVFEM). FEM has geometrical flexibility, CVM has physical intuition and CVFEM is a powerful combination of these two methods [19] .

Methodology
In the present work, a new approach is used for conducting numerical simulations in this field which provides reliable and realistic binding cycles. Results have been validated with experimental data in the companion publication by the present authors [15] . The main challenge in numerical simulation is achieving a realistic modelling of the transport of bio-species and their reactions on the functionalized surface [10] . In order to overcome this challenge, in the current model, full Navier-Stokes equations coupled with convection, diffusion and reaction of targeted molecules Eqs. (1) -( (3) ) are solved implicitly with control-volume based finite-element method (CVFEM). The coupled system of equations is discretized using high order discretization methods such as the improved skewed upwind differencing (SUD) scheme. Details of the numerical domain and boundary conditions that are implemented in this model are summarized in Fig. 2 and Table 2 . For concentration ( c ) of bio-species on the sensor and on the walls, 'Homogeneous Neumann' and 'Neumann' conditions are prescribed, respectively. In the present code, unsteady viscous fluid flow in a channel is solved and the flow is assumed to be laminar, incompressible and 2D. The flow is solved using the finite element-based control volume method.

Governing equations
The proposed numerical platform solves the Navier-Stokes equations Eqs. (1) -(3) to simulate the buffer fluid flow and Langmuir-Hinshelwood model ( Eq. (11) ) to simulate the kinetics of adsorption of molecules. The fluid is considered as continuum and incompressible [20] , due to negligible changes in density.
ρ ∂v ∂t where u and v are the velocity in the x and y direction, respectively, while ρ is the density and μ is the molecular viscosity. Continuity equation has the role of controlling pressure field in the domain, but as it is shown, there is no pressure term in continuity equation. In order to relate the pressure field effects in continuity equations, the Karimian-Schneider [5] method ( Eqs. (4) to (6) ) is used.
where BL refers to the values in boundary layer, s is the area of each element, L is the reference length and V j is the velocity of the main element. In the convection term of momentum, u and v are obtained using the physical influence scheme (PIS) method ( Eqs. (7) and (8) ).
For simulating the kinetics of adsorption of molecules on a functionalized solid surface (e.g. DNA hybridization and virus detection etc.), the Langmuir-Hinshelwood model [1] is used. Since targets are constantly captured by ligands and dissociate at a smaller rate, such reactions are weakly reversible. Using Fick's law ( Eq. (9) ) and evaluating the mass balance and effect of velocity field would give us the complete equation system for modelling the concentration ( Eq. (10) ). The reaction rates can be written as Eq. (11) . F = −D ∇c (9) ∂c where c is the concentration of the targeted molecules, c 0 is the inlet concentration, b is the surface concentration of bound analyte [14] , b max is the density of binding sites on the sensor, D is the diffusion coefficient and S is the source/sink term.

Discretization
Since the control-volume based finite-element method is used in the present computational code, equations are discretized in a conservative form in order to satisfy the conservations through all surfaces. For implementing the control-volume based finite-element method, first the numerical domain is discretized into nodes, elements and sub-elements, as it is shown in Fig. 3 . Each element consists of four sub-elements.
For discretization, the integration of the main equations Eqs. (1) -(11) on each sub-element (control volume) is calculated. To solve the above-mentioned equations, velocity and pressure are needed at integration points (see Fig. 3 ). Different interpolation methods are used for generating values in integration points, since interpolation does not take the physics of the flow into account and it is not reliable for the entire range of Peclet numbers. This case is more evident in flows with high Peclet numbers. In this code, the variables in integration points are calculated with the PIS upwind method. For this method, the required variables in upstream points of flow are generated. With the values of velocity in x and y directions and coordinates of upwind points, the equation of the streamline passing by these points can be generated ( Eq. (12) ).
where a and b represent the portions of an element split by the streamline, as it is shown in Fig. 3 and subscripts 1 and 4 indicate the number of sub-element. Subsequently, a 4 × 4 matrix is generated for each element ( Eq. (14) ), which connects the value of upwind point to 4 nodes of each element. With these values, flow can be solved by the upwind method.
The area of each element and sub-element ( ds x and ds y ) needed in the calculations of integrations over control volume are generated within the subroutine 'Area'. At the end, for solving the main equations, the value of each element is obtained by a bilinear interpolation between the 4 subelements. Details of these procedures are described in the next section.

Structure of the code
The proposed numerical model can solve the microfluidic-integrated biosensor with different geometries, sensor, position of the sensor, buffer fluid and meshes. The list of these parameters is presented in Table 3 . Fig. 4 and Table 4 show the structure of the proposed numerical model.

Method validation
This model has been validated with experimental data provided by Berthier and Silberzan [1] . Experiments were done in a micro chamber with the details provide in Table 5 for different DNA strands. The buffer fluid enters the chamber while carrying the targeted molecules. This experiment takes 10 5 s. The results are shown in Fig. 5 , where good agreement is obtained, demonstrating a success in using the present methodology for modelling convective diffusive and reactive targeted molecules inside microfluidic integrated biosensors. The developed numerical model is suitable for modelling flow inside microfluidic-integrated biosensors, which is adaptable to various Table 4 The structure of the numerical model.  geometries and operating conditions. Although this version of the simulation platform has several limitations, based on the boundary conditions and assumptions defined in the code, it is capable of improvements can be made by through implementing more applicable numerical models for future studies. For instance, the concentration of the targeted molecules is assumed to be high enough for a bulk concentration, although this assumption cannot be valid for samples with low analyte concentrations. The performance of this numerical platform in other favourable applications in analysis of miniaturised biological devices such as indirect measurement of chemical detection [9] , could be investigated and enhanced in future versions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.