Efficient Simulation of Arbitrary Multicomponent First-Order Binding Kinetics for Improved Assay Design and Molecular Assembly.

Traditional enzyme-linked immunosorbent assay (ELISA), long the workhorse for specific target protein detection using microplate wells, is nearing its fundamental limit of sensitivity. New opportunities in health care call for in vitro diagnostic tests with ultrahigh sensitivity. Magnetic bead-based sandwich immunoassay formats have been developed that can reach unprecedented sensitivities, orders of magnitude better than are allowed for by the rate constants for a single ligand-receptor interaction. However, these ultrahigh sensitivity assays are vulnerable to a host of confounding factors, including nonspecific binding from background molecules and loss of low-abundance target to tube walls and during wash steps. Moreover, the optimization of workflow is often time-consuming and expensive. In this work, we present a simulation tool that allows users to graphically define arbitrary binding assays, including fully reversible first-order binding kinetics, timed addition of extra components, and timed wash steps. The tool is freely available as a user-friendly webapp. The framework is lightweight and fast, allowing for inexpensive simulation and visualization of arbitrarily complex assay schemes, including but not limited to digital immunoassays, DNA hybridization, and enzyme kinetics, for validation and optimization of assay designs without requiring any programming knowledge from the user. We demonstrate some of these capabilities and provide practical guidance on assay simulation design.

Introduction

The ability to sensitively detect specific biomarkers in a clinical sample containing a mixture of off-target components is a cornerstone of diagnostic medicine. The most common method by which this is achieved for proteins is the enzyme-linked immunosorbent assay (ELISA),[1] in which target proteins are captured either by direct adsorption to the surface of a plate or through precoated “capture” antibodies and subsequently detected by a labeled secondary “detector” antibody and optically read via colorimetric, fluorescence, or chemiluminescence detection strategies.[2−4] This technique is ubiquitous in life sciences and medicine to detect and quantify a specific protein in a complex mixture and is driving everything from pregnancy tests to cancer detection. However, as our understanding of the human proteome advances, there is a growing need for the detection of target proteins in the femtomolar concentration range, and standard ELISA is limited to target concentrations within a few logs of the dissociation constant (KD) for the receptor–ligand interaction on which it is based, which practically limits it in most cases to the picomolar to nanomolar range.

Recently, impressive work has been done using the so-called digital sandwich immunoassay schemes, in which just a few copies of a target protein can be detected and counted directly.[5−7] Digital counting methods overcome measurement uncertainty associated with the integration of an analog optical signal and is, in principle, only limited by Poisson counting noise, making femtomolar and in some cases even attomolar concentrations detectable.[8,9]

One of the best-known examples of this approach is the SiMoA technology pioneered by Walt and colleagues and now commercialized by Quanterix.[5−7,10] In their scheme, paramagnetic beads are coated with capture antibodies specific to the target protein and are mixed with a clinical sample. Because each bead has on the order of N ∼ 105 antibodies, the effective KD for a bead is N times smaller than the KD for any individual antibody alone, permitting efficient capture in 3D of very low concentration targets. Coupled with digital detection of the beads with or without a target protein bound, the SiMoA technology has improved the sensitivity of standard ELISA by >1000x.[11]

While this approach has made possible the quantification of low-abundance biomarkers from complex biofluids,[12] the development of an assay for a particular target remains an extremely laborious and expensive task, involving weeks of optimization, varying the assay steps, component concentrations, incubation times, number of washes, etc., to maximize the assay performance. To assist in these tedious experimental tasks, there is a need for simulation tools to shift the burden of optimization away from the expensive and time-consuming empirical framework.

A few studies have attempted to address this knowledge gap. Chang et al.[11] presented a model that involves a multistep, multicomponent framework for simulation of their assay workflows, which performs well in the low-concentration regime of target as compared to capture antibody (i.e., 1 or 0 target per bead).[11] However, their model assumes irreversible binding between additions of new components to the mixture, which is overly optimistic about the final signal generated. More recently, a Python library that allows the definition of arbitrary first-order coupled kinetics was introduced, which has similar capabilities similar to the tool presented in this work but relies on the user being able to program their own assays in Python.[13] Other computational resources with a slightly different specialization can be readily found,[14] but applying these tools can be a laborious task.

