from pathlib import Path from pandas import DataFrame from biodivine_aeon import *

original_model = BooleanNetwork.from_sbml(Path('butanol-pathway.sbml').read_text()) model = BooleanNetwork.from_aeon(Path('butanol-pathway-relaxed.aeon').read_text()) # Try to build a symbolic asynchronous graph from the original model. # This should list problems in the model structure. try: graph = SymbolicAsyncGraph(original_model) except Exception as e: print("Cannot create graph from the original model.") print(e) # Now print regulations that changed between the two models. original_rg = original_model.graph() rg = model.graph() print("==== spoIIE -> spoIIAB updated from activation to inhibition ====") print(original_rg.find_regulation("spoIIE", "spoIIAB")) print(rg.find_regulation("spoIIE", "spoIIAB")) print("==== sigA -> spo0A_p observability is removed ====") print(original_rg.find_regulation("sigA", "spo0A_p")) print(rg.find_regulation("sigA", "spo0A_p")) print("==== sporulation -> glucose___PTS observability is removed ====") print(original_rg.find_regulation("sporulation", "glucose___PTS")) print(rg.find_regulation("sporulation", "glucose___PTS")) # And we can succssfully create a symbolic graph from this relaxed model: graph = SymbolicAsyncGraph(model)

attractors = find_attractors(graph) assert(len(attractors) == 2) assert(attractors[0].colors().intersect(attractors[1].colors()).is_empty()) attractor = attractors[0].union(attractors[1])

# A simple helper function to fix the value of a variable in the given set def select_value(states, variable, value): return states.intersect(graph.fix_variable(variable, value)) # Classify the parametrisations (colors) depending on whether a variable can be always true, # always false, or change its value (i.e. the set states contains both instances where the # variable is true and false). def check_stability(states, variable): on_colors = select_value(states, variable, True).colors() off_colors = select_value(states, variable, False).colors() both = on_colors.intersect(off_colors) return { "on": on_colors.minus(both), "off": off_colors.minus(both), "both": both } # Separate the attractor set into three subsets corresponding # to the three modes in the decision tree: # Mode I sigA_off = select_value(attractor, "sigA", False) sigA_on = select_value(attractor, "sigA", True) # Mode II glucose_off = select_value(sigA_on, "glucose", False) # Mode III glucose_on = select_value(sigA_on, "glucose", True) print(check_stability(glucose_off, "butanol")) print(check_stability(glucose_on, "butanol"))

table = [] def process_stability(data): size_all = data["on"].cardinality() + data["off"].cardinality() + data["both"].cardinality() result = {} if not data["on"].is_empty(): result["on"] = round((data["on"].cardinality() / size_all) * 100.0) if not data["off"].is_empty(): result["off"] = round((data["off"].cardinality() / size_all) * 100.0) if not data["both"].is_empty(): result["both"] = round((data["both"].cardinality() / size_all) * 100.0) return result def add_row(table, variable): c1 = process_stability(check_stability(sigA_off, variable)) c2 = process_stability(check_stability(glucose_off, variable)) c3 = process_stability(check_stability(glucose_on, variable)) table.append([variable, c1, c2, c3]) for variable in ["acetic_acid", "butyric_acid", "lactic_acid", "butanol", "ethanol", "acetone", "sporulation"]: add_row(table, variable) DataFrame(table, columns=["variable", "sigA=off (I)", "glucose=off (II)", "glucose=on (III)"])

gene_table = [] for variable in ["ack", "adc", "bdhAB", "buk1", "crt", "hbd", "pta"]: add_row(gene_table, variable) DataFrame(gene_table, columns=["variable", "sigA=off (I)", "glucose=off (II)", "glucose=on (III)"])

model = BooleanNetwork.from_sbml(Path('butanol-pathway.sbml').read_text()) # Declare uninterpreted functions: model.add_parameter({ "name": "f_1", "arity": 2 }) model.add_parameter({ "name": "f_2", "arity": 2 }) model.add_parameter({ "name": "f_3", "arity": 3 }) # Change original functions and confirm that the change has been successful: f_spo0A_p = "spo0A & phosphorylation & f_1(sigA, sporulation)" f_glucose_PTS = "glucose & PTS & f_2(cell_membrane, sporulation)" f_spoIIAB = "spo0A_p & f_3(sigH, spoIIAA, spoIIE)" model.set_update_function("spo0A_p", f_spo0A_p) model.set_update_function("glucose___PTS", f_glucose_PTS) model.set_update_function("spoIIAB", f_spoIIAB) model.get_update_function("spo0A_p")

Path('butanol-pathway-f1-f2-f3.aeon').write_text(model.to_aeon())

graph = SymbolicAsyncGraph(model) attractors = find_attractors(graph) assert(len(attractors) == 2) assert(attractors[0].colors().intersect(attractors[1].colors()).is_empty()) attractor = attractors[0].union(attractors[1])

# A simple helper function to fix the value of a variable in the given set def select_value(states, variable, value): return states.intersect(graph.fix_variable(variable, value)) # Classify the parametrisations (colors) depending on whether a variable can be always true, # always false, or change its value (i.e. the set states contains both instances where the # variable is true and false). def check_stability(states, variable): on_colors = select_value(states, variable, True).colors() off_colors = select_value(states, variable, False).colors() both = on_colors.intersect(off_colors) return { "on": on_colors.minus(both), "off": off_colors.minus(both), "both": both } # Separate the attractor set into three subsets corresponding # to the three modes in the decision tree: # Mode I sigA_off = select_value(attractor, "sigA", False) sigA_on = select_value(attractor, "sigA", True) # Mode II glucose_off = select_value(sigA_on, "glucose", False) # Mode III glucose_on = select_value(sigA_on, "glucose", True)

