Team:Vilnius-Lithuania/Model
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.
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.
Video 1. The velocity estimation experiment. Different buffer concentrations (50, 60, 70, 80, 90, 100, 110 $\mu L$) were applied onto strip membranes and set next to a ruler. The experiment was recorded with a video camera. The result shows a clear exponential fit. The velocity was robust to changes in volume in the range 60-100 $\mu L$, but started increasing afterwards. 50 $\mu L$ test was discarded because the solution never permeated into the strip membrane.
Video 2. The second velocity estimation experiment. Buffer concentrations up to 200 $\mu L$ were considered. At volumes 150 $\mu L$ and larger the buffer started overflowing. We conclude that valid volumes range 60-125 $\mu L$.
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 literature2 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$.
Figure 1. Fitted time-series data for the lateral flow velocity. The result shows a clear exponential fit.
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.
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$.
Figure 2. Simulation of the complex formation for the 9 wells ($a_0$ = 500, 250, 125, 62.5, 31.25, 15.63, 7.81, 3.91, 1.95 nM) at $p_0$ = 11.4. The upper threshold concentration $c_{upper}$ marks the point where the signal was still visible, and the $c_{lower}$ where it was not. The complex saturates at around $a_0$ = 31.25 nM. $t_{ref}$ is the time the complex reached 0.8$c_{max}$ concentration.
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:
- There exists a minimum threshold concentration $c_{min}$ of the complex that allows for the signal $S$ to be observable.
- 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.
- 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 3 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.
Figure 3. Simulation of the complex formation when gold nano particle concentration is below the lower threshold at $p_0$ = 6.51 nM. The analyte concentrations are at $a_0$ = 62.5, 31.25, 15.63, 7.81, 3.91, 0 nM. None of the results should have been visible, but the first two were. The lower threshold hypothesis gives inconsistent results. The upper and lower threshold concentrations are set the same as in the first control experiment.
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.
Figure 4. Simulation of the complex formation when gold nano particle concentration is above the upper threshold at $p_0$ = 46 nM. The analyte concentrations are at $a_0$ = 50, 25, 12.5, 6.25, 3.175, 0 nM. Three of the results should have been visible, but only the first two were. This suggests a different reaction model.
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.
Figure 5. The control experiment ($p_0 = 11.4 nM$) simulation with redefined upper and lower thresholds to match 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.
New model. Thus we form several potential ODE models that could describe the formation of a complex, including mass-action kinetics of highter order, Hill kinetics, and a cascade of differential equations, giving rise to several predictions on the effect of the initial analyte and gold nano particle concentrations. We calibrated these models using control experiments at $p_0$ = 11.4 nM where the signal was not observed for $a_0$ = 3.9 nM with $c_{thresh}$ corresponding to 3.7 nM. And we fix $c_{thresh}$ as a concentration below which we assume no signal is visible after 30 minutes.
We have found empirically that the equation \begin{equation}\label{empirical} \dot{c} = k_1 \frac{(a_0-c)^n (p_0 - c)}{(k p_0)^n + (a_0 - c)^n} - k_{-1} c \end{equation} explains experimental results, with \(k_1 = 10^{-3}\), $k_{-1} = 10^{-4}$, $n = 2$, $k=2$.
The formula predicts the experimental results by separating them into three classes: visible (above upper threshold), not-visible (below lower threshold) and indeterminate/random (intermediate values).
Experiment. From this phenomenological model we made predictions that gold nano particle concentration is not main the determining factor for the visibility of the test line, but the concentration of analyte. For very low concentrations of gold nano particles (3.7 nM = $c_{thresh}$ and 6 nM = $c_{lower}$ from first experiment), we saw the result for 120, 60, 30 and 15 nM of analyte, and with higher concentrations (11.4, 23 nM), adding 40, 20, 10, 5 nM of analyte, we saw the results only when $a_0 > \frac{1}{2} p_0$.
Figure 6. The model predictions for the visibility of the signal with an uncertainty region in the middle. $p_0$ = 11.4 nM.
Figure 7. The model predictions for the visibility of the signal with an uncertainty region in the middle. $p_0$ = 23 nM.
Figure 8. The model predictions for the visibility of the signal with an uncertainty region in the middle. $p_0$ = 6 nM.
Figure 9. The model predictions for the visibility of the signal with an uncertainty region in the middle. $p_0$ = 3.7 nM.
Conclusion. The model incorporates both the existence of a lower threshold along with the condition that there must be enough analyte in comparison to gold nano particles, suggesting that one analyte per gold nano particle is not enough to bind to a probe with high probability. We can use this model to evaluate if a signal will be visible from the initial concentrations.
We chose to implement the model using the Python-friendly computational finite element method platform FEniCS5,6. This was done because it is open-source, has a Python interface, and its back-end is implemented in C++, as a result it does not compromise speed.
FEniCS requires that the PDE is expressed in terms of a weak or variational formulation: \begin{align} \int_{\Omega} \partial_t [A]v_{A} + D_{A} \nabla [A] \cdot \nabla v_{A} + V(t) \cdot \nabla [A] v_{A} \, dx \\ + \int_{\Omega} \partial_t [P]v_{P} + D_{P} \nabla [P] \cdot \nabla v_{P} + V(t) \cdot \nabla [P] v_{P} \, dx \\ + \int_{\Omega} \partial_t[C]v_{C} + D_{C} \nabla [C] \cdot \nabla v_{C} + V(t) \cdot \nabla [C] v_{C} \, dx \\ + \int_{\Omega} \Big(k_1[A][P]v_{A} - k_{-1}[C]v_{A} + k_1[A][P]v_{P} - k_{-1}[C]v_{P} \\ - k_1[A][P]v_{C} + k_{-1}[C]v_{C} + k_2[C][R]v_{R} - k_{-2}[S]v_{R} \\ - k_2[C][R]v_{S} + k_{-2}[S]v_{S}\Big)\, dx = 0. \end{align} The functions $[A](t), [P](t), [C](t), [R](t), [S](t)$ are taken to be second-order Lagrange finite elements because the equation has quadratic non-linearities. Functions $v_\alpha$ are the corresponding test functions equal to zero on the boundary $\partial \Omega$. Neumann boundary conditions force any surface integrals to vanish. Time-discretization is done by Crank-Nicholson finite-difference scheme, because it leads to quick simulations, yet is sufficiently accurate and stable. Higher-order schemes cause the simulations to become too time-consuming. The computation of finite elements is done by Newton's conjugate gradient method.
To make the app freely available and fast, we require to use only open-source, freely-distributed software tools that achieve fast computations. Open-source condition is satisfied by Python, and we opted to use it. However, specialized PDE-constrained optimization tools cannot be distributed freely, and are, hence, not implementable as a back-end in an open-source software. Hence, we needed an alternative approach.
We reason as follows: the rate of formation of the signal S is directly proportional to the available concentration of complex C at the test line. Hence we argue that the location of the test line ought to be where the complex C concentration is maximized, which is a trade-off between the rate of complex formation and the lateral flow velocity.
The forward simulation of the advection-diffusion-reaction system runs with given initial conditions until the complex C reaches a threshold value at some point x, and this point is output as the optimal test line location. The threshold value is computed from equation (\ref{maximum}) as a fraction of a theoretical maximum.
Video 3. The simulation of the advection-diffusion-reaction system. The simulation estimates the optimal test line location. The simulation is run until the complex reaches a threshold value 0.9 $c_{max}$. The top simulation is for the complex, the middle one for analyte and lower one for gold nano particles. $a_0$ = 25 nM, $p_0$ = 11.4 nM. The velocity profile is taken to correspond to 100 $\mu L$ of analyte.
The second simulation runs the entire system and predicts the time it takes to obtain a signal.
Video 4. Simulation of the full ADR system which predicts the time a sufficient concentration is achieved on the test line.
To check if the model correctly predicts the optimal location, we simulated the system for a 25 nM analyte solution at 90% threshold with gold nano particle concentration 11.4 nM, and achieved the optimal line at 10mm. Having obtained the simulation results, we conducted corresponding experiments. We added different analyte concentrations (100 nM, 25 nM, 6.5 nM) onto strips with test lines at different locations: 5mm, 10mm (optimal for 25 nM), 12mm.
The hypothesis is that the middle test line location (simulated from onFlow) would obtain a more intense signal than the closest test line, and the most distant test line signal will not differ much from optimal location determined from software.
Video 5. To make sure that this location is the optimal one, we examined 2 more locations that are 0.5 cm further and 0.5 cm closer to the location obtained from the simulation. 100 nM, 25 nM, 6.25 nM solutions were dipped into three strips with test line locations at 0.5 cm, 1 cm (optimal for 25 nM), and 1.5 cm. Here the strip tests are in groups of three. The first group's strip is dipped into analyte solution where the concentration is 100 nM, second - 25 nM, third - 6.25 nM, and fourth - 0 nM. The first strip in each group has a test line at 1.5 cm, the second strip test line is at 1 cm and the third strip is at 0.5 cm. The hypothesis is that the middle test line location (simulated from onFlow) would obtain a more intense signal than the closest test line, and the most distant test line signal will not differ much from the optimal location determined from software.
References
- 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).
- Berli, C. L. A. & Kler, P. A. A quantitative model for lateral flow assays. Microfluidics and Nanofluidics 20, 104 (2016).
- 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).
- 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).
- 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).
- Logg, A. & Wells, G. N. DOLFIN: Automated finite element computing. ACM Transactions on Mathematical Software 37, 20:1–20:28 (2010).
contact us
follow us
