5  使用FBA进行模拟

使用磁通平衡分析的仿真可以使用 Model.optimize() 进行求解。这将最大化或最小化(最大化是默认设置)通过客观反应的通量。

from cobra.io import load_model
model = load_model("textbook")

5.1 运行 FBA

solution = model.optimize()
print(solution)
<Solution 0.874 at 0x195d0067050>

Model.optimize()方法会给出一个解决方案(Solution)对象。这个对象包含了几个关键属性:

  • objective_value:表示目标函数的优化值,也就是模型追求最大化或最小化的那个数值结果。
  • status:显示线性规划求解器给出的解决状态,告诉你问题是否得到解决、有没有特殊情况需要留意。
  • fluxes:这是一个按照反应标识符来组织的pandas序列,里面记录了每个反应的通量值。通量是指正向反应速率与逆向反应速率之差,直观反映了物质在该反应中的净流动情况。
  • shadow_prices:同样是一个pandas序列,不过这次是依据代谢物标识符来排列的,展示了每个代谢物的影子价格。影子价格揭示了在保持目标函数最优的前提下,如果系统可利用的某种代谢物增多或减少一个单位,目标函数值的变化量,有助于理解代谢物的重要性及限制条件。

例如,在最后一次调用“model.optimize()”之后,如果优化成功,则其状态将为最佳。如果模型不可行,则会引发错误。

solution.objective_value
0.8739215069684279

可与 cobrapy 一起使用的求解器速度非常快,以至于对于许多中小型模型来说,计算求解的速度甚至比从求解器收集值并将 python 对象转换为它们所需的速度还要快。使用“model.optimize”,我们收集所有反应和代谢物的值,如果重复进行,这可能需要大量时间。如果我们只对单个反应或目标的通量值感兴趣,那么使用“model.slim_optimize”会更快,它只执行优化并返回目标值,由你来获取你可能需要的其他值。

%%time
model.optimize().objective_value
CPU times: total: 0 ns
Wall time: 1e+03 µs
0.8739215069684305
%%time
model.slim_optimize()
CPU times: total: 0 ns
Wall time: 1e+03 µs
0.8739215069684305

5.1.1 分析 FBA solutions

使用FBA求解的模型可以通过使用摘要方法进一步分析,这些方法输出打印文本以快速表示模型行为。对整个模型调用 summary 方法将显示有关模型的输入和输出行为的信息,以及优化的目标。

model.summary()

Objective

1.0 Biomass_Ecoli_core = 0.8739215069684297

Uptake

Metabolite Reaction Flux C-Number C-Flux
glc__D_e EX_glc__D_e 10 6 100.00%
nh4_e EX_nh4_e 4.765 0 0.00%
o2_e EX_o2_e 21.8 0 0.00%
pi_e EX_pi_e 3.215 0 0.00%

5.1.1.1 Secretion

Metabolite Reaction Flux C-Number C-Flux
co2_e EX_co2_e -22.81 1 100.00%
h2o_e EX_h2o_e -29.18 0 0.00%
h_e EX_h_e -17.53 0 0.00%

此外,还可以使用汇总方法检查单个代谢物的输入输出行为。例如,以下命令可用于检查模型的整体氧化还原平衡

model.metabolites.nadh_c.summary()

nadh_c

C21H27N7O14P2

Producing Reactions

Percent Flux Reaction Definition
13.14% 5.064 AKGDH akg_c + coa_c + nad_c --> co2_c + nadh_c + succoa_c
8.04% 3.1 Biomass_Ecoli_core 1.496 3pg_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 + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 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 + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c
41.58% 16.02 GAPD g3p_c + nad_c + pi_c <=> 13dpg_c + h_c + nadh_c
13.14% 5.064 MDH mal__L_c + nad_c <=> h_c + nadh_c + oaa_c
24.09% 9.283 PDH coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c

5.1.1.2 Consuming Reactions

Percent Flux Reaction Definition
100.00% -38.53 NADH16 4.0 h_c + nadh_c + q8_c --> 3.0 h_e + nad_c + q8h2_c

或者了解主要的能源生产和消费反应

model.metabolites.atp_c.summary()

atp_c

C10H12N5O13P3

Producing Reactions

Percent Flux Reaction Definition
66.58% 45.51 ATPS4r adp_c + 4.0 h_e + pi_c <=> atp_c + h2o_c + 3.0 h_c
23.44% 16.02 PGK 3pg_c + atp_c <=> 13dpg_c + adp_c
2.57% 1.758 PYK adp_c + h_c + pep_c --> atp_c + pyr_c
7.41% 5.064 SUCOAS atp_c + coa_c + succ_c <=> adp_c + pi_c + succoa_c

5.1.1.3 Consuming Reactions

