A computational simulation platform for designing real-time monitoring systems with application to COVID-19.

With the aim of contributing to the fight against the coronavirus disease 2019 (COVID-19), numerous strategies have been proposed. While developing an effective vaccine can take months up to years, detection of infected patients seems like one of the best ideas for controlling the situation. The role of biosensors in containing highly pathogenic viruses, saving lives and economy is evident. A new competitive numerical platform specifically for designing microfluidic-integrated biosensors is developed and presented in this work. Properties of the biosensor, sample, buffer fluid and even the microfluidic channel can be modified in this model. This feature provides the scientific community with the ability to design a specific biosensor for requested point-of-care (POC) applications. First, the validation of the presented numerical platform against experimental data and then results and discussion, highlighting the important role of the design parameters on the performance of the biosensor is presented. For the latter, the baseline case has been set on the previous studies on the biosensors suitable for SARS-CoV, which has the highest similarity to the 2019 nCoV. Subsequently, the effects of concentration of the targeted molecules in the sample, installation position and properties of the biosensor on its performance were investigated in 11 case studies. The presented numerical framework provides an insight into understanding of the virus reaction in the design process of the biosensor and enhances our preparation for any future outbreaks. Furthermore, the integration of biosensors with different devices for accelerating the process of defeating the pandemic is proposed.

Introduction

COVID-19

The current pandemic, the coronavirus disease is known as “COVID-19” (World Health Organization, 2020). According to the Coronaviridae Study Group (CSG) of the International committee on Taxonomy of viruses' studies, this virus is as forming a sister clade to the severe acute respiratory syndrome coronaviruses (“SARS-CoVs”). It is named as “SARS-CoV-2” by the CSG (Gorbalenya et al., 2020) (Previously named “2019 nCoV”). Fig. 1 (step 1) presents illustration of the structure of the coronavirus. It uses its spikes protein to enter the cell and delays the immune system response so that when immune system responses, the infection is progressed enough that fighting it would be hard (Guevara-Carrion et al., 2011). The three spill-over of animal to human coronaviruses are Middle East respiratory syndrome coronavirus (“MERS-CoV”), severe acute respiratory syndrome coronavirus (“SARS-CoV”) and 2019 nCoV (“SARS-CoV-2’).

Fig. 1

Schematic view of the four stages of designing a biosensor with the use of computational fluid dynamics. In stage 1, the virus is studied, and a conceptual design is prepared in stage 2. In stage 3, numerical analysis takes place and the design is enhanced in stage 4 based on the results of the simulation.

Since these viruses have the same origin, scientists are using the studies on MERS-CoV and SARS-CoV in order to reach to a faster solution for the SARS-CoV-2. SARS-CoV and SARS-CoV-2 both use angiotensin converting enzyme II (ACE2) as a cellular entry receptor and facilitate entry into the same panel of cell lines (WHO, 2020; Hoffmann et al., 2020; Seo et al., 2020). Antibodies, vaccines, and drug candidates are being developed throughout the world with the aim of preventing virus from entering the cell, viral replication or delaying the immune system.

The quantitative real-time polymerase chain reaction (qRT-PCR) test has been used broadly to detect SARS-CoV-2, SARS and MERS guides has been employed as a help for this test (WHO, 2020). Although the conventional technologies such as quantitative.

-real-time polymerase chain reaction (qRT-PCR) are cost, labour, and time consuming, they have a higher sensitivity and specificity, hence they are the gold standard for biosensor design (Choi, 2020).

Numerical simulations

Biosensors are sensitive, fast, low-cost, and easy-to-use pathogenic virus detection systems which are one of the best solutions for the reduction of death in pandemics (Amritsar et al., 2020). Microfluidic integrated biosensors have been designed and developed for these means, in which a microfluidic channel is linked to electrical chips as micro-scale laboratories for analyzing reaction results. In Fig. 1 (step 2), a schematic view of microfluidic channel equipped with a biosensor is illustrated. Targeted molecules flow into the channel and have chemical reactions with biological recognition element on the surface of the biosensor. Designing a reliable, efficient, affordable and quick biosensor for point-of-care (POC) testing still remains a challenge (Liu et al., 2020).

The most important part in designing a procedure is related to the transport of targeted molecules to the functionalized surface of the biosensor (Luka et al., 2015; Nair and Alam, 2006). Robust simulation (fast, reliable, and stable numerical modelling) of these systems and their chemical reactions would be a huge help in designing efficient biosensor. Computational fluid dynamics (CFD) is one of the best solutions for this matter. In recent years, CFD has successfully been applied to numerous biomedical-related projects involving design, validation and proof-of-concept; some examples from the present authors include (Deyranlou et al., 2020; McElroy and Keshmiri, 2018; Ruiz-Soler et al., 2017; Swanson et al., 2020).

