partial fix kinetic models

Author: PHvanLent

scripts/build_pathway_modelB.py

@@ -76,7 +76,7 @@
     name="v_biomass",
     species=["c_A","c_atp","c_biomass","co2"],
     stoichiometry = [-0.175, -1, 1,  0.05],
-    compartments = ['c','c','c', 'c','c'],
+    compartments = ['c','c','c', 'c'],
     mechanism=mech_biomass, )
 
 mech_maintenance = mechanisms.Jax_MM_Irrev_Uni(substrate='c_atp',

scripts/build_pathway_modelC.py

@@ -90,7 +90,7 @@
 v_maintenance = build.Reaction(
     name="v_maintenance",
     species=["c_A", "c_atp", "co2",'c_adp', 'c_biomass'],
-    stoichiometry=[-1/6, -1, 1,
+    stoichiometry=[0, -1, 1,
                    1, 0],
     compartments=['c', 'c', 'c', 'c', 'c'],
     mechanism=mech_maintenance, )
@@ -105,9 +105,9 @@
 
 v_respiration = build.Reaction(
     name="v_respiration",
-    species=["c_A", "c_adp","co2", "c_atp", 'c_biomass'],
-    stoichiometry=[-1,-28, 6, 28, 0],
-    compartments=['c', 'c', 'c', 'c', 'c'],
+    species=["c_A","co2","c_atp", 'c_biomass'],
+    stoichiometry=[-1, 6, 1, 0],
+    compartments=['c', 'c', 'c', 'c'],
     mechanism=mech_respiration,
 )
 
@@ -435,9 +435,9 @@
 
 v_nadh_regeneration = build.Reaction(
     name="v_atpsynthase",
-    species= ['c_nadh','c_nad','c_atp','c_biomass'],
-    compartments = ['c', 'c', 'c', 'c'],
-    stoichiometry=[-1, 1, 3, 0],
+    species= ['c_nadh','c_nad','c_biomass'],
+    compartments = ['c', 'c', 'c'],
+    stoichiometry=[-1, 1, 0],
     mechanism=mech_nadh_regeneration,
 )
 
@@ -570,13 +570,13 @@
 
 protein_load = build.AssignmentRule(
     rule=("0.3*"
-          "(((kcat_1 *enzyme_1) + (kcat_2 * enzyme_2) + (kcat_3 * enzyme_3) + (kcat_4*enzyme_4) "
-          "+ (kcat_5* enzyme_5) +"
-          "(kcat_6 * enzyme_6) + (kcat_7 * enzyme_7) + (kcat_8 * enzyme_8) + "
-          "(kcat_9 * enzyme_9) + (kcat_10 * enzyme_10) + (kcat_11 * enzyme_11) +"
-          "(kcat_12 * enzyme_12) + (kcat_13 * enzyme_13)+ (kcat_14 * enzyme_14)+"
-          "(kcat_15 * enzyme_15) + (kcat_16 * enzyme_16) + (kcat_17 * enzyme_17)+"
-          "+90)/100)"))
+          "((enzyme_1 + enzyme_2 + enzyme_3 + enzyme_4"
+          "+  enzyme_5 +"
+          "enzyme_6+  enzyme_7 +  enzyme_8 + "
+          "enzyme_9 + enzyme_10 + enzyme_11 +"
+          "enzyme_12 +  enzyme_13+  enzyme_14+"
+          "enzyme_15 +  enzyme_16 +  enzyme_17+"
+          "+1)/17)"))
 
 kmodel.add_parameter_assignment_rule(variable="vmax_maintenance",
                                      rule=protein_load)
@@ -751,7 +751,7 @@
 
 for p in increasing_protein_load:
     newparameters = model.parameters.copy()
-    newparameters['enzyme_12'] = p
+    newparameters['enzyme_8'] = p
 
     ys = JaxKmodel(ts=ts, y0=model.y0, params=newparameters)
     ys = pd.DataFrame(ys, columns=S.index)

scripts/build_pathway_modelD.py

