3D simulation of microfluidic biosensor for SARS-CoV-2 S protein binding kinetics using new reaction surface design.

In this study, we performed 3D finite element simulations on the binding reaction kinetics of SARS-CoV-2 S protein (target analyte) and its corresponding immobilized antibody (ligand) in a heterogeneous microfluidic immunoassay. Two types of biosensors with two different shapes and geometries of the reaction surface and electrodes were studied. Alternating current electrothermal (ACET) force was applied to improve the binding efficiency of the biomolecular pairs by accelerating the transport of analytes to the binding surface. The ACET force stirs the flow field, thereby reducing the thickness of the diffusion boundary layer, often developed on the reaction surface due to the slow flow velocity, low analyte diffusion coefficient, and surface reaction high rate. The results showed that the detection time of one of the biosensors can be improved by 69% under an applied voltage of 10 Vrms and an operating frequency of 100 kHz. Certain control factors such as the thermal boundary conditions as well as the electrical conductivity of the buffer solution were analyzed in order to find the appropriate values to improve the efficiency of the biosensor.

Introduction

The severe acute respiratory pandemic of coronavirus 2 (SARS-CoV-2) which was reported in 2019 in Wuhan, China [1], causes the coronavirus disease (COVID-19) responsible for the death of many people across the world. Present information indicates that SARS-CoV-2 is more contagious than the former severe acute respiratory syndrome coronavirus SARS-CoV [2]. SARS-CoV-2, like SARS-CoV, uses the angiotensin converting enzyme (ACE2) to bind to human cell through its structural S-spike glycoprotein (S protein) [3].

Among the methods of diagnosing COVID-19 infection, the real-time polymerase chain reaction (real-time PCR) is known to be an efficient and sensitive technique [4] despite several false-positive or false-negative results that can sometimes occur, especially in the early stages of the infection. On the other hand, laboratories using real-time PCR methods need sample preparation, expensive laboratory tools, heating sample to several temperatures for denaturation, hybridization, and extension which limits their applications [5, 6].

For later stages, 5 to 15 days after infection with the virus, antibody-based techniques, named serological assays, can be used such as enzyme-linked immunosorbent assay (ELISA). However, false-positive results can also happen due to cross-reactivity between specific antibodies with antibodies versus other coronavirus epitopes [7]. Although the serological assay is fast and needs minimal equipment, its effectiveness may be restricted only in detecting an acute COVID-19 infection.

Despite some drawbacks, point-of-care (POC) diagnostic tools, based on optical, colorimetric, electrochemical, and magnetic methods, remain promising methods for decentralized, rapid, sensitive, and inexpensive diagnosis of COVID-19 infection [8].

To eliminate the complicated steps in sample preparation and diminish the possibility of false positives and false negatives, much research has focused on innovative detection methods such as biosensors. In the last few years, there has been great concern in biosensors, which are transportable devices, made of biological molecules with a sensor [9].

In the biomedical field, biosensors are used to detect target analytes such as proteins or biological molecules using the sensitivity of antibody-antigen interactions. Two types of immunoassays can be used, the homogeneous immunoassay, where antibodies and antigens take place in solution in the microfluidic channel, and the heterogeneous immunoassays where the interaction takes place between the antibodies immobilized on the sensitive surface of the microfluidic channel and the antigens in solution just at the boundary layer of the sensitive surface. The binding reaction between the target analyte (antigen) and the fixed ligand (antibody) leads to the formation of the analyte/ligand complex. The surface concentration of this formed complex has a crucial role for the detection process [10, 11].

The generation of fluid flow in these microfluidic devices is a major challenge due to the dominance of surface forces over inertial forces (i.e., low Reynold number). The ratio of the diffusion rate of the analytes to the surface reaction is comparatively small, and a diffusion boundary layer is developed and affects the efficiency of the biosensor [12, 13].

Recently, the number of improved microfluidic biosensor devices has increased significantly [14]. Saad et al. used a magnetic field to improve the binding reaction [15]. Alternating current electrokinetic technology (ACEK) is typically used to manipulate fluids in microfluidic biosensor devices [16-18]. Among the phenomena of alternating current electrokinetics, alternating current electrothermal force (ACET) is found to be the popular technique used to mix fluids and improve the binding efficiency of diffusion-limited molecules in microfluidic devices [19-22].

