Although there is a vast number of cancer medication which is developed, but the most common problem is the resistance of the cancer cell to a single treatment. The future of the anticancer treatment will be direct more to the drug combination between two or more of the anticancer medications, but it’s not easy to test every possible drug combination experimentally (1,2,3).
Using the knowledge available in different literatures and databases about cancer drugs targets and the ability to create a computer model that represent the different cancer cell lines and make it easy and a promising field for developing new drugs combinations without spending a long time in testing all the possible combinations (1).
MEK, AKT, TAK and PI3K are some of the targets that participate in more than one pathway and used as a drug target for treatment of different cancer types (4,5,6,7), but the single therapy can lead to drug resistance in addition to the side effects due to the use of high dose of one drug (8). In contrast, the use of a combination between two or more therapeutics lead to decreasing the expected side effects due to decreasing the drug concentration, having better therapeutic effect because of targeting several pathways and decreasing the drug resistance (9).
Our model based on Flobak et al. 2015 gastric adenocarcinoma cell line (AGS) logical model. The original model created using GINsim depending on the prior biological knowledge of the AGS intracelleular signaling pathways, and the final network consists of 75 components and 149 direct interactions. The model used to predect the synergy between seven inhibitors each inhibitor targets a specific node in the network which leads to 21 pairwise combinations. The result was that the model confirm the synergy between the previously discovered combinations MEK-AKT or MEK-PI3K inhibitions in addition to discovering a new synergy between TAK1-AKT or TAK1-PI3K inhibitions combinations (1).
The aims of the project is to verify the results from Flobak et al. 2015 using bioLQM perturpation and MaBoSS Stochastic simulations, in addition to using pint and the different tools to find a new drug targets that have synergy in the model.
GINsim is a free software used to facilitates the definition of multi-valued logical models, modeling, and analysis of regulatory networks (11).
# First GINsim was used to load and show the original model.
import ginsim
lrg = ginsim.load("/tmp/colomotoxcpb9yzb_Flobak_FullModel_S2_Dataset.zginml")
ginsim.show(lrg)
Figure 1. The original model graph of 75 components represented by the gray and orange (targets) colors and the two outputs colored with pink for antisurvival and green for prosurvival and 149 direct interactions (the green arrows represent activatory interactions and red T arrows represent inhibitory interactions).
Using bioLQM for conversion of the multi-valued model to booleanized model to make it easy to use the boolean models tools, in addition, identificcation of stable states, the states where there is no updates happend to the network components and each component can continue to operate at its current level (10).
import biolqm
# Conversion of our ginsim model to bioLQM model to be able to compute the stable states of the model using Java library bioLQM.
lqm = ginsim.to_biolqm(lrg)
fixpoints = biolqm.fixpoints(lqm)
# To show how the stable states in our model.
print(len(fixpoints), "fixpoints")
# To display the list of stable states
from colomoto_jupyter import tabulate
tabulate(fixpoints)
1 fixpoints
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival | cMYC | TCF | NLK | DKK1gene | CCND1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 3 | 1 | 1 | 1 | 0 | 2 |
Figure 2. Shows the model stable state with the result of antisurvival=0 (inactive) and prosurvival=3 (active) and we can notice that MEK,TAK,AKT and PI3K nodes which used as targets in the model are active in the stable state.
# GINsim can be used to visualize The stable state
ginsim.show(lrg,fixpoints[0])
Figure 3. GINsim regulatory graph showing the result of stable state. The white node represent the value 0 means inactive while the blue nodes represent value 1 means active.
The aim of the next section is to confirm the Flobak et al.2015 results which show that there are four combinations have synergies predected by the model, two combinationes was known MEK-AKT or MEK-PI3K inhibitions, in addition to two new combinations TAK1-AKT or TAK1-PI3K inhibitions (1). We will confirm these results using both bioLQM perturbation and MaBoSS stochastic simulation.
BioLQM perturbation (mutation) used to knockout a specific component in the model by fixing the activity level of this component to 0. This type of perturbation has a role in evaluating the importance of the component in the model (10).
pert_AKT = biolqm.perturbation(lqm, "AKT%0")
# define stable state
fps_AKT = biolqm.fixpoints(pert_AKT)
tabulate(fps_AKT)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival | cMYC | TCF | NLK | DKK1gene | CCND1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 3 | 1 | 1 | 1 | 0 | 2 |
The result from the single perturbation of AKT node showed a growth 3-1 = 2
pert_TAK1 = biolqm.perturbation(lqm, "TAK1%0")
fps_TAK1 = biolqm.fixpoints(pert_TAK1)
tabulate(fps_TAK1)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival | cMYC | TCF | NLK | DKK1gene | CCND1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 3 | 1 | 1 | 0 | 0 | 2 |
The result from the single perturbation of TAK1 node showed a growth 3-0 = 3
# By Knoking out both AKT and TAK1 nodes we can notice the effect on the model stable state.
pert_AKT_TAK1 = biolqm.perturbation(lqm, "AKT%0 TAK1%0")
# define stable state
fps_AKT_TAK1 = biolqm.fixpoints(pert_AKT_TAK1)
tabulate(fps_AKT_TAK1)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival | cMYC | TCF | NLK | DKK1gene | CCND1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 2 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 3 | 1 | 1 | 0 | 0 | 2 |
Figure 4. The stable state of the model after knoking out both AKT and TAK1 shows Antisurvival = 2 and Prosurvival = 3
import maboss
mbs_AKT_TAK1 = biolqm.to_maboss(lqm)
# To check that the two nodes are the desired we simulate mutations.
mbs_AKT_TAK1.mutate("AKT", "OFF")
mbs_AKT_TAK1.mutate("TAK1", "OFF")
# Spicify the outputs
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_AKT_TAK1,output)
# Setting the parameters
mbs_AKT_TAK1.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
# Run and visualize simulation results
AKT_TAK1_sim = mbs_AKT_TAK1.run()
AKT_TAK1_sim.plot_trajectory()
AKT_TAK1_sim.plot_piechart(autopct=True)
Figure 5. Shows the result of knocking out both AKT and TAK1 nodes and the two outputs are multi-value nodes so in this state we noticed that Antisurvival_b1, b2 are active (means 1) so the value = 2, while prosurvival_b1, b2, b3 are active so the value = 3.
pert_AKT_MEK = biolqm.perturbation(lqm, "AKT%0 MEK%0")
# define stable state
fps_pert_AKT_MEK = biolqm.fixpoints(pert_AKT_MEK)
tabulate(fps_pert_AKT_MEK)
# Trapspaces
traps_pert_AKT_MEK = biolqm.trapspace(pert_AKT_MEK)
tabulate(traps_pert_AKT_MEK)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 254 | 1 | 1 | 1 | 1 | 1 | 0 | 254 | 254 | 254 | 1 | 1 | 254 | 254 | 254 | 1 | 254 | 1 | 1 | 1 | 254 | 254 | 254 | 254 | 254 | 0 | 0 | 0 | 1 | 0 | 254 | 254 | 0 | 1 | 254 | 254 | 1 | 0 | 0 | 254 | 0 | 254 | 1 | 0 | 254 | 254 | 254 | 254 | 1 | 1 | 1 | 254 | 0 | 1 | 1 | 254 | 254 | 0 | 0 | 254 | 254 | 254 | 0 | 0 | 1 | 0 | 0 | 1 | 254 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 254 | 1 | 1 | 254 | 0 | 1 | 254 |
Figure 6. Perturbation of AKT and MEK nodes give cyclic attractors so we run a trapspace and in this case we found that Antisurvival_b1 is active with value 1, while both Antisurvival_b2 and 3 are oscillate between 0 and 1 so we gave the value 0.5 the same with Prosurvival_b3 so the final values are:
mbs_AKT_MEK = biolqm.to_maboss(lqm)
# To check that the two nodes are the desired we simulate mutations.
mbs_AKT_MEK.mutate("AKT", "OFF")
mbs_AKT_MEK.mutate("MEK", "OFF")
# Spicify the outputs
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_AKT_MEK,output)
# Setting the parameters
mbs_AKT_MEK.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
# Run and visualize simulation results
AKT_MEK_sim = mbs_AKT_MEK.run()
AKT_MEK_sim.plot_trajectory()
AKT_MEK_sim.plot_piechart(autopct=True)
Figure 7. shows the results of AKT and MEK mutation by visualizing the chart for example the orange part 33.2% have Antisurvival_b1 as active while b2 and b3 are inactive so this gives Antisurvival value = 1 while the Prosurvival_b1,b2,b3 are active in the same state so it has the value 3 so the growth predicted to be 2 and 30 % for growth = 1 .
pert_PI3K_MEK = biolqm.perturbation(lqm, "PI3K%0 MEK%0")
# define stable state
fps_pert_PI3K_MEK = biolqm.fixpoints(pert_PI3K_MEK)
tabulate(fps_pert_PI3K_MEK)
# Trapspaces
traps_pert_PI3K_MEK = biolqm.trapspace(pert_PI3K_MEK)
tabulate(traps_pert_PI3K_MEK)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 254 | 1 | 1 | 1 | 1 | 1 | 0 | 254 | 254 | 254 | 1 | 1 | 254 | 254 | 254 | 1 | 254 | 1 | 1 | 1 | 254 | 254 | 254 | 254 | 254 | 0 | 0 | 254 | 1 | 0 | 254 | 254 | 0 | 1 | 0 | 254 | 1 | 0 | 0 | 1 | 0 | 254 | 0 | 254 | 254 | 254 | 254 | 254 | 0 | 0 | 1 | 254 | 0 | 1 | 1 | 254 | 254 | 0 | 0 | 254 | 0 | 254 | 0 | 0 | 1 | 0 | 0 | 1 | 254 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 254 | 0 | 1 | 0 |
Figure 8. The same happened with PI3K and MEK perturbation using bioLQM and having cyclic attractor and the results are:
mbs_PI3K_MEK = biolqm.to_maboss(lqm)
# To check that the two nodes are the desired we simulate mutations.
mbs_PI3K_MEK.mutate("PI3K", "OFF")
mbs_PI3K_MEK.mutate("MEK", "OFF")
# Spicify the outputs
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_PI3K_MEK,output)
# Setting the parameters
mbs_PI3K_MEK.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
# Run and visualize simulation results
PI3K_MEK_sim = mbs_PI3K_MEK.run()
PI3K_MEK_sim.plot_trajectory()
PI3K_MEK_sim.plot_piechart(autopct=True)
Figure 9. Shows the stochastic simulation of PI3K and MEK.
pert_PI3K_TAK1 = biolqm.perturbation(lqm, "PI3K%0 TAK1%0")
# define stable state
fps_pert_PI3K_TAK1 = biolqm.fixpoints(pert_PI3K_TAK1)
tabulate(fps_pert_PI3K_TAK1)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival | cMYC | TCF | NLK | DKK1gene | CCND1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 2 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 2 | 1 | 1 | 0 | 0 | 1 |
mbs_PI3K_TAK1 = biolqm.to_maboss(lqm)
# To check that the two nodes are the desired we simulate mutations.
mbs_PI3K_TAK1.mutate("PI3K", "OFF")
mbs_PI3K_TAK1.mutate("TAK1", "OFF")
# Spicify the outputs
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_PI3K_TAK1,output)
# Setting the parameters
mbs_PI3K_TAK1.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
# Run and visualize simulation results
PI3K_TAK1_sim = mbs_PI3K_TAK1.run()
PI3K_TAK1_sim.plot_trajectory()
PI3K_TAK1_sim.plot_piechart(autopct=True)
Figure 11. Shows the stochastic simulation of PI3K and TAK1 and the results were Antisurvival 2 and Prosurvival 2.
The aim of model reduction is to facilitate the analysis of the large model by creating a reduced model which contains fewer components with the similar dynamic characteristics (10). In the reduction step we used GINsim software to reduce the model, and we decided to use the reduced model from Flobak et al.2015. In addition keeping mTORC1 node as a new target as it has no direct connection to the outputs or directly linked to either another drug target In our reduced model DUSP6 was added as it's a self-loop node in the model.
red_model = ginsim.load("/tmp/colomotox_ff8atw_CASCADE_reduction.zginml")
ginsim.show(red_model)
Figure 12. Represent the reduced model after selecting mTORK1 to be a new node to be added to Flobak et al 2015 reduced model as it has no direct connection to the outputs nodes and in this model DUSP6 was kept because its a self regulated node.
# First we checked if the reduced model has the same stable state as the original model using bioLQM fixpoints.
lqm_red = ginsim.to_biolqm(red_model) # Conversion of reduced model from GINsim to bioLQM
fps_red = biolqm.fixpoints(lqm_red)
print(len(fps_red), "fixpoints")
tabulate(fps_red)
1 fixpoints
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 3 |
ginsim.show(red_model, fps_red[0])
Figure 13. GINsim regulatory graph showing the result of stable state. The white node represent the value 0 means inactive while the blue nodes represent value 1 means active.
BioLQM and MaBoSS simulation used to compare the perturbation of each single node and the combination between mTORC1 and each node in the model to see if there is a synergy between any of these combinations.
import sys
import itertools
perturbations_single = ["MEK%0","TAK1%0","p38alpha%0","AKT%0","PI3K%0","GSK3%0", "betacatenin%0"]
'''Define dictionaries to store attractors for each perturbation'''
fixpointlist = {}
'''trapspacelist is used for perturbations where no stable state is found
to find eventual cyclic attractors.'''
trapspacelist = {}
for p in perturbations_single:
'''No stable state found, so we look for cycles'''
trapspace = biolqm.trapspace(biolqm.perturbation(lqm_red, p))
trapspacelist[p] = trapspace
for k, v in trapspacelist.items():
print(k)
sys.displayhook(tabulate(v))
print()
MEK%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 254 | 254 | 1 | 1 | 1 | 254 | 0 | 0 | 0 | 1 | 1 | 1 | 0 |
TAK1%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
p38alpha%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
AKT%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
PI3K%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 |
GSK3%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
betacatenin%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
Figure 13. The results of bioLQM single perturbation of the reduced model.
# One example of MaBoSS single mutation using MEK.
mbs_red_MEK = biolqm.to_maboss(lqm_red)
mbs_red_MEK.mutate("MEK","OFF")
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_red_MEK,output)
mbs_red_MEK.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
MEK_sim = mbs_red_MEK.run()
MEK_sim.plot_trajectory()
MEK_sim.plot_piechart(autopct=True)
Figure 14. MaBoSS simulation of MEK node using the reduced model showed that ~ 50% propability that the model will reach growth of 1 due to activation of (Antisurvival b1 ) and (prosurvival_b1 and b2).
To confirm the efficacy of our reduced model by comparing the perturbation results with the original model.
perturbations_single = ["MEK%0","TAK1%0","p38alpha%0","AKT%0","PI3K%0","GSK3%0", "betacatenin%0"]
'''Define dictionaries to store attractors for each perturbation'''
fixpointlist = {}
'''trapspacelist is used for perturbations where no stable state is found
to find eventual cyclic attractors.'''
trapspacelist = {}
for p in perturbations_single:
'''No stable state found, so we look for cycles'''
trapspace = biolqm.trapspace(biolqm.perturbation(lqm, p))
trapspacelist[p] = trapspace
for k, v in trapspacelist.items():
print(k)
sys.displayhook(tabulate(v))
print()
MEK%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 254 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 254 | 254 | 0 | 254 | 254 | 254 | 1 | 254 | 1 | 1 | 1 | 254 | 254 | 1 | 1 | 254 | 1 | 1 | 0 | 254 | 0 | 0 | 0 | 0 | 0 | 1 | 254 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 254 | 1 | 254 | 1 | 1 | 1 | 254 | 0 | 1 | 254 | 0 | 0 | 0 | 0 | 0 | 0 | 254 | 254 | 0 | 1 | 0 | 0 | 254 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 254 | 0 | 1 | 0 |
TAK1%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 |
p38alpha%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 |
AKT%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 |
PI3K%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 0 |
GSK3%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 |
betacatenin%0
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
Figure 14. The single perturbation of the original model to compare the results with our reduced model.
# Results from the single perturbation of the nodes in the reduced and the original model are the same.
#Reduced model original model
# MEK growth = 1.5 1.5
# TAK1 growth = 3 3
# p38alpha growth = 3 3
# AKT growth = 2 2
# PI3K growth = 1 1
# GSK growth = 3 3
# betacatenin growth = -1 -1
To see if there is any synergy between our new node (mTORC1) in the reduced model and any of the other nodes in the model we run a double perturbation between mTORC1 and all the nodes in the reduced model using bioLQM and MaBoSS.
perturbations = ["MEK%0","TAK1%0","p38alpha%0","AKT%0","PI3K%0","GSK3%0", "betacatenin%0"]
# Perform double perturbations - All nodes against one
per = ["mTORC1%0"] # Node that should be in combinatorial perturbations
Node_perturbations = list(itertools.product(per, perturbations))
Node_perturbations = str(Node_perturbations).replace(',', ' ')
Node_perturbations
"[('mTORC1%0' 'MEK%0') ('mTORC1%0' 'TAK1%0') ('mTORC1%0' 'p38alpha%0') ('mTORC1%0' 'AKT%0') ('mTORC1%0' 'PI3K%0') ('mTORC1%0' 'GSK3%0') ('mTORC1%0' 'betacatenin%0')]"
'''Define dictionaries to store attractors for each perturbation'''
fixpointlist = {}
'''trapspacelist is used for perturbations where no stable state is found
to find eventual cyclic attractors.'''
trapspacelist = {}
for p in perturbations:
fixpoints = biolqm.fixpoints(biolqm.perturbation(lqm_red, p))
if(fixpoints):
fixpointlist[p] = fixpoints
else:
'''No stable state found, so we look for cycles'''
trapspace = biolqm.trapspace(biolqm.perturbation(lqm_red, p))
trapspacelist[p] = trapspace
#see attractors in tabulated tables for each perturbation
for k, v in fixpointlist.items():
print(k)
sys.displayhook(tabulate(v))
print()
for k, v in trapspacelist.items():
print(k)
sys.displayhook(tabulate(v))
print()
TAK1%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 3 |
p38alpha%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 3 |
AKT%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 3 |
PI3K%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 2 |
GSK3%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 3 |
betacatenin%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival | GSK3 | betacatenin | Prosurvival | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
MEK%0
MEK | DUSP6 | TAK1 | p38alpha | AKT | mTORC1 | PI3K | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | GSK3 | betacatenin | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 1 | 254 | 254 | 1 | 1 | 1 | 254 | 0 | 0 | 0 | 1 | 1 | 1 | 0 |
Figure 15. Results of double perturbation between mTORC1 and all the nodes in the reduced model, and by comparing the results with the sigle perturbation showed that there are no synergys between any of these combinations.
In order not to repeat the same simulation twice we run bioLQM double perturpation in the reduced model and MaBoSS simulation in the original model also to see if there a synergy between mTORC1 and any of the nodes in the original model.
# One example of double perturbation between mTORC1 and MEK.
mbs_mTORC1_MEK = biolqm.to_maboss(lqm)
# To check that the two nodes are the desired we simulate mutations.
mbs_mTORC1_MEK.mutate("mTORC1", "OFF")
mbs_mTORC1_MEK.mutate("MEK", "OFF")
# Spicify the outputs
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_mTORC1_MEK,output)
# Setting the parameters
mbs_mTORC1_MEK.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
# Run and visualize simulation results
mTORC1_MEK_sim = mbs_mTORC1_MEK.run()
mTORC1_MEK_sim.plot_trajectory()
mTORC1_MEK_sim.plot_piechart(autopct=True)
Figure 16. MaBoSS simulation between mTORC1 and the nodes in the original model shows that ~ 40% propability that the model will reach growth of 3 activation of (prosurvival_b1,b2,b3).
Although there are synergys between the nodes in the reduced model of Flobak et al 2015, but we couldn't find any synergy between mTORC1 and any of these nodes so we tried to use Pint to predict a new possibilities of combinations in the original model that may have synergies.
import pypint
pint_model = pypint.load("/tmp/colomototoaptmmw_Flobak_FullModel_S2_Dataset.zginml")
1 state(s) have been registered: initState_
# Using (having) method to Define the new initial state to our model here for exampe MEK and PI3K
model_initial = pint_model.having({"MEK": 0, "PI3K": 0})
# Then we specified the goal.
from pypint import Goal
alt = Goal("Antisurvival=2", "Prosurvival=3") | Goal("Antisurvival=1", "Prosurvival=2") | Goal("Antisurvival=0", "Prosurvival=1")
# To see if the goal is reachable in the original model initial state or not
model_initial.reachability(alt)
True
# To make sure that our modified initial state reach the goal or not
model_initial.having(MEK=0 , PI3K=0).reachability(alt)
True
# Pint reachability analysis can be used to identify mutations that control goal reachability.
mutations = model_initial.oneshot_mutations_for_cut(alt,maxsize=2)
mutations
This computation is an under-approximation: returned mutations are all valid, but they may be non-minimal, and some solutions may be missed.
Limiting solutions to mutations of at most 2 automata. Use maxsize
argument to change.
[{'ERK': 0, 'TCF': 0}, {'PTEN': 1, 'TCF': 0}, {'PI3K': 0, 'TCF': 0}, {'PDK1': 0, 'TCF': 0}, {'RSK': 0, 'TCF': 0}, {'cMYC': 0, 'CCND1': 0}, {'TCF': 0, 'CCND1': 0}]
First we run perturbation for TCF alone to identify the growth( output)
mbs_TCF = biolqm.to_maboss(lqm)
mbs_TCF.mutate("TCF","OFF")
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_TCF,output)
mbs_TCF.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
TCF_sim = mbs_TCF.run()
TCF_sim.plot_trajectory()
TCF_sim.plot_piechart(autopct=True)
Figure 17. Represent the single perturbation result of TCF in the original model
MaBoSS simulatio used to confirm the predicted results from Pint
mbs_TCF_PI3K = biolqm.to_maboss(lqm)
mbs_TCF_PI3K.mutate("TCF","OFF")
mbs_TCF_PI3K.mutate("PI3K","OFF")
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_TCF_PI3K,output)
mbs_TCF_PI3K.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
TCF_PI3K_sim = mbs_TCF_PI3K.run()
TCF_PI3K_sim.plot_trajectory()
TCF_PI3K_sim.plot_piechart()
Figure 18. MaBoSS simulation result of double mutation of TCF and PI3K.
Triple mutation was used to see if there is a synergy if we knocked out MEK with the double mutation of TCF and PI3K. I used both BioLQm to identify the cyclic attractor and MaBoSS simulation.
# BioLQM triple perturbation
pert_MEK_PI3K_TCF = biolqm.perturbation(lqm, "MEK%0 PI3K%0 TCF%0")
# define stable state
fps_pert_MEK_PI3K_TCF = biolqm.fixpoints(pert_MEK_PI3K_TCF)
tabulate(fps_pert_MEK_PI3K_TCF)
# Trapspaces
traps_pert_MEK_PI3K_TCF = biolqm.trapspace(pert_MEK_PI3K_TCF)
tabulate(traps_pert_MEK_PI3K_TCF)
RTPK | SHC1 | GRB2 | SOS | Ras | Raf | MEK | ERK | DUSP6 | TAK1 | MKK3 | ASK1 | TAB | p38alpha | DUSP1 | MKK4 | MKK7 | MEKK4 | MLK3 | Rac | JNK | GRAP2 | SHP2 | GAB | IKKB | IKKA | AKT | IRS1 | Caspase8 | PTENgene | Caspase9 | FOXO | TSC | pras40 | BCL2 | p53 | Rheb | mTORC1 | PTEN | BAD | mTORC2 | MDM2gene | S6K | MDM2 | CytochromeC | BAX | NFkB | MAP3K8 | PI3K | PDK1 | LRP | ITCH | DKK1 | Fz | Antisurvival_b1 | Antisurvival_b2 | Antisurvival_b3 | Axin | GSK3 | RTPKgene | RSK | MSK | CFLAR | CK1 | Dvl | SFRP1 | SFRP1gene | Caspase37_b1 | Caspase37_b2 | betacatenin | betaTrCP | MMP | Egr1 | LEF | Prosurvival_b1 | Prosurvival_b2 | Prosurvival_b3 | cMYC | TCF | NLK | DKK1gene | CCND1_b1 | CCND1_b2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 254 | 254 | 254 | 254 | 254 | 254 | 0 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 1 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 254 | 0 | 0 | 254 | 254 | 0 | 0 | 254 | 254 | 254 | 254 | 254 | 254 | 0 | 254 | 254 | 254 | 0 | 1 | 1 | 254 | 254 | 254 | 254 | 254 | 1 | 254 | 0 | 0 | 0 | 0 | 0 | 254 | 0 | 0 | 0 |
Figure 19. Shows the results of trible perturbation using bioLQM. Antisurvival = 1.5, prosurvival = 0 so the growth will be - 1.5 which is less than the growth of each node perturbation alone so there is synergy in triple perturbation.
# MaBoSS simulation was used to confirm the results.
mbs_MEK_PI3K_TCF = biolqm.to_maboss(lqm)
mbs_MEK_PI3K_TCF.mutate("MEK","OFF")
mbs_MEK_PI3K_TCF.mutate("PI3K","OFF")
mbs_MEK_PI3K_TCF.mutate("TCF","OFF")
output = ['Prosurvival_b1','Prosurvival_b2','Prosurvival_b3', 'Antisurvival_b1', 'Antisurvival_b2','Antisurvival_b3']
maboss.set_output(mbs_MEK_PI3K_TCF,output)
mbs_MEK_PI3K_TCF.update_parameters(sample_count=1000,time_tick=1.0, max_time=100)
MEK_PI3K_TCF_sim = mbs_MEK_PI3K_TCF.run()
MEK_PI3K_TCF_sim.plot_trajectory()
MEK_PI3K_TCF_sim.plot_piechart(autopct=True)
Figure 20. MaBoSS simulation for triple perturbation shows that ~70% propability that the model will reach growth of -1 activation of (Antisurvival_b1) so here we can predict that TCF is a promsining drug target.
By using bioLQM perturbation and MaBoSS Stochastic simulation we confirmed the synergy between the combination of MEK or TAK1 inhibitors with PI3K or AKT inhibitors (results of Flobak et al. 2015), and this gave us a good indication about the efficacy of these tools.
In the reduced model we decided to choose a new target (mTORC1) and add it to the Flobak reduced model to see if there any synergy between this node and any of the other nodes in the reduced model, but the results after running single and double perturbation showed that there is no synergy. Pint was used to predict new combinations in the model which can lead to a synergy so by using oneshot_mutations_for_cut method and PI3K = 0 and MEK = 0 as initial state, TCF node was suggested to have a combination with PI3K and this mutation affect goal reachability, so we run double perturbation using biLQM and stochastic simulation using MaBoSS to confirm the synergy. MaBoSS simulation and bioLQm perturbation used for double mutation between TCF and PI3K but the result showed that there is no inhibition of the growth better than using TCF inhibitor alone so there is no synergy in double mutation, while in triple mutation of TCF, PI3K and MEK showed ~70% propability that the model will reach growth of -1 activation of (Antisurvival_b1) which means that there is a synergy in the triple mutation case and give an prediction that TCF could be a promising drug target.
As we mentiond that the first aim of the project is to confirm the predicted results from Flobak et al. 2015, so using both bioLQM and MaBoSS stochastic simulation we confirmed the four drug combinations. The second aim of the project is to identify a new drug targets that may have synergy in AGS model, in this step we used pint reachability analysis which suggested TCF as a new target that have triple synergy with MEK and PI3K
According to the results from Pint triple mutation, TCF inhibitor considered as a promising drug target and further analysis can be done to find mor triple combinations or by changing the initial state we can find more nodes which have a good synergy in the model.
1- Flobak Å, Baudot A, Remy E, Thommesen L, Thieffry D, et al. (2015) Discovery of Drug Synergies in Gastric Cancer Cells Predicted by Logical Modeling. PLOS Computational Biology 11(8): e1004426. https://doi.org/10.1371/journal.pcbi.1004426
2- Flobak, Å, Niederdorfer, B., Nakstad, V. T., Thommesen, L., Klinkenberg, G., and Lægreid, A. (2019). A high-throughput drug combination screen of targeted small molecule inhibitors in cancer cell lines. Sci. Data 6:237. doi: 10.1038/ s41597-019-0255-7
3- Niederdorfer, B., Touré, V., Vazquez, M., Thommesen, L., Kuiper, M., Lægreid, A., & Flobak, Å. (2020). Strategies to Enhance Logic Modeling-Based Cell Line-Specific Drug Synergy Prediction [Original Research]. Frontiers in Physiology, 11. https://doi.org/10.3389/fphys.2020.00862
4- Degirmenci U, Wang M, Hu J. Targeting Aberrant RAS/RAF/MEK/ERK Signaling for Cancer Therapy. Cells. 2020 Jan 13;9(1):198. doi: 10.3390/cells9010198. PMID: 31941155; PMCID: PMC7017232.
5- Hinz N, Jücker M. Distinct functions of AKT isoforms in breast cancer: a comprehensive review. Cell Commun Signal. 2019 Nov 21;17(1):154. doi: 10.1186/s12964-019-0450-3. PMID: 31752925; PMCID: PMC6873690.
6- Mosele F, Stefanovska B, Lusque A, Tran Dien A, Garberis I, Droin N, Le Tourneau C, Sablin MP, Lacroix L, Enrico D, Miran I, Jovelet C, Bièche I, Soria JC, Bertucci F, Bonnefoi H, Campone M, Dalenc F, Bachelot T, Jacquet A, Jimenez M, André F. Outcome and molecular landscape of patients with PIK3CA-mutated metastatic breast cancer. Ann Oncol. 2020 Mar;31(3):377-386. doi: 10.1016/j.annonc.2019.11.006. Epub 2020 Jan 24. PMID: 32067679.
7- Huang Z, Tang B, Yang Y, Yang Z, Shi L, Bai Y, Yan B, Karnes RJ, Zhang J, Jimenez R, Wang L, Wei Q, Yang J, Xu W, Jia Z, Huang H. MAP3K7-IKK Inflammatory Signaling Modulates AR Protein Degradation and Prostate Cancer Progression. Cancer Res. 2021 Sep 1;81(17):4471-4484. doi: 10.1158/0008-5472.CAN-20-4194. Epub 2021 Jun 22. PMID: 34158377.
8- Menden, M.P., Wang, D., Mason, M.J. et al. (2019) Community assessment to advance computational prediction of cancer drug combinations in a pharmacogenomic screen. Nat Commun 10, 2674. https://doi.org/10.1038/s41467-019-09799-2
9- Gregory, J. V., Vogus, D. R., Barajas, A., Cadena, M. A., Mitragotri, S., Lahann, J., Programmable Delivery of Synergistic Cancer Drug Combinations Using Bicompartmental Nanoparticles. Adv. Healthcare Mater. 2020, 9, 2000564. https://doi.org/10.1002/adhm.202000564
10- Naldi, A. (2018). bioLQM: a java library for the manipulation and conversion of Logical Qualitative Models of biological networks. bioRxiv, 287011. https://doi.org/10.1101/287011
11- Naldi A, Hernandez C, Abou-Jaoudé W, Monteiro PT, Chaouiya C and Thieffry D (2018) Logical Modeling and Analysis of Cellular Regulatory Networks With GINsim 3.0. Front. Physiol. 9:646. doi: 10.3389/fphys.2018.00646