@@ -138,10 +138,9 @@
 
 v_maintenance = build.Reaction(
     name="v_maintenance",
-    species=["c_A", "c_atp", "co2",'c_adp', 'c_biomass'],
-    stoichiometry=[-1/6, -1, 1,
-                   1, 0],
-    compartments=['c', 'c', 'c', 'c', 'c'],
+    species=["c_A", "c_atp", "co2",'c_biomass'],
+    stoichiometry=[0, -1, 1, 0],
+    compartments=['c', 'c', 'c', 'c'],
     mechanism=mech_maintenance, )
 
 mech_respiration = mechanisms.Jax_MM_Irrev_Bi(substrate1='c_A',
@@ -154,12 +153,11 @@
 
 v_respiration = build.Reaction(
     name="v_respiration",
-    species=["c_A", "c_adp",
-             "co2", "c_atp",
+    species=["c_A",
+             "co2",
+             "c_atp",
              'c_biomass'],
-    stoichiometry=[-1, -28,
-                   6, 28,
-                   0],
+    stoichiometry=[-1, 6, 1, 0],
     compartments=['c', 'c', 'c', 'c', 'c'],
     mechanism=mech_respiration,
 )
@@ -555,7 +553,7 @@
 v_nadh_regeneration = build.Reaction(
     name="v_atpsynthase",
     species= ['c_nad','c_nadh','c_biomass'],
-    compartments = ['c', 'c', 'c', 'c'],
+    compartments = ['c', 'c', 'c'],
     stoichiometry=[-1, 1,  0],
     mechanism=mech_nadh_regeneration,
 )
@@ -704,17 +702,15 @@
 #we sum all vmaxes as a protein load, divide this by the sum of vmaxes (which
 #for the baseline returns one, and then we multiply by 0.2 to get the vmax maintenance.
 
-
 protein_load = build.AssignmentRule(
     rule=("0.3*"
-          "(((kcat_1 *enzyme_1) + (kcat_2 * enzyme_2) + (kcat_3 * enzyme_3) + (kcat_4*enzyme_4) "
-          "+ (kcat_5* enzyme_5) +"
-          "(kcat_6 * enzyme_6) + (kcat_7 * enzyme_7) + (kcat_8 * enzyme_8) + "
-          "(kcat_9 * enzyme_9) + (kcat_10 * enzyme_10) + (kcat_11 * enzyme_11) +"
-          "(kcat_12 * enzyme_12) + (kcat_13 * enzyme_13)+ (kcat_14 * enzyme_14)+"
-          "(kcat_15 * enzyme_15) + (kcat_16 * enzyme_16) + (kcat_17 * enzyme_17)+(kcat_18 *enzyme_18)+"
-          "(kcat_19 * enzyme_19) + (kcat_20 * enzyme_20)+ (kcat_21 * enzyme_21)+"
-          "+1)/100)"))
+          "((enzyme_1 + enzyme_2 + enzyme_3 + enzyme_4 "
+          "+ enzyme_5 + enzyme_6 + enzyme_7 + enzyme_8 + "
+          "enzyme_9 + enzyme_10 + enzyme_11 + enzyme_12 + "
+          "enzyme_13 + enzyme_14 + enzyme_15 + enzyme_16 + "
+          "enzyme_17 + enzyme_18 + enzyme_19 + enzyme_20 + enzyme_21 "
+          "+1)/100)")
+)
 
 kmodel.add_parameter_assignment_rule(variable="vmax_maintenance",
                                      rule=protein_load)

scripts/build_pathway_modelE.py

@@ -64,20 +64,23 @@
 
 
 # now biomass, from the internal carbon source species
-mech_biomass = Jax_MM_Irrev_Uni_e(substrate="c_A",
+mech_biomass = Jax_MM_Irrev_Bi_e(substrate1="c_A",
+                                 substrate2="c_atp",
                                   kcat="kcat_1",
                                   enzyme="enzyme_1",
-                                  km_substrate="km_v1_gluc")
+                                  km_substrate1="km_v1_gluc",
+                                 km_substrate2="km_v1_atp")
 mech_biomass.add_modifier(modifier.BasicDivision(symbol="a", scaling="scale"))
 mech_biomass.add_modifier(modifier.BasicMultiplication(symbol="c_biomass", scaling="scale"))
-mech_biomass.add_modifier(modifier.SimpleInhibitor(inhibitor="toxic_byproduct",k_I="ki_toxicbyproduct"))
 v_biomass = build.Reaction(
     name="v_biomass",
-    species=["c_A","c_biomass","co2"],
-    stoichiometry = [-0.175, 1,  0.05],
+    species=["c_A","c_atp","c_biomass","co2"],
+    stoichiometry = [-0.175, -1, 1,  0.05],
     compartments = ['c','c','c', 'c','c'],
     mechanism=mech_biomass, )
 
