In this notebook is reported the code for evaluation of the predictive power of FLT3-ITD logic models.

In each cell line, we simulated:

We used the Boolean model with the lowest error between experimental and simulated data in the two cell lines (best model).

Computation of steady-state in tumor and FLT3 inhibition condition

for(cellline in c('JMD', 'TKD')){
 
  for (condition in c('tumor', 'FLT3i')){
    
    # Load the subfamily of 100 models
    load(paste0('./results/', cellline, '_subfamily_threshold.RData'))
    
    # If the condition is tumor set IGFR, TNF and FLT3 ON
    if(condition == 'tumor'){
      ## -- TUMOR SIMULATION -- ##
      
      # Restrict cues to the nodes of the simulation
      CNOlist@cues[1,1:3] <- rep(1, 3) 
      CNOlist@cues <- t(as.matrix(CNOlist@cues[1,1:3])) 
      
      # Remove the inhibitors of FLT3 to set the ON
      CNOlist@inhibitors[1,1] <- 0 
      CNOlist@inhibitors <- t(as.matrix(CNOlist@inhibitors[1,1]))
      
      # Add the stimuli of FLT3, TNF and IGF1R to set the ON
      v <- c(1,1,1)
      names(v) <- c('FLT3', 'TNF', 'IGF1R')
      CNOlist@stimuli <- t(as.matrix(v))
      
    }else{ 
      ## -- FLT3 INHIBITION -- ##
      
      # Restrict cues to the nodes of the simulation
      CNOlist@cues[1,1:3] <- rep(1, 3)
      CNOlist@cues <- t(as.matrix(CNOlist@cues[1,1:3]))
      
      # Put the inhibitor of FLT3 to set it OFF
      CNOlist@inhibitors 
      CNOlist@inhibitors[1,1] <- 1
      CNOlist@inhibitors <- t(as.matrix(CNOlist@inhibitors[1,1]))
      
      # Put the stimuli of IGFR1 and TNF to set them ON
      CNOlist@stimuli[1,1:2] <- 1
      CNOlist@stimuli <-  t(as.matrix(CNOlist@stimuli[1,1:2]))
    }
    
    ### -- MANUAL SIMULATION -- ##
    # Creating object for simulation
    simList <- prep4sim(modelT0) 
    
    matrix_sim_q <- list()
    for(i in c(1)){ #select only the best model
      
      bString <- as.vector(unlist(best_100_models[i,]))
      #print(bString)
      # cutting expanded gates to optimal bString
      modelCut <- cutModel(modelT0, list(bString)[[1]])
      
      # creating a new simList from cutted model
      newSimList <- cutSimList(simList, list(bString)[[1]])
      
      # creating index list
      indexList <- indexFinder(CNOlist, modelT0)
      
      ## ** SIMULATION A T0 (before treatment) ** 
      simResT0 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList,
                              indexList=indexList, mode = 0)
      colnames(simResT0) <- modelT0$namesSpecies
      
      ## ** SIMULATION A T1 (after treatment) **
      simResT1 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList, 
                              indexList=indexList, mode = 1)
      colnames(simResT1) <- modelT0$namesSpecies
      
      ### -- STORING RESULTS in data objects -- ###
      
      # list object
      simRes <- list(T0 = simResT0, T1 = simResT1)
      
      if(sum(is.na(simRes$T0)) > 0 | sum(is.na(simRes$T1)) > 0){
        matrix_sim_q[[i]] <- list()
      }else{
        matrix_sim_q[[i]] <- list()
        matrix_sim_q[[i]] <- simRes
      }
    }
    
    matrix_sim_t0 <- lapply(matrix_sim_q, function(x){x$T0})
    matrix_sim_t0 <- Filter(Negate(is.null), matrix_sim_t0)
    averaget0 <-  Reduce("+", matrix_sim_t0) / length(matrix_sim_t0)
    
    matrix_sim_t1 <- lapply(matrix_sim_q, function(x){x$T1})
    matrix_sim_t1 <- Filter(Negate(is.null), matrix_sim_t1)
    averaget1 <-  Reduce("+", matrix_sim_t1) / length(matrix_sim_t1)
    
    # matrix object
    sim_matrix <- matrix(nrow = length(simRes$T0), ncol = 2)
    
    sim_matrix[,1] <- t(as.matrix(averaget0)) 
    sim_matrix[,2] <- t(as.matrix(averaget1))
    rownames(sim_matrix) <- colnames(simResT1)
    colnames(sim_matrix) <- c('T0', 'T1')
    
    
    if(condition == 'tumor'){
      sim_matrix_new <- sim_matrix
    }else{
      sim_matrix_new <- cbind(sim_matrix_new, sim_matrix[,'T1'])
      colnames(sim_matrix_new) <- c('T0', 'T10', 'T90')
    }
  }
  
  sim_matrix_new <- sim_matrix_new[,c('T10', 'T90')] 
  colnames(sim_matrix_new) <- c('tumor', 'FLT3i')
  ### -- PLOTTING RESULTS AS A HEATMAP -- ###
  paletteLength <-100
  myColor <- colorRampPalette(c('darkblue', "red3"))(paletteLength)
  # simulation's result as heatmap
  hp <- pheatmap(t(sim_matrix_new),
                 color = myColor,
                 #display_numbers = TRUE)
                 cellwidth = 20,
                 cellheight = 20,
                 cluster_rows = FALSE,
                 cluster_cols = FALSE,
                 border_color = 'white',
                 legend_breaks = c(0,0.5, 1),
                 fontsize = 12,
                 width = 5,
                 heigth = 30,
                 #labels_col = c('T10', 'T90'),
                 main = paste0('FLT3 ITD-', cellline),
                 angle_col = c('90'),
                 silent = TRUE)
  
  if(cellline == 'TKD'){
    TKD_hp <- hp
    sim_matrix_TKD <- sim_matrix_new
  }else{
    JMD_hp <- hp
    sim_matrix_JMD <- sim_matrix_new
  }
}

