Solving Problems with Precision

Modeling is a vital part of the design process for any project because it allows variables that cannot be feasibly tested in Wet Lab to be simulated. Project modeling is more than analyzing important variables; it accelerates the engineering design cycle and saves resources.
For our project, modeling has provided a streamlined process for wet lab experiments by testing variables that are in question and optimizing the conditions. We used modeling to predict the 3D structures of antibodies that we are working with, which would normally take several months to determine experimentally. We also focused on LFA optimization to predetermine which parameters would produce the best results in the lab.

model program
References

1. Holvoet, P. et al. Circulating Oxidized LDL Is a Useful Marker for Identifying Patients With Coronary Artery Disease. Arteriosclerosis, Thrombosis, and Vascular Biology 21, 844–848 (2001).
2. Protein Molecular Weight. https://www.bioinformatics.org/sms/prot_mw.html.
3. Gasperino, D. Improving Lateral Flow Assay Performance Using Computational Modeling | Annual Review of Analytical Chemistry. https://www.annualreviews.org/doi/10.1146/annurev-anchem-061417-125737 (2018).
4. Mass action kinetics - Mathematics of Reaction Networks. https://reaction-networks.net/wiki/Mass_action_kinetics.
5. Landry, J. P., Ke, Y., Yu, G.-L. & Zhu, X. D. Measuring Affinity Constants of 1,450 Monoclonal Antibodies to Peptide Targets with a Microarray-based Label-Free Assay Platform. J Immunol Methods 417, 86–96 (2015).
6. Liu, Z. et al. The effect of report particle properties on lateral flow assays: A mathematical model. Sensors and Actuators B: Chemical 248, 699–707 (2017).
7. Glinowiecka-Cox, M. Analytic Solution of 1D Diffusion-Convection Equation with Varying Boundary Conditions. (Portland State University, 2022). doi:10.15760/honors.1224.
8. T--Rochester--lfam.pdf.

Protein Modeling

The original aim for protein modeling was to model the sandwich complex between the labeled antibody and the test-line antibody, which indicates a positive test. The plan was to predict the epitopes on the biomarker OxLDL that our labeled full-length IgG antibody and single-chain variable fragments (scFv) would bind to. By knowing these epitopes, we could determine which scFv fragment out of the variety of candidates would not interfere with the binding of the labeled antibody.

We first generated the structures of the three antibodies and apolipoprotein B-100 (apoB-100), which is the carrier protein on OxLDL. We successfully created Protein Data Bank (PDB) structure files for each antibody using AlphaFold. However, because of the size and complexity of apoB-100 due to the lipid composition of LDL, we were unable to model the complete structure of the biomarker. This meant we were unable to utilize docking softwares to predict the binding of the biomarker with our lateral flow assay components.

However, with further reading and meeting with researchers in the field, the modeling team was able to determine an ideal sequence for our scFv. The protein models of the lateral flow assay components were ultimately utilized in determining which strain of E. coli to synthesize the scFv and nanobody in. Since each of the selected scFv fragments and the nanobody contained disulfide bonds, the fragments must be expressed in SHuffle E. coli cells instead of another strain of E. coli.

McPC603 scFv Ribbon and Surface Model

This is a model generated from AlphaFold of the McPC603 antibody fragment. We generated models like these for our other antibodies to visualize them. Scroll to zoom in and out and click/drag to move the model around!

Lateral Flow Assay Modeling: Introduction

The Goal

Modeling works closely with wet lab to make major design decisions in the construction of the lateral flow assay (LFA). In this project, the goal for modeling is to determine the optimal (i.e. the minimum) concentration of detector units (full length IgG antibodies conjugated to a gold nanoparticle) that need to be placed on the sample pad in order to produce a visible signal at the test and control lines when at least 0.800 mg/dL (or 800 mg/L) of oxLDL is present in the blood sample.1 This is the threshold for oxLDL circulating in the bloodstream that indicates early development of atherosclerosis. Our team determined that the optimum concentration of the detector units is the minimum possible concentration because gold nanoparticles are the most costly portion in building our LFA pilots. Thus, the work in modeling necessarily reduced the time and cost of producing the physical test strip.

Constraints

1. The test must produce a signal with 0.800 mg/dL of oxLDL in the blood sample. We aim for our test to be sensitive to this stage of atherosclerotic development because it marks the stage at which users are at high risk for atherosclerosis but still have time to halt or even reverse the development process (through healthy diets, increased exercise, etc).

