Team:Vilnius-Lithuania/Model

Introduction

Consider a permeable strip with length X greater than width Y. On one end of the strip we have a sample pad onto which we drop analytes A that are the oligos we want to test. Beneath the sample pad is the conjugate pad that holds gold nanoparticles P. Analytes and nanoparticles diffuse and react to form a complex C. We also assume that this end of the strip is dipped into a volume (or "bath") of water that starts permeating into the strip with flow velocity $V(t)$. Somewhere along the strip at points $L \leq x \leq U$ is a test line with fixed probes R that are specific to the oligo sequence we want to test. Probes R react with the complex C to form the signal S. \begin{align} \begin{split} \ce{A + P <=>[\ce{k_1}][\ce{k_{-1}}] C}, \\ \\ \ce{C + R <=>[\ce{k_2}][\ce{k_{-2}}] S}. \end{split} \end{align} The goal. The goal is to provide a mathematical basis for a fast and open-source app that helps in the design of a lateral flow assay test.

Advection-diffusion-reaction model of the lateral flow assay

Corresponding to the description we form the model: \begin{equation} \max_{L \in [0, X]} J(L) = S(t_{end}) \end{equation} subject to the advection-diffusion-reaction (ADR) system1 \begin{equation} \begin{cases} \frac{\partial [A]}{\partial t} = D_A \Delta [A] -V(t) \cdot \nabla [A] - k_1[A][P] + k_{-1}[C],\\ \frac{\partial [P]}{\partial t} = D_P \Delta [P] -V(t) \cdot \nabla [P] - k_1[A][P] + k_{-1}[C],\\ \frac{\partial [C]}{\partial t} = D_{C} \Delta [C] -V(t) \cdot \nabla [C] + k_1[A][P] - k_{-1}[C],\\ \frac{\partial [R]}{\partial t} = - k_2[R][C] + k_{-2}[S],\\ \frac{\partial [S]}{\partial t} = + k_2[R][C] - k_{-2}[S], \end{cases} \end{equation} with the initial conditions $A_0(\vec{x}) = P_0(\vec{x}) = C_0(\vec{x}) = S_0(\vec{x}) = 0$ on the strip and \begin{align} R_0(\vec{x}) = \begin{cases} R_0, \ L \leq x \leq U, \\ 0, \ \text{elsewhere}. \end{cases} \end{align} We consider the case where the inflow boundary satisfies \begin{equation} \begin{cases} [A](t,x=0) = A_0 H(T-t), \\ [P](t,x=0) = P_0 H(T-t), \\ [C](t,x=0) = 0, \\ \end{cases} \end{equation} $H(t)$ is the Heaviside function and $A_0, P_0$ are some positive constants.

These boundary conditions correspond to the case where the sample and conjugate pads are partially dipped into a water bath and that there has not been substantial mixing of analytes A and gold nanoparticles P in the water bath before time $t=0$. Time $t=0$ corresponds to the moment the sample and conjugate pads are saturated and the solution starts moving laterally into the rest of the strip through the boundary x=0 for a period of time $T$.

The non-mixing assumption is relevant to the optimization because we may consider the lower extreme case of the initial complex concentration since it monotonically increases with time. This allows us to mathematically ignore the initial diffusion of a drop of analytes through the sample and conjugate pads.

Everywhere else we assume that the particles cannot move through the boundaries, corresponding to the Neumann conditions $\nabla[A] \cdot \vec{n} = 0, \nabla[P] \cdot \vec{n} = 0, \nabla[C] \cdot \vec{n} = 0, \nabla[R] \cdot \vec{n} = 0, \nabla[S] \cdot \vec{n} = 0$.

