March 22, 2019

Are you ready?

What do you need to know?

  • Computational Thinking
    • Decomposition
    • Pattern Recognition
    • Abstraction
    • Algorithm Design
  • Not memory
  • Problem solving

Key ideas

  • Loops: for(){}, while(){}
  • Conditional: if(){}, if(){} else {}
  • Recursive functions:
    • fac(N) = N * fac(N-1)
  • Cumulative values: today = yesterday + change

Practical tools

  • Turtle graphics
  • Normal graphics: plot(), segments()
  • Reading DNA sequences read.fasta()

The most important part

Do you want to pass this class?

Are you lazy?

Are you a hard worker?

If you are here, you are smart

Do you have the will to succeed?

Did you do all the homework?

Google is your friend

Did you find new exercises?

Systems

Today we will see definitions

Be sure of understand them

What is life?

Biology is the study of life

Life is a complex phenomenon

To understand life, biologists invented Systems Theory

Systems Theory

Systems Theory was invented in 1948 by the biologist Ludwig von Bertalanffy

Today it is used by engineers, psychologists, managers, and other scientists

It is going to be one of the important ways of thinking in this century

Systems

A system is a group of parts that are interconnected and affect each other

  • The parts are separated and well defined
  • Instead of parts, we can say items or components
  • They have interactions, for example a process that transforms some items into others

Example

Last class we saw a system with two items

  • Cells
  • Food

and one process

  • Eating, which also was cell duplication

Drawing Systems

We will represent any system with a network with two kinds of nodes

  • Circles will represent the items
  • Boxes will represent their interactions

Circles are connected to boxes with arrows and vice-versa

Arrows

  • Circles are connected to boxes
  • Boxes are connected to circles

  • There cannot be arrows between circles or between boxes
    (only between nodes of different shape)

  • There can be more than one arrow between a box and a circle, in any direction

The technical name of this kind of network is directed bipartite multigraph

Example: Molecular Biology’s Dogma

  • Which are the items?
  • Which are the processes?

How to draw the system

The arrows going into a box show which items are required for the process

  • The items that point towards a process are the inputs that will be consumed by the process
    • When a item is consumed by a process, its amount gets smaller
  • The arrows going out of a box show which are the products of the process
    • The amount of these items will increase with the process.

Exercise: eating-dupication system

We want to separate eating from duplication

There are two processes and three items

  • food
  • cells
  • fat cells

Draw this model with pen and paper

Each item has two values

At any given moment, there is some amount of each item in the system

  • for example, the number of cells

Also, we can know how much did the amount changed between the previous moment and this moment

  • for example, the growth of cells

Representing items in R

Each item has two values, that change with time

We will represent these values with vectors

  • The amount of item will be called item
  • The change of item will be called d_item
    • We read it “delta-item”

These are dynamical systems

Dynamical system means that the system changes with time

We will represent time by a positive integer number

Time i represents “now”. The units can be seconds, hours, days, years, or whatever is best to understand the system

If we use “days”, then i-1 means “yesterday”. It can also be “last year” or “one second ago”

Systems are in a state

  1. The particular condition that someone or something is in at a specific time
    • water in a liquid state.
  2. The amounts of all items (circles) in a system.
    • cells[i]
    • food[i]

The state of the system changes with time

In any fixed time, the state is all the values of items’ quantities

  • one value for each item
  • Examples:
    • cells[i], food[i]
    • Concentration of H, O, Water
    • Number of foxes and rabbits
    • Number of primers and DNA molecules

How items change

The “boxes” of the system represent processes

They do not change with time

The have one constant value called rate

For each process we have a value process_rate

The state today depends only on the past states

The new state depends only on the past states, nothing else

If we know the initial state and rate constants, we can calculate everything

Finding the formulas for boxes

For each box in the graph we get only one term

We multiply the rate constant and each of the amount variables of the circles connected by incoming arrows

We use the index i-1. Processes depend on the previous values

If there are several incoming arrows from the same circle, then the variable is multiplied several times

The outgoing arrows are not important in this part

Example

There are two arrows coming into the eating box

The formula for this box is

eating_rate * cells[i-1] * food[i-1]

We use the index i-1 because we do not know yet the value of cells[i] or food[i]

We will calculate them later

At the begin of “today” we only know “yesterday”

Formulas for circles

To get the formula for each circle

  • sum all the terms of boxes connected with incoming arrows,
  • minus all the terms of boxes connected with outgoing arrows

This value is assigned to the delta variable of the circle.

In our example the formula for delta_food is

d_food[i] <-  -eating_rate * cells[i-1] * food[i-1]

since the food circle has only one outgoing arrow

Formula for the circle (cont…)

The cell circle has two incoming arrows
from eat and one outgoing arrow to eat. Therefore

d_cells[i] <-   eating_rate * cells[i-1] * food[i-1] 
              + eating_rate * cells[i-1] * food[i-1] 
              - eating_rate * cells[i-1] * food[i-1]

which, after simplification, is just

d_cells[i] <- eating_rate * cells[i-1] * food[i-1]

resulting on a final result of one positive incoming arrow

