Transport Effects on Multiple-Component Reactions in Optical Biosensors
Ryan M. Evans David A. Edwards










R. M. Evans Applied and Computational Mathematics Division

	       Information and Technology Laboratory
 National Institute of Standards and Technology,  Gaithersburg, MD 20899, USA
             
 ryan.evans@nist.gov      
           D. A. Edwards Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA

	    dedwards@udel.edu







Transport Effects on Multiple-Component Reactions in Optical Biosensors
This work was done with the support of the National Science Foundation under award number NSF-DMS 1312529. The first author was also partially supported by the National Research Council through an NRC postdoctoral fellowship.
    Ryan M. Evans        David A. Edwards 

    Received: date / Accepted: date
=========================================================================================================================================================================================================================================================================================================




Optical biosensors are often used to measure kinetic rate constants associated with chemical reactions.  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 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.




















§ 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_min,x_max] to produce measurements of the form

    B(t)=1/x_max-x_min∫_x_min^x_maxB(x,t) dx.

See Figure <ref> 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:  <cit.>.  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 <cit.>.  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 δ, and polymerase η.  Moreover, in order for a successful DNA translesion synthesis event to occur  polymerase η must bind with the PCNA molecule.  A central question surrounding DNA translesion synthesis is whether the polymerase η and PCNA complex forms through direct binding, or through a  catalysis-type ligand switching process <cit.>. 




The former scenario is depicted in Figure <ref>, where we have shown polymerase η directly binding with a PCNA molecule, 
i.e. the reaction:


    P_1: E+L_2[_2k_d]_2
    k_a  E L_2.

Here, we have denoted the PCNA molecule and polymerase η as E and L_2 respectively.  Additionally, _2k_ denotes the rate at which L_2 binds with an empty receptor E, and _2k_ denotes the rate at which L_2 dissociates from a receptor E.  We will refer to this as pathway one, or simply P_1 as in (<ref>).



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

    P_2: E+L_1 [_1k_d]_1k_a E L_1, 
    EL_1+L_2 [^1_2k_d]^1_2k_a E L_1 L_2 [^2_1k_a]^2_1k_d E L_2+L_1 , EL_2_2k_dE+L_2.

In (<ref>) and Figure <ref> we have denoted polymerase δ 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_1L_2; then L_1 dissociates from EL_1L_2, leaving EL_2; finally, L_2 dissociates from EL_2.  Furthermore, in (<ref>) and Figure 
<ref> the rate constants _1k_ and _1k_ denote the rates at which L_1 binds and unbinds with a receptor E,  ^j_ik_ denotes the rate at which ligand L_i binds with the product EL_j, and ^j_ik_ denotes the rate at which L_i dissociates from the product EL_1L_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 P_2 as in (<ref>). 






Though Zhuang et al. provided indirect evidence of the ligand switch in <cit.>,  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 (<ref>).  By measuring the rate constants in (<ref>), one could determine whether EL_1L_2  forms through direct binding, or the catalysis-type ligand switching process.  We note that the latter manifests itself mathematically with _2k_a=0, while the former with _1^2k_a= _1^2k_d= _2^1k_a= _2^1k_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_1L_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

    𝒮(t)=s_1B_1(t)+(s_1+s_2)B_12(t)+s_2B_2(t)

are produced.  In (<ref>)

    B_i(t)=1/x_max-x_min∫_x_min^x_maxB_i(x,t) dx

denotes the average reacting species concentration, for i=1, 12, 2, and s_i denotes the molecular weight of L_i. The lumped signal (<ref>) raises uniqueness concerns, since more than one set of rate constants may possibly correspond to the same signal (<ref>).  Fortunately, through varying the uniform in-flow concentrations of the ligands,  C_1(0,y,t)=C_1,u 
and C_2(0,y,t)=C_2,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 <cit.>.

The presence of multiple interacting species and the lumped signal (<ref>) complicate parameter estimation even for systems accurately descibed by the well-stirred kinetics approximation.  However in <cit.>, Edwards has shown that transport dynamics affect ligand binding in a thin boundary layer near the sensor surface.  Hence, we begin in Section <ref> by summarizing the relevant boundary layer equations, which take the form of a set of nonlinear Integrodifferential Equations (IDEs).   In Section <ref>, 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 <ref>.  Convergence properties are examined in Subsection <ref>, and in Section <ref> the accuracy of our ODE approximation is verified by comparing results of our ODE model with results from our numerical method described in Section <ref>.  Conclusions and plans for future work are discussed in Section <ref>.