Percent Flux Reaction Definition
12.27% -8.39 ATPM atp_c + h2o_c --> adp_c + h_c + pi_c
76.46% -52.27 Biomass_Ecoli_core 1.496 3pg_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 + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 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 + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c
0.33% -0.2235 GLNS atp_c + glu__L_c + nh4_c --> adp_c + gln__L_c + h_c + pi_c
10.94% -7.477 PFK atp_c + f6p_c --> adp_c + fdp_c + h_c

5.2 更改目标

目标函数由客观反应的objective_coefficient属性确定。通常,使用描述构成细胞的代谢物组成的“生物量”函数。

biomass_rxn = model.reactions.get_by_id("Biomass_Ecoli_core")

目前在模型中,目标中只有一个反应(生物质反应),线性系数为1。

from cobra.util.solver import linear_reaction_coefficients
linear_reaction_coefficients(model)
{<Reaction Biomass_Ecoli_core at 0x195cff33490>: 1.0}

可以通过分配 Model.objective 来更改目标函数,Model.objective 可以是反应对象(或只是它的名称),也可以是 {Reaction: objective_coefficient} 的字典。

# 将目标改为 ATPM
model.objective = "ATPM"

# 上界应该是1000,这样我们才能得到实际的最佳值
model.reactions.get_by_id("ATPM").upper_bound = 1000.
linear_reaction_coefficients(model)
{<Reaction ATPM at 0x195cff09450>: 1.0}
model.optimize().objective_value
175.00000000000006

我们还可以有更复杂的目标,比如二次项。

5.3 运行 FVA

FBA不会总是给出独特的解决方案,因为多个通量状态可以达到相同的最佳状态。FVA(或通量变异性分析)找到每个代谢通量处于最佳状态的范围。

import cobra
from cobra.flux_analysis import flux_variability_analysis
flux_variability_analysis(model, model.reactions[:10])
minimum maximum
ACALD -9.375513e-15 0.000000e+00
ACALDt -1.198454e-14 0.000000e+00
ACKr -1.832929e-14 1.039260e-14
ACONTa 2.000000e+01 2.000000e+01
ACONTb 2.000000e+01 2.000000e+01
ACt2r -3.615752e-15 0.000000e+00
ADK1 0.000000e+00 1.924093e-13
AKGDH 2.000000e+01 2.000000e+01
AKGt2r -9.821844e-15 0.000000e+00
ALCD2x 9.994214e-15 0.000000e+00

设置参数“fraction_of_optimium=0.90”将给出 90% 最佳反应的通量范围。

cobra.flux_analysis.flux_variability_analysis(
    model, model.reactions[:10], fraction_of_optimum=0.9)
minimum maximum
ACALD -2.692308 0.000000e+00
ACALDt -2.692308 0.000000e+00
ACKr -4.117647 1.039260e-14
ACONTa 8.461538 2.000000e+01
ACONTb 8.461538 2.000000e+01
ACt2r -4.117647 0.000000e+00
ADK1 0.000000 1.750000e+01
AKGDH 2.500000 2.000000e+01
AKGt2r -1.489362 0.000000e+00
ALCD2x -2.333333 0.000000e+00

标准FVA可能包含环路,即高绝对通量值,只有在允许它们参与环路时才能很高(一种在体内不会发生的数学伪影)。使用“loopless”参数来避免此类循环。接下来,我们可以在环中看到 FRD7 和 SUCDi 反应,与此同时可使用无环 FVA 来避免这种情况。

loop_reactions = [model.reactions.FRD7, model.reactions.SUCDi]
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=False)
minimum maximum
FRD7 0.0 980.0
SUCDi 20.0 1000.0
flux_variability_analysis(model, reaction_list=loop_reactions, loopless=True)
minimum maximum
FRD7 0.0 0.0
SUCDi 20.0 20.0

5.3.1 在摘要方法中运行 FVA

通量变异性分析也可以嵌入到对摘要方法的调用中。例如,可以通过以下方式快速找到基材消耗和产物形成的预期变化

model.optimize()
model.summary(fva=0.95)

Objective

1.0 ATPM = 175.0000000000001

Uptake

Metabolite Reaction Flux Range C-Number C-Flux
glc__D_e EX_glc__D_e 10 [9.5; 10] 6 100.00%
o2_e EX_o2_e 60 [55.88; 60] 0 0.00%

5.3.1.1 Secretion

