Simulation of patients

We tested whether we could exploit the cell-derived Boolean models to identify novel personalized combinatorial treatments. To these aims, each patient’s mutational profile (Fig. S6A and D, Table S5) was turned into the initial activity of a personalized Boolean model (Fig. 5A, panel c).

Next, for each patient, we performed a simulation of the following conditions in silico:

Importantly, our approach enabled us to obtain patient-specific predictive Boolean models able to describe the drug-induced signaling rewiring (Fig. 5F and Fig. S7) and to quantify apoptosis inhibition and proliferation levels (Fig. 5A, panel d and e, Fig. 5E and Fig. S6E).

# Read the mutation file of the patients
mutations <- read_tsv('./input/FLT3_ITD_patients_mutational_profile.tsv')
mutations[mutations == 0] <- 3 #Set NA mutation to a non Boolean value
mutations[mutations == 1] <- 0 #Set LOF mutation to Boolean value 0
mutations[mutations == 2] <- 1 #Set GOF mutation to Boolean value 1

# Create an empty list for each patient' plots
patients_plot_list <- list()

# Counter for plot list

plot_idx = 1
# For each patient
for(p_idx in 1:nrow(mutations)){
  
  # Patient name and mutation vector
  patient <- mutations$Patient[p_idx]
  mutations_vector <- unlist(mutations[p_idx,2:6])
  
  # Get the two cell-derived models subfamily
  JMD_session <- load(paste0('./results/JMD_subfamily_threshold.RData'))
  JMD_best <- as.vector(unlist(best_100_models[1,]))
  JMD_model <- modelT0
  JMDsimList <- simList
  
  TKD_session <- load(paste0('./results/TKD_subfamily_threshold.RData'))
  TKD_best <- as.vector(unlist(best_100_models[1,]))
  TKD_model <- modelT0
  TKDsimList <- simList
  
  # set the conditions we want to simulate 
  # 1) patients' mutational profile + active TNF IGF1 receptors
  # 2) condition 1 + FLT3i
  # 3) condition1 + FLT3i + druggable node
  
  # To create the collection of nodes to inhibit
  # we start from the node inhibited in the perturbation experiment
  inhibited <- colnames(CNOlist@inhibitors)
  
  # We expand the 'inhibited' condition with
  # the 'mutations' and four additive kinases 'KRAS', 'ERK1/2' and 'PDPK1'
  conditions <- c('mutations', inhibited,
                  'ERK1_2', 'KRAS', 'PDPK1')
  
  # For each condition 
  for (condition in conditions){
    
    # SET THE INITIAL STATE OF EACH CONDITION
    if(condition == 'mutations'){
      
      # set stimuli conditions as a named vector of 1
      # and assign it to CNOlist@stimuli
      
      # --> TNF and IGF1 receptors
      d <- c(1,1)
      names(d) <- c('TNF', 'IGF1R')
      
      # --> GOF mutations
      s <- mutations_vector[mutations_vector == 1]
      sp <- c(s,d)
      sf <- rep(1, length(sp))
      names(sf) <- names(sp)
      CNOlist@stimuli <-  t(as.matrix(sf))
      
      # Set inhibitory conditions as a named vector of 1
      # and assign it to CNOlist@inhibitors
      # --> LOF mutations
      ini <- mutations_vector[mutations_vector == 0]
      i_new <- rep(1, length(ini))
      names(i_new) <- names(ini)
      
      # Modify the built-in object
      CNOlist@inhibitors <- t(as.matrix(i_new))
      
    }else if(condition == 'FLT3'){ # set stimuli conditions
      # --> TNF and IGF1 receptors
      d <- c(1,1)
      names(d) <- c('TNF', 'IGF1R') 
      
      # --> all the GOF 
      # BUT we exclude FLT3 since we inhibit it!
      s <- mutations_vector[mutations_vector == 1]
      sp <- c(s,d)
      sf <- rep(1, length(sp))
      names(sf) <- names(sp)
      sf <- sf[names(sf) != 'FLT3'] # exclude FLT3
      
      # Modify the built-in object
      CNOlist@stimuli <-  t(as.matrix(sf))
      
      # Set inhibitory condition 
      # --> FLT3i 
      in_val <- c(1)
      names(in_val) <- condition
      
      # --> LOF mutations
      i <- mutations_vector[mutations_vector == 0]
      i_new <- rep(1, length(i))
      names(i_new) <- names(i)
      
      # combine the two
      i_trans <- c(i_new, in_val)
      names(i_trans) <- c(names(i_new), names(in_val))
      
      # modify the built-in object
      CNOlist@inhibitors <- t(as.matrix(i_trans))}else{
        ## MUTATIONS + FLT3i + NODES ##
        # set stimuli conditions
        # --> TNF and IGF1 receptor
        d <- c(1,1)
        names(d) <- c('TNF', 'IGF1R') 
        
        # --> GOF mutations excluding FLT3
        s <- mutations_vector[mutations_vector == 1]
        sp <- c(s,d)
        sf <- rep(1, length(sp))
        names(sf) <- names(sp)
        sf <- sf[names(sf) != 'FLT3'] #exclude FLT3
        CNOlist@stimuli <-  t(as.matrix(sf)) # modify the built-in object
        
        # set inhibitory conditions
        # --> LOFs mutations 
        i <- mutations_vector[mutations_vector == 0]
        i_new <- rep(1, length(i))
        names(i_new) <- names(i)
        
        # --> FLT3i and new druggable node specified in condition 
        in_val <- c(1, 1)
        names(in_val) <- c('FLT3', condition)
        
        i_trans <- c(i_new, in_val) #combine them
        names(i_trans) <- c(names(i_new), names(in_val))
        
        CNOlist@inhibitors <- t(as.matrix(i_trans))
      }
    
    # RUN THE SIMULATION
    
    matrix_sim_q <- list()
    
    # Select the cellular model corresponding to the patient mutation type
    if(!grepl('^TKD', patient)){
      #print('selecting JMD model')
      bString <- JMD_best
      modelT0 <- JMD_model
      simList <- JMDsimList
    }else{
      #print('Selecting TKD model')
      bString <- TKD_best
      modelT0 <- TKD_model
      simList <- TKDsimList
    }
    
    modelCut <- cutModel(modelT0, list(bString)[[1]])  # cutting expanded gates to optimal bString
    newSimList <- cutSimList(simList, list(bString)[[1]])  # creating a new simList from cutted model
    indexList <- indexFinder(CNOlist, modelT0)  # creating index list
    
    # T0 simulation (before treatment)
    simResT0 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList,
                            indexList=indexList, mode = 0)
    colnames(simResT0) <- modelT0$namesSpecies
    
    # T1 simulation (after treatment)
    simResT1 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList, 
                            indexList=indexList, mode = 1)
    colnames(simResT1) <- modelT0$namesSpecies
    
    simRes <- list(T0 = simResT0, T1 = simResT1) # store results in a list
    
    # If there is some oscillatory node, set to 50
    simRes$T0[is.na(simRes$T0)] <- 50
    simRes$T1[is.na(simRes$T1)] <- 50
    
    # Create a matrix of the simulation result
    sim_matrix <- tibble(nodes = colnames(simRes$T0), T0 = simRes$T0[1,], T1 = simRes$T1[1,]) 
    
    colnames(sim_matrix) <- c('nodes', 'T0', paste0('FLT3i+',condition,'i'))
    
    sim_matrix <- sim_matrix %>%
      column_to_rownames('nodes') %>% as.matrix
    
    if(condition == 'mutations'){
      sim_matrix_new <- sim_matrix
      conditions_names <- c('T0', 'FLT3-ITD+TNF+IGF')
      colnames(sim_matrix_new) <- conditions_names
    }else{
      conditions_names <- colnames(sim_matrix_new)
      
      sim_matrix_new <- cbind(sim_matrix_new, sim_matrix[,2])
      if(condition == 'FLT3'){
        colnames(sim_matrix_new) <- c(colnames(sim_matrix_new)[1:length(colnames(sim_matrix_new))-1], 
                                      paste0(condition,'i'))
      }else{
        colnames(sim_matrix_new) <- c(colnames(sim_matrix_new)[1:length(colnames(sim_matrix_new))-1], 
                                      paste0('FLT3i+',condition,'i'))
      }
    }
  }
  
  # Remove the T0 column (not important)
  sim_matrix_new <- sim_matrix_new[,-1] 
  
  # Store simulation result in a tibble
  sim_tibble_p <- tibble(rownames_to_column(data.frame(sim_matrix_new), 'node'))
  
  # Remove combinatorial treatments with oscillatory outcome
  sim_tibble_p <- sim_tibble_p[,colSums(sim_tibble_p == 50) == 0]
  
  # Create a file with model nodes' states for each patient
  write_tsv(sim_tibble_p, 
            paste0('./results/patients_simulations/', patient, '_steady_states.tsv'))
  
  # PHENOTYPES INFERENCE ACTIVITY
  # Infer proliferation and apoptosis activity in all conditions
  phens <- c('APOPTOSIS', 'PROLIFERATION')
  
  for(phen_i in 1:length(phens)){
    phen <- phens[phen_i]
    
    # read ProxPath output 
    apoptosis_annotation <- read_tsv('./results/phenotypes_analysis/paths_to_apoptosis_proliferation.txt') %>%
      select(node = QueryNode, phenotype = EndNode, apoptosis = Final_Effect)
    
    # select the phenotype and remove duplicated nodes
    apoptosis_annotation <- apoptosis_annotation %>% 
      filter(phenotype == phen) %>% 
      select(node, apoptosis)  %>%
      distinct()

    # Join nodes activity in simulation with the regulatory role
    # on the phenotype 
    
    if(grepl('JMD', patient)){
      annotated <- left_join(sim_tibble_p, apoptosis_annotation %>% 
                               select(node, apoptosis), by = 'node')
      
      # keep only terminal nodes in the signal cascade
      if(phen == 'APOPTOSIS'){
        annotated <- annotated %>% 
          filter(node %in% c('MTOR', 'STAT5A'))
      }else{
        annotated <- annotated %>% 
          filter(node %in% c('PI3K', 'MAPK14', 'STAT5A', 'MTOR'))
      }
    }else{
      annotated <- left_join(sim_tibble_p, apoptosis_annotation %>% 
                               select(node, apoptosis), by = 'node')
      
      if(phen == 'APOPTOSIS'){
        annotated <- annotated %>%
          filter(node %in% c('STAT5A', 'AKT'))
      }else{
        annotated <- annotated %>%
          filter(node %in% c('PI3K', 'KRAS', 'CREB1', 'STAT5A'))
      }
    }
    
    # modify the table for the phenotype barplot
    annotated <- annotated %>% column_to_rownames('node')
    annotated <- annotated*annotated$apoptosis
    annotated$apoptosis <- NULL
    
    pheno_df <- tibble(condition = str_replace_all(names(colSums(annotated)), '\\.', '+'), 
                       value = colSums(annotated))
    
    #revert the APOPTOSIS sign to create 'apoptosis inhibition'
    if(phen == 'APOPTOSIS'){ 
      pheno_df$value <- pheno_df$value * -1
    }
    
    
    barplot_apoptosis <- ggplot(pheno_df, aes(x = condition, y = value, fill)) +
      geom_bar(stat = 'identity', position = 'dodge', width = 0.8) +
      #labs(title = str_replace(conditions_sim, '\\.', '+'))+
      xlab('') + ylab(phen)+
      theme_classic() +
      ylim(0,5)+
      ggtitle(paste0(patient, ' ', phen))+
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 10),
            axis.text.y = element_text(size = 10), 
            axis.title = element_text(size = 12))
    
    patients_plot_list[[plot_idx]] <- barplot_apoptosis
    plot_idx = plot_idx + 1
    print(plot_idx)
    
    # create a table for proliferation and apoptosis phenotype for heatmap
    if(p_idx == 1){
      if(phen == 'APOPTOSIS'){
        apoptosis_df <- pheno_df 
        colnames(apoptosis_df) <- c('condition', patient)
      }else{
        proliferation_df <- pheno_df
        colnames(proliferation_df) <- c('condition', patient)
      }
    }else{
      if(phen == 'APOPTOSIS'){
        colnames(pheno_df) <- c('condition', patient)
        apoptosis_df <- full_join(apoptosis_df, pheno_df, by = 'condition') 
      }else{
        colnames(pheno_df) <- c('condition', patient)
        proliferation_df <- full_join(proliferation_df, pheno_df, by = 'condition') 
      }
    }
  }
  patients_plot_list[[plot_idx]] <- plot_list
}

