# NUMERICAL MODELING by hcj

VIEWS: 7 PAGES: 5

• pg 1
```									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 1jA,
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 1jA,
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 1jA,
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 1jA,

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 [Cexp(-(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)

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)=Cexp(-(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.

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=h2h. 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/2R=30 (i.e. 10m/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 n150.
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.

```
To top