In [51]:
%reload_ext autoreload
%aimport mentos
%autoreload 1
from IPython.display import Image, display
import sys
import pandas as pd
Image('Metabolic-network.JPG')
Out[51]:
The chemical equations for the ABC model are:
$$ \begin{eqnarray} A_{ext} & \overset{\overset{+1}{\rightharpoonup}}{\underset{-1}{\leftharpoondown}} & A \\ A & \overset{\overset{+2}{\rightharpoonup}}{\underset{-2}{\leftharpoondown}} & B \\ A & \overset{\overset{+3}{\rightharpoonup}}{\underset{-3}{\leftharpoondown}} & C \\ B + E & \overset{\overset{+4}{\rightharpoonup}}{\underset{-4}{\leftharpoondown}}& 2D\\ E_{ext} & \overset{\overset{+5}{\rightharpoonup}}{\underset{-5}{\leftharpoondown}} & E \\ 2B & \overset{\overset{+6}{\rightharpoonup}}{\underset{-6}{\leftharpoondown}} & C + F \\ C &\overset{\overset{+7}{\rightharpoonup}}{\underset{-7}{\leftharpoondown}} & D \\ D & \overset{\overset{+8}{\rightharpoonup}}{\underset{-8}{\leftharpoondown}}& D_{ext} \\ F &\overset{\overset{+9}{\rightharpoonup}}{\underset{-9}{\leftharpoondown}} & F_{ext} \\ \end{eqnarray} $$which can be represented as a Stoichiometric matrix $S_{full}$:
$$ S_{full} = \left[ {\begin{array}{cccccccccc} & R_1 & R_2 & R_3 & R_4 & R_5 & R_6 & R_7 & R_8 & R_9 \\ A & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ B & 0 & 1 & 0 & -1 & 0 & -2 & 0 & 0 & 0 \\ C & 0 & 0 & 1 & 0 & 0 & 1 & -1 & 0 & 0 \\ D & 0 & 0 & 0 & 2 & 0 & 0 & 1 & -1 & 0 \\ E & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\ F & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -1 \\ A_{ext} & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ E_{ext} & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ D_{ext} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ F_{ext} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array}}\right] $$
In [52]:
import escher,os
import cobra
abc_model = cobra.Model('ABC_model')
M = {}
A_ext = cobra.Metabolite('A_ext', compartment='ext')
E_ext = cobra.Metabolite('E_ext', compartment= 'ext')
D_ext = cobra.Metabolite('D_ext', compartment='ext')
F_ext = cobra.Metabolite('F_ext', compartment='ext')
A = cobra.Metabolite('A',compartment='c')
B = cobra.Metabolite('B',compartment='c')
C = cobra.Metabolite('C',compartment='c')
D = cobra.Metabolite('D',compartment='c')
E = cobra.Metabolite('E',compartment='c')
F = cobra.Metabolite('F',compartment='c')
R_1 = cobra.Reaction('R_1')
R_2 = cobra.Reaction('R_2')
R_3 = cobra.Reaction('R_3')
R_4 = cobra.Reaction('R_4')
R_5 = cobra.Reaction('R_5')
R_6 = cobra.Reaction('R_6')
R_7 = cobra.Reaction('R_7')
R_8 = cobra.Reaction('R_8')
R_9 = cobra.Reaction('R_9')
R_1.add_metabolites({A_ext: -1, A: 1})
R_2.add_metabolites({A:-1, B:1})
R_3.add_metabolites({A:-1, C:1})
R_4.add_metabolites({B:-1, E:-1, D:2})
R_5.add_metabolites({E_ext:-1, E:1})
R_6.add_metabolites({B:-2, C:1,F:1})
R_7.add_metabolites({C:-1,D:1})
R_8.add_metabolites({D:-1, D_ext:1})
R_9.add_metabolites({F:-1, F_ext: 1})
abc_model.add_reactions([R_1, R_2, R_3, R_4, R_5, R_6, R_7, R_8, R_9])
cobra.io.save_json_model(abc_model, 'Mentos/ABC_model.json')
cobra.io.write_sbml_model(abc_model, 'Mentos/ABC_model.sbml')
abc_json = cobra.io.to_json(abc_model)
In [53]:
abc_array = abc_model.to_array_based_model()
abc_array.S
Out[53]:
In [54]:
sys.path.append('')
mets = ['A','B','C','D','E','F','A_ext', 'D_ext', 'E_ext','F_ext']
internal_mets = [m for m in mets if 'ext' not in m]
rxns = ['R_{}'.format(i) for i in range(1,10)]
data = {'R_1': pd.Series({'A':1, 'A_ext': -1}),
'R_2': pd.Series({'A':-1,'B':1}),
'R_3': pd.Series({'A':-1,'C':1}),
'R_4': pd.Series({'B':-1, 'D': 2, 'E': -1}),
'R_5': pd.Series({'E':1, 'E_ext':-1}),
'R_6': pd.Series({'B':-2, 'C':1, 'F': 1}),
'R_7': pd.Series({'C':-1, 'D':1}),
'R_8': pd.Series({'D': -1, 'D_ext': 1}),
'R_9': pd.Series({'F':-1, 'F_ext':1})}
fullS = pd.DataFrame(data, columns=rxns, index=mets,dtype='int64').fillna(0)
biomass = fullS.columns.get_loc('R_8')
A_uptake = 0
E_uptake = 4
R = 8.3144598/1000.0 # ideal gas constant
T = 298.15 # standard temperature
n_A = 6.022e23 # Avogadro's number
V = 1e-15 # volume of cell in Liters
c_L = 1e-8 # lower bound of metabolite concentrations
c_U = 1e-3 # upper bound of metabolite concentrations
v_L = -100
v_U = 100
A_ext = 6 # Index of A_ext
E_ext = 8 # Index of E_ext
D_ext = 7 # Index of D_ext
F_ext = 9 # Index of F_ext
A,B,C,D,E,F = range(6)
lambda_x = 0.5
S = fullS.loc[internal_mets].as_matrix()
m,n = fullS.shape
metab = {}
reactions = {}
true_metab, true_reactions = {},{}
display(fullS)
How do we predict metabolite concentrations? First we need to to obtain the chemical potential for each metabolite.
Chemical potential $\vec{\mu}^0$ is the Gibbs energy of formation $\Delta G^0_+$ for each metabolite: $$ \vec{\mu}^0 = \Delta \vec{G^0_+} = \left[ {\begin{array}{cc} A & 4 \\ B & 2 \\ C & 2 \\ D & 0 \\ E & 2 \\ F & 1 \\ \end{array}}\right] $$
Which implies a Change in Standard Gibbs free energy $S^T\vec{\mu}^0 = \Delta G^0$ for each reaction: $$ \Delta \vec{G}^0 = \left[ {\begin{array}{cc} R_1 & 0 \\ R_2 & -2 \\ R_3 & -2 \\ R_4 & -4 \\ R_5 & 0 \\ R_6 & -1 \\ R_7 & -2 \\ R_8 & 0 \\ R_9 & 0 \\ \end{array}}\right] $$
In [67]:
mu0 = pd.Series([4.0,2,2,0,2,1,4,0,2,1], index=mets,dtype='float64')
display(pd.DataFrame(mu0,columns=['$\mu^0$']))
deltaG0 = fullS.T.dot(mu0)
display(pd.DataFrame(deltaG0,columns=['$\Delta G^0$']))
In [57]:
import pyOpt
from scipy.stats import entropy
import numpy as np
import cvxpy as cvx
import pandas as pd
import numpy as np
import sys, cobra
from IPython.display import display, HTML
#import pandas.core.format as fmt
pd.options.display.float_format = '{:.3g}'.format
#pd.options.display.float_+ormat = '{:.0f}'.format
pd.options.display.float_format = '{:.3g}'.format
Steady-state fluxes can be predicted from flux balance analysis:
$$\begin{array}{ll} \underset{\vec{\bf v}}{\mbox{maximize}} & v_{biomass} \\ \mbox{subject to} & S\cdot v = 0 \\ & 0 \leq v \leq v_{upper} \\ \end{array}$$
In [58]:
v = cvx.Variable(n)
obj = cvx.Maximize(v[biomass])
constraints = [S*v == 0,
0 <= v,
v <= v_U]
prob = cvx.Problem(obj, constraints)
prob.solve(verbose=True)
v.value
Out[58]:
At steady state, the log likelihood is equal to the difference of the logs of the equilibrium constant $K$ and the mass action ratio $Q$, subject to the boundary conditions
$$\begin{array}{ll} \underset{\log\vec{\bf L}_+, \log\vec{\bf c}}{\mbox{minimize}} & \|\log L_+\|_{2} \\ \mbox{subject to} & \log L_+ = \log K_{eq} - \log Q_r \\ & \log [A_{ext}] = \log c_{upper} \\ & \log [E_{ext}] = \log c_{upper} \\ & \log [D_{ext}] = \log c_{lower} \\ & \log [F_{ext}] = \log c_{lower} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_A \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_B \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_C \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_D \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_E \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_F \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ \end{array}$$Where,
$$\begin{eqnarray} L_+ & = & \frac{r_+}{r_-} \\ \log K_{eq} & = & -\frac{1}{RT}S^T\cdot\mu^0 \\ \log Q_r & = & S^T\cdot\log c \\ \mu & = & RT\log c + \mu^0 \\ \end{eqnarray}$$
In [59]:
log_c = cvx.Variable(m) # log of the concentrations
log_likelihood = cvx.Variable(n)
log_Q = fullS.as_matrix().T*log_c # log of the Reaction quotient
log_K = cvx.Constant(-1.0/(R*T)*deltaG0.as_matrix()) # log of the Equilibrium constnat
mu = R*T*log_c + mu0.values
#obj = cvx.Maximize(cvx.sum_entries(cvx.entr(log_likelihood))) # entropy of the log likelihoods
obj = cvx.Minimize(cvx.norm2(log_likelihood))
constraints = [
log_c[A_ext] == cvx.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] == cvx.log(c_U),
log_c[D_ext] == cvx.log(c_L),
log_c[F_ext] == cvx.log(c_L),
log_likelihood == log_K - log_Q, # Thermodynamic constraint
mu[A] >= mu[D_ext] + mu[F_ext], # Energy sink constraint
mu[B] >= mu[D_ext] + mu[F_ext],
mu[C] >= mu[D_ext] + mu[F_ext],
mu[D] >= mu[D_ext] + mu[F_ext],
mu[E] >= mu[D_ext] + mu[F_ext],
mu[F] >= mu[D_ext] + mu[F_ext],
mu[A] <= mu[A_ext] + mu[E_ext], # Energy barrier constraint
mu[B] <= mu[A_ext] + mu[E_ext],
mu[C] <= mu[A_ext] + mu[E_ext],
mu[D] <= mu[A_ext] + mu[E_ext],
mu[E] <= mu[A_ext] + mu[E_ext],
mu[F] <= mu[A_ext] + mu[E_ext],
]
prob = cvx.Problem(obj, constraints)
prob.solve( verbose=True)
display(pd.DataFrame({'$\log c$':np.squeeze(np.asarray(log_c.value))},index=mets))
display(pd.DataFrame({'$\log L_+$':np.squeeze(np.asarray(log_likelihood.value))},index=rxns))
In [60]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = dict([(rxns[i], R*T*log_likelihood[i].value) for i in range(len(rxns))])
metabolite_data = (mu0 + R*T*np.squeeze(np.asarray(log_c.value))).to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data, metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[60]:
At steady state, the log likelihood is equal to the difference of the logs of the equilibrium constant $K$ and the mass action ratio $Q$, subject to the boundary conditions
$$\begin{array}{ll} \underset{\log\vec{\bf L}_+, \log\vec{\bf c}}{\mbox{minimize}} & \sum_j \text{entropy}(\log L_+) \\ \mbox{subject to} & \log L_+ = \log K_{eq} - \log Q_r \\ & \log [A_{ext}] = \log c_{upper} \\ & \log [E_{ext}] = \log c_{upper} \\ & \log [D_{ext}] = \log c_{lower} \\ & \log [F_{ext}] = \log c_{lower} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_A \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_B \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_C \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_D \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_E \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ & \mu_{D_{ext}} + \mu_{F_{ext}} \leq \mu_F \leq \mu_{A_{ext}} + \mu_{E_{ext}} \\ \end{array}$$Where,
$$\begin{eqnarray} L_+ & = & \frac{r_+}{r_-} \\ \log K_{eq} & = & -\frac{1}{RT}S^T\cdot\mu^0 \\ \log Q_r & = & S^T\cdot\log c \\ \mu & = & RT\log c + \mu^0 \\ \end{eqnarray}$$
In [61]:
log_c = cvx.Variable(m) # log of the concentrations
log_likelihood = cvx.Variable(n)
log_Q = fullS.as_matrix().T*log_c # log of the Reaction quotient
log_K = cvx.Constant(-1.0/(R*T)*deltaG0.as_matrix()) # log of the Equilibrium constnat
mu = R*T*log_c + mu0.values
obj = cvx.Maximize(cvx.sum_entries(cvx.entr(log_likelihood))) # entropy of the log likelihoods
#obj = cvx.Minimize(cvx.norm1(log_likelihood))
constraints = [
log_c[A_ext] == cvx.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] == cvx.log(c_U),
log_c[D_ext] == cvx.log(c_L),
log_c[F_ext] == cvx.log(c_L),
log_likelihood == log_K - log_Q, # Thermodynamic constraint
mu[A] >= mu[D_ext] + mu[F_ext], # Energy sink constraint
mu[B] >= mu[D_ext] + mu[F_ext],
mu[C] >= mu[D_ext] + mu[F_ext],
mu[D] >= mu[D_ext] + mu[F_ext],
mu[E] >= mu[D_ext] + mu[F_ext],
mu[F] >= mu[D_ext] + mu[F_ext],
mu[A] <= mu[A_ext] + mu[E_ext], # Energy barrier constraint
mu[B] <= mu[A_ext] + mu[E_ext],
mu[C] <= mu[A_ext] + mu[E_ext],
mu[D] <= mu[A_ext] + mu[E_ext],
mu[E] <= mu[A_ext] + mu[E_ext],
mu[F] <= mu[A_ext] + mu[E_ext],
]
prob = cvx.Problem(obj, constraints)
prob.solve( verbose=True)
display(pd.DataFrame({'$\log c$':np.squeeze(np.asarray(log_c.value))},index=mets))
display(pd.DataFrame({'$\log L_+$':np.squeeze(np.asarray(log_likelihood.value))},index=rxns))
In [62]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = dict([(rxns[i], R*T*log_likelihood[i].value) for i in range(len(rxns))])
metabolite_data = (mu0 + R*T*np.squeeze(np.asarray(log_c.value))).to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data, metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[62]:
In [63]:
v_L, v_U = mentos.get_rxn_bounds_from_log_likelihood( np.squeeze(np.asarray(log_likelihood.value)))
v_L, v_U
Out[63]:
In [64]:
v = cvx.Variable(n)
obj = cvx.Maximize(v[biomass])
constraints = [S*v == 0,
v_L <= v,
v <= v_U]
prob = cvx.Problem(obj, constraints)
prob.solve(verbose=True)
v.value
Out[64]:
In [ ]:
From the definition of forward likelihood and net flux:
$$\begin{eqnarray} L_+ &=& \frac{r_{+}}{r_{-}} \\ v &=& r_{+} - r_{-} \\ \end{eqnarray}$$We can solve for the forward reaction rates $r_{+}$ and $r_{-}$: $$\begin{eqnarray} r_{-} & = & \frac{r_{+}}{L_+} \\ r_{+} - r_{-} & = & v\\ r_{+} - \frac{r_{+}}{L_+} & = & v \\ r_{+}(1 -\frac{1}{L_+}) & = & v \\ r_{+} & = & \frac{v}{1-\frac{1}{L_+}} \\ & = & \frac{L_+v}{L_+ - 1} \\ r_{-} & = & \frac{v}{L_+ - 1} \\ \end{eqnarray}$$
In [65]:
net_flux = np.squeeze(np.asarray(v.value))
forward_likelihood = np.exp(np.squeeze(np.asarray(log_likelihood.value)))
net_flux, forward_likelihood
forward_rate = net_flux*forward_likelihood/(forward_likelihood - 1)
backward_rate = net_flux/(forward_likelihood - 1)
forward_rate - backward_rate
Out[65]:
Let $r_+ = v(N+1)$ and $r_- = vN$. By definition, the thermodynamic likelihood is equal to ratio of the forward and backward rates, so
$$\begin{eqnarray} \frac{r_+}{r_-} & = & \frac{v(N+1)}{vN} \\ & = & \frac{N+1}{N} \\ & = & 1 + \frac{1}{N} \end{eqnarray}$$The thermodynamic constraint at steady state is
$$\begin{eqnarray} -\frac{1}{RT}S^T\mu^0 -S^T\log a & = & \log \frac{r_+}{r_-} \\ & = & \log \frac{N+1}{N} \\ \end{eqnarray}$$
In [86]:
reload( mentos)
metab['convex_mentos'] = mentos.generate_metabolite_report( np.squeeze(np.asarray(log_c.value))
, forward_rate, backward_rate, S, mets, internal_mets,
mu0).astype(np.float64)
metab['convex_mentos']['Chemical potential'].to_csv('Mentos/ABC_chemical_potentials.csv',header=True)
true_metab['convex_mentos']= pd.read_csv('ABC_metab_convex_mentos.csv',index_col=0).astype(np.float64)
pd.testing.assert_frame_equal(true_metab['convex_mentos'], metab['convex_mentos'])
mentos.compare_frames(true_metab['convex_mentos'].astype(np.float64),metab['convex_mentos'])
Out[86]:
In [76]:
reactions['convex_mentos'] = mentos.generate_rxn_report(mets, log_c, log_Q.value, log_K.value,
forward_rate, backward_rate, rxns, deltaG0,'R_8')
reactions['convex_mentos']['Delta G'].to_csv(
'Mentos/ABC_DeltaG.csv',header=True)
reactions['convex_mentos']['Net flux'].to_csv(
'Mentos/ABC_net_flux.csv',header=True)
reactions['convex_mentos']['Macroscopic Reaction Entropy Production Net Flux'].to_csv(
'Mentos/ABC_maxent_prod_net_flux.csv',header=True)
true_reactions['convex_mentos']= pd.read_csv('ABC_reactions_convex_mentos.csv',index_col=0)
pd.testing.assert_frame_equal(true_reactions['convex_mentos'], reactions['convex_mentos'])
mentos.compare_frames(true_reactions['convex_mentos'], reactions['convex_mentos'])
Out[76]:
In [87]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions['convex_mentos']['Net flux'].to_dict()
metabolite_data = metab['convex_mentos']['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[87]:
Where
Decision variables are emphasized in bold.
In [ ]:
np.minimum(0, driving_force)
In [89]:
import mentos
from IPython.display import HTML
obj = 'maxent_production'
def make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U):
"""Maximum entropy of ABC model"""
def maxent_production(x):
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = entropy(np.append(forward_probability, backward_probability ) )
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return maxent_production
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
true_metab[obj] = pd.read_csv('ABC_metab_{obj}.csv'.format(obj=obj), index_col=0).T
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
true_reactions[obj] = pd.read_csv('ABC_reactions_{obj}.csv'.format(obj=obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[89]:
In [65]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[65]:
Where
Decision variables are emphasized in bold.
In [91]:
from IPython.display import HTML
obj = 'minent_production'
def make_minent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U):
"""Minimum entropy of ABC model"""
def minent_production(x):
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = entropy(np.append(forward_probability, backward_probability ) )
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return f, g, fail
return minent_production
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj].astype(np.float64),metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[91]:
In [92]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[92]:
Where
\begin{itemize} \item $\vec{\bf r}_+,\vec{\bf r}_-$ are decision variables representing the forward and backward rates, respectively. \item ${\mathscr P}_{+i},{\mathscr P}_{-i}$ are the normalized forward and backward thermodynamic driving forces $\frac{r_{+i}}{r_{-i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ and $\frac{r_{-i}}{r_{+i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ \item $\vec{\bf c}$ is a decision variable representing the chemical concentrations of each metabolite \item $S$ is the $m\times n$ stoichiometric matrix of representing $m$ metabolites and $n$ reactions of the model \item $\vec{\mu}^0$ is the vector of standard chemical potentials \item $L_+ = \frac{r_+}{r_-}$ \item $L_- = \frac{r_-}{r_+}$ \item $v_{lower}(L) =\min\left(-1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) + 1$ \item $v_{upper}(L) = \max\left(1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) - 1$ \end{itemize}Decision variables are emphasized in bold.
In [93]:
obj = 'micro_maxent_production_net_flux_diff'
def make_micro_maxent_production_net_flux_diff_objective(fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U):
def micro_maxent_production_net_flux_diff(x):
"""Maximize product of microscopic entropy production and net_flux of the ABC model"""
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = -np.sum( np.log(forward_probability)*forward_probability*net_flux
- np.log(backward_probability)*backward_probability*net_flux)
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return micro_maxent_production_net_flux_diff
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
#metab[obj].to_csv('ABC_metab_{}.csv'.format(obj))
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
#display(metab[obj].T)
#HTML(reactions[obj].T.to_html())
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0).astype(np.float64)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[93]:
In [22]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[22]:
Where
\begin{itemize} \item $\vec{\bf r}_+,\vec{\bf r}_-$ are decision variables representing the forward and backward rates, respectively. \item ${\mathscr P}_{+i},{\mathscr P}_{-i}$ are the normalized forward and backward thermodynamic driving forces $\frac{r_{+i}}{r_{-i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ and $\frac{r_{-i}}{r_{+i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ \item $\vec{\bf c}$ is a decision variable representing the chemical concentrations of each metabolite \item $S$ is the $m\times n$ stoichiometric matrix of representing $m$ metabolites and $n$ reactions of the model \item $\vec{\mu}^0$ is the vector of standard chemical potentials \item $L_+ = \frac{r_+}{r_-}$ \item $L_- = \frac{r_-}{r_+}$ \item $v_{lower}(L) =\min\left(-1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) + 1$ \item $v_{upper}(L) = \max\left(1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) - 1$ \end{itemize}Decision variables are emphasized in bold.
In [94]:
obj = 'micro_maxent_production_net_flux'
def make_micro_maxent_production_net_flux_objective(fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U):
def micro_maxent_production_net_flux(x):
"""Maximize product of microscopic entropy production and net_flux of the ABC model"""
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = -np.sum( np.log(forward_probability)*forward_probability*net_flux
+ np.log(backward_probability)*backward_probability*net_flux)
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return micro_maxent_production_net_flux
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
#metab[obj].to_csv('ABC_metab_{}.csv'.format(obj))
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
#display(metab[obj].T)
#HTML(reactions[obj].T.to_html())
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0).astype(np.float64)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[94]:
In [71]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[71]:
Where
\begin{itemize} \item $\vec{\bf r}_+,\vec{\bf r}_-$ are decision variables representing the forward and backward rates, respectively. \item ${\mathscr P}_{+i},{\mathscr P}_{-i}$ are the normalized forward and backward thermodynamic driving forces $\frac{r_{+i}}{r_{-i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ and $\frac{r_{-i}}{r_{+i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ \item $\vec{\bf c}$ is a decision variable representing the chemical concentrations of each metabolite \item $S$ is the $m\times n$ stoichiometric matrix of representing $m$ metabolites and $n$ reactions of the model \item $\vec{\mu}^0$ is the vector of standard chemical potentials \item $L_+ = \frac{r_+}{r_-}$ \item $L_- = \frac{r_-}{r_+}$ \item $v_{lower}(L) =\min\left(-1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) + 1$ \item $v_{upper}(L) = \max\left(1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) - 1$ \end{itemize}Decision variables are emphasized in bold.
In [95]:
obj = 'micro_maxent_production_rate'
def make_micro_maxent_production_rate_objective(fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U):
def micro_maxent_production_rate(x):
"""Maximize product of microscopic entropy production and rate of the ABC model"""
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = -np.sum( np.log(forward_probability)*forward_probability*forward_rate
+ np.log(backward_probability)*backward_probability*backward_rate)
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return micro_maxent_production_rate
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
#metab[obj].to_csv('ABC_metab_{}.csv'.format(obj))
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
#display(metab[obj].T)
#HTML(reactions[obj].T.to_html())
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0).astype(np.float64)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[95]:
In [73]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[73]:
Where
\begin{itemize} \item $\vec{\bf r}_+,\vec{\bf r}_-$ are decision variables representing the forward and backward rates, respectively. \item ${\mathscr P}_{+i},{\mathscr P}_{-i}$ are the normalized forward and backward thermodynamic driving forces $\frac{r_{+i}}{r_{-i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ and $\frac{r_{-i}}{r_{+i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ \item $\vec{\bf c}$ is a decision variable representing the chemical concentrations of each metabolite \item $S$ is the $m\times n$ stoichiometric matrix of representing $m$ metabolites and $n$ reactions of the model \item $\vec{\mu}^0$ is the vector of standard chemical potentials \item $L_+ = \frac{r_+}{r_-}$ \item $L_- = \frac{r_-}{r_+}$ \item $v_{lower}(L) =\min\left(-1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) + 1$ \item $v_{upper}(L) = \max\left(1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) - 1$ \end{itemize}Decision variables are emphasized in bold.
In [96]:
obj = 'macro_maxent_production_net_flux'
def make_macro_maxent_production_net_flux_objective(fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U, biomass):
def macro_maxent_production_net_flux(x):
"""Maximize product of macroscopic entropy production and net_flux of the ABC model"""
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = -np.sum( np.log(forward_probability)*forward_probability*net_flux[biomass]
+ np.log(backward_probability)*backward_probability*net_flux[biomass])
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return macro_maxent_production_net_flux
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
#metab[obj].to_csv('ABC_metab_{}.csv'.format(obj))
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
#display(metab[obj].T)
#HTML(reactions[obj].T.to_html())
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0).astype(np.float64)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[96]:
In [76]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[76]:
Where
\begin{itemize} \item $\vec{\bf r}_+,\vec{\bf r}_-$ are decision variables representing the forward and backward rates, respectively. \item ${\mathscr P}_{+i},{\mathscr P}_{-i}$ are the normalized forward and backward thermodynamic driving forces $\frac{r_{+i}}{r_{-i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ and $\frac{r_{-i}}{r_{+i}}\left(\sum_j\frac{r_{+j}}{r_{-j}} + \frac{r_{-j}}{r_{+j}}\right)^{-1}$ \item $\vec{\bf c}$ is a decision variable representing the chemical concentrations of each metabolite \item $S$ is the $m\times n$ stoichiometric matrix of representing $m$ metabolites and $n$ reactions of the model \item $\vec{\mu}^0$ is the vector of standard chemical potentials \item $L_+ = \frac{r_+}{r_-}$ \item $L_- = \frac{r_-}{r_+}$ \item $v_{lower}(L) =\min\left(-1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) + 1$ \item $v_{upper}(L) = \max\left(1,\text{sign}(\log L)\cdot L^{\text{sign}(\log L)}\right) - 1$ \end{itemize}Decision variables are emphasized in bold.
In [97]:
obj = 'macro_maxent_production_net_flux'
def make_macro_maxent_production_net_flux_diff_objective(fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U, biomass):
def macro_maxent_production_net_flux_diff(x):
"""Maximize product of macroscopic entropy production and net_flux difference of the ABC model"""
log_c, \
forward_rate, backward_rate, \
log_Q, log_K, \
forward_likelihood, backward_likelihood, \
forward_probability, backward_probability, mu, \
thermodynamic_driving_force, \
net_flux = mentos.make_variables( x, fullS, mu0, deltaG0, R, T )
f = -np.sum( np.log(forward_probability)*forward_probability*net_flux[biomass]
- np.log(backward_probability)*backward_probability*net_flux[biomass])
g = np.concatenate((np.dot(S,net_flux),
np.log(forward_rate) - np.log(backward_rate) + log_Q - log_K,
[ log_c[A_ext] - np.log(c_U), # Creating a concentration gradient between A, E and D,F
log_c[E_ext] - np.log(c_U),
log_c[D_ext] - np.log(c_L),
log_c[F_ext] - np.log(c_L)],
[-mu[A] + mu[D_ext] + mu[F_ext], # Energy sink constraint
-mu[B] + mu[D_ext] + mu[F_ext],
-mu[C] + mu[D_ext] + mu[F_ext],
-mu[D] + mu[D_ext] + mu[F_ext],
-mu[E] + mu[D_ext] + mu[F_ext],
-mu[F] + mu[D_ext] + mu[F_ext]],
[ mu[A] - mu[A_ext] - mu[E_ext], # Energy barrier constraint
mu[B] - mu[A_ext] - mu[E_ext],
mu[C] - mu[A_ext] - mu[E_ext],
mu[D] - mu[A_ext] - mu[E_ext],
mu[E] - mu[A_ext] - mu[E_ext],
mu[F] - mu[A_ext] - mu[E_ext]],
np.minimum(-1, thermodynamic_driving_force)+ 1 -net_flux,
net_flux - np.maximum(1, thermodynamic_driving_force ) - 1))
fail = 0
return -f, g, fail
return macro_maxent_production_net_flux_diff
metab[obj], reactions[obj] = mentos.run_mentos(fullS, S, internal_mets, deltaG0, mu0,
c_L, c_U, v_L, v_U,
initial_log_c=metab['convex_mentos']['Log Concentrations'].values,
initial_forward_rate=reactions['convex_mentos']['Forward rate'].values,
initial_backward_rate=reactions['convex_mentos']['Backward rate'].values,
obj=make_maxent_production_objective( fullS, S, mu0, deltaG0, R, T, c_L, c_U,v_L, v_U),
biomass_rxn='R_8')
#metab[obj].to_csv('ABC_metab_{}.csv'.format(obj))
#reactions[obj].to_csv('ABC_reactions_{obj}.csv'.format(obj=obj))
#display(metab[obj].T)
#HTML(reactions[obj].T.to_html())
true_metab[obj] = pd.read_csv('ABC_metab_{}.csv'.format(obj), index_col=0).astype(np.float64)
true_reactions[obj] = pd.read_csv('ABC_reactions_{}.csv'.format(obj),index_col=0)
display(mentos.compare_frames(true_metab[obj], metab[obj]))
pd.testing.assert_frame_equal(true_metab[obj],metab[obj].astype(np.float64))
pd.testing.assert_frame_equal(true_reactions[obj],reactions[obj])
HTML(mentos.compare_frames(true_reactions[obj],reactions[obj]).to_html())
Out[97]:
In [98]:
#import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions[obj]['Net flux'].to_dict()
metabolite_data = metab[obj]['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[98]:
In [113]:
import os, re
nonchar = re.compile('[^a-zA-Z0-9]')
for obj in reactions:
print(obj)
if not os.path.exists(obj):
os.makedirs(obj)
for col in reactions[obj].columns:
filename = nonchar.sub('_',col)
print("\t{}".format(filename))
reactions[obj][col].to_csv('{obj}/reaction_{col}.tsv'.format(obj=obj,col=filename), sep='\t', header=True)
In [110]:
pd.Series.to_csv?
In [112]:
import os, re
nonchar = re.compile('[^a-zA-Z0-9]')
for obj in metab:
print(obj)
if not os.path.exists(obj):
os.makedirs(obj)
for col in metab[obj].columns:
filename = nonchar.sub('_',col)
print("\t{}".format(filename))
metab[obj][col].to_csv('{obj}/metab_{col}.tsv'.format(obj=obj,col=filename), sep='\t',header=True)
In [80]:
pmetabolites = pd.Panel.from_dict({'micro_maxent_production_net_flux':metab['micro_maxent_production_net_flux'],
'macro_maxent_production_net_flux':metab['macro_maxent_production_net_flux']},orient='minor')
preactions = pd.Panel.from_dict({'micro_maxent_production_net_flux':reactions['micro_maxent_production_net_flux'],
'macro_maxent_production_net_flux':reactions['macro_maxent_production_net_flux']},orient='minor')
HTML(preactions.swapaxes(1,0).to_frame(filter_observations=False).xs(
'Net flux',level='major', drop_level=False).to_html())
Out[80]:
In [81]:
reactions = pd.Panel.from_dict(reactions,orient='minor')
HTML(preactions.swapaxes(1,0).to_frame(filter_observations=False).to_html())
Out[81]:
In [82]:
HTML(pmetabolites.swapaxes(1,0).to_frame(filter_observations=False).to_html())
Out[82]:
Note that
In [158]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions['micro_maxent_production_rate']['Net flux'].to_dict()
metabolite_data = metab['micro_maxent_production_rate']['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[158]:
In [231]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions['analytic_mentos']['Net flux'].to_dict()
metabolite_data = metab['convex_mentos']['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[231]:
In [168]:
import escher,os
abc_dir = 'Mentos/ABC'
map_name = os.path.join(abc_dir,'ABC_map.json')
model_name = os.path.join(abc_dir,'ABC_model.json')
reaction_data = reactions['macro_maxent_production_net_flux']['Delta G'].to_dict()
metabolite_data = metab['macro_maxent_production_net_flux']['Chemical potential'].to_dict()
builder = escher.Builder(map_json=map_name, model_json=model_name,
reaction_data=reaction_data,metabolite_data=metabolite_data)
builder.display_in_notebook(menu='zoom')
Out[168]:
In [242]:
pmetabolites = pd.Panel.from_dict(metabolites,orient='minor')
preactions = pd.Panel.from_dict(reactions, orient='minor')
display(preactions['Forward likelihoods'].add_prefix('L_f '))
display(preactions['Backward likelihoods'].add_prefix('L_b '))
display(preactions['Entropy'].add_prefix('Entropy '))
display(pmetabolites['Concentrations'].add_prefix('Conc '))
display(pmetabolites['Counts'].add_prefix('Count '))
In [151]:
pmetabolites = pd.Panel.from_dict(metab,orient='minor')
preactions = pd.Panel.from_dict(reactions,orient='minor')
HTML(preactions.swapaxes(1,0).to_frame(filter_observations=False).to_html())
Out[151]:
When the network is in steady state, the concentrations of the internal metabolites do not change. Therefore $$ \frac{d\vec{c}}{dt} = S\cdot\vec{v} = 0$$ where $\vec{c}$ are the concentrations of the internal metabolites, and $\vec{v}$ are the steady state fluxes, and $S$ is the stoichiometric matrix that only includes internal metabolites.
$$ S\cdot\vec{v} = \left[ {\begin{array}{cccccccccc} & R_1 & R_2 & R_3 & R_4 & R_5 & R_6 & R_7 & R_8 & R_9 \\ A & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ B & 0 & 1 & 0 & -1 & 0 & -2 & 0 & 0 & 0 \\ C & 0 & 0 & 1 & 0 & 0 & 1 & -1 & 0 & 0 \\ D & 0 & 0 & 0 & 2 & 0 & 0 & 1 & -1 & 0 \\ E & 0 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\ F & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -1 \\ \end{array}}\right]\cdot\left[ {\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \\ v_6 \\ v_7 \\ v_8 \\ v_9 \\ \end{array}}\right] = 0 $$Given a set of initial conditions, we can find the optimal growth rate by solving the following optimization problem: $$\max_v( v_{growth} )$$ $$\mbox{Subject to:}$$ $$\begin{eqnarray} & S\cdot v & = 0 \\ \|v\|_1 \leq 1\\ \end{eqnarray}$$
In [115]:
from cvxpy import *
S = fullS.loc[internal_mets].as_matrix()
m,n = S.shape
v = Variable(n)
growth = v[7]
prob = Problem(
Maximize( growth ),
[ S*v == 0,
norm1(v) <= 1])
prob.solve()
print("status: {}".format( prob.status))
print("Optimal value: {}".format( prob.value))
print("Optimal var: \n", )
reactions['Fluxes'] = pd.DataFrame(v.value, index=rxns)
display(reactions[['Net probabilities', 'Fluxes']])
reactions['Net probabilities'].sum()
Out[115]:
In [ ]: