CASCADE 1.0 MODEL ANALYSIS¶

1. Introduction.¶

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.

2. Materials and Methods¶

2.1. GINsim¶

GINsim is a free software used to facilitates the definition of multi-valued logical models, modeling, and analysis of regulatory networks (11).

In [2]:
# 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)
Out[2]:

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).

2.2. bioLQM¶

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).

In [3]:
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)

2.1.1 Identification of the stable state (fixepoints)¶

In [4]:
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
Out[4]:
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.

In [5]:
#  GINsim can be used to visualize The stable state 
ginsim.show(lrg,fixpoints[0])
Out[5]:

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.

3. Confirmation of the synergy in Flobak et al. 2015 results.¶

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.

3.A1. bioLQM confirmation of AKT & TAK1 synergy.¶

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).

Single perturbation of the synergy nodes.¶
In [6]:
 pert_AKT = biolqm.perturbation(lqm, "AKT%0")
# define stable state 
fps_AKT = biolqm.fixpoints(pert_AKT)
tabulate(fps_AKT)
Out[6]:
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
Using the Growth equation = Prosurvival - Antisurvival from flobake et al.2015¶

The result from the single perturbation of AKT node showed a growth 3-1 = 2

In [7]:
pert_TAK1 = biolqm.perturbation(lqm, "TAK1%0")
fps_TAK1 = biolqm.fixpoints(pert_TAK1)
tabulate(fps_TAK1)
Out[7]:
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

Double perturbation of the nodes AKT and TAK1 to cofirm synergy.¶
In [8]:
# 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)
Out[8]:
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

Using the equation Growth = Prosurvival - Antisurvival from flobake et al.2015, so the Growth is 3 - 2 = 1¶
This equation used to confirm the synergy ---growth (perturbation1 & perturbation2) < Min (growth (perturbation1), growth (perturbation2)).¶
So in our result here the combination growth is 1 which is lower than the growth of AKT or TAK1 perturbation alone which is 2 and 3 respectively and this confirm synergy.¶

3.A2. MaBoSS Stochastic simulation confirmation of AKT & TAK1 synergy.¶

MaBoSS library used for stochastic simulation of boolean networks, and here used to confirm the results by mutating both AKT and TAK1 nodes in the model. In the code we fixed the activity of both AKT and TAK1 "OFF" during the simulation and focused on the output nodes and finaly visulaized the results as plot and piechart.¶

In [9]:
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.

Using the equation (Growth)=Prosurvival - Antisurvival 3 - 2 = 1 this value is less than the growth of each node perturbation alone (2 for AKT and 3 for TAK1) so there is a synergy between these two nodes.¶

3.B1. bioLQM confirmation of AKT & MEK synergy¶

In [10]:
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)
Out[10]:
In [11]:
# Trapspaces

traps_pert_AKT_MEK = biolqm.trapspace(pert_AKT_MEK)
tabulate(traps_pert_AKT_MEK)
Out[11]:
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:

Antisurvival = 2¶
Prosurvival = 2.5¶
(Growth) 2.5 - 2 = 0.5 which is less than the growth of each node perturbation alone (2 for AKT and 1.5 for MEK) so there is a synergy between these two nodes.¶

3.B2. MaBoSS confirmation of AKT & MEK synergy¶

In [12]:
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 .

3.C1. bioLQM perturbation of PI3K & MEK synergy¶

In [13]:
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)
Out[13]:
In [14]:
# Trapspaces

traps_pert_PI3K_MEK = biolqm.trapspace(pert_PI3K_MEK)
tabulate(traps_pert_PI3K_MEK)
Out[14]:
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:

Antisurvival = 2¶
Prosurvival = 2¶
(Growth) 2 - 2 = 0 In this perturbation we found that in Flobak et al. 2015 the growth is 0.5 but our results also still less than the growth of each node alone (1 for PI3K and 1.5 for MEK) so there is a synergy between these two nodes.¶

3.C2. MaBoSS confirmation of PI3K & MEK synergy¶

In [15]:
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.

~ 70% probability of the growth of 1 and ~ 25% for the growth to be 0¶

3.D1. bioLQM confirmation of PI3K & TAK1 synergy.¶

In [15]:
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)
Out[15]:
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

Figure 10. shows the biolqm double perturbation of PI3K and TAK1 nodes and the results are:

Antisurvival = 2¶
Prosurvival = 2¶
(Growth) 2 - 2 = 0 The same result and less than the growth of each node perturbation alone (1 for PI3K and 3 for TAK1) so there is a synergy between these two nodes.¶

3.D2. MaBoSS confirmation of PI3K & TAK1 synergy.¶

In [16]:
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.