In this work, we present a user-friendly tool that allows the definition of arbitrarily complex multicomponent systems of first-order binding interactions, which allows for the timed addition of components in the middle of the workflow, and timed wash steps. The effects of the nonspecific binding of background molecules and the influence of the timing and duration of wash steps can easily be simulated to inform assay design before undertaking complex experimental work. We demonstrate its utility through the simulation of a full bead-based digital sandwich immunoassay, complete with nonspecific binding and wash steps, and discuss best practices for digital immunoassay design that arise from the conclusions of the model. As another use case, we employ this tool to model a hybridization reaction involving the assembly of small DNA nanostructures used in our previous work,[15,16] with the goal of aiding in explaining previously ambiguous experimental results.

Materials and Methods

Experimental

Solid-state nanopore data is reused from our previous work for the purposes of comparison to this tool. Briefly, nanopore fabrication is performed using the controlled breakdown method[17−20] on SiN membranes and used to sense DNA nanostructures, the design of which has been previously published.[21] Experiments are performed in 3.2 M LiCl, pH 8, at 100 mV using a 12 nm pore, and 8-pole low-pass Bessel filtered at 200 kHz for analysis.

Theoretical

A simple reversible first-order receptor–ligand binding or DNA hybridization reaction between arbitrary components A and B is governed by the following system of equationsandwhere [°] denotes the concentration, and kon and koff are the association and dissociation constants, respectively, which together define the dissociation constant . This can be readily extended to a multicomponent system with a vector of component concentrations using Einstein summation notation aswhere β is a rank-3 tensor of on-rates in such a way that β = β, for i ≠ j and i ≠ k, is the rate at which c is produced by binding between c and c and β for i ≠ k is the rate at which c is depleted as it forms various complexes with c. α is a matrix of off-rates such that α, for i ≠ j, is the rate at which c breaks into c and another component, and α is the total rate at which c is enriched by the breakup of all other components. To represent a physically valid system of first-order binding kinetics, α and β must satisfy conservation of mass requirements, which can be expressed asand

Prior to simulation, all variables are internally normalized to ensure numerical stability as follows. We define two calculated normalization constants specific to a given system of interactions:is the on-rate normalization constant, whileis the off-rate normalization constant. One could equivalently use the spectral norm of α to normalize, but we use the Frobenius norm here since it is easier to generalize to higher-order tensors. From these, we define the population normalization, which can be thought of as a generalized analog of KD for the system, as

Defining and , and , the normalized system of equations becomes

Note that all variables are rescaled back to the same units in which they are configured by the user prior to the actual output.

While most nontrivial reactions do not have analytical solutions, this system can readily be numerically solved using an explicit Runge–Kutta integration scheme.[22] The tool presented here uses the standard RK4 to propagate eq through time, with a time step of Δτ that is calculated at runtime to ensure numerical stability, under the assumption that the system will be farthest from equilibrium when the simulation starts. The model outputs a snapshot of the current concentration vector 10 times per τ. Since τ is by construction shorter than the shortest dynamic time scale in the system, this ensures that no interesting features are lost to sampling issues.

The model makes several key assumptions. We are considering only the evolution of concentration over time, which means that we are assuming that the probability of two species encountering one another is proportional to the ratio of their concentrations, or equivalently, that all of the components are perfectly mixed. We do not explicitly model diffusion nor are steric interactions considered. As such, it is important that this tool only be applied to the simulation of assays that occur in the bulk phase (i.e., in 3D). If antibodies are immobilized on a 2D surface, for example, it is not to be expected that this model will give correct results.

For the purposes of numerical simulation, we define a system at equilibrium to be one that satisfies

It is assumed that different components in the mixture will respond differently to washes. The use of magnetic beads, which can be pelleted with a magnet to pipette out the supernatant to remove all unbound molecules in solution, is one way that this can be achieved practically. In the simulation, wash steps are modeled by multiplying the concentration of all components by a component-specific “wash efficiency” factor between 0 and 1, which represents the fraction of the component concentration that is removed during the wash step.

The full user interface and a guide to using the tool for practical assay system and workflow definition are given in Supporting Information Section S1. The webapp is freely available at https://tcossalab.net/binding-assay/.