The applied voltage, which leads to a rise in the temperature of the fluid, has an upper limit, beyond which the electrochemical effects can damage the nature of the biofluids [23-25]. As the electrothermal force and therefore the velocity of the fluid are extremely dependent on the applied voltage (i.e., quartic dependence on Vrms), the increase in the velocity of electro-thermally based biomedical microdevices is consequently hampered.

To overcome this problem, several studies have proposed and developed innovative designs of microfluidic biosensors to modify the electric field and consequently flow topology by changing the geometry of the electrodes as well as the shape of the reaction surface [26, 27]. In this context, two different geometries of electrodes and reaction surface have been proposed to improve the performance of biosensors for the detection of specific proteins of the coronavirus.

Kaziz et al. [22] analyzed the ACET effect on the binding kinetics of the SARS-CoV-2 S protein. The location of the electrodes on the walls of the microchannel and the appropriate choice of the adsorption and desorption constants of the binding reaction, whereas conserving the same affinity, are discussed in order to reduce the response time of the biosensor and thus improve its performance. However, these numerical simulations have been done for two-dimensional configurations which is somewhat constraining for spatial factors and practical design, especially the sensing area. In the present work, three-dimensional simulations are performed. The effect of the sensing area shape and the geometrical location of the electrodes on the microfluidic biosensor performance are investigated. Several parameters such as thermal boundary conditions and the effect of electrical conductivity were also tested in this work.

Theoretical formulation

Geometrical configurations

Before giving the equations of the model, we start by specifying the geometry of the investigated biosensors. We have considered two types of biosensors with distinct designs of reaction surface and electrodes as shown in Fig. 1. The length, width, and height of the microchannel used for the two types of biosensors are, respectively, 250, 60, and 40 µm.

Fig. 1

Design of two types of biosensor models. a 3D microchannel; b type 1 biosensor; c type 2 biosensor

For the first type of biosensor, the reaction surface was a disk (reaction surface 1), while for the second type it was a ring (reaction surface 2). For both biosensor types, the thickness of the reaction surface is 3 μm and that of the electrodes is neglected. The reaction surface for both types of biosensors is in the bottom wall and at a distance X = 100 µm from the microchannel inlet. All the geometric parameters of the type 1 biosensor and the type 2 biosensor are given in Tables 1 and 2, respectively [26].

Table 1

Geometrical parameters for the type 1 biosensor (the reaction surface is a disk)

Parameters	Value
External radius of electrode Rext (µm) Internal radius of electrode Rint (µm) Radius of reaction surface RS (µm) Bottom electrode angle α (°) Top electrode angle θ (°) Distance from the inlet X (µm)	25 20 15.96 160 160 100

Table 2

Geometrical parameters for the type 2 biosensor (the reaction surface is a ring)

Parameters	Value
Positive electrode radius REl1 (µm)	10
Internal radius of reaction surface RSint(µm)	15
External radius of reaction surface RSext(µm)	22
Internal radius of the negative electrode REl2int(µm)	25
External radius of the negative electrode REl2ext(µm)	30
Distance from the inlet X (µm)	100

The simulation was performed by calculating the electric field, the temperature field, the flow field, and the analyte concentration field within the microfluidic channel and calculating the reaction kinetics of the SARS-CoV-2 S protein/antibody pairs formed on the reaction surface of the biosensor. All calculation details are presented in the next sections.

Electric field calculation

The electric field E was calculated using the Poisson equation (Eq. 1) where V is the electric potential due to the external potential applied across the electrodes.

Temperature field calculation

The temperature field in the microchannel fluid was calculated by means of the following thermal energy equation (Eq. 2):where ρ, λ, σ, C and are the density, the thermal conductivity, the electrical conductivity, the specific heat at constant pressure, and the velocity field of the fluid, respectively.

The non-uniform alternating current (AC) electric field applied on the fluid induces the variations of the electrical conductivity σ and permittivity ε of the solution due to the temperature gradient generated from the inhomogeneous Joule heating serving as the heat source of the fluid, and it is defined as .

