In [67]:
# define the names of the nodes (molecules) in the model
nodes <- c("Cell Size",
"Cln3",
"MBF",
"Clb5,6",
"Mcm1/SFF",
"Swi5",
"Sic1",
"Clb1,2",
"Cdc20&Cdc14",
"Cdh1",
"Cln1,2",
"SBF")
# count the number of nodes in the network
N <- length(nodes)
# define an NxN matrix of dependencies (these are the "edges" in the network)
# the row (first index) identifies the "regulator", and the column
# (second index) identifies the target molecule.
# NOTE: I define aij as the transpose of how Li et al. define theirs!
a <- matrix(nrow=N, ncol=N)
# label the rows and columns of the matrix by the node names
rownames(a) <- nodes
colnames(a) <- nodes
# initialize all entries of the matrix to zero
a[,] <- 0
# define the activation edges (green arrows)
a["Cell Size", "Cln3"] <- 1
a["Cln3", "MBF"] <- 1
a["Cln3", "SBF"] <- 1
a["SBF", "Cln1,2"] <- 1
a["MBF", "Clb5,6"] <- 1
a["Cdc20&Cdc14", "Sic1"] <- 1
a["Clb5,6", "Mcm1/SFF"] <- 1
a["Clb5,6", "Clb1,2"] <- 1
a["Mcm1/SFF", "Clb1,2"] <- 1
a["Clb1,2", "Mcm1/SFF"] <- 1
a["Mcm1/SFF", "Cdc20&Cdc14"] <- 1
a["Mcm1/SFF", "Swi5"] <- 1
a["Cdc20&Cdc14", "Swi5"] <- 1
a["Cdc20&Cdc14", "Cdh1"] <- 1
a["Clb1,2", "Cdc20&Cdc14"] <- 1
a["Swi5", "Sic1"] <- 1
# define the repression edges (red arrows)
a["Clb1,2", "SBF"] <- -1
a["Clb1,2", "MBF"] <- -1
a["Clb1,2", "Cdh1"] <- -1
a["Clb1,2", "Swi5"] <- -1
a["Clb1,2", "Sic1"] <- -1
a["Cln1,2", "Sic1"] <- -1
a["Clb5,6", "Sic1"] <- -1
a["Sic1", "Clb5,6"] <- -1
a["Cln1,2", "Cdh1"] <- -1
a["Clb5,6", "Cdh1"] <- -1
a["Cdc20&Cdc14", "Clb5,6"] <- -1
a["Cdc20&Cdc14", "Clb1,2"] <- -1
a["Sic1", "Clb1,2"] <- -1
a["Cdh1", "Clb1,2"] <- -1
# define the autorepression edges (yellow loops)
a["Cln3", "Cln3"] <- -1
a["Swi5", "Swi5"] <- -1
a["Cdc20&Cdc14", "Cdc20&Cdc14"] <- -1
a["Cln1,2", "Cln1,2"] <- -1
a["Mcm1/SFF", "Mcm1/SFF"] <- -1
# define the fixed points of the model
num.fp <- 7
fixed.points <- matrix(nrow=num.fp,
ncol=N)
fixed.points[,] <- 0
colnames(fixed.points) <- nodes
# define the nonzero entries of the fixed-point matrix
fixed.points[1,"Cdh1"] <- 1
fixed.points[1,"Sic1"] <- 1
fixed.points[2,"SBF"] <- 1
fixed.points[2,"Cln1,2"] <- 1
fixed.points[3,"MBF"] <- 1
fixed.points[3,"Cdh1"] <- 1
fixed.points[3,"Sic1"] <- 1
fixed.points[4,"Sic1"] <- 1
fixed.points[5,"MBF"] <- 1
fixed.points[5,"Sic1"] <- 1
fixed.points[7,"Cdh1"] <- 1
# set up the vector for recording the number of times you end up at each
# of the different fixed points
basin.counts <- vector(length=num.fp, mode="integer")
basin.counts[] <- 0
Define a function hamming.dist that gives the hamming distance between two states of the Boolean network
In [68]:
hamming.dist <- function(x1,x2) {
# return the Hamming distance between state vectors x1 and x2
return(sum(abs(x1-x2)))
}
Define a function evolve that takes the network from one Boolean vector state to another Boolean vector state
In [69]:
evolve <- function(state) {
# "state" is a vector of length N, containing zeros or ones.
# A zero for a given node means it is in the "off" state; a
# one means that the node is in the "on" state. Update for one
# time-step according to the instructions in Eq. 1 of the Li et al. paper,
# and return the new state
# compute the influence function as product of state and the a matrix
result <- state %*% a
# positive entries get set to 1
result[result > 0] <- 1
# for any node for which the influence function is zero, use the current state
result[result == 0] <- state[result == 0]
# negative entries get set to 0
result[result < 0] <- 0
return(result)
}
Write a function that runs 10,000 simulations of the network. In each simulation, the procedure is:
state (make sure the first position of the state vector is set to zero)state to evolve and then updating state with the return value from evolvestate changes in the last call to evolve; if it does not, then you have reached a fixed point; stop iteratingfixed.points; for the unique row i for which you find a match, increment the element in position i of basin.countsbasin.counts
In [70]:
basin_ids <- replicate(10000, {
state <- sample.int(2, N, replace=TRUE)-1
state[1] <- 0
state.new <- NULL
while(1) {
state.new <- evolve(state)
if (0 == hamming.dist(state, state.new)) {
break()
}
state <- state.new
}
fpi <- which(0 == apply(fixed.points, 1, function(fixedpoint_vec) {hamming.dist(fixedpoint_vec, state.new)}))
})
table(basin_ids)