In the previous works, different software packages and codes were utilized to model the flow inside the microfluidic channel. Discretization method is the key component in these numerical solutions. There are three different discretization methods; finite element method (FEM), finite volume method (FVM) 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 (Tombarevic et al., 2013). In some of them, the control-volume method is employed in the microfluidic flow simulation using software packages such as CFD-ACE+, Ansys-CFX, Ansys-Fluent and Flow3D, however, to avoid numerical diffusion (when algorithm aggregates numerical errors through the simulation) the actual diffusion constant from the experiments has not been used (Glatzel et al., 2008). Another common numerical study reported for this application is based on using an analytical velocity profile (assuming fully developed Poiseuille flow) in the convective-diffusive transport equation instead of using the converged velocity from Navier-stokes equations. Gelarkin finite-element method (Selmi et al., 2017) and COMSOL commercial code (Hu et al., 2007; Liu et al., 2019), (including FEMLAB, MATLAB subsidiary of COMSOL) are in this category, in which it has been reported that numerical results has more than 2 orders of magnitude difference with the experimental data (Squires et al., 2008).

Aims and objective

It is, therefore, the aim of the present work to introduce a reliable and stable numerical model that has been developed to solve full Navier-Stokes equations coupled with convection, diffusion and reaction of targeted molecules based on CVFEM. The model can simulate microfluidic-integrated biosensor without the need to modify the input values (properties of the real experimental samples) and the results will have good agreement with already existing experimental data from the literature.

The next section presents the description of the methods, governing equations, and boundary conditions. As was alluded to earlier, Fig. 1 demonstrates the four stages of designing a microfluidic-integrated biosensor with the proposed numerical platform. In step 1, the virus and its origins is studied, which is critical for step 2, for preparing a conceptual design of a microfluidic-integrated biosensor. In the next step, the numerical model simulates the buffer fluid flow, convection, diffusion and reaction of the targeted molecules in the biosensor. Finally, in step 4, design optimization takes place, based on the results and analysis. As a result, the present study will focus on studying the effects of varying the most design parameters on the performance of the designed biosensor.

Methods

In the proposed new approach for modelling virus detection in a microfluidic-integrated biosensor, full Navier-Stokes equations (Equation (1)–(3)) coupled with convection, diffusion and reaction of targeted molecules (Equation (4)–(6)) are solved implicitly with the CVFEM. Step 3 in Fig. 1 demonstrates the numerical domain and boundaries. Since variation of velocity in the span wise direction of this channel is negligible (due to negligible side wall effects on sensor surface), 2D simulation is valid for this case (Liu et al., 2019).

In the micro-scale, any parameter proportionate to the scale of the surface area becomes more important and viscous forces dominate over inertial forces. In order to avoid numerical diffusion or suppression of small numbers during numerical loops, high-order discretization methods are needed (Glatzel et al., 2008), which have been implemented in the proposed model. In addition, diffusion and reaction are coupled and solved fully implicitly using CVFEM (Karimian and Schneider, 1995). Furthermore, instead of using an analytical velocity profile inside the channel, full Navier-Stokes equations are solved for a comprehensive numerical model.

Governing equations

The proposed numerical platform solves the two-dimensional Navier-Stokes equations (Equations (1)–(3)) to simulate the buffer fluid flow. where u and v are the velocity in the x and y direction, respectively, while ρ is the density and μ is the molecular viscosity. The fluid is considered as continuum and in-compressible, because of its negligible changes in density. In the present test cases, the concentration of biological species is significantly low, and they do not have any effect on density or viscosity of the carrier fluid (Berthier and Silberzan, 2001). By using Fick's law (Equation (4)) and effect of velocity field a comprehensive equation for modelling kinetics of biological species can be generated (Equation (5)): where c is the concentration of the targeted molecules, D is the diffusion coefficient and S is the source/sink term which is assumed zero in this study. Kinetics of heterogeneous biosensors (which are used for DNA hybridization, virus detection etc.) is modelled using Langmuir-Hinshelwood mechanism (Berthier and Silberzan, 2001). Ligands constantly trap targeted bio-species and they dissociate at a smaller rate (Equation (6)):where k is the adsorption rate, k is the dissociation rate of the sensor, b is the surface concentration of bound analyte (Sevenler et al., 2019), c is the inlet concentration and b is the density of binding sites on the sensor. For discretizing the convection term, the upwind differencing scheme is used, while for the diffusion term, improved skewed upwind differencing (SUD) is used (Karimian and Schneider, 1995), which helps generate accurate and robust results. Upwind differencing scheme (UDS) is first order accurate and gives high level of false diffusion. SUD scheme uses flow direction which provides more realistic results. In the improved SUD scheme, effects of the diffusion and source term are included in the differencing procedure.

