4 - benchmark-simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark - Simulations -"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import biolqm\n",
    "import subprocess\n",
    "import pandas as pd\n",
    "from glob import glob\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "from colomoto.minibn import BooleanNetwork\n",
    "from merrin.bnet2flexflux import bn2flexflux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Files\n",
    "metabolic_network_file = f'data/covert/metabolic_network.xml'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Benchmark\n",
    "benchmark_dir = './benchmark/results/bn_index.csv'\n",
    "\n",
    "# Bnet\n",
    "bnet_dir = './benchmark/results/bnet/*.bnet'\n",
    "\n",
    "# Simulations\n",
    "simulations_dir = './benchmark/results/simulations/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_folder(folder):\n",
    "    if not os.path.exists(folder):\n",
    "        os.makedirs(folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulations parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regulators = {'RPO2', 'RPb', 'RPcl', 'RPh'}\n",
    "\n",
    "constraints = \"\"\"\\\n",
    "obj : MAX(Growth)\n",
    "\n",
    "\n",
    "Tc1 0 10.5\n",
    "Tc2 0 10.5\n",
    "Td 0 12\n",
    "Te 0 12\n",
    "Tf 0 5\n",
    "Th 0 5\n",
    "To2 0 15\n",
    "\"\"\"\n",
    "\n",
    "initParams = {\n",
    "    \"biomassReaction\": \"Growth\",\n",
    "    \"timeStep\": \"0.01\",\n",
    "    \"nSteps\": 300,\n",
    "    \"initBiomass\": \"0.01\"\n",
    "}\n",
    "\n",
    "initConcentrations = {\n",
    "    \"Carbon1\": 10,\n",
    "    \"Carbon2\": 10,\n",
    "    \"Hext\": 0,\n",
    "    \"Dext\": 0,\n",
    "    \"Eext\": 0,\n",
    "    \"Fext\": 0,\n",
    "    \"Oxygen\": 100,\n",
    "    \"Biomass\": 0\n",
    "}\n",
    "\n",
    "initRegStates = {\n",
    "    \"RPcl\": 0,\n",
    "    \"RPO2\": 0,\n",
    "    \"RPb\": 0,\n",
    "    \"RPh\": 0\n",
    "}\n",
    "\n",
    "default = {\n",
    "    \"constraints\": constraints,\n",
    "    \"initParams\": initParams,\n",
    "    \"initConcentrations\": initConcentrations,\n",
    "    \"initRegStates\": initRegStates\n",
    "}\n",
    "\n",
    "simulations = {\n",
    "    \"fig5\": default | {\"initRegStates\": initRegStates | {\n",
    "        \"RPcl\": 1\n",
    "    }}, \n",
    "    \"fig6\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Oxygen\": 5,\n",
    "        \"Hext\":0\n",
    "    }} | {\"initRegStates\": initRegStates},\n",
    "    \"fig7\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Hext\": 2,\n",
    "    }} | {\"initRegStates\": initRegStates | {\n",
    "        \"RPh\": 1\n",
    "    }},\n",
    "    \"fig8\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 0,\n",
    "        \"Carbon2\": 5,\n",
    "        \"Hext\": 10,\n",
    "    }} | {\"initRegStates\": initRegStates | {\n",
    "        \"RPh\": 1\n",
    "    }},\n",
    "    \"fig9\": default | {\"initConcentrations\": initConcentrations | {\n",
    "        \"Carbon1\": 1,\n",
    "        \"Hext\": 5,\n",
    "        \"Fext\": 0.1,\n",
    "    }} | {\"initRegStates\": initRegStates | {\n",
    "        \"RPcl\": 1,\n",
    "        \"RPh\": 1\n",
    "    }},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_constraints(cfg, filename):\n",
    "    with open(filename, \"w\") as fp:\n",
    "        fp.write(cfg[\"constraints\"])\n",
    "\n",
    "def _write_state(fp, sect, state):\n",
    "    fp.write(f\"{sect}\\n\")\n",
    "    for node, value in state.items():\n",
    "        fp.write(f\"{node}\\t{value}\\n\")\n",
    "        \n",
    "def write_initial(cfg, filename):\n",
    "    with open(filename, \"w\") as fp:\n",
    "        for sect, value in cfg[\"initParams\"].items():\n",
    "            fp.write(f\"{sect}\\n{value}\\n\\n\")\n",
    "        sect = \"initConcentrations\"\n",
    "        _write_state(fp, sect, cfg[sect])\n",
    "        sect = \"initRegStates\"\n",
    "        if sect in cfg:\n",
    "            fp.write(\"\\n\")\n",
    "            _write_state(fp, sect, cfg[sect])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def met4j_dynamicRFBA(mn, rn, cons, init, output):\n",
    "    subprocess.run([\"java\", \"-Dlog4j.configuration=\",\n",
    "        \"-jar\", \"/opt/met4j-dynamicRFBA.jar\",\n",
    "        \"-sol\", \"GLPK\", \n",
    "        \"-s\", mn, \"-reg\", rn,\n",
    "        \"-cons\", cons, \"-init\", init, \"-out\", output], check=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_simulation(out):\n",
    "    df = pd.read_csv(out, sep=\"\\t\")\n",
    "    df.set_index(\"Time\", inplace=True)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulations(mn, bnet_file, rn, output_dir, output_tmpl):\n",
    "    output = {}\n",
    "    \n",
    "    bnet = BooleanNetwork.load(bnet_file)\n",
    "    bn2flexflux(bnet, rn, regulators)\n",
    "\n",
    "    for name, cfg in tqdm(simulations.items()):\n",
    "        out = os.path.join(output_dir, output_tmpl.format(name))\n",
    "        cons = os.path.join(output_dir, f\"constraints-{name}.txt\")\n",
    "        init = os.path.join(output_dir, f\"initialParameters-{name}.txt\")\n",
    "        write_constraints(cfg, cons)\n",
    "        write_initial(cfg, init)\n",
    "        met4j_dynamicRFBA(mn, rn, cons, init, out)\n",
    "        output[name] = read_simulation(out)\n",
    "        \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "benchmark_data = []\n",
    "\n",
    "for bnet_file in sorted(glob(bnet_dir)):\n",
    "    solution_name = bnet_file.split('/')[-1].split('.')[0]\n",
    "    sim_path = simulations_dir + solution_name\n",
    "    create_folder(sim_path)\n",
    "\n",
    "    rn = sim_path + '/regulatory_network.sbml'\n",
    "    template = 'output-{}.csv'\n",
    "\n",
    "    run_simulations(metabolic_network_file, bnet_file, rn, sim_path, template)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}