Source code for crnt4sbml.c_graph

import networkx
import numpy
import sympy
import scipy.optimize
import sys
import matplotlib.cbook
import warnings
from sys import platform as sys_pf
if sys_pf == 'darwin':
    import matplotlib
    matplotlib.use("TkAgg")
    import matplotlib.pyplot as plt
else:
    import matplotlib.pyplot as plt

warnings.filterwarnings("ignore", category=matplotlib.cbook.mplDeprecation)


[docs]class Cgraph: """ Class for constructing core CRNT values and C-graph of the network. """
[docs] def __init__(self, model): """ Initialization of Cgraph class. See also --------- crnt4sbml.CRNT.get_c_graph() """ self.__g = networkx.DiGraph() # list of nodes and edges in the order they were added self.__g_nodes = [] self.__g_edges = [] self.__dic_species_id_bcs = {} for i in model.getListOfSpecies(): self.__dic_species_id_bcs.update({i.getId(): i.getBoundaryCondition()}) self.__species = [] self.__parse_reactions(model) self.__assign_reaction_type() # self.__species = [i.getId() for i in model.getListOfSpecies() if not i.getBoundaryCondition()] # temporary makes testing of higher deficiency theory easier ###################################################################################### #self.__species_name = [i.getName() for i in model.getListOfSpecies()] #self.__dic_species_id_name = {} #for i in range(len(self.__species)): # self.__dic_species_id_name.update({self.__species[i]: self.__species_name[i]}) ###################################################################################### self.__linkage_classes = [self.__g.subgraph(c) for c in networkx.weakly_connected_components(self.__g)] self.__complexes = [self.__g.nodes[n]['label'] for n in self.__g_nodes] self.__reactions = [self.__g.edges[e]['label'] for e in self.__g_edges] # self.__check_weak_reversibility_of_linkage_classes() # todo: candidate for deletion # core calculations self.__create_y_matrix() self.__create_s_matrix() self.__create_a_matrix() self.__create_mass_action_monomials() if all([i <= 1 for i in self.get_number_of_terminal_strong_lc_per_lc()]) and \ self.get_dim_equilibrium_manifold() > 0: self.__create_BT_matrix() self.__create_lambda_matrix()
# self.__calculate_deficiency() # todo: candidate for deletion # self.__calculate_number_of_conservation_relations() # todo: candidate for deletion # self.__classify_dimensionality() # todo: candidate for deletion def __parse_reactants_and_products(self, reaction): # extracting reactants reactants = [i.getSpecies() for i in reaction.getListOfReactants()] bound_cond = any([self.__dic_species_id_bcs.get(n, n) for n in reactants]) if bound_cond: color = 'red' else: color = 'green' stoichiometries = [int(i.getStoichiometry()) for i in reaction.getListOfReactants()] temp = zip([repr(i) if i != 1 else '' for i in stoichiometries], reactants) reactant_complex_name = '+'.join(['*'.join(filter(lambda x: x != '', i)) for i in temp]) if self.__g.number_of_nodes() == 0 or reactant_complex_name not in [self.__g.nodes[v]['label'] for v in self.__g.nodes]: self.__g.add_node(reactant_complex_name, label=reactant_complex_name, species=reactants, stoichiometries=stoichiometries, species_bc=bound_cond, color=color) self.__g_nodes.append(reactant_complex_name) # extracting products products = [i.getSpecies() for i in reaction.getListOfProducts()] bound_cond = any([self.__dic_species_id_bcs.get(n, n) for n in products]) if bound_cond: color = 'red' else: color = 'green' stoichiometries = [int(i.getStoichiometry()) for i in reaction.getListOfProducts()] temp = zip([repr(i) if i != 1 else '' for i in stoichiometries], products) product_complex_name = '+'.join(['*'.join(filter(lambda x: x != '', i)) for i in temp]) if product_complex_name not in [self.__g.nodes[v]['label'] for v in self.__g.nodes]: self.__g.add_node(product_complex_name, label=product_complex_name, species=products, stoichiometries=stoichiometries, species_bc=bound_cond, color=color) self.__g_nodes.append(product_complex_name) [self.__species.append(i) for i in reactants if (i not in self.__species) and (not self.__dic_species_id_bcs[i])] [self.__species.append(i) for i in products if (i not in self.__species) and (not self.__dic_species_id_bcs[i])] return [reactant_complex_name, product_complex_name] def __parse_reactants_products_and_modifiers(self, reaction): modifiers = [i.getSpecies() for i in reaction.getListOfModifiers()] # extracting reactants reactants = [i.getSpecies() for i in reaction.getListOfReactants()] reactant_complex_name = self.__create_reactant_or_product_complex(reaction, reactants, modifiers, "reactant") # extracting products products = [i.getSpecies() for i in reaction.getListOfProducts()] product_complex_name = self.__create_reactant_or_product_complex(reaction, products, modifiers, "product") # extracting modifiers bound_cond = any([self.__dic_species_id_bcs.get(n, n) for n in modifiers]) if bound_cond: color = 'red' else: color = 'green' stoichiometries = [int(1.0)] modifiers_complex_name = "".join(modifiers) + "".join(reactants) if modifiers_complex_name not in [self.__g.nodes[v]['label'] for v in self.__g.nodes]: self.__g.add_node(modifiers_complex_name, label=modifiers_complex_name, species=[modifiers_complex_name], stoichiometries=stoichiometries, species_bc=bound_cond, color=color) self.__g_nodes.append(modifiers_complex_name) [self.__species.append(i) for i in reactants if (i not in self.__species) and (not self.__dic_species_id_bcs[i])] [self.__species.append(i) for i in products if (i not in self.__species) and (not self.__dic_species_id_bcs[i])] [self.__species.append(i) for i in modifiers if (i not in self.__species) and (not self.__dic_species_id_bcs[i])] [self.__species.append(i) for i in [modifiers_complex_name] if (i not in self.__species)] return [reactant_complex_name, product_complex_name, modifiers_complex_name] def __create_reactant_or_product_complex(self, reaction, non_modifier, modifiers, comp): bound_cond = any([self.__dic_species_id_bcs.get(n, n) for n in non_modifier + modifiers]) if bound_cond: color = 'red' else: color = 'green' # this statement assumes one species in modifier indices = [i for i in range(len(non_modifier)) if modifiers[0] == non_modifier[i]] if comp == 'reactant': reactants_stoich = [int(i.getStoichiometry()) for i in reaction.getListOfReactants()] # assumes only one modifier modifiers_sotich = int(1.0) if indices: stoich_sum = sum([reactants_stoich[i] for i in indices]) + modifiers_sotich stoichiometries = [reactants_stoich[i] for i in range(len(reactants_stoich)) if i not in indices] + \ [stoich_sum] species_list = [non_modifier[i] for i in range(len(non_modifier)) if i not in indices] + modifiers else: stoichiometries = reactants_stoich + [modifiers_sotich] species_list = non_modifier + modifiers else: products_stoich = [int(i.getStoichiometry()) for i in reaction.getListOfProducts()] # assumes only one modifier modifiers_sotich = int(1.0) if indices: stoich_sum = sum([products_stoich[i] for i in indices]) + modifiers_sotich stoichiometries = [products_stoich[i] for i in range(len(products_stoich)) if i not in indices] + \ [stoich_sum] species_list = [non_modifier[i] for i in range(len(non_modifier)) if i not in indices] + modifiers else: stoichiometries = products_stoich + [modifiers_sotich] species_list = non_modifier + modifiers temp = zip([repr(i) if i != 1 else '' for i in stoichiometries], species_list) complex_name = '+'.join(['*'.join(filter(lambda x: x != '', i)) for i in temp]) if self.__g.number_of_nodes() == 0 or complex_name not in [self.__g.nodes[v]['label'] for v in self.__g.nodes]: self.__g.add_node(complex_name, label=complex_name, species=species_list, stoichiometries=stoichiometries, species_bc=bound_cond, color=color) self.__g_nodes.append(complex_name) return complex_name def __extract_direct_reaction(self, reaction): r, p = self.__parse_reactants_and_products(reaction) if (r, p) not in [i for i in self.__g.edges]: self.__g.add_edge(r, p, label=reaction.getId(), k=None, sbml_label=reaction.getId(), type=None) self.__g_edges.append((r, p)) def __extract_reverse_reaction(self, reaction): r, p = self.__parse_reactants_and_products(reaction) if (p, r) not in [i for i in self.__g.edges]: self.__g.add_edge(p, r, label=reaction.getId() + 'r', k=None, sbml_label=reaction.getId(), type=None) self.__g_edges.append((p, r)) def __extract_catalysis_reaction(self, reaction): modifier = self.__get_celldesigner_reaction_type(reaction) if modifier != 'CATALYSIS': exception_statement = "For reaction " + reaction.getId() + " there is a modifier that cannot be " \ "identified, there are multiple modifiers," \ " or is not catalysis. \nThis may be" \ " a result of not constructing the SBML file " \ "using CellDesigner. \nPlease make the appropriate" \ " changes to this reaction to continue." raise Exception(exception_statement) r, p, m = self.__parse_reactants_products_and_modifiers(reaction) # add complex formation if (r, m) not in [i for i in self.__g.edges]: self.__g.add_edge(r, m, label=reaction.getId() + 'f', k=None, sbml_label=reaction.getId(), type="complex formation") self.__g_edges.append((r, m)) # add complex dissociation if (m, r) not in [i for i in self.__g.edges]: self.__g.add_edge(m, r, label=reaction.getId() + 'd', k=None, sbml_label=reaction.getId(), type="complex dissociation") self.__g_edges.append((m, r)) # add catalysis if (m, p) not in [i for i in self.__g.edges]: self.__g.add_edge(m, p, label=reaction.getId() + 'c', k=None, sbml_label=reaction.getId(), type="catalysis") self.__g_edges.append((m, p)) def __parse_reactions(self, model): for i in model.getListOfReactions(): if i.getNumModifiers() != 0: self.__extract_catalysis_reaction(i) else: self.__extract_direct_reaction(i) if i.getReversible(): self.__extract_reverse_reaction(i) def __get_celldesigner_reaction_type(self, i): i_initial = i.getAnnotation().getChild(0) child_number = self.__get_child_number(i_initial, 'reactionType') if child_number is not None: if i.getNumModifiers() == 1: child_number1 = self.__get_child_number(i_initial, 'listOfModification') if child_number1 is not None: child_number2 = self.__get_child_number(i_initial.getChild(child_number1), 'modification') ii = i_initial.getChild(child_number1).getChild(child_number2) modification = self.__get_modification_type(ii) return modification else: return None else: return None else: return None def __get_child_number(self, i, string_value): num_children = i.getNumChildren() children_strings = [i.getChild(ii).getName() for ii in range(num_children)] if string_value in children_strings: child_number = children_strings.index(string_value) else: child_number = None return child_number def __get_modification_type(self, ii): attribute_names = [ii.getAttrName(iii) for iii in range(ii.getAttributesLength())] if 'type' in attribute_names: indx_type = attribute_names.index('type') attribute_values = [ii.getAttrValue(iii) for iii in attribute_names] return attribute_values[indx_type] else: return None def __assign_reaction_type(self): for i in self.get_g_edges(): current_edge = self.__g.edges[i] if current_edge['type'] is None: reactant_node = self.__g.nodes[i[0]] product_node = self.__g.nodes[i[1]] if sum(reactant_node['stoichiometries']) == 2 and sum(product_node['stoichiometries']) == 1: current_edge['type'] = "complex formation" elif sum(reactant_node['stoichiometries']) == 1 and sum(product_node['stoichiometries']) == 2: if (i[1], i[0]) in self.get_g_edges(): current_edge['type'] = "complex dissociation" else: current_edge['type'] = "catalysis" elif sum(reactant_node['stoichiometries']) == 1 and sum(product_node['stoichiometries']) == 1: current_edge['type'] = "catalysis" elif sum(reactant_node['stoichiometries']) == 1 and sum(product_node['stoichiometries']) == 3: current_edge['type'] = "catalysis" else: output_statement = "The reaction type of reaction " + current_edge['label'] + " could not be " +\ "identified as it does not fit any biological criteria established. \n" +\ "Automatic creation of bounds for optimization routines cannot be completed. \n" print(output_statement) # core CRNT calcs def __create_y_matrix(self): self.__Y = sympy.zeros(len(self.__species), len(self.__complexes)) count = 0 for i in self.__g_nodes: v = self.__g.nodes[i] stoich = dict(zip(v['species'], v['stoichiometries'])) self.__Y[:, count] = [stoich[i] if i in stoich.keys() else 0 for i in self.__species] count += 1 def __create_s_matrix(self): self.__S = sympy.zeros(len(self.__species), len(self.__reactions)) count = 0 for i in self.__g_edges: target = self.__complexes.index(i[1]) source = self.__complexes.index(i[0]) self.__S[:, count] = self.__Y[:, target] - self.__Y[:, source] count += 1 def __create_lambda_matrix(self): self.__Lambda = sympy.zeros(len(self.__complexes), len(self.__linkage_classes)) for li in range(len(self.__linkage_classes)): l_complexes = [i for i in self.__linkage_classes[li].nodes] self.__Lambda[:, li] = [1 if i in l_complexes else 0 for i in self.__complexes] def __create_a_matrix(self): self.__A = sympy.zeros(len(self.__complexes), len(self.__complexes)) inc_list = [list(self.__g.out_edges(n, data='label')) for n in self.__g_nodes] for i in range(len(self.__complexes)): for ii in inc_list[i]: self.__A[i, i] = self.__A[i, i] - sympy.Symbol(ii[2], positive=True) for i in range(len(self.__complexes)): for j in range(len(self.__complexes)): if self.__g.has_edge(self.__complexes[j], self.__complexes[i]): self.__A[i, j] = self.__A[i, j] + sympy.Symbol(self.__g.get_edge_data(self.__complexes[j], self.__complexes[i])['label'], positive=True) def __create_BT_matrix(self): # (bt = NullSpace[Transpose[Y.A]]) // MatrixForm the_null_space = (self.__Y * self.__A).T.nullspace() # getting the number of columns given in nullspace computation sizes = len(the_null_space) # simplifying the entries given in the columns # produced by the nullspace for i in range(0, sizes): the_null_space[i] = sympy.simplify(the_null_space[i]) bt_temp = self.__create_non_negative_b_matrix(the_null_space, sizes) # taking out zero rows if they are in unique vectors bt = sympy.Matrix(bt_temp[~numpy.all(bt_temp == 0, axis=1)]).evalf(10) # making the dimension of bt be lambda by N, just incase # we have more vectors than we should bt = bt[0:self.get_dim_equilibrium_manifold(), :] # putting in a check to make sure that the rows # of bt are linearly independent _, vals = bt.T.rref() if len(vals) != self.get_dim_equilibrium_manifold(): print(" ") print("Error: Not all rows of B^T are linearly independent!") print(" ") sys.exit() self.__B = bt def __simplify_nullspace(self, the_null_space): # Flag that says a species or reaction was found in one of the column vectors of the nullspace denom_flag = False print("beginning null space") #print(the_null_space) print("") for i in range(0, len(the_null_space)): #print(the_null_space[i]) the_null_space[i] = sympy.simplify(the_null_space[i]) print(the_null_space) # find those columns that have a common denominator denom_list = [] for ii in the_null_space: # getting common denominator denom_in_column = [sympy.denom(i) for i in ii] denom_list.append(denom_in_column) print(f"denom_list {denom_list}") like_columns = [] for i in denom_list: same_denoms = [ii for ii in range(len(denom_list)) if denom_list[ii] == i] if same_denoms not in like_columns: like_columns.append(same_denoms) print(like_columns) temp_null = [sympy.zeros(len(self.__species), 1)] * (len(self.__species) - self.__S.rank()) count = 0 for ii in like_columns: for i in ii: temp_null[count] += the_null_space[i] temp_null[count] = sympy.simplify(temp_null[count]) count += 1 print(temp_null) sys.exit() # collecting nullspace vectors with common denominators all_denom = list(set(denom_list)) temp_null = [sympy.zeros(len(self.__species), 1)] * (len(self.__species) - self.__S.rank()) try: for ii in range(len(all_denom)): indices = [i for i, x in enumerate(denom_list) if x == all_denom[ii]] print(f"indices {indices}") for i in indices: temp_null[ii] += the_null_space[i] temp_null[ii] = sympy.simplify(temp_null[ii]) except Exception as e: print("hi") print(f"final temp null {temp_null}") sys.exit() # checking to see if a reaction or species is in the reduced nullspace for i in temp_null: atoms_of_col_vec = [ii.atoms() for ii in i[:, 0]] for ii in atoms_of_col_vec: if any([iii in self.__reactions + self.__species for iii in ii]): denom_flag = True if denom_flag: raise Exception("Nullspace calculation from S contains SymPy variables and this could not be resolved.") return temp_null def __create_non_negative_b_matrix(self, the_null_space, sizes): ''' Idea inspired by https://r.789695.n4.nabble.com/convex-nonnegative-basis-vectors-in-nullspace-of-matrix-td4548822.html. ''' a_null = numpy.zeros((len(self.__species), sizes)) for i in range(len(self.__species)): temp_vec = [] for j in range(sizes): temp_vec.append(the_null_space[j][i]) a_null[i, :] = temp_vec a_eq = numpy.array([numpy.sum(a_null, axis=0)]) # must multiply by negative one because in optimization we have the inequality <= 0.0 a_ub = -1.0*a_null b_ub = numpy.zeros(len(self.__species)) b_eq = numpy.array([1.0]) # defining the number of solutions to simulate num_sols = ((a_ub.shape[0]+1)*(a_ub.shape[1]+1))*10 # a matrix to hold the different solutions of the method sol = numpy.zeros((num_sols, len(self.__species))) # searching the nullspace for nonnegative vectors for i in range(num_sols): # generating a multivariate normal random variates crit = numpy.random.normal(0, 1.0, a_ub.shape[1]) # Solving the linear programming problem out = scipy.optimize.linprog(crit, A_eq=a_eq, b_eq=b_eq, A_ub=a_ub, b_ub=b_ub, bounds=(None, None), method='simplex') # multiplying our nullspace vectors with the minimized value # to create a nonnegative vector sol[i, :] = numpy.dot(a_null, out.x) sol[i, [numpy.abs(sol[i, :]) < numpy.finfo(float).eps][0]] = numpy.float64(0.0) # getting the unique vectors unique_vecs = numpy.unique(sol, axis=0) num_rows = numpy.size(unique_vecs, 0) # taking the smallest nonzero entry of the unique vectors # and dividing by it this will hopefully create nice # looking vectors that are whole numbers for i in range(num_rows): minval = numpy.min(unique_vecs[i, numpy.nonzero(unique_vecs[i, :])]) unique_vecs[i, :] = unique_vecs[i, :]/numpy.float64(minval) unique_vecs = numpy.unique(unique_vecs.round(decimals=10), axis=0) return unique_vecs def __create_mass_action_monomials(self): self.__psi = sympy.ones(len(self.__complexes), 1) for j in range(len(self.__complexes)): for i in range(len(self.__species)): self.__psi[j] *= sympy.Symbol(self.__species[i], positive=True)**self.__Y[i, j] # TODO: possibly make concetrations positive # public getters def get_dict_id_name(self): return self.__dic_species_id_name
[docs] def get_ode_system(self): """ Returns SymPy matrix representing the ODE system. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_ode_system()) ⎡ -re₁⋅s₁⋅s₂ + re1r⋅s₃ + re₄⋅s₁₆ - re₅⋅s₁⋅s₆ + re5r⋅s₁₅ ⎤ ⎢ ⎥ ⎢ -re₁⋅s₁⋅s₂ + s₃⋅(re1r + re₂) ⎥ ⎢ ⎥ ⎢ re₁⋅s₁⋅s₂ + s₃⋅(-re1r - re₂) ⎥ ⎢ ⎥ ⎢re₂⋅s₃ - re₃⋅s₆⋅s₇ + re3r⋅s₁₆ - re₅⋅s₁⋅s₆ + s₁₅⋅(re5r + 2⋅re₆)⎥ ⎢ ⎥ ⎢ -re₃⋅s₆⋅s₇ + s₁₆⋅(re3r + re₄) ⎥ ⎢ ⎥ ⎢ re₃⋅s₆⋅s₇ + s₁₆⋅(-re3r - re₄) ⎥ ⎢ ⎥ ⎣ re₅⋅s₁⋅s₆ + s₁₅⋅(-re5r - re₆) ⎦ """ return self.__Y*self.__A*self.__psi
[docs] def get_graph(self): """ Returns the NetworkX DiGraph representation of the network. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_graph() """ return self.__g
[docs] def get_species(self): """ Returns Python list of strings representing the species of the network. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> print(network.get_c_graph().get_species()) ['s1', 's2', 's3', 's6', 's7', 's16', 's15'] """ return self.__species
[docs] def get_complexes(self): """ Returns Python list of strings representing the complexes of the network. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> print(network.get_c_graph().get_complexes()) ['s1+s2', 's3', 's6+s2', 's6+s7', 's16', 's7+s1', 's1+s6', 's15', '2*s6'] """ return self.__complexes
[docs] def get_reactions(self): """ Returns Python list of strings representing the reactions of the network. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> print(network.get_c_graph().get_reactions()) ['re1', 're1r', 're2', 're3', 're3r', 're4', 're5', 're5r', 're6'] """ return self.__reactions
[docs] def get_a(self): """ Returns SymPy matrix representing the kinetic constant matrix, :math:`A`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_a()) ⎡-re₁ re1r 0 0 0 0 0 0 0⎤ ⎢ ⎥ ⎢re₁ -re1r - re₂ 0 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 re₂ 0 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 0 0 -re₃ re3r 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 0 0 re₃ -re3r - re₄ 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 0 0 0 re₄ 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 0 0 0 0 0 -re₅ re5r 0⎥ ⎢ ⎥ ⎢ 0 0 0 0 0 0 re₅ -re5r - re₆ 0⎥ ⎢ ⎥ ⎣ 0 0 0 0 0 0 0 re₆ 0⎦ """ return self.__A
[docs] def get_y(self): """ Returns SymPy matrix representing the molecularity matrix, :math:`Y`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_y()) ⎡1 0 0 0 0 1 1 0 0⎤ ⎢ ⎥ ⎢1 0 1 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢0 1 0 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢0 0 1 1 0 0 1 0 2⎥ ⎢ ⎥ ⎢0 0 0 1 0 1 0 0 0⎥ ⎢ ⎥ ⎢0 0 0 0 1 0 0 0 0⎥ ⎢ ⎥ ⎣0 0 0 0 0 0 0 1 0⎦ """ return self.__Y
[docs] def get_s(self): """ Returns SymPy matrix representing the stoichiometric matrix, :math:`S`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_s()) ⎡-1 1 0 0 0 1 -1 1 0 ⎤ ⎢ ⎥ ⎢-1 1 1 0 0 0 0 0 0 ⎥ ⎢ ⎥ ⎢1 -1 -1 0 0 0 0 0 0 ⎥ ⎢ ⎥ ⎢0 0 1 -1 1 0 -1 1 2 ⎥ ⎢ ⎥ ⎢0 0 0 -1 1 1 0 0 0 ⎥ ⎢ ⎥ ⎢0 0 0 1 -1 -1 0 0 0 ⎥ ⎢ ⎥ ⎣0 0 0 0 0 0 1 -1 -1⎦ """ return self.__S
[docs] def get_b(self): """ Returns SymPy matrix representing the mass conservation matrix, :math:`B`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_b()) ⎡ 0 0 0 0 1.0 1.0 0 ⎤ ⎢ ⎥ ⎢ 0 1.0 1.0 0 0 0 0 ⎥ ⎢ ⎥ ⎣1.0 0 1.0 1.0 0 1.0 2.0⎦ """ return self.__B
[docs] def get_lambda(self): """ Returns SymPy matrix representing the linkage class matrix, :math:`\Lambda`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_lambda()) ⎡1 0 0⎤ ⎢ ⎥ ⎢1 0 0⎥ ⎢ ⎥ ⎢1 0 0⎥ ⎢ ⎥ ⎢0 1 0⎥ ⎢ ⎥ ⎢0 1 0⎥ ⎢ ⎥ ⎢0 1 0⎥ ⎢ ⎥ ⎢0 0 1⎥ ⎢ ⎥ ⎢0 0 1⎥ ⎢ ⎥ ⎣0 0 1⎦ """ return self.__Lambda
[docs] def get_psi(self): """ Returns SymPy matrix representing the mass action monomials, :math:`\psi`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> import sympy >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> sympy.pprint(network.get_c_graph().get_psi()) ⎡s₁⋅s₂⎤ ⎢ ⎥ ⎢ s₃ ⎥ ⎢ ⎥ ⎢s₂⋅s₆⎥ ⎢ ⎥ ⎢s₆⋅s₇⎥ ⎢ ⎥ ⎢ s₁₆ ⎥ ⎢ ⎥ ⎢s₁⋅s₇⎥ ⎢ ⎥ ⎢s₁⋅s₆⎥ ⎢ ⎥ ⎢ s₁₅ ⎥ ⎢ ⎥ ⎢ 2 ⎥ ⎣ s₆ ⎦ """ return self.__psi
[docs] def get_g_nodes(self): """ Returns a list of strings that represent the order of the nodes of the NetworkX DiGraph. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_g_nodes() """ return self.__g_nodes
[docs] def get_g_edges(self): """ Returns a list of tuples of strings that represent the order of the edges of the NetworkX DiGraph. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_g_edges() """ return self.__g_edges
[docs] def get_weak_reversibility_of_linkage_classes(self): """ Returns list of Python boolean types for the weak reversibility of each linkage class. If the linkage class is weakly reversible then the entry in the list is True, False otherwise with order as defined by :func:`crnt4sbml.Cgraph.get_linkage_classes`. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_weak_reversibility_of_linkage_classes() """ return [networkx.is_strongly_connected(i) for i in self.__linkage_classes]
[docs] def get_if_cgraph_weakly_reversible(self): """ Returns weak reversibility of the network. If the network is weakly reversible True is returned, False otherwise. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_if_cgraph_weakly_reversible() """ return all(self.get_weak_reversibility_of_linkage_classes())
[docs] def get_dim_equilibrium_manifold(self): """ Returns integer value representing the dimension of the equilibrium manifold, :math:`\lambda`. This value is the number of mass conservation relationships. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> print(network.get_c_graph().get_dim_equilibrium_manifold()) 3 """ return len(self.__species) - self.__S.rank()
[docs] def get_deficiency(self): """ Returns integer value representing the deficiency of the network, :math:`\delta`. :download:`Fig1Ci.xml <../../sbml_files/Fig1Ci.xml>` for the provided example. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/Fig1Ci.xml") >>> print(network.get_c_graph().get_deficiency()) 2 """ return len(self.__complexes) - len(self.__linkage_classes) - self.__S.rank()
[docs] def get_linkage_classes(self): """ Returns list of NetworkX subgraphs representing the linkage classes. Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_linkage_classes() """ return self.__linkage_classes
[docs] def get_network_dimensionality_classification(self): """ Returns a two element list specifying the dimensionality of the network. Possible output: ["over-dimensioned",0] or ["proper",1] or ["under-dimensioned",2] or ["NOT DEFINED!",3] Example -------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_network_dimensionality_classification() """ temp = self.get_dim_equilibrium_manifold() - self.get_deficiency() if temp < 0: network_class = "over-dimensioned" classification = 0 elif temp == 0: network_class = "proper" classification = 1 elif temp > 0: network_class = "under-dimensioned" classification = 2 else: network_class = "NOT DEFINED!" classification = 3 return [network_class, classification]
[docs] def get_linkage_classes_deficiencies(self): """ Returns an interger list of each linkage class deficiency. Here, the first element corresponds to the first linkage class with order as defined by :func:`crnt4sbml.Cgraph.get_linkage_classes`. Example --------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_linkage_classes_deficiencies() """ deficiencies = [] for l in self.__linkage_classes: temp = [self.__g.edges[e]['label'] for e in l.edges] reactions_idx = [self.__reactions.index(i) for i in temp] deficiencies.append(networkx.number_of_nodes(l) - 1 - self.__S[:, reactions_idx].rank()) return deficiencies
[docs] def get_number_of_terminal_strong_lc_per_lc(self): """ Returns an integer list stating the number of terminally strong linkage classes per linkage class. Here, the first element corresponds to the first linkage class with order as defined by :func:`crnt4sbml.Cgraph.get_linkage_classes`. Example --------- >>> import crnt4sbml >>> network = crnt4sbml.CRNT("path/to/sbml_file.xml") >>> network.get_c_graph().get_number_of_terminal_strong_lc_per_lc() """ number_of_terminal_strong_lc_per_lc = [] for lc_i in self.__linkage_classes: num_term_strong = 0 strong_lcs = [lc_i.subgraph(c).copy() for c in networkx.strongly_connected_components(lc_i)] for slc_j in strong_lcs: slc_j_lbls = [i for i in slc_j.nodes] x = [lc_i.successors(v_i) for v_i in slc_j_lbls] connected_out_vertices = [item for sublist in x for item in sublist] d = set(connected_out_vertices) - set(slc_j_lbls) if d == set(): num_term_strong += 1 number_of_terminal_strong_lc_per_lc.append(num_term_strong) return number_of_terminal_strong_lc_per_lc
# public methods
[docs] def print(self): """ Prints edges and nodes of NetworkX DiGraph. See also --------- crnt4sbml.CRNT.print_c_graph """ print("") print("Reaction graph of the form") print("reaction -- reaction label:") for e in self.__g_edges: print(e[0] + ' -> ' + e[1] + ' -- ' + self.__g.edges[e]['label']) print("")
[docs] def plot(self): """ Plots NetworkX DiGraph. See also --------- crnt4sbml.CRNT.plot_c_graph """ pos = networkx.circular_layout(self.__g, scale=1.5) pos = networkx.kamada_kawai_layout(self.__g, pos=pos, scale=1.5) node_colors = [self.__g.nodes[n]['color'] for n in self.__g_nodes] networkx.draw(self.__g, pos, node_color=node_colors) node_labels = networkx.get_node_attributes(self.__g, 'label') networkx.draw_networkx_labels(self.__g, pos, labels=node_labels) edge_labels = networkx.get_edge_attributes(self.__g, 'label') networkx.draw_networkx_edge_labels(self.__g, pos, edge_labels, label_pos=0.3) plt.show()
[docs] def plot_save(self): """ Saves the plot of the NetworkX DiGraph. See also --------- crnt4sbml.CRNT.plot_save_c_graph """ pos = networkx.circular_layout(self.__g, scale=1.5) pos = networkx.kamada_kawai_layout(self.__g, pos=pos, scale=1.5) node_colors = [self.__g.nodes[n]['color'] for n in self.__g_nodes] networkx.draw(self.__g, pos, node_color=node_colors) node_labels = networkx.get_node_attributes(self.__g, 'label') networkx.draw_networkx_labels(self.__g, pos, labels=node_labels) edge_labels = networkx.get_edge_attributes(self.__g, 'label') networkx.draw_networkx_edge_labels(self.__g, pos, edge_labels, label_pos=0.3) plt.savefig('network_cgraph.png')