Transport Effects on Multiple-Component Reactions in Optical Biosensors

TL;DR

This paper develops a simplified ODE model for analyzing multi-component reactions in optical biosensors, validated through numerical simulations, to improve the measurement of kinetic rate constants in DNA repair processes.

Contribution

It introduces a reduced ODE model for complex IDE systems in biosensor experiments, enabling easier parameter estimation and analysis.

Findings

ODE model agrees with IDE simulations

Simplifies analysis of multi-component biosensor reactions

Validates the model for DNA translesion synthesis studies

Abstract

Optical biosensors are often used to measure kinetic rate constants associated with chemical reactions. Such instruments operate in the \textit{surface-volume} configuration, in which ligand molecules are convected through a fluid-filled volume over a surface to which receptors are confined. Currently, scientists are using optical biosenors to measure the kinetic rate constants associated with DNA translesion synthesis--a process critical to DNA damage repair. Biosensor experiments to study this process involve multiple interacting components on the sensor surface. This multiple-component biosensor experiment is modeled with a set of nonlinear integrodifferential equations (IDEs). It is shown that in physically relevant asymptotic limits these equations reduce to a much simpler set of Ordinary Differential Equations (ODEs). To verify the validity of our ODE approximation, a numerical method for the IDE system is developed and studied. Results from the ODE model agree with simulations of the IDE model, rendering our ODE model useful for parameter estimation.

Tables

Table 1. Table 4.1: Convergence results for the reacting species. Here Da = .01 , 1 , 10 Da .01 1 10 \mathrm{Da}=.01,\ 1,10 .

	$\mathrm{Da}\ll 1$	$\mathrm{Da}=O\left(1\right)$	$\mathrm{Da}\gg 1$
${\overline{B}}_1$	$O\left(\Delta x^{2.09}\right)$	$O\left(\Delta x^{1.33}\right)$	$O\left(\Delta x^{1.48}\right)$
${\overline{B}}_{12}$	$O\left(\Delta x^{2.04}\right)$	$O\left(\Delta x^{1.42}\right)$	$O\left(\Delta x^{1.53}\right)$
${\overline{B}}_2$	$O\left(\Delta x^{2.09}\right)$	$O\left(\Delta x^{1.33}\right)$	$O\left(\Delta x^{1.46}\right)$

Table 2. Table 5.1: Maximum difference between our injection phase ERC approximation ( 3.7 ) and our finite difference solution. All rate constants take equal to 1 in both cases.

	${\overline{B}}_1$	${\overline{B}}_{12}$	${\overline{B}}_2$	$𝒮$
$\mathrm{Da}=.1$	7.81$\times {10}^{-5}$	$5.04\times {10}^{-5}$	$7.81\times {10}^{-5}$	$1.47\times {10}^{-4}$
$\mathrm{Da}=.45$	$1.00\times {10}^{-3}$	$6.93\times {10}^{-4}$	$1.00\times {10}^{-3}$	$2.00\times {10}^{-3}$

Table 3. Table 5.2: Maximum difference between our wash phase ERC approximation ( 3.8 ) and our finite difference solution. All rate constants take equal to 1 in both cases.

	${\overline{B}}_1$	${\overline{B}}_{12}$	${\overline{B}}_2$	$𝒮$
$\mathrm{Da}=.1$	3.578$\times {10}^{-5}$	$3.40\times {10}^{-5}$	$3.58\times {10}^{-5}$	$7.37\times {10}^{-5}$
$\mathrm{Da}=.45$	$4.33\times {10}^{-4}$	$4.62\times {10}^{-4}$	$4.33\times {10}^{-4}$	$9.45\times {10}^{-4}$

Table 4. Table A.1: Dimensional parameter ranges, taken from references (de la Torre et al., 2000 ) , (Gen, 2013 ) , (Rich et al., 2008 ) , (Yarmush et al., 1996 ) .

Parameter	Rich (2008)	Yarmush	Biacore T200	Torre
$k_a$ $\left({10}^8{\mathrm{cm}}^3/\left(\mathrm{mol}\cdot s\right)\right)$	${10}^{-4}$–${10}^{-2}$	$.5$–$5\times 10$	${10}^{-5}$–$3\times 10$
$k_d$ $\left({10}^{-3}{s}^{-1}\right)$	$1$	$8.9$	${10}^{-2}$–${10}^3$
$D_1$ $\left({10}^{-7}{\mathrm{cm}}^2/s\right)$				4.0
$D_2$ $\left({10}^{-7}{\mathrm{cm}}^2/s\right)$				6.88
$H$ $\left(1\mathrm{cm}\right)$	$.05$	$.04$
$L$ $\left(1\mathrm{cm}\right)$	$2$
$W$ $\left(1\mathrm{cm}\right)$	$1.3$
$R_T$ $\left({10}^{-12}\mathrm{mol}/{\mathrm{cm}}^2\right)$	$1.11\times {10}^{-1}$–$2.33\times {10}^1$	$2.5$–4
$Q$ $\left(1\mu L/\min \right)$	$100$–$1500$		$1$–$100$
$V$ $\left(1\mathrm{cm}/s\right)$	$.153$–$2.88$	$.36$–$.6$	$.001$–$1.92$
$C_{i,u}$ $\left({10}^{-11}\mathrm{mol}/{\mathrm{cm}}^3\right)$	$2.96\times {10}^{-1}$–$2\times {10}^1$

Table 5. Table A.2: Dimensionless parameters.

Parameter	Bound
$ϵ$	$0.02$–$0.025$
Re	$8.00\times {10}^{-5}$–$0.36$
Pe	$0.12$–$523.26$
$\mathrm{Da}$	$9.29\times {10}^{-8}$–$1.49\times {10}^3$
${}_2{}^{}K_a^{}$	$2.96\times {10}^{-9}$–$3.38\times {10}^8$
${}_2{}^{1}K_a^{}$	$2.96\times {10}^{-9}$–$3.38\times {10}^8$
${}_1{}^{2}K_a^{}$	$2\times {10}^{-7}$–$5\times {10}^6$
${}_i{}^{}K_d^{}$	$1\times {10}^{-5}$–$3.38\times {10}^8$
${}_i{}^{j}K_d^{}$	$1\times {10}^{-5}$–$3.38\times {10}^8$
$D_r$	$0.58$
$F_r$	$0.01$–$39.28$

Equations

$$\overline{B}(t)=\frac{1}{x_{\max }-x_{\min }}{\int }_{x_{\min }}^{x_{\max }}B(x,t)\mathrm{d}x.$$

$$\displaystyle\mathrm{P}_{1}:E+L_{2}\xrightleftharpoons[{}_{2}k_{\mathrm{d}}]{{}_{2}k_{\mathrm{a}}}EL_{2}.$$

$$\displaystyle\mathrm{P}_{2}:E+L_{1}\xrightleftharpoons[{}_{1}k_{\mathrm{d}}]{{}_{1}k_{\mathrm{a}}}EL_{1},\ EL_{1}+L_{2}\xrightleftharpoons[{}^{1}_{2}k_{\mathrm{d}}]{{}^{1}_{2}k_{\mathrm{a}}}EL_{1}L_{2}\xrightleftharpoons[{}^{2}_{1}k_{\mathrm{a}}]{{}^{2}_{1}k_{\mathrm{d}}}EL_{2}+L_{1},\ EL_{2}\xrightharpoonup{{}_{2}k_{\mathrm{d}}}E+L_{2}.$$

