Effects of diffusion boundary layer on reaction kinetics of immunoassay
in a biosensor
Chih-Kai Yang, Jeng-Shian Chang,a兲 Sheng D. Chao,b兲and Kuang-Chong Wu
Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan
共Received 28 June 2007; accepted 27 February 2008; published online 28 April 2008兲
Specific binding reaction is a natural characteristic that is applied to design biosensors. This work simulates the binding reaction kinetics of two commonly used proteins, C-reactive protein and immunoglobulin G, in a reaction chamber 共microchannel兲 of a biosensor. For a diffusion-limited protein, the diffusion boundary layer on the reaction surface of the biosensor would hinder the binding reaction from association and dissociation. Several crucial factors, which influence the binding reaction curves in the simulation, are discussed, including the concentration of analyte, the inlet flow velocity, the channel height, and the length of the reaction surface. A higher channel causes the diffusive transport of the analyte to take longer time to reach the reaction surface, which in turn decreases the reaction rate of the protein pairs. The length of the reaction surface plays an important role in the formation of the boundary layer. For longer reaction surface, it takes more time to allow diffusion to overcome the larger zone of the diffusion boundary layer, resulting in a slower binding rate and a longer time to reach saturation. The presented data of simulation are useful in designing the biosensors. © 2008 American Institute of Physics.关DOI:10.1063/1.2909980兴
I. INTRODUCTION
In recent years, the study of biosensors has become an active research field. It brings many interesting applications of nanotechnology to the area of biochemical and biological systems. The development of immunoassay places special emphasis on the high sensitivity and real-time detecting abil-ity for various biosensors. The three most commonly used devices in detecting biomolecules are the microcantilever-based biosensor, the surface plasmon resonance 共SPR兲 sen-sor, and the quartz crystal microbalance 共QCM兲 sensor. Al-though the basic principles of detection for these sensors are different, they all involve the same kinetics of specific rec-ognition of analytes and immobilized ligands, in which the concentration of the binding complex of analytes and ligands on the reaction surface plays a key role.
It is important to monitor the patient who suffers from a high risk epidemic such as cardiovascular disease. Human C-reactive protein共CRP兲, first discovered in 1930,1is a typi-cal acute-phase protein as a clinitypi-cal indicator of many in-flammatory conditions. CRP is usually present in human se-rum with a concentration ⬍1g/ml. However, CRP concentration level can increase up to tens or even hundreds times when the inflammation occurs. Human CRP is a ho-mogeneous protein 共molecular weight of ⬃115 kDa兲 which consists of five identical noncovalently bonded monomer subunits共molecular weight of ⬃23 kDa兲 linking in the form of a cyclic pentamer. In an immunoassay experiment, the common protein pair of antigen and antibody is human CRP and its monoclonal antihuman CRP. It is known that diffu-sion and binding ability are important factors in the reaction kinetics. These two factors diversely differ for different pro-tein pairs. Thus for the purpose of comparison, not only the
protein pair of human CRP and antihuman CRP but another common protein pair of mouse Immunoglobulin G共IgG兲 and antimouse IgG are simulated for the study of stoichiometry in this work.
The specific recognition of analytes and immobilized ligands occurs at the reaction surface of a biosensor, which is a solid-liquid interface. The reaction kinetics can be de-scribed as a two-step process.2
共1兲 Mass-transport process: The analyte is transported by diffusion from the bulk solution toward the reaction sur-face.
关A兴bulk 关A兴surface.
共2兲 Chemical reaction process: The binding of the protein pair takes place.
关A兴surface+关B兴 kd ka
关AB兴,
where关A兴bulkis the concentration of the analyte in the bulk, 关A兴surface is the analyte concentration at the reaction surface, 关B兴 is the ligand concentration, 关AB兴 is the analyte-ligand complex concentration, ka is the association rate constant,
and kdis the dissociation rate constant.
When the analyte takes a longer time to transport 共by convection and diffusion兲 to the reaction surface than chemi-cal reaction, the whole reaction is restrained by mass-transport process. This usually causes the formation of a dif-fusion boundary layer.3For the biosensors operating in fluid environment, the diffusion velocity of many biomolecules is relatively slow compared to the reaction velocity, as mea-sured by the Damkoehler number.4Thus on the reaction sur-face, where the flow is static under the assumption of nonslip flow condition, the diffusion boundary is formed. The size of the zone of the diffusion boundary layer depends on the ratio a兲Electronic mail: jschang@spring.iam.ntu.edu.tw.
b兲Electronic mail: sdchao@spring.iam.ntu.edu.tw.
of reaction rate to the diffusion rate and the flow velocity adjacent to the reaction surface.5 Besides, the size of the diffusion boundary layer will grow along the reaction surface in the flow direction. As the whole reaction 共association or dissociation兲 is complete, the diffusion boundary layer will disappear. Thus the size of the diffusion boundary layer is a good measure of reaction velocity.
In this paper we simulate two dimensional共2D兲 full time scale processes of association and dissociation in a biosensor immunoassay with the samples of CRP pairs and IgG pairs, using the finite element analysis software, COMSOL MULTIPHYSICS™.6Calculations of the surface concentration of the analyte-ligand complex versus time during the entire biochemical reaction are performed to study the effects of several related parameters on the detection of biomolecules with the biosensor. The plan of the paper is as follows. In Sec. II, the equations governing the flow field, the concen-tration field, and the biochemical reaction are described. De-tailed geometry, flow properties, diffusion constant, and binding constants, initial and boundary conditions that are required for simulation are described in Sec. III. Numerical calculations are then performed to solve the velocity field of the fluid flow, the concentration field of the analyte, in the bulk, the surface concentration of the analyte and the surface concentration of the analyte-ligand complex. Finally the re-sults are presented in Sec. IV. As for the biosensor with a narrower reaction surface, such as a microcantilever beam sensor, 2D simulations could show large deviation from three dimensional共3D兲 simulations. In Sec. IV F, we will demon-strate the difference of the binding reaction curves based on the 3D simulations and the 2D simulations for the reaction surface with various aspect ratios 共length to width兲.
II. THEORY
Consider a 2D model of a biosensor in a microchannel, as shown in Fig.1. The dimension of the biosensor is 40m in length and 3m in thickness and that of the microchannel is 500m in length and 150m in height. It is assumed that the reaction surface has the same length as that of the bio-sensor. The x-axis is set to be along the flow direction and
y-axis is the direction of the channel height. The left corner
of the bottom side is chosen as the origin. The buffer solution mixed with the analyte flows from the left to the right along the x-direction. On the reaction surface 共230艋x艋270,y = 3兲 of the biosensor are the immobilized ligands. The center of the biosensor is located in the middle of the bottom side, 共250, 1.5兲, where the unit is micrometers. Notice that the
The equations governing the flow field, the concentration field, and the biochemical reaction are given below.
A. The flow field
In this work it is assumed that the fluid is incompressible so that
u
x+
v
y= 0, 共1兲
where u and v are, respectively, the x and the y velocity
components. The equations of motion are as follows:
u t +
冉
u u x+v u y冊
−ⵜ 2u +p x= 0, 共2兲 v t +冉
u v x+v v y冊
−ⵜ 2v +p y= 0, 共3兲whereis the dynamic viscosity of fluid,is the density of fluid, and p is the pressure, respectively. In this work, it is assumed that the density and viscosity of the modeled incompressible fluid are constant independent of the tem-perature and the concentration.
B. The concentration field
Transport of analytes to and from the surface is assumed to be described by the Fick second law with convective terms: 关A兴 t + u 关A兴 x +v 关A兴 y = D
冉
2关A兴 x2 + 2关A兴 y2冊
, 共4兲where关A兴 共=关A兴bulk兲 is the bulk concentration of analyte and D is the diffusion coefficient of analyte.
C. The reaction surface
The reaction between immobilized ligand and analyte is assumed to follow the first order Langmuir adsorption model.3,7 During the reaction, the concentration of the analyte-ligand complex关AB兴 increases as a function of time according to the reaction rate equation
关AB兴
t = ka关A兴surface兵关B0兴 − 关AB兴其 − kd关AB兴, 共5兲
where关A兴surfaceis the concentration of the analyte at the re-action surface by mass transport, 关B0兴 is the initial surface concentration of the ligand,关AB兴 is the surface concentration of the analyte-ligand complex, kais the association rate
con-stant, and kdis the dissociation rate constant.
III. SIMULATION DETAILS
An unstructured mesh consisting of triangular elements is generated for the calculation. Figure 2 shows this mesh where the region near the biosensor is refined for a better resolution.
FIG. 1.共Color online兲 Sketch of the 2D model. The size of the biosensor is 40m in length and 3m in thickness. The channel size is 500m in length and 150m in height.
A. The flow field configuration
In an immunoassay experiment, phosphate buffer saline 共PBS兲 is usually used to be a neutral buffer solution 共pH = 7.2兲, which is mixed with analytes as a carrier fluid. The physical properties of such a fluid are similar to water. The dynamic viscosity is 10−3Pa s. Since the flow in the mi-crochannel is in a low Reynolds number condition, it is as-sumed to be a laminar flow. The average velocity of the parabolic profile is set as u = 10−4 m/s at the inlet of the microchannel. Boundary conditions are p = 0 at the outlet, and nonslip elsewhere. The fluid is assumed to be at rest initially. The flow field can be determined by solving Eqs. 共1兲–共3兲 in conjunction with the boundary and initial condi-tions.
B. The concentration field configuration and kinetics of the specific binding
The diffusion coefficients of human CRP are taken as 2.175⫻10−11m2/s and that of IgG as 5⫻10−11m2/s.8,9
Five inlet concentrations of the analyte关A兴, 0.64, 1.92, 6.4, 19.2, and 64 nM, are considered. The initial surface concentration 关B0兴 of the immobilized ligand is assumed as 1.4 ⫻10−8mole/m2.10
At the reaction surface, the diffusive flux is balanced against the reaction rate
− D
冉
关A兴y
冊
surface= ka关A兴surface兵关B0兴 − 关AB兴其 − kd关AB兴.
共6兲 The association rate constant kaand the dissociation rate
constant kdof the protein pairs can be found in Ref.11. The
values of ka and kd for CRP-anti-CRP binding interactions
are 1.0⫻107M–1s−1 and 2.6⫻10−2s–1, respectively. The
values of kaand kdfor IgG-anti-IgG binding interactions are
2.5⫻105M–1s−1 and 3⫻10−4s–1, respectively.
The initial conditions for both the concentration of the analyte in the bulk,关A兴, and the concentration of the analyte-ligand complex on the reaction surface,关AB兴, are all zero.
Once the flow field is determined, Eqs.共4兲 and 共5兲can be solved with the initial and boundary conditions described above to yield关A兴bulk,关A兴surface, and关AB兴 versus time during the entire reaction. In particular, the area average of 关AB兴 over the reaction surface is performed and presented in the simulation results below.
The dissociation phase of the binding reaction of the two protein pairs is simulated by terminating the supply of the analyte at a time after the binding reaction is saturated.
It is noted in our examples that the association reaction rates 共ka关B0兴兲 for CRP and IgG in the PBS solution are,
respectively, 140 and 3.5m/s. The diffusion rates 共D/h兲 for CRP molecule and IgG in the PBS solution are, respec-tively, 0.14 and 0.33m/s when h is picked as 150m. Since on the reaction surface the association rate is faster than the diffusion rate for both CRP and IgG, the diffusion boundary layer will form during the binding reaction. IV. RESULTS
We have performed the 2D finite element simulations to evaluate the binding curves during the association and disso-ciation and presented the results in Secs. IV A–IV E. The results are examined to be mesh independent by means of a convergence test. Effects of some crucial parameters are also discussed in this section. It is expected that for the reaction surface having small aspect ratio共length to width兲, namely, wider reaction surface, the 2D simulation should be accurate enough. However, for the system with relatively narrower reaction surface, 3D calculations are expected to be required. This is discussed in Sec. IV F
A. Binding kinetics of CRP and IgG
The simulated binding reaction curves for CRP and IgG are shown in Figs.3and4, respectively, corresponding to the different concentrations of the analyte, namely, 0.64, 1.92, 6.4, 19.2, and 64 nM. Although the dissociation equilibrium constants KDof the two proteins CRP and IgG are
compa-rable as medium affinity 共CRP:2.6⫻10−9M, IgG: 1.2
⫻10−9M兲, the characteristic behaviors in the association phase or the dissociation phase are different. The response of CRP is apparently faster than that of IgG. The main reason is that both values of kaand kdof CRP are greater than those of
IgG.
The analyte concentration 关A兴 versus the complex con-centration关AB兴satcan be obtained through Eq.共5兲by requir-ing关AB兴/t = 0共i.e., at the equilibrium of the association兲:
关A兴 ⬇ 关A兴surface=
kd/ka
关B0兴/关AB兴sat− 1
= KD
关B0兴/关AB兴sat− 1 , 共7兲 where KD⬅kd/ka is the dissociation equilibrium constant.
Equation共7兲is plotted in Fig.5, in which the squares denote the simulated data 共关A兴,关AB兴sat兲 picked from Figs. 3 and4. Consistency of the simulated data with Eq. 共7兲 serves to demonstrate the correctness of our simulation.
Equation 共7兲 is not convenient to use from the experi-mental point of view since the measured response R of a biosensor, such as QCM or SPR, is not the concentration of the analyte-ligand complex itself but rather a quantity that is assumed to be proportional to the concentration of the analyte-ligand complex 共i.e., R=␣关AB兴, where ␣ is a con-FIG. 2. 2D unstructured mesh with triangular elements.
stant兲. Thus a simple revision of Eq.共7兲 is desirable. In the experiment with a very dense concentration of the analyte 关Aˆ兴 such that KD/关Aˆ兴Ⰶ1, the initial concentration 关B0兴 of the ligand can be obtained as the maximum equilibrium response
Rˆeqdivided by the constant␣: 关B0兴 = 关ABˆ 兴sat
冉
1 +KD
关Aˆ兴surface冊⬇ 关ABˆ 兴sat
冉
1 +KD
关Aˆ兴
冊
⬇ 关ABˆ 兴sat= Rˆeq/␣, 共8兲
where关ABˆ 兴satis the maximum concentration of the analyte-ligand complex. When 关B0兴 is known, KDcan be estimated
as
KD⬇ 关A兴*共Rˆeq/Req*− 1兲, 共9兲
where Req* is the response at equilibrium in an experiment, in which the concentration of the analyte is 关A兴*. Finally, the relationship关A兴 versus Reqis given as
关A兴 ⬇ KD
Rˆeq/Req− 1
. 共10兲
In other words, by conducting only two experiments, one with very dense concentration of the analyte and the other one with a medium concentration of the analyte, one can determine Rˆeq and KD. Thus the concentration 关A兴 of the
analyte in the follow-up experiments can be predicted ac-cording to Eq.共10兲 when the response of the association at equilibrium is measured.
B. The diffusion boundary layer
Under the assumption of nonslip flow boundary condi-tion, the consumption of关A兴surfacein the association phase is faster than the supply from the bulk on the reaction surface, if the diffusion rate of analyte is slower than the reaction rate. Thus a small diffusion boundary layer of the analyte is formed near the reaction surface. Within the boundary layer, there is a shortage of the analyte. The diffusion boundary FIG. 3.共Color online兲 The average surface concentra-tion of CRP-anti-CRP complex along the surface as a function of time for different CRP bulk concentrations.
FIG. 4.共Color online兲 The average surface concentra-tion of IgG-anti-IgG complex along the surface as a function of time for different IgG bulk concentrations.
FIG. 5. 共Color online兲 The saturated concentration of analyte-ligand complex as a function of the concentra-tion of analyte in the bulk.
FIG. 6. 共Color online兲 The development of the diffusion boundary layer of the CRP binding reaction. The biosensor is located at 共250, 1.5兲 and the concentration of CRP is 6.4 nM. The left three illustrations are in association phase at times of 500, 1000, and 1500 s, and the right three illustrations are in dissociation phase at times of 3000, 4000, and 5000 s. Notice that the density scales are different to increase the plot visibility.
FIG. 7. 共Color online兲 The development of the diffusion boundary layer of the IgG binding reaction. The biosensor is located at 共250, 1.5兲 and the inlet concentration of IgG is 6.4 nM. The left three illustrations are in association phase at times of 1000, 2000, and 3000 s, and the right three illustrations are in dissociation phase at times 6000, 7000, and 8000 s. Notice that the density scales are different to increase the plot visibility.
FIG. 8.共Color online兲 Influence of raising the inlet flow velocity on the curve of the surface concentration of CRP complex vs time. Notice that the biosensor is lo-cated at共250, 1.5兲 and the inlet concentration of CRP is 6.4 nM.
layer also occurs in the dissociation phase when the supply of analyte in the flow is terminated. In the latter case the situation is reversed; namely, within the boundary layer, there is a denser concentration of the analyte than that in the bulk. Figures6and7show the time evolutions of the diffu-sion boundary layers of the CRP and IgG protein pairs, re-spectively, from the association phase共the left panels兲 to the dissociation phase共the right panels兲. Notice that the density scales are different to increase the plot visibility. With a con-tinuous supply of the analyte, the boundary layer is progres-sively eliminated during the association phase. Similarly, the boundary layer also gradually vanishes in the dissociation phase due to the termination of the supply of the analyte in the bulk. Obviously, any means to limit the growth of this boundary layer will raise the rates of the association and dissociation, as will be discussed in the next few subsections. C. Effect of inlet flow velocity
Raising flow velocity is an effective means to reduce the thickness of the diffusion boundary layer. Figure8 presents
the effect of raising the inlet flow velocity on the CRP-anti-CRP binding kinetics. As expected, the faster the inlet flow velocity is, the faster association and dissociation rates are.
D. Effect of height of the reaction channel
In this subsection the height of the reaction channel is allowed to vary. It is also noted that in each case of various channel heights, the flow rate is different but the average velocity共=100m/s兲 through the channel cross sections re-mains the same. In addition, the center of the biosensor is set at the position共250, 1.5兲 with the inlet analyte concentration of 6.4 nM. Figures 9 and 10 show the initial slopes of the binding curves for average surface concentration of the CRP complex and IgG complex along the surface, respectively. The height of the channel varies from 50 to 250m. Larger slopes are found for lower channel height for both CRP and IgG protein pairs. This can be understood as faster diffusion in transporting the analytes for the channel with a lower height.
FIG. 9.共Color online兲 The initial slope of CRP binding reaction as a function of channel’s height. The inlet concentration of CRP is 6.4 nM. The channel’s inlet flow velocity is 100m/s.
FIG. 10.共Color online兲 The initial slope of IgG binding reaction as a function of channel’s height. The inlet concentration of IgG is 6.4 nM. The channel’s inlet flow velocity is 100m/s.
E. Effect of the length of the reaction surface
To study the effect of the length of the reaction surface on the reaction rates, we vary the length of the reaction sur-face from 20 to 200m. The inlet concentration of the ana-lyte is constant 共=6.4 nM兲, and the supply of the analyte is sustained for 3000 s. The binding reaction curves for CRP are shown in Fig. 11, and the formation of the diffusion boundary layer on the reaction surface at time= 500 s is shown in Fig.12. As expected, the required time to reach a steady state gets longer when the span of the boundary layer is longer. The results for IgG are similar to those for CRP and are not shown here. Although shortening the length of the reaction surface could restrain the expansion of the boundary layer, the detecting techniques of a biosensor may require a sufficiently long reaction surface to satisfy its limit of detection. Consequently, some kind of compensation should be considered in the biosensor design.
F. Effect of width: 3D simulation
In some biosensors, the reaction surface may occupy the whole or half of the bottom surface of the microchannel such
as SPR and QCM. As for the microcantilever beam sensor, the size of reaction surface, especially the width, is rather small compared to that of the microchannel. Thus it is ap-propriate to compare the difference between the binding curve predicted by 2D simulation and that by 3D simulation. Figure 13 shows the sketch of a 500m共length兲 ⫻500m共width兲⫻150m共height兲 microchannel with a biosensor put on the middle of the bottom surface. There are four sizes of biosensor chosen to perform 3D simulation, namely, 40⫻500⫻3 m3, 40⫻200⫻3m3, 40⫻50 ⫻3m3, and 40⫻20⫻3m3. Results of 3D simulation to-gether with that of 2D simulation共the sensor with 40m in length and 3 m in height兲 for the average concentration of CRP complex along the length direction of the biosensor are shown in Fig. 14, where the average inlet flow velocity of 100m/s and the inlet bulk concentration of CRP of 6.4M are used. It is seen that the binding curve of the biosensor of width 500m, which is the same as the width of the micro-channel, predicted by 3D simulation is virtually the same as that predicted by 2D simulation. Besides, it is the trend that 3D simulation yields the faster binding reaction than 2D FIG. 11.共Color online兲 The average surface concentra-tion of CRP complex along the surface as a funcconcentra-tion of time for various lengths of the reaction surface. The concentration of CRP is 6.4 nM. The channel’s inlet flow velocity is 100m/s. The size of the channel is shown in Fig.1.
FIG. 12. 共Color online兲 The expansion of the diffusion boundary layer for varying lengths of the reaction surface. It is noticed that the time=500 s for 共a兲 20m,共b兲 30m,共c兲 40m,共d兲 80m,共e兲 120m, and共f兲 200m, respectively.
simulation does as the width of the biosensor gets smaller. Hence, it is suggested that 3D simulation is necessary for high aspect ratio共ratio of length to width兲 biosensor to yield more accurate prediction and that the biosensor with high aspect ratio will have better performance.
V. CONCLUSION
This paper presents 2D numerical simulations on the im-munoassay in a biosensor by using the finite element analysis software, COMSOL MULTIPHYSICS™. Two commonly used proteins, CRP and IgG, are used as the analytes for the analy-sis of the binding kinetics.
The existence of the diffusion boundary layer restrains the association and dissociation reaction rates for the slow diffusion protein pairs and nonslip flow boundary condition, which affects the performance of the biosensors. Several cru-cial factors have been discussed, namely, the concentration of the analyte, the inlet flow velocity, the height of the
channel, and the length of the reaction surface. It is expected that the denser concentration of the analyte yields the faster association and dissociation rates. Raising flow velocity can effectively reduce the thickness of the diffusion boundary layer. A lower channel causes the diffusion to take less time to transport the analyte to the reaction surface and hence increases the association and dissociation rates of the protein pairs. Moreover, the shorter the reaction surface along the flow direction is, the smaller the diffusion boundary layer is. Thus shortening the length of the reaction surface can also raise both the association and dissociation rates. In summary, the faster inlet flow velocity, the lower channel, and the shorter reacting surface can reduce the growth of the diffu-sion boundary layer and raise the performance of the biosen-sor. Finally, if the width of the biosensor is relatively small compared to that of the microchannel, 3D simulations are necessary to predict more accurate results. These results pre-sented in this work should be very useful for designing bio-sensors.
FIG. 13. 共Color online兲 Sketch of the 3D model. The channel size is 500m in length, 500m in width, and 150m in height. The biosensor, 40m in length, 3m in height, and variable width, is put on the center of the bottom surface.
FIG. 14.共Color online兲 The average surface concentra-tion of CRP-anti-CRP complex along the surface as a function of time for the reaction surface with various widths.
Council in Taiwan through NSC 96-2120-M-002-014. We thank the NCHC for providing computing resources.
1W. S. Tillet and T. Francis,J. Exp. Med.52, 561共1930兲. 2N. Camillone,Langmuir20, 1199共2004兲.
3D. B. Hibbert and J. J. Gooding,Langmuir18, 1770共2002兲.
4W. M. Deen, Analysis of Transport Phenomena共Oxford University Press,
7I. Langmuir,J. Am. Chem. Soc.40, 1361共1918兲.
8Y. Hokama, M. K. Coleman, and R. F. Riley, J. Immunol. 95, 156共1965兲. 9H. A. Leddy and F. Guilak,Ann. Biomed. Eng.31, 753共2003兲. 10E. Behravesh, V. I. Sikavitsas, and A. G. Mikos,Biomaterials24, 4365
共2003兲.
11C. Chou, H. Y. Hsu, H. T. Wu, K. Y. Tseng, A. Chiou, C. J. Yu, Z. Y. Lee, and T. S. Chan, J. Biomed. Opt. 12, 0240251共2007兲.
