<jupyter_start><jupyter_text>_E. coli_ core model with some changes:	* Corrected the transhydrogenase (THD2) to one proton translocation;* Deactivated pyruvate formate lyase (PFL) under aerobic coditions.<jupyter_code>model = cobra.io.load_json_model(r'..\0_ecoli_models\e_coli_core.json')

model.reactions.THD2.subtract_metabolites({"h_e": -1, "h_c": 1})
KORxn_base = ['PFL']
KORxn(model, KORxn_base)
model.reactions.EX_glc__D_e.bounds = (0,0)

model_wt = model.copy()<jupyter_output><empty_output><jupyter_text>Model curations: Anticipated results of a model knocked out the _lpd_ gene:  * no growth on acetate only;  * grow on acetate + formate (with formate dehydrogenase), and where formate serves only as electron source<jupyter_code># knock out lpdA
model.genes.b0116.knock_out()
# on acetate only
with model as m: 
    m.reactions.EX_ac_e.lower_bound = -10  # arbitrary 
    cobra.flux_analysis.pfba(m)
    m.summary()
    flux2file(m,'biomass_1','lpd_core')<jupyter_output>IN FLUXES      OUT FLUXES    OBJECTIVES
-------------  ------------  ----------------------
o2_e   15      h2o_e  15.9   BIOMASS_Ecol...  0.115
ac_e   10      co2_e  15.1
h_e     7.69
nh4_e   0.627
pi_e    0.423<jupyter_text>The model uses the oxidative pentose phosphate pathway generating NADPH, which is then used by NAD transhydrogenase (NADTRHD, sthA) transferring the electrons.  Try to block it by KO GND.<jupyter_code>with model_wt as m:
    m.reactions.EX_ac_e.lower_bound = -10  
    cobra.flux_analysis.pfba(m)
    m.summary()
    
    print('\n','* - * - *')
    m.reactions.GND.knock_out()
    cobra.flux_analysis.pfba(m)
    m.summary()
    
# KO GND doesn't change the results on WT, so block it in the model.
with model as m: 
    m.reactions.EX_ac_e.lower_bound = -10
    KORxn(m, ['GND'])
    cobra.flux_analysis.pfba(m)
    m.summary()
    flux2file(m,'biomass_2','lpd_core')

# becoming infeasible
# not growing on acetate now, try to add formate
# adding FDH reaction to the model
model_FDH = model.copy()
FDH = cobra.Reaction(id="FDH", name="formate dehydrogenase", lower_bound=0, upper_bound=1000) 
model_FDH.add_reaction(FDH)
FDH.add_metabolites({"for_c": -1, "nad_c": -1, "co2_c": 1, "nadh_c": 1})

with model_FDH as m: 
    m.reactions.EX_ac_e.lower_bound = -10
    m.reactions.EX_for_e.lower_bound = -100
    KORxn(m, ['GND'])
    cobra.flux_analysis.pfba(m)
    m.summary()
    flux2file(m,'biomass_for_1','lpd_core')<jupyter_output>IN FLUXES     OUT FLUXES    OBJECTIVES
------------  ------------  ----------------------
h_e    37.6   h2o_e  41.6   BIOMASS_Ecol...  0.391
for_e  35.4   co2_e  38.8
o2_e   20.6
ac_e   10
nh4_e   2.13
pi_e    1.44<jupyter_text>The `model` now cannot use acetate as a sole carbon source.  Formate (with FDH) can be used as electron source. And it serves only as electron donor.<jupyter_code>KORxn(model_wt,['GND'])
KORxn(model, ['GND'])
KORxn(model_FDH,['GND'])<jupyter_output><empty_output><jupyter_text>Calculating the substrate uptake rates: WT grown on aceate has doubling time **2.9 h**, i.e., **0.24 h-1**.<jupyter_code>DT = Q_(2.9,'hr')   # WT grow on 20 mM acetate -> DT = 2.9 h
gr = np.log(2) / DT  
print(f'growth rate is {gr: ~.2f}')
with model_wt as m:
    m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
    m.reactions.EX_glc__D_e.bounds = (0,0)
    m.reactions.EX_ac_e.bounds = (-1000,0)
    m.objective = {m.reactions.EX_ac_e: -1}
    m.objective_direction = 'min'
    cobra.flux_analysis.pfba(m)
    m.summary()

    AcUp = Q_(m.objective.value,'mmol/gCDW/hr')
    print(f'acetate uptake rate is {AcUp:.2f~P}')
