Polymerase Chain Reaction

The Polymerase Chain Reaction (PCR) is a method used to synthesize millions of copies of a given DNA sequence.

A typical PCR reaction consists of series of cycles:

• template DNA denaturation,
• primer annealing, and
• extension of the annealed primers by DNA polymerase.

This loop is repeated between 25 and 30 times

PCR system

We simplify and we forget about the polymerase and the dNTP. They both will be represented by primers.

We have a system with two parts, DNA and primers and one process, the thermal cycle.

The system is represented by this diagram:

We have seen this before

This system was in Homework 8

Many of you noticed that this system has the same shape of cell-food system

Same shape means same behavior

Simulate several PCR cycles

Let’s write a function to simulate the PCR reaction.

The function name should be pcr and it must take four inputs:

• the number of cycles N,
• the initial DNA concentration DNA_ini,
• the initial primer concentration primer_ini, and
• the reaction rate cycle_rate.

The function must return a data frame with cycle, DNA and primer.

Simulate several PCR cycles

pcr <- function(N, DNA_ini, primer_ini, cycle_rate) {
    DNA    <- d_DNA    <- rep(NA, N)
    primer <- d_primer <- rep(NA, N)
    DNA[1] <- DNA_ini
    primer[1] <- primer_ini
    d_DNA[1] <- d_primer[1] <- 0
    for(i in 2:N) {
        cycle <-  cycle_rate * primer[i-1] * DNA[i-1]
        d_DNA[i]    <-  cycle
        d_primer[i] <- -cycle
        DNA[i]    <-    DNA[i-1] + d_DNA[i]
        primer[i] <- primer[i-1] + d_primer[i]
    }
    return(data.frame(cycle=1:N, DNA, primer, d_DNA, d_primer))
}

Simulate several PCR cycles

pcr(N=6, DNA_ini=1e6, primer_ini=1e8, cycle_rate=1e-8)

  cycle      DNA    primer    d_DNA  d_primer
1     1  1000000 100000000        0         0
2     2  2000000  99000000  1000000  -1000000
3     3  3980000  97020000  1980000  -1980000
4     4  7841396  93158604  3861396  -3861396
5     5 15146331  85853669  7304935  -7304935
6     6 28150012  72849988 13003681 -13003681

Simulate several PCR cycles

exprmnt <- pcr(N=20, DNA_ini=1e6, primer_ini=1e8,
                               cycle_rate=1e-8)
plot(DNA ~ cycle, exprmnt, type="o")

Question

Does final DNA depends on

• initial DNA concentration?

• initial primer concentration?

• PCR reaction rate?

PCR depends on initial concentration

The PCR reaction curve depends on the initial concentration of DNA.

We want to understand this dependency for the following values of initial DNA concentration:

initial_DNA <- 10**(0:6)
initial_DNA

[1] 1e+00 1e+01 1e+02 1e+03 1e+04 1e+05 1e+06

PCR depends on initial concentration

for(i in 1:6) {
  exprmnt <- pcr(N=30, DNA_ini=initial_DNA[i], primer_ini=1e8, cycle_rate=1e-8)
  plot(DNA ~ cycle, exprmnt, main=initial_DNA[i])
}

PCR depends on initial concentration

Quantitative PCR

• Long term behavior is the same in all cases
  • DNA concentration is saturated
  • All primers are used
• Initial DNA concentration changes when the saturation happens
  • More initial DNA means faster reaction

Initial DNA concentration affects reaction time

• We need to find when DNA concentration increases a lot
• One easy way to do that is to find in which cycle does the DNA concentration crosses the 50% line
• We need a function to find position_of_half()

How would you write that function?

position_of_half() function

one possible solution

position_of_half <- function(x) {
    half_value <- (max(x)+min(x))/2
    for(i in 1:length(x)) {
        if(x[i]>= half_value) {
            return(i)
        }
    }
}

Example

DNA concentration crosses 50% at 21 cycles

Cycle Threshold

CT <- rep(NA, length(initial_DNA))
for(i in 1:length(initial_DNA)) {
    exprmnt <- pcr(N=30, DNA_ini=initial_DNA[i],
                    primer_ini=1e8, cycle_rate=1e-8)
    CT[i] <- position_of_half(exprmnt$DNA)
}

or, better

CT <- sapply(1:length(initial_DNA), function (i) {
    exprmnt <- pcr(N=30, DNA_ini=initial_DNA[i],
                    primer_ini=1e8, cycle_rate=1e-8)
    return(position_of_half(exprmnt$DNA))
})

Result

We conclude that the CT value can be used to predict the initial DNA concentration.

Bidirectional Reactions

We can see that

• The speed of reaction r1 depends on the reaction rate, the number of oxygen atoms, and the square of the number of hydrogen atoms

  r1_rate *  H * H * O

• The speed of r2 depends on the reaction rate and the number of water molecules

  r1_rate *  W

• The actual reaction rates depend on the temperature and other conditions

Long term behavior: equilibrium

mix <- water_formation(N=200, r1_rate=0.01, r2_rate=0.01,
     H_ini=2, O_ini=1, W_ini=0)

Equilibrium values depends on r1_rate

rates_of_r1 <- 10^seq(from=-3, to=-1, length=5)
r2_rate=0.01, H_ini=2, O_ini=1, W_ini=0

We can predict final outcome for each rate

rates_of_r1 <- 10^seq(from=-3, to=-1, length=5)
r2_rate=0.01, H_ini=2, O_ini=1, W_ini=0

Same with more r1_rate values

rates_of_r1 <- seq(from=0.001, to=0.1, length=50)
r2_rate=0.01, H_ini=2, O_ini=1, W_ini=0

This is the same pattern as CT on qPCR

final_Water <- rep(NA, length(rates_of_r1))
for(i in 1:length(rates_of_r1)) {
    mix <- water_formation(N=100, r1_rate=rates_of_r1[i],
                r2_rate=0.01, H_ini=2, O_ini=1, W_ini=0)
    final_Water[i] <- mix$Water[100]
}
plot(rates_of_r1, final_Water)

What was the reaction rate?

We can measure how much water we got.
If at the end we have 0.5 mol of water, then …

… with the curve we can find what was the reaction rate

Other Long-Time Behavior

The two cases we saw had the same long-term behavior: equilibrium

Other systems (such as predator-prey model) can have other behavior: cycle

Key Questions

• What is the long time behavior
  • Equilibrium
  • Cycle
  • Chaos
• What is the long-term effect of
  • initial values
  • process rates
• If we know the long-term behavior, what was
  • the initial values
  • the process rates