Although the electrothermal force applied causes heating of the fluid, the dependence of the specific heat at constant pressure and of the thermal conductivity of the fluid as a function of the temperature T can be neglected, since the conduction of the ambient heat allows this received heat to dissipate, and therefore, the growth of the temperature in the fluid mostly stands low [24].

Here, the viscous dissipation term is smaller than the Joule effect [18], so it has been neglected.

Flow field calculation

The laminar flow velocity field of the incompressible fluid in the microchannel was calculated using the three-dimensional stationary Navier–Stokes equations (Eqs. 3 and 4). defines the electrothermal force.

Here, p is the pressure and the dynamic viscosity of the fluid. The applied non-uniform electric field causes non-uniform Joule heating of the fluid inside the microfluidic channel which gives rise to a temperature gradient and therefore to gradients of permittivity and electrical conductivity, thus generating the electrothermal force responsible of the agitation of the fluid. The expression of the electrothermal force is given by (Eq. 5):where is the fluid charge relaxation time and is the angular frequency of the electric field . The used transporter fluid is water at a temperature close to 293 K, which allowed to consider, according to Green et al. [23], the following equations (Eqs. 6 and 7):

Therefore, the electrothermal force is given by (Eq. 8):where is the permittivity of free space and is the relative permittivity of the fluid.

Analyte concentration field calculation

The reaction kinetics between specific analytes and attached ligands to the sensitive surface of the biosensor can be described by two processes [28]: the process of mass transport of analytes by convection and diffusion to the reaction surface and the process of the chemical reaction of adsorption of analytes to the ligands to form analyte/ligand complexes on the reaction surface. Since the concentration of target analytes in the carrier fluid is too low, it therefore has no effect on viscosity [29]. By using Fick's second law and the impact of the velocity field, a full convection–diffusion equation was obtained to model the transport of targeted molecules (Eq. 9):where indicates the concentration of the analyte and is the diffusion constant. Obviously, the above equation does not contain any generation term. The binding reaction is confined on the reaction surface, ant it is involved only on the boundary condition.

Reaction kinetics calculation.

For heterogeneous biosensors, chemical kinetics are calculated using the Langmuir–Hinshelwood first-order adsorption model [29, 30]. On the reaction surface, antibody ligands (B) trap target analytes (A) to create analyte/ligand complexes (AB), and then, these complexes dissociate at a slower rate according to the following equation (Eq. 10):where is the adsorption rate constant, is the dissociation rate constant, and and are, respectively, the concentration of analytes at the binding surface and the concentration of ligands immobilized on this surface. is the concentration of the analyte/ligand complexes formed at the binding surface.

Boundary and initial conditions

Table 3 summarizes all the boundary conditions applied in this model. An electric potential of and a temperature equal to that of the ambient have been applied to the electrodes.

Table 3

Boundary conditions—electric potential, temperature, velocity and analyte concentration for walls, reaction surface, electrodes, inlet, and outlet for both biosensor types

Type	Electric potential (V)	Temperature (T)	Velocity (u)	Concentration ([A])
Microchannel walls	∂V∂n=0	∂T∂n=0	u = 0	∂A∂n=0
Reaction surface	∂V∂n=0	∂T∂n=0	u = 0	∂A∂n=-1D∂AB∂t
Electrodes	±V_rms	T_0	u = 0	∂A∂n=0
Inlet	∂V∂n=0	n→.λ∇T=0	uave	A_0
Outlet	∂V∂n=0	n→.λ∇T=0	∂u∂n=0	n→.(D∇A)=0

is the unit normal vector to the surface

The walls, reaction surface, inlet, and outlet are assumed to be electrically insulated. In the matter of the thermal boundary conditions, the walls and reaction surface are presumed to be thermally insulated. No heat flow is exchanged at the inlet and outlet of the channel.

