#Supplementary information

#R Code used for the simulation of the Boolean model:

#Dependencies
library(BoolNet)

#Boolean rules exportation from text file
boolnetwork <- loadNetwork('file.net')

#Search of synchronous and asynchronous attractors
synch_attractors <- getAttractors(boolnetwork, 'synchronous', method = "exhaustive")
asynch_attractors <- getAttractors(boolnetwork, 'asynchronous', method = "exhaustive")


########################################################################################3
#R Code used for the simulation of the network as a continuous model:

#Dependencies
library(BoolNet)
library(deSolve)
library(BoolNetPerturb)

#Generation and simulation of the continuous model 
cont_simulation <- function(boolnetwork, LOGIC="Zadeh", EQ="Villarreal", indx_nodes_1=c(1,2,19,20,23,26), b=10, timeSteps=20){
  ODEnetwork <- booleanToODE(boolnetwork, logic=LOGIC, eq = EQ)
  ODEnetwork$parameters[1] <- b
  ODEnetwork$state[1:length(boolnetwork$genes)] <- 0
  ODEnetwork$state[indx_nodes_1] <- 1
  sim = ode(func = ODEnetwork$func, 
          parms = ODEnetwork$parameters, 
          y = ODEnetwork$state,
          times = seq(0, timeSteps, 0.2))
  return(sim)
}
o=cont_simulation(network)

#Mutant simulation, example for NFkB overactivation
network_mut <-fixGenes(network, 13, 1)
o=cont_simulation(network_TF_diff_mut)


########################################################################################3
#R Code used for the simulation of the continuous model and the Jaccard index calculation:

#Dependencies
library(BoolNet)
library(deSolve)
library(BoolNetPerturb)
library(philentropy)
library(ggplot2)

#Jaccard index calculation between two vectors
Jaccard = function (x, y) {
  M.11 = sum(x == y)
  M.10 = length(x)
  M.01 = length(y)
  return (M.11 / (M.10 + M.01 - M.11))
}

#Jaccard index calculation between two simulations with n vectors, where n is equal to the number of the simulation time-steps
jaccardTS <- function(WTsim, MUTsim, ind_geneMut=0){
  Jacc_TS <- c()
  for(r in 2:nrow(WTsim)){
    xy <- rbind(WTsim[r,2:ncol(WTsim)], MUTsim[r,2:ncol(MUTsim)])
    xy <-round(xy,1)
    if(ind_geneMut>0){
      #print("Deleting mutated gene")
      xy <- xy[,-ind_geneMut]
    }
    #print(Jaccard(xy[1,],xy[2,]))
    Jacc_TS <- c(Jacc_TS, Jaccard(xy[1,],xy[2,]))
  }
  return(Jacc_TS)
}

#Distance to the wild type simulation for all the mutant continuous simulations
muts_Jaccard <-function(network, TimeSteps=20){
  wt=round(cont_simulation(network, timeSteps = TimeSteps),1)
  result_table <- as.data.frame(matrix(0, ncol=nrow(wt)-1, nrow=((length(network$genes))*2)+1))
  rownames(result_table) <- c("wt", paste0(network$genes, "_0"), paste(network$genes, "_1"))
  count=1
  result_table[count,] <- jaccardTS(wt, wt)
  for(b in 0:1){
    for(a in 1:length(network$genes)){
      mutNet <- fixGenes(network, a, b)
      mut <- round(cont_simulation(mutNet, timeSteps = TimeSteps),1)
      count=count+1
      print(rownames(result_table)[count])
      result_table[count,] <- jaccardTS(wt, mut, a)
    }
  }
  return(result_table)
}

#Figure 3 distance heatmap generation
muts_Jaccarad_result <- muts_Jaccard(network)
result <- as.matrix(muts_Jaccarad_result)
result <- result[-1,]
my_palette <- colorRampPalette(c("firebrick3", "lightgoldenrodyellow", "dodgerblue4"))(n = 299)
heatmap(mm, Colv = NA, scale="none", labCol = c(1:ncol(result)), col=my_palette)