# Write steady-state results
sim_tibble_JMD <- tibble(rownames_to_column(data.frame(sim_matrix_JMD), 'node'))
#write_tsv(sim_tibble_JMD, './results/JMD_quizartinib_simulation_steady_state.tsv')

sim_tibble_TKD <- tibble(rownames_to_column(data.frame(sim_matrix_TKD), 'node'))
#write_tsv(sim_tibble_TKD, './results/TKD_quizartinib_simulation_steady_state.tsv')

The results of the quizartinib simulation are shown below:

grid.arrange(grobs = list(JMD_hp[[4]], TKD_hp[[4]]), ncol = 1, nrow = 2)

Inference of Apoptosis and Proliferation

To functionally interpret the results of the simulations, for each network, we first extracted key regulators of ‘apoptosis’ and ‘proliferation’ hallmarks from SIGNOR. To this aim, we applied our recently-developed ProxPath algorithm, a graph-based method able to retrieve significant paths linking nodes of our two optimized models to proliferation and apoptosis phenotypes.

Read table with regulators of Proliferation and Apoptosis of ProxPath

## # A tibble: 6 × 3
##   node   phenotype       apoptosis
##   <chr>  <chr>               <dbl>
## 1 STAT5A APOPTOSIS              -1
## 2 MTOR   APOPTOSIS              -1
## 3 JNK    APOPTOSIS               1
## 4 AKT    APOPTOSIS              -1
## 5 AKT    APOPTOSIS              -1
## 6 ERK1_2 DIFFERENTIATION         1

Compute Apoptosis and Proliferation state