§ 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,u and C_2(x,y,t)=C_2,u.  Injection continues until the signal (<ref>) 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.



 §.§ Injection Phase


To present our governing equations we introduce the dimensionless variables:

    x=x/L, y=y/H, t=_1k_aC_1,ut, B_i(x,t)=B_i(x,t)/R_T, C_i(x,y,t)=C_i(x,y,t)/C_i,u,
       _i^jK_=C_i,u·_i^jk_/C_1,u·_1k_, _i^jK_=k_/C_1,u·_1k_, F_r=C_rD_r, C_r=C_1,u/C_2,u, D_r=D_1/D_2.
 
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^jK_a and _i^jK_d are the dimensionless analogs of _i^jk_ and _i^jk_.  In the latter expressions the index i=1, 2, whereas j=1, 2, or can be blank.  Furthermore, F_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

    S(t)=𝒮(t)/R_T· s_1=B_1(t)+(1+s_2/s_1)B_12(t)+s_2/s_1B_2(t).

 Moreover, we may use (<ref>) 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 (<ref>) gives the kinetics equations:


    B_1/ t   = (1-B_Σ)C_1(x,0,t)- _1K_B_1 - ^1_2K_B_1C_2(x,0,t) + ^1_2K_B_12, 
        B_12/ t   =_2^1K_aB_1C_2(x,0,t)-_2^1K_dB_12+_1^2K_aB_2C_1(x,0,t)-_1^2K_dB_12,  
        B_2/ t   = _2K_(1-B_Σ)C_2(x,0,t)- _2K_B_2+ _1^2K_B_12  - ^2_1K_B_2C_1(x,0,t),  
       𝐁(x,0)=0,


which hold on the reacting surface when y=0 and x∈[0,1].  In (<ref>), 𝐁=(B_1, B_12, B_2)^T is a vector in ℝ^3 whose components contain the three bound state concentrations.  In addition, the terms in equations (<ref>)–(<ref>) have been ordered in accordance with Figures <ref> and <ref>.

Edwards has shown <cit.> 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


    D_r^2 C_1/η^2=η C_1/ x, 
       ^2 C_2/η^2=η C_2/ x.

In (<ref>)–(<ref>): η=^1/3y is the boundary layer variable,  =VH^2/(LD_2)≫ 1 is the Péclet number, and V is the characteristic velocity associated with our flow.  

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:

    C_1/η(x,0,t)=/F_r(   B_1/ t+ B_12/ t),
        C_2/η(x,0,t)=(   B_12/ t+ B_2/ t).

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

    C_i(0,η,t)=1, 
    lim_η→∞C_i(x,η,t)=1,


the solution to (<ref>) is given by


    C_1(x,0,t)=1-D_r^1/3/F_rΓ(2/3)3^1/3∫_0^x( B_1/ t+ B_12/ t)(x-ν, t)ν/ν^2/3, 
       C_2(x,0,t)=1-/Γ(2/3)3^1/3∫_0^x( B_12/ t+ B_2/ t)(x-ν, t)ν/ν^2/3.

See <cit.> for details of a similar calculation.  During the injection phase, the bound state concentration is then governed by (<ref>) using (<ref>).  

In (<ref>)–(<ref>) and (<ref>)

    =_1k_aR_T(HL)^1/3/(VD^2)^1/3
 
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 ≪ 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 (<ref>) shows that the unbound concentration at the surface is only a perturbation away from uniform inlet concentration.  When → 0 in (<ref>) using (<ref>), one recovers the well-stirred approximation in which transport into the surface  completely decouples from reaction.


On the other hand, when =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 (<ref>).  When x≪ 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 <ref>, where we have shown the results of our numerical simulations described in Section <ref>. 

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_1L_2 to form.  Thus, in this case EL_1L_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 <ref>–this reflects the weakly singular nature of the functions which we are attempting to approximate.  When x≪ 1, one may show 𝐁 has the perturbation expansion

    𝐁(x,t)=^0𝐁(t)+ x^1/3·^1𝐁(t)+ O(^2 x^2/3)