table = [] def process_stability(data): size_all = data["on"].cardinality() + data["off"].cardinality() + data["both"].cardinality() result = {} if not data["on"].is_empty(): result["on"] = round((data["on"].cardinality() / size_all) * 100.0) if not data["off"].is_empty(): result["off"] = round((data["off"].cardinality() / size_all) * 100.0) if not data["both"].is_empty(): result["both"] = round((data["both"].cardinality() / size_all) * 100.0) return result def add_row(table, variable): c1 = process_stability(check_stability(sigA_off, variable)) c2 = process_stability(check_stability(glucose_off, variable)) c3 = process_stability(check_stability(glucose_on, variable)) table.append([variable, c1, c2, c3]) for variable in ["acetic_acid", "butyric_acid", "lactic_acid", "butanol", "ethanol", "acetone", "sporulation"]: add_row(table, variable) DataFrame(table, columns=["variable", "sigA=off", "glucose=off", "glucose=on"])

bdd_variables = graph.bdd_variables() bdd = check_stability(glucose_on, "acetic_acid")["off"].to_bdd() # We leave the BDD visualisation commented out because the graph is large # graphviz.Source(bdd.to_dot(bdd_variables))

# Reload model again model = BooleanNetwork.from_sbml(Path('butanol-pathway.sbml').read_text()) # Do not add f_1 any more model.add_parameter({ "name": "f_2", "arity": 2 }) model.add_parameter({ "name": "f_3", "arity": 3 }) # Re-declare update functions, using the new fixed f_1 f_spo0A_p = "spo0A & phosphorylation & sigA & !sporulation" f_glucose_PTS = "glucose & PTS & f_2(cell_membrane, sporulation)" f_spoIIAB = "spo0A_p & f_3(sigH, spoIIAA, spoIIE)" model.set_update_function("spo0A_p", f_spo0A_p) model.set_update_function("glucose___PTS", f_glucose_PTS) model.set_update_function("spoIIAB", f_spoIIAB) model.get_update_function("spo0A_p")

Path('butanol-pathway-f2-f3.aeon').write_text(model.to_aeon())

graph = SymbolicAsyncGraph(model) attractors = find_attractors(graph) assert(len(attractors) == 2) assert(attractors[0].colors().intersect(attractors[1].colors()).is_empty()) attractor = attractors[0].union(attractors[1])

# Recompute the sets corresponding to different modes: sigA_off = select_value(attractor, "sigA", False) sigA_on = select_value(attractor, "sigA", True) glucose_off = select_value(sigA_on, "glucose", False) glucose_on = select_value(sigA_on, "glucose", True) table = [] for variable in ["acetic_acid", "butyric_acid", "lactic_acid", "butanol", "ethanol", "acetone", "sporulation"]: add_row(table, variable) DataFrame(table, columns=["variable", "sigA=off (I)", "glucose=off (II)", "glucose=on (III)"])

gene_table = [] for variable in ["ack", "adc", "bdhAB", "buk1", "crt", "hbd", "pta"]: add_row(gene_table, variable) DataFrame(gene_table, columns=["variable", "sigA=off", "glucose=off", "glucose=on"])

bdd_variables = graph.bdd_variables() def check_attractor_size(table_row, attractor): bdd_param_true = bdd_variables.mk_literal(row, True) attractor_true = attractor.to_bdd().and_not(bdd_variables.mk_literal(row, False)) attractor_false = attractor.to_bdd().and_not(bdd_variables.mk_literal(row, True)) true_vertices = graph.unit_colored_vertices().copy_with(attractor_true).vertices() false_vertices = graph.unit_colored_vertices().copy_with(attractor_false).vertices() only_true = true_vertices.minus(false_vertices) only_false = false_vertices.minus(true_vertices) all_card = only_true.cardinality() + only_false.cardinality() if not all_card == 0.0: true_percent = 100.0 * only_true.cardinality() / all_card false_percent = 100.0 * only_false.cardinality() / all_card impact = 100.0 * all_card / attractor.vertices().cardinality() print(table_row, { "true": true_percent, "false": false_percent, "impact": impact }) else: print(table_row, " - no influence") for row in ["f_2[0,0]", "f_2[0,1]", "f_2[1,0]", "f_2[1,1]"]: check_attractor_size(row, glucose_on) print("=====") for row in ["f_3[0,0,0]", "f_3[0,0,1]", "f_3[0,1,0]", "f_3[0,1,1]", "f_3[1,0,0]", "f_3[1,0,1]", "f_3[1,1,0]", "f_3[1,1,1]"]: check_attractor_size(row, glucose_on)

# Reload model again model = BooleanNetwork.from_sbml(Path('butanol-pathway.sbml').read_text()) # Fix update functions using the repaired values f_spo0A_p = "spo0A & phosphorylation & sigA & !sporulation" f_glucose_PTS = "glucose & PTS & cell_membrane & !sporulation" f_spoIIAB = "spo0A_p & sigH & (spoIIE | !spoIIAA)" model.set_update_function("spo0A_p", f_spo0A_p) model.set_update_function("glucose___PTS", f_glucose_PTS) model.set_update_function("spoIIAB", f_spoIIAB) Path('butanol-pathway-repaired.sbml').write_text(model.to_sbml())