+
+
 mech_maintenance = mechanisms.Jax_MM_Irrev_Bi(substrate1='c_A',
                                               substrate2= 'c_atp',
                                                vmax="vmax_maintenance",
@@ -90,28 +93,22 @@
 v_maintenance = build.Reaction(
     name="v_maintenance",
     species=["c_A", "c_atp", "co2",'c_adp', 'c_biomass'],
-    stoichiometry=[-1/6, -28, 1,
-                   28, 0],
-    compartments=['c', 'c', 'c', 'c', 'c'],
+    stoichiometry=[0, -1, 1, 0],
+    compartments=['c', 'c', 'c', 'c'],
     mechanism=mech_maintenance, )
 
-mech_respiration = mechanisms.Jax_MM_Irrev_Bi(substrate1='c_A',
-                                              substrate2= "c_adp",
+mech_respiration = mechanisms.Jax_MM_Irrev_Uni(substrate='c_A',
                                                vmax="vmax_respiration",
-                                               km_substrate1="km_respiration_A",
-                                              km_substrate2="km_respiration_adp",)
+                                               km_substrate="km_respiration_A")
 mech_respiration.add_modifier(modifier.BasicDivision(symbol="a", scaling="scale"))
 mech_respiration.add_modifier(modifier.BasicMultiplication(symbol="c_biomass", scaling="scale"))
 
 v_respiration = build.Reaction(
     name="v_respiration",
-    species=["c_A", "c_adp",
-             "co2", "c_atp",
+    species=["c_A","co2", "c_atp",
              'c_biomass'],
-    stoichiometry=[-1, -28,
-                   6, 28,
-                   0],
-    compartments=['c', 'c', 'c', 'c', 'c'],
+    stoichiometry=[-1, 6,1,0],
+    compartments=['c', 'c', 'c', 'c'],
     mechanism=mech_respiration,
 )
 
@@ -555,9 +552,9 @@
 
 v_nadh_regeneration = build.Reaction(
     name="v_atpsynthase",
-    species= ['c_nadh','c_adp','c_nad', 'c_atp', 'c_biomass'],
-    compartments = ['c', 'c', 'c', 'c'],
-    stoichiometry=[-1, -1, 1, 1, 0],
+    species= ['c_nadh','c_nad', 'c_biomass'],
+    compartments = ['c', 'c', 'c'],
+    stoichiometry=[-1, 1,  0],
     mechanism=mech_nadh_regeneration,
 )
 
@@ -569,17 +566,17 @@
     "kcat_1": 0.3,
     "enzyme_1": 1,
     "km_v1_gluc": 0.01,
-
+    "km_v1_atp": 0.005,
     "a": 1.6,
     "scale": 1,  # is the broth volume
     # "respiration_flux": 0.02,
     "vmax_respiration": 0.35,
     "km_respiration_A": 0.6,
-    "km_respiration_adp":0.1,
     "vmax_maintenance": 0, #has assignment rule
     "km_vmaintenance_atp":0.1,
     "km_vmaintenance_gluc": 0.1,
     "km_v2_A": 0.01,
+
     "km_v2_adp":0.01,
     "kcat_2": 0.1,
     "enzyme_2": 1.0,
@@ -706,35 +703,23 @@
 #for the baseline returns one, and then we multiply by 0.2 to get the vmax maintenance.
 
 
-# protein_load = build.AssignmentRule(
-#     rule=("0.3*"
-#           "(((kcat_1 *enzyme_1) + (kcat_2 * enzyme_2) + (kcat_3 * enzyme_3) + (kcat_4*enzyme_4)"
-#           "+(kcat_5* enzyme_5) +"
-#           "(kcat_6 * enzyme_6) + (kcat_7 * enzyme_7) + (kcat_8 * enzyme_8) +"
-#           "(kcat_9 * enzyme_9) + (kcat_10 * enzyme_10) + (kcat_11 * enzyme_11) +"
-#           "(kcat_12 * enzyme_12) + (kcat_13 * enzyme_13)+ (kcat_14 * enzyme_14)+"
-#           "(kcat_15 * enzyme_15) + (kcat_16 * enzyme_16) + (kcat_17 * enzyme_17) + (kcat_18 * enzyme_18) +"
-#           "(kcat_19 * enzyme_19) + (kcat_20 * enzyme_20) + (kcat_21 * enzyme_21) +"
-#           "(kcat_22 * enzyme_22)+ (kcat_23 enzyme_23)+ (kcat_24 * enzyme_24)+ (kcat_25 * enzyme_25)"
-#           "+75)/100)"))
 
 protein_load = build.AssignmentRule(
-    rule=("0.2*"
-          "(((kcat_1 *enzyme_1) + (kcat_2 * enzyme_2) + (kcat_3 * enzyme_3) + (kcat_4*enzyme_4) "
-          "+ (kcat_5* enzyme_5) +"
-          "(kcat_6 * enzyme_6) + (kcat_7 * enzyme_7) + (kcat_8 * enzyme_8) + "
-          "(kcat_9 * enzyme_9) + (kcat_10 * enzyme_10) + (kcat_11 * enzyme_11) +"
-          "(kcat_12 * enzyme_12) + (kcat_13 * enzyme_13)+ (kcat_14 * enzyme_14)+"
-          "(kcat_15 * enzyme_15) + (kcat_16 * enzyme_16) + (kcat_17 * enzyme_17)+(kcat_18 *enzyme_18)+"
-          "(kcat_22 * enzyme_22)+ (kcat_23* enzyme_23)+ (kcat_24 * enzyme_24)+ (kcat_25 * enzyme_25)+"
-          "+49.4991)/55)"))
+    rule=("0.3*"
+          "((enzyme_1 + enzyme_2 + enzyme_3 + enzyme_4 "
+          "+ enzyme_5 + enzyme_6 + enzyme_7 + enzyme_8 + "
+          "enzyme_9 + enzyme_10 + enzyme_11 + enzyme_12 + "
+          "enzyme_13 + enzyme_14 + enzyme_15 + enzyme_16 + "
+          "enzyme_17 + enzyme_18 + enzyme_22 + enzyme_23 + "
+          "enzyme_24 + enzyme_25 + 12)/25)")
+)
 
 kmodel.add_parameter_assignment_rule(variable="vmax_maintenance",
                                      rule=protein_load)
 
 kmodel_sim = build.NeuralODEBuild(kmodel)
 
-ts = jnp.linspace(0, 50, 1000)
+ts = jnp.linspace(0, 90, 1000)
 
 print(kmodel_sim.species_names)  #check unit consistency (cmol/L) later
 print('number of parameters', len(parameters))
@@ -910,26 +895,26 @@
 norm = colors.Normalize(vmin=increasing_protein_load.min(), vmax=increasing_protein_load.max())
 cmap = plt.cm.Reds  # or any other colormap you like
 #
-# for p in increasing_protein_load:
-#     print(p)
-#     newparameters = model.parameters.copy()
-#     newparameters['enzyme_5'] = p
-#
-#     ys = JaxKmodel(ts=ts, y0=model.y0, params=newparameters)
-#     ys = pd.DataFrame(ys, columns=S.index)
+for p in increasing_protein_load:
+    print(p)
+    newparameters = model.parameters.copy()
+    newparameters['enzyme_5'] = p
+
+    ys = JaxKmodel(ts=ts, y0=model.y0, params=newparameters)
+    ys = pd.DataFrame(ys, columns=S.index)
+
+    color = cmap(norm(p))  # map protein load to a color
+    ax.plot(ts, ys['product_E'], color=color, label=f'p={p:.2f}')
 #
-#     color = cmap(norm(p))  # map protein load to a color
-#     ax.plot(ts, ys['product_E'], color=color, label=f'p={p:.2f}')
-# #
-# # # Optional: add colorbar
-# sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
-# sm.set_array([])  # Needed for ScalarMappable in some versions
-# cbar = plt.colorbar(sm, ax=ax)
-# cbar.set_label('Protein Load (vmax_maintenance)')
+# # Optional: add colorbar
+sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
+sm.set_array([])  # Needed for ScalarMappable in some versions
+cbar = plt.colorbar(sm, ax=ax)
+cbar.set_label('Protein Load (vmax_maintenance)')
+
+ax.set_title("Protein load (0-100 vmax3)")
+ax.set_xlabel('Time')
+ax.set_ylabel('Product')
+plt.tight_layout()
 #
-# ax.set_title("Protein load (0-100 vmax3)")
-# ax.set_xlabel('Time')
-# ax.set_ylabel('Product')
-# plt.tight_layout()
-# #
-# fig.savefig(f"{output_dir}/protein_load_effect.png", bbox_inches='tight')
+fig.savefig(f"{output_dir}/protein_load_effect.png", bbox_inches='tight')