Boundary conditions

Boundary conditions implemented in this model are summarized in Table 1 and illustrated in Fig. 1. For walls, since there is no mass flux from the concentration field to or through them, the homogeneous Neumann condition is used. On the surface of the biosensor, where chemical reactions govern mass flux to the wall, the Neumann condition is used. The sequence of our developed model in this study is provided in Table 2 .

Table 1

Boundary condition – velocity and concentration for walls, sensor, inlet and outlet of the channel.

Type	Velocity (u)	Concentration (c)
Interior	Navier-Stokes equations	Convection-diffusion-reaction
Walls	No slip	Homogeneous Neumann ( ∂c∂n=0)
Sensor	No slip	Neumann ( ∂c∂n=−1D∂b∂t)
Inlet	u = u₀	c = c₀
Outlet	Zero gradient	n→.(D∇c)=0

Table 2

Graphical representation of the Algorithm behind the numerical model proposed in the present work.

Algorithm 1: Numerical model developed in this study.
	Data: geometry and grid resolution
	Input: initial condition and boundary condition (Table 1)
	Result: same for output data
1	initialization;
2	Whilet < tmaxdo
3		set Δt;
4		set calculation matrices to zero;
5		forallCVFEM nodes/volumesdo
6			calculate upwind vectors;
7			calculate center points and distance to edges for the improved SUD scheme;
8			forallsub elementsdo
9				generate coefficients of the conservation of mass (Equation (1));
10				generate coefficients of transient, convection, diffusion, pressure and source term of the momentum equation (Equations (2), (3)));
11				generate coefficients of convection, diffusion and chemical reaction (Equations (4), (5), (6));
12			end
13			build the global matrix based on the coefficients;
14		end
15		apply the velocity and pressure boundary conditions of the buffer fluid (Table 1) on the global matrix;
16		apply the boundary conditions of the sensor (Table 1) on the global matrix;
17		generate the band matrix of the global matrix and solve it for velocity, pressure, and concentration;
18		foralldomaindo
19			assign the results to the variables and define old values;
20		end
21		t←t + Δt;
22	end

Results

In this study, the proposed model is initially validated against the existing data from the literature (Berthier and Silberzan, 2001). The validation is conducted for a case where the buffer-fluid flows through a microfluidic channel (10−3 m in height and 10−2 m in width) with a flow rate of 10−6 m 3 /s, diffusion ratio of 7 × 10−11 m 2 /s and an input concentration of targeted molecules of 2.5 × 10−6 Mol/m 3. The density of binding sites, adsorption and dissociation rate of the sensor are 1.668 × 10−8 Mol/m 2, 75 m3 /Mol·s and 10−2 1/s, respectively. Fig. 2 demonstrates that the results of the current numerical study is in excellent agreement with the experimental data.

Fig. 2

Validation of the current numerical model with the experimental results (Berthier and Silberzan, 2001), the normalized surface concentration over time. The surface concentration is normalized to the maximum density of binding sites on the sensor.

Following a successful validation exercise, a base setup, suitable for coronavirus detection, is chosen for studying the effects of three functional parameters on the performance of the microfluidic-integrated biosensor. The previous study by Qi et al. (2006) shows that the affinity (k  = k /k ) of antibodies used in biosensors for SARS-CoV (b1 and h12) is in the order of 10−6 Mol/m 3. Specifically, the order of magnitude of adsorption rate is 103 m3 /Mol·s and that for dissociation rate is 10−3 1/s. In the numerical setup, sample flows into microfluidic channel with input concentration (c ), Diffusion constant (D) and Flow rate (Q) of 1 × 10−11 Mol/m 3, 1 × 10−11 m 2 /s and 1.7 × 10−10 m 3 /s, respectively. The density of binding on the sensor (b ) is 3.3 × 10−8 Mol/m and it is installed at the middle of the wall of the microfluidic channel, as Fig. 1 (step 3) presents. The length of the microfluidic channel (L) is 10-3 m and its height is 10-4 m, suitable for fabrication of biosensors especially designed for SARS-CoV-2 (Seo et al., 2020).

