Tutorial of Boolean Network Analysis of Time Series Data
Written by Xiao Yang
last updated Aug, 2023
Overview
Dynamical system methods have been instrumental in explaining mechanisms of change and predicting future trajectories and states of psychological systems (Boker & Graham, 1998; Chow, Ferrer, & Hsieh, 2010). A natural next step is to use these methods to design control strategies and guide system trajectories (Molenaar & Nesselroade, 2015; Wang et al., 2014).
In this tutorial, we introduce and forward a Boolean network method because it can give an intuitive explanation of psychological dynamics that are useful for theorists and suggestions for network control that are useful for practitioners. We demonstrate the utility of the Boolean network through analysis of multivariate binary time-series data simulated to mimic intensive longitudinal data (e.g., ecological momentary assessment).
Boolean network
Boolean network (Kauffman, 1969, 1993) method is a specific instantiation of discrete dynamical system model. The special characteristic of the Boolean network method is to use the Boolean operators, including AND (&), OR (|), or NOT (!), to describe and model temporal dynamics between variables.
The Boolean network method is like a vector-autoregression (VAR) model, whereas the temporal relations are expressed using the Boolean operators, and the variables are binary.
Data type that can use the Boolean network method: The Boolean netowrk method can be applied on intensive longitudinal data (e.g., Ecological Momentary Assessment, physiological data, physical activity) to model the temporal dynamics and facilitate control design. The method can work on binary time-series, and continuous-scale time-series, which need to be binarized (see Step 0).
We are using the R package BoolNet (Mussel et al., 2010) to conduct the data analysis in this tutorial.
This tutorial is to demonstrate the following 6 steps of Boolean network approach, including:
Binarization of continuous data
Simulation of binary data
Inference of Boolean functions
Plotting state space transition graph
Extraction of attractors
Design of network control
library(ggplot2) # for data visualization
library(BoolNet) # for Boolean network modeling
library(dplyr) # for data organization ##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(igraph) # for network plotting ##
## Attaching package: 'igraph'## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionlibrary(stringdist) # for computing distance between strings/sequences
library(stringr) # for computing distance between strings/sequences
library(reshape2) # for data organization
library(metafolio) # for color specification
set.seed(12345)Step 0: Binarization of continuous data
Binarization is to convert the continuous-scale variables to binary variables, so that it can be fitted into the Boolean model. This is done using the “binarizeTimeSeries” function in BoolNet package, and note here the binarization is done by rows, instead of by columns. Since we usually use the long-format to organize intensive longitudinal data, we take the transpose here to transform the data into a format that each row is a variable.
The “binarizeTimeSeries” function is using k-means clustering and k (number of groups) = 2.
Here we show an example of how binarization works on continuous-scale data. You can plug in your continuous-scale data in this step, to get binarized data.
n.obs <- 20 # set the length of the time-series
# simulate a 3-variate continuous-scale time-series using a uniform distribution between 0 and 1.
continuous.time.series <- cbind(runif(n.obs,0, 1),
runif(n.obs,0, 1),
runif(n.obs,0, 1))If you have continuous-scale time-series data, start with this step to convert your data
# transpose to wide-format
continuous.time.series <- t(continuous.time.series)
# check to see the continuous-scale
print(continuous.time.series[,1:5]) # print the first 5 columns in the binarized data## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.7209039 0.8757732 0.7609823 0.8861246 0.4564810
## [2,] 0.4537281 0.3267524 0.9654153 0.7074819 0.6445426
## [3,] 0.7821933 0.4291988 0.9272740 0.7732432 0.2596812# binarize the continuous-scale variables
bin.continuous.time.series <- binarizeTimeSeries(continuous.time.series)
# show the threshold of the 3 variables
bin.continuous.time.series$thresholds## [1] 0.5600568 0.5183026 0.5204997# print the first five time points
print(bin.continuous.time.series$binarizedMeasurements[,1:5]) # print the first 5 columns in the binarized data## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 1 1 1 0
## [2,] 0 0 1 1 1
## [3,] 1 0 1 1 0Here, we can see the function “binarizeTimeSeries” split the continuous-scale data into two groups using the threshold, e.g., any value that is above 0.5600568 in the first variable (first row) will be converted to 1, otherwise it will be converted to 0.
Step 1: Simulation of binary time-series data
Why simulate?
This step is to generate a binary-scale multivariate time-series which allow us to look at how the model works without empirical data. It has the advantage that we know what is the underlying temporal dynamics with the simulated data, and examine how the method recovered the temporal dynamics in estimates.
How is the data simulated?
Data are simulated for a 3-node network. Time series was simulated based on the temporal relations and process noise. Our hypohetical 3-node newtork was using simulated data based on a pre-defined Boolean functions and process noise. The Boolean functions are expressed in the following equations:
x1(t+1) = x3(t)
x2(t+1) = x1(t) OR x3(t)
x3(t+1) = NOT x1(t) AND x2(t)
The process noise is to flip the binary variable at each time step with a probability of 0.05.
The three-node network is a hypothetical example used for the purpose of illustration. The Boolean functions and noise probability were arbitrarily chosen.
We simulate 20 steps to generate a relatively short time-series.
p <- 3 # number of variables
# randomly generated number for each time-step, so that the noise can be added
zeta <- cbind(runif(n.obs,0, 1),
runif(n.obs,0, 1),
runif(n.obs,0, 1)) # 3 time-series of noise for 3 variables
zeta <- t(zeta) # adjust zeta to a wide-format
time.series <- matrix(rep(0, p * n.obs), nrow =p , ncol = n.obs)
time.series[,1] <- runif(p,0, 1) # initial value of the three variables were randomly chosen
chance.to.flip <- 0.05 # probability of adding noise at each time step
time.series[,1] <-ifelse(time.series[,1] >= .5, TRUE,FALSE)
for (col in 2:n.obs)
{
# simulate each time step based on the previous time step and Boolean functions
time.series[1,col] <- time.series[3,col-1] # x1 = x3
time.series[2,col] <- time.series[1,col-1] | time.series[3,col-1] # x2 = x1 OR x3
time.series[3,col] <- !time.series[1,col-1] & time.series[2, col -1] # x3 = NOT x1 AND x2
# add noise, if the random number in zeta is less than 0.05, then taking the opposite value
time.series[1,col] <- ifelse(zeta[1,col] < chance.to.flip,
!time.series[1,col],
time.series[1,col] )
time.series[2,col] <- ifelse(zeta[2,col] < chance.to.flip,
!time.series[2,col],
time.series[2,col] )
time.series[3,col] <- ifelse(zeta[3,col] < chance.to.flip,
!time.series[3,col],
time.series[3,col] )
}
# adjust the time-series to long-format
time.series <- t(time.series)
time.series <- data.frame(time.series)
names(time.series) <- c("x1", "x2", "x3")Plot the binary time-series
Here we plot the converted time-series, and now it is in the binary format. Any colored area indicates the variable is 1, and white area indicates the variable is 0.
# specify the colors of the time-series
n <- ncol(time.series)
cols <- gg_color_hue(n)
# for plotting purposes, convert the binary time-series back to long-format
bin.data.plot <- time.series
bin.data.plot <- data.frame(bin.data.plot)
bin.data.plot$index <- 1:nrow(bin.data.plot)
bin.data.plot.melt <- melt(bin.data.plot,id = "index")
bin.data.plot.melt$value.variable <- ifelse(bin.data.plot.melt$value==0,"0", as.character(bin.data.plot.melt$variable))
bin.data.plot.melt$value.variable <- factor(bin.data.plot.melt$value.variable,
levels = c("0","x1",
"x2",
"x3"))
ggplot(data = bin.data.plot.melt)+
geom_rect(aes(xmin = index - .5, xmax = index + .5,
ymin = 0, ymax = 1, fill = factor( value.variable)))+
facet_wrap(~variable, ncol = 1) +
scale_fill_manual(values = c("#FFFFFF", cols))+
theme(
strip.background = element_blank(),
panel.background = element_blank(),
legend.title = element_blank(),
legend.key = element_blank(),
legend.position = "none",
axis.text.y=element_text(color="black",size=10),
axis.text.x=element_text(color="black",size=10),
axis.title.y=element_text(color="black",size=10),
axis.title.x=element_text(color="black",size=10),
axis.line = element_line(color = 'black'))+
ylim(0,1)+
xlab("Time") 
Step 2: Inference of Boolean functions
net.data <- t(time.series)
network.size<- p
booleannet <- reconstructNetwork(net.data,
method = "bestfit",
maxK = 2,
readableFunctions=T,
returnPBN = F)
print(booleannet)## Probabilistic Boolean network with 3 genes
##
## Involved genes:
## x1 x2 x3
##
## Transition functions:
##
## Alternative transition functions for gene x1:
## x1 = (x3) (error: 1)
##
## Alternative transition functions for gene x2:
## x2 = (x3) | (x1) (error: 1)
##
## Alternative transition functions for gene x3:
## x3 = (!x1 & x2) (error: 1)# saveNetwork(net, paste(outputpath, id,"-network.txt", sep = "")) # can't save a collection of BNWe can see the Boolean functions we simulated the data from are recovered here by the inference, which are:
x1(t+1) = x3(t)
x2(t+1) = x1(t) OR x3(t)
x3(t+1) = NOT x1(t) AND x2(t)
Note that this is only to showcase how the Boolean network inference works, future work can execute a comprehensive simulation study, which requires robust testing of multiple factors, including a variety of forms of Boolean funtions, observation length, noise.
How Boolean functions are inferred?
The details of how the Boolean functions are inferred can be found in Akutsu et al. (2000). The essence of the inference is to compare the outcome varible’s time-series at time t+1 and the time-series of the input variables at time t, after appyling the Boolean operator AND, OR, or NOT; and find which predictors and Boolean operators are the closest match of the outcome variable.
What is error?
Error here is the sum of false negative (predicting the outcome as 0 when it was actually 1) and false positive (predicting the outcome as 1 when it was actually 0). Note the total number of time points is 20 in the simulated time-series, and error = 1 indicates out of 20 time points 1 was predicted incorrectly using the “bestfit” algorithm.
How to interpret the Boolean functions?
Take function x2 = (x1 | x3) as an example:
This indicates x2(t+1) is predicted by x1(t) OR x3(t), which means only when either x1(t) or x3(t) are 1, x2(t+1) is 1; otherwise, x2(t+1) will be 0. OR rule is similar to an additive rule, so whenever there is a 1 in an OR rule, the outcome will be 1.
Another way to talk about the interpretation of this Boolean function is x2 will be turned ON at time t+1 when either x1 or x3 was turned ON at time t.
What are the possible applications in psychology?
In psychology, when two variables have to work together to activate another variable, it is an AND rule, e.g., when both group member A and group member B are self-disclosing, group member B will continue to self-disclose the next time (maybe member B needed to feel supported from member A). A Boolean function that looks like B(t+1) = A(t) & B(t) describes such a dynamic.
Here, we select the fisrt network if more than 1 rule inferred per node, so that the attractors can be extracted.
# note here we only chose the first function of all possible functions
singlenet <- chooseNetwork(booleannet,
functionIndices = rep(1,network.size),
dontCareValues=rep(1,network.size),
readableFunctions=T)
print(singlenet)## Boolean network with 3 genes
##
## Involved genes:
## x1 x2 x3
##
## Transition functions:
## x1 = (x3)
## x2 = (x3) | (x1)
## x3 = (!x1 & x2)Why extract attractors first?
Note here, we extract attractors first (Step 4) because state transition graph needs the attractor object as an input. But we will explain what the attractors mean later, as the state transition graph (Step 3) will help you understand the concept of attractor.
ga <- getAttractors(network = singlenet,
returnTable = TRUE)
# print(ga, activeOnly=TRUE)Step 3: Plot state space transition graph
p<-plotStateGraph(ga,
piecewise=TRUE,
drawLabels = T,
plotIt = F,
colorsAlpha = c(colorBasinsNodeAlpha = 1,
colorBasinsEdgeAlpha = 1,
colorAttractorNodeAlpha = 1,
colorAttractorEdgeAlpha = 1))
plot.igraph(p,
label.cex = 1.2,
vertex.label.color="black",
vertex.label.dist=3,
remove.loops = T,
edge.arrow.size=.4)
What does state transition graph mean?
Plot the state transition graph, where each node is a state of the system (a vector of state of x, y, z), and the arrows/edges are state-to-state transition. E.g., state 100 means (x = 1, y = 0, z = 0); and the state 100 will transition to state 000 in attractor 1. The state transition is computed based on the Boolean functions inferred in the first step. e.g., because the x(t+1) = x(t) & z(t), so when x(t) = 1 and z(t) = 0, x(t+1) = 1 & 0 = 0, Similarly, y(t+1) = y(t) = 0, z(t+1) = z(t) = 0, so the next state of the system is (0,0,0), hence state 100 transitions to state 000.
What is an attractor and how can I tell from the state transition graph?
When a state transitions back to itself, it is an attractor - as the system gravitates toward an attractor. Here in each attractor, there is 1 or multiple attractor state(s) that the system tends to stay for a long run.
What are the possible applications in psychology?
In psychology, using the same example of group-level self-disclosure dynamic, both member A and member B can be stuck in an attractor where neither of them self-discloses - A = 0, and B = 0 - depending on the Boolean functions. So discovering where the attractors are will lead to understanding of the problems in the system; and this also alludes to the network control part in Step 6 because naturally we would want to intervene when problematic or undesirable attractors are found.
Step 4: Extraction of attractors
print(ga, activeOnly=TRUE)## Attractor 1 is a simple attractor consisting of 1 state(s) and has a basin of 1 state(s).
## Active genes in the attractor state(s):
## State 1: --
##
## Attractor 2 is a simple attractor consisting of 3 state(s) and has a basin of 7 state(s).
## Active genes in the attractor state(s):
## State 1: x2
## State 2: x3
## State 3: x1, x2What do we know about the attractors?
We have 2 attractor here. The first one is a fixed point attractor, where all three nodes are OFF. The second one is a complex attractor, where the system transition among 3 states, which are x2 is ON corresponding system state (0,1,0), x3 is ON corresponding system state (0,0,1), and both x1 and x2 are ON corresponding state (1,1,0).
Step 6: Network control
Why network control is useful?
After extracting the attractors, we can define its desirability based on practical concerns. Using the parent-child dynamic example, we might want to direct the parent and/or child out of a yelling attractor. That is why network control is useful.
How to design network control?
Network control can be designed based on a search algorithm to find the node to flip so that the system can move from the undesirable attractors to the desirable attractor, as an intervention. For example, if we consider 000 is a desirable attractor, and 101 is an undesirable attractor, we can flip the third node, so that 101 will become 100, and then go into 000 in 1 step.
What are the possible applications in psychology?
For the example of group-level self-disclosure, there could be a node flip to direct the system out of an undesirable state such as neither member A or B is self-disclosing to a desirable state such as both are self-disclosing. A node flip is when a member’s behavior changes from ON (disclose) to OFF (not disclose) or from OFF to ON, and this can serve as an intervention strategy to move the group into a much more desirable attractor, given the right conditions.
A function to find the element index of the minimal value in a vector
which.is.min <- function (x)
{
y <- seq_along(x)[x == min(x)]
return(y)
}Search algorithm for control strategy - flipping one node
The search algorithm here is to compute Hamming distance between two states, where the first state is a state in an undesirable attractor, and the second state is a state in a desirable attractor basin.
Then the shortest Hamming distance can help identify how to perturb specific nodes to move the system from an undesirable attractor to a desirable attractor.
The specific control strategy as such: Suppose node x1 is the node that differs in the two states in the undesirable attractor and the desirable attractor basin (denoted su, sd respectively) that have the shortest Hamming distance, then the node to perturb is x1, the state to perturb to is the corresponding state of x1 in sd.
There are nested loops in the following code, here are a couple reasons why nested loops are used:
- in two loops we searched for all the states in undesirable attractor and all the states in a desirable attractor basin
- It becomes complicated to compute the Hamming distance when the undesirable attractor is a complex attractor (or limit cycle).
control.solution <- NULL
ga.undesirable <- NULL # undesirable attractor list
ga.desirable <- NULL # desirable attractor listManually assign the desirablity to (the index of) attractors
note here the index used is according to the index in the code as follows (executed previously):
ga <- getAttractors(network = singlenet,returnTable = TRUE)
You can check the index of attractors by running the following function
print(ga, activeOnly=TRUE)
ga.undesirable <- 1 # hypothetically desirable
ga.desirable <- 2 # hypothetically undesirable# state space transition table - sstt
sstt <- getTransitionTable(ga)
# att.no <- unique(sstt$attractorAssignment)
# node.list <- singlenet$genes
# the states where transition starts
transition.from <- NULL
for (col in 1:network.size)
{
transition.from <- paste(transition.from, sstt[,col], sep = "")
}
# the states where transition ends
transition.to <- NULL
for (col in (network.size + 1):(2 * network.size))
{
transition.to <- paste(transition.to, sstt[,col], sep = "")
}
# transition.table will be used to compute distance from undesirable states to desirable states
transition.table <- data.frame(cbind(transition.from, transition.to))
names(transition.table) <-c("from","to")
transition.table$attractor.assignment <- ga$stateInfo$attractorAssignment
transition.table$stepstoa <-ga$stateInfo$stepsToAttractor
if (length(ga.desirable) > 0 &length(ga.undesirable)>0 )
{
for (ua in 1:length(ga.undesirable))
{
for (da in 1:length(ga.desirable))
{
undesirable.states <-transition.table[transition.table$attractor.assignment ==
ga.undesirable[ua] & transition.table$stepstoa ==0 , ]$from
# "from" are the states in the desirable attractor basin
desirable.states<- transition.table[transition.table$attractor.assignment ==
ga.desirable[da], ]$from
# "to" is the state of the desirable attractor
desirable.attractor <- transition.table[transition.table$attractor.assignment ==
ga.desirable[da]& transition.table$stepstoa ==0 , ]$to
# distance from undesirable states to desirable states
distance.mat <- stringdistmatrix(undesirable.states, desirable.states, method = "hamming")
# identify the node to flip, if there is a simple control strategy by flipping one node
if(min(distance.mat)==1)
{
# multiple undesirable states (complex attractor)
for (row in 1:nrow(distance.mat))
{
# find the states in attractor basin that has distance of 1
state.index <- which.is.min(distance.mat[row,])
# multiple minimal distance states
for(index in 1:length(state.index))
{
# column index of the minimal distance state in the desirable attractor basin
col.index <- state.index[index]
# loop of search which character has the distance of 1, that will be the node to flip
for (ichar in 1:nchar(as.character(desirable.states[col.index])))
{
# find the character that differs between desirable state and undesirable state/attractor
if (substr(as.character(desirable.states[col.index]), ichar, ichar) !=
substr(as.character(undesirable.states[row]), ichar, ichar))
{
flip.node.index <- ichar
}# end of the if condition to find character that differs
}# end of the loop of search which character has the distance of 1, that will be the node to flip
control.solution <- rbind(control.solution,
c(as.character(undesirable.states[row]), # undesirable attractor
as.character(desirable.states[col.index]), # a state in desirable attractor basin
flip.node.index, # node to flip
as.character(desirable.attractor) # desirable attractor
))
} # end of the loop of multiple minimal distance states
}# end of loop of multiple undesirable states (complex attractor)
} # end of if condition of a simple control strategy by flipping one node
else{
flip.node.index <- "no.bitflip.found"
# multiple undesirable states (complex attractor)
for (row in 1:nrow(distance.mat))
{
# record each undesirable state, which node can you flip
control.solution <- rbind(control.solution,
c(as.character(undesirable.states[row]), # undesirable attractor
as.character(desirable.states[col.index]), # a state in desirable attractor basin
flip.node.index, # node to flip
as.character(desirable.attractor) # desirable attractor
))
}# end of loop of multiple undesirable attracotrs
}# end of else condition there is a bitflip
}# end of da loop
}# end of ua loop
}# end of if condition at least 1 desirable and undesirable attractor
control.solution <- data.frame(control.solution)
names(control.solution) <- c("undesirable.attractor", "state.in.desirable.basin", "node.to.flip","desirable.attractor")Print out the control strategy
print(control.solution)## undesirable.attractor state.in.desirable.basin node.to.flip
## 1 000 100 1
## 2 000 010 2
## 3 000 001 3
## desirable.attractor NA NA
## 1 001 010 110
## 2 001 010 110
## 3 001 010 110Now you have the control strategy identified! There are 3 strategies you can use to move the system from the undesirable attractor (000) to the desirable attractor (a limit cycle attractor with 3 states 001, 010, 110).
The first control strategy is to turn node 1 ON (in the node to flip column) when the system is in the undesirable state 000, then the system will move into the state 100, which is in the desirable attractor basin, and the system will eventually move and stay in the desirable attractor once it is in the basin of the desirable attractor (see the previous state transition graph for a visual represetation of the desirable attractor basin).
Similarly, two other control strategies are shown in row 2 and 3.
Summary
Now you have finished following the tutorial, and let me reiterate the key steps to model and control Boolean network:
- Inference of the Boolean functions
- Extraction of attractors and assign desirability to the attractors
- Design network control
If you have questions or suggestions, feel free to drop me an email xfy5031[at]psu.edu! I will be very excited to hear what you plan or have done with your own data using the Boolean network method!
References
Akutsu, T., Miyano, S., Kuhara, S. (2000). Algorithms for identifying Boolean networks and related biological networks based on matrix multiplication and fingerprint function. Journal of Computational Biology, 7(3), 331-343.
Boker, S., & Graham, J. (1998). A dynamical systems analysis of adolescent substance use. Multivariate Behavioral Research, 33(4), 479-507.
Chow, S., Ferrer, E., & Hsieh, F. (2010). Statistical Methods for Modeling Human Dynamics: An Interdisciplinary Dialogue (Notre Dame Series on Quantitative Methodology, Vol 4). New York, NY: Taylor & Francis.
Kauffman, S. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology, 22, 437-467.
Kauffman, S. (1993). The Origins of Order: Self-organization and Selection in Evolution. Oxford, UK: Oxford University Press.
Lähdesmäki, H., Shmulevich, I., & Yli-Harja, O. (2003). On learning gene regulatory networks under the Boolean network model. Machine Learning, 52(1), 147-167.
Molenaar, P.C.M., & Nesselroade, J.R. (2015). Systems methods for developmental research. In R.M. Lerner (Ed) Handbook of Child Psychology and Developmental Science (pp 1-31). New York, NY: Wiley.
Müssel, C. Hopfensitz, M., & Kestler, H. (2010). BoolNet – An R package for generation, reconstruction, and analysis of Boolean networks. Bioinformatics, 26(10), 1378-1380.
Wang, Q., Molenaar, P.C.M., Harsh, S., Freeman, K., Xie, J., Gold, C., Rovine, M., & Ulbrecht, J. (2014). Personalized state-space modeling of glucose dynamics for Type I diabetes using continuously monitored glucose, insulin dose, and meal intake: An extended Kalman filleter approach. Journal of Diabetes Science and Technology, 8(2), 331-345.