Concerning the laminar flow, the inlet fluid flows with a parabolic profile such that the average velocity, u, was set at 100 μm/s, and the outlet was open to the atmosphere. A no-slip condition is applied at the walls, reaction surface, and electrodes. For the transport of the analytes, a constant analyte concentration and the condition were imposed at the inlet and at the outlet of the channel, respectively. Since there was no interaction between the target analytes and the electrodes or the other walls, the homogeneous Neumann condition was used [31]. In other terms, these surfaces are assumed to be impermeable. On the reaction surface, the diffusive flow condition was then imposed because there was a loss of analyte concentration due to the analyte/ligand adsorption. In other words, at the reaction surface, the diffuse flow is balanced temporal flow rate. This equation is valid since the Langmuir model is a monolayer model and the stoichiometric is equal to one: .

Concerning the initial conditions related to Eqs. (9) and (10), the analyte concentration and the complex concentration have been set to zero: and .

Numerical simulations

The system of equations describing this model was solved numerically with the Galerkin finite element method [32]. The whole fluid domain has been meshed with 9337 tetrahedral elements with a refinement in the vicinity of the reaction surface area. All surfaces of the microchannel have been meshed with 2210 triangular elements with a finer mesh for the electrodes. The velocity components are approximated using second-order interpolation and first order for the pressure in each element, whereas the analyte concentration is approximated using quadratic interpolation in each element.

First, the electric field has been calculated by means of the electrostatic equation, and then, the fields of temperature, pressure, and velocity are obtained by solving together the Navier–Stokes equations and the thermal energy. Finally, the two equations depending on the time that are the transport of analytes and the chemical binding reaction equations are calculated simultaneously to determine the concentrations of the target analytes, [A], in the microchannel and the analyte/ligand complex concentration ([AB]) on the sensitive surface.

Here, the analyte is the SARS-CoV-2 S protein, and the ligand is its corresponding antibody (b1 and h12) [31, 33]. The total concentration of the protein/antibody complex formed was obtained by integrating its local concentration over the entire binding surface, and the simulation of the dissociation step was obtained by stopping the supply of analyte after the saturation time. The normalized surface concentration of SARS-CoV-2 S protein/antibody,, is the surface concentration, which is normalized to the total binding sites on the sensor surface, .

All the physical and biological (binding) parameters of SARS-CoV-2 S protein/antibody [31, 33] used for this numerical study are illustrated in Table 4.

Table 4

Physical and binding parameters

Constant	Name	Value
λ (W/(K·m))	Thermal conductivity	0.6
ρ (kg/m2)	Fluid density	1000
µ (Pa.s)	Dynamic viscosity	1.08 × 10–3
Cp (kJ/(kg·K))	Specific heat	4.184
σ (S/m)	Electrical conductivity	5.75 × 10–2
εr	Relative permittivity	80.2
f (kHz)	Frequency	100
k_on(m3/ (Mol.s))	Adsorption rate constant	103
k_off(s−1)	Desorption rate constant	10–3
D (m2/s)	Diffusion constant	10–11
B_max(Mol/m2)	Surface ligand concentration	3.3 × 10–8
A_0(Mol/m3)	Analyte input concentration	10–9

Results and discussion

Numerical model validation

Initially, the model was proven compared to experimental data point of Berthier and Silberzan [29], without applying AC electric field (Fig. 2). The used microfluidic channel is 1 cm wide and 1 mm high. The flow rate of the carrier fluid is 10–6 m3/s. The concentration of the target molecules at the entrance of the channel and the diffusion constant are, respectively, 2.5 × 10–6 Mol/m3 and 7 × 10–11 m2/s. The density of binding locations, association, and dissociation constant are 1.668 × 10–8 Mol/m2, 75 m3/Mol⋅s, and 10–2 1/s, respectively. We noticed that our results were in good agreement with the experimental data. The difference between the numerical and experimental results comes from numerical model, but also experimental results contain some errors. However, this difference is relatively small, and we can consider that the model is validated.

Fig. 2

Validation of the current numerical model compared to the experimental study of Berthier and Silberzan [29], without electrothermal force

Following the successful experimental validation, we also validate the model with the numerical results of C-reactive protein (CRP) complex concentration calculated by Echouchene et al. [26] for the two biosensor types without (Fig. 3a) and with (Fig. 3b) applying AC electric field. A well concordance between the two results was also observed especially during the adsorption phase. A small error is observed at the end of desorption phase. Overall, the average error is smaller than 5%. We can consider that the model is validated, and it may be used for further applications.

Fig. 3