# save df for apoptosis and proliferation levels

# arrange all the plots in a grid
combined_plots <- plot_grid(plotlist = patients_plot_list, labels = "AUTO", ncol=8,nrow = 4)
write_rds(combined_plots, './results/patients_simulations/barplot_combinations.rds')
write_tsv(apoptosis_df, './results/patients_simulations/apoptosis_combination_patients.tsv')
write_tsv(proliferation_df, './results/patients_simulations/proliferation_combination_patients.tsv')

Phenotypes activity of each patient model

combined_plots <- readRDS('./results/patients_simulations/barplot_combinations.rds')
combined_plots

Figure S6E Each couple of bar plots show the in-silico apoptosis inhibition (right) and proliferation activation (left) levels in control and FLT3 inhibition conditions in combination with the knock-out of each of 10 kinases in FLT3 ITD-JMD and FLT3 ITD-JMD + TKD and FLT3 ITD-TKD patients.

Heatmap of phenotypes

i_plot = 1
plot_list <- list()
for(phenotype in c('apoptosis', 'proliferation')){
  
  pheno_tb <- read_tsv(paste0('./results/patients_simulations/', phenotype, '_combination_patients.tsv'))
  
  pheno_matrix <-  pheno_tb %>% column_to_rownames('condition') 
  
  if(phenotype == 'apoptosis'){
    myColor <- c('#557887', '#FFE2E2', 'red3')
  }else{
     myColor <- colorRampPalette(c("#FFE2E2", "red3"))(4)
  }
  
  p <- pheatmap(as.matrix(pheno_matrix),
              color = myColor,
              cellwidth = 20,
              cellheight = 20,
              cluster_rows = F,
              cluster_cols = T,
              border_color = 'white',
              fontsize = 13,
              width = 5,
              heigth = 30,
              silent = TRUE,
              main = paste0('in silico ', phenotype, 
                            ifelse(phenotype == 'apoptosis', 
                                   ' inhibition ', ' '), 
                            'levels'),
              angle_col = c('90'))
  
  plot_list[[i_plot]] <- p[[4]]
  i_plot = i_plot + 1
}

plot_grid(plotlist = plot_list, ncol = 2)

Figure 5E Heatmap representing patient-specific in silico apoptosis inhibition (left panel) and proliferation levels (right panel) upon each simulation condition.