Results and Discussion

Model Validation

To validate our model, we first begin with a demonstration of a reaction that can be solved analytically, namely, a simple ligand–receptor system defined by the reaction A + B ↔ AB. We then reproduce results from Chang et al.[11] for the time evolution of the concentration of a generic first-order binding system, matching values of the concentrations, on- and off-rates, of all components for that particular system. The results from these simulations are shown in Supporting Information Section S2 and reveal a perfect agreement of our model with the existing method of Chang et al. in the low-concentration regime.

Idealized Digital Immunoassay

Digital immunoassay design is becoming an increasingly important challenge. As the limits of detection are pushed further down toward single-molecule copies, the details of the sample preparation steps leading up to the final detection step become critical in enabling greater sensitivity. To better understand the impact of upstream sample preparation and biochemical reactions leading up to the final measurement step, we simulate a variety of model systems, including a full assay workflow used in the SiMoA technology, to better understand the influence of the timing and duration of each step of a typical digital immunoassay.

The simulation setup for this model consists of four primary components: capture antibodies bound to magnetic beads (A), a target protein (T), a detector antibody (D), and a labeling molecule (L), such as an enzyme[6] or a DNA strand.[16] Assuming no cross-reactivity between species and no nonspecific background molecules, these four components can eventually bind together into an ATDL complex, forming every permutation of subcomplexes along the way. Our idealized assay consists of the 10 possible species representing the set of subcomponents: four base components (A, T, D, and L), three 2-component subcomplexes (AT, TD, and DL), two 3-component subcomplexes (ATD and TDL), and the full ATDL complex. These subcomponents can form or break up at any point in accordance with eq . We will assume for simplicity that the binding of a subsection of the complex does not alter the binding constants for subsequent assembly steps, though the simulation framework can accommodate the changing binding constants in response to partial complexation. This system is diagramatically shown in Figure a,b.

Figure 1

(a) Schematic diagram of a typical 4-component bead-based digital immunoassay, in which a capture antibody A bound to a paramagnetic bead captures targets T from solution, which are then labeled with a combination of a detector antibody D and a label L that is eventually used for downstream detection. On the way to assemble the full complex, all possible subpermutations of 2- and 3-components will form, leading to a total of 10 possible interacting species in solution. (b) Network diagram of interactions between the four primary components of a bead-based digital immunoassay and complexes thereof. Paired binding pathways are color-coded to highlight interaction partners in the network. (c) Typical assay workflow. Assay begins at the “Mix” step and follows the arrows in the labeled order.

Note that we are not explicitly simulating beads in this case, but instead simulate the capture antibodies as though they are uniformly distributed through the reaction volume. This is a subtlety that is nonetheless important: when a bead is coated with N capture antibodies, to first order, it acts as a single antibody with a kon that is N times larger than the kon for a single antibody alone, while koff remains unchanged. In the case where many of those binding sites are occupied, the effective on-rate reduces proportionately, that is, a bead with N binding sites of which n are occupied has an effective on-rate of (N–n) kon with respect to the next binding event, which means that it is not possible to simulate these beads directly unless we consider the time-dependent on-rate that results. Because these beads are spread through the three-dimensional volume of the reaction chamber, however, it is reasonable from a simulation standpoint to ignore those spatial correlations and to simply simulate the capture antibodies as though they were freely diffusing on their own, which circumvents the need to account for partial binding affecting the on-rate of the beads themselves. This is reasonable if the ratio of the capture antibodies to target molecules is large and steric effects are negligible, which will always be true for any real assay attempting to measure the concentration of a low-abundance target. As an example, in the SiMoA platform, the ratio of target to bead is <1 for a digital readout, while each bead contains hundreds of thousands of antibodies.[5,11] Nevertheless, care should be taken when using this model to simulate a situation in which target concentrations are comparable to the capture-antibody concentrations, as in this case, these steric effects may be important and this model would be expected to overestimate binding.

A typical paramagnetic microbead-based digital assay workflow generally follows the steps below. This workflow is shown in a flowchart in Figure c.

Capture-antibody-coated beads are mixed with the target and allowed to equilibrate.

