COMSOL 6.4 - Modeling an LFA Rapid Detection Test

The COVID-19 pandemic has made the world painfully aware of the importance of contact tracing. Knowing who is infected, as well as when and where that person was infected, is required in order to be able to contain an outbreak of a virus infection. The lack of this knowledge has forced most of the countries in the world to implement restrictions and lockdowns, which have led to severe social disruptions and a global economic recession.

The development of rapid detection test devices has allowed for contact tracing in the areas of society that are critical and cannot shut down, such as healthcare and food supply. Access to inexpensive tests has also allowed for people to self-diagnose and then self-isolate, if infected. This has also helped in slowing down the spreading of the virus. South Korea and Germany have successfully used testing to slow down the pandemic on a broader scale, in all areas of society (Ref. 1).

One possible rapid detection test for COVID-19 is based on lateral flow assay (LFA), also called lateral flow immunoassay (LFIA) or immunochromatographic tests. The standard pregnancy test that can be bought at supermarkets is probably the most familiar application of LFA. LFA offers an inexpensive, relatively reliable, robust, and easy-to-use test for many diseases and conditions (Ref. 2).

The test consists of the following parts (shown in Figure 1): The sample pad, the conjugation pad, the membrane, and the absorption pad, also called a wick pad. The assembly of the different components (pads and membrane) is often referred to as the test strip. The test strip is protected by a plastic housing (Ref. 3).

Figure 1: Schematic drawing of a rapid detection test device. Here, the thickness of the membranes and pads is exaggerated by a factor of 5 in order to show the structure. The assembly of the pads and the membrane can be referred to as a test strip.

Test Functionality

When the body is infected with the coronavirus, SARS-CoV-2, the immune system quickly reacts by forming antibodies. Dendritic cells may present virus antigens to allow recognition by T-cells. T-cells may activate B-cells to secrete antibodies that target the antigen (Ref. 1). Among the first to form are the so-called IgM antibodies. These antibodies attach to the surface of so-called antigens on the virus particles as soon as they come close. In the case of the coronavirus, the antigens can be the spike proteins at the surface of the virus (S antigen). Once attached to the antigens, the antibodies block the spike proteins of the virus, impeding them to attach and infect human cells. This neutralizes the virus, since it is not able to replicate outside of an infected cell. There are a number of different antibodies that may target different antigens. Note that there are also other mechanisms for fighting the infection. In addition, T-cells that recognize the virus may also target infected cells directly. They can instruct the cells to self-destruct or they may kill the infected cells, thus also neutralizing the virus.

The IgM antibodies manufactured by the immune system attach, for example, to the spike antigen (S antigen) of the SARS-CoV-2 virus particle, thus neutralizing the particle. The neutralized virus cannot enter a human cell and therefore cannot replicate. The virus particle is destroyed.

The IgM antibodies patrol the body in groups of five, forming small particles (or large molecules), attaching to every virus particle they encounter. Later in the infection, the immune system also forms other antibodies, for example, IgG, that patrol the body on their own and fasten to every virus particle they see. The IgG antibodies take longer to manufacture by the body, but they also last longer and grant immunity as long as they are present.

Some of the LFA rapid detection tests for COVID-19 are based on the detection of both IgM and IgG antibodies. These are the tests that are investigated in the model featured here.

Figure 2: Schematic drawing of a rapid detection test device. The sample contains the human antibodies IgM and IgG, created due to COVID-19 infection. An animal antibody (AA) is also added to the sample liquid with the buffer solution. The test lines have immobilized antibody detectors on the three different areas. Note that the test lines are not yet visible.

The patient’s blood (or saliva) may be applied to the sample well followed by an application of the buffer, by dripping a couple of droplets of buffer in the sample well. The buffer contains an animal antibody AA to test whether the test is working correctly.

The sample is transported by the capillary forces to the conjugate pad. Here, the IgM and IgG antibodies attach to a conjugate label. A conjugate label may be a gold nanoparticle that has SARS-CoV-2 antigens on its surface. Two different complexes are then formed: IgM with the conjugate (=IgM-C) and IgG with the conjugate (=IgG-C). Also, the animal virus antibody picks up its respective particle (AA-C). These complexes are dissolved in the sample liquid.

The sample is then transported to the membrane by the capillary forces. In the first test line, IgM detectors are immobilized at the surface of the membrane during the manufacturing of the membrane. These IgM detectors capture the IgM-C complex and keep these complexes stationary in the area of this test line. The nanoparticles accumulate and reveal the presence of the complex on the test line by coloring the test line red.

In the same way, the IgG-C complex and the AA-C complex react with immobilized detectors in the other test lines. The coloring of the control test line reveals that the sample has passed the membrane area.

