NUMERICAL MODELING Basic definitions Let f(x,y) and g(x,y) denote the amount of the free PCR product and extended solid phase forward primer, respectively, at the position (x,y) (Figure I). Let u(x,y) denotes amount of annealed PCR product at position (x,y) and v(x,y) denotes the amount of annealed reverse primers at position (x,y). Free reverse primers Free PCR product, fi(x,y) gi(x,y) Annealed PCR ui(x,y) product, ui(x,y) vi(x,y) Extended forward primers, gi(x,y) Annealed reverse primer, vi(x,y) Solid phase forward primers Melting Annealing Elongation Figure I. Scheme of experiment and the basic disinganions. At the melting stage free and immobilised primers are totally separated from the extended forward primers (gi(x,y)) and free PCR products (fi(x,y)), respectively. At the annealing stage free reverse primers and free PCR products are annealed to their immobilised counterparts (denoted in numerical modelling as vi(x,y) and ui(x,y), respectively). At the elongation stage all the annealed primers are extended to form full length double stranded products. As this model accounts for all three stages of PCR, it requires some indexes for f, g, u and v. The first numerical index denotes the number of PCR cycle. The second index is a letter m, a, or e, denoting melting, annealing, or elongation stage of PCR cycle, respectively. As at the annealing stage the hybridisation and diffusion occure in parallel, this stage requires a separate iterative process and an additional numerical index for denotation of these iterations. E.g., fiaj(x,y) specifies the amount of free PCR product at i-th cycle at j-th iteration (or j-th second) of annealing stage at position (x,y); gie(x,y) specifies the amount of extended forward primers at i-th cycle of PCR at the end of elongation stage at position (x,y). Melting stage As no primers are hybridised to their targets at the melting stage, uim(x,y)=vim(x,y)0. The amount of extended forward primers remains unchanged during the melting stage: g0m(x,y)0 gim(x,y)=gi-1e(x,y). The free PCR product diffuses during the melting stage: f0m(x,y)1 fim(x,y)=DmMfi-1e(x,y). Here Dm is the operator transforming the initial distribution of PCR product to its distribution after one second of diffusion; Dm is a convolution of f with fundamental solution of the diffusion equation, Dmf(x,y)=exp(-(x2+y2)/rm2)*f(x,y) where rm is the diffusion distance per one second of melting stage, M is the duration of melting stage (in seconds). Annealing stage The initial conditions of annealing stage correspond to the end of preceding melting stage: fia0(x,y)=fim(x,y), gia0(x,y)=gim(x,y), via0(x,y)=uia0(x,y)0, At the annealing stage the diffusion of free PCR product occures in parrallel with its hybridisation. We represent this by means of alternating steps of diffusion and hybridisation during the annealing stage: first we calculate f after one second diffusion, Daf(x,y); and then accoun for hybridisation of free PCR product. That is, fiaj(x,y)=(1-h1(x,y))Dafiaj-1(x,y) for 1jA, where A is the duration of annealing stage (in seconds); (x,y) equals 1 if the point (x,y) belongs to the gel pad (microarray element), and equals 0 otherwise; h1[0;1] is the hybridisation parameter (a portion of free PCR product hybridised to solid phase primers for one second); together (1-h1(x,y)) comprise the hybridisation operator. The operator Da is similar to Dm: Daf(x,y)exp(-(x2+y2)/ra2)*f(x,y), where ra is the diffusion distance per one second of annealing stage. The annealing of free PCR product to solid phase forward primers increases the parameter u (by the value equal to decrease of f): uiaj(x,y)=uiaj-1(x,y)+h1(x,y)Dafiaj-1(x,y) for 1jA, Similarly, free reverse primer anneals to extended forward primer resulting in decrease of g: giaj(x,y)=(1-h2(x,y))giaj-1(x,y) for 1jA, where h2[0;1] is the portion of extended forward primer hybridised per one second of annealing stage. It is assumed here that the consumption of reverse primer at exponential phase of PCR is insignificant in comparison with its total amount. The annealing of reverse primer increases the parameter v: viaj(x,y)=viaj-1(x,y)+h2(x,y)gim(x,y) for 1jA, Elongation stage Following the elongation stage, the amount of exteded forward primers (g) grows again due to (i) the extention of solid phase forward primers with annealed free PCR products (contribution of parameter u) and (ii) the release of extended reverse primers from immobilised PCR product (contribution of parameter v, see Figure I, Elongation stage): gie(x,y)=giaA(x,y)+uiaA(x,y)+viaA(x,y). The amount of free PCR product changes similarly, however, here the diffusion should be taken into account: fie(x,y)=uiaA(x,y)+viaA(x,y)+DeEfiaA(x,y) , where E is the duration of elongation stage, De is defined as follows: Def(x,y)=exp(-(x2+y2)/re2)*f(x,y), where re is the diffusion distance per one second of elongation stage. Capillary microarray (one dimensional model) For one dimensional model the f, g, u, v, and are the functions of one variable (x only); D, Dm, Da, and De are the operators of one dimensional convolution. Calculation of convolution The calculation of convolution with fundamental solution of the diffusion equation is outlined at Figure II. The f(x,y) and fundamental solution [Cexp(-(x2+y2)/r2)] are presented as (2n+1)(2n+1) and (2m+1)(2m+1) arrays of floating point numbers, respectively (C is a normalisation constant). Here the elements (n+1,n+1) and (m+1,m+1) correspond to (x,y)=(0,0). A B f(x,y) gel pad wall Cexp(-(x2+y2)/r2) D C f1(x,0) f1(x,y)=f1((x2+y2)1/2,0) Figure II. Scheme of calculation of convolution. For calculation of convolution f1(x,y)=Cexp(-(x2+y2)/r2)*f(x,y) it was assumed that the reaction volume was surrounded with impenetrable wall of (2n+1) diameter (Figure II A, red line). The (2n+1)(2n+1) array was duplicated producing the (4n+2)(2n+1) array (Figure II B). Then (1+n+m) subarrays of (2m+1)(2m+1) size (Figure II B, red dashed square frames) centered on (x,0) line (Figure II B, blue frame) were consequitively selected. The elements of these subarrays were multiplied with corresponsing elements of fundamental solution array and the sums of these products were taken. These sums form 1(1+n+m) array corresponding to f1(x,0) (Figure II C). For other positions (x,y) the convolution was calculated as follows: f1(x,y)=f1((x2+y2)1/2,0) (Figure II D). This procedure assumes circular simmetry of f(x,y) (including circular simmetry of (x,y) as a prerequisite) and the condition of n >> m. ADDITIONAL RESULTS Estimation of diffusion and hybridisation parameters For comparison with experimental data the two dimensional model was executed with M=20, A=45 (standard condition), M=20, A=22 (short annealing condition), and M=40, A=45 (long melting condition). For the reasons of simplicity it was assumed that h1=h2h. For every condition the model was executed with the same combinations of rm, ra and h; n=150, m=15, (x,y)=1 if (x2+y2)1/2R=30 (i.e. 10m/pixel) and (x,y)=0 otherwise. It was found that for ra=0 the profile of extended solid phase primer amount across the gel pad (microarray element) gie(x,0) was identical to that of simple model described in the main text (Figure 4A therein). However, for ra>0.48÷0.5 (depending on h) the maximum of gie was achieved at the border of gel pad (Figure III A) rather than at its centre (as it was shown experimentally; main text, Figure 1A). For concordance with experimental data the amplification factor should be from 1.51 to 1.67; the CT difference between the standard and the short annealing conditions should be from 3 to 5. For ra=0 this corresponds to rm between 5 and 7.5 and h between 0.08 and 0.1 (Figure III B). A 1,3 ra=0 . B 6,5 1,2 ra=0.5. 6 ra=0.6. h=0.07 5,5 1,1 5 1 signal . CT 4,5 D CT 0,9 h=0.08 4 0,8 3,5 h=0.09 0,7 3 0,6 h=0.1 2,5 rm= 7.5 7.0 6.5 6.0 5.5 5.0 0,5 2 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 1,5 1,52 1,54 1,56 1,58 1,6 1,62 1,64 1,66 1,68 1,7 position amplification factor C 5,5 5 max in the centre max at the border 4,5 h=0.08 . 4 D CT CT 3,5 h=0.09 3 h=0.1 ra=0.5 ra=0.6 2,5 ra=0 ra=0.48 2 1,56 1,565 1,57 1,575 1,58 1,585 1,59 amplification factor Figure III. Estimation of diffusion and hybridisation parameters rm, ra, and h. A. The cross-section of gel pad for rm=6.5, h=0.08. Black squares - ra=0, red squares - ra=0.5, green diamonds - ra=0.6. The latter line achieves its maximum at the border rather than in the centre of gel pad. B. The difference of threshold cycle value (CT) between reactions with 22s and 45s annealing time versus the amplification factor with 45s annealing time for ra=0 and variable rm and h. Blue lines are the curves of constant h, red lines are the curves of constant rm. Dashed frame outlines the area of experimentally obtained values of CT and apmlification factor. C. The same as the pannel B, but for constant rm=6.5 and variable h and ra (blue and red lines, respectively). Solid squares correspond to the conditions where the signal maximum is achieved in the centre of gel pad, open squares correspond to the conditions where the signal maximum is achieved at the border of gel pad (like green diamond line at pannel A). To study how these conditions could be affected by ra, the numerical simulation data are plotted for rm=6.5 in Figure III C. The increase of ra from 0 to 0.5 leads to insignificant (within 0.01) increase of amplification factor and small (within 0.3) increase of CT difference between standard and short annealing conditions. That is, for concordance with experimental data one should assume that rm ranges from 5 to 7.5, h ranges from 0.08 to 0.1, ra ranges from 0 to 0.5. The effect of reaction chamber radius One can assume that the small size of reaction chamber could lead to the hindered diffusion of free PCR product during the melting stage which in turn could lead to the faster growth of PCR curve with the lower CT. To verify this hypothesis the simple model (accounting for the melting stage only) was executed with n ranging from 300 to 30, m=2 (to keep the condition n>>m), rm=1, M=900, R=30.15. These combinations of parameters were ajusted to obtain the amplification factor of 1.59 (equals to the experimental value, see main text) for n150. A B The wall R=30 R=30.15 Figure IV. The comparison of R=30 (A), and R=30.15 (B) conditions. The reason of choosing R=30.15 instead of R=30 is outlined on Figure IV. On this figure the pixels belonging to the gel pad and its mirror image in the duplicated (2n+1)(2n+1) array are marked with black and grey colour, respectively for R=30 (Figure IV A) and R=30.15 (Figure IV B). The pixels used for convolution calculation for the position close to the reaction chamber wall are surrounded with the red frame. In the Figure IV B all the pixels in the red frame are either balck or grey, i.e. belong either to the gel pad or to its image, while in Figure IV A there are some white pixels in the red frame. That is, for n around 30 the Figure IV B provides more natural environment than the Figure IV A. 31 30 2D 29 1D 28 27 CT 26 25 24 23 22 21 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 R/n Figure V. Threshold cycle value for different ratios of gel pad radius (R) to the dictance from the centre of gel pad to the wall of hybridisation chamber (n). Solid circles – two dimensional model, open triangles – two dimensional model. The plot of CT versus R/n ratio is given in Figure V. For 0<R/n<1/2 the CT value changes insignificantly for the both one dimensional and two dimensional models, while for 1/2<R/n<1 it decreases rapidly to its minimal value. That is, for significant effect on PCR the size of reaction chamber should be reduced below the double size of gel pad. This finding could be important for optimisation of emulsion PCR where the sepharose bead with immobilised primers is located in the droplet of PCR solution immersed in the mineral oil.
Pages to are hidden for
"NUMERICAL MODELING"Please download to view full document