Beads are immobilized and washed to remove the unbound target and any nonspecific background molecules. Note that some captured targets will dissociate during and after this step, although this is minimized by the enhanced kon if the number of capture antibodies is high relative to the number of targets.

Beads are resuspended, and detection antibodies and labeling enzymes are added and allowed to bind and label the captured targets remaining on the beads.

Beads are immobilized and washed to remove the unbound detector and label. Note that some captured target–label complexes, as well as labels themselves, will dissociate during and after this step.

Postprocessing done to count the fraction of labeled beads.

As an illustration of the capabilities of this software, we first approximate a generic bead-based immunoassay with a workflow similar to that which is used in the SiMoA technology[6,10] or to the nanopore electrical readout, which we conducted recently.[16] Since accurate values of the on- and off- rates are not always available for the capture antibody and detector antibody pairings in those digital assays, we instead make some approximations. To do so while still gaining useful physical insights, we normalize parameters by the KD value of the capture-antibody pairing. We thus set on- and off-rates for AT and TD binding pairs to 1 (arbitrary units, which translates to matching the KD for all pairings to that of the capture antibody) and on- and off-rates for DL to 10 and 0.01, respectively, to approximate the much stronger biotin–streptavidin interaction usually used to bind the label to the detector. Target concentration is set to 0.001, which can be thought of as being equal to 0.001KD for the target–capture-antibody pairing, consistent with the typical ultrasensitive digital immunoassays that use very high affinity picomolar KD pairings to detect femtomolar target concentrations. Capture-antibody concentration is set to 103 to simulate a ratio of the capture antibody to target of 106, which is typical of ultrasensitive assays when operated near the limit of sensitivity, while detector antibody and label, when added, are at a concentration of 102. These concentration ratios are typical of digital immunoassay workflows. To make the example concrete, this could correspond to using a very high sensitivity capture antibody with a KD of 1 pM using a 1 nM concentration of capture antibody to detect a 1 fM concentration of target using biotin–streptavidin labeling with a KD of 1 fM, which is near the limit of sensitivity for the best assays currently available.

Wash steps are modeled by setting the free concentration of all species that are not bound to a bead (in this case, complexes that do not contain an A) to zero, though the general framework allows for the effects of imperfect washing to be simulated as well by setting the wash efficiency to a number between 0 and 1. Three washes are performed after allowing the initial equilibration of A and T. Following the last wash, D and L are added, and an additional three washes are performed. Wash timings are set to 1 normalized time unit after each component addition event or to run to equilibrium in all cases, whichever is shorter, which simulates the typical digital immunoassay workflow of having fixed wash timings. Note that this means that equilibrium is not necessarily established before the washes during this stage, which is typical of experimental workflows. Figure c shows the full workflow being simulated.

Note that the normalization of units (or use of relative units) is not required in general. The webapp framework lets users choose units for time and concentration and enter all parameters in those units. Results are also reported in the chosen units. The use of relative concentrations for this example is to make clear the generality of the tool.

Figure presents the time evolution of the relative concentrations of the various interacting species in such a typical 4-component bead-based sandwich immunoassay. Figure a is a zoom into step 1, showing the time evolution of the capture of targets by antibody-coated beads. Equilibration of AT complexes occurs very fast due to the excess of capture antibodies, and there is essentially no loss of target to the wash steps since any targets that dissociate are immediately recaptured (note that a concentration of ∼10–6 represents aM in a typical SiMoA assay with a capture antibody with a KD of 1 pM). These results clearly show the minimal effect that washing has on the target at this stage, and it is during this part of the assay that washes are least disruptive to the overall downstream assay performance. Figure b then shows the addition of the detector antibody, D, and label, L, and the time evolution of all complexes that are not removed during the washes (anything containing an A, i.e., stuck to the bead). The concentration of the labeled sandwiched target on bead, ATDL, very rapidly approaches the initial target concentration of 0.001, implying that all targets are captured and labeled. Subsequent washes to remove excess free-floating detector antibodies and labels, once equilibrium is reached, lead to a marked decrease of that full ATDL complex concentration with time. In contrast, the concentration of AT complexes (just the target captured on bead) gradually returns to the initial concentration of targets as detector antibodies dissociate. On the other hand, the concentration of free targets in solution, T, remains very low (<10–6), and the concentration of targets sandwiched with an unlabeled detector antibody, ATD, is also low (∼10–5) and does not change much due to the high affinity between the detector and label. Concentrations in this part of the assay respond much more strongly to washes, and equilibrium takes orders of magnitude longer to establish at this stage than during the initial target–capture step. Figure c shows the total amount of T and L in the system (summing over all complexes containing those components), demonstrating the loss of signal that occurs as a function of wash timing and steps. While essentially no target is lost throughout the assay, the amount of label available to indicate its initial presence in the sample is highly dependent on the details of the wash steps following the addition of the detector antibody, D, and label, L.