2. We will focus on optimizing the amount of detector units we place onto the test strip. So, the other parameters that are required to create a LFA are fixed as follows:

a. The distance between the end of the conjugate pad and the test line is 20.0 mm
b. The distance between the test line and the control line is 10.0 mm
c. The concentration of receptor antibodies (the antibodies at the control line and the test line) are 10 nM

Remarks

We will consider all concentrations in terms of mg/L. This is due to the unknown molar mass of oxLDL. Currently, the concentration of oxLDL in the blood which correlates to a high risk of developing atherosclerosis is only reported in terms of mg/dL. Thus, we are unable to convert this threshold into molar concentrations. So, we used an online protein calculator to calculate the weight of each component (listed below) of the LFA.2

The Problem Setup

A standard LFA, including our design, is composed of four different pads: the sample pad, the conjugate pad, the nitrocellulose membrane, and the absorbent pad. See the design page for a walkthrough of how an LFA works.

Scenario 1: At least 0.800 mg/dL of oxLDL is present in the sample

oxLDL will bind to the gold-conjugated full-sized IgG on the conjugate pad. Then, the oxLDL-IgG-gold nanoparticle complex will diffuse and advect from the conjugate pad onto the nitrocellulose membrane. The complexes will diffuse resulting in a net movement from an area of high concentration to an area of low concentration. As these particles move in a net direction during diffusion, the liquid sample will further be transported by the momentum of the fluid’s bulk motion, also known as advection, thereby moving the complexes over the test strip. These complexes will bind at both the test and control lines. Specifically, oxLDL will bind to both types of single chain variable fragments (scFv) present at the test line. The gold-conjugated full-sized IgG, also known as the detector unit, will bind to the control line receptor nanobody. This results in a visible line at both the test and control lines (i.e. a positive test).

Scenario 2: Less than 0.800 mg/dL of oxLDL is present in the sample

The sample will diffuse into the conjugate pad but the gold-conjugated full sized IgGs will not form a new complex with the sample. Then, these gold-conjugated full sized IgGs, the detector unit, will both diffuse and advect onto the nitrocellulose membrane. These detector units will not bind to test line receptor scFvs. The detector units will instead continue to move to the control line, where the constant region of the full length IgG of the detector unit will bind to the control line receptor nanobody.

Components and terms

Variable Term Molecule
A Analyte oxLDL (biomarker)
D Detector Unit Gold nanoparticles conjugated to full sized IgG antibodies
AD Analyte detector unit complex oxLDL-conjugated-full-sized-IgG complex
T1 IK17 scFv Test line receptor antibody
T2 McPC603 scFv Test line receptor antibody 2
T1AD Test line receptor antibody 1 analyte detector unit complex IK17-detector-unit complex
T2AD Test line receptor antibody 2 analyte detector unit complex McPC603-detector-unit complex
C Control line receptor antibody VHH (nanobody)
CAD Control line receptor antibody analyte detector unit complex VHH-oxLDL-detector unit complex

Remarks

We will consider all concentrations in terms of mg/L. This is due to the unknown molar mass of oxLDL. Currently, the concentration of oxLDL in the blood which correlates to a high risk of developing atherosclerosis is only reported in terms of mg/dL. Thus, we are unable to convert this threshold into molar concentrations. So, we used an online protein calculator to calculate the weight of each component (listed below) of the LFA.2

The Problem Set-Up: Developing and Using Differential Equations and Systems of Equations

Assumptions

Before using useful equations to model the concentration of each component, we need to make a few assumptions that both simplify the problem and allow us to combine certain equations.

  • The analyte, the detector antibody, and the receptor antibodies are dissolved in a well mixed solution formed by the sample wetting the test strip.
  • There is negligible absorption of the analyte, detector, and the analyte-detector complex into the nitrocellulose membrane pores (i.e. no loss in concentration in A, D, nor AD due to the nitrocellulose).
  • Each component reacts with their target at a constant rate (i.e. kinetic rate constants are constant with respect to concentration).
  • All concentrations are monotonic (i.e. the concentrations only increase or they only decrease over time and space). We assumed before that each reaction’s rate is constant, which means the rate is monotonic until equilibrium. After equilibrium is reached, the concentrations of all products and reactants remain constant.

Using the Advection Dispersion Reaction Equation (ADRE)