# lpd on glucose + acetate
with model as m:
    m.reactions.EX_glc__D_e.lower_bound = -10.5  # in mmol/gCDW/hr
    m.reactions.EX_ac_e.lower_bound = - AcUp.magnitude        # in mmol/gCDW/hr
    cobra.flux_analysis.pfba(m)
    m.summary()

    sol = Q_(m.objective.value,'1/hr')

print('\n', f'The doubling time is {(np.log(2) / sol): ~.2f}')<jupyter_output>IN FLUXES        OUT FLUXES    OBJECTIVES
---------------  ------------  ----------------------
o2_e      15     h2o_e  29.7   BIOMASS_Ecol...  0.446
ac_e      13     pyr_e  21.8
glc__D_e  10.5   h_e    17.7
nh4_e      2.43  co2_e   4.65
pi_e       1.64

 The doubling time is 1.56 hr<jupyter_text>The maximum acetate uptake rate is **13.0 mmol/gCDW/h**, which is within the range of 11 - 15.5, reported in https://doi.org/10.1038/84379.   Also, on glc + Ac, the doubling time is colse to the experiment value *1.6 h* (Fig. 2). Here the glc uptake rate is set to 10.5 according to [<NAME>, 1994](https://pubmed.ncbi.nlm.nih.gov/7986045/).Therefore, the model is good: _lpd_ (b0116, affecting PDH & AKGDH), GND and PFL are knocked out. The yield can be calculated as the ratio between the growth rate and the substrate uptake/utilization rate:   $$\dfrac{1/h}{mmol/gCDW/h} = \dfrac{gCDW}{mmol} = \dfrac{gCDW/L}{mmol/L}$$Using BionNumber [109836](https://bionumbers.hms.harvard.edu/bionumber.aspx?s=n&v=4&id=109836), one unit of OD600 _E. coli_ corresponds to a cell dry weight of _0.39 gCDW/L_, biomass yield can be converted to OD600 (increase).  **Assuming all the limiting substrates are consumed in the experiments**<jupyter_code>OD2Y = Q_(0.39,'gCDW/L')  # per unit of OD600 =1 culture<jupyter_output><empty_output><jupyter_text>_$\Delta$lpd_ growing on acetate + formate:  For the _$\Delta$lpd_ strain grows on 20 mM acetate and formate:  The growth rates of the $\Delta$_lpd_ strain on acetate and formate are lower than WT on acetate. Especially when formate concentration below 13 mM, the growth rate increased along with the formate concentration increase, this could indicate that formate uptake (rate) or formate oxidation is the limit under low concentrations.  Also, the maximum ODs are all below the WT on acetate alone. The maximum ODs increased along with the increase of formate concentration. These indicates that acetate is not the limit but formate (energy source). For the $\Delta$_lpd_ strain on acetate and formate, formate is the major energy source, acetate serves mainly as carbon source, the strain should have higher biomass yield than WT on soly acetate, where acetate is both carbon and energy sources.<jupyter_code># Block all exchange reactions to avoid secretion
for rxn in model_FDH.exchanges:
    if "C" in rxn.check_mass_balance(): # knock out all other carbon-related transporters
        rxn.bounds = (0, 0)
model_FDH.reactions.EX_co2_e.bounds = (-1000, 1000)

with model_FDH as m: 
    gr = np.log(2) / Q_(5.2,'hr')  # 20 mM ac + 44 mM formate, the highest growth rate
    m.reactions.EX_ac_e.lower_bound = - AcUp.magnitude
    m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
    m.reactions.EX_for_e.lower_bound = -1000
    m.objective = {m.reactions.EX_for_e: -1}
    m.objective_direction = 'min'
    cobra.flux_analysis.pfba(m)
    m.summary()
    for_up = Q_(m.objective.value,'mmol/gCDW/hr')
    flux2file(m,'biomass_for_4','lpd_core')
    