If there are several arrows

In the last slide we had several arrows between the cells circle and the eating box

It is easy to see that we only care about the resulting number and direction of arrows

Two input - One output = One input

In summary

Now we can write the formulas
for all the delta variables

d_food[i]  <- -eating_rate * cells[i-1] * food[i-1]
d_cells[i] <-  eating_rate * cells[i-1] * food[i-1]

Finally, the amount variables have to be updated

Each amount variable is the cumulative sum of the delta variables

food[i]  <- food[i-1]  + d_food[i]
cells[i] <- cells[i-1] + d_cells[i]

Exercise

Write the formulas for the eating-dupication system

  • First write the terms for each box (process)
  • Then write the formulas for each circle (item)

Initial conditions

The last missing piece are the initial values of the circles’ variables.

The value of cells[i] depends on the value of cells[i-1], and we can only calculate that for i >= 2

We will use the variables cells_ini and food_ini to define the values used in cells[1] and food[1]

For the delta variables, we can assume that they are initially zero.

Last class we wrote this

ncells <- growth      <- rep(NA, 100)
food   <- food_change <- rep(NA, 100)
ncells[1]      <- 1
growth[1]      <- 0
food[1]        <- 20
food_change[1] <- 0
eat_rate <- 0.003
for(i in 2:100) {
  growth[i]      <-  ncells[i-1]*food[i-1]*eat_rate
  food_change[i] <- -ncells[i-1]*food[i-1]*eat_rate
  ncells[i] <- ncells[i-1] + growth[i]
  food[i]   <-   food[i-1] + food_change[i]
}

Using the new names

N <- 100
cells   <- d_cells <- rep(NA, N)
food    <- d_food  <- rep(NA, N)
cells[1] <- 1
food[1]  <- 20
d_cells[1] <- d_food[1]  <- 0

for(i in 2:N) {
  d_cells[i] <- cells[i-1]*food[i-1]*eating_rate
  d_food[i]  <- -cells[i-1]*food[i-1]*eating_rate
  cells[i]  <- cells[i-1] + d_cells[i]
  food[i]   <- food[i-1]  + d_food[i]
}

Algorithm design: write as function

cell_culture <- function(N, eating_rate, cells_ini, food_ini) {
  cells   <- d_cells <- rep(NA, N)
  food    <- d_food  <- rep(NA, N)
  cells[1] <- cells_ini
  food[1]  <- food_ini
  d_cells[1] <- d_food[1]  <- 0
  for(i in 2:N) {
    d_cells[i] <- cells[i-1]*food[i-1]*eating_rate
    d_food[i]  <- -cells[i-1]*food[i-1]*eating_rate
    cells[i]  <- cells[i-1] + d_cells[i]
    food[i]   <- food[i-1]  + d_food[i]
  }
  return(data.frame(cells, food))
}

Exercise

Write the code to simulate the eating-dupication system

How does the behavior change with different initial values?

How does the behavior change with different process rates?

Another system: Water

A chemical reaction combines hydrogen and oxygen to produce water

To make it simple we will assume that the reaction is this: \[2H+O\leftrightarrow H_2O\]

Our system has 3 items:

  • Hydrogen
  • Oxygen
  • Water

The simple system

Reactions are bidirectional

The items of the system

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

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

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

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

More rules

  • The speed of reaction1 depends on the number of oxygen atoms and the square of the number of hydrogen atoms \[r_1 = k_1 H^2O\]
  • The speed of reaction2 depends on the number of water molecules \[r_2=k_2 W\]
  • The constants \(k_1\) and \(k_2\) depend on temperature and other conditions

More rules

  • The number water molecules increase with the first reaction and decreases with the second reaction \[d\_W = r1 - r2\]
  • The number of oxygen atoms reduces with the first reaction and increases with the second \[d\_O = -r1 + r2\]
  • The number of hydrogen atoms reduces two times with the first reaction and increases two times with the second \[d\_H = 2*(-r1 + r2)\]

The reactions

Then the main part of the simulation is

for(i in 2:N) {
    d_W[i] <-    k1*H[i-1]*H[i-1]*O[i-1] - k2*W[i-1]
    d_O[i] <-   -k1*H[i-1]*H[i-1]*O[i-1] + k2*W[i-1]
    d_H[i] <- -2*k1*H[i-1]*H[i-1]*O[i-1] + 2*k2*W[i-1]

    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, default 1
  • O_ini: initial amount of oxygen, default 1
  • W_ini: initial amount of water, default 0
  • k1: reaction 1 rate, default 0.1
  • k2: reaction 2 rate, default 0.1

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

References

Ludwig von Bertalanffy, “General System Theory: Foundations, Development, Applications” New York (1968)

Koch, Ina. “Petri Nets - A Mathematical Formalism to Analyze Chemical Reaction Networks.” Molecular Informatics 29, no. 12 (December 17, 2010): 838–43. doi:10.1002/minf.201000086.

Baez, John C., and Jacob Biamonte. “Quantum Techniques for Stochastic Mechanics,” September 17, 2012. http://arxiv.org/abs/1209.3632.