Time T can be estimated from the conservation of mass. \begin{equation} n = \int_0^T c A |\vec{V}(t)| \, dt, \end{equation} n is the molar mass of a substance, c is the concentration, A is the cross-sectional area of the inflow boundary, $|\vec{V}(t)|$ is the speed of the fluid flow through the boundary. The solution for time $T$ is \begin{equation} T = x^{-1}(\frac{n}{cA}). \end{equation} Initial concentrations. If we drop a volume $V_1$ of an analyte with concentration $c_1$ onto a sample or conjugate pad, then the initial concentrations of the analytes, $a_0$, can also be calculated by the conservation of mass: \begin{equation} c_1 V_1 = c_2 V_2, \end{equation} where $V_2$ is the volume of a water bath, and $c_2$ is the resulting initial concentration.
The initial concentration of gold nano-particles $p_0$ is calculated similarly, except that the volume $V_2$ is taken with respect to the volume of the conjugate pad.
Velocity. To calculate the velocity in the flow direction, we propose several simple models for the wetting front velocity that were found in the literature2: \begin{equation} V(t) = C, \ V(t) = Ce^{-rt}, \ V(x) = \frac{C}{x}, \end{equation} corresponding to different wetting front profiles: \begin{equation} x(t) = Ct, \ x(t) = \frac{C}{r}(1 - e^{-rt}), \ x(t) = \sqrt{2Ct}. \end{equation} Velocity estimation. Time-series data of a video recording of lateral flow experiments was used to fit the models. In the experiments we put two groups of 8 flow assays next to a ruler under a video camera and added different volumes of buffer solution. Dependence on y coordinate along the width of the strip was masked by random noise, and thus we neglected it.

The experiment showed clearly that the velocity is non-constant, and for 60-100$\mu L$ range the profile was nearly identical, whereas for volumes larger than 100$\mu L$ the velocities started increasing. The exponential velocity profile had the closest fit. The radical profile that was proposed in the literature\sup{berli_quantitative_2016} appeared to be inaccurate.

At volumes less than 60 $\mu L$ the solution never left the sample pad, hence there was no flow through the boundary at $x=0$. At 60 $\mu L$ the solution never reached the control line at the other end of the strip. At initial volumes larger than 125 $\mu L$ the solute started overflowing. Thus, in addition, we conclude that reasonable analyte volumes are between 60 and 125 $\mu L$.

Diffusion coefficient. The Stokes-Einstein equation \begin{equation} D = \frac{k_B T}{2 \pi \eta r} \end{equation} is used to calculate the effective diffusion coefficient \begin{equation} D_e = D\frac{\epsilon_t \delta}{\tau} \end{equation} in a porous medium.
$k_B$ is Boltzmann's constant, T is absolute temperature, $\eta$ is the viscosity of fluid in a medium, and r is the radius of a particle. $\epsilon_t$ is the porosity, $\delta$ is the restrictivity, and $\tau$ is the tortuosity of the medium.

Medium. Buffer solution is mainly water and its viscosity is 8.90 $\cdot 10^{-4}$ $Pa \cdot s$. Porosity for nitrocellulose is at 78.9%3. Since the medium is homogeneous along its length, we assume that the constrictivity $\delta$ is 1. For roughly cylindrical shapes tortuosity $\tau$ can be estimated as $\frac{1}{\epsilon_t}$4.