Validation of the current numerical model versus the numerical data of Echouchene et al. [26] with and without electrothermal force for both biosensor types. a Without electrothermal force and b with electrothermal force

Applied voltage effect on binding reaction kinetics

In this section, different AC voltages are applied to determine which voltage produces the greatest permitted temperature rise, which is about 3 K, for biochemical process conditions [18]. Figure 4 shows the evolution over time of the normalized surface concentration, , for the two types of biosensors as a function of different values of applied voltage. A strong improvement in the binding reaction was spotted in the case where the electrothermal force was applied for both types of biosensors. This is the essential proof of the efficacy of ACET in improving immunoassays of biosensors. Certainly, the electrothermal force generates vortices, thus stirring the flow, which enhances the transportation of analyte molecules to the binding surface and thus reduces the diffusion boundary layer thickness. As can be seen, the limit of the applied voltage is 20 V for the first type of biosensor and 10 V for the second type because the temperature rises exceed 3 k beyond these limit voltages.

Fig. 4

Normalized surface concentration as a function of time, with and without applying voltage for the two biosensor types. a Type 1 biosensor; b type 2 biosensor

To estimate the ACET impact on the kinetics of the binding reaction of SARS-CoV-2 S protein, we registered the detection time of the biosensor, drop percentage, maximum velocity, and temperature rise for both types of biosensors when applying various voltage values for both studies. The detection time is the biosensor response time when the surface concentration achieved 95% of its maximum value. Results are listed in Tables 5 and 6 for type 1 and type 2 biosensors, respectively. We can see the decrease in the detection time of both biosensors with the application of the electric field. Consequently, the ACET force enhances the kinetics of SARS-CoV-2 S protein/antibody ligand binding reaction and, therefore, the performance of the biosensor. We also see that the highest temperature rise is smaller than 3 K for applied voltages less than or equal to 20 V for the type 1 biosensor and 10 V for the type 2 biosensor, which limits their use beyond these optimal voltage values. The temperature rise recorded for the type 1 biosensor is 7.06 K for an electrical excitation of 25 V and that of the type 2 biosensor is 5.76 for 15 V applied voltage.

Table 5

Type 1 biosensor response time, drop percentage, maximum velocity, and temperature rise for various applied voltages

Applied Voltage (V)	Detection Time (s)	Drop percentage (%)	Maximum velocity (mm/s)	Temperature rise (K)
0 V	28,975	–	0.19	–
5 V	27,400	5	0.28	0.28
10 V	27,350	6	1.59	0.29
15 V	23,200	20	5.24	1.13
20 V	17,450	40	13.48	2.55

Table 6

Type 2 biosensor response time, drop percentage, maximum velocity, and temperature rise for various applied voltages

Applied voltage (V)	Detection time (s)	Drop percentage (%)	Maximum velocity (mm/s)	Temperature rise (K)
0 V	21,875	–	0.20	–
5 V	15,925	27	1.51	0.64
10 V	6875	69	16.80	2.57

Figure 5 shows the velocity field in both types of biosensors with the application of electrothermal flow. As expected, it is observed that the velocity profile is no longer parabolic in the vicinity of the reaction surface and the electrodes, giving rise to a flow structure which is not perfectly unidirectional. The fluid velocity increases near the reaction surface, in both types of biosensors, and reaches the highest values recorded in Tables 5 and 6. This certainly indicates the presence of eddy structures, near the reaction surface and the electrodes, responsible for the disturbance and agitation of the fluid and therefore for improving the speed of the biosensor binding reaction.

Fig. 5

Velocity field with ACET effect for both biosensor types. The plotted plane is x—y for z = 0 μm

For better examination of changes in biosensor performance, detection time versus applied voltage is investigated for both biosensors. As shown in Fig. 6, the detection time decreases when the applied voltage increases. Variations in voltages occur at the same rate, while variations in binding cycle and detection times are not at the same rate. The relationship between detection time and applied voltage is nonlinear.

Fig. 6

Polynomial adjustment of the response time as a function of the voltage for the two types of biosensors

Equations (11) and (12) show the best-fitting polynomial line across the numerical results for type 1 and type 2 biosensor, respectively:where τ is the detection time and V is the applied voltage.