phens <- c('APOPTOSIS', 'PROLIFERATION')
plot_list <- list()
for(i_phen in c(1:length(phens))){
  i_phen
  phen <- phens[i_phen]
  
  # Select the regulators of the i- phenotype
  apoptosis_annotation <- apoptosis_df %>% 
    filter(phenotype == phen) %>% 
    select(node, apoptosis)  %>%
    distinct()
  
  # Annotate with protein regulatory role on phenotype each model protein
  annotated_JMD <- left_join(sim_tibble_JMD, apoptosis_annotation %>% 
                               select(node, apoptosis), by = 'node') %>%
    filter(!is.na(apoptosis))
  
  annotated_TKD <- left_join(sim_tibble_TKD, apoptosis_annotation %>% 
                               select(node, apoptosis), by = 'node') %>%
    filter(!is.na(apoptosis)) 
  
  # Select only the regulators in each cell line that are endpoint proteins
  # in high-confidence signaling axes (edge frequency 0.4). As such, 
  # if two regulators of the same phenotype were linked in the same axis, 
  # we considered only the one at the end of the cascade.
  if(phen == 'APOPTOSIS'){
    annotated_JMD <- annotated_JMD %>% 
      filter(node %in% c('MTOR', 'STAT5A'))
    
    annotated_TKD <- annotated_TKD %>%
      filter(node %in% c('STAT5A', 'AKT'))
  }else{
    annotated_JMD <- annotated_JMD %>% 
      filter(node %in% c('PI3K', 'MAPK14', 'STAT5A', 'MTOR'))
    
    annotated_TKD <- annotated_TKD %>%
      filter(node %in% c('PI3K', 'KRAS', 'CREB1', 'STAT5A'))
  }
  
  for(cellline in c('JMD', 'TKD')){
    #cellline = 'JMD'
    
    if(cellline == 'JMD'){
      annotated <- annotated_JMD
    }else{
      annotated <- annotated_TKD
    }
    
    # Multiply the role of proteins (apoptosis) and their activation status (T10/T90)
    # and integrate them with an OR logic (equal contribution with the mean)
    apo <- annotated %>% 
      filter(!is.na(apoptosis)) %>%
      mutate(tumor_apo = tumor*apoptosis,
             FLT3i_apo = FLT3i*apoptosis) %>% 
      group_by(apoptosis) %>% 
      summarise(tumor = sum(tumor_apo),
                FLT3i = sum(FLT3i_apo))
    
    apo <- apo %>% summarise(
      tumor = sum(tumor),
      FLT3i = sum(FLT3i))
    
    if(cellline == 'JMD'){
      apo$cellline <- 'JMD'
      apo_JMD <- apo
    }else{
      apo$cellline <- 'TKD'
      apo_TKD <- apo
    }
  }
  
  apo_tb <- bind_rows(apo_JMD, apo_TKD)
  apo_matrix <- as.matrix(apo_tb %>% column_to_rownames('cellline'))
  
  a <- data.frame(t(apo_tb)) %>% rownames_to_column()
  JMD <- a[-3, c('rowname', 'X1')]
  colnames(JMD) <- c('time', 'apoptosis')
  JMD$cellline <- 'JMD'
  
  TKD <- a[-3, c('rowname', 'X2')]
  colnames(TKD) <- c('time', 'apoptosis')
  TKD$cellline <- 'TKD'
  
  barplot_tb <- bind_rows(JMD, TKD) %>% 
    mutate_at('apoptosis', as.numeric)
  
  # To make the plot more readble set opposite the Apoptosis
  if(phen == 'APOPTOSIS'){
    barplot_tb$apoptosis <- barplot_tb$apoptosis * -1
  }
  
  # DRAW A BARPLOT OF EACH PHENOTYPE IN EACH CELL LINE
  barplot_apoptosis <- ggplot(barplot_tb, aes(x = time, y = apoptosis, fill = cellline)) +
    geom_bar(stat = 'identity', position = 'dodge') +
    #labs(title = str_replace(conditions_sim, '\\.', '+'))+
    xlab('') + ylab(phen)+
    scale_fill_manual(values = c('#36437F','#E1BD6A'))+
    theme_classic() +
    ylim(0,6)+
    scale_x_discrete(labels = c('tumor' = 'TNF + IGF \n + FLT3-ITD',
                                'FLT3i' = 'FLT3i'))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 8),
          axis.text.y = element_text(size = 8), 
          axis.title = element_text(size = 8), 
          legend.text = element_text(size=8))
  
  plot_list[[i_phen]] <- barplot_apoptosis
  
  
  # DRAW AN HEATMAP OF THE STATUS OF PHENOTYPE REGULATORS 
  colnames(annotated_JMD) <- c('node', paste0('JMD', c('', '_FLT3i')), 'apoptosis')
  colnames(annotated_TKD) <- c('node', paste0('TKD', c('', '_FLT3i')), 'apoptosis')
  
  cell_simulation_tb <- full_join(annotated_JMD, annotated_TKD)
  
  apoptosis_phen <- cell_simulation_tb %>% column_to_rownames('node')
  
 
  apoptosis_phen <- apoptosis_phen*apoptosis_phen$apoptosis
  apoptosis_phen$apoptosis <- NULL
  
  if(phen == 'APOPTOSIS'){
    apoptosis_phen = apoptosis_phen * -1
  }
  
  myColor <- colorRampPalette(c("darkblue", "white", "red3"))(3)
  
  hp <- pheatmap(t(apoptosis_phen),
                 color = myColor,
                 #display_numbers = TRUE)
                 cellwidth = 10,
                 cellheight = 10,
                 cluster_rows = FALSE,
                 cluster_cols = FALSE,
                 border_color = 'grey',
                 #breaks = c(-1, 0, 1),
                 fontsize = 8,
                 width = 5,
                 heigth = 30,
                 #labels_col = c('T10', 'T90'),
                 #main = paste0('quizartinib simulation in ', cellline),
                 angle_col = c('45'),
                 silent = TRUE)

  plot_list[[i_phen+2]] <- hp[[4]]
}

n <- length(plot_list)
nCol <- floor(sqrt(n))
p <- do.call("grid.arrange", c(plot_list, ncol=2))

Figure 3E. Heatmaps (bottom) report the activation level of positive and negative phenotype regulators present in the two Boolean models. Bar plots (upper) showing the proliferation activation and apoptosis inhibition levels in untreated and FLT3i conditions in the steady states of FLT3 ITD-JMD (yellow) and FLT3 ITD-TKD (blue) Boolean models.