The flow continues to the absorption pad (wick pad). The pore volume of this pad determines the sample volume that can flow through the test strip. Once the absorption path is full, the flow in the test strip stops. The only possibility to restart the flow is to evaporate some of the sample liquid in the absorption pad.

Model Definition

It could be shown that once passed the sample pad, the sample fluid flow quickly forms a flat velocity profile (Ref. 4). This means that it flows uniformly along the width of the test strip. This implies that a 2D model is enough to understand the challenges and functionality of the rapid detection test device, as long as the sample pad does its job in distributing the flow uniformly. We therefore used 2D models to investigate the transport and reactions in the test strip. The 2D models allow us to use a higher resolution of the mesh along both the length and thickness of the test strip.

The model combines the Richards’ Equation interface for porous media flow and the Transport of Diluted Species in Porous Medium interface in COMSOL Multiphysics.

Here, the model was simplified assuming that the reaction sites are spread all across the thickness of the conjugate pad and, for the test line, across the membrane thickness. Therefore, the reactants do not have to be transported to the surface of the test strip and the reaction between the conjugate label and the antibodies are maximized. In this way, the conjugate label picks up as much antibodies as possible.

The reactions that form the IgM-C, IgG-C, and AA-C complex species are defined by the Domain ODEs and DAEs interface.

Fluid Flow

Richards’ equation describes the unsaturated–saturated flow of water in porous media. In this problem, Richards’ equation can only be used for the water since the air in the test strip is supposed to be at atmospheric pressure. The governing equation for the model is

(1)

The pressure, p, is the dependent variable, C denotes the specific moisture capacity (m−1), Se is the effective saturation, Sp is a storage coefficient (m−1), t is time, κs denotes the permeability at saturation (m2), and κr is the relative permeability, which is calculated, after Brooks and Corey, as

(2)

The material properties are assumed to be constant within the whole domain, regardless of the different pads or membranes used in the test.

Chemical Reactions

On the conjugate pad, two processes take place. First, IgM and IgG antibodies from the saliva sample bind to gold nanoparticles SCoAu, creating complexes IgMC and IgGC.

Second, animal antibodies are released into the liquid.

On the first test line, the IgMC complex reacts with the adsorbed IgMd detector protein to form the adsorbed IgMPos surface species. The formation of the immobilized IgMPos gives the first test line its detection color.

Analogously to above, in the second test line the IgGC complex reacts with the adsorbed IgGd detector protein. The adsorbed species IgGPos gives the second test line its red color.

On the third test the antigen complex AAC binds to immobile detector molecules AAd. The adsorbed species AAPos gives the control test line its red color.

In this homogenized version, the chemical reactions are provided in the Domain ODEs and DAEs interface which can be imported from an external file using the option.

Transport of Antibodies

The species transport is simulated using the Transport of Diluted Species in Porous Media interface. The values used for diffusion and reaction rates are listed in Table 1.

Table 1: Transport properties.
Name	Value	Description
D1	3·10-11 m2/s	Diffusion coefficient within all domains
k1	10 m3/(s·mol)	Rate constant, conjugate pad
k2	1 m3/(s·mol)	Rate constant, test line
c0	2.3·10-4 mol/m3	Inlet antibody concentration
c0ads	5·10-7 mol/m2	Adsorbed species surface concentration,conjugate pad

Results and Discussion

Figure 3 and Figure 4 show how the pressure front and the effective saturation front proceed with time. At about 260 s of simulated time the test strip is fully saturated.

Figure 3: Pressure distribution at different output times (t = 20 s, t = 60 s, t = 210 s, t = 280 s, and t = 410 s).

Figure 4: Effective saturation in the test strip at 20 s, 60 s, 210 s, 280 s, and 410 s (from top to bottom).

Figure 5: Concentration field of IgMC after 20 s (top), 60 s, 210 s, 280 s, and 410 s (bottom).

Figure 5, Figure 6, and Figure 7 illustrate the concentration fields of IgM-C, IgG-C, and AA-C complex species at different simulation times. After 410 s, there is no longer any flow through the test strip. We get a zone with high concentrations in the conjugate pad and a zone with low concentrations below the test line for the corresponding antibody.

Figure 6: Concentration field of IgGC after 20 s (top), 60 s, 210 s, 280 s, and 410 s (bottom).

Figure 7: Concentration field of AAC after 20 s (top), 60 s, 210 s, 280 s, and 410 s (bottom).

