{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Network Analysis\n",
"\n",
"StrainDesign provides canonical functions for maximizing and minimizing metabolic fluxes in network. The output format is identical to the ones of COBRApy's functions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Flux optimization (FBA/pFBA)\n",
"\n",
"Flux Balance Analysis (FBA) is a linear program that optimizes a flux rate under the given steady-state network constraints. For the case of a growth-rate maximization, the problem is written as:\n",
"$$\\text{maximize}\\:\\:\\: v_{growth} \\\\ \\text{ subject to } \\\\ \\mathbf{S\\, v = 0 } \\\\ \\mathbf{ lb \\le v \\le ub}$$\n",
"\n",
"Where $\\mathbf{S}$ is the stoichiometric matrix of the metabolic model, $\\mathbf{v}$ is the vector of metabolic flux rates and $\\mathbf{lb}$ and $\\mathbf{ub}$ are the physiological lower and upper bounds of each flux rate which also define whether a reaction can run in the reverse direction ($\\mathbf{lb} < 0$) or not ($\\mathbf{lb} \\ge 0$). $\\mathbf{S\\, v = 0}$ represents all steady-state constraints, and $\\mathbf{lb \\le v \\le ub}$ the allowed flux ranges.\n",
"\n",
"\n",
"\n",
"**Note:**\n",
"\n",
"All of the following computation examples will require the COBRApy and the StrainDesign package.\n",
"\n",
"
\n",
"\n",
"Here, we load both packages and the e_coli_core model from BiGG:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Set parameter Username\n",
"Academic license - for non-commercial use only - expires 2023-07-20\n"
]
}
],
"source": [
"import cobra\n",
"import straindesign as sd\n",
"model = cobra.io.load_model('e_coli_core')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An FBA is launched by a single function call. By default the model's objective function is optimized. The function returns a solution object, in which the objective value and the fluxes are stored in `solution.objective_value` and `solution.fluxes`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Maximum growth: 0.8739215069684305.\n"
]
}
],
"source": [
"solution = sd.fba(model)\n",
"\n",
"print(f\"Maximum growth: {solution.objective_value}.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may also use a custom objective (in form of a linear expression) and change the optimization sense. Here we minimize the Glucose uptake rate through the PTS:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Minimum flux through GLCpts: 0.4794285714285715.\n"
]
}
],
"source": [
"solution = sd.fba(model,obj='GLCpts',obj_sense='minimize')\n",
"\n",
"print(f\"Minimum flux through GLCpts: {solution.objective_value}.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also consider custom constraints (in this case limited oxygen uptake and an increased fixed ATP maintenance demand):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Maximum growth at limited oxygen uptake and high ATP maintenance: 0.26305573292588313.\n"
]
}
],
"source": [
"solution = sd.fba(model,constraints=['-EX_o2_e <= 5', 'ATPM = 20'])\n",
"\n",
"print(f\"Maximum growth at limited oxygen uptake and high ATP maintenance: {solution.objective_value}.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Parsimonious FBA (pFBA)\n",
"\n",
"Parsimonious flux balance analysis optimizes a flux rate under the given steady-state network constraints, but also minimizes the sum of absolute fluxes to achieve this optimum. One can write:\n",
"$$\\text{minimize}\\:\\:\\:\\Sigma | v_i | \\\\ \\text{subject to} \\\\$$ \n",
"$$\\text{maximize}\\:\\:\\: v_{growth} \\\\ \\text{subject to} \\\\ \\mathbf{S\\, v = 0} \\\\ \\mathbf{lb \\le v \\le ub}$$\n",
"\n",
"pFBA is the simpleset (although very rough) approach of emulating a cell's enzyme cost minimization (after an assumed growth maximization).\n",
"\n",
"\n",
"\n",
"**Note:**\n",
"\n",
"pFBA solutions are more often \"unique\" than pure FBA solutions, since the outer minimization leaves fewer degrees of freedom in the solution space.\n",
"\n",
"
\n",
"\n",
"StrainDesign computes pFBA solutions when you pass the 'pFBA'-argument with a value of 1. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The sum of fluxes of the regular FBA: 2508.293334194643 is usually higher than of the parsimoniuos FBA: 518.4220855176071\n"
]
}
],
"source": [
"fba_sol = sd.fba(model)\n",
"pfba1_sol = sd.fba(model,pfba=1)\n",
"print(f\"The sum of fluxes of the regular FBA: {sum([abs(v) for v in fba_sol.fluxes.values()])} \"+\\\n",
" f\"is usually higher than of the parsimoniuos FBA: {sum([abs(v) for v in pfba1_sol.fluxes.values()])}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Likewise, it is possible to minimize the number of active reactions to attain an optimal flux distribution. We therefore use pFBA mode 2. Most of the times modes 1 and 2 yield the same results."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The number of active reactions in pFBA1: 48, and pFBA2: 48, is often identical.\n"
]
}
],
"source": [
"pfba1_sol = sd.fba(model,pfba=1)\n",
"pfba2_sol = sd.fba(model,pfba=2)\n",
"print(f\"The number of active reactions in pFBA1: {sum([v!=0 for v in pfba1_sol.fluxes.values()])}, \"+\\\n",
" f\"and pFBA2: {sum([v!=0 for v in pfba2_sol.fluxes.values()])}, is often identical.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Flux variability analysis (FVA)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The fva function determines the possible maximal and minimal flux ranges under the given model constraints:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" minimum maximum\n",
"PFK 0.0 176.61\n",
"PFL 0.0 40.00\n",
"PGI -50.0 10.00\n",
"PGK -20.0 -0.00\n",
"PGL 0.0 60.00\n",
"... ... ...\n",
"NADH16 0.0 120.00\n",
"NADTRHD 0.0 378.22\n",
"NH4t 0.0 10.00\n",
"O2t 0.0 60.00\n",
"PDH 0.0 40.00\n",
"\n",
"[95 rows x 2 columns]\n"
]
}
],
"source": [
"solution = sd.fva(model)\n",
"print(solution)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parameters 'solver' and 'constraints' can also be used in the FVA function call. As an example, we determine the flux ranges for the case that the flux sum of PDH and PFL is smaller than 8."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" minimum maximum\n",
"PFK 0.0 147.610000\n",
"PFL 0.0 8.000000\n",
"PGI -50.0 10.000000\n",
"PGK -20.0 -0.000000\n",
"PGL 0.0 60.000000\n",
"... ... ...\n",
"NADH16 0.0 120.000000\n",
"NADTRHD 0.0 375.220000\n",
"NH4t 0.0 8.300263\n",
"O2t 0.0 60.000000\n",
"PDH 0.0 8.000000\n",
"\n",
"[95 rows x 2 columns]\n"
]
}
],
"source": [
"solution = sd.fva(model, constraints='PDH + PFL <= 8')\n",
"print(solution)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Yield optimization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yield optmization aims to maximize a given flux expression (e.g., a product's synthesis rate) divided by another such expression (e.g., a substrate's uptake rate).\n",
"\n",
"\n",
"\n",
"**Warning:**\n",
"\n",
"Yield optimization and Flux Balance Analysis are two different methods that produce distinct optimal values and also distinct optimal flux vectors.\n",
"\n",
"
\n",
"\n",
"Consider the following example of maximizing growth and biomass yield under limited oxygen uptake."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Maximum yield (FBA): 0.031965425062067315.\n",
"Maximum yield (yOpt): 0.03629426243040193.\n"
]
}
],
"source": [
"numerator = 'BIOMASS_Ecoli_core_w_GAM'\n",
"denominator = '-EX_glc__D_e'\n",
"constraint = '-EX_o2_e <= 3'\n",
"\n",
"sol_fba = sd.fba(model,obj='BIOMASS_Ecoli_core_w_GAM',obj_sense='maximize',constraints=constraint)\n",
"fba_yield = sol_fba.fluxes[numerator] / -sol_fba.fluxes['EX_glc__D_e'] \n",
"\n",
"sol_yopt = sd.yopt(model,obj_num=numerator,obj_den=denominator,obj_sense='maximize',constraints=constraint)\n",
"yopt_yield = sol_yopt.objective_value\n",
"\n",
"print(f\"Maximum yield (FBA): {fba_yield}.\")\n",
"print(f\"Maximum yield (yOpt): {yopt_yield}.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The best biomass yield is achieved when only respiration is used. The best growth rates use respiration and *additionally* overflow metabolism with poorer biomass yield."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plotting chapter shows how relationships between yields and rates can be visualized."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Mathematical background\n",
"\n",
"Constraint-based models ($\\mathbf{S\\, v = 0}$, $\\mathbf{lb \\le v \\le ub}$) can be rewritten in a single matrix-inequality term ($\\mathbf{A \\, x \\le b}$).\n",
"\n",
"With this notation, the yield optimization is a linear fractional program (LFP):\n",
"$$\n",
"\\begin{array}{l c}\n",
"\t\\text{maximize} \\; &\\mathbf{\\tfrac{c^\\intercal x}{d^\\intercal x}} \\\\\n",
"\t\\text{subject to} &\\mathbf{A \\, x \\le b}.\n",
"\\end{array}\n",
"$$\n",
"Under the condition that the denominator term is strictly positive ($\\mathbf{A\\,x\\le b}\\, \\Rightarrow\\, \\mathbf{d}^\\intercal \\mathbf{x}>\\,0$), the LFP may be rewritten as an LP problem, using the Charnes-Cooper transformation. The formerly fixed boundaries $\\mathbf{b}$ of the problem are then scaled by the auxiliary variable $e$ while the variable $y = \\mathbf{\\tfrac{c^\\intercal x}{d^\\intercal x}}$ expresses the original objective function:\n",
"$$\n",
"\\begin{array}{l c}\n",
"\t\\text{maximize} \\; & y \\\\\\\\\n",
"\t\\text{subject to} &\n",
"\t\\begin{bmatrix}\n",
"\t\t\\mathbf{A~} & \\mathbf{-b} & \\mathbf{0} \\\\\n",
"\t\t\\mathbf{d}^\\intercal & 0 \t\t & 0 \\\\\n",
"\t\t\\mathbf{c}^\\intercal & 0 \t\t & -1 \\\\\n",
"\t\\end{bmatrix} \n",
"\t\\begin{bmatrix}\\mathbf{\\tilde{x}} \\\\ e \\\\ y \\end{bmatrix} \n",
"\t\\begin{matrix} \\le \\\\ = \\\\ = \\end{matrix} \n",
"\t\\begin{bmatrix}\\mathbf{0} \\\\ 1 \\\\ 0 \\end{bmatrix}\\\\\\\\\n",
"\t&e \\ge 0.\n",
"\\end{array}\n",
"$$\n",
"\n",
"Solutions of $\\mathbf{x}$ to the LFP (first problem) can be retrieved from a solution of the LP through $\\mathbf{x}=\\frac{\\mathbf{\\tilde{x}}}{e}$.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.7 ('cnapy-dev')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "50a36e699ca3834a05ca1bd86fc5b7db4829f1cbbb3e27f51544575d88525899"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}