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 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} \min_{L \in [0, X]} J(L) = -S(t_{end}) \end{equation} subject to the advection-diffusion-reaction (ADR) system \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 literature: \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 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%. 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}$.
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 $(# 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}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 in turn: 1. 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.
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.
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.
We have found empirically that the equation \begin{equation} \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^{-5}$, $n = 2$, $k=2$.
From this phenomenological model we made predictions that with very low gold nanoparticle concentrations (3.7 nM = $c_{thresh}$ and 6 nM = $c_{lower}$ from first experiment), adding 12, 6, 3 and 1.5 nM of analyte would produce a result for the first two analyte concentrations, and with higher concentrations (11.4, 23 nM), adding 40, 20, 10, 5 nM of analyte, we would only see results when $a_0 > \frac{1}{2} p_0$.
The results were consistent with the predictions of the model. 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.

Numerical PDE-constrained optimization software tools are, in general, not freely-distributed, even those that are adapted for Python, and hence cannot be used for an open-source application. And for such solvers the number of computational steps scales with the number of PDE simulations and optimization iterations.
To make the app freely available and fast, we take 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 is a trade-off between the rate of formation of the complex C and the velocity V(t) of the lateral flow.
We chose to implement the model using the computational finite element method platform FEniCS. This was done because it is open-source, has a Python interface, and its back-end is implemented in C++, thus 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.