The impact of the flow is clear if we plot the particle concentration of the detection species IgMPos, IgGPos, and AAPos at their respective test line area (Figure 8). The increase in particle concentrations starts with a delay given by when flow reaches the respective test line. The concentration grows almost linearly as IgMC is adsorbed at the test line area and the detector species IgMPos is formed. IgMPos and IgGPos compete for the same conjugate label and therefore AAPos grows slightly faster. The linear increase means a constant adsorption rate. When the flow stops, due to the saturation of the absorption pad, the rate of formation of IgMPos continues about the same rate. This implies that, in this case, the formation of IgMPos is determined by the adsorption kinetics; that is, it is determined by the adsorption rate of the gold nanoparticles in IgMC to the adsorption site. If it was mass transport controlled, we would see a change of the slope of the curve and the grow would be slowed down when the flow stops. This could, of course, change if we change the rate constant for the adsorption kinetics. The visibility of the test lines starts at about 1 × 108 particles/mm2.

Figure 8: Particles per unit surface (t = 20 s, t = 60 s, t = 210 s, t = 280 s, and t = 410 s).

References

1. T. Kilic, R. Weissleder, and H. Lee, “Molecular and Immunological Diagnostic Tests of COVID-19: Current Status and Challenges,” iScience, vol. 23, no. 8, 101406, 2020; doi.org/10.1016/j.isci.2020.101406.

2. B.G. Andryukov, “Six decades of lateral flow immunoassay: from determining metabolic markers to diagnosing COVID-19,” AIMS Microbiology, vol. 6, no. 3, pp. 280–304, 2020; pubmed.ncbi.nlm.nih.gov/33134745/.

3. “Rapid Lateral Flow Test Strips, Considerations for Product Development,” Merck Millipore, 2013 EMD Millipore Corporation, Billerica, MA, USA.

4. www.comsol.com/blogs/modeling-a-rapid-detection-test-in-comsol-multiphysics/

Porous_Media_Flow_Module/Solute_Transport/rapid_detection_test_2d

Modeling Instructions

Start setting up the physics for flow in unsaturated porous media using the Richards’ Equation interface. The physics for the chemical reactions and the transport mechanisms are added later.

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Import the geometry sequence from a separate file.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file rapid_detection_test_2d_geom_sequence.mph.

In the toolbar, click

Click the

button in the toolbar.

Definitions

As the geometry is long and very flat, the view is scaled in the y direction to make it easier to select different domains within the geometry.

In the window, expand the > node.

Axis

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the text field, type 10.

Click

Click the

button in the toolbar.

Global Definitions

The geometry has been set up using parameters. Add more parameters relevant for the material, chemistry and operation conditions. For convenience, all these parameters can be imported from external files.

Parameters: Geometry

In the window, under click .

In the window for , type Parameters: Geometry in the text field.

Parameters: Material

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file rapid_detection_test_2d_parameters_material.txt.

In the text field, type Parameters: Material.

Parameters: Operating Conditions

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file rapid_detection_test_2d_parameters_operating_conditions.txt.

In the text field, type Parameters: Operating Conditions.

Richards’ Equation (dl)

Next, set up the physics for unsaturated porous media flow.

In the window, under click .

In the window for , locate the section.

Clear the checkbox as the geometry is so flat that gravity does not take any effect.

To improve the convergence of the model, set the shape function to linear.

Click to expand the section. From the list, choose .

Unsaturated Porous Medium 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the dp text field, type d_p.

Locate the section. From the list, choose .

In the α text field, type alpha.

In the l text field, type 1.

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the p text field, type p0.

Materials

Now add the materials needed.

Fluid

In the window, under right-click and choose .

In the window for , type Fluid in the text field.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho_l	kg/m³	Basic
Dynamic viscosity	mu	mu_l	Pa·s	Basic

Porous Material 1 (pmat1)

Right-click and choose > .

In the window for , locate the section.

From the list, choose .

Locate the section. In the εp text field, type por.

Locate the section. Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Global Definitions

A pressure boundary condition is prescribed at the inlet. To improve stability, the inlet pressure is ramped up smoothly over a short initial time interval rather than being applied instantaneously. Therefore, create a step function first.

Step 1 (step1)

In the toolbar, click

and choose > .

In the window for , locate the section.

In the text field, type 0.5[s].

Click to expand the section. In the text field, type 1.

Richards’ Equation (dl)

Pressure 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

In the p0 text field, type p0+(pw-p0)*step1(t).

Mesh 1

To capture the saturation front correctly, a high mesh resolution is needed. In this case, a mapped mesh is chosen as this combines efficiency and good resolution.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Mapped 1

In the toolbar, click

Size

In the window, click .

In the window for , click to expand the section.

In the text field, type 1e-4.

