Changes default calculations to symbolic outputs

As an additional test for this branch, I have consolidated all the code from the KDA manuscript and executed it on this branch to ensure it does not break the current code.

KDA Manuscript Code

```python def main_paper_code(): import numpy as np import networkx as nx from kda import graph_utils # define connectivity matrix K = np.array([ [0, 1, 0, 1], [1, 0, 1, 1], [0, 1, 0, 1], [1, 1, 1, 0], ]) # initialize empty graph object G = nx.MultiDiGraph() # populate the edge data using KDA utility graph_utils.generate_edges(G, K) from kda import calculations as calcs prob_strs = calcs.calc_state_probs(G, output_strings=True) p1 = prob_strs[0] print(p1) from sympy import symbols # create symbols for current variables (k12, k21, k23, k32, k34, k43, k41, k14, k24, k42) = symbols("k12, k21, k23, k32, k34, k43, k41, k14, k24, k42") # create symbols for substitution (R_on, R_off, L_on, L_off, k_leak, R_int, R_ext, L_int, L_ext) = symbols("R_on, R_off, L_on, L_off, k_leak, R_int, R_ext, L_int, L_ext") # define variable substitution mapping model = {k12:R_off, k21:R_on*R_int, k23:L_on*L_int, k32:L_off, k34:L_off, k43:L_on*L_ext, k41:R_on*R_ext, k14:R_off, k24:k_leak, k42:k_leak} # perform variable subs, simplify p1 = p1.subs(model).simplify() # generate the expression for the # net cycle flux of cycle a J_a = calcs.calc_net_cycle_flux(G, cycle=[0, 1, 2, 3], order=[0, 1], key="name", output_strings=True) # perform variable substitutions J_a = J_a.subs(model).simplify() # create state probability expressions p1, p2, p3, p4 = calcs.calc_state_probs(G, output_strings=True) # simplify probability expressions p1 = p1.subs(model) p2 = p2.subs(model) p3 = p3.subs(model) # create operational flux expressions J_R = (p1*model[k12] - p2*model[k21]) J_L = (p2*model[k23] - p3*model[k32]) def antiporter_comparison_code(): import numpy as np import networkx as nx from sympy import symbols from kda import graph_utils, calculations as calcs # define connectivity matrix K = np.array( [ [0, 1, 0, 1, 0, 1], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [1, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [1, 0, 0, 0, 1, 0], ] ) # create the 6-state kinetic diagram G_leak = nx.MultiDiGraph() # populate edge data using KDA utility graph_utils.generate_edges(G_leak, K) # create symbols for current variables (k12, k21, k23, k32, k34, k43, k45, k54, k56, k65, k61, k16, k14, k41) = symbols( "k12 k21 k23 k32 k34 k43 k45 k54 k56 k65 k61 k16 k14, k41") # create variables for substitution (H_on, H_off, Na_on, Na_off, k_conf, k_leak, H_in, H_out, c_Na) = symbols( "H_on, H_off, Na_on, Na_off, k_conf, k_leak, H_in, H_out, c_Na") # define variable substitution mapping model = { k12:H_on*H_out, k21:H_off, k23:k_conf, k32:k_conf, k34:H_off, k43:H_on*H_in, k45:Na_on*c_Na, k54:Na_off, k56:k_conf, k65:k_conf, k61:Na_off, k16:Na_on*c_Na, k14:k_leak, k41:k_leak, } # collect the cycles using KDA cycles = graph_utils.find_all_unique_cycles(G_leak) # label cycles according to figure cycle_labels = ["a", "c", "b"] # set the cycle directions (CW) cycle_orders = [[0, 1], [0, 1], [5, 0]] # generate net cycle flux expressions G_leak_cycles = {} for label, _cycle, _order in zip(cycle_labels, cycles, cycle_orders): func = calcs.calc_net_cycle_flux( G_leak, cycle=_cycle, order=_order, key='name', output_strings=True) # perform variable substitutions func = func.subs(model).simplify() G_leak_cycles[label] = {"cycle": _cycle, "order": _order, "func": func} # assign net cycle flux expressions J_a = G_leak_cycles["a"]["func"] J_b = G_leak_cycles["b"]["func"] J_c = G_leak_cycles["c"]["func"] # generate state probability expressions prob_strs = calcs.calc_state_probs( G_leak, key='name', output_strings=True) p1, p2, p3, p4, p5, p6 = prob_strs # calculate operational fluxes using cycles J_H_cycle = (J_a + J_c).simplify() J_Na_cycle = (J_b + J_c).simplify() # calculate operational fluxes using transition fluxes # J_H = p1 * k12 - p2 * k21 = J12 J_H_trans = (p1.subs(model) * model[k12] - p2.subs(model) * model[k21]).simplify() # J_Na = p6 * k61 - p1 * k16 = J61 J_Na_trans = (p6.subs(model) * model[k61] - p1.subs(model) * model[k16]).simplify() # compare operational flux expressions # created using cycles and transitions print((J_H_cycle - J_H_trans).simplify() == 0) print((J_Na_cycle - J_Na_trans).simplify() == 0) def emre_comparison_code(): import sys # raise the recursion limit so SymPy doesn't raise a RecursionError # when parsing the EmrE normalization expression (sigma) sys.setrecursionlimit(5000) import numpy as np import networkx as nx from sympy import symbols from sympy.parsing.sympy_parser import parse_expr from kda import graph_utils, diagrams, calculations as calcs # define connectivity matrix K = np.array( [ [0, 1, 1, 0, 0, 0, 1, 0], [1, 0, 0, 1, 0, 0, 0, 1], [1, 0, 0, 1, 1, 0, 0, 0], [0, 1, 1, 0, 0, 1, 0, 0], [0, 0, 1, 0, 0, 1, 1, 0], [0, 0, 0, 1, 1, 0, 0, 1], [1, 0, 0, 0, 1, 0, 0, 1], [0, 1, 0, 0, 0, 1, 1, 0], ] ) # create the 8-state kinetic diagram G = nx.MultiDiGraph() # populate edge data using KDA utility graph_utils.generate_edges(G, K) # create symbols for current variables (k31, k13, k57, k75, k42, k24, k68, k86, k34, k43, k56, k65, k12, k21, k78, k87, k71, k17, k53, k35, k64, k46, k82, k28) = symbols( "k31 k13 k57 k75 k42 k24 k68 k86 k34 k43 k56 k65 k12 k21 k78 k87 k71 k17 k53 k35 k64 k46 k82 k28") # create variables for substitution (H_on, H_off, D_on, D_off, H_in, H_out, D_in, D_out, k_AA_anti, k_AA_sym) = symbols("H_on H_off D_on D_off H_in H_out D_in D_out k_AA_anti k_AA_sym") # define variable substitution mapping model = { k31: H_on * H_out, k13: H_off, k57: H_on * H_in, k75: H_off, k42: H_on * H_out, k24: H_off, k68: H_on * H_in, k86: H_off, k34: D_on * D_out, k43: D_off, k56: D_on * D_in, k65: D_off, k12: D_on * D_out, k21: D_off, k78: D_on * D_in, k87: D_off, k71: k_AA_anti, k17: k_AA_anti, k53: k_AA_sym, k35: k_AA_sym, k64: k_AA_anti, k46: k_AA_anti, k82: k_AA_sym, k28: k_AA_sym, } # collect the directional edges beforehand dir_edges = diagrams.generate_directional_diagrams(G, return_edges=True) # create the normalization factor for net cycle fluxes sigma = calcs.calc_sigma(G, dir_edges, key="name", output_strings=True) sigma = parse_expr(sigma) sigma = sigma.subs(model).simplify() # define cycle info EmrE_cycles = { 1: {"cycle": [0, 1, 3, 2, 4, 5, 7, 6], "order": [6, 0]}, 2: {"cycle": [0, 6, 7, 5, 4, 2], "order": [6, 0]}, 3: {"cycle": [0, 6, 7, 5, 3, 2], "order": [6, 0]}, 4: {"cycle": [0, 6, 7, 5, 3, 1], "order": [6, 0]}, 6: {"cycle": [0, 6, 7, 1, 3, 2], "order": [6, 0]}, 7: {"cycle": [0, 6, 7, 1], "order": [6, 0]}, 8: {"cycle": [0, 2, 3, 1, 7, 5, 4, 6], "order": [6, 0]}, 9: {"cycle": [0, 6, 4, 5, 7, 1], "order": [6, 0]}, 10: {"cycle": [0, 6, 4, 5, 3, 2], "order": [6, 0]}, 11: {"cycle": [0, 6, 4, 5, 3, 1], "order": [6, 0]}, 13: {"cycle": [0, 6, 4, 2, 3, 1], "order": [6, 0]}, 14: {"cycle": [0, 6, 4, 2], "order": [6, 0]}, 15: {"cycle": [0, 1, 3, 5, 7, 6, 4, 2], "order": [0, 2]}, 16: {"cycle": [0, 2, 4, 6, 7, 1], "order": [0, 2]}, 17: {"cycle": [0, 2, 4, 5, 7, 1], "order": [0, 2]}, 18: {"cycle": [0, 2, 4, 5, 3, 1], "order": [0, 2]}, 19: {"cycle": [0, 2, 3, 5, 7, 1], "order": [0, 2]}, 20: {"cycle": [0, 1, 7, 6, 4, 5, 3, 2], "order": [0, 2]}, 22: {"cycle": [1, 7, 6, 4, 5, 3], "order": [1, 3]}, 23: {"cycle": [1, 7, 6, 4, 2, 3], "order": [2, 4]}, 24: {"cycle": [1, 7, 5, 4, 2, 3], "order": [2, 4]}, 25: {"cycle": [1, 7, 5, 3], "order": [1, 3]}, 26: {"cycle": [2, 3, 5, 7, 6, 4], "order": [2, 4]}, 27: {"cycle": [2, 3, 5, 4], "order": [2, 4]}, } # generate net cycle flux expression numerators for idx, cycle_info in EmrE_cycles.items(): pi_diff = calcs.calc_pi_difference(G, cycle=cycle_info["cycle"], order=cycle_info["order"], key="name", output_strings=True) flux_diags = diagrams.generate_flux_diagrams(G, cycle=cycle_info["cycle"]) sigma_K = calcs.calc_sigma_K(G, cycle=cycle_info["cycle"], flux_diags=flux_diags, key="name", output_strings=True) if sigma_K == 1: func = parse_expr(pi_diff) else: func = parse_expr(pi_diff) * parse_expr(sigma_K) # perform variable substitutions func = func.subs(model).simplify() EmrE_cycles[idx]["func"] = func # define contributing cycles for both ligands H_cycles = [1, 2, 3, 4, 10, 11, 13, 14, 16, 17, 19, 20, 22, 23, 24, 25] D_cycles = [3, 4, 6, 7, 8, 9, 10, 11, 15, 16, 17, 18, 23, 24, 26, 27] J_H_cycle = 0 for idx in H_cycles: J_H_cycle += EmrE_cycles[idx]["func"] J_D_cycle = 0 for idx in D_cycles: # negative contributors if idx in (3, 4, 6, 7, 8, 9, 10, 11): J_D_cycle -= EmrE_cycles[idx]["func"] else: J_D_cycle += EmrE_cycles[idx]["func"] # normalize the operational fluxes (cycles) J_H_cycle = (J_H_cycle / sigma).simplify() J_D_cycle = (J_D_cycle / sigma).simplify() # generate state probability expressions prob_strs = calcs.calc_state_probs(G, key='name', output_strings=True) p1, p2, p3, p4, p5, p6, p7, p8 = prob_strs # calculate operational fluxes using transition fluxes J_13 = p1.subs(model) * model[k13] - p3.subs(model) * model[k31] J_24 = p2.subs(model) * model[k24] - p4.subs(model) * model[k42] J_21 = p2.subs(model) * model[k21] - p1.subs(model) * model[k12] J_43 = p4.subs(model) * model[k43] - p3.subs(model) * model[k34] J_H_trans = (J_13 + J_24).simplify() J_D_trans = (J_21 + J_43).simplify() # compare operational flux expressions # created using cycles and transitions print((J_H_cycle - J_H_trans).simplify() == 0) print((J_D_cycle - J_D_trans).simplify() == 0) if __name__ == "__main__": main_paper_code() antiporter_comparison_code() emre_comparison_code() ```

The above script executes and runs to completion on this branch. However, it actually does not run on the current master due to a missing key="name"... in one of the very first lines of code in the paper... 🤦‍♂️

So I guess this is serving multiple purposes 🦝

Becksteinlab / kda

Changes default calculations to symbolic outputs #88

Description

Changes