$$\mathcal{S}(t)=s_1{\overline{B}}_1(t)+\left(s_1+s_2\right){\overline{B}}_{12}(t)+s_2{\overline{B}}_2(t)$$

$${\overline{B}}_i(t)=\frac{1}{x_{\max }-x_{\min }}{\int }_{x_{\min }}^{x_{\max }}{B}_i(x,t)\mathrm{dx}$$

$$\begin{array}{rl}& \tilde{x}=\frac{x}{L},\tilde{y}=\frac{y}{H},\tilde{t}={}_1{k}_{\mathrm{a}}{C}_{1,\mathrm{u}}t,{\tilde{B}}_i(x,t)=\frac{B_i(x,t)}{R_{\mathrm{T}}},{\tilde{C}}_i(x,y,t)=\frac{C_i(x,y,t)}{C_{i,\mathrm{u}}},\\ & {}_i^j{\tilde{K}}_{\mathrm{a}}=\frac{C_{i,\mathrm{u}}\cdot {}_i^j{k}_{\mathrm{a}}}{C_{1,\mathrm{u}}\cdot {}_1{k}_{\mathrm{a}}},{}_i^j{\tilde{K}}_{\mathrm{d}}=\frac{k_{\mathrm{d}}}{C_{1,\mathrm{u}}\cdot {}_1{k}_{\mathrm{a}}},{\tilde{F}}_{\mathrm{r}}={\tilde{C}}_{\mathrm{r}}{\tilde{D}}_{\mathrm{r}},{\tilde{C}}_{\mathrm{r}}=\frac{C_{1,\mathrm{u}}}{C_{2,\mathrm{u}}},{\tilde{D}}_{\mathrm{r}}=\frac{D_1}{D_2}.\end{array}$$

$$S(t)=\frac{\mathcal{S}(t)}{R_{\mathrm{T}}\cdot s_{1}}=\overline{B}_{1}(t)+\left(1+\frac{s_{2}}{s_{1}}\right)\overline{B}_{12}(t)+\frac{s_{2}}{s_{1}}\overline{B}_{2}(t){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}.}$$

$$\begin{array}{rl}\frac{\partial B_1}{\partial t}& =(1-B_{\Sigma })C_1(x,0,t)-{}_1{K}_{\mathrm{d}}{B}_1-{}_2^1{K}_{\mathrm{a}}{B}_1{C}_2(x,0,t)+{}_2^1{K}_{\mathrm{d}}{B}_{12},\end{array}$$

$$\begin{array}{rl}\frac{\partial B_{12}}{\partial t}& ={}_2^1{K}_{\mathrm{a}}{B}_1{C}_2(x,0,t)-{}_2^1{K}_{\mathrm{d}}{B}_{12}+{}_1^2{K}_{\mathrm{a}}{B}_2{C}_1(x,0,t)-{}_1^2{K}_{\mathrm{d}}{B}_{12},\end{array}$$

$$\begin{array}{rl}\frac{\partial B_2}{\partial t}& ={}_2{K}_{\mathrm{a}}(1-B_{\Sigma })C_2(x,0,t)-{}_2{K}_{\mathrm{d}}{B}_2+{}_1^2{K}_{\mathrm{d}}{B}_{12}-{}_1^2{K}_{\mathrm{a}}{B}_2{C}_1(x,0,t),\end{array}$$

$$\mathbf{B}(x,0)=0,$$

$$D_{\mathrm{r}}\frac{{\partial }^2{C}_1}{\partial {\eta }^2}=\eta \frac{\partial C_1}{\partial x},$$

$$\frac{{\partial }^2{C}_2}{\partial {\eta }^2}=\eta \frac{\partial C_2}{\partial x}.$$

$$\frac{\partial C_1}{\partial \eta }(x,0,t)=\frac{\mathrm{Da}}{F_{\mathrm{r}}}\left(\frac{\partial B_1}{\partial t}+\frac{\partial B_{12}}{\partial t}\right),$$

$$\frac{\partial C_2}{\partial \eta }(x,0,t)=\mathrm{Da}\left(\frac{\partial B_{12}}{\partial t}+\frac{\partial B_2}{\partial t}\right).$$

$$C_i(0,\eta ,t)=1,$$

$$\underset{\eta \to \infty }{\lim }{C}_i(x,\eta ,t)=1,$$

$$C_1(x,0,t)=1-\frac{D_{\mathrm{r}}^{1/3}\mathrm{Da}}{F_{\mathrm{r}}\Gamma (2/3)3^{1/3}}{\int }_0^x\left(\frac{\partial B_1}{\partial t}+\frac{\partial B_{12}}{\partial t}\right)(x-\nu ,t)\frac{\mathrm{d}\nu }{{\nu }^{2/3}},$$

$$C_2(x,0,t)=1-\frac{\mathrm{Da}}{\Gamma (2/3)3^{1/3}}{\int }_0^x\left(\frac{\partial B_{12}}{\partial t}+\frac{\partial B_2}{\partial t}\right)(x-\nu ,t)\frac{\mathrm{d}\nu }{{\nu }^{2/3}}.$$

$$\mathrm{Da}=\frac{{}_1{k}_{\mathrm{a}}{R}_{\mathrm{T}}(HL{)}^{1/3}}{(V{D}^2{)}^{1/3}}$$

$$\mathbf{B}(x,t)={}^0\mathbf{B}(t)+\mathrm{Da}{x}^{1/3}\cdot {}^1\mathbf{B}(t)+O({\mathrm{Da}}^2{x}^{2/3})$$

$$\frac{\partial \mathbf{B}}{\partial x}(x,t)=\frac{\mathrm{Da}{}^1\mathbf{B}(t)}{3x^{2/3}}+O\left(\frac{{\mathrm{Da}}^2}{x^{1/3}}\right).$$

$$\mathbf{B}(x,0)=A^{-1}\mathbf{f},$$

$$A=\left(\begin{array}{ccc}(1+\text{prescript}1K_{\mathrm{d}}+\text{prescript}12K_{\mathrm{a}})& 1-\text{prescript}12K_{\mathrm{d}}& 1\\ -\text{prescript}12K_{\mathrm{a}}& (\text{prescript}12K_{\mathrm{d}}+\text{prescript}21K_{\mathrm{d}})& -\text{prescript}21K_{\mathrm{a}}\\ \text{prescript}2K_{\mathrm{a}}& \text{prescript}2K_{\mathrm{a}}-\text{prescript}21K_{\mathrm{d}}& (\text{prescript}2K_{\mathrm{a}}+\text{prescript}2K_{\mathrm{d}}+\text{prescript}21K_{\mathrm{a}})\end{array}\right),$$

$$\mathbf{f}=\left(\begin{array}{c}1\\ 0\\ \text{prescript}2K_{\mathrm{a}}\end{array}\right).$$

$$C_i(0,\eta ,t)=0,$$

$$\underset{\eta \to \infty }{\lim }{C}_i(x,\eta ,t)=0.$$

$$C_1(x,0,t)=-\frac{D_{\mathrm{r}}^{1/3}\mathrm{Da}}{F_{\mathrm{r}}\Gamma (2/3)3^{1/3}}{\int }_0^x\left(\frac{\partial B_1}{\partial t}+\frac{\partial B_{12}}{\partial t}\right)(x-\nu ,t)\frac{\mathrm{d}\nu }{x^{2/3}},$$