Mapped 1

In the window, click .

In the window for , click to expand the section.

Select the checkbox.

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 1 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 3e-5.

Distribution 1

Right-click and choose .

Select Boundary 66 only.

In the window for , locate the section.

In the text field, type 3.

Distribution 2

Right-click and choose .

Select Boundaries 54 and 56 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 200.

In the text field, type 5.

Select the checkbox.

This expression creates a gradually increasing mesh size, so the elements become larger toward the right end.

Distribution 3

Right-click and choose .

In the window for , locate the section.

Click

In the dialog, type 6, 8, 25, 27, 53, 55 in the text field.

Click .

In the window for , locate the section.

In the text field, type 20.

Click

Zoom in and make sure the mesh size corresponds to that in the image below.

Definitions (comp1)

Next, set up the study. As the wetting front advances and the strip approaches full saturation, capillary suction decreases and no longer provides a stabilizing gradient, which may cause the solver to have convergence issues. To prevent this, use an that reinitializes the pressure to be in equilibrium with the applied pressure when the strip is almost fully saturated, i.e., when the effective saturation exceeds 0.8 at position 0.105 m. After that, the pressure remains constant and the Darcy velocity becomes zero.

Domain Point Probe 1

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to 0.105.

In row , set to th/2.

Point Probe Expression 1 (ppb1)

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type dl.Se.

Add Physics

In the toolbar, click

and choose .

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Events (ev)

Indicator States 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	g(v,vt,vtt,t)	Initial value (u0)	Description
saturated	ppb1-0.8	0

Implicit Event 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type saturated>0 as the event is triggered when the condition changes sign.

Locate the section. In the table, enter the following settings:


	Expression
p	pw

Locate the section. Clear the checkbox.

Study 1

This is a highly nonlinear problem, so solver settings need some adjustments. At the start, changes happen quickly, so choose time steps that capture the moving liquid front accurately. Use a stable damped Newton method to manage the nonlinearity. While convergence is the top priority, settings should also keep computation time reasonable. For more details of each setting, see the manual.

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type 5*(exp(range(0,1/19,1)*4)-1)/(exp(4)-1) range(10,5,420).

The first five seconds are divided into 20 time steps, which gradually increase in size. Keep in mind that the software uses extra time steps if needed for convergence.

Click to expand the section. From the list, choose .

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

Find the subsection. From the list, choose .

From the list, choose .

In the window, expand the > > > node, then click .

In the window for , click to expand the section.

In the text field, type 0.6.

In the text field, type 30.

In the window, under > > > click .

In the window for , locate the section.

From the list, choose .

In the window, click .

In the window for , locate the section.

Clear the checkbox.

In the toolbar, click

Results

A pressure plot and a plot of the effective saturation can be loaded from the menu. Follow the steps below to reproduce Figure 3 and Figure 4.

Result Templates

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > and > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Results

Pressure (dl)

In the window, expand the > node, then click .

In the window for , click to expand the section.

From the list, choose .

Click to expand the section. From the list, choose .

From the list, choose .

In the text field, type -1.25e-3.

Solution Array 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type 20 60 200 420.

In the toolbar, click

Click the

button in the toolbar.

Pressure (dl)

In the toolbar, click

and choose .

Table Annotation 1

In the window for , locate the section.

From the list, choose .

In the table, enter the following settings:


x-coordinate	y-coordinate	Annotation
0	0	20s
0	-1.25e-3	60s
0	-2.5e-3	200s
0	-3.75e-3	420s

Locate the section. Clear the checkbox.

Pressure (dl)

In the window, click .

In the toolbar, click

Effective Saturation (dl)

In the window, expand the > node, then click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

In the text field, type -1.25e-3.

Solution Array 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type 20 60 200 420.

Table Annotation 1

In the window, under > right-click and choose .

Table Annotation 1

In the window, right-click and choose .

In the toolbar, click

Click the

button in the toolbar.

Component 1 (comp1)

Next, set up the chemical reactions together with mass transport. If available in your license, the interface provides a convenient way to define reactions. In this model, the interface is used. For convenience, a fully set-up interface is imported from an external file.

In the toolbar, click and choose .

In the dialog, click

Browse to the model’s Application Libraries folder and double-click the file rapid_detection_test_2d_geom_sequence.mph.

Click .

Inspect the node (the Domain ODEs and DAEs interface). Eight species and hence eight dependent variables are defined. For each section of the detection test, a domain feature specifies the reactions described in Chemical Reactions. Some variables are still undefined at this stage, namely the species that are transported. These will be introduced in the next step, when setting up the transport interface.

Add Physics

In the toolbar, click