March 22, 2019
for(){}
, while(){}
if(){}
, if(){} else {}
fac(N) = N * fac(N-1)
plot()
, segments()
read.fasta()
Biology is the study of life
Life is a complex phenomenon
To understand life, biologists invented 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
A system is a group of parts that are interconnected and affect each other
Last class we saw a system with two items
and one process
We will represent any system with a network with two kinds of nodes
Circles are connected to boxes with arrows and vice-versa
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
The arrows going into a box show which items are required for the process
We want to separate eating from duplication
There are two processes and three items
Draw this model with pen and paper
At any given moment, there is some amount of each item in the system
Also, we can know how much did the amount changed between the previous moment and this moment
Each item has two values, that change with time
We will represent these values with vectors
item
d_item
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”
cells[i]
food[i]
In any fixed time, the state is all the values of items’ quantities
cells[i]
, food[i]
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 new state depends only on the past states, nothing else
If we know the initial state and rate constants, we can calculate everything
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
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”
To get the formula for each circle
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
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
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
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]
Write the formulas for the eating-dupication system
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.
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] }
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] }
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)) }
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?
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:
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
There are two reactions at the same time:
from hydrogen and oxygen into water (that is, left to right),
from water into hydrogen and oxygen (that is, right to left).
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] }
Create a water_formation
function.
Inputs are:
N
: Number of steps in the simulationH_ini
: initial amount of hydrogen, default 1O_ini
: initial amount of oxygen, default 1W_ini
: initial amount of water, default 0k1
: reaction 1 rate, default 0.1k2
: reaction 2 rate, default 0.1One 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
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.