Figure 2

Time evolution of the relative concentrations of the various interacting species in a typical 4-component bead-based sandwich immunoassay. (a) Fast equilibration of the target molecule T, with an initial concentration of 0.001, in the presence of an extreme excess of capture antibody A, forming AT complexes (∼99.99% of targets are captured very rapidly). D and L are added to the mixture between panels (a) and (b). (b) Time evolution of complexed species that do not get washed away during wash steps after the addition of D and L at a normalized concentration of 100. (c) Time evolution of the total target T and total label L available in the system, showing a strong loss of ATDL complexes to dissociation after washes that occur late in the workflow. Vertical dashed lines indicate wash steps.

From the results of Figure , a few things are immediately apparent. The first is that not all wash steps are equal: very little target is lost in the first round of wash steps (Figure a) to remove the nonspecific background molecules from the complex biofluids, while a majority of detector–label complex is lost in the second round (Figure b) used to remove excess, unbound labels. While it is possible to calibrate the assay to account for these losses, the signal will get progressively weaker with the number of wash steps that occur after the label is added, implying that wash steps should be front-loaded to the extent possible and that minimal washing is desirable after adding the label (to the extent that free-floating labels are removed satisfactorily to minimize false positives below the detection limit). Also noteworthy is that whereas the capture-antibody target equilibrium is established very quickly, equilibrium takes a much longer time to establish for the system after the addition of the detector and label. The reason for the asymmetry is simple. Due to the enormous excess of capture antibodies, any targets that dissociate almost immediately rebind, whereas there are a very few sites for the detector antibody–label complex to bind, meaning that an unbinding event between T and D or any complexes containing both components are usually final at this stage. It is interesting to note that most of the losses occur on longer time scales and that immediately after a wash everything remains bound, approaching equilibrium via a stretched exponential process. This confirms the intuitive understanding that the timing and duration of wash steps are critical considerations in designing an effective digital immunoassay workflow. The inherent nonequilibrium nature of the final measurement means that accurate comparison between assay runs requires that timings be perfectly matched, meaning that automation of the measurement process is necessary for reproducible, consistent results.

It is interesting to note that, in Figure c, the full ATDL complex makes up only a minority of the total L available in the system at any given time, with most of the L free-floating and not attached to an antibody once equilibrium is established. This is of significant consequence for downstream readout mechanisms. In the SiMoA assay model, only ATDL complexes give rise to downstream signal detection since beads need to be confined to microwells for optical readout, whereas in the nanopore assay, all sources of L leftover at the end of the washes will contribute, whether or not they are still bound to the bead. This is not necessarily a problem: after the first wash, any label that dissociates and remains in the solution was at one point bound to a target and does represent target–capture; however, care must be taken when constructing calibration curves with the downstream detector to ensure that the source of the signal, and losses thereto, is properly accounted for in the upstream assay design.

When on- and off-rates are all known for the assay components involved, the tool can also be used to estimate the time required to achieve equilibrium at each step, a critical piece of information that can inform consistent assay design without requiring complex and expensive experimentation with varied timing.

Results are consistent with the intuitive conclusion that fewer wash steps will always be better from the perspective of the total available labeled target at the end of the assay, but it should be noted carefully that the details of the readout scheme may affect what is optimal from a signal-to-background ratio perspective. For example, in the nanopore readout case,[16] because any free-floating label that remains after washing was the result of dissociation from a complex involving the capture antibody and target, it still reports accurately on the target concentration even if it is no longer complexified during readout.

DNA Nanostructure Binding

We next consider a different binding model that we recently explored experimentally:[16,21] the binding together of two star-shaped DNA nanostructures using a linker strand into a dumbbell shape, as shown in Figure a. We simulate this using a two-component system initially—stars and linkers. DNA stars have a single-stranded region at the end of a double-stranded tail that is complementary to half as of ssDNA linker. Once a linker binds to one star, the star-linker complex can bind with an additional star to form the dumbbell. We estimate kon ∼ 106 s–1 M–1,[23] and assume that binding is irreversible so that all off-rates are 0. As we will show, this simple model fails to fully capture the experimental results, and we complexify the model from this simple starting point, considering the effect first of misassembled DNA star structures that cannot properly bind as well as stars that can bind into a pseudo-dumbbell in the absence of a linker strand, and finally, including the presence of an additional population of misformed stars that can bind into duplexes with a slower on-rate. These three cases and comparison of the resulting simulations to experiment are shown in Figure .

Figure 3

(a) Comparison between an idealized simulation (blue triangles) and experimental results from He et al.[16] (black squares). (b) Illustration of the effects of nonbinding misassembled probes and of the effects of probes that can bind in the absence of a linker strand. (c) Illustration of the asymmetry introduced when a subset of probes has weaker and competitive binding kinetics with the linker. For the experimental data, the linker strand concentration ranges from 200 pM to 400 nM and shooting star probes are fixed at 20 nM, with ∼1100 single-molecule events at each concentration.

In its simplest form (Figure a,b), one intuitively expects that the fraction of dumbbells f formed should be equal to the ratio of linkers to stars x (assuming equal concentrations of both star-halves) or the inverse of that ratio, whichever is smaller, that is

This is a consequence of irreversibility of binding since when an excess of linker strands is present, the star nanostructures will get capped by the linker strand and be unable to bind further to another star, since two linker strands cannot bind together. This supports the notion that we previously presented that any practical application of these schemes must operate in the regime where linker strands are the limiting reagent, which is practically the case for any real assay.

In the experimental case, shown as the black curve in Figure b, eq only holds true for values of x ≪ 1. The peak value at x = 1 is smaller experimentally than the theoretical prediction. Below f ∼ 10–2, the dumbbell fraction reaches a minimum (noise floor) that persists even in the absence of linker strands. Finally, there exists an asymmetry showing more binding than expected for x > 1. We hypothesized that the reduced binding near x = 1 is likely due to a fraction of misassembled stars that cannot bind properly, while the noise floor is likely due to misassemblies that occur in such a way as to allow the stars to bind together in the absence of a linker. The simulation tool allows us to validate these hypotheses. If we include both types of misassembled products in the simulation, we indeed recover both the reduced maximum value of f and the false-positive behavior at both ends of the spectrum, as shown in Figure b.

It is interesting to note the discrepancy between even the modified simulation and experiment for x > 1, where simulation still underestimates the binding, and an asymmetry is present in the experimental results that is not predicted by eq . A hypothesis to explain this discrepancy is that what we are modeling as nonbinding misassemblies instead binds less strongly with linkers, and when linker strands are present in extreme excess, these weakly binding species begin to matter for the kinetics. This scenario could arise, for example, if the single-stranded tail meant to interact with the linker instead weakly binds to a misassembled single-stranded arm through some weak base-pairing that would then require a strand displacement reaction for the linker to bind properly. If this was the case, one would expect that strand displacement only to occur in the case where there is an excess of linkers leftover after binding the unaffected probes.

To test this hypothesis, we simulated this as well using an on-rate for the strand displacement component that is 10x smaller than the on-rate for the properly assembled DNA stars. The introduction of these binding elements reproduces the asymmetry observed experimentally, as the linkers will bind the misassemblies in significant quantities only in the case where there is an excess beyond that required to saturate the properly assembled probe molecules. This agreement can be seen in Figure c.

Note that the total concentration of dumbbells remains unchanged when introducing misassemblies, but is broken up between proper assemblies and misassemblies, and that the reported dumbbell fraction included in the numerator the sum of all structures that would look like a dumbbell to a nanopore (two stars bound together by a linker strand) and in the denominator the sum of the concentrations of all structures that would look like stars to the nanopore.