(this is simply (<ref>) for x≪ 1).  It therefore follows that

    𝐁/ x
    (x,t)= ^1𝐁(t)/3x^2/3+O( ^2/x^1/3).

Hence, although the function 𝐁 is well-defined and continuous near x=0, it has a vertical tangent at x=0.  The weakly-singular nature of 𝐁 is magnified since =2.  To resolve this region, one may think to 
adaptively change Δ x with the magnitude of 𝐁/ x.  However, because the sensogram reading S(t) is computed over the region [x_min,x_max], we are not concerned with resolving this region and  a uniform step size is sufficient.  Moreover, our convergence results in Subsection <ref> demonstrate that a lack of resolution at x=0 does not affect our results in the region of interest [x_min,x_max].




 §.§ 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
 <cit.>.  This implies that because the bound ligand concentration evolves on a much slower time scale than the unbound ligand concentration <cit.>, 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 (<ref>), with (<ref>) replaced by the steady solution to (<ref>) during the injection phase:
 

    𝐁(x,0)=A^-1𝐟,
       A=[ (1++)   1 -     1;     -   (+)     -;           - ( ++) ] ,
       =[ 1; 0;   ].
 

Equations similar to (<ref>) hold:


    C_i(0,η,t)=0, 
    lim_η→∞C_i(x,η,t)=0.
Equation (<ref>) is the inflow condition, and (<ref>) expresses the requirement that the concentration in the boundary layer must match the concentration C_i(x,y,t)=0 in the outer region.  Moreover, as in the injection phase one can use (<ref>)–(<ref>) together with (<ref>) to show:
 

    C_1(x,0,t)=-D_r^1/3/F_rΓ(2/3)3^1/3∫_0^x( B_1/ t+ B_12/ t)(x-ν, t)ν/x^2/3, 
       C_2(x,0,t)=-/Γ(2/3)3^1/3∫_0^x( B_12/ t+ B_2/ t)(x-ν, t)ν/ν^2/3.

Thus, during the wash phase the bound state evolution is governed by the (<ref>)–(<ref>), (<ref>), and (<ref>).




§ EFFECTIVE RATE CONSTANT APPROXIMATION
 