Our primary interest for this project is to examine the behavior of the concentration of each component (listed in the above table) over the length of the test strip and over time. These relationships can be described using the the Advection Dispersion Reaction Equation. These equations have been cited by other researchers as innovative methods to optimize the performance and design of an LFA. These articles, as well as our team, are interested in how the transport of the detector unit AD to the test and control lines affect the concentrations of AD bound at the test and control lines.2,2 The ADRE is derived from the conservation of mass for a transported chemical, which means these equations provide the insight needed to analyze the transport of AD over the test strip. For example, a question we could consider using the ADRE is if the reaction between A and D occurs before they reach the test line; if this reaction happens too slowly, the test could result in a false negative.

Using Mass Action Kinetics Equations

Assumptions 1 and 3 allow us to use mass action kinetics to define the concentration of the “fixed” components of our test. This means we have equations to describe the concentrations of the receptor antibodies anchored at the test and control lines over time. Additionally, we observe that these concentrations over time vary only due to the association and dissociation of the receptors and their respective targets.4 Furthermore, we accounted for the association and dissociation of AD as the complex moves along the test strip. That is, the concentration of AD varies over the length of the test strip not only due to the advection of the bulk fluid and the diffusion/dispersion of the complex molecules, but also because AD has a possibility of dissociating into A and D, or A and D in the solution associating to form AD. We accounted for this cause of AD concentration fluctuation by adding/subtracting appropriate expressions derived from mass action kinetics equations in the ADRE.

However, we could not find the appropriate rate constants for each binding reaction from literature. So, we used an approximation from which found the rate constants for antibodies with similar equilibrium constants as our antibodies/antibody fragments.5 So, we take the rate constant for the forward binding reactions to be 10.05 M/s and the rate constant for the backward dissociating reactions to be 10.0-5 M/s.

Reaction Chemical Equations

The following chemical reaction equations are needed to describe the mass action kinetics equations. The forward/backward rate constants are also denoted by their respective equations.

kfAD denotes the forward rate constant and kbADdenotes the backwards rate constant for this reaction.kfT1AD denotes the forward rate constant and kbT1AD denotes the backwards rate constant for this reaction. kfT2AD denotes the forward rate constant and kbT2AD denotes the backwards rate constant for this reaction. kfCAD denotes the forward rate constant and kbCAD denotes the backwards rate constant for this reaction.

Breaking the Lateral Flow Assay (LFA) into Three Systems

We chose to separate the LFA into three discrete systems of partial differential equations. We separated these systems based on the placement of the conjugate pad, the test line, and the control line (shown in the figure below). In doing so, we can define the behavior at the system boundaries (described in subsequent paragraphs), which contains variables we wanted to solve for, more easily. We could also characterize the sudden decrease in AD concentration after AD passes over the test lines this way. Specifically, we break the LFA into the following three parts:

  • Leftmost edge of the conjugate pad to the rightmost edge of the conjugate pad
  • Rightmost edge of the conjugate pad to the test line
  • The test line to the control line

Three systems the LFA was divided into. 1. Leftmost edge of the conjugate pad to the rightmost edge of the conjugate pad, 2. Rightmost edge of the conjugate pad to the test line, 3. The test line to the control line

Partial Differential Equation (PDE) System

We will write one partial differential equation for each component of the test strip for each system of differential equations. Conceptually, our goal is to define the concentration over each of the three sections of the LFA. These equations allow us to characterize the concentration as a function of both position on the test strip and time passed since loading the sample onto the sample pad. We can also begin to identify the phenomena that affect the concentration of each component.

The ADRE requires the diffusivity constant for the analyte (DA), detector unit (DD), and the analyte-detector-unit complex (DA +DD). We take these constants to be the same from Lui et al,.6 This means DA = 10.0-4 mm/s2 and DD = 10.0-6 mm/s2. Since we assumed the diffusivity of the analyte-detector-unit complex is the sum of their respective diffusivity constants, we let DA+DD = 10.0-4 + 10.0-6mm/s2.

The Role of Mass Action Kinetics and Rate Constants

As shown above, mass action kinetics are leveraged to predict the concentrations of the components fixed to the LFA (i.e. the receptor antibodies). For example, take equation 1 System 1. This equation considers the concentration of the analyte as it moves into the conjugate pad from the sample pad. While the concentration differs over time, as well as over the pad, due to the movement of the analyte, the concentration must also differ due to the analyte reacting with the detector unit. The rate of decrease in analyte concentration due to the analyte reacting with the detector unit is exactly the forward reaction rate constant multiplied by the concentration of the analyte multiplied by the concentration of the detector unit due to mass action kinetics. Similarly, the rate of increase in analyte concentration due to the analyte dissociating with the detector unit is also governed by mass action kinetics. By taking these reactions into account for each equation, we are able to more accurately determine the behavior of concentration over space and time for each component.