Analyte radius. To calculate the radius of an irregular shape like an oligonucleotide, we consider it as an ellipse with semi-axes a and b, and consider its radius r to be $\sqrt{ab}$. The width of a nucleotide is 1 nm and the length of an oligonucleotide is $0.33 \cdot (\# \text{of nucleotides})$ nm. Thus $r = \sqrt{0.33 \cdot (\# \text{of nucleotides})}$ nm, which for an oligonucleotide with 104 monomers is roughly 5.86 nm.

Gold nano particle radius. The radius of a gold nano particle is 6.5 nm, and it has 40 nucleotide tails in each direction, increasing the radius by $\sqrt{0.33 \cdot 40}$ = 3.6 nm to a total radius of 10.1 nm.

Complex radius. We assume that the radius of a complex is the radius of a gold nano particle + the attached analyte with extra 20 nucleotides of an unattached tail, so it is 6.5 nm + 6.4 nm = 12.9 nm.

Effective diffusion coefficients. Hence, the effective diffusion coefficients are $D_{A} = 26.1 \cdot 10^{-12}$ and $D_{P} = 15.1 \cdot 10^{-12}$, $D_{C} = 11.8$ $m^2 s^{-1}$. The rest of the particles are fixed.

\section{Development of the model} First, consider the ideal case without diffusion or advection. Then we postulate that the formation of complex is described by the equation \begin{equation} \dot{c} = k_1 (a_0 - c)(p_0 - c) - k_{-1} c,\\ \end{equation} where \(a_0\) and \(p_0\) are the initial conditions of the analyte and nanoparticles, respectively.

The solution \begin{equation}\label{ODE solution} c = c_1 c_2 \frac{1-e^{k_1 (c_2 - c_1) t}}{c_1 - c_2 e^{k_1 (c_2 - c_1)t}} \end{equation} is a monotonically increasing function that does not reach a maximum in finite time. $c_1$ and $c_2$ are the two roots of the equation \begin{equation}\label{complex relation} (a_0 - c)(p_0 - c) = K c,\end{equation} \(c_1 < c_2\). And the theoretical maximum is \begin{equation}\label{maximum}c_1 = \frac{1}{2}(a_0 + p_0 + K - \sqrt{(a_0 - p_0)^2 + 2K(a_0 + p_0) + K^2}) < min(a_0,p_0),\end{equation} which is an increasing function in both $a_0$ and $p_0$. $K$ is $\frac{k_{-1}}{k_1}$.

Equation (\ref{ODE solution}) shows that the rate of complex formation is directly proportional to the exponent $k_1 (c_2 - c_1)$ which achieves minimum values along the line $a_0 - p_0 + K = 0$.

Experiment. 9 lateral flow assay strip tests were done. Gold nanoparticles were added on a conjugate pad of volume $2.46 \cdot 10^{-8} m^3$ at fixed initial concentration $p_0 = 11.4nM$. $100 \mu L$ of analyte were added at concentrations starting from $a_0 = 500nM$ and diluted by a factor of 2 down to a concentration of $a_0 = 1.95nM$. After 20 minutes we observed the results.

Result. It was found that at some cut-off threshold of the analyte concentration no signal was visible (at $a_0$ = 15nM the signal is still visible, but at $a_0$ = 7.8nM no signal was visible).
There could be several explanations:

  1. There exists a minimum threshold concentration $c_{min}$ of the complex that allows for the signal $S$ to be observable.
  2. The rate of formation of the complex was too slow and thus the location of the test line was not far enough for the complex to form in sufficient amounts.
  3. Each gold nanoparticle can bind to several analyte molecules, and maybe several (more than one, on average) analytes have to bind to a nanoparticle before the resulting complex has a significantly large probability to bind to a probe $R$.

We consider each. Model predictions. Simulations of the advection-diffusion-reaction equations show that diffusion and advection act to reduce the "effective" concentration of the substrates. Hence, the ideal solution overestimates the real concentration. Since the experimental conditions relative to diffusion and advection are the same in each experiment, we use equation (\ref{ODE solution}) to calculate the reference concentration values $c_{ref}$. We take $k_1 = 10^{6} s^{-1}$ and $k_{-1} = 10^{-3} s^{-1}$.
First, we note that there is no need to add more than 30 nM of the analyte because at this point the complex formation saturates at the initial gold nanoparticle concentration. Second, if a threshold complex concentration exists, then we postulate that it is between $c_{lower} = 6.6 nM$ and $c_{upper} = 9.7 nM$.
The time $t_{ref}$ in FIGURE refers to the moment the complex reaches $0.8 c_{max}$ to check if the slowest complex formation rate is achieved near the theoretically computed $a_0 = p_0 - K = 10.4$ nM.

Experiments. To see if the visual threshold is real, we conducted two new tests. First, we fixed $p_0$ at 6.51 nM $ = c_{lower}$ on 6 strips, and added 100 $\mu L$ of $a_0 =$ 62.5, 31.25, 15.62, 7.8, 3.9 and 0 nM on each, respectively. Theoretically, at this point we should not see any visible signal at all, because the complex concentration cannot exceed the gold nanoparticle concentration, which is set below the visibility threshold.

Second, we fixed $p_0$ at 4 times the initial amount $= 46$ nM. Then we invert the equation (\ref{complex relation}) into \begin{equation}\label{a from complex} a_0 = c(1-\frac{K}{c-p_0}) \end{equation} and compute the analyte concentrations for $c=c_{lower}$ and $c_{upper}$ that yield $a_0 = 6.7$ nM and $a_0 = 10.1$ nM. Then we set an experiment with 100 $\mu L$ of analyte added on 6 strips at concentrations 50, 25, 12.5, 6.75, 3.37 and 0 nM. By hypothesis, we should not see the last three test lines.

Lastly, we added 6 more strips with the $p_0 = 11.4 nM$ and $a_0 = $ 62.5, 31.25, 15.62, 7.8, 3.9 and 0 nM as control. After 20 minutes we looked at the results.

Results. This time the control produced two clearly visible lines at $a_0 = 62.5$, 31.75 nM and a barely visible result for $a_0 = 7.81$ nM. This shows that a potential threshold is lower, but we cannot reject the hypothesis yet. Based on this result we reset $c_{upper}$ and $c_{lower}$ as 6.6 nM and 3.47 nM.
For $p_0$ = 6.6 nM we saw two barely visible lines for $a_0$ = 62.5 nM and 31.5 nM. This result is consistent with the estimate that $c_{thresh}$ is between 6.6 nM and 3.47 nM, putting it at around 6 nM of complex concentration.
For $p_0$ = 46 nM we saw two barely visible lines for $a_0$ = 50, 25 nM, but by hypothesis with redefined thresholds we should have seen at least faint signals for $a_0$ = 12.5, 7.8 nM, which were never observed. The control line, however, had a very high expression, which agrees with having a high gold nanoparticle concentration. From this we conclude that it is not necessary to have a large concentration of nano particles, as long as the concentration of the analytes is large enough.

Conclusion. We cannot reject the minimum threshold hypothesis, but this result is inconsistent with the model and shows that the threshold hypothesis cannot be the only explanation for why we do not see a signal below some point. After calculating the results, we note that the signal was not seen whenever $a_0 < \frac{1}{2} p_0$. This gives credibility to the third hypothesis that gold nano particles need to bind several analytes, on average, before the complex can bind probes with high probability. Consequently, we conclude that we need a new model.

References

  1. Ms, R. & Cm, A. A mathematical model to predict the optimal test linelocation and sample volume for lateral flow immunoassays. Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Annual International Conference 2012, 2408–2411 (2012).
  2. Berli, C. L. A. & Kler, P. A. A quantitative model for lateral flow assays. Microfluidics and Nanofluidics 20, 104 (2016).
  3. Das, D. & Panigrahi, P. CFD simulations for paper-based DNA amplifi-cation reaction (LAMP) of Mycobacterium tuberculosis—point-of-care diagnostic perspective. Medical & Biological Engineering & Computing 58, 271–289 (2020).
  4. Tjaden, B., Cooper, S. J., Brett, D. J., Kramer, D. & Shearing, P. R. Onthe origin and application of the Bruggeman correlation for analysing transport phenomena in electrochemical systems. Current Opinion in Chemical Engineering. Nanotechnology / Separation Engineering 12, 44–51 (2016).
  5. Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A.,Richardson, C., Ring, J., Rognes, M. E. & Wells, G. N. The FEniCS ProjectVersion 1.5. Archive of Numerical Software 3. Number: 100 (2015).
  6. Logg, A. & Wells, G. N. DOLFIN: Automated finite element computing. ACM Transactions on Mathematical Software 37, 20:1–20:28 (2010).