Metabolite Reaction Flux Range C-Number C-Flux
ac_e EX_ac_e 0 [-2.059; 0] 2 0.00%
acald_e EX_acald_e 0 [-1.346; 0] 2 0.00%
akg_e EX_akg_e 0 [-0.7447; 0] 5 0.00%
co2_e EX_co2_e -60 [-60; -54.17] 1 100.00%
etoh_e EX_etoh_e 0 [-1.167; 0] 2 0.00%
for_e EX_for_e 0 [-5.833; 0] 1 0.00%
glu__L_e EX_glu__L_e 0 [-0.6731; 0] 5 0.00%
h2o_e EX_h2o_e -60 [-60; -54.17] 0 0.00%
h_e EX_h_e 0 [-5.833; 0] 0 0.00%
lac__D_e EX_lac__D_e 0 [-1.129; 0] 3 0.00%
nh4_e EX_nh4_e 0 [0; 0.6731] 0 0.00%
pi_e EX_pi_e 0 [0; 0.171] 0 0.00%
pyr_e EX_pyr_e 0 [-1.346; 0] 3 0.00%
succ_e EX_succ_e 0 [-0.875; 0] 4 0.00%

同样,代谢物质量平衡的变异性也可以通过通量变异性分析来检查。

model.metabolites.pyr_c.summary(fva=0.95)

pyr_c

C3H3O3

Producing Reactions

Percent Flux Range Reaction Definition
50.00% 10 [9.5; 10] GLCpts glc__D_e + pep_c --> g6p_c + pyr_c
0.00% 0 [-1.129; 0] LDH_D lac__D_c + nad_c <=> h_c + nadh_c + pyr_c
0.00% 0 [0; 8.75] ME1 mal__L_c + nad_c --> co2_c + nadh_c + pyr_c
0.00% 0 [0; 8.75] ME2 mal__L_c + nadp_c --> co2_c + nadph_c + pyr_c
50.00% 10 [1.25; 18.75] PYK adp_c + h_c + pep_c --> atp_c + pyr_c
0.00% 0 [-1.346; 0] PYRt2 h_e + pyr_e <=> h_c + pyr_c

5.3.1.2 Consuming Reactions

Percent Flux Range Reaction Definition
0.00% 0 [-0.1316; 0] Biomass_Ecoli_core 1.496 3pg_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 + 4.9414 glu__L_c + 59.81 h2o_c + 3.547 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 + 4.1182 akg_c + 3.7478 coa_c + 59.81 h_c + 3.547 nadh_c + 13.0279 nadp_c + 59.81 pi_c
100.00% -20 [-28.75; -13] PDH coa_c + nad_c + pyr_c --> accoa_c + co2_c + nadh_c
0.00% 0 [-5.833; 0] PFL coa_c + pyr_c --> accoa_c + for_c
0.00% 0 [-8.75; 0] PPS atp_c + h2o_c + pyr_c --> amp_c + 2.0 h_c + pep_c + pi_c

在这些汇总方法中,这些值报告为中心点 +/- FVA 解的范围,由最大值和最小值计算得出。

5.4 运行 pFBA

简约的 FBA(通常写成 pFBA)找到一个通量分布,该分布可提供最佳增长率,但使通量总和最小化。这涉及求解两个顺序线性程序,但由 cobrapy 透明地处理。有关pFBA的更多详细信息,请参阅 Lewis et al. (2010).

model.objective = 'Biomass_Ecoli_core'
fba_solution = model.optimize()
pfba_solution = cobra.flux_analysis.pfba(model)

这些函数可能会给出完全不同的目标值,因为 pFBA 显示的目标值定义为 ‘sum(abs(pfba_solution.fluxes.values))’,而标准 FBA 的目标值定义为通过优化反应的加权通量(例如 ‘fba_solution.fluxes[“Biomass_Ecoli_core”]’)。

pFBA 和 FBA 都应在求解器容差范围内为要优化的目标返回相同的结果。例如,对于使反应“Biomass_Ecoli_core”最大化的 FBA 问题:

abs(fba_solution.fluxes["Biomass_Ecoli_core"] - pfba_solution.fluxes[
    "Biomass_Ecoli_core"])
5.551115123125783e-16
import numpy as np
np.isclose(
    fba_solution.fluxes["Biomass_Ecoli_core"], 
    pfba_solution.fluxes["Biomass_Ecoli_core"]
)
True

5.5 运行几何FBA

几何 FBA 找到一个独特的最佳通量分布,这是可能通量范围的核心。有关几何FBA的更多详细信息,请参阅 K Smallbone, E Simeonidis (2009).

geometric_fba_sol = cobra.flux_analysis.geometric_fba(model)
geometric_fba_sol
Optimal solution with objective value 0.000
fluxes reduced_costs
ACALD 0.000000e+00 0.0
ACALDt 0.000000e+00 0.0
ACKr 1.260427e-13 0.0
ACONTa 6.007250e+00 0.0
ACONTb 6.007250e+00 0.0
... ... ...
TALA 1.496984e+00 0.0
THD2 0.000000e+00 0.0
TKT1 1.496984e+00 0.0
TKT2 1.181498e+00 0.0
TPI 7.477382e+00 0.0

95 rows × 2 columns