gba_sflux_problem variable names aligned with gba_problem

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