$$C_2(x,0,t)=-\frac{\mathrm{Da}}{\Gamma (2/3)3^{1/3}}{\int }_0^x\left(\frac{\partial B_{12}}{\partial t}+\frac{\partial B_2}{\partial t}\right)(x-\nu ,t)\frac{\mathrm{d}\nu }{{\nu }^{2/3}}.$$

$$\overline{B_1{C}_2}=\overline{B_1\left(1-\frac{\mathrm{Da}}{3^{1/3}\Gamma (2/3)}{\int }_0^x\left(\frac{\partial B_{12}}{\partial t}+\frac{\partial B_2}{\partial t}\right)\frac{\mathrm{d}\nu }{(x-\nu {)}^{2/3}}\right)},$$

$$\mathbf{B}(x,t)={}^0\mathbf{B}(x,t)+O(\mathrm{Da}).$$

$$\frac{\mathrm{d}{}^0\mathbf{B}}{\mathrm{d}t}=-A{}^0\mathbf{B}+\mathbf{f},$$

$$\text{prescript}0\mathbf{B}(t)=A^{-1}(I-{\mathrm{e}}^{-At})\mathbf{f}$$

$$\mathbf{B}(x,t)={}^0\mathbf{B}(t)+\mathrm{Da}{x}^{1/3}\cdot {}^1\mathbf{B}(t)+O({\mathrm{Da}}^2).$$

$$\mathrm{Da}{B}_i(x,t)=\mathrm{Da}\text{prescript}0B_i(t)+O({\mathrm{Da}}^2),$$

$${\overline{B}}_1-\mathrm{Da}\overline{h}\cdot \text{prescript}0B_1\left(\frac{\mathrm{d}\text{prescript}0B_{12}}{\mathrm{d}t}+\frac{\mathrm{d}\text{prescript}0B_2}{\mathrm{d}t}\right)+O({\mathrm{Da}}^2),$$

Full text

††journal: Bull. Math. Bio.

1 Introduction

Note: this manuscript now appears in the Bulletin of Mathematical Biology, and may be found through the following reference: Evans, R.M. & Edwards, D.A. Bull Math Biol (2017) 79: 2215. https://doi.org/10.1007/s11538-017-0327-9

Kinetic rate constants associated with chemical reactions are often measured using optical biosensors. Such instruments operate in the surface-volume configuration in which ligand molecules are convected through a fluid-filled volume, over a surface to which receptors are immobilized. Ligand molecules are transported through the fluid onto the surface to bind with available receptor sites, creating bound ligand molecules at concentration $B(x,t)$. Mass changes on the surface due to ligand binding are averaged over a portion of the channel floor $[x_{\mathrm{min}},x_{\mathrm{max}}]$ to produce measurements of the form

[TABLE]

See Figure 1.1 for a schematic of one such biosensor experiment.

Measuring kinetic rate constants with optical biosensors requires an accurate model of this process, and models have been successfully proposed and progessively refined throughout the years: (Edwards, 1999, 2000, 2001, 2006, 2011; Edwards et al., 1999; Lebedev et al., 2006; Zumbrum and Edwards, 2014, 2015). Although such models are typically limited to reactions involving only a single molecule or a single step, chemists are currently using biosensor technology to measure rate constants associated with reactions involving multiple interacting components. In particular, chemists are now using biosensor experiments to elucidate how cells cope with DNA damage. Harmful DNA lesions can impair a cell’s ability to replicate DNA, and its ability to survive. One way a cell may respond to a DNA lesion is through DNA translesion synthesis (Friedberg, 2005; Lehmann et al., 2007; Plosky and Woodgate, 2004). For a description of this process we refer the interested reader to the references included herein; however, for our purposes it is sufficient to know that DNA translesion synthesis involves three interacting components: a Proliferating Cell Nuclear Antigen (PCNA) molecule, polymerase $\delta$, and polymerase $\eta$. Moreover, in order for a successful DNA translesion synthesis event to occur polymerase $\eta$ must bind with the PCNA molecule. A central question surrounding DNA translesion synthesis is whether the polymerase $\eta$ and PCNA complex forms through direct binding, or through a catalysis-type ligand switching process (Zhuang et al., 2008).

The former scenario is depicted in Figure 1.2, where we have shown polymerase $\eta$ directly binding with a PCNA molecule, i.e. the reaction:

[TABLE]

The catalysis-type ligand switching process is depticted in Figure 1.3 and stated precisely as:

[TABLE]

In (1.2b) and Figure 1.3 we have denoted polymerase $\delta$ as $L_{1}$. This process is summarized as follows: first $L_{1}$ binds with an available receptor $E$; next $L_{2}$ associates with $EL_{1}$ to create the product $EL_{1}L_{2}$; then $L_{1}$ dissociates from $EL_{1}L_{2}$, leaving $EL_{2}$; finally, $L_{2}$ dissociates from $EL_{2}$. Furthermore, in (1.2b) and Figure 1.3 the rate constants ${}_{1}k_{\mathrm{a}}$ and ${}_{1}k_{\mathrm{d}}$ denote the rates at which $L_{1}$ binds and unbinds with a receptor $E$, ${}^{j}_{i}k_{\mathrm{a}}$ denotes the rate at which ligand $L_{i}$ binds with the product $EL_{j}$, and ${}^{j}_{i}k_{\mathrm{d}}$ denotes the rate at which $L_{i}$ dissociates from the product $EL_{1}L_{2}$. In the latter two expressions the indices $i$ and $j$ can equal one or two. We shall refer to this pathway two, or simply $\mathrm{P}_{2}$ as in (1.2b).

Though Zhuang et al. provided indirect evidence of the ligand switch in (Zhuang et al., 2008), a direct demonstration of this process has not been possible with conventional techniques such as fluorescence microscopy, since such techniques introduce the possibility of modifying protein activity. Hence, scientists are using label-free optical biosensors to measure the rate constants in (1). By measuring the rate constants in (1), one could determine whether EL1L2 forms through direct binding, or the catalysis-type ligand switching process. We note that the latter manifests itself mathematically with ${}_{2}k_{\mathrm{a}}=0$, while the former with ${}_{1}^{2}k_{\mathrm{a}}={}_{1}^{2}k_{\mathrm{d}}={}_{2}^{1}k_{\mathrm{a}}={}_{2}^{1}k_{\mathrm{d}}=0$.

However, the presence of multiple intereacting components on the sensor surface complicates parameter estimation. In the present scenario there are three species $EL_{1}$, $EL_{1}L_{2}$, and $EL_{2}$ at concentrations $B_{1}(x,t)$, $B_{12}(x,t)$, and $B_{2}(x,t)$, and since optical biosensors typically measure only mass changes at the surface, lumped measurements of the form

[TABLE]

are produced. In (1.3)

[TABLE]

