{
"cells": [
{
"cell_type": "markdown",
"id": "f2d7b3c9",
"metadata": {},
"source": [
"# Sparse FBA"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "a6747aa4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Set parameter Username\n",
"Academic license - for non-commercial use only - expires 2024-11-12\n"
]
},
{
"data": {
"text/plain": [
"(72, 95)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from cobra.io import load_model\n",
"\n",
"model = load_model(\"textbook\")\n",
"len(model.metabolites), len(model.reactions)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "a1684e94",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"
\n",
" \n",
" \n",
" \n",
" | \n",
" \n",
" \n",
" Installed version: | v0.9.2-beta.1 (latest: v0.9.1-alpha.6) | Available backends: | CVXPY v1.4.1, PICOS v2.4.17 | Default backend (corneto.K): | CVXPY | Installed solvers: | CBC, CLARABEL, COPT, CPLEX, CVXOPT, DIFFCP, ECOS, ECOS_BB, GLPK, GLPK_MI, GUROBI, MOSEK, OSQP, PROXQP, SCIP, SCIPY, SCS | Graphviz version: | v0.20.1 | Repository: | https://github.com/saezlab/corneto | \n",
" \n",
" | \n",
"
\n",
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import corneto as cn\n",
"\n",
"cn.info()"
]
},
{
"cell_type": "markdown",
"id": "d5189703",
"metadata": {},
"source": [
"## Manual implementation of sparse FBA"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7454c898",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'_flow': Variable((95,), _flow), 'flow': Variable((95,), _flow)}"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from corneto.methods.metabolism.fba import fba_problem\n",
"\n",
"G = cn.Graph.from_cobra_model(model)\n",
"P = fba_problem(G)\n",
"P.expr"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "e0b2542a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.8739215069684304\n"
]
}
],
"source": [
"rid = list(G.get_edges_by_attr(\"id\", \"Biomass_Ecoli_core\"))[0]\n",
"biomass = P.expr.flow[rid]\n",
"P.add_objectives(biomass, weights=-1)\n",
"P.solve()\n",
"print(biomass.value)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "5aa68f5a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.5"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Add custom constraints to the problem\n",
"P += biomass <= 0.5\n",
"P.solve(solver=\"SCIPY\")\n",
"biomass.value"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "98608b05",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Warning: node 'acald_c', graph '%3' size too small for label\n",
"Warning: node 'accoa_c', graph '%3' size too small for label\n",
"Warning: node 'nadh_c', graph '%3' size too small for label\n",
"Warning: node 'acald_e', graph '%3' size too small for label\n",
"Warning: node 'acon_C_c', graph '%3' size too small for label\n",
"Warning: node 'amp_c', graph '%3' size too small for label\n",
"Warning: node 'succoa_c', graph '%3' size too small for label\n",
"Warning: node 'etoh_c', graph '%3' size too small for label\n",
"Warning: node 'nadph_c', graph '%3' size too small for label\n",
"Warning: node 'glu__L_c', graph '%3' size too small for label\n",
"Warning: node 'gln__L_c', graph '%3' size too small for label\n",
"Warning: node 'nadp_c', graph '%3' size too small for label\n",
"Warning: node 'q8h2_c', graph '%3' size too small for label\n",
"Warning: node 'lac__D_e', graph '%3' size too small for label\n",
"Warning: node 'lac__D_c', graph '%3' size too small for label\n",
"Warning: node 'etoh_e', graph '%3' size too small for label\n",
"Warning: node 'glc__D_e', graph '%3' size too small for label\n",
"Warning: node 'gln__L_e', graph '%3' size too small for label\n",
"Warning: node 'glu__L_e', graph '%3' size too small for label\n",
"Warning: node 'mal__L_e', graph '%3' size too small for label\n",
"Warning: node 'succ_e', graph '%3' size too small for label\n",
"Warning: node 'dhap_c', graph '%3' size too small for label\n",
"Warning: node 'succ_c', graph '%3' size too small for label\n",
"Warning: node 'mal__L_c', graph '%3' size too small for label\n",
"Warning: node '6pgl_c', graph '%3' size too small for label\n",
"Warning: node '13dpg_c', graph '%3' size too small for label\n",
"Warning: node '6pgc_c', graph '%3' size too small for label\n",
"Warning: node 'ru5p__D_c', graph '%3' size too small for label\n",
"Warning: node 'xu5p__D_c', graph '%3' size too small for label\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"G.plot(custom_edge_attr=cn.pl.flow_style(P), layout=\"fdp\")"
]
},
{
"cell_type": "markdown",
"id": "81c714ae",
"metadata": {},
"source": [
"Now we are going to see how to use specialized building blocks from CORNETO to automatically add indicators for the presence of reactions in a model. Indicators are binary variables associated to a continuous variable: when the indicator is 0, the continuous variable is 0; and when the indicator is 1, the continuous variable is unconstrained. \n",
"\n",
"By creating indicators for flow, we can try for example to minimize the number of reactions with flux in the model, by minimizing the total number of ones in the indicator vector."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e0c07efb",
"metadata": {},
"outputs": [],
"source": [
"from corneto.backend._base import Indicator\n",
"\n",
"# Automatically create a new indicator variable for the flow\n",
"# If the indicator is 0, the flow is blocked, if it is 1, the flow is unblocked (can take any value within bounds, including 0)\n",
"P += Indicator(\"unblocked_flow\")\n",
"P.expr"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "03c7efa4",
"metadata": {},
"outputs": [],
"source": [
"# We try to minimize the number of 1s (i.e., maximize reactions with blocked flux)\n",
"P.add_objectives(sum(P.expr.unblocked_flow))\n",
"P.solve()\n",
"print(biomass.value)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fb5e3254",
"metadata": {},
"outputs": [],
"source": [
"# How many reactions are unblocked?\n",
"sum(P.expr.unblocked_flow.value)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "294c77d5",
"metadata": {},
"outputs": [],
"source": [
"G.plot(custom_edge_attr=cn.pl.flow_style(P), layout=\"fdp\")"
]
},
{
"cell_type": "markdown",
"id": "432c7910",
"metadata": {},
"source": [
"## Sparse FBA with min. flux"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2cbe27ec",
"metadata": {},
"outputs": [],
"source": [
"P = fba_problem(G)\n",
"P += P.expr.flow[rid] >= 0.2\n",
"P += Indicator(\"unblocked_flow\")\n",
"P.add_objectives(sum(P.expr.unblocked_flow))\n",
"P.solve(solver=\"SCIPY\")"
]
},
{
"cell_type": "markdown",
"id": "2a846d74",
"metadata": {},
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "601a7432",
"metadata": {},
"outputs": [],
"source": [
"sum(P.expr.unblocked_flow.value)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "55060e52",
"metadata": {},
"outputs": [],
"source": [
"sum(abs(P.expr.flow.value) >= 1e-6)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.10"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}