May 5, 2020

Our course so far

  • Solve simple problems
    • Decomposition
    • Pattern matching
    • Abstraction
  • Systems Analysis
    • Identify the system parts and processes
    • Use R to simulate a system
    • See what happens in the long term

Long term

One important question that we want to answer is

What happens in the long term

In may cases it is not easy to guess. Our intuition is bad

We need to simulate and see what happens

Polymerase Chain Reaction

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
    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))

Simulate several PCR cycles

pcr(N=6, DNA_ini=1e6, primer_ini=1e8, cycle_rate=1e-8)
  cycle      DNA    primer
1     1  1000000 100000000
2     2  2000000  99000000
3     3  3980000  97020000
4     4  7841396  93158604
5     5 15146331  85853669
6     6 28150012  72849988

Simulate several PCR cycles

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

Long term behavior

  • The system reaches an equilibrium
  • Once the primers are finished, the reaction stops
  • We can see the final DNA and primer concentrations
  • Question: Does final DNA depends on
    • initial DNA concentration?
    • initial primer concentration?
    • PCR reaction rate?

Effect of initial values

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)
[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) {
    experiment <- pcr(N=30, DNA_ini=initial_DNA[i],
                      primer_ini=1e8, cycle_rate=1e-8)
    plot(DNA ~ cycle, experiment, main=initial_DNA[i])

Result in next slide

PCR depends on initial concentration

All in one plot

Remember that each plot() command creates a new plot

We want to have only one plot

We will plot() the first result, the we will add lines() for the rest

PCR depends on initial concentration

experiment <- pcr(N=30, DNA_ini=initial_DNA[1],
                  primer_ini=1e8, cycle_rate=1e-8)
plot(DNA ~ cycle, data=experiment, type="o",
            main="Effect of initial DNA on PCR")
for(i in 2:length(initial_DNA)) {
    experiment <- pcr(N=30, DNA_ini=initial_DNA[i],
                      primer_ini=1e8, cycle_rate=1e-8)
    lines(DNA ~ cycle, experiment, pch=i, type="b")
legend("topleft", legend=initial_DNA, pch=1:7)

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 measure 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 can use the function position_of_half() from Homework 6.

position_of_half() function

one solution for Homework 6

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


DNA concentration crosses 50% at 21 cycles

Cycle Threshold

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


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

Water formation


A chemical reaction combines hydrogen and oxygen to produce water

To make it simple we will assume that the reaction is this:

2 H + OH2O

Our system has 3 items:

  • Hydrogen
  • Oxygen
  • Water

Reactions are bidirectional

The items of the system

We represent hydrogen by H, oxygen by O, and water by W.

Each vector has the number of atoms or molecules for different times

We will measure the number of atoms and molecules in moles

The initial quantities are:

  • 1 mole oxygen,
  • 1 mole hydrogen, and
  • 0 mole water

The rules are the following:

There are two reactions at the same time:

  1. from hydrogen and oxygen into water
    (that is, left to right),

  2. from water into hydrogen and oxygen
    (that is, right to left).

We also see that

  • The number water molecules increase with the first reaction and decreases with the second reaction
  • The number of oxygen atoms reduces with the first reaction and increases with the second
  • The number of hydrogen atoms reduces two times with the first reaction and increases two times with the second

Simplifying notation

According to our rules, we have two rates:

  • r1_rate
  • r2_rate

In chemistry, these rates are usually written as k1 and k2

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
  • actual reaction rates depend on temperature and other conditions

Main part of simulation

for(i in 2:N) {
    r1 <- r1_rate*H[i-1]*H[i-1]*O[i-1]
    r2 <- r2_rate*W[i-1]
    d_W[i] <-   r1 - r2
    d_O[i] <-  -r1 + r2
    d_H[i] <- -2*r1 + 2*r2
    W[i] <- W[i-1] + d_W[i]
    O[i] <- O[i-1] + d_O[i]
    H[i] <- H[i-1] + d_H[i]

Exercise: finish it

Create a water_formation function.

Inputs are:

  • N: Number of steps in the simulation
  • H_ini: initial amount of hydrogen
  • O_ini: initial amount of oxygen
  • W_ini: initial amount of water
  • r1_rate: reaction 1 rate
  • r2_rate: reaction 2 rate

Exercise 2: conservation of mass

One atom of oxygen has 16 times the mass of one atom of hydrogen

Therefore one molecule of water has mass 18

Show that the total mass never changes

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

… 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

Extra material