The interest of the fittings (11) and (12) linking detection times to voltage is to predict approximately the detection time of each biosensor for any given applied voltage (within the permitted voltage range).

Temperature boundary conditions effect

In this subsection, type 2 biosensor was used since it has the best detection time for a lower applied voltage. To show the influence of thermal boundary conditions on the fluid temperature field and on the efficiency of the binding reaction, four cases of adjustment of the thermal limit applied to the electrodes and to the walls of the microchannel are considered. For case 1, a temperature of 300 K was applied for the two electrodes only and the upper and lower walls of the microchannel were thermally insulated. This was the case used in previous simulations. For case 2, the electrodes and the top wall were kept at 300 K, but the lower wall was thermally insulated. For case 3, the electrodes and the lower wall were kept at 300 K and the upper wall was insulated, and finally, for case 4, the two walls, upper and lower, as well as the electrodes, were kept at 300 K. The used voltage was 10 V, and the average input velocity was 100 μm/s.

Figure 7 shows the temperature contours for the type 2 biosensor with the four various temperature control cases. It is noted that setting the temperature limits of Case 1 produces the greater temperature increase and gradient nearby the reaction surface than the other three cases and it is also noted, as shown in Fig. 8 and Table 7, that case 1 gives the best performance in terms of improving the detection time and the kinetics of the biosensor binding reaction compared to the other three cases of thermal adjustment.

Fig. 7

Temperature distribution for four thermal boundary conditions cases of the type 2 biosensor. The outlines are plotted on the x–y plane for z = 0 μm

Fig. 8

Time normalized surface concentration for type 2 biosensor with four distinct instances of temperature boundary conditions

Table 7

Detection time of type 2 biosensor under 10 V applied voltage for the four cases of thermal limits setting

Temperature boundary conditions	Detection time (s)
Case 1 (electrodes are at 300 K)	6950
Case 2 (top wall and electrodes are at 300 K)	7225
Case 3 (lower wall and electrodes are at 300 K)	9525
Case 4 (top wall, bottom wall, and electrodes are at 300 K)	9750

Electrical conductivity effect

From the expression of the electrothermal force, we can see that it depends on the electrical conductivity σ. As the biochemical solutions used as carrier fluid cover a wide range of electrical conductivity, we thought, here, study the effect of these different conductivities on the SARS-CoV-2 S protein/antibody binding reaction kinetics. The electrical conductivities used in this section are, respectively, σ = 0.01, 0.03, 0.0575, and 0.1 S / m. Type 2 biosensor model was used. The electrodes were kept at T0, and all the walls were insulated (Case1 thermal boundary condition). Figures 9 and 10 show the simulation findings of the binding reaction and the temperature jumps for the fluid with different electrical conductivities, respectively, when the applied voltage is 10 V. The results of the two figures show that the growth in temperature of the fluid caused by the Joule heating is nearly linear to the electrical conductivity but not for the reaction rate reinforced by the electrothermal effect. Consequently, care would be taken to use the high conductivity buffer solution as a means of decreasing the time detection biosensor for biological applications.

Fig. 9

Normalized surface concentrations versus time for type 2 biosensor with various electrical conductivities

Fig. 10

Temperature rise in the microfluidic channel of the type 2 biosensor as a function of the conductivity of the fluid with an applied voltage of 10 and 15 V

Conclusion

This work aimed to find new designs of biosensors for reliable, economical, and sensitive detection, especially for the early detection of infectious diseases such as COVID-19. We studied the reaction kinetics of the SARS-Cov-2 S protein binding under electrothermal flow by performing a 3D simulation on two heterogeneous immunoassay models having different shapes and geometric locations of the reaction surface and electrodes. The results show that the electrothermal force improves the binding efficiency in both types of biosensors, especially in the second type which has an annular surface reaction, where the enhancement was about 69% with an applied voltage of 10 V. Four cases of thermal boundary conditions are studied to show their influences on the fluid temperature field and on the efficiency of the binding reaction. Although increasing the electrical conductivity of the buffer solution can lead to a bigger electrothermal flux and aid to rise the binding reaction rate, care should be taken as this also rapidly increases the temperature of the buffer solution.