denotes the average reacting species concentration, for $i=1,\ 12,\ 2$, and $s_{i}$ denotes the molecular weight of $L_{i}$. The lumped signal (1.3) raises uniqueness concerns, since more than one set of rate constants may possibly correspond to the same signal (1.3). Fortunately, through varying the uniform in-flow concentrations of the ligands, $C_{1}(0,y,t)=C_{1,\mathrm{u}}$ and $C_{2}(0,y,t)=C_{2,\mathrm{u}}$, one may resolve this ill-posedness in certain physically relevant scenarios (Evans, R. M. and Edwards, D. A. and Li, W., submitted). This approach to identifying the correct set of rate constants in the presence of ambiguous data is related to the “global analysis” technique in biological literature (Karlsson and Fält, 1997; Morton et al., 1995).

The presence of multiple interacting species and the lumped signal (1.3) complicate parameter estimation even for systems accurately descibed by the well-stirred kinetics approximation. However in (Edwards, 1999), Edwards has shown that transport dynamics affect ligand binding in a thin boundary layer near the sensor surface. Hence, we begin in Section 2 by summarizing the relevant boundary layer equations, which take the form of a set of nonlinear Integrodifferential Equations (IDEs). In Section 3, it is shown that in experimentally relevant asymptotic limits our IDE model reduces to a much simpler set of Ordinary Differential Equations (ODEs) which can be used for parameter estimation. To verify the accuracy of our ODE approximation, a numerical method is developed in Subsection 4.1. Convergence properties are examined in Subsection 4.2, and in Section 5 the accuracy of our ODE approximation is verified by comparing results of our ODE model with results from our numerical method described in Section 4. Conclusions and plans for future work are discussed in Section 6.

2 Governing Equations

For our purposes, biosensor experiments are partitioned into two phases: an injection phase, and a wash phase. During the injection phase $L_{1}$ and $L_{2}$ are injected into the biosensor via a buffer fluid at the uniform concentrations $C_{1}(x,y,0)=C_{1,\mathrm{u}}$ and $C_{2}(x,y,t)=C_{2,\mathrm{u}}$. Injection continues until the signal (1.3) reaches a steady-state, at which point the biosensor is washed with the buffer fluid–this is the wash phase of the experiment. Only pure buffer is flowing through biosensor during the wash phase, not buffer containing ligand molecules. This causes all bound ligand molecules at the surface to dissociate and flow out of the biosensor, thereby preparing the device for another experiment. We first summarize the governing equations for the injection phase.

2.1 Injection Phase

To present our governing equations we introduce the dimensionless variables:

[TABLE]

We have scaled the spatial variables with the instrument’s dimensions, time with the association rate of $L_{1}$ onto an empty receptor, the bound ligand concentrations $B_{i}$ with the initial free receptor concentration, and the unbound ligand concentrations with their respective uniform inflow concentrations. The rate constants ${}_{i}^{j}\widetilde{K}_{\mathrm{a}}$ and ${}_{i}^{j}\widetilde{K}_{\mathrm{d}}$ are the dimensionless analogs of ${}_{i}^{j}k_{\mathrm{a}}$ and ${}_{i}^{j}k_{\mathrm{d}}$. In the latter expressions the index $i=1,\ 2$, whereas $j=1,\ 2,$ or can be blank. Furthermore, $\widetilde{F}_{\mathrm{r}}$ measures the diffusion strength of each reacting species, as characterized by the product of the input concentrations and the diffusion coefficients. Henceforth, we shall drop the tildes on our dimensionless variables for simplicity. In particular, we denote the dimensionless sensogram reading as

[TABLE]

Moreover, we may use (1.4) to denote the dimensionless average concentration, as it is of the same form in both the dimensionless and dimensional contexts.

Applying the law of mass action to (1) gives the kinetics equations:

[TABLE]

which hold on the reacting surface when $y=0$ and $x\in[0,1]$. In (4.12d), $\mathbf{B}=(B_{1},\ B_{12},\ B_{2})^{T}$ is a vector in $\mathbb{R}^{3}$ whose components contain the three bound state concentrations. In addition, the terms in equations (2.3a)–(2.3c) have been ordered in accordance with Figures 1.2 and 1.3.

Edwards has shown (Edwards, 1999) that transport effects dominate in a thin boundary layer near the reacting surface where diffusion and convection balance. Hence the governing equations for $C_{i}$ are

[TABLE]

Since $C_{1}$ is used up in the production of $B_{1}$ and $B_{12}$, and $C_{2}$ is used up in the production of $B_{12}$ and $B_{2}$, we have the diffusive flux conditions:

[TABLE]

Equations (2.4a)–(2.4d) reflect the fact that in the boundary layer $C$ is in a quasi-steady-state where change is driven solely by the surface reactions (2.4c)–(2.4d). Then, given the inflow and matching conditions

[TABLE]

the solution to (2.1) is given by

[TABLE]

See (Edwards, 1999) for details of a similar calculation. During the injection phase, the bound state concentration is then governed by (2.3) using (2.5).

In (2.4c)–(2.4d) and (2.5)

[TABLE]

is the Damköhler number–a key dimensionless parameter which measures the speed of reaction relative to the transport into the surface. In the experimentally relevant parameter regime of $\mathrm{Da}\ll 1$, the time scale for transport into the surface is much faster than the time scale for reaction. In this case there is a only a weak coupling between the two processes, and (2.5) shows that the unbound concentration at the surface is only a perturbation away from uniform inlet concentration. When $\mathrm{Da}\to 0$ in (2.3) using (2.5), one recovers the well-stirred approximation in which transport into the surface completely decouples from reaction.

On the other hand, when $\mathrm{Da}=O(1)$ the two processes occur on the same time scale, and ligand depletion effects become more evident. This is a phenomenon in which ligand molecules are transported into the surface to bind with receptor sites upstream, before they bind with receptor sites downstream. Mathematically, this is reflected in the convolution integrals in (2.5). When $x\ll 1$ the convolution integral influences the unbound concentration at the surface less than when $x$ is larger.

A sample space-time curve for each of the reacting species concentrations $B_{i}(x,t)$ is depicted in Figure 2.1, where we have shown the results of our numerical simulations described in Section 4.

The $x$-axis represents the sensor, and $t$-axis represents time. Injection begins at $t=0$, and ligand molecules bind with receptor sites as they are transported into the surface. Binding proceeds as the injection continues; finally each of the concentrations achieve a chemical equilibrium in which there is a balance between association and dissociation. Observe the spatial heterogeneity present in each of the bound state concentrations–the reaction proceeds faster near the inlet at $x=0$ than the rest of the surface. This is precisely the ligand depletion phenomenon described in the above paragraph, and is particularly evident in the surface plot of $B_{12}$. This is because in this simulation we have taken all of the rate constants equal to one, and either $EL_{1}$ or $EL_{2}$ must be present in order for $EL_{1}L_{2}$ to form. Thus, in this case $EL_{1}L_{2}$ experiences effectively twice the ligand depletion of the other reacting species.

Furthermore, one may notice an apparent discontinuity in each of the surface plots depicted in Figure 2.1–this reflects the weakly singular nature of the functions which we are attempting to approximate. When $x\ll 1$, one may show $\mathbf{B}$ has the perturbation expansion

[TABLE]

(this is simply (3.4) for $x\ll 1$). It therefore follows that

[TABLE]

