1. Deriving Alpha Table from FBA

import os
import pandas as pd

from scarcc.preparation.metabolic_model import BasicModel
from scarcc.preparation.find_directory import find_directory
from scarcc.preparation.target_gene.gene_combination_handler import get_DG_list, get_SG_list

# get file directory
model_directory = find_directory('models', os.path.abspath(''))
data_directory = find_directory('Data', os.path.abspath(''))

# initialize model
E0, S0, _ = BasicModel(model_directory=model_directory, flux_weighting=True).load_ES_models()
DG_list = get_DG_list(os.path.join(data_directory, 'GeneCombos.csv'), n_combos=None)
SG_list = get_SG_list(DG_list)

1.1. Search for \(\alpha\) that corresponds to IC50

Class MonocultureAlphaFinder is responsible for finding the alpha values the corresponds to the IC50 dosage(default) of antibiotics, given a target gene and the metabolic model. Which is based on the reduction in biomaffss flux as the objective function to be evaluated.

Steps for deriving the appropriate alpha:

1. Assign medium for each species

2. Program will iteratively search for the \(\alpha\) such that the normalized biomass flux is within threshold around target_normalized_biomass = 0.5

3. Conditions for stop searching

   • Convergence of biomass flux around target_normalized_biomass

   • Reached desired \(\alpha\) precision

   • Maximum iteration

NOTE : Non-essential genes are assigned with values 100000, which ensures no flux will went through that reaction even if it is non-essential, this measures allows for synergism through blockage of pathway and change in metabolic flux pattern

Minimum arguments for construction of MonocultureAlphaFinder: model and current_gene

They are automatically handled in get_alpha_biomass_df when given model_list and desired gene lists, see section “Production of alpha_table”

from scarcc.preparation.perturbation.alpha_finder.monoculture import MonocultureAlphaFinder
ma = MonocultureAlphaFinder(model=E0, current_gene='folA')
alpha, norm_bm, summary = ma.find_feasible_alpha()
print('alpha:', alpha, 'Normalized biomass', norm_bm, 'Flux summary', summary.T)

alpha: 0    62.28375
Name: E0_alpha, dtype: float64 Normalized biomass 0    0.500011
Name: E0_Normalized_biomass, dtype: float64                                                                          0
E0_alpha                                                          62.28375
E0_biomass                                                         0.21007
E0_Normalized_biomass                                             0.500011
E0_forward_reactions_id                                      DHFR, DHFR_v1
E0_forward_reactions     DHFR: dhf_c + h_c + nadph_c --> 0.016055552210...
E0_forward_flux_values                                           [0.35111]
E0_opposite_reactions    DHFR_v1: 0.016055552210648805 dhf_c + 0.016055...
E0_opposite_flux_values                                              [0.0]
E0_Net_Flux_Boolean                                               Net Flux
E0_Net_Flux                                                        0.35111
E0_is_growth_switch                                                  False
Gene_inhibition                                                       folA

Arguments Customization

maf_kwargs = {
    'model': E0,
    'current_gene': 'folA', # target_gene
    'target_normalized_biomass': 0.5, # desired normalized biomass flux
    'exp_leap': 2, # raise alpha to the power of exp_leap in next search step
    'precision': 2, # precision of optimal alpha, stop search when reaching this precision
    'acceptance_threshold_lower': 1, # lower bound of desired normalized biomass flux: target_normalized_biomass*acceptance_threshold_lower
    'acceptance_threshold_upper': 1.1, # upper bound of desired normalized biomass flux: target_normalized_biomass*acceptance_threshold_upper
    'search_alpha': 1.02, # initial alpha to start search
    'iter_max': 25 # maximum number of iterations before stopping search
    }
maf = MonocultureAlphaFinder(**maf_kwargs)
alpha, norm_bm, summary = maf.find_feasible_alpha()

1.2. Production of alpha_table

Integration of results with given list of models and potential target genes provided, get_alpha_biomass_df automatically pass each to MonocultureAlphaFinder and produce a complete alpha_table_m1

Note : Keys model and current_genes should not be included in maf_kwargs

from scarcc.preparation.perturbation.alpha_finder.monoculture import get_alpha_biomass_df
maf_kwargs = {
            'target_normalized_biomass': 0.5,
            'potential_genes': SG_list[:5],
            'precision': 2,
            'acceptance_threshold_lower': 1}
# default detailed_alpha_table=False
alpha_table = get_alpha_biomass_df(model_list = [E0, S0], data_directory=data_directory, **maf_kwargs)
alpha_table.head()

alpha_table_m1.csv is saved in d:\scarccpy\Data

	E0_alpha	E0_Normalized_biomass	S0_alpha	S0_Normalized_biomass	E0_response	S0_response
