Acyclic boolean models of cell signaling (experimental)#

This guide shows the basics of modeling intracellular signaling using acyclic boolean models for multi-perturbation experiments. For this, we will use an extension of the CellNopt ILP method implemented with CORNETO.

  • Author: Attila Gabor (Heidelberg University, SaezLab)

  • Reviewers: Pablo Rodriguez Mier (Heidelberg University, SaezLab)

NOTE: This notebook is still experimental.

import corneto as cn
import corneto.methods.signal.cellnopt_ilp as cno

cn.info()
Installed version:v1.0.0.dev3 (latest stable: v1.0.0-alpha)
Available backends:CVXPY v1.6.0, PICOS v2.5.1
Default backend (corneto.opt):CVXPY
Installed solvers:CVXOPT, GLPK, GLPK_MI, HIGHS, SCIP, SCIPY
Graphviz version:v0.20.3
Installed path:/home/runner/work/corneto/corneto/corneto
Repository:https://github.com/saezlab/corneto

Introduction to a Boolean model#

As a first step we define start with a small network and corresponding datasets just to illustrate the concepts. The prior knowledge network of ‘G1’ contains three nodes, EGF, TNFa and Ras. Edges show the potential activation of the RAS protein.

# Definition of the prior knowledge network
G1 = cn.Graph.from_sif_tuples(
    [
        ("EGF", 1, "AND1"),
        ("TNFa", 1, "AND1"),
        ("AND1", 1, "Ras"),
        ("EGF", 1, "Ras"),
        ("TNFa", 1, "Ras"),
    ]
)
G1.plot()
../../_images/4e46465ca2e4f77c764f290c442ebba66ba6d1e724855d58478655f7aa98c482.svg

Now, we organize some in-silico experimental data into a format that can be used by the CNO toolbox.

First, we define 4 conditions (exp0 … exp3). In each condition, the activated nodes (EGF or TNFa) are set to 1 and the non-activated nodes to 0. We also define the values of the corresponding outputs (Ras).

Note that the output, Ras is only active (1.0) in the 4th condition, when both of the inputs are also activated. This corresponds to an AND relationship: RAS is active only if both EGF and TNFa are activated. We expect that the optimization finds the correct subgraph of the above network, which contains only the 3 nodes with the AND gates, while individual edges from EGF and TNFa to RAS are removed.

# RAS is only active iff both EGF and TNFa are active -> we need to identify the AND gate
exp_list_G1_and = {
    "exp0": {"input": {"EGF": 0, "TNFa": 0}, "output": {"Ras": 0}},
    "exp1": {"input": {"EGF": 1, "TNFa": 0}, "output": {"Ras": 0}},
    "exp2": {"input": {"EGF": 0, "TNFa": 1}, "output": {"Ras": 0}},
    "exp3": {"input": {"EGF": 1, "TNFa": 1}, "output": {"Ras": 1}},
}
G1_prep = cno.expand_graph_for_flows(G1, exp_list_G1_and)
G1_prep.plot()
../../_images/b79eb8beaa21c67ef98dc11e8713b8280d5e7def532eaad4e12009a795ed48ed.svg
cno.plot_data(exp_list_G1_and)
../../_images/be9d2266e948081639e284b6f0a37d084b530339435290e40d10f028978692fe.png
P = cno.cellnoptILP(G1_prep, exp_list_G1_and, verbose=False)
cno.plot_fitness(G1, exp_list_G1_and, P, measured_only=True)
../../_images/e1484e218055cd7792d31c0d7f3283ee501311430190dea602f356a70c32e896.png
G1_prep.plot(
    custom_edge_attr=cno.cno_style(P, flow_name="edge_activates", scale=None, iexp=3)
)
../../_images/c61d3607bda18fa3e977dc4136911ddbee35e8fcef97a32cc8634abb595da2cd.svg

Tasks:

  1. modify the experiments so that it corresponds to a scenario where only EGF could activate RAS, but not TNFa.

  2. modify the experiments so that it corresponds to a scenario where both TNFa or EGF could activate RAS.

Inhibition: !EGF activates RAS#

Let’s define the same network but now we will add an inhibition of EGF on RAS. This means that if EGF is active, RAS should be inactive. Note the difference in the definition of the prior knowledge network.

We also change the experiments so that the output RAS is active only when TNFa is active.

