diff --git a/xbanalysis/problems/gba_sflux_problem.py b/xbanalysis/problems/gba_sflux_problem.py index e431909bc9bf62d5176d69af7ce8fd69ae7be98e..a8f985728fb1a09086d05b41241aaa5b72405eda 100644 --- a/xbanalysis/problems/gba_sflux_problem.py +++ b/xbanalysis/problems/gba_sflux_problem.py @@ -35,46 +35,49 @@ class GbaSfluxProblem: # while XBA model contains enzymes for all metabolic reactions with related synthesis reactions, # in scaled flux GBA optimization one a single protein (here we use the ribosome) is used metab_sids = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000247')] - self.sids_e = [sid for sid in metab_sids if xba_model.species[sid].constant is True] - sids_s = [sid for sid in metab_sids if xba_model.species[sid].constant is False] - sid_p = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000250')][0] - self.sids_sp = sids_s + [sid_p] - self.sids_esp = self.sids_e + self.sids_sp - - rids_be = [r.id for r in xba_model.reactions.values() if r.sboterm.is_in_branch('SBO:0000167')] - rid_r = xba_model.species[sid_p].ps_rid - self.rids_ber = rids_be + [rid_r] - - self.mw_esp = np.array([xba_model.species[sid].mw for sid in self.sids_esp]) - self.mw_ber = np.array([xba_model.species[xba_model.reactions[rid].enzyme].mw - for rid in self.rids_ber]) - # mw_ber is used for scaling of kcats, weight for protein synthesis needs to be rescaled - s = xba_model.species[sid_p] + self.const_metab_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is True] + self.var_metab_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is False] + sid_ribo = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000250')][0] + self.var_sids = self.var_metab_sids + [sid_ribo] + self.sids = self.const_metab_sids + self.var_sids + self.sid2idx = {msid: idx for idx, msid in enumerate(self.sids)} + + self.mr_rids = [r.id for r in xba_model.reactions.values() + if r.sboterm.is_in_branch('SBO:0000167')] + rid_ps_ribo = xba_model.species[sid_ribo].ps_rid + self.rids = self.mr_rids + [rid_ps_ribo] + self.rid2idx = {mrid: idx for idx, mrid in enumerate(self.rids)} + + self.mw_sids = np.array([xba_model.species[sid].mw for sid in self.sids]) + self.mw_enz = np.array([xba_model.species[xba_model.reactions[rid].enzyme].mw + for rid in self.rids]) + # mw_enz is used for scaling of kcats, weight for protein synthesis needs to be rescaled + s = xba_model.species[sid_ribo] if hasattr(s, 'rba_process_costs'): p_len = s.rba_process_costs['translation'] else: - r = xba_model.reactions[rid_r] + r = xba_model.reactions[rid_ps_ribo] p_len = min(r.reactants.values()) - self.mw_ber[-1] *= p_len + self.mw_enz[-1] *= p_len # retrieve S-matrix: internal metabos + total protein and metabolic reaction + ribosome synthesis - sm_espxber = xba_model.get_stoic_matrix(self.sids_esp, self.rids_ber) + sids_x_rids = xba_model.get_stoic_matrix(self.sids, self.rids) # calulate mass turnover in gram reactants/products per mol reactions (g/mol) - sm_espxber_react = np.where(sm_espxber > 0, 0, sm_espxber) - sm_espxber_prod = np.where(sm_espxber < 0, 0, sm_espxber) - mass_turnover_react = -sm_espxber_react.T.dot(self.mw_esp) - mass_turnover_prod = sm_espxber_prod.T.dot(self.mw_esp) + sids_x_rids_react = np.where(sids_x_rids > 0, 0, sids_x_rids) + sids_x_rids_prod = np.where(sids_x_rids < 0, 0, sids_x_rids) + mass_turnover_react = -sids_x_rids_react.T.dot(self.mw_sids) + mass_turnover_prod = sids_x_rids_prod.T.dot(self.mw_sids) self.mass_turnover = np.max((mass_turnover_react, mass_turnover_prod), axis=0) # mass normalized stoichiometric matrix - self.smmn_espxber = (sm_espxber * self.mw_esp.reshape(-1, 1)) / self.mass_turnover + self.sidx_x_rids_mn = (sids_x_rids * self.mw_sids.reshape(-1, 1)) / self.mass_turnover self.rho = xba_model.parameters['rho'].value - self.cemn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw - for sid in self.sids_e]) - self.cimn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw - for sid in sids_s]) + self.ext_conc_mn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw + for sid in self.const_metab_sids]) + self.var_initial_mn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw + for sid in self.var_metab_sids]) # identify unit ids for kinetic parameters units_per_s = [XbaUnit({'kind': 'second', 'exp': -1.0, 'scale': 0, 'mult': 1.0})] @@ -163,37 +166,35 @@ class GbaSfluxProblem: :return: gba model mass normalized as per Hugo Dourado :rtype: dict of pandas DataFrames """ - protein_mn = self.rho - sum(self.cimn) + protein_mn = self.rho - sum(self.var_initial_mn) - kcats = np.array([self.get_kcats(rid) for rid in self.rids_ber]).T - skcats = kcats * 3600.0 * self.mass_turnover / self.mw_ber + kcats = np.array([self.get_kcats(rid) for rid in self.rids]).T + skcats = kcats * 3600.0 * self.mass_turnover / self.mw_enz - km = np.zeros_like(self.smmn_espxber) - ki = np.zeros_like(self.smmn_espxber) - rid2idx = {mrid: idx for idx, mrid in enumerate(self.rids_ber)} - sid2idx = {msid: idx for idx, msid in enumerate(self.sids_esp)} - for rid in self.rids_ber: + km = np.zeros_like(self.sidx_x_rids_mn) + ki = np.zeros_like(self.sidx_x_rids_mn) + for rid in self.rids: kms, kis = self.get_kms(rid) for sid, val in kms.items(): - km[sid2idx[sid], rid2idx[rid]] = val + km[self.sid2idx[sid], self.rid2idx[rid]] = val for sid, val in kis.items(): - ki[sid2idx[sid], rid2idx[rid]] = val + ki[self.sid2idx[sid], self.rid2idx[rid]] = val # scale Michaelis constants with molecluar weights - skm = km * self.mw_esp.reshape(-1, 1) - ski = ki * self.mw_esp.reshape(-1, 1) - ska = np.zeros_like(self.smmn_espxber) + skm = km * self.mw_sids.reshape(-1, 1) + ski = ki * self.mw_sids.reshape(-1, 1) + ska = np.zeros_like(self.sidx_x_rids_mn) mn_model = { - 'N': pd.DataFrame(self.smmn_espxber, index=self.sids_esp, columns=self.rids_ber), - 'KM': pd.DataFrame(skm, index=self.sids_esp, columns=self.rids_ber), - 'KI': pd.DataFrame(ski, index=self.sids_esp, columns=self.rids_ber), - 'KA': pd.DataFrame(ska, index=self.sids_esp, columns=self.rids_ber), - 'kcat': pd.DataFrame(skcats, index=['kcat_f', 'kcat_b'], columns=self.rids_ber), - 'conditions': pd.DataFrame(np.hstack((self.rho, self.cemn)), - index=['rho'] + self.sids_e, columns=[1]), - 'lower_c': pd.DataFrame(np.zeros(len(self.sids_sp)), - index=self.sids_sp, columns=['lower']), - 'initial': pd.DataFrame(np.hstack((self.cimn, protein_mn)), - index=self.sids_sp, columns=['initial']), + 'N': pd.DataFrame(self.sidx_x_rids_mn, index=self.sids, columns=self.rids), + 'KM': pd.DataFrame(skm, index=self.sids, columns=self.rids), + 'KI': pd.DataFrame(ski, index=self.sids, columns=self.rids), + 'KA': pd.DataFrame(ska, index=self.sids, columns=self.rids), + 'kcat': pd.DataFrame(skcats, index=['kcat_f', 'kcat_b'], columns=self.rids), + 'conditions': pd.DataFrame(np.hstack((self.rho, self.ext_conc_mn)), + index=['rho'] + self.const_metab_sids, columns=[1]), + 'lower_c': pd.DataFrame(np.zeros(len(self.var_sids)), + index=self.var_sids, columns=['lower']), + 'initial': pd.DataFrame(np.hstack((self.var_initial_mn, protein_mn)), + index=self.var_sids, columns=['initial']), } return mn_model