Then for enhancing the performance of the biosensor and studying the effect of important parameters, 11 cases has been studied. Details of these numerical simulations is demonstrated, and color scaled in Table 3 . In the first study (group A) the installation position of the functionalized surface is changed. Groups B and C present results of the study on affinity of the biosensor and concentration of the targeted molecules, respectively. The center position of the biosensor for the first case is located 312.5 μm away from the inlet. The subsequent cases are moved at the intervals of 125 μm towards the outlet. Fig. 3 (a) provides the binding cycle for these five cases and their saturation time. In stage 1 of the binding cycle, targeted molecules flow into the microfluidic channel equipped with functionalized zone on the surface of the wall with ligands, which will have affinity with targeted molecules. At the start of stage 2, targeted molecules are injected into the buffer fluid and they start to be adsorbed and immobilized temporarily by the available sites of the sensor. In this stage, the detection of surface concentration will begin. When concentration reaches its asymptotic value (i.e. stage 3), an analyte free fluid passes through the channel and the desorption of ligands takes place (i.e. stage 4).

Table 3

Details of the numerical simulations. Group A (cases 1 to 5), which is for studying the effects of the position of the biosensor and case 3 is the base case. Group B (cases 6 to 8) for studying effect of varying dissociation and adsorption rates (but keeping the affinity the same) of the biosensor and case 7 is the base case. Group C (cases 9 to 11) for studying the effect of concentration of the targeted molecules in the sample and case 10 is the base case in this group.

Case	xs [μm]	kon [m³/Mol.s]	koff [1/s]	c0 × 10⁻11 [Mol/m³ ]
A-1	312.5	1000	0.001	1
A-2	437.5	1000	0.001	1
A-3	562.5	1000	0.001	1
A-4	687.5	1000	0.001	1
A-5	812.5	1000	0.001	1
B-6	562.5	100	0.0001	1
B-7	562.5	1000	0.001	1
B-8	562.5	10000	0.01	1
C-9	562.5	1000	0.001	1
C-10	562.5	1000	0.001	100
C-11	562.5	1000	0.001	1000000

Fig. 3

(a) Binding cycle for different installation position, the normalized surface concentration () over time. The surface concentration is normalized to the maximum density of binding sites on the sensor. The other properties of these five cases (cases 1 to 5) are constant and similar to the base case. (b) Variation of the saturation time with sensor position for cases 1 to 5. (c) Binding cycle of three sensors with the same affinities but different magnitudes of dissociation and adsorption rate (cases 6 to 8). (d) Binding cycle of a sensor for different inlet concentrations (cases 9 to 11). (e) Graph of variation of the saturation time with sensor position, adsorption rate, and dissociation rate and inlet concentration. Cases 1 to 11 of Table 3. Values presented in this graph are normalized.

For a better analysis of the changes in the performance of the biosensor, saturation time versus the position of the sensor is illustrated in Fig. 3(b). As the biosensor transfers from the inlet, its saturation time increases As Fig. 3(b) shows, the variation of the saturation time with respect to the location of the sensor reduces as the biosensor gets far from the inlet or outlet. Displacements happen at the same rate, although changes in binding cycle and saturation times are not at the same rate. As Equation (7) presents, the relation between the saturation time (t ) and the position of the sensor from inlet plane (x ) is non-linear. Equation (7) presents the best-fit polynomial line through the numerical results:where t represents the saturation time and x is the position of the sensor from the inlet boundary.

Next, we study the effect of functionalized surface with the same affinities but different magnitudes of dissociation and adsorption rate on the total performance of the biosensor. Kinetics in heterogeneous reactions depend on reaction rates and diffusion of species towards the functionalized surface. These values can change by modifying the ligands used on the functionalized surface that effects the efficiency of the whole biosensor. Fig. 3(c) provides the results for this study. As reaction rates increase, the total time over which hybridization reaches its saturation value (stage 3) decreases. In addition, these time differences show a significant increase in cycle time when adsorption and dissociation rates reach orders of 102 m3/Mol·s and 10−4 1/s, respectively (case 6). For instance, saturation time would increase 40% as both adsorption and dissociation rates are reduced by a factor of 10.

Fig. 3(d) demonstrates the binding cycle of cases 9 to 11 for different concentration of targeted molecules in the sample. As the input concentration increases, the saturation time decreases. Changes in the saturation time is non-linear. Moreover, the shape of association period (stage 4) changes significantly for different concentrations of the bio-species. The concentration of antibody is critically important for the biosensor performance. Various concentration of antibody is tested here in order to achieve the best performance for the biosensor.

Overall, as Fig. 3(e) demonstrates, the three design parameters tested here have shown to have significant impact on the saturation time of the biosensor. In this design process, the saturation time can reduce from several hours to approx. one hour (in case 11) with only a small change in the imported sample. It is important to note the fact that the saturation time is needed for measuring the amount of infectious virus in the sample. For detecting the virus in the sample, it is not necessary to reach the saturation point.

The fundamental rule of biosensors in fighting against pandemic is proven in previous studies (Lee et al., 2020; Morales-Narváez and Dincer, 2020; Zhu et al., 2020). As it is shown above, design and optimization of an effective, low cost and reliable detection system in a short period of time aimed for fighting pandemics can be possible with an accurate numerical simulation. This platform is specifically designed for microfluidic integrated biosensors, since these devices reduce the amount of reagent, energy consumption, waste and cost due to their nature; they are an integration of chemical and bio-logical process in a single platform and can sense small volumes of analytes. These results prove the previous statement regarding the better accuracy of this numerical model in comparison to the other previous numerical studies using other software packages and codes (provided in the introduction section).

The above analysis elucidates the design controls imposed by convection, diffusion and reaction. First to determine the optimal position of the biosensor in microfluidic channel for better approach of the targeted molecules to the binding surface. Then to enhance the binding reactions, affinity of the biosensors are studied. At last, concentration of the targeted bio-species that is an important factor in the improvement of the mass transport. The performance of the biosensor can be improved based on an appropriate choice of these three design parameters.

Conclusion

Time is of paramount importance in the fight against a pandemic. The crucial role of biosensors in containing highly pathogenic viruses in saving lives and the economy is evident. The challenge is to design a biosensor for each specific for point-of-care (POC) applications.

To address this challenge, a new numerical solver for competitive biosensors is developed and presented in this work. The model is based on the control-volume finite-element method (CVFEM) and all the convective-diffusive-Langmuir equations are solved fully coupled together with a high order discretization scheme.

The first part of the study presented the validation test of the preliminary results against the experimental data. Unlike previously published numerical studies, the present model makes use of the exact flow parameters, adsorption/dissociation rate and diffusion constant, that shows the reliability and robustness of the developed approach.

The next part of the study showed the effect of 3 design parameters on the performance of the biosensor. It was found that the installation position of the biosensor has a significant effect on its efficiency; it can reduce the detection time by over 50% with only 500-μm displacement. It was also found that the reaction rates would lead to a considerable increase in the total time of the binding cycle when they reach a specific order of magnitude, thus it would be better to avoid using such ligands in reaction surfaces. Furthermore, changing the inlet concentration was also shown to have a significant impact on the saturation time and shape of binding cycle stages (which demonstrates the behavior of targeted molecules through time). These variations take place in a non-linear manner.

Since May 2020 until the time of writing this paper (July 2020), more than 20 diagnostic tests have been approved by FDA by receiving the emergency use authorization (EUA) (BioCentury, 2020). The sample-to-answer tests are chip-based biosensors (nucleic acid test) and paper-based (antibody test) biosensors. Because of complicated procedure of nucleic acid tests in comparison to the antibody test, recent studies show that antibody tests are playing an important role in reducing the overall diagnostic time and help in a fast decision making. (Choi, 2020; Du et al., 2020; Li et al., 2020; Sheridan, 2020).

Future work

The presented computational model would provide a valuable insight into understanding of the virus reaction in the designed biosensor in an everchanging inlet conditions and biosensor itself. The unique platform presented here would enhance our readiness for the current and future outbreaks and would help to develop innovative designs. The integration of biosensors with different devices would help in reducing the risk of the potential next waves of COVID-19. It would also help in containing any other virus-related disasters in future.

The current tests have the potential to produce wrong negative results and they have different sensitivity (BioCentury, 2020). Most of the FDA approved tests are qRT-PCR, which has lower limit of detection (LoD), hence higher level of sensitivity. Limit of detection (LoD) is one of the crucial performance characteristics that describes the lowest concentration of the target that can be reliably measured by the test (Lavín et al., 2018). This numerical platform can in understanding operational behavior of designed biosensors and reducing their limit of detection.

CRediT authorship contribution statement

Fatemeh Shahbazi: Software, Methodology, Validation, Investigation, Visualization, Conceptualization, Writing - original draft. Masoud Jabbari: Supervision, Conceptualization, Methodology. Mohammad Nasr Esfahani: Supervision, Writing - review & editing. Amir Keshmiri: Conceptualization, Resources, Writing - review & editing.

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.