During both phases of the experiment, the bound state concentration 𝐁(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 𝐁(t), rather than the spatially-dependent function 𝐁(x,t), since from 𝐁(t) we can construct the sensogram signal (<ref>) (the quantity of interest).  Thus, we seek to find an approximation to 𝐁(t), and begin by
 finding one during the injection phase.  We first average each side of (<ref>), with C_1(x,0,t) and C_2(x,0,t) given by (<ref>), in the sense of (<ref>).  Immediately, we are
 confronted with terms such as

    B_1C_2=B_1(1-/3^1/3Γ(2/3)∫_0^x( B_12/ t+ B_2/ t)ν/(x-ν)^2/3),

on the right hand side of (<ref>).   In the experimentally relevant case of small , we are motivated to expand B(x,t) in a perturbation series:

    B(x,t)=^0 B(x,t)+O().

In this limit, the leading order of (<ref>) is just C_i=1.  Using this result in (<ref>), we have that the governing equation for ^0 B is independent of x:

    d ^0 Bdt=-A ^0 B+𝐟,

where A is given by (<ref>) and 𝐟 by (<ref>).  Hence the leading-order approximation

    0𝐁(t)=A^-1(I-^-At)𝐟

is independent of space. Substituting (<ref>) into (<ref>), 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 (<ref>) at O(); hence we may write

    𝐁(x,t)=^0𝐁(t)+ x^1/3·^1𝐁(t)+O(^2).

As a result of (<ref>) we have the relation

    B_i(x,t)=0B_i(t)+O(^2),

which may be used to show  the right hand side of (<ref>) is equal to

    B_1- h·0B_1(0B_12/ t+0B_2/ t)+O(^2),
       h(x)=3^2/3x^1/3/Γ(2/3).

We then average (<ref>), and use the resulting relation in (<ref>) to show the right hand side of (<ref>) reduces to:


    B_1 C_2=B_1[ 1- h(B_12/ t+B_2/ t)]+O(^2).

In this manner, we can derive a set of nonlinear ODEs for 𝐁(t) of the form:


    𝐁/ t=M^-1(𝐁)(-A𝐁+)+O(^2), 
       𝐁(0)=0,

where

    M(𝐁)=I+ N(𝐁),
       N(𝐁)=[       D_r^1/3h/F_r (1-) D_r^1/3h/F_r (1-)  - h·                     -h·;                      h·       h·+(D_r^1/3h/F_r)         ( D_r^1/3h/F_r);        -( D_r^1/3h/F_r)   -(D_r^1/3h/F_r)+h(1-)                   h(1-) ].

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


    𝐁/ t=M^-1(𝐁)(-𝒟𝐁)+O(^2), 
       𝐁(0)=A^-1𝐟, 
       𝒟=[   - 0; 0 + 0; 0 -   ],

where M(𝐁) is as in (<ref>)  


Following <cit.>, we refer to the Ordinary Differential Equation (ODE) systems  (<ref>) and (<ref>) 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 (<ref>) or (<ref>), 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 (<ref>).  Since experimental data is still forthcoming, we do not employ a regression algorithm to fit the rate constants in (<ref>) and (<ref>) 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 <ref>.

First consider the solutions depicted on the left.  Here the injection phase (<ref>) has been run from t=0 to t=5 and the wash phase (<ref>) has been run from t=5 to t=10.  Furthermore, all rate constants were taken equal to one and the Damköhler number was =0.1.  During the injection phase it is seen that B_1 and B_2 reach equilibrium after approximately one second, while 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_1L_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 𝐁 during the injection phase is given by (<ref>), and one can readily verify that A^-1𝐟=(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  B_12 decays to zero faster than the other two species: EL_1L_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 <ref>.  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^1K_a, which was taken equal to _2^1K_a=10. During the injection phase it is seen that B_1 quickly reaches a local maximum, and then decreases to steady-state. Since _2^1K_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_1L_2 depicted on the right in Figure <ref>.  From these observations it is clear why the steady-state value of B_12 is larger than the previous case.  However, it may be counterintuitive to observe that  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_1L_2, the increase in EL_1L_2 also increases the concentration of empty receptor sites. The continuous injection of L_2 therefore drives the average concentration B_2 to a larger steady-state value.  During the wash phase, it is seen that B_1 reaches a global maximum after approximately t≈5.75 seconds.  The increase in B_1 during the wash phase is a direct consequence of L_2 molecules dissociating from EL_1L_2.  Since only pure buffer is flowing through the biosensor during the wash phase, it is seen in Figure <ref> that each of the average concentrations B_i decay to zero.



















§ NUMERICS
 

To verify the O(^2) accuracy of our ERC approximation derived in Section <ref>, we now develop a numerical approximation to the IDE system (<ref>), where C_1(x,0,t) and C_2(x,0,t) are given by (<ref>).  We focus on the injection phase, since the wash phase is similar. Our approach is based on the numerical method described in <cit.>.  Semi-implicit methods have been previously used with great
success to solve reaction-diffusion equations <cit.>, 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 <ref>;  however,  we first turn our attention to deriving our numerical method in Subsection <ref>.






 §.§ Semi-implicit finite difference algorithm



We discretize the spatial interval [0,1] by choosing N+1 equally spaced discretization nodes x_i=iΔ x, for i=0, …, N, and discretize time by setting t_n=nΔ t, for n=0, …. Having chosen our discretization nodes and time steps, we seek to discretize (<ref>), where C_1(x,0,t) and C_2(x,0,t) are given by (<ref>).  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 (<ref>), however the integrand of C_1(x,0,t) is singular when ν=0.  To handle the singularity we subtract and add

    ( B_1/ t+ B_12/ t)(x-ν,t)|_ν=0

from the integrand.  Doing so yields

    C_1(x,0,t)=   1-D_r^1/3/F_r 3^1/3Γ(2/3) {∫_0^x[( B_1/ t+ B_12/ t)(x-ν,t)
        - ( B_1/ t+ B_12/ t)(x,t)]dν/ν^2/3+3 x^1/3( B_1/ t+ B_12/ t)(x,t)},

where we have used the fact that (<ref>) is independent of ν.  Then choosing a discretization node x=x_i and a time step t=t_n, we apply the trapezoidal rule to (<ref>) to obtain

    C_1(x_i,0,t_n)=   1-D_r^1/3/F_r 3^1/3Γ(2/3) {0·/2+∑_j=1^i-1[( B_1/ t+ B_12/ t)(x_i-x_j,t_n)
        - ( B_1/ t+ B_12/ t)(x_i,t_n)]Δ x/x_j^2/3+[( B_1/ t+ B_12/ t)(0,t_n)
        -( B_1/ t+ B_12/ t)(x_i,t_n) ]/2 x_i^2/3+3 x^1/3_i( B_1/ t+ B_12/ t)(x_i,t_n)},

when x_i>0; simply evaluating (<ref>) 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 <cit.> we have

    lim_ν→ 0( B_k/ t(x-ν, t_n) - B_k/ t(x,t))1/ν^2/3=lim_ν→ 0ν^1/3( B_k/ t(x-ν, t_n) - B_k/ t(x,t))1/ν,

which implies

    lim_ν→ 0( B_k/ t(x-ν, t_n) - B_k/ t(x,t))1/ν^2/3=lim_ν→ 0ν^1/3^2 B_k/ x t(x,t)=0,

for k=1, 12, or 2.  The last equality follows since we expect B_i/ t to be sufficiently regular for fixed x>0.  The expansion (<ref>) shows that this is certainly true when ≪ 1, however when =O(1) or larger the nonlinearity in (<ref>) renders any analytic approximation to B_i beyond reach.  Our results in Subsection <ref> show that our method indeed converges when =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 

    B_j(x_i,t_n)≈ B^j_i,n,

and approximate the time derivatives through the formula

    B_j/ t(x_i,t_n)≈B^j_i,n-B^j_i,n-1/Δ t:=Δ B^j_i,n/Δ t.

Our approximation (<ref>) 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 Δ 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 (<ref>), the fully-discretized version of C_1(x,0,t) is given by substituting (<ref>) into (<ref>):

    C^1_i,n=   1-D_r^1/3/F_r 3^1/3Γ(2/3) {∑_j=1^i-1[(Δ B^1_i-j,n/Δ t+Δ B^12_i-j,n/Δ t)- (Δ B^1_i,n/Δ t+Δ B^12_i,n/Δ t)]Δ x/x_j^2/3
        +[(Δ B^1_0,n/Δ t+Δ B^12_0,n/Δ t) -(Δ B^1_i,n/Δ t+Δ B^12_i,n/Δ t)]/2 x_i^2/3+3 x^1/3_i(Δ B^1_i,n/Δ t+Δ B^12_i,n/Δ t)},

for i>0, and C^1_0,n=1.  The function C_2(x,0,t) has a similar discretization which we denote as C_2(x_i,0,t_n)≈ C^2_i,n.  Thus, our numerical method takes the form:

    Δ B^1_i,n+1/=(1-B^Σ_i,n)C^1_i,n+1- B^1_i,n- B^1_i,nC^2_i,n+1+ B^12_i,n,
       Δ B^12_i,n+1/= B^1_i,nC^2_i,n+1- B^12_i,n+ B^2_i,nC^1_i,n+1- B^12_i,n,
       Δ B^2_i,n+1/=(1-B^Σ_i,n)C_i,n+1^2- B^2_i,n+ B^12_i,n- B^2_i,nC^1_i,n+1.



We enforce the initial condition (<ref>) at our N+1 discretization nodes through the condition B^j_i,0=0 for j=1, 12, 2, and i=1,…, N.  Observe that our method (<ref>) is semi-implicit rather than fully-implicit.  This renders (<ref>) linear in Δ B^j_i,n+1, and as a result we can write


    Δ𝐁_i,n+1/=M^-1_i,n(𝐁_i,n)(A^-1_i,n+1𝐁_i,n+𝐟_i,n+1),

where 𝐁_i,n=(B^1_i,n, B^12_i,n, B^2_i,n)^T.  Hence, by using a method which is only semi-implicit in time we  avoid the expense and complication of solving a nonlinear system at each time step.  Having solved for Δ𝐁_i,n+1 using (<ref>), we march forward in time at a given node x_i through the formula

    𝐁_i,n+1=𝐁_i,n+1/2(3Δ𝐁_i,n+1-Δ𝐁_i,n),



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 (<ref>) we see  C_1(x,0,t) and C_2(x,0,t) depend on B_j(ν,t) only for ν≤ 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 (<ref>) the matrix M^-1_i,n(𝐁_i,n) depends only upon 𝐁_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 𝐟_i,n+1 depend upon 𝐁_l,n+1 for l<i.  Thus, at each time step n+1 we first determine 𝐁_0,n+1.  Next, we increment i and use the value of 𝐁_0,n+1 in (<ref>) to determine 𝐁_1,N+1.  We proceed by iteratively marching our way downstream from x_2 to x_N to determine 𝐁_2,n+1, …, 𝐁_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 (<ref>) was initialized with one step of Euler's method.  

Furthermore, with our finite difference approximation to 𝐁(x,t), we can determine the average quantity 

    B(t)=(B_1(t), B_12(t), B_2(t))^T

with the trapezoidal rule

    𝐁(t_n)≈1/x_max-x_min(/2𝐁_m,n+∑_i=m+1^M-1𝐁_i,n+/2𝐁_M,n).

In (<ref>), the indices i=m and i=M correspond to x_min=m and x_max=M.  Our nodes were chosen to align with x_min and x_max to avoid interpolation error.



 §.§ Convergence study
 






  §.§.§ Spatial Convergence


We now examine the spatial rate of convergence of our numerical method.  Since from 𝐁 we can compute the quantity of interest (<ref>), we derive estimates for the rate at which our numerical approximation converges to 𝐁.  Furthermore, because the system (<ref>), (<ref>) is nonlinear, our analysis will focus on the experimentally relevant case of ≪ 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 (<ref>), (<ref>).  Averaging (<ref>) in the sense of (<ref>) gives:


    dB_1/d t   = (1-B_Σ)C_1(x,0,t)- _1K_B_1  - ^1_2K_B_1C_2(x,0,t)+ ^1_2K_B_12 ,
    	  
       dB_12/d t   = ^1_2K_B_1C_2(x,0,t)- ^1_2K_B_12 + ^2_1K_B_2C_1(x,0,t) 
      - _1^2K_B_12,  
       dB_2/d t   = _2K_(1-B_Σ)C_2(x,0,t)- _2K_B_2+ _1^2K_B_12  - ^2_1K_B_2C_1(x,0,t).  
       𝐁(0)=0.


As in Subsection <ref>, we handle the singularity in (<ref>) by adding and subtracting (<ref>) from the integrand of (<ref>) to write C_1(x,0,t) as in (<ref>).  The unbound ligand concentration C_2(x,0,t) has a representation analogous to (<ref>).  In the following analysis we limit our attention to (<ref>), since the analysis for equations (<ref>)–(<ref>) is nearly identical.

We proceed by anaylzing each of the terms in (<ref>):


    -B_1,
       B_12,
       C_1(x,0,t),
       -B_ΣC_1(x,0,t),
       -B_1 C_2(x,0,t).


Upon inspecting (<ref>) and using linearity of the averaging operator, we see that three terms contribute to (<ref>):

    1
       -D_r^1/3/F_r3^1/3Γ(2/3)∫_0^x( B_1/ t+ B_12/ t)(x-ν,t)-( B_1/ t+ B_12/ t)(x,t)dν/ν^2/3,
       3x^1/3( B_1/ t+ B_12/ t).

In a similar manner, (<ref>) and (<ref>) each imply that we incur error from the terms:

    - B_Σ
       D_r^1/3 B_Σ/F_r3^1/3Γ(2/3)∫_0^x( B_1/ t+ B_12/ t)(x-ν,t)-( B_1/ t+ B_12/ t)(x,t)dν/ν^2/3,
       -3x^1/3B_Σ( B_1/ t+ B_12/ t),
       - B_1
        B_1/3^1/3Γ(2/3)∫_0^x( B_12/ t+ B_2/ t)(x-ν,t)-( B_12/ t+ B_2/ t)(x,t)dν/ν^2/3,
       -3 x^1/3B_1( B_12/ t+ B_2/ t).



Let us denote the trapezoidal rule of a function f(x) over the interval [a,b] by 𝒯(f(x),[a,b]).  Then since 𝒯(1,[x_min,x_max]) is exact, the term (<ref>) does not contribute to the spatial discritization error.

Next we decompose the expansion (<ref>) into its individual components to obtain

    B_j(x,t)=^0B_j(t)+ x^1/3·^1B_j(t)+O(x^2/3^2),
  
for j=1, 12, 2.  Substituting (<ref>) into (<ref>), (<ref>), (<ref>) (<ref>),  and using the fact that 𝒯(x^1/3,[x_min,x_max]) converges at a rate O(^2) shows that each of these terms converge at a rate of O(^2).  Similarly, one can substitute (<ref>) into (<ref>), (<ref>), and (<ref>), and use the fact that 𝒯(x^1/3,[x_min,x_max]) converges at a rate of O(^2), to show that each of these terms converge at a rate of of O(^2).  

It remains to determine the error associated with (<ref>), (<ref>), and (<ref>), so we turn our attention to (<ref>) and substitute (<ref>) into (<ref>) to obtain

    -D_r^1/3/F_r3^1/3Γ(2/3)(x_max-x_min)(d^1B_1/d t+d^1B_12/dt)(t)∫_x_min^x_max∫_0^x(x-ν)^1/3-x^1/2)dν/ν^2/3,

where we have used the definition of our averaging operator (<ref>).  In writing (<ref>), neglected higher-order terms which do not contribute to the leading-order spatial discretization error.  Since the coefficient of the integral in (<ref>) 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 (<ref>), we calculate the error  associated with applying the trapezoidal rule to the double integral

    ∫_x_min^x_max∫_0^x[(x-ν)^1/3-x^1/3]dν/ν^2/3 dx.

Treating the inner integral as a function of x we define

    f(x)=∫_0^x ( (x-ν)^1/3-x^1/3)ν^-2/3 ν,

whose closed form is given by

    f( x)= x^2/3/2(2^1/3√(π)Γ(1/3)/Γ(5/6)-6).

Towards applying the trapezoidal rule to (<ref>), we first note 𝒯(f,[0,x_i]) converges at a rate of O(^4/3).  This is seen by first rewriting (<ref>) as

    ∫_0^Δ x [(x-ν)^1/3-x^1/3]ν^-2/3 ν+∫_Δ x^ x_i [(x-ν)^1/3-x^1/3]ν^-2/3 ν.

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


Applying the trapezoidal rule to (<ref>) then gives 

    ^2∫_^∫_0^x((x-ν)^1/3-x^1/3) ν^-2/3 ν  x
          =^2/2𝒯(f(x),[0,x_m])+∑_i=m+1^M-1^2𝒯(f(x),[0,x_i]) 
           +^2/2𝒯(f(x),[0,x_M])+O(^2^2),

where we have let x_min=x_m=m, and x_max=x_M=M.  Since 𝒯(f,[0,x_i]) converges at a rate of O(^4/3), the right hand side of the above is 

    (^2/2f(x_m)+O(^2^7/3))+∑_i=m+1^M-1(^2 f(x_i)+O(^2^7/3))
         +(^2/2f(x_M)+O(^2^7/3))+O(^2^2).

To compute our results in Section <ref>, we took =0.2, =0.8, in accordance with the literature <cit.>.  Hence, in the above sum there are approximately 0.6N=0.6^-1 terms on the order of O(^2^7/3), and the above sum reduces to

    ^2 (/2f(x_m)+∑_i=m+1^M-1 f(x_i)
     +/2f(x_M))+O(^2^4/3)
         +O(^2^5/3).

The dominant error in (<ref>) is O(^2^4/3), thus the spatial discretization error associated with (<ref>) is O(^2^4/3).  When measuring convergence we used values of =.25, =.75 to facilitate progressive
 grid refinement; however, it is clear that these values of  and  do not change our argument.  A similar argument shows the spatial discretization error associated with the nonlinear terms (<ref>) and (<ref>) is O(^2^4/3).   

We have depicted our spatial convergence measurements for B_1 in Figure <ref> and tabulated them in Table <ref>.  To obtain these results, we first computed a reference solution, with ==1/512.  We then created a series of test solutions with mesh width =1/2^j, for j=2,…, 7, keeping =1/512 constant. Next, we computed 𝐁 by averaging our reference solution and test solutions at each time step with the trapezoidal rule as in (<ref>).  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(^2) when ≪ 1, O(^4/3) when =O(1), and O(^3/2) when ≫ 1.  The reduction in convergence when  increases from small to moderate may be attributed to the O(^2^4/3) contributions from (<ref>), (<ref>), and (<ref>).  There are
 two competing magnitudes of error in (<ref>): one of O( ^2) (from terms  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>)), and one of O(^2 ^4/3) (from the integral terms (<ref>), (<ref>), and (<ref>)).  When ^2^4/3<  ^2, or <^2/3, the former is larger.  Conversely, when ^2/3< the latter is larger. 

 