Solving the Systems of Equations (Using MatLab)

To solve partial differential equations, two important conditions are needed for every equation.

Initial Conditions

Initial conditions are the solution functions defined at time zero. Conceptually, these conditions are known concentrations at all positions of the system starting at time to or t = 0.00 sec. We define to to be the moment the sample of blood reaches the conjugate pad from the sample pad. An example of an initial condition is the concentration of D in system 1; our objective is to find concentration of D, which will be uniformly distributed in the conjugate pad. We defined the initial condition of every component of every system to be 0.00 nM = 0.00 mg/dL.

We define the initial conditions of every component in every system to be 0 mg/mL except for [D] in system 1 because no component except for [D] in system 1 is uniformly distributed throughout the entire length of the system. For example, the receptor antibodies only occupy a line within their respective systems, not the entire portion of the nitrocellulose membrane.

Boundary Conditions

Boundary conditions are the system’s behavior at edge points at any given time t. In this case, the boundary conditions are the known concentrations (or known behavior of concentrations) of specific components defined at the leftmost and rightmost edges of each system.

Left Edge Boundary Condition

The left edge boundary condition for A of system 1 would be the concentration of oxLDL at the leftmost edge of the conjugate pad at any given time. We let xo = 0.00 mm denote the leftmost edge for each system and we let xf= L mm denote the rightmost edge for each system; L is the length of the corresponding pad of the LFA.

Dirichlet Boundary Conditions

Notice, we used the Dirichlet Boundary Conditions for components with known concentrations at each boundary.

Neumann Boundary Conditions

Neumann boundary conditions were used for components with unknown concentrations at each boundary.7

We define the boundary conditions as follows.

System o mm L mm
System 1 [A] = 8.00, [AD] = 0.00 [AD] is the same as thee concentration of AD at the rightmost edge of System 2
System 2 [A] = 0.00, [AD] = 0.00
[AD] is the concentration of AD that must be present from system 1 to produce enough [T1] = 2.67*10-3
[T1AD] = 0.00, [T2] = 7.18 * 10-3, [T2AD] = 0.00 		[AD] is the same as the right edge boundary of System 3
System 3 [A] = 0.00, [D] = 0.00
[AD] is the variable we are solving for in this system (e.g. how much AD should be present at the beginning of the system to produce a visible signal
[C] = 1.50*10-3, [CAD] = 0.00 		[AD] = 0.00

We set the boundary conditions of [AD] (the concentration of AD) to be 0 at the right edge of system 3 because we want to minimize the excess AD (i.e. unnecessary AD) that an LFA must be constructed with. Also, we consider the left edge boundary conditions of [A] = 0.00 and [D] = 0.00 for systems 2 and 3 as a simplifying assumption; instead of defining an interpolated function, we decided to assess the approximated reaction rates and determine that the decrease in concentration due to dissociation for the following complexes are negligible (too close to zero to significantly impact concentration): T1AD, T2AD, and CAD.

Finally, since we placed the right edges of systems 2 and 3 at areas of the test strip that contain receptor antibodies, we can define the corresponding boundary functions: the concentrations of the respective antibodies is 10.0 nM (which was then converted to mg/L). For example, the right edge of system 3 is the control line. So, we know that the boundary condition for [C] in system 3 is the initial concentration of control line receptor antibodies we would place at the control line.

Solving the PDE System

We employed MatLab’s pdepe function to solve the three systems. Based on the results of the Rochester 2020 iGEM team model, we chose 6.40 nM to be the minimum concentration of gold nanoparticles that needs to bind at the control line to produce a visible signal.8 By using stoichiometric conversions, we found this minimum concentration to be equivalent to 1.26 * 10(-3) mg/L. Using this threshold, we implemented our MatLab model to iteratively solve System 3 using different starting concentrations of the AD until the concentration of CAD met the threshold. In other words, we solved for ADc or the concentration of AD that needed to reach the control line in order to produce a visible signal.

Then, we used this solved value ADc as the right edge boundary condition for system 2. This means we needed a concentration of ADc to be left after the complexes T1AD and T2AD formed at the test line. In MatLab, we iteratively solved for the initial concentration of AD of system 2, or ADT, that produced a solution satisfying the right boundary condition. So, ADT is the concentration of AD needed at the start of the nitrocellulose membrane in order to produce a signal at both the test and control lines. We ensured the solver found the solution that produces a visible solution at the test line by setting the right boundary conditions of T1AD and T2AD to be 6.40 nM, or the concentration of gold nanoparticles that needs to bind at the test line to produce a visible signal.

Finally, we assigned ADT to be the right boundary condition of ADi or the concentration of AD that must be present after A and D bind on the conjugate pad. We also specified the right boundary conditions of A and D to be zero (once the complexes run to the absorbent pad), in order to find the optimal initial concentrations of A and D at the start of the test. We then used these boundary conditions to iteratively solve for the initial concentration of D; we fix the initial concentration of A to be 8.00 mg/L or the minimum concentration of oxLDL we expect to be on the test and produce a visible signal.

By solving the three systems of PDE in this way, we are able to mathematically find the initial concentration of D to produce visible signals at the test and control lines if the expected minimum threshold of oxLDL is present in the initial blood sample.

Results

Informing Physical Design of AtheroSHuffle

After solving the three systems in the manner as described above, we were able to conclude that the ideal minimum concentration of detector unit is 5.61*10-4 mg/L. In other words, the model predicts that this concentration of the detector unit loaded onto the conjugate pad of the LFA is sufficient to generate a visible signal on the test and control lines when at least 8 mg/L of oxLDL is present in the bloodstream.

Model Applications

In addition to predicting the minimum concentration of detector unit to generate a visible signal, the model is also informative of the behavior of each LFA component when 5.61*10-4 mg/L of the detector unit is used on the conjugate pad. For example, we can examine the behavior of the CAD complex at this optimal detector unit concentration.

Fig 2. This figure depicts the concentration of CAD or the complex-analyte-detector-unit complex as a function vs. time.

As shown in the figure above, we observe that the CAD concentration will rapidly begin to increase at around 30.0 seconds. Then, due to mass action kinetics, the concentration will peak at over 1.20 * 10-3 mg/L, which means a visible signal will be produced. Furthermore, we can see that this concentration will last a very short amount of time. So, the model determines that the concentration of CAD will indeed reach the minimum concentration of gold nanoparticles bound to the control line in order to produce a visible signal but will not remain at that minimum concentration for more than one second. So, one model improvement would be to find the minimum concentration of D such that the concentration of CAD remains above the minimum threshold for at least 10 seconds.

Model Improvements

One important parameter that the model cannot predict is the initial volume of the sample that travels through the LFA. However, the minimum volume threshold can be determined in Wet Lab without excessive reagent or resource usage. Specifically, we can test this in Wet Lab using a few drops of synthetic blood (a glycerol-water solution) and a formulated chase buffer. Our team prioritizes the user’s comfort in using this LFA in an at-home setting, which means we prioritize minimizing the amount of blood needed to be drawn for this test to work. So, we initially hypothesized the initial volume of the sample would impact the ability of the detector units to move up the test strip. We were concerned that the detector units would remain at the start of the test strip with an insufficient volume of fluid to allow for the diffusion/dispersion/advection of the detector unit and analyte. See the Results page for more details on these experiments. Briefly, we observed that AD complexes were able to flow down the entire length of the membrane when they were followed with a few drops of chase buffer. From a user perspective, a solution we can provide to this issue includes equipping the LFA kit with an aliquot of the chase buffer we formulated in the aforementioned experiments. The users would be instructed to add a few drops of chase buffer after adding their blood samples to aid in the movement of the protein components of the LFA. However, with future quantitative development, we will find pre-existing models and/or create new models to account for the initial volume of sample. Specifically, we should examine the effect of varying the initial amount of oxLDL on the amount of detector unit required to produce a visible signal at the test and control lines, rather than their respective concentrations.

Another improvement of the model includes generating data for the constants utilized. Specifically, using the reaction rate constants specific to the antibodies and antibody fragments used on the test strip as well as the diffusivity constants of oxLDL and the gold-conjugated full sized IgG. Since oxLDL is one of the largest macromolecules known to date, the diffusivity of oxLDL through the nitrocellulose membrane may be smaller than the variable found from literature. Additionally, the rate constants were estimated using the equilibrium constants; that is, the rate constants of antibodies with similar equilibrium constants were used in place of the rate constants that correspond to the components of the test strip. Since there is no correlation between thermodynamic and kinetic constants, these values may not estimate the rate constants of the components of the lateral flow assay. These refinements would allow the model to more accurately predict the concentration as a function of time and space.