So the (Growth) 2 -2 = 0 the same result and confirm synergy.¶

4. Model reduction.¶

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.

In [18]:
red_model = ginsim.load("/tmp/colomotox_ff8atw_CASCADE_reduction.zginml")
ginsim.show(red_model)
Out[18]:

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.

4.1. compute the reduced model stable state.¶

In [19]:
# 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
Out[19]:
MEK DUSP6 TAK1 p38alpha AKT mTORC1 PI3K Antisurvival GSK3 betacatenin Prosurvival
0 1 1 1 0 1 1 1 0 0 1 3
The reduced model after adding mTORC1 as new target has the same stable state as the original model with MEK, TAK1, AKT and PI3K active. And the Two outputs Antisurvival = 0 (Inactive) and Prosurvival = 3 (active)¶
In [20]:
ginsim.show(red_model, fps_red[0])
Out[20]:

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.

5. Perturbation.¶

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.

5.1 Single perturbations¶

5.1a. Single perturbations in the reduced model using bioLQM.¶

In [16]:
import sys
import itertools

perturbations_single = ["MEK%0","TAK1%0","p38alpha%0","AKT%0","PI3K%0","GSK3%0", "betacatenin%0"]
In [22]:
'''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
Out[22]:
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
Out[22]:
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
Out[22]:
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
Out[22]:
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
Out[22]:
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
Out[22]:
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
Out[22]:
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.

5.1b. Single perturbations of the reduced model using MaBoSS.¶

In [23]:
# One example of MaBoSS single mutation using MEK.
mbs_red_MEK = biolqm.to_maboss(lqm_red)
In [24]:
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).

5.1c. Single perturbation of the original model.¶

To confirm the efficacy of our reduced model by comparing the perturbation results with the original model.

In [17]:
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
Out[17]:
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
Out[17]:
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
Out[17]:
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
Out[17]:
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
Out[17]:
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
Out[17]:
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
Out[17]:
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.

In [26]:
# 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

5.2. Double perturbation.¶

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.

5.2.a. Double perturbation (mTORC1) in the reduced model ( bioLQM ).¶

In [27]:
perturbations = ["MEK%0","TAK1%0","p38alpha%0","AKT%0","PI3K%0","GSK3%0", "betacatenin%0"]
In [28]:
# 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
Out[28]:
"[('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')]"
In [29]:
'''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
Out[29]:
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
Out[29]:
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
Out[29]:
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
Out[29]:
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
Out[29]:
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
Out[29]:
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
Out[29]:
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.

5.2.b. Double perturbation of (mTORC1) the original model ( MaBoSS ).¶

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.

In [30]:
# 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).

The results after calculating the growth of the double perturbation between mTORC1 and all the other nodes in both the reduced and the original models and comparing it with the growth of the single perturbation showed that there is no synergy between any of the combinations.¶

6. Pint.¶

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.

In [33]:
import pypint

pint_model = pypint.load("/tmp/colomototoaptmmw_Flobak_FullModel_S2_Dataset.zginml")

1 state(s) have been registered: initState_

In [34]:
# 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")
In [35]:
# To see if the goal is reachable in the original model initial state or not 
model_initial.reachability(alt) 
Out[35]:
True
In [36]:
# To make sure that our modified initial state reach the goal or not
model_initial.having(MEK=0 , PI3K=0).reachability(alt)
Out[36]:
True
In [37]:
# 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.

Out[37]:
[{'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}]
The list of mutations predected using pint showed that TCF with PI3K which is a drug target in the original model affect goal reachability. And to see if there is a synergy between them we run MaBoSS stochastic simulation to compare the growth.¶

6.1. TCF single perturbation (MaBoSS simulation).¶

First we run perturbation for TCF alone to identify the growth( output)

In [38]:
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

(Growth) prosurvival (0) - Antisurvival (1) = -1¶

6.2.Double mutation (TCF&PI3K)¶

MaBoSS simulatio used to confirm the predicted results from Pint

In [39]:
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.

Growth = Prosurvival (0) - Antisurvival (1) = -1 Which is the same as the growth of TCF alone (-1) so there is no synergy in double mutation¶

6.3. MEK,PI3K&TCF (Triple mutation).¶

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.

In [21]:
# 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)
Out[21]:
In [22]:
# Trapspaces

traps_pert_MEK_PI3K_TCF = biolqm.trapspace(pert_MEK_PI3K_TCF)
tabulate(traps_pert_MEK_PI3K_TCF)
Out[22]:
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.

In [20]:
# 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.

7. Results and Discussion.¶

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.

8. Conclusion and further analysis.¶

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.

9. Sources.¶

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

In [ ]: