{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3ef581f4-2a5b-4070-8a47-cf798b876941",
   "metadata": {},
   "source": [
    "# PyCoMo discretized growth #\n",
    "This tutorial show-cases the use of discretized growth for optimization across the full solution space. \n",
    "With this, linearizing the model with fixed abundance or growth rate is no longer necessary. \n",
    "However, this approach ustilizes binary variable, using a MILP to solve. \n",
    "This is faster than non-linear solvers, but slower than LPs.\n",
    "\n",
    "The implementation is based on an idea by Salvy and Hatzimanikatis:\n",
    "\n",
    "Salvy, P., Hatzimanikatis, V. \n",
    "The ETFL formulation allows multi-omics integration in thermodynamics-compliant metabolism and expression models. \n",
    "Nat Commun 11, 30 (2020). \n",
    "https://doi.org/10.1038/s41467-019-13818-7\n",
    "\n",
    "The expected runtime for this notebook is less than 10 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bd639378-75b0-4d05-8781-69e82451bc78",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2026-03-15 15:39:37,627 - PyCoMo - INFO - Logger initialized.\n"
     ]
    }
   ],
   "source": [
    "import cobra\n",
    "import pycomo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e835ead6-7fc3-4272-934e-125d687de369",
   "metadata": {},
   "outputs": [],
   "source": [
    "cobra.Configuration().solver=\"glpk\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2641bff9-26df-47e5-994c-79a26e1296d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pycomo.configure_logger(level=\"info\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d308cb07-ef8d-44bd-9896-bf2cfa9bc860",
   "metadata": {},
   "source": [
    "## Base Model ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "65127ba9-3ce4-4207-a761-e9e8a266142e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_single_org_choice(is_one=True):\n",
    "    model_id = 1 if is_one else 2\n",
    "    other_id = 2 if is_one else 1\n",
    "    model = cobra.Model()\n",
    "\n",
    "    model.name = f\"Org{model_id}\"\n",
    "    \n",
    "    model.add_metabolites([cobra.Metabolite(i) for i in [\n",
    "        \"S\", \n",
    "        f\"X{model_id}\", \n",
    "        \"Y1\",\n",
    "        \"Y2\",\n",
    "        \"A\", \n",
    "        \"P\", \n",
    "        \"BM\"]])\n",
    "    model.add_reactions([cobra.Reaction(i) for i in [\n",
    "        f\"EX_X{model_id}\", \n",
    "        \"EX_Y1\", \n",
    "        \"EX_Y2\", \n",
    "        f\"TP_X{model_id}\", \n",
    "        \"TP_Y1\", \n",
    "        \"TP_Y2\", \n",
    "        \"v4\", \n",
    "        \"v5\", \n",
    "        \"bio\"]])\n",
    "\n",
    "    for internal_met in [\n",
    "        \"S\", \n",
    "        \"A\", \n",
    "        \"P\", \n",
    "        \"BM\"\n",
    "    ]:\n",
    "        model.metabolites.get_by_id(internal_met).compartment = \"i\"\n",
    "\n",
    "    for external_met in [\n",
    "        f\"X{model_id}\", \n",
    "        \"Y1\",\n",
    "        \"Y2\",\n",
    "    ]:\n",
    "        model.metabolites.get_by_id(external_met).compartment = \"e\"\n",
    "\n",
    "    for met in model.metabolites:\n",
    "        met.name = met.id\n",
    "\n",
    "    for rxn in model.reactions:\n",
    "        rxn.name = rxn.id\n",
    "    \n",
    "    model.reactions.get_by_id(f\"EX_X{model_id}\").add_metabolites({f\"X{model_id}\": -1})\n",
    "    model.reactions.get_by_id(f\"TP_X{model_id}\").add_metabolites({f\"X{model_id}\": -1, \"S\": 1})\n",
    "\n",
    "    model.reactions.get_by_id(f\"EX_Y{model_id}\").add_metabolites({f\"Y{model_id}\": -1})\n",
    "    model.reactions.get_by_id(f\"EX_Y{other_id}\").add_metabolites({f\"Y{other_id}\": -1})\n",
    "    model.reactions.get_by_id(f\"TP_Y{model_id}\").add_metabolites({f\"Y{model_id}\": -1, \"A\": 1})\n",
    "    model.reactions.get_by_id(f\"TP_Y{other_id}\").add_metabolites({f\"Y{other_id}\": -1, \"P\": 1})\n",
    "    \n",
    "    model.reactions.get_by_id(\"v4\").add_metabolites({\"S\": -1, \"BM\": 1, \"P\": 1})\n",
    "    model.reactions.get_by_id(\"v5\").add_metabolites({\"A\": -1, \"BM\": 1})\n",
    "    \n",
    "    model.reactions.get_by_id(\"bio\").add_metabolites({\"BM\": -1})\n",
    "\n",
    "    model.reactions.get_by_id(f\"EX_X{model_id}\").bounds = (-1000.,0.)\n",
    "    model.reactions.get_by_id(f\"TP_X{model_id}\").bounds = (0.,1.)\n",
    "\n",
    "    model.reactions.get_by_id(f\"EX_Y{model_id}\").bounds = (0.,1000.)\n",
    "    model.reactions.get_by_id(f\"EX_Y{other_id}\").bounds = (0.,1000.)\n",
    "    model.reactions.get_by_id(f\"TP_Y{model_id}\").bounds = (0.,1.)\n",
    "    model.reactions.get_by_id(f\"TP_Y{other_id}\").bounds = (-1.,0.)\n",
    "    \n",
    "    model.reactions.get_by_id(\"v4\").bounds = (0.,1000.)\n",
    "    model.reactions.get_by_id(\"v5\").bounds = (0.,1000.)\n",
    "    \n",
    "    model.reactions.get_by_id(\"bio\").bounds = (0.,1000.)\n",
    "\n",
    "    model.objective = \"bio\"\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bb759890-99e6-42a6-80f1-96c09ab4ebe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_toy_pycomo_model():\n",
    "    model1 = construct_single_org_choice(is_one=True)\n",
    "    model2 = construct_single_org_choice(is_one=False)\n",
    "    \n",
    "    single_org1 = pycomo.SingleOrganismModel(model1, name=\"org1\")\n",
    "    single_org2 = pycomo.SingleOrganismModel(model2, name=\"org2\")\n",
    "    \n",
    "    com_model = pycomo.CommunityModel([single_org1, single_org2], name=\"SymmetricToy\")\n",
    "    \n",
    "    com_model.medium = {'EX_X1_medium': 1000.0,\n",
    "     'EX_Y1_medium': 0.,\n",
    "     'EX_Y2_medium': 0.,\n",
    "     'EX_X2_medium': 1000.0}\n",
    "    com_model.apply_medium()\n",
    "    return com_model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db87718-44e4-4135-a993-ba581d15c6ca",
   "metadata": {},
   "source": [
    "## Find maximum growth rate ##\n",
    "The maximum growth rate is 2.0 at 50%/50% abundance of the two organisms. It can be found quickly by the discretized growth rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8588935b-9cd1-40ef-b201-8cc0760ef9b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2026-03-15 15:39:37,924 - PyCoMo - INFO - No community model generated yet. Generating now:\n",
      "2026-03-15 15:39:37,930 - PyCoMo - INFO - Identified biomass reaction from objective: bio\n",
      "2026-03-15 15:39:37,930 - PyCoMo - INFO - Note: no products in the objective function, adding biomass to it.\n",
      "2026-03-15 15:39:37,931 - PyCoMo - INFO - Biomass reaction bio is a boundary reaction.\n",
      "2026-03-15 15:39:37,931 - PyCoMo - INFO - Moving reactant BM out of the boundary compartment into compartment 'bio', to avoid further boundary reactions being added to it.\n",
      "2026-03-15 15:39:37,956 - PyCoMo - INFO - Identified biomass reaction from objective: bio\n",
      "2026-03-15 15:39:37,956 - PyCoMo - INFO - Note: no products in the objective function, adding biomass to it.\n",
      "2026-03-15 15:39:37,957 - PyCoMo - INFO - Biomass reaction bio is a boundary reaction.\n",
      "2026-03-15 15:39:37,958 - PyCoMo - INFO - Moving reactant BM out of the boundary compartment into compartment 'bio', to avoid further boundary reactions being added to it.\n",
      "2026-03-15 15:39:37,994 - PyCoMo - WARNING - No annotation overlap found for matching several metabolites (2). Please make sure that the matched metabolites are indeed representing the same substance in all models! The list of metaboliteswithout annotation overlap can be accessed via 'model.no_annotation_overlap'\n",
      "2026-03-15 15:39:38,007 - PyCoMo - WARNING - There are 18 metabolites without elements in the model. Mass balance checks may be unreliable.\n",
      "2026-03-15 15:39:38,008 - PyCoMo - INFO - Generated community model.\n"
     ]
    }
   ],
   "source": [
    "com_model = create_toy_pycomo_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ee50d5b6-c051-499d-b2de-18757c959607",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 <= 1 <= 1, 0 <= 2 <= 1, 0 <= 4 <= 1, 0 <= 8 <= 1, 0 <= 16 <= 1, 0 <= 32 <= 1, 0 <= 64 <= 1, 0 <= 128 <= 1, 0 <= 256 <= 1, 0 <= 512 <= 1]\n"
     ]
    }
   ],
   "source": [
    "pycomo.helper.discretized_growth.relax_model_linearisation(com_model)\n",
    "f_vars = pycomo.helper.discretized_growth.init_f_variables(com_model)\n",
    "pycomo.helper.discretized_growth.discretize_growth_of_model(com_model, f_vars, mu_lb=0., mu_ub=100., n_mu_bins=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d0afb67d-e50f-477f-9eb7-56286297bf22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong><em>Optimal</em> solution with objective value 2.000</strong><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fluxes</th>\n",
       "      <th>reduced_costs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>org1_TF_X1_org1_e</th>\n",
       "      <td>-0.50</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TF_Y1_org1_e</th>\n",
       "      <td>-0.50</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TF_Y2_org1_e</th>\n",
       "      <td>0.50</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TP_X1</th>\n",
       "      <td>0.50</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TP_Y1</th>\n",
       "      <td>0.50</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SK_org2_bio_ub</th>\n",
       "      <td>4.99</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SK_org2_to_community_biomass_ub</th>\n",
       "      <td>4.99</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>f_final</th>\n",
       "      <td>1.00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>abundance_reaction</th>\n",
       "      <td>0.00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>community_biomass</th>\n",
       "      <td>2.00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>55 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "<Solution 2.000 at 0x7f0352a94070>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "com_model.model.objective = \"community_biomass\"\n",
    "sol = com_model.model.optimize(\"maximize\")\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44610f49-2b3d-42e6-9810-c1168796f8f0",
   "metadata": {},
   "source": [
    "## Find maximum of flux ##\n",
    "flux v4 of organism 1 can reach its maximum flux of 1.0 at a growth rate 1.0 - not at maximum growth. See if the optimum is found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "241dbf93-922e-4bc4-8a77-9745cbeb70c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong><em>Optimal</em> solution with objective value 1.000</strong><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fluxes</th>\n",
       "      <th>reduced_costs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>org1_TF_X1_org1_e</th>\n",
       "      <td>-1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TF_Y1_org1_e</th>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TF_Y2_org1_e</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TP_X1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>org1_TP_Y1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SK_org2_bio_ub</th>\n",
       "      <td>-0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SK_org2_to_community_biomass_ub</th>\n",
       "      <td>-0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>f_final</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>abundance_reaction</th>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>community_biomass</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>55 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "<Solution 1.000 at 0x7f0352a94ac0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "com_model.model.objective = \"org1_v4\"\n",
    "sol = com_model.model.optimize(\"maximize\")\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c94bfba-cc74-41d8-9037-40fada1b2526",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pycomo_efmtool_dev_cplex_3.10",
   "language": "python",
   "name": "pycomo_efmtool_dev_cplex_3.10"
  },
  "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.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}