print(f'Yield is {(gr/for_up/OD2Y).magnitude:.1e} unit OD600 per mM formate')<jupyter_output>IN FLUXES      OUT FLUXES    OBJECTIVES
-------------  ------------  --------------
h_e    17.2    h2o_e  18.6   EX_for_e  16.5
for_e  16.5    co2_e  17.6
o2_e    9.23
ac_e    3.41
nh4_e   0.727
pi_e    0.49
Yield is 2.1e-02 unit OD600 per mM formate<jupyter_text>Formate is the limiting as shown in the experimental results.  Here acetate uptake rate is _3.41 mmol/gCDW/h_, (the maximum) formate uptake rate (or FDH) is calculated to be **16.5 mmol/gCDW/h**. * The ratio of formate/acetate is **4.8**.  * Biomass yield is **$ 2 \times{10 ^ {-2}}$ unit OD600 per mM formate**.<jupyter_code>growth_for = pd.DataFrame({
    'formate_conc': PA_(np.logspace(0,5,6,base=1/1.5)*100, dtype='mM'),
    'max_OD600': PA_([0.43,0.363,0.273,0.211,0.163,0.122],dtype=''),
    'doubling_time': PA_([6.1,5.6,5.2,5.2,5.5,5.8], dtype='hr')
})

growth_for['Yield'] = growth_for['max_OD600'] / growth_for['formate_conc']    # unit OD600 per mM formate
# growth_for['growth_rate'] = PA_(np.full(len(growth_for),np.log(2)),dtype='') / growth_for['doubling_time']
growth_for['growth_rate'] = PA_(np.full(len(growth_for),np.log(2)),dtype='') / growth_for['doubling_time']
growth_for.pint.dequantify()<jupyter_output><empty_output><jupyter_text>Since formate (uptake or oxidation) is the limit for growth, it determines the growth rate.<jupyter_code>formate = []
ac = []

with model_FDH as m: 
    m.reactions.EX_ac_e.lower_bound = - AcUp.magnitude  # uptake rate, in mmol/gCDW/h
    m.reactions.EX_for_e.lower_bound = -1000
    m.objective = {m.reactions.EX_for_e: -1}
    m.objective_direction = 'min'
    for i in range(len(growth_for['growth_rate'])):
        gr = growth_for.loc[i, 'growth_rate']
        m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
        cobra.flux_analysis.pfba(m)
        formate.append(abs(m.reactions.EX_for_e.flux))
        ac.append(abs(m.reactions.EX_ac_e.flux))

growth_for['EX_formate'] = PA_(formate, dtype='mmol/gCDW/hr')
growth_for['EX_acetate'] = PA_(ac, dtype='mmol/gCDW/hr')

growth_for['Yield_cal'] = growth_for['growth_rate'] / growth_for['EX_formate'] / OD2Y  # unit OD600 per mM formate
growth_for['for/ac'] = growth_for['EX_formate'] / growth_for['EX_acetate']

growth_for.pint.dequantify()<jupyter_output><empty_output><jupyter_text>_$\Delta$lpd_ growing on acetate + methanol:<jupyter_code># building the model
model_MDH = model.copy()

meoh_c = cobra.Metabolite(id='meoh_c',
                          formula='CH4O',
                          name='methanol',
                          compartment='c')
fald_c = cobra.Metabolite(id='fald_c',
                          formula='CH2O',
                          name='Formaldehyde',
                          compartment='c')

EX_meoh_e = cobra.Reaction(id="EX_meoh_e", name="methanol exchange", lower_bound=0, upper_bound=1000) 
model_MDH.add_reaction(EX_meoh_e)
EX_meoh_e.add_metabolites({meoh_c: -1})

MeDH = cobra.Reaction(id="MeDH", name="methanol dehydrogenase", lower_bound=0, upper_bound=1000) 
model_MDH.add_reaction(MeDH)
MeDH.add_metabolites({meoh_c: -1, "nad_c": -1, fald_c: 1, "nadh_c": 1, 'h_c': 1})

