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