Gene_inhibition
folA	62.283750	0.500011	19.698750	0.500092	Essential or Phenotype	Essential or Phenotype
argD	11.793750	0.499911	8.494688	0.499914	Essential or Phenotype	Essential or Phenotype
tktA	27.540000	0.499759	17.276250	0.500377	Essential or Phenotype	Essential or Phenotype
pheA	6.916875	0.500100	6.371016	0.500089	Essential or Phenotype	Essential or Phenotype
rffG	100000.000000	1.000000	119.658750	0.500014	Non_essential	Essential or Phenotype

rffG is non_essential and alpha is set to 100000

df_detailed = get_alpha_biomass_df(model_list=[E0, S0], potential_genes=['gltD'], detailed_alpha_table=False, **maf_kwargs)

	E0_alpha	E0_Normalized_biomass	S0_alpha	S0_Normalized_biomass	E0_response	S0_response
Gene_inhibition
gltD	1.1475	0.500403	100000.0	1.0	Essential or Phenotype	Non_essential

alpha_table with only alpha and the normalized biomass

1.3. Passing detailed dataframe as argument to retrieve alpha_table

alpha_table2 = get_alpha_biomass_df(alpha_biomass_df=df_detailed,
                                    detailed_alpha_table=False)
alpha_table2.head()

	E0_alpha	E0_Normalized_biomass	S0_alpha	S0_Normalized_biomass	E0_response	S0_response	E0_response	S0_response
Gene_inhibition
gltD	1.1475	0.500403	100000.0	1.0	Essential or Phenotype	Non_essential	Essential or Phenotype	Non_essential

1.4. Essentiality difference between models

Some encoded reactions corresponding to the target gene are nonessential for one species while being essential for the other

1. Flux dropped to 0 even without reduction in biomass yield

2. Reduction in flux through reactions lead to reduction in biomass production).

alpha_table2.loc[['gltD']]

	E0_alpha	E0_Normalized_biomass	S0_alpha	S0_Normalized_biomass	E0_response	S0_response	E0_response	S0_response
Gene_inhibition
gltD	1.1475	0.500403	100000.0	1.0	Essential or Phenotype	Non_essential	Essential or Phenotype	Non_essential

1.5. Heterogeneity of \(\alpha\) across medium property

Value of \(\alpha\) is dependent on the type and concentration of nutrients and carbon source(lactose vs methionine limiting condition), which also affects the normalization denominator under each environmental condition.

Note : Default medium is set in BasicModel

E.coli medium under limiting lactose condition

print(E0.medium)
print('biomass flux =', round(E0.slim_optimize()), 'under EX_lcts_e uptake limit: =', E0.medium['EX_lcts_e'],)

{'EX_pi_e': 10, 'EX_met__L_e': 0.0647, 'EX_fe3_e': 10, 'EX_mn2_e': 10, 'EX_fe2_e': 10, 'EX_zn2_e': 10, 'EX_mg2_e': 10, 'EX_lcts_e': 10, 'EX_ca2_e': 10, 'EX_ni2_e': 10, 'EX_cu2_e': 10, 'EX_cobalt2_e': 10, 'EX_mobd_e': 10, 'EX_so4_e': 10, 'EX_nh4_e': 10, 'EX_k_e': 10, 'EX_cl_e': 10, 'EX_o2_e': 10}
biomass flux = 0 under EX_lcts_e uptake limit: = 10

Try setting medium with limited methionine and unlimited lactose

E0.medium = {
    'EX_pi_e': 10,
    'EX_met__L_e': 0.0647,
    # 'EX_met__L_e': 1,
    'EX_fe3_e': 10,
    'EX_mn2_e': 10,
    'EX_fe2_e': 10,
    'EX_zn2_e': 10,
    'EX_mg2_e': 10,
    'EX_lcts_e': 10,
    'EX_ca2_e': 10,
    'EX_ni2_e': 10,
    'EX_cu2_e': 10,
    'EX_cobalt2_e': 10,
    'EX_mobd_e': 10,
    'EX_so4_e': 10,
    'EX_nh4_e': 10,
    'EX_k_e': 10,
    'EX_cl_e': 10,
    'EX_o2_e': 10}
print('biomass flux =', round(E0.slim_optimize(),3), 'under EX_met__L_e uptake limit: =', E0.medium['EX_met__L_e'],)

biomass flux = 0.42 under EX_met__L_e uptake limit: = 0.0647

Warning : Medium is a dictionary assigned to the metabolic model, changing medium inside context manager WILL NOT revert the assignment of medium

exaggerated_alpha = get_alpha_biomass_df(model_list=[E0], potential_genes=['folA'], detailed_alpha_table=False, **maf_kwargs)
# df_detailed = get_alpha_biomass_df(model_list=[E0], potential_genes=['folA'], detailed_alpha_table=True, **maf_kwargs)

# exaggerated_alpha.columns
position = 'folA', 'E0_alpha'

print('Methionie Limiting IC50 alpha: ', exaggerated_alpha.loc[position])
print('Lactose Limiting IC50 alpha:', df_detailed.loc[position])

Methionie Limiting IC50 alpha:  99.32249999999999
Lactose Limiting IC50 alpha: 52.785000000000004