FRM = cobra.Reaction(id="FRM", name="formaldehyde oxidation, overall", lower_bound=0, upper_bound=1000) 
model_MDH.add_reaction(FRM)
FRM.add_metabolites({fald_c: -1, "nad_c": -1, 'h2o_c': -1, "for_c": 1, "nadh_c": 1, 'h_c': 2})
# Block exchange reactions to avoid secretion, expect formate
for rxn in model_MDH.exchanges:
    if "C" in rxn.check_mass_balance(): # knock out all other carbon-related transporters
        rxn.bounds = (0, 0)
model_MDH.reactions.EX_co2_e.bounds = (-1000, 1000)
model_MDH.reactions.EX_for_e.bounds = (0,1000)  # no FDH overexpression

with model_MDH as m:
    gr = np.log(2) / Q_(32,'hr')  # 20 mM ac + 667 mM MeOH, the highest growth rate
    m.reactions.EX_ac_e.lower_bound = -AcUp.magnitude
    m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
    m.reactions.EX_meoh_e.lower_bound = -1000
    m.objective = {m.reactions.EX_meoh_e: -1}
    m.objective_direction = 'min'
    cobra.flux_analysis.pfba(m)
    m.summary()
    flux2file(m,'biomass_meoh_2','lpd_core')<jupyter_output>IN FLUXES      OUT FLUXES    OBJECTIVES
-------------  ------------  ---------------
o2_e   4.77    h2o_e  4.96   EX_meoh_e  4.61
ac_e   0.554   for_e  4.61
nh4_e  0.118   h_e    4.49
pi_e   0.0797  co2_e  0.186<jupyter_text>Assuming the acetate uptake rate to be the maximum as WT on acetate, _13 mmol/gCDW/h_, the maximum MeOH uptake is calculated to be **4.61 mmol/gCDW/h**. Similar to the formate.<jupyter_code>growth_meoh = pd.DataFrame({
    'MeOH_conc': PA_(np.logspace(0,3,4,base=1/1.5)*1000, dtype='mM'),
    'max_OD600': PA_([0.278,0.504,0.341,0.234],dtype=''),
    'doubling_time': PA_([38,32,37,37],dtype='hr')
})

growth_meoh['Yield'] = growth_meoh['max_OD600'] / growth_meoh['MeOH_conc']  # unit OD600 per mM methanol
growth_meoh['growth_rate'] = PA_(np.full(len(growth_meoh),np.log(2)),dtype='')  / growth_meoh['doubling_time']

meoh = []
ac = []
with model_MDH as m:
    m.reactions.EX_ac_e.lower_bound = -AcUp.magnitude
    m.reactions.EX_meoh_e.lower_bound = -1000
    m.objective = {m.reactions.EX_meoh_e: -1}
    m.objective_direction = 'min'
   
    for i in range(len(growth_meoh['growth_rate'])):
        gr = growth_meoh.loc[i, 'growth_rate']
        m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
        cobra.flux_analysis.pfba(m)
        flux2file(m,f'meoh_{i}','lpd_for')
        meoh.append(abs(m.reactions.EX_meoh_e.flux))
        ac.append(abs(m.reactions.EX_ac_e.flux))

growth_meoh['EX_meoh'] = PA_(meoh, dtype='mmol/gCDW/hr')
growth_meoh['EX_acetate'] = PA_(ac, dtype='mmol/gCDW/hr')

growth_meoh['Yield_cal'] = growth_meoh['growth_rate'] / growth_meoh['EX_meoh'] / OD2Y  # unit OD600 per mM methanol
growth_meoh['meoh/ac'] = growth_meoh['EX_meoh'] / growth_meoh['EX_acetate']

growth_meoh.pint.dequantify()<jupyter_output><empty_output><jupyter_text>The gC1M-gC2M strain growing on glucose + methanol<jupyter_code># growth data without formate
growth_MDH = pd.DataFrame({
    'doubling_time': PA_([2.7,3,3,10],dtype='hr')
},index=['CnMDH', 'BmMDH*','BsMDH','CgMDH'])

growth_MDH['growth_rate'] = PA_(np.full(len(growth_MDH),np.log(2)),dtype='') / growth_MDH['doubling_time']
growth_MDH.pint.dequantify()
# growth data with foramte
growth_MDH_for = pd.DataFrame({
    'doubling_time': PA_([2.2,1.7,2.1,1.6],dtype='hr')
},index=['CnMDH', 'BmMDH*','BsMDH','CgMDH'])