Hence, although the function $\mathbf{B}$ is well-defined and continuous near $x=0$, it has a vertical tangent at $x=0$. The weakly-singular nature of $\mathbf{B}$ is magnified since $\mathrm{Da}=2$. To resolve this region, one may think to adaptively change $\Delta x$ with the magnitude of $\partial\mathbf{B}/\partial x$. However, because the sensogram reading $S(t)$ is computed over the region $[x_{\mathrm{min}},x_{\mathrm{max}}]$, we are not concerned with resolving this region and a uniform step size is sufficient. Moreover, our convergence results in Subsection 4.2 demonstrate that a lack of resolution at $x=0$ does not affect our results in the region of interest $[x_{\mathrm{min}},x_{\mathrm{max}}]$.

2.2 Wash Phase

We now summarize the relevant equations for the wash phase. In practice the injection phase is run until the bound state concentration reaches a steady-state (Rich and Myszka, 2009). This implies that because the bound ligand concentration evolves on a much slower time scale than the unbound ligand concentration (Edwards, 1999), the unbound ligand concentration will have also reached steady-state by the time the wash phase begins. In particular, the unbound concentration on the surface will be uniform by the time the wash phase starts–i.e., $C_{i}(x,0,0)=1$. Thus, the kinetics equations are given by (2.3), with (4.12d) replaced by the steady solution to (2.3) during the injection phase:

[TABLE]

Equations similar to (2.1) hold:

[TABLE]

[TABLE]

Thus, during the wash phase the bound state evolution is governed by the (2.3a)–(2.3c), (2.9), and (2.11).

3 Effective Rate Constant Approximation

During both phases of the experiment, the bound state concentration $\mathbf{B}(x,t)$ obeys a nonlinear set of IDEs which is hopeless to solve in closed form. However, we are ultimately interested in the average concentration $\overline{\mathbf{B}}(t)$, rather than the spatially-dependent function $\mathbf{B}(x,t)$, since from $\overline{\mathbf{B}}(t)$ we can construct the sensogram signal (2.2) (the quantity of interest). Thus, we seek to find an approximation to $\overline{\mathbf{B}}(t)$, and begin by finding one during the injection phase. We first average each side of (2.3), with $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ given by (2.5), in the sense of (1.4). Immediately, we are confronted with terms such as

[TABLE]

on the right hand side of (2.3a). In the experimentally relevant case of small $\mathrm{Da}$, we are motivated to expand ${\bf B}(x,t)$ in a perturbation series:

[TABLE]

In this limit, the leading order of (2.5) is just $C_{i}=1$. Using this result in (2.3), we have that the governing equation for ${}^{0}{\bf B}$ is independent of $x$:

[TABLE]

where $A$ is given by (2.9b) and $\mathbf{f}$ by (2.9c). Hence the leading-order approximation

[TABLE]

is independent of space. Substituting (3.3) into (2.5), the time-dependent terms may be factored out of the integrand, leaving the spatial dependence of $C_{j}$ varying as $x^{1/3}$. This is the only spatial variation in (2.3) at $O(\mathrm{Da})$; hence we may write

[TABLE]

As a result of (3.4) we have the relation

[TABLE]

which may be used to show the right hand side of (3.1) is equal to

[TABLE]

We then average (3.5), and use the resulting relation in (3.6) to show the right hand side of (3.1) reduces to:

[TABLE]

In this manner, we can derive a set of nonlinear ODEs for $\overline{\mathbf{B}}(t)$ of the form:

[TABLE]

We have also derived a set of ERC equations for the wash phase, they take take the form:

[TABLE]

Following (Edwards and Jackson, 2002), we refer to the Ordinary Differential Equation (ODE) systems (3.7) and (3.8) as our Effective Rate Constant (ERC) Equations. A significant advantage of our ERC equations is that these ODEs are far easier to solve numerically than their IDE counterparts. To solve (3.7) or (3.8), one may simply apply their linear multistage or multistep formula of choice. This feature renders our ERC equations attractive for data analysis, since they can be readily implemented into a regression algorithm when attempting to determine the rate constants associated with the reactions (1). Since experimental data is still forthcoming, we do not employ a regression algorithm to fit the rate constants in (3.7) and (3.8) to biosensor data. Synthetic data for the kinetic rate constants was used in our numerical simulations.

Solutions of our ERC equations for different parameter values are depicted in Figure 3.1.

First consider the solutions depicted on the left. Here the injection phase (3.7) has been run from $t=0$ to $t=5$ and the wash phase (3.8) has been run from $t=5$ to $t=10$. Furthermore, all rate constants were taken equal to one and the Damköhler number was $\mathrm{Da}=0.1$. During the injection phase it is seen that $\overline{B}_{1}$ and $\overline{B}_{2}$ reach equilibrium after approximately one second, while $\overline{B}_{12}$ takes approximately two seconds. This is not a surprise: we are injecting equal amounts of both ligands, all the rate constants are the same, and either $EL_{1}$ or $EL_{2}$ must already be present in order for $EL_{1}L_{2}$ to form. The equality of the rate constants is also the reason why all three species attain the same steady-state. Mathematically, the steady-state of $\mathbf{B}$ during the injection phase is given by (2.9a), and one can readily verify that $A^{-1}\mathbf{f}=(1/4,\ 1/4,\ 1/4)^{T}$ when all of the rate constants are equal to one. Physically, each of the species ultimately achieves the same balance between association and dissociation. Furthermore, the fact that all of the rate constants are the same is the reason why $\overline{B}_{12}$ decays to zero faster than the other two species: $EL_{1}L_{2}$ transitions to either $EL_{1}$ or $EL_{2}$ at the same rate as the latter two species transition into an empty receptor $E$.

Now consider the solutions depicted on the right in Figure 3.1. As with the previous case the injection phase has been run from $t=0$ to $t=5$ and the wash phase has been run from $t=5$ to $t=10$. However this time, all the rate constants have been taken equal to one except ${}_{2}^{1}K_{\mathrm{a}}$, which was taken equal to ${}_{2}^{1}K_{\mathrm{a}}=10$. During the injection phase it is seen that $\overline{B}_{1}$ quickly reaches a local maximum, and then decreases to steady-state. Since ${}_{2}^{1}K_{\mathrm{a}}$ is an order of magnitude larger than the other rate constants, after a short period of time $L_{2}$ molecules bind with $EL_{1}$ at a faster rate than $L_{1}$ molecules bind with empty receptor sites. This results in the chemical equilibrium between $EL_{1}$ and $EL_{1}L_{2}$ depicted on the right in Figure 3.1. From these observations it is clear why the steady-state value of $\overline{B}_{12}$ is larger than the previous case. However, it may be counterintuitive to observe that $\overline{B}_{2}$ reaches a larger steady-state value in the solutions depicted on the right than the solutions depicted to the left. Although one may think the vast majority $L_{2}$ molecules should be used in forming $EL_{1}L_{2}$, the increase in $EL_{1}L_{2}$ also increases the concentration of empty receptor sites. The continuous injection of $L_{2}$ therefore drives the average concentration $\overline{B}_{2}$ to a larger steady-state value. During the wash phase, it is seen that $\overline{B}_{1}$ reaches a global maximum after approximately $t\approx 5.75$ seconds. The increase in $\overline{B}_{1}$ during the wash phase is a direct consequence of $L_{2}$ molecules dissociating from $EL_{1}L_{2}$. Since only pure buffer is flowing through the biosensor during the wash phase, it is seen in Figure 3.1 that each of the average concentrations $\overline{B}_{i}$ decay to zero.