# Test graph
G2 = cn.Graph.from_sif_tuples(
    [
        ("EGF", -1, "AND1"),
        ("TNFa", 1, "AND1"),
        ("AND1", 1, "Ras"),
        ("EGF", -1, "Ras"),
        ("TNFa", 1, "Ras"),
    ]
)

# RAS is only active iff both EGF and TNFa are active -> we need to idetify the AND gate
exp_list_G2_egf = {
    "exp0": {"input": {"EGF": 0, "TNFa": 0}, "output": {"Ras": 1}},
    "exp1": {"input": {"EGF": 1, "TNFa": 0}, "output": {"Ras": 0}},
    "exp2": {"input": {"EGF": 0, "TNFa": 1}, "output": {"Ras": 1}},
    "exp3": {"input": {"EGF": 1, "TNFa": 1}, "output": {"Ras": 0}},
}
G2_prep = cno.expand_graph_for_flows(G2, exp_list_G2_egf)
G2_prep.plot()
../../_images/015ceeaeea34d4052d32e66820424fe5446084f93d597ca8cc81e22fb11a9e41.svg
P2 = cno.cellnoptILP(G2_prep, exp_list_G2_egf, verbose=False)
cno.plot_fitness(G2_prep, exp_list_G2_egf, P2, measured_only=True)
../../_images/75a13bb968797193660f9eb77cede1d05c6f26a06227507e2a8f04565845e4bc.png
G2_prep.plot(
    custom_edge_attr=cno.cno_style(P2, flow_name="edge_activates", scale=None, iexp=0)
)
../../_images/4ce7860609960047e739d838509ef98935a8904cd8f2aa979eaab392052342cf.svg

Inhibition with AND gates#

Please remember, inhibitions are added to the incoming edge of the AND node. The outgoing edges of the AND node are always activating by convention.

# Test graph
G2 = cn.Graph.from_sif_tuples(
    [
        ("EGF", -1, "AND1"),
        ("TNFa", 1, "AND1"),
        ("AND1", 1, "Ras"),
        ("EGF", -1, "Ras"),
        ("TNFa", 1, "Ras"),
    ]
)

# RAS is only active iff both EGF and TNFa are active -> we need to idetify the AND gate
exp_list_G2_and = {
    "exp0": {"input": {"EGF": 0, "TNFa": 0}, "output": {"Ras": 0}},
    "exp1": {"input": {"EGF": 1, "TNFa": 0}, "output": {"Ras": 0}},
    "exp2": {"input": {"EGF": 0, "TNFa": 1}, "output": {"Ras": 1}},
    "exp3": {"input": {"EGF": 1, "TNFa": 1}, "output": {"Ras": 0}},
}
G2_prep = cno.expand_graph_for_flows(G2, exp_list_G2_and)
G2_prep.plot()
../../_images/015ceeaeea34d4052d32e66820424fe5446084f93d597ca8cc81e22fb11a9e41.svg
P2 = cno.cellnoptILP(G2_prep, exp_list_G2_and, verbose=False)
G2_prep.plot(
    custom_edge_attr=cno.cno_style(P2, flow_name="edge_activates", scale=None, iexp=2)
)
../../_images/2ed3796501afc91e03c6e1583309ce55fc40f84d6e4b51e9ee525c1413413a7b.svg
cno.plot_fitness(G2_prep, exp_list_G2_and, P2, measured_only=True)
../../_images/1c87fef36311dd592ce07e7c7534ab6df2b353208392db2528da4df29fc898dc.png

Complex case study with inhibition#

Let’s look at a more realistic case study which contains a larger network and more experiments.

Also note that the we define inhibition in the experimental description. This could be, for example, the effect of small molecule kinase inhibitor.

G3 = cn.Graph.from_sif_tuples(
    [
        ("EGF", 1, "Ras"),
        ("EGF", 1, "PI3K"),
        ("TNFa", 1, "PI3K"),
        ("TNFa", 1, "TRAF6"),
        ("TRAF6", 1, "p38"),
        ("TRAF6", 1, "Jnk"),
        ("TRAF6", 1, "NFkB"),
        ("Jnk", 1, "cJun"),
        ("p38", 1, "Hsp27"),
        ("PI3K", 1, "Akt"),
        ("Ras", 1, "Raf"),
        ("Raf", 1, "Mek"),
        ("Akt", -1, "Mek"),
        ("Mek", 1, "p90RSK"),
        ("Mek", 1, "Erk"),
        ("Erk", 1, "Hsp27"),
    ]
)
# G3.add_edge((), 'EGF')
# G3.add_edge((), 'TNFa')
# G3.add_edge('Akt', ())
# G3.add_edge('cJun', ())
# G3.add_edge('Hsp27', ())
# G3.add_edge('NFkB', ())
# G3.add_edge('Erk', ())
# G3.add_edge('p90RSK', ())
# G3.add_edge('Jnk', ())