Because the actual rate constants are unknown and simply chosen to show qualitative effects of having different kinds of misassemblies present, the concentration breakdown is arbitrarily chosen to match the experimental data. It should be carefully noted at this point that the fact that a model accurately reproduces experimental behavior is not definitive proof that the model is an accurate representation of the physics, nor are the rate constants or concentrations chosen expected to be an accurate reproduction of the real breakdown of misassemblies. This tool is not meant to prove that a system of interaction is the underlying physical system. Rather, it allows rapid exploration of the downstream expectation if certain interactions are included so that physical intuition and hypotheses can be rapidly tested for consistency without needing to invest upfront experimental time and resources.

The second issue that arises in such first-order binding systems is that the two-step binding required to make a full DNA nanostructure greatly extends the time required to reach equilibrium when x ∼ 1, experimentally taking more than 24 h at parity (see the Supporting Information from He et al.[16]). Figure presents a comparison of the time response of the concentration of dumbbells normalized to the equilibrium concentration ratios of x = 1 and 0.1. Equilibrium takes much longer to establish at parity (red curve), even though the concentration of DNA stars in both cases is the same.

Figure 4

Comparison of time-to-equilibrium for two different ssDNA linkers to dsDNA star ratios, demonstrating that equilibrium takes much longer to establish when binary components are present in equal concentration. Note that time is log-scaled here to show the difference more clearly.

This can be understood, under the assumption of irreversible binding, by the fact that at parity the concentration of available reagents to bind is depleted as the reaction progresses, leading to progressively slower kinetics. The need for two binding events to form the positive molecular signal comes at the cost of the increased assay time, which is clearly seen in Figure and suggests that multistep bindings are to be avoided to minimize assay time when required reagents are present in equal concentrations. Future digital assays with nanopore readout should take this into consideration.

In addition to avoiding multistep bindings to minimize equilibration time, the results make clear that for DNA binding systems, the purity of the sample is of paramount importance. Misassemblies are not only able to prevent binding where desired but can introduce secondary binding modes through the exposure of partially complementary sequences. While the explanation we have provided here for the interpretation of our experimental results is not proven by the simulation, it provides a clear next step in optimizing our DNA binding assay through optimization of the purity of the DNA origami assemblies used to construct the final molecular complex.

Conclusions

We have presented a computational tool that allows for simulation of arbitrarily complex first-order binding kinetics-based assay workflows from start to finish without requiring any programming knowledge on the part of the user. We have demonstrated its utility for providing insights into assay performance and validating hypotheses through two case studies. We compared our experimental results with those of the simulation tool, finding good agreement. A detailed description of the user interface along with a basic user guide is included in the Supporting Information. The webapp is freely available at https://tcossalab.net/binding-assay/.

While methods exist to measure on- and off-rates, many ligand–receptor pairs do not have separate on- and off-rate values listed in the literature and only have the equilibrium constant for the reaction available, that is, . If this is the case for the experimental system to be simulated, an educated guess must be made as to how this quantity splits into on- and off-rates individually, as we did in this paper for one of our illustrative examples. By varying these parameters and simulating over the possible splits, the simulation tool can provide some insight into the effect of these assumptions and allow for an understanding of how these parameters affect output uncertainty. Even in the absence of exact on- and off- rate values to provide quantitative results, the tool offers a method by which assay designs are validated generally, giving clear insight into the effect on the downstream signal of wash steps and nonspecific binding of background molecules, allowing full assay workflows to be validated in minutes prior to conducting expensive and time-consuming experimental work that can take weeks. This application of our simulation tool was demonstrated in our first case study. It is also a tool with which to test ideas relating to experimental workflow failure, making it simple to test hypotheses about any discrepancies between prediction and experiment by defining potential unintended component interactions, as we demonstrated in our second case study.

Finally, the results presented so far have provided valuable insight into the bead-based digital immunoassay design generally: wash steps should be conducted as early in the assay workflow as practicable and must be carefully timed to be consistent between calibration and experiment if not allowed to equilibrate between washes. Readout scheme is of critical importance in terms of assay calibration, since depending on the nature of the scheme, it may be the case that the calibration will be sensitive to both dissociated and complexified labels. A thorough understanding of the source of the signal being measured will aid greatly in calibrating real biomolecular assays.