When ≫ 1, the bound state evolves on a longer time scale of the form <cit.>

    t_w=t/.

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  ≫ 1 as the transport-limited regime.  Substituting (<ref>) into (<ref>), (<ref>), one may find the leading-order approximation to the resulting system for ≫ 1 by neglecting the left hand side of (<ref>)–(<ref>).  Doing so one finds that even a leading-order approximation to (<ref>)–(<ref>) is nonlinear, rendering any error estimates in the transport-limited regime beyond reach. Nonetheless, our results in Figure <ref> and Table <ref> show that convergence is not an issue when ≫ 1.



  §.§.§ Temporal Convergence


Since our time stepping method (<ref>) is analogous to a second-order Adams-Bashforth formula, we expect our method to achieve second-order accuracy in time.  Figure <ref> shows that this is indeed the case when =.01, and the rate constants are =1, =1, =1, =1/2, =2, =2, and =1/2 (as in Subsection <ref> when measuring spatial convergence).  Temporal convergence was measured in an analogous manner to spatial convergence.

 However, we note that measuring temporal convergence when =O(1) is computationally prohibitive, since spatial convergence is O(^2^4/3) in this case, so in order for the spatial and temporal errors to balance one must have O(^2^4/3)=O(^2) or =^3/2.  
 
 Nonetheless, our results from Section <ref> 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 =O(1) or larger.






  §.§.§ Stability Remarks