growth_MDH_for['growth_rate'] = PA_(np.full(len(growth_MDH_for),np.log(2)),dtype='') / growth_MDH_for['doubling_time']  # 1/h
growth_MDH_for.pint.dequantify()
# build the model
model_C1 = model_wt.copy()
full_model = cobra.io.load_json_model(r'..\0_ecoli_models\iML1515.json')

newRxns = [
    model_MDH.reactions.EX_meoh_e,
    model_MDH.reactions.MeDH,
    model_MDH.reactions.FRM,
    full_model.reactions.FTHFLi,  # formate-THF ligase
    full_model.reactions.MTHFC,   # Methenyltetrahydrofolate cyclohydrolase
    full_model.reactions.MTHFD,   # Methylenetetrahydrofolate dehydrogenase (NADP)
    full_model.reactions.GLYCL,   # GCV
    full_model.reactions.GHMT2r   # GlyA
]

# to check the reactions
for i in newRxns:
    print(i)
newRxns[6].bounds = (-1000,1000)  # change GCV to be reversible
model_C1.add_reactions(newRxns)   # add the new reactions to the model

# update the biomass function by replacing 3pg to ser, gly and C1
# from DOI 10.1021/acssynbio.8b00093
model_C1.reactions.BIOMASS_Ecoli_core_w_GAM.build_reaction_from_string("0.462 ser__L_c + 1.033 gly_c + 1.033 mlthf_c + 3.7478 accoa_c + 59.81 atp_c + 0.361 e4p_c + 0.0709 f6p_c + 0.129 g3p_c + 0.205 g6p_c + 0.2557 gln__L_c + 3.4454 glu__L_c + 59.35 h2o_c + 2.051 nad_c + 13.0279 nadph_c + 1.7867 oaa_c + 0.5191 pep_c + 2.8328 pyr_c + 0.8977 r5p_c -> 59.81 adp_c + 2.6222 akg_c + 3.7478 coa_c + 58.27 h_c + 2.051 nadh_c + 13.0279 nadp_c + 58.31 pi_c + 1.033 thf_c")
print(model_C1.reactions.BIOMASS_Ecoli_core_w_GAM)
# assuming glc uptake rate is -10.5 as WT
MeDH = []
with model_C1 as m:
    m.reactions.EX_glc__D_e.lower_bound = -10.5
    m.reactions.EX_meoh_e.lower_bound = -1000
    m.objective = 'MeDH'
    m.objective_direction = 'min'
    for i in range(len(growth_MDH)):
        gr = growth_MDH['growth_rate'].iloc[i]
        m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
        cobra.flux_analysis.pfba(m)
        MeDH.append(m.reactions.MeDH.flux)

growth_MDH['MeDH_flux'] = PA_(MeDH,dtype='mmol/gCDW/hr')

growth_MDH.pint.dequantify()
# While the maximum growth the strain can reach is as with foramte
MeDH = []
with model_C1 as m:
    m.reactions.EX_glc__D_e.lower_bound = -10.5
    m.reactions.EX_meoh_e.lower_bound = -1000
    m.objective = 'MeDH'
    m.objective_direction = 'min'
    for i in range(len(growth_MDH_for)):
        gr = growth_MDH_for['growth_rate'].iloc[i]
        m.reactions.BIOMASS_Ecoli_core_w_GAM.bounds = (gr.magnitude, gr.magnitude)
        cobra.flux_analysis.pfba(m)
        MeDH.append(m.reactions.MeDH.flux)
growth_MDH_for['MeDH_flux'] = PA_(MeDH,dtype='mmol/gCDW/hr')
growth_MDH_for.pint.dequantify()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14362325094002343, 'fraction_numerical': 0.045429328730814275, 'mean_word_length': 3.7203956939191154, 'pattern_counts': {'":': 11, '<': 45, '<?xml version=': 0, '>': 47, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '427788', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 7201, 'n_tokens_neox': 6438, 'n_words': 1355, 'path': '2020_Energy_Auxotroph/LPD_FBA.ipynb', 'repo_name': 'he-hai/PubSuppl', 'size': '16232'}