exp_list_toy_full = {
    "exp0": {
        "input": {"EGF": 0, "TNFa": 0},
        "inhibition": {},
        "output": {
            "Akt": 0,
            "Hsp27": 0,
            "NFkB": 0,
            "Erk": 0,
            "p90RSK": 0,
            "Jnk": 0,
            "cJun": 0,
        },
    },
    "exp1": {
        "input": {"EGF": 1, "TNFa": 0},
        "inhibition": {},
        "output": {
            "Akt": 0.91,
            "Hsp27": 0,
            "NFkB": 0.86,
            "Erk": 0.8,
            "p90RSK": 0.88,
            "Jnk": 0,
            "cJun": 0,
        },
    },
    "exp2": {
        "input": {"EGF": 0, "TNFa": 1},
        "inhibition": {},
        "output": {
            "Akt": 0.82,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.0,
            "p90RSK": 0.0,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
    "exp3": {
        "input": {"EGF": 1, "TNFa": 1},
        "inhibition": {},
        "output": {
            "Akt": 0.91,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.8,
            "p90RSK": 0.88,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
    "exp4": {
        "input": {"EGF": 1, "TNFa": 0},
        "inhibition": {"Raf": 1},
        "output": {
            "Akt": 0.91,
            "Hsp27": 0,
            "NFkB": 0.86,
            "Erk": 0.0,
            "p90RSK": 0.0,
            "Jnk": 0,
            "cJun": 0,
        },
    },
    "exp5": {
        "input": {"EGF": 0, "TNFa": 1},
        "inhibition": {"Raf": 1},
        "output": {
            "Akt": 0.82,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.0,
            "p90RSK": 0.0,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
    "exp6": {
        "input": {"EGF": 1, "TNFa": 1},
        "inhibition": {"Raf": 1},
        "output": {
            "Akt": 0.91,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.0,
            "p90RSK": 0.0,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
    "exp7": {
        "input": {"EGF": 1, "TNFa": 0},
        "inhibition": {"PI3K": 1},
        "output": {
            "Akt": 0.0,
            "Hsp27": 0,
            "NFkB": 0.0,
            "Erk": 0.8,
            "p90RSK": 0.88,
            "Jnk": 0,
            "cJun": 0,
        },
    },
    "exp8": {
        "input": {"EGF": 0, "TNFa": 1},
        "inhibition": {"PI3K": 1},
        "output": {
            "Akt": 0.0,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.0,
            "p90RSK": 0.0,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
    "exp9": {
        "input": {"EGF": 1, "TNFa": 1},
        "inhibition": {"PI3K": 1},
        "output": {
            "Akt": 0.0,
            "Hsp27": 0.7,
            "NFkB": 0.90,
            "Erk": 0.8,
            "p90RSK": 0.88,
            "Jnk": 0.25,
            "cJun": 0.4,
        },
    },
}


G3.plot()
../../_images/a7b2152f0be063ecdd4d073baa74125df189602839b0ffcfd773741e4eb18dcf.svg
cno.plot_data(exp_list_toy_full)
../../_images/b43ae7053011795c9a529bcee260656ec412695d580e057931338a775733af2d.png
G3_prep = cno.expand_graph_for_flows(G3, exp_list_toy_full)
G3_prep.plot()
../../_images/7b5cab6502738b247db3cad4e00793f65f28273b9b9fd7746dc69057850c2f9b.svg
P_G3 = cno.cellnoptILP(G3_prep, exp_list_toy_full, verbose=False)
G3_prep.plot(
    custom_edge_attr=cno.cno_style(P_G3, flow_name="edge_activates", scale=None, iexp=3)
)
../../_images/e2f5b6047330ace76c8a3dd1ae1c5d854857d9f7c898c8f9e6b56c4647fbaa83.svg
cno.plot_fitness(G3_prep, exp_list_toy_full, P_G3)
../../_images/c95638a37e5f9424ee76f1a3d530d5207b9ed22dc81d01dc62952dbdea95403d.png