4 Numerics

To verify the $O(\mathrm{Da}^{2})$ accuracy of our ERC approximation derived in Section 3, we now develop a numerical approximation to the IDE system (2.3), where $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ are given by (2.5). We focus on the injection phase, since the wash phase is similar. Our approach is based on the numerical method described in (Edwards and Jackson, 2002). Semi-implicit methods have been previously used with great success to solve reaction-diffusion equations (Nie et al., 2006), as they are typically robust, efficient, and accurate. Similarly, in our problem we exploit the structure of the integrodifferential operator, which naturally suggests a semi-implicit method in time. Moreover, since our method is semi-implicit in time we avoid the expense and complication of solving a nonlinear system at each time step. Convergence properties and remarks concerning stability, are discussed in Subsection 4.2; however, we first turn our attention to deriving our numerical method in Subsection 4.1.

4.1 Semi-implicit finite difference algorithm

We discretize the spatial interval $[0,1]$ by choosing $N+1$ equally spaced discretization nodes $x_{i}=i\Delta x$, for $i=0,\ \ldots,N$, and discretize time by setting $t_{n}=n\Delta t$, for $n=0,\ \ldots$. Having chosen our discretization nodes and time steps, we seek to discretize (2.3), where $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ are given by (2.5). Note that this requires discretizing both the time derivatives and the convolution integrals; we first turn our attention to the latter, and focus on $C_{1}(x,0,t)$. We would like to apply the trapezoidal rule to spatially discretize (2.5a), however the integrand of $C_{1}(x,0,t)$ is singular when $\nu=0$. To handle the singularity we subtract and add

[TABLE]

from the integrand. Doing so yields

[TABLE]

where we have used the fact that (4.1) is independent of $\nu$. Then choosing a discretization node $x=x_{i}$ and a time step $t=t_{n}$, we apply the trapezoidal rule to (4.2) to obtain

[TABLE]

when $x_{i}>0$; simply evaluating (4.2) at $x=x_{0}$ gives $C(x_{0},0,t_{n})=1$. The first term in the sum is zero, because in a similar manner to Appendix B of (Zumbrum, 2013) we have

[TABLE]

which implies

[TABLE]

for $k=1,\ 12,$ or $2$. The last equality follows since we expect $\partial B_{i}/\partial t$ to be sufficiently regular for fixed $x>0$. The expansion (3.4) shows that this is certainly true when $\mathrm{Da}\ll 1$, however when $\mathrm{Da}=O(1)$ or larger the nonlinearity in (2.3) renders any analytic approximation to $B_{i}$ beyond reach. Our results in Subsection 4.2 show that our method indeed converges when $\mathrm{Da}=O(1)$ or larger.

We now turn our attention to discretizing the time derivatives. We denote our approximation to $B_{j}(x_{i},t_{n})$ by

[TABLE]

and approximate the time derivatives through the formula

[TABLE]

Our approximation (4.7) holds for all reacting species $j=1,\ 12,\ 2$, each of our discretization nodes $x_{i}$, and each time step $t_{n}$. As we shall show below, we treat $\Delta B^{j}_{i,n}$ as separate variable used to update $B^{j}_{i,n}$ at each iteration of our algorithm.

With our time derivatives discretized as (4.7), the fully-discretized version of $C_{1}(x,0,t)$ is given by substituting (4.7) into (4.3):

[TABLE]

We enforce the initial condition (4.12d) at our $N+1$ discretization nodes through the condition $B^{j}_{i,0}=0$ for $j=1,\ 12,\ 2$, and $i=1,\ldots,N$. Observe that our method (4.8) is semi-implicit rather than fully-implicit. This renders (4.8) linear in $\Delta B^{j}_{i,n+1}$, and as a result we can write

[TABLE]

which is analogous to a second-order Adams-Bashforth formula.

In addition, we chose a method that is implicit in $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ also due to the form of the convolution integrals. From (2.5) we see $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ depend on $B_{j}(\nu,t)$ only for $\nu\leq x$. Thus by choosing a method that is implicit in $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$, we are able to use the updated values of $B_{j}(x,t)$ in the convolution integrals by first computing the solution at $x=0$, and marching our way downstream at each time step.

To make this notion more precise we note that in (4.9a) the matrix $M^{-1}_{i,n}(\mathbf{B}_{i,n})$ depends only upon $\mathbf{B}_{i,n}$, however because of the convolution integrals $C^{1}_{i,n+1}$ and $C^{2}_{i,n+1}$, the matrix $A_{i,n+1}$ and vector $\mathbf{f}_{i,n+1}$ depend upon $\mathbf{B}_{l,n+1}$ for $l<i$. Thus, at each time step $n+1$ we first determine $\mathbf{B}_{0,n+1}$. Next, we increment $i$ and use the value of $\mathbf{B}_{0,n+1}$ in (4.9b) to determine $\mathbf{B}_{1,N+1}$. We proceed by iteratively marching our way downstream from $x_{2}$ to $x_{N}$ to determine $\mathbf{B}_{2,n+1},\ \ldots,\ \mathbf{B}_{N,n+1}$. Intuitively, the updated information from the convolution integral flows downstream from left to right at each time step. We may repeat this procedure for as many time steps as we wish. In addition, we remark that the formula (4.8) was initialized with one step of Euler’s method.

Furthermore, with our finite difference approximation to $\mathbf{B}(x,t)$, we can determine the average quantity

[TABLE]

with the trapezoidal rule

[TABLE]

In (4.11), the indices $i=m$ and $i=M$ correspond to $x_{\mathrm{min}}=m\Delta x$ and $x_{\mathrm{max}}=M\Delta x$. Our nodes were chosen to align with $x_{\mathrm{min}}$ and $x_{\mathrm{max}}$ to avoid interpolation error.

4.2 Convergence study

4.2.1 Spatial Convergence

We now examine the spatial rate of convergence of our numerical method. Since from $\overline{\mathbf{B}}$ we can compute the quantity of interest (2.2), we derive estimates for the rate at which our numerical approximation converges to $\overline{\mathbf{B}}$. Furthermore, because the system (2.3), (2.5) is nonlinear, our analysis will focus on the experimentally relevant case of $\mathrm{Da}\ll 1$. In addition, we will derive estimates only for the injection phase of the experiment, since the wash phase is similar.

To proceed, we consider the average variant of (2.3), (2.5). Averaging (2.3) in the sense of (1.4) gives:

[TABLE]

As in Subsection 4.1, we handle the singularity in (2.5a) by adding and subtracting (4.1) from the integrand of (2.5a) to write $C_{1}(x,0,t)$ as in (4.2). The unbound ligand concentration $C_{2}(x,0,t)$ has a representation analogous to (4.2). In the following analysis we limit our attention to (4.12a), since the analysis for equations (4.12b)–(4.12c) is nearly identical.

We proceed by anaylzing each of the terms in (4.12a):

[TABLE]

Upon inspecting (4.2) and using linearity of the averaging operator, we see that three terms contribute to (4.13c):

[TABLE]

In a similar manner, (4.13d) and (4.13e) each imply that we incur error from the terms:

[TABLE]

Let us denote the trapezoidal rule of a function $f(x)$ over the interval $[a,b]$ by $\mathcal{T}(f(x),[a,b])$. Then since $\mathcal{T}(1,[x_{\mathrm{min}},x_{\mathrm{max}}])$ is exact, the term (4.14a) does not contribute to the spatial discritization error.

Next we decompose the expansion (3.4) into its individual components to obtain

[TABLE]

for $j=1,\ 12,\ 2$. Substituting (4.16) into (4.13b), (4.13a), (4.15a) (4.15d), and using the fact that $\mathcal{T}(x^{1/3},[x_{\mathrm{min}},x_{\mathrm{max}}])$ converges at a rate $O(\Delta x^{2})$ shows that each of these terms converge at a rate of $O(\mathrm{Da}\Delta x^{2})$. Similarly, one can substitute (4.16) into (4.14c), (4.15c), and (4.15f), and use the fact that $\mathcal{T}(x^{1/3},[x_{\mathrm{min}},x_{\mathrm{max}}])$ converges at a rate of $O(\Delta x^{2})$, to show that each of these terms converge at a rate of of $O(\mathrm{Da}\Delta x^{2})$.

It remains to determine the error associated with (4.14b), (4.15b), and (4.15e), so we turn our attention to (4.14b) and substitute (4.16) into (4.14b) to obtain

[TABLE]

where we have used the definition of our averaging operator (1.4). In writing (4.17), neglected higher-order terms which do not contribute to the leading-order spatial discretization error. Since the coefficient of the integral in (4.17) is a function of time alone, this coefficient does not contribute to the leading-order spatial discretization error and we neglect it in our analysis. Hence, to compute the spatial discretization error associated with (4.17), we calculate the error associated with applying the trapezoidal rule to the double integral

[TABLE]

Treating the inner integral as a function of $x$ we define

[TABLE]

whose closed form is given by

[TABLE]

Towards applying the trapezoidal rule to (4.18), we first note $\mathcal{T}(f,[0,x_{i}])$ converges at a rate of $O(\Delta x^{4/3})$. This is seen by first rewriting (4.19) as

[TABLE]

The term on the right converges at a rate of $O(\Delta x^{2})$, and the term on the left converges at a rate of $O(\Delta x^{4/3})$, which follows from expanding $(x-\nu)^{1/3}$ about $\nu=0$, and using the definition of the trapezoidal rule.

Applying the trapezoidal rule to (4.18) then gives

[TABLE]

where we have let $x_{\min}=x_{m}=m\Delta x$, and $x_{\max}=x_{M}=M\Delta x$. Since $\mathcal{T}(f,[0,x_{i}])$ converges at a rate of $O(\Delta x^{4/3})$, the right hand side of the above is

[TABLE]

To compute our results in Section 5, we took $x_{\mathrm{min}}=0.2,\ x_{\mathrm{max}}=0.8$, in accordance with the literature (Edwards and Jackson, 2002). Hence, in the above sum there are approximately $0.6N=0.6\Delta x^{-1}$ terms on the order of $O(\mathrm{Da}^{2}\Delta x^{7/3})$, and the above sum reduces to

[TABLE]

The dominant error in (4.24) is $O(\mathrm{Da}^{2}\Delta x^{4/3})$, thus the spatial discretization error associated with (4.14b) is $O(\mathrm{Da}^{2}\Delta x^{4/3})$. When measuring convergence we used values of $x_{\mathrm{min}}=.25,\ x_{\mathrm{max}}=.75$ to facilitate progressive grid refinement; however, it is clear that these values of $x_{\mathrm{min}}$ and $x_{\mathrm{max}}$ do not change our argument. A similar argument shows the spatial discretization error associated with the nonlinear terms (4.15b) and (4.15f) is $O(\mathrm{Da}^{2}\Delta x^{4/3})$.

We have depicted our spatial convergence measurements for $\overline{B}_{1}$ in Figure 4.1 and tabulated them in Table 4.1. To obtain these results, we first computed a reference solution, with $\Delta x=\Delta t=1/512$. We then created a series of test solutions with mesh width $\Delta x=1/2^{j}$, for $j=2,\ldots,7$, keeping $\Delta t=1/512$ constant. Next, we computed $\overline{\mathbf{B}}$ by averaging our reference solution and test solutions at each time step with the trapezoidal rule as in (4.11). We then computed the error between each test solution and the reference solution by taking the maximum difference of the two over all time steps.

From our results, we see that our method converges at a rate of $O(\Delta x^{2})$ when $\mathrm{Da}\ll 1$, $O(\Delta x^{4/3})$ when $\mathrm{Da}=O(1)$, and $O(\Delta x^{3/2})$ when $\mathrm{Da}\gg 1$. The reduction in convergence when $\mathrm{Da}$ increases from small to moderate may be attributed to the $O(\mathrm{Da}^{2}\Delta x^{4/3})$ contributions from (4.14b), (4.15b), and (4.15e). There are two competing magnitudes of error in (4.12a): one of $O(\mathrm{Da}\ \Delta x^{2})$ (from terms (4.13b), (4.13a), (4.14c), (4.15a), (4.15c), (4.15d), and (4.15f)), and one of $O(\mathrm{Da}^{2}\Delta x^{4/3})$ (from the integral terms (4.14b), (4.15b), and (4.15e)). When $\mathrm{Da}^{2}\Delta x^{4/3}<\mathrm{Da}\ \Delta x^{2}$, or $\mathrm{Da}<\Delta x^{2/3}$, the former is larger. Conversely, when $\Delta x^{2/3}<\mathrm{Da}$ the latter is larger.

When $\mathrm{Da}\gg 1$, the bound state evolves on a longer time scale of the form (Edwards, 1999)

[TABLE]

In this case, the characteristic time scale for reaction for reaction is much faster than the characteristic time scale for transport into the surface, and one typically refers to $\mathrm{Da}\gg 1$ as the transport-limited regime. Substituting (4.25) into (2.3), (2.5), one may find the leading-order approximation to the resulting system for $\mathrm{Da}\gg 1$ by neglecting the left hand side of (2.3a)–(2.3c). Doing so one finds that even a leading-order approximation to (2.3a)–(2.3c) is nonlinear, rendering any error estimates in the transport-limited regime beyond reach. Nonetheless, our results in Figure 4.1 and Table 4.1 show that convergence is not an issue when $\mathrm{Da}\gg 1$.

4.2.2 Temporal Convergence

Since our time stepping method (4.9b) is analogous to a second-order Adams-Bashforth formula, we expect our method to achieve second-order accuracy in time. Figure 4.2 shows that this is indeed the case when $\mathrm{Da}=.01$, and the rate constants are $\prescript{}{1}{K_{\mathrm{d}}}=1,\ \prescript{}{2}{K_{\mathrm{a}}}=1,\ \prescript{}{2}{K_{\mathrm{d}}}=1,\ \prescript{2}{1}{K_{\mathrm{a}}}=1/2,\ \prescript{2}{1}{K_{\mathrm{d}}}=2,\ \prescript{1}{2}{K_{\mathrm{a}}}=2,$ and $\prescript{1}{2}{K_{\mathrm{d}}}=1/2$ (as in Subsection 4.2.1 when measuring spatial convergence). Temporal convergence was measured in an analogous manner to spatial convergence.

However, we note that measuring temporal convergence when $\mathrm{Da}=O(1)$ is computationally prohibitive, since spatial convergence is $O(\mathrm{Da}^{2}\Delta x^{4/3})$ in this case, so in order for the spatial and temporal errors to balance one must have $O(\mathrm{Da}^{2}\Delta x^{4/3})=O(\Delta t^{2})$ or $\Delta x=\Delta t^{3/2}$. Nonetheless, our results from Section 5 demonstrate that our finite difference approximation agrees with our ERC approximation for a wide parameter range, so we are not concerned with temporal convergence of our method when $\mathrm{Da}=O(1)$ or larger.

4.2.3 Stability Remarks

We now make brief remarks concerning the stability of our method. Recall from Subsection 4.1 we first determine the value of $B_{i}(x,t)$ upstream at $x=0$, and iteratively march our way downstream to $x=1$ at each time step. Therefore, we expect any instabilities at $x=0$ to propagate downstream. Requiring that there no instabilities at $x=0$ is equivalent to asking that our time stepping method (4.9b) is stable for the ODE system found by replacing $C_{1}(x,0,t)$ and $C_{2}(x,0,t)$ with the constant function 1 in (2.3). Though we do not have precise stability estimates for this system, numerical experimentation has shown that our time steps need to be sufficiently small in order to ensure that our numerical approximation is well behaved.

5 Effective Rate Constant Approximation Verification

With our numerical method in hand, we are now in a position to verify the accuracy of our ERC approximations (3.7) and (3.8). We tested the accuracy of our ERC equations when $\mathrm{Da}=0.1$, and $\mathrm{Da}=0.45$; the results are below in Figure 5.1 and Tables 5.1–5.2.

From these results, it is evident that our ERC equations accurately characterize $\overline{\mathbf{B}}$ and the sensogram reading (2.2) not only for small $\mathrm{Da}$, but for moderate $\mathrm{Da}$ as well. Motivated by (Edwards and Jackson, 2002), we ran a series of simulations for different values of $\mathrm{Da}$, ranging from $\mathrm{Da}\approx 0.02$ to $\mathrm{Da}=150.$ We measured the maximum absolute error for each value of $\mathrm{Da}$, and created the curves shown in Figure 5.2. The error starts off small as expected, and increases at rates which compares favorably with our $O(\mathrm{Da}^{2})$ prediction, and finally reaches an asymptote corresponding to roughly two percent absolute error. Thus, although our ERC approximations (3.7) and (3.8) are formally valid for only small values of $\mathrm{Da}$, their solutions agree with our finite difference approximation for moderate and large values of $\mathrm{Da}$.

6 Conclusions

Scientists are attempting to determine whether the polymerase $\eta$ and PCNA complex (denoted $EL_{2}$ throughout) which results from DNA translesion synthesis forms through direct binding (1.2a), or through a catalysis-type ligand switching process (1.2b). Since fluorescent labeling techniques may modify protein behavior, label-free optical biosensor experiments are used. Interpreting experimental data relies on a mathematical model, and modeling multiple-component biosensor experiments results in a complicated and unwieldy set of equations. We have shown that in experimentally relevant limits this model reduces to a much simpler set of ODEs (our ERC equations), which can be used to fit rate constants using biosensor data. In contrast with the standard well-stirred kinetics approximation, our ERC equations accurately characterize binding when mass transport effects are significant. This renders our ERC equations a flexible tool for estimating the rate constants in (1). In turn, estimates for the rate constants in (1) will reveal whether the polymerase $\eta$ and PCNA complex forms via direct binding (1.2a), or the catalysis-type ligand switching process (1.2b).

Furthermore, the consideration of both direct binding (1.2a) and the ligand switching process (1.2b) has several mathematical and physical consequences. First, due the form of (1), the species are directly coupled through the kinetics equations. This is true even in the well-stirred limit in which $\mathrm{Da}\to 0$, and (2.5) reduces to $C_{1}(x,0,t)=C_{2}(x,0,t)=1$. However, transport effects manifested in (2.5) nonlinearly couple the reacting species. So we see in Figure 2.1 that there is a more pronounced depletion region in $B_{12}$ than in either of the other two species. Physically, this is a consequence of the fact that either $EL_{1}$ or $EL_{2}$ must be present in order for $EL_{1}L_{2}$ to form, thus the latter is affected by depletion of the former two species. Additionally, the multiple-component reactions (1) alter the form of the sensogram reading to the lumped signal (2.2), thereby complicating parameter estimation.

In addition to establishing a firm foundation for studying the inverse problem of estimating the rate constants in (1), the present work also opens the door for future work on modeling and simulating multiple-component biosensor experiments. This includes considering other physical effects like cross-diffusion, or steric hinderance; and comparing the finite difference method described herein to the method of lines algorithm discussed in (Zumbrum, 2013).

Appendix A Parameter Values

Parameter values from the literature are tabulated below.

The variables $W,\ Q$, represent the dimensional width, and flow rate; the other dimensional variables are as in Section 2. The flow rate is related to the velocity through the formula (Edwards, 2011)

[TABLE]

Using the dimensional values above, we calculated the following extremal bounds on the dimensionless variables.

Here $\epsilon=H/L$ is the aspect ratio, and $\mathrm{Re}=VH^{2}/(\nu L)$ is the appropriate Reynolds number associated with our system.

The authors wish to emphasize that the bounds in Table A.2 are naïve extremal bounds calculated by using minimum and maximum values for the dimensional parameters in Table A.1. In particular, the values for the dimensionless rate constants in Table A.2 are not estimates of their true values; they are minimum and maximum values calculated using combinations of extremal values for the parameters in Table A.2. A large variation in the dimensionless rate constants is highly unlikely, since this scenario corresponds to one in which one of the association rate constants is very large, and another association rate constant very small. We would also like to note that a large variation in some of the parameters, such as the kinetic rate constants or $\mathrm{Da}$, would necessitate very small values for either or both of $\Delta t$ and $\Delta x$ in our numerical method.

Furthermore, one may be concerned about the upper bound on the Reynolds number, the lower bound on the Péclet number, and the upper bound on the Damköhler number. All of these extremal bounds were calculated using a flow rate of $1\ \mu\mathrm{L}/\mathrm{min}$–the slowest flow rate possible on the BIAcore T200 (Gen, 2013). Even with the fastest reactions, one can still design experiments to minimize transport effects by increasing the flow rate $Q$ (thus the velocity), decreasing the initial empty receptor concentration $R_{\mathrm{T}}$, and decreasing the ligand inflow concentrations $C_{1,\mathrm{u}}$ and $C_{2,\mathrm{u}}$. In the case of the fastest reaction ${}_{1}k_{\mathrm{a}}=3\times 10^{9}\ \mathrm{cm}^{3}/(\mathrm{mol}\cdot\mathrm{s})$, we can take: $Q=390\ \mu\mathrm{L}/\mathrm{min},\ V=.75\ \mathrm{cm}/\mathrm{s},\ R_{\mathrm{T}}=7.76\times 10^{-13}\ \mathrm{mol}/\mathrm{cm}^{2},\ C_{1,\mathrm{u}}=C_{2,\mathrm{u}}=2.96\times 10^{-12}\ \mathrm{mol}/\mathrm{cm}^{3}.$ These choices yield the dimensionless parameters $\mathrm{Re}=0.09,\ \mathrm{Pe}=136.26,\ \mathrm{Da}=5.16;$ these values are perfectly in line with our analysis, and the validity of our ERC equations.