We now make brief remarks concerning the stability of our method.  Recall from Subsection <ref> 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 (<ref>) 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 (<ref>).    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. 







§ 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 (<ref>) and (<ref>).    We tested the accuracy of our ERC equations when =0.1, and =0.45;  the results are below in Figure <ref> and Tables <ref>–<ref>.














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







§ CONCLUSIONS



Scientists are attempting to determine whether the polymerase η and PCNA complex (denoted EL_2 throughout) which results from DNA translesion synthesis forms through direct binding (<ref>), or through a catalysis-type ligand switching process (<ref>).  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 (<ref>). In turn, estimates for the rate constants in (<ref>) will reveal whether the polymerase η and PCNA complex forms via direct binding (<ref>), or the catalysis-type ligand switching process (<ref>).



 

Furthermore, the consideration of both direct binding (<ref>) and the ligand switching process (<ref>) has several mathematical and physical consequences.  First, due the form of (<ref>), the species are directly coupled through the kinetics equations.  This is true even in the well-stirred limit in which → 0, and (<ref>) reduces to C_1(x,0,t)=C_2(x,0,t)=1.  However, transport effects manifested in (<ref>) nonlinearly couple the reacting species.  So we see in Figure <ref> 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_1L_2 to form, thus the latter is affected by depletion of the former two species.  Additionally, the multiple-component reactions (<ref>) alter the form of the sensogram reading to the lumped signal (<ref>), thereby complicating parameter estimation.



