4. Stoichiometry Scaling

import pandas as pd
import os
from scarcc.preparation.find_directory import find_directory
from scarcc.preparation.metabolic_model import BasicModel

data_directory = find_directory('Data', os.path.abspath(''))
model_directory = find_directory('models', os.path.abspath(''))
sim_chamber_directory = find_directory('SimChamber', os.path.abspath(''))

E0, S0, all_components = BasicModel(flux_weighting=True, model_directory=model_directory).load_ES_models()
alpha_table = pd.read_csv(os.path.join(data_directory, 'alpha_table_m1.csv'), index_col=0)

This implementation reduce the enzymatic activities controlled by the target gene by \(\alpha\), regardless of duplicated genes controlling the reactions

• Irreversible reaction:

\begin{align} A \rightarrow B \Rightarrow A \rightarrow 1/\alpha \times B \end{align}

• Reversible reaction:

\begin{align} A \leftrightarrow B \Rightarrow \begin{cases} A \rightarrow \frac{1}{\alpha} \times B\\ \frac{1}{\alpha} \times A \leftarrow B \end{cases} \end{align}

Note 1: Perform scaling within context manager to recover to normal state for new set of gene perturbation

Note 2: This implementation is because of the discrepancy between gene knock out and gene inhibition. In gene knock out, the reaction will remain unaffected in case of gene redundancy, whereas to quantify the inhibitory effect of one antibiotics(e.g. IC50), the gene redundancy is ignored for simplicity

alpha_table = pd.read_csv(os.path.join(data_directory, 'alpha_table_m1.csv'), index_col=0)

from scarcc.preparation.metabolic_model.core.component import get_gene_id
from scarcc.preparation.perturbation import alter_Sij
current_gene = 'folA'
with E0, S0:
    alter_Sij(E0, alphas=alpha_table.loc[current_gene, 'E0'], genes=current_gene, ko=False)
    alter_Sij(S0, alphas=alpha_table.loc[current_gene, 'S0'], genes=current_gene, ko=False)

    print('After: ')
    print('E0: ', [str(rxn) for rxn in E0.genes.get_by_id(get_gene_id(E0, 'folA')).reactions])
    print('S0: ', [str(rxn) for rxn in S0.genes.get_by_id(get_gene_id(S0, 'folA')).reactions])
print('Before: ')
print('E0: ', [str(rxn) for rxn in E0.genes.get_by_id(get_gene_id(E0, 'folA')).reactions])
print('S0: ', [str(rxn) for rxn in S0.genes.get_by_id(get_gene_id(S0, 'folA')).reactions])

After:
E0:  ['DHFR_v1: 0.01896481730065458 dhf_c + 0.01896481730065458 h_c + 0.01896481730065458 nadph_c <-- nadp_c + thf_c', 'DHFR: dhf_c + h_c + nadph_c --> 0.01896481730065458 nadp_c + 0.01896481730065458 thf_c']
S0:  ['DHFR: dhf_c + h_c + nadph_c --> 0.05612262794205336 nadp_c + 0.05612262794205336 thf_c', 'DHFR_v1: 0.05612262794205336 dhf_c + 0.05612262794205336 h_c + 0.05612262794205336 nadph_c <-- nadp_c + thf_c']
Before:
E0:  ['DHFR: dhf_c + h_c + nadph_c <=> nadp_c + thf_c']
S0:  ['DHFR: dhf_c + h_c + nadph_c <=> nadp_c + thf_c']

If both target gene in gene combination leads to the inhibition of the same reaction, the inhibition effect to that particular reaction is set to the maximum one rather than scaled twice

Example:

both sucA and sucB gene encode for AKGDH

• Case 1:

  • sucA with \(\alpha\) = 10, sucB with \(\alpha\) = 50

• Case 2:

  • sucA with \(\alpha\) = 50, sucB with \(\alpha\) = 50

with E0, S0:
    alter_Sij(E0, alphas=[10,50], genes=['sucA', 'sucB'], ko=False)

    message1 = 'E0: '+ str([str(rxn) for rxn in E0.genes.get_by_id(get_gene_id(E0, 'sucA')).reactions])
    print('Case 1: \n', message1)

with E0, S0:
    alter_Sij(E0, alphas=[50,50], genes=['sucA', 'sucB'], ko=False)

    message2 = 'E0: '+ str([str(rxn) for rxn in E0.genes.get_by_id(get_gene_id(E0, 'sucA')).reactions])
    print('Case 2: \n', message2)

if message1 == message2:
    print(f'Output of Case 1 and Case 2 are the same')

Case 1:
 E0: ['AKGDH: akg_c + coa_c + nad_c --> 0.020000000000000018 co2_c + 0.020000000000000018 nadh_c + 0.020000000000000018 succoa_c']
Case 2:
 E0: ['AKGDH: akg_c + coa_c + nad_c --> 0.020000000000000018 co2_c + 0.020000000000000018 nadh_c + 0.020000000000000018 succoa_c']
Output of Case 1 and Case 2 are the same