In addition to establishing a firm foundation for studying the inverse problem of estimating the rate constants in (<ref>), 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 <cit.>.







§ 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 <ref>.  The flow rate is related to the velocity through the formula <cit.>

    V=6Q/WH.


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

Here ϵ=H/L is the aspect ratio, and Re=VH^2/(ν L) is the appropriate Reynolds number associated with our system.

The authors wish to emphasize that the bounds in Table <ref> are naïve extremal bounds calculated by using minimum and maximum values for the dimensional parameters in Table <ref>.  In particular, the values for the dimensionless rate constants in Table <ref> are not estimates of their true values; they are minimum and maximum values calculated using combinations of extremal values for the parameters in Table <ref>.  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 , would necessitate very small values for either or both of Δ t and Δ 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 μL/min–the slowest flow rate possible on the BIAcore T200 <cit.>.  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_T, and decreasing the ligand inflow concentrations C_1,u and C_2,u.  In the case of the fastest reaction _1k_=3×10^9 cm^3/(mol·s), we can take:
Q=390 μ L/min,V=.75 cm/s,R_T=7.76× 10^-13 mol/cm^2,C_1,u=C_2,u=2.96×10^-12 mol/cm^3 .
These choices yield the dimensionless parameters Re=0.09, Pe=136.26, Da=5.16; these values are perfectly in line with our analysis, and the validity of our ERC equations. 






plainnat



