publication_figures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cecb389d-12a4-43ef-8f23-d89d40b98697",
   "metadata": {},
   "source": [
    "Copyright by Paul Rudolph\n",
    "\n",
    "Research Group Applied Systems Biology - Head: Prof. Dr. Marc Thilo Figge\n",
    "\n",
    "https://www.leibniz-hki.de/en/applied-systems-biology.html\n",
    "\n",
    "HKI-Center for Systems Biology of Infection\n",
    "\n",
    "Leibniz Institute for Natural Product Research and Infection Biology - Hans Knöll Insitute (HKI)\n",
    "\n",
    "Adolf-Reichwein-Straße 23, 07745 Jena, Germany\n",
    "\n",
    "This code is licensed under BSD 2-Clause\n",
    "\n",
    "See the LICENSE file provided with this code for the full license."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2efc4c77-2abc-40bd-857b-e1cce9f4d46d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import seaborn as sns\n",
    "import plotly.express as px\n",
    "import plotly.io as pio\n",
    "import plotly.graph_objects as go\n",
    "import matplotlib\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import multiprocessing,functools\n",
    "from tqdm import tqdm\n",
    "import sys\n",
    "from scipy.stats import chi2\n",
    "from SALib.sample import saltelli\n",
    "from SALib.analyze import sobol\n",
    "\n",
    "from model import system_size,adjust_CaL,time_scaling,ldh_model,cal_model,cal_nb_model,candida_model,cal_nb_with_binding_model\n",
    "from fitting import simulate, validation_interval,distance,distance_max,log_distance, get_parameters\n",
    "from model_utils import objective_LDH,objective_CAL,objective_Nb,generate_data_LDH,generate_data_CAL,generate_data_Nb,prepare_data_LDH,prepare_data_CaL,prepare_data_Nb,objective_Candida,prepare_data_Candida\n",
    "from preprocessing import read_data, preprocess_data\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "RESULT_FOLDER = os.path.abspath(f\"../PUBLICATION/SIMULATION/\")\n",
    "synthetic_nb = \"CAL1-H1\"\n",
    "native_nb = \"CAL1-F1\"\n",
    "PLOT_FIT_FOLDER = \"../PUBLICATION/PLOTS/\"\n",
    "os.makedirs(PLOT_FIT_FOLDER,exist_ok=True)\n",
    "DATA_FOLDER = os.path.abspath(\"../DATA/\")\n",
    "file_name = os.path.join(DATA_FOLDER,\"Experimental Data For Modelling.xlsx\")\n",
    "df_ex5,df_ex3,df_ex3_candida,df_ex1,df_fig_cal,df_fig_candida = read_data(file_name)\n",
    "df_data_LDH,df_data_sim,df_data_pre,df_data_pre_co,df_data_post,df_data_candida_sim,df_data_candida_pre,df_data_candida_post = preprocess_data(file_name,time_scaling, adjust_CaL=adjust_CaL, adjust_LDH=True,system_size=system_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df281e7c-6ac4-4d28-a393-3faa823c1fb2",
   "metadata": {},
   "source": [
    "## Set style for plotly plots\n",
    "  * Font: 12px(16pt), Arial\n",
    "  * Ticklabels: 11px(14.66pt),bold,black\n",
    "  * Horizontal legend at the bottom of the figure\n",
    "  * Bold axis lines and bold ticks\n",
    "  * Custom function for line plots with error bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db97b18-3b12-4a7c-9872-47a697a0f23a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"whitegrid\")\n",
    "plotly_template = pio.templates[\"simple_white\"]\n",
    "plotly_template.layout.update(font={\"size\":14+2/3, \"family\":\"Arial\"},\n",
    "                              legend_title=dict(font=dict(size=16)))\n",
    "plotly_template.layout.legend.update(font=dict(size=16),\n",
    "                                     orientation='h',\n",
    "                                     xanchor='center',\n",
    "                                     yanchor='bottom',\n",
    "                                     x=0.5,\n",
    "                                     y=-0.15,\n",
    "                                     itemwidth=40,\n",
    "                                     itemsizing=\"constant\",\n",
    "                                     entrywidth=40)\n",
    "plotly_template.layout.update(xaxis=dict(linecolor=\"Black\",tickprefix=\"<b>\",ticksuffix=\"</b>\",tickangle=0, linewidth=2,tickcolor=\"black\",tickwidth=2),\n",
    "                              font=dict(size=14.66),\n",
    "                              titlefont=dict(size=16,family=\"Arial\"))\n",
    "plotly_template.layout.update(yaxis=dict(linecolor=\"Black\",tickprefix=\"<b>\",ticksuffix=\"</b>\",tickangle=0, linewidth=2,tickcolor=\"black\",tickwidth=2),\n",
    "                              font=dict(size=14.66),\n",
    "                              titlefont=dict(size=16,family=\"Arial\"))\n",
    "pio.templates[\"plotly_custom\"] = plotly_template\n",
    "px.defaults.template = \"plotly_custom\"\n",
    "px.defaults.width = 600\n",
    "px.defaults.height = 400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab955dc6-1de8-4e60-9e3d-657cd462cd77",
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_with_bands(**kwargs):\n",
    "    '''Plotly line plot with error bands'''\n",
    "    figure_with_errok_bars = px.line(**kwargs)\n",
    "    fig = px.line(**{arg: val for arg,val in kwargs.items() if arg != 'error_y'})\n",
    "    for data in figure_with_errok_bars.data:\n",
    "        x = list(data['x'])\n",
    "        y_upper = list(data['y'] + data['error_y']['array'])\n",
    "        y_lower = list(np.maximum(data['y'] - data['error_y']['array'],0.0) if data['error_y']['arrayminus'] is None else data['y'] - data['error_y']['arrayminus'])\n",
    "        color = f\"rgb{matplotlib.colors.to_rgb(data['line']['color'])}\".replace(\"rgb\",\"rgba\").replace(\")\",\",.3)\")\n",
    "        fig.add_trace(\n",
    "            go.Scatter(\n",
    "                x = x+x[::-1],\n",
    "                y = y_upper+y_lower[::-1],\n",
    "                fill = 'toself',\n",
    "                fillcolor = color,\n",
    "                line = dict(\n",
    "                    color = 'rgba(255,255,255,0)'\n",
    "                ),\n",
    "                hoverinfo = \"skip\",\n",
    "                showlegend = False,\n",
    "                legendgroup = data['legendgroup'],\n",
    "                xaxis = data['xaxis'],\n",
    "                yaxis = data['yaxis'],\n",
    "            )\n",
    "        )\n",
    "    reordered_data = []\n",
    "    for i in range(int(len(fig.data)/2)):\n",
    "        reordered_data.append(fig.data[i+int(len(fig.data)/2)])\n",
    "        reordered_data.append(fig.data[i])\n",
    "    fig.data = tuple(reordered_data)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce015dcc-0d41-4d77-9b91-6fdb7d11fd5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_df(RESULT_FOLDER):\n",
    "    '''Load in the fitting results (always the best evaluation)'''\n",
    "    df_list = list()\n",
    "    for result in os.listdir(RESULT_FOLDER):\n",
    "            result_path = os.path.join(RESULT_FOLDER, result)\n",
    "            if os.path.isdir(result_path):\n",
    "                df = {**pd.read_feather(result_path+\"/LDH.feather\").assign(A= int(result)).iloc[0].rename({\"Fun\":\"LDH_ERROR\"}),\n",
    "                                **pd.read_feather(result_path+\"/CAL.feather\").iloc[0].rename({\"Fun\":\"CAL_ERROR\"}),\n",
    "                                **pd.read_feather(result_path+f\"/NB_CAL.feather\").iloc[0].rename({\"Fun\":\"NB_ERROR\"}),\n",
    "                     **pd.read_feather(result_path+f\"/CANDIDA.feather\").iloc[0].rename({\"Fun\":\"CANDIDA_ERROR\"})}\n",
    "                df_list.append(df)\n",
    "    return pd.DataFrame(df_list).sort_values(\"A\").assign(k_c=0.0)\n",
    "\n",
    "def simulate_data_for_dosing(model,args):\n",
    "    '''Simulate model for a specific CaL and Nb concentration'''\n",
    "    initials,kwargs,CaL,Nb = args\n",
    "    column_names = [\"CaL\",\"Nb\",\"CaL8\",\"E_Alive\",\"E_Dead\",\"LDH\"]\n",
    "    time_max = 100*time_scaling\n",
    "    time_points = np.linspace(0,time_max,time_max+1,dtype=\"float64\")\n",
    "    funcptr = model\n",
    "    order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta']\n",
    "    parameters = np.array([kwargs[op] for op in order_parameters])\n",
    "    return pd.DataFrame(simulate(parameters, initials, time_points,funcptr).T,columns = column_names).assign(Time = time_points, CaL = CaL, Nb = Nb, Agg =kwargs[\"A\"])\n",
    "\n",
    "def simulate_dosing(model, df):\n",
    "    '''Screen over nb and cal concentration to find perfect dosing'''\n",
    "    screening_list = list()\n",
    "    for row in df.to_dict(orient='records'):\n",
    "        for CaL in np.linspace(0,120,121):\n",
    "            for Nb in np.linspace(0,70,701):\n",
    "                time_max = 100*time_scaling\n",
    "                time_points = np.linspace(0,time_max,time_max+1,dtype=\"float64\")\n",
    "                initials = np.array([CaL/row[\"A\"]*system_size*1e-3,Nb*system_size*1e-3,0.0,1.0,0.0,0.0])\n",
    "                screening_list.append((initials,row,CaL, Nb))\n",
    "\n",
    "\n",
    "\n",
    "    results = list()\n",
    "    constant_args = [model]\n",
    "    partial_parallelize = functools.partial(simulate_data_for_dosing, *constant_args)\n",
    "    with multiprocessing.Pool(np.minimum(len(screening_list),250)) as pool:\n",
    "        results = list(tqdm(pool.imap(partial_parallelize, screening_list), total=len(screening_list), file=sys.stdout))\n",
    "\n",
    "    df_effect = (pd.concat(results)[[\"E_Dead\",\"Time\",\"CaL\", \"Nb\",\"Agg\"]])\n",
    "\n",
    "    return df_effect\n",
    "\n",
    "def generate_dosing_figure(df):\n",
    "    '''Helper function for generating the final plot'''\n",
    "    data = df.groupby([\"CaL\"])[\"Nb_max\"].agg([\"mean\",\"std\"]).reset_index().assign(std= lambda df_:df_[\"std\"].where(df_[\"std\"]!=0, 0.01)).sort_values(\"CaL\",ascending=False)\n",
    "    fig = line_with_bands(data_frame=data,\n",
    "               x=\"CaL\",y=\"mean\", \n",
    "               labels={\"Effect\":\" Dead VEC[%]\",\"Time\":\"Time [h]\",\"mean\":f\"<b>{synthetic_nb} (µM)</b>\",\"CaL\":\"<b>CaL (µM)</b>\",\"A\":\"Aggregate Size\"}, \n",
    "               error_y=\"std\")\n",
    "    fig.update_yaxes(range=[0,61])\n",
    "    fig.update_xaxes(range=[0,111],tickmode=\"array\",tickvals=[0,10,20,30,40,50,60,70,80,90,100,110])\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "675d4dd6-4737-47be-b1ba-234ecf9bd2dc",
   "metadata": {},
   "source": [
    "# Fig 4a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d49e8ae2-95b1-452b-a18b-b742fce2f973",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_dosing = simulate_dosing(cal_nb_model.address, generate_df(RESULT_FOLDER))\n",
    "df = (pd.concat([df_dosing])\n",
    "     .query(\"E_Dead != 0 and Time == 24\")\n",
    "     .assign(E_Dead = lambda df_:df_.E_Dead.where(df_.E_Dead == 0,df_.eval(\"E_Dead*100\")))\n",
    "     .groupby([\"CaL\",\"Nb\",\"Agg\"])[\"E_Dead\"].agg(\"max\").reset_index()\n",
    "      .query(\"E_Dead < 10 \")\n",
    "     .groupby([\"CaL\",\"Agg\"],group_keys=False).apply(lambda grp_: grp_.assign(Nb_max = grp_.query(f\"Nb == {grp_.Nb.min()}\").sort_values(\"Nb\",ascending=True).Nb.iloc[0])) \n",
    "      [[\"Nb_max\",\"CaL\",\"Agg\"]].drop_duplicates().reset_index(drop=True)\n",
    "    )\n",
    "fig = generate_dosing_figure(df)\n",
    "fig.add_hline(y=50.63,opacity=1.,visible=True,line=dict(width=2,dash=\"longdash\",color=\"rgb(0,0,0)\"),x1=107/111)\n",
    "fig.add_vline(x=107,opacity=1.,visible=True,line=dict(width=2,dash=\"longdash\",color=\"rgb(0,0,0)\"),y1=50.8/61)\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4A.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4A.png\"), format=\"png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9ba6129-84e4-4642-8b33-9306590c2419",
   "metadata": {},
   "source": [
    "# Fig 4b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5477800-9b4c-46cf-be47-aa448cc3f0ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_data(funcptr,initials, time_points,**kwargs):\n",
    "    '''Simulate model with a certain parmaater configuration'''\n",
    "    order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta',\"k_s\"]\n",
    "    parameters = np.array([kwargs[op] for op in order_parameters])\n",
    "    data = simulate(parameters, initials, time_points,funcptr)\n",
    "    return data\n",
    "\n",
    "df_list = list()\n",
    "column_names = [\"CaL\",\"Nb\",\"CaLA\",\"E_Alive\",\"E_Dead\",\"LDH\",\"Candida\",\"F_NI\",\"F_I\"]\n",
    "for funccptr,df in [(candida_model.address,generate_df(RESULT_FOLDER))]:\n",
    "    for row in df.to_dict(orient='records'):\n",
    "        for delay in np.arange(0,24):\n",
    "            row.update(delay=delay)\n",
    "            #row.update(k_b = row[\"k_b_c\"])\n",
    "            candida_coverage= df_data_candida_sim.query(\"Time==72\").Readout.mean()/df.beta.mean() \n",
    "            time_delay = delay*time_scaling\n",
    "            Y0 = 2.0e4\n",
    "            initials = np.array([0.0,0,0.0,candida_coverage,0.0,0.0,Y0,0,0])\n",
    "            if delay != 0.0:\n",
    "                time_points = np.linspace(0,time_delay,time_delay+1,dtype=\"float64\")\n",
    "                post_initials = np.array([0.0,0,0.0,candida_coverage,0.0,0.0,Y0,0,0])\n",
    "                initials = simulate_data(funccptr,post_initials, time_points, **row)[:,-1]\n",
    "                df_list.append((pd.DataFrame(simulate_data(funccptr,post_initials, time_points, **row).T,columns = column_names)\n",
    "                .assign(Delay = delay, \n",
    "                  Time = time_points,  \n",
    "                  initNb = Nb, \n",
    "                  Agg= row[\"A\"],\n",
    "                  E_Alive= lambda df_:df_.E_Alive +(1-candida_coverage),\n",
    "                  E_Dead=lambda df_:df_.E_Dead*100)))\n",
    "            for Nb in [0,4,8,16]:\n",
    "                setup_initials = np.copy(initials)\n",
    "                time_max = 72*time_scaling\n",
    "                time_points = np.linspace(0,time_max,time_max+1,dtype=\"float64\")\n",
    "                setup_initials[1] = Nb*system_size*1e-3\n",
    "                df_list.append(pd.DataFrame(simulate_data(funccptr,setup_initials, time_points, **row).T,columns = column_names)\n",
    "                          .assign(Delay = delay, \n",
    "                                  Time = time_points+delay,\n",
    "                                  secretion= lambda df_: df_.eval(f\"E_Alive*F_I*{row['k_s']}*{row['A']}\"),\n",
    "                                  initNb = Nb, \n",
    "                                  Agg= row[\"A\"],\n",
    "                                  E_Alive= lambda df_:df_.E_Alive +(1-candida_coverage),\n",
    "                                  E_Dead=lambda df_:df_.E_Dead*100)\n",
    "                 )\n",
    "fig = line_with_bands(data_frame=(\n",
    "    pd.concat(df_list).assign(secretion = lambda df_: df_.secretion*1e6)\n",
    "    .query(\"Delay == 0 and Time <= 72 and initNb==0\")\n",
    "    .sort_values(\"Time\")\n",
    "    .groupby([\"initNb\",\"Delay\",\"Agg\"], group_keys=False).apply(lambda grp_: grp_.assign(secretion=lambda df_:df_.secretion.cumsum()))\n",
    "    .groupby([\"initNb\",\"Time\",\"Delay\"])[\"secretion\"].agg([\"mean\",\"std\"]).reset_index()\n",
    "), \n",
    "           x=\"Time\",y=\"mean\",  error_y=\"std\", labels={\"mean\":'<b>\"Effective\" CaL (µM)<br>during infection</b>',\"Time\":\"<b>Time (h)</b>\",\"initNb\":\"Initial CAL1-F1 [µM]\", \"Delay\":\"Post Treatment [h]\"})\n",
    "fig.update_xaxes(\n",
    "    range=[0,24.3],\n",
    "    tickmode='array',\n",
    "    tickvals=[0,8,16,24,32,40,48],  # Use all time points as tick values\n",
    ")\n",
    "fig.update_yaxes(range=[0,123])\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4B.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4B.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "450b552a-9f00-4f80-81f0-ce9336929816",
   "metadata": {},
   "source": [
    "# Fig 4c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4317446b-3710-498d-b526-6513f8cabae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_list = list()\n",
    "column_names = [\"CaL\",\"Nb\",\"CaLA\",\"E_Alive\",\"E_Dead\",\"LDH\",\"Candida\",\"F_NI\",\"F_I\"]\n",
    "for funcptr,df in [(candida_model.address,generate_df(RESULT_FOLDER))]:\n",
    "    for row in df.to_dict(orient='records'):\n",
    "        for delay in np.arange(0,24):\n",
    "            row.update(delay=delay)\n",
    "            candida_coverage= df_data_candida_sim.query(\"Time==72\").Readout.mean()/df.beta.mean()\n",
    "            time_delay = delay*time_scaling\n",
    "            Y0 = 2.0e4\n",
    "            initials = np.array([0.0,0,0.0,candida_coverage,0.0,0.0,Y0,0,0])\n",
    "\n",
    "            for Nb in [0,4,8,16]:\n",
    "                if delay != 0.0:\n",
    "                    time_points = np.linspace(0,time_delay,time_delay*10+1,dtype=\"float64\")\n",
    "                    post_initials = np.array([0.0,0,0.0,candida_coverage,0.0,0.0,Y0,0,0])\n",
    "                    initials = simulate_data(funcptr,post_initials, time_points, **row)[:,-1]\n",
    "                    df_list.append((pd.DataFrame(simulate_data(funcptr,post_initials, time_points, **row).T,columns = column_names)\n",
    "                    .assign(Delay = delay, \n",
    "                      Time = time_points,  \n",
    "                      initNb = Nb, \n",
    "                      Agg= row[\"A\"],\n",
    "                      E_Alive= lambda df_:df_.E_Alive +(1-candida_coverage),\n",
    "                      E_Dead=lambda df_:df_.E_Dead*100)))\n",
    "                setup_initials = np.copy(initials)\n",
    "                time_max = 72*time_scaling\n",
    "                time_points = np.linspace(0,time_max,time_max*10+1,dtype=\"float64\")\n",
    "                setup_initials[1] = Nb*system_size*1e-3\n",
    "                df_list.append(pd.DataFrame(simulate_data(funcptr,setup_initials, time_points, **row).T,columns = column_names)\n",
    "                          .assign(Delay = delay, \n",
    "                                  Time = time_points+delay,  \n",
    "                                  initNb = Nb, \n",
    "                                  Agg= row[\"A\"],\n",
    "                                  E_Alive= lambda df_:df_.E_Alive +(1-candida_coverage),\n",
    "                                  E_Dead=lambda df_:df_.E_Dead*100)\n",
    "                 )\n",
    "times = np.arange(0,24)\n",
    "delay_array= [0,3,9,12]\n",
    "\n",
    "data = pd.concat(df_list).assign(LDH = lambda df_:df_.LDH *1e3).query(\"Time <= 48\").groupby([\"initNb\",\"Time\",\"Delay\"])[\"LDH\"].agg([\"mean\",\"std\"]).reset_index()\n",
    "time_10 = data.query('mean <= 10 and Delay ==0 and initNb==0').sort_values('Time').iloc[-1].Time             \n",
    "fig = line_with_bands(data_frame=data.query(f\"Delay in {delay_array}\").assign(initNb= lambda df_:\"<b>\"+df_.initNb.map(str)+\" µM</b>\\t\"), x=\"Time\",y=\"mean\",color=\"initNb\", error_y=\"std\",\n",
    "                      width=1200,height=800, facet_col=\"Delay\",\n",
    "                      facet_col_wrap=2,facet_row_spacing=.2,facet_col_spacing=.12,\n",
    "           labels={\"mean\":\" LDH (ng/ml)\",\"Time\":\"Time (h)\",\"initNb\":f\"<b>{native_nb}:\\t</b>\",\"Delay\":\"post-treatment\"})\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \" \".join(anno.text.__add__(\" h\").split(\"=\")[::-1])))\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"Simultaneous addition\" if anno.text.split(\" \")[0] == \"0\" else anno.text))\n",
    "\n",
    "fig.add_hline(y=200,opacity=1.,visible=True,line=dict(width=3,dash=\"dash\",color=\"rgba(255,0,0)\"))\n",
    "fig.update_xaxes(\n",
    "    range=[0,24.3],showticklabels=True,\n",
    "    tickmode='array',\n",
    "    tickvals=np.arange(0,48,2)  # Use all time points as tick values\n",
    ")\n",
    "\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = \"<b>LDH (ng/ml)</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Time (h)</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+\"</b>\"))\n",
    "\n",
    "fig.update_yaxes(range=[0,404],showticklabels=True)\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4C.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_4C.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f713fa-9c89-41fa-b1e7-ff8cf8d9d441",
   "metadata": {},
   "source": [
    "# Fig S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbfe8f6c-84c5-4ef1-9813-a5b3f188351e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = generate_df(RESULT_FOLDER)\n",
    "df_liste = list()\n",
    "df_scatter = list()\n",
    "column_names = [\"CaL\",\"Nb\",\"CaLA\",\"E_Alive\",\"E_Dead\",\"LDH\"]\n",
    "def simulate_data(initials, time_points,**kwargs):\n",
    "    funcptr = cal_nb_model.address\n",
    "    order_parameters = ['k_a', 'k_c', 'k_n','k_l',\"k_b\", 'k_d','A', 'alpha', 'beta']\n",
    "    parameters = np.array([kwargs[op] for op in order_parameters])\n",
    "    return simulate(parameters, initials, time_points,funcptr)\n",
    "for row in df.to_dict(orient='records'):\n",
    "    for CaL in [1,10,35,70]:\n",
    "        for Nb in [0,4,8,16]:\n",
    "            time_max = 72*time_scaling\n",
    "            time_points = np.linspace(0,time_max,time_max+1,dtype=\"float64\")\n",
    "            initials = np.array([CaL/row[\"A\"]*system_size*1e-3,Nb*system_size*1e-3,0.0,1.0,0.0,0.0])\n",
    "            df_liste.append(pd.DataFrame(simulate_data(initials, time_points, **row).T,columns = column_names)\n",
    "                      .assign(Time = time_points, CaL = CaL, Nb = Nb, Agg =row[\"A\"])\n",
    "                 )\n",
    "column_names = [\"CaL\",\"Nb\",\"CaLA\",\"E_Alive\",\"E_Dead\",\"LDH\",\"_1\",\"_2\",\"_3\"]\n",
    "def simulate_data(initials, time_points,**kwargs):\n",
    "    funcptr = candida_model.address\n",
    "    order_parameters = ['k_a', 'k_c', 'k_n','k_l',\"k_b\", 'k_d','A', 'alpha', 'beta',\"k_s\"]\n",
    "    parameters = np.array([kwargs[op] for op in order_parameters])\n",
    "    return simulate(parameters, initials, time_points,funcptr)\n",
    "for row in df.to_dict(orient='records'): \n",
    "    for Nb in [0,4,8,16]:\n",
    "        time_max = 72*time_scaling\n",
    "        time_points = np.linspace(0,time_max,time_max+1,dtype=\"float64\")\n",
    "        initials = np.array([0,Nb*system_size*1e-3,0.0,df_data_candida_sim.query(\"Time==72\").Readout.mean()/df.beta.mean(),0.0,0.0,20000,0,0])\n",
    "        df_liste.append(pd.DataFrame(simulate_data(initials, time_points, **row).T,columns = column_names)\n",
    "                  .assign(Time = time_points, CaL = \"Candida\", Nb = Nb, Agg =row[\"A\"])\n",
    "             )            \n",
    "\n",
    "data_list=list()\n",
    "for  key, row_s in df_data_sim.query(\"Nb == 0\").groupby([\"CaL\",\"Nb\"],group_keys=False):\n",
    "    initial_CaL = key[0]*adjust_CaL*1/system_size*1e3\n",
    "    data_list.append(row_s.reset_index().drop(columns=[\"CaL\",\"Nb\"]).assign(CaL = initial_CaL, Nb = 0))\n",
    "for  key, row_s in df_data_candida_sim.query(\"Nb == 0\").groupby([\"Nb\"],group_keys=False):\n",
    "    data_list.append(row_s.reset_index().drop(columns=[\"Nb\"]).assign(CaL = \"Candida\", Nb = 0))\n",
    "df_scatter = pd.concat(data_list).sort_values([\"CaL\",\"Nb\"])\n",
    "df_effect = (pd.concat(df_liste)\n",
    "             .assign(LDH = lambda df_:df_.LDH*1e3)\n",
    "             .groupby([\"CaL\",\"Nb\",\"Time\"])[\"LDH\"].agg([\"mean\",\"std\"])\n",
    "             .reset_index()\n",
    "            )\n",
    "fig = line_with_bands(data_frame=df_effect.assign(Nb= lambda df_:\"<b>\"+df_.Nb.map(str)+\" µM</b>\\t\"), x=\"Time\",y=\"mean\",facet_col=\"CaL\",color=\"Nb\", error_y=\"std\",width=1200,height=800, \n",
    "           labels={\"mean\":\"LDH (ng/ml)\",\"Time\":\"Time (h)\",\"Nb\":f\"<b>{synthetic_nb}/{native_nb} concentration:</b>\\t\"},facet_col_spacing=0.1,facet_row_spacing=0.2,facet_col_wrap=4)\n",
    "fig.update_xaxes(\n",
    "    matches=None,\n",
    "    showticklabels=True,\n",
    "    range=[0,74],\n",
    "    tickmode='array',\n",
    "    tickvals=[0,24,48,72],  # Use all time points as tick values\n",
    ")\n",
    "\n",
    "fig.update_yaxes(\n",
    "    matches=None,\n",
    "    showticklabels=True,\n",
    "    tickmode='array',\n",
    "    tickvals=[0,500,1000,1500],  range=[0,1800]\n",
    ")\n",
    "\n",
    "\n",
    "scatter= px.scatter(df_scatter.assign(LDH = lambda df_:df_.Readout *1e3),x=\"Time\",y=\"LDH\",facet_col=\"CaL\",facet_col_wrap=4)\n",
    "scatter.for_each_trace(lambda trace:trace.update(name=\"DATA\",marker=dict(size=10,color=\"black\"),showlegend=False,legendgroup=\"DATA\"))\n",
    "for i in range(len(scatter.data)):\n",
    "    fig.add_trace(scatter.data[i])\n",
    "    \n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = \"<b>LDH (ng/ml)</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Time (h)</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \" \".join(anno.text.__add__(\" µM\").split(\"=\")[::-1])  if \"Candida\" != anno.text.split(\"=\")[1] else \"<i>C. albicans</i>-secreted CaL\"))\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+\"</b>\"))\n",
    "fig.update_layout(legend=dict(y=.4,x=.6))\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_S3.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"FIG_S3.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2ed687-1a67-4d5f-b603-9c3289758542",
   "metadata": {},
   "source": [
    "## Give parameter values by aggregate size in correct units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4e3e657-f689-4e08-b2cd-c598cfcc8bd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_dependent_variables = ['k_b', 'k_a', 'k_d', 'k_l', 'k_s']\n",
    "\n",
    "parameter_table = (generate_df(RESULT_FOLDER)\n",
    " [time_dependent_variables+[\"A\",\"alpha\",\"beta\"]]\n",
    ".melt(id_vars=\"A\", var_name=\"Parameter\")\n",
    " .assign(value=lambda df_:df_.value.where(df_.Parameter.apply(lambda x:  x not in time_dependent_variables),df_.value/3600*time_scaling))\n",
    " .rename(columns={\"value\":\"Parameter Value\",\"A\":\"Aggregate Size\"})\n",
    " .pivot(index='Parameter', columns='Aggregate Size')\n",
    ")\n",
    "display(parameter_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42a0d136-0ec1-41bb-a5da-3505dd5aaa4c",
   "metadata": {},
   "source": [
    "# Identifiability analysis using Profile Likelihood Methode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a696e51-d054-44f7-98c8-0b89f65b002e",
   "metadata": {},
   "outputs": [],
   "source": [
    "AGG_SIZE = 8\n",
    "all_parameters = generate_df(RESULT_FOLDER).query(f\"A == {AGG_SIZE}\").iloc[0].to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5845859b-8112-4ab4-8817-469557a52013",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "CRITICAL_VALUE = chi2.ppf(0.95, 1)\n",
    "time_dependent_variables = ['k_b', 'k_a', 'k_d', 'k_l', 'k_s']\n",
    "AGG = \"mean\"\n",
    "order_parameters = ['k_l', 'beta']\n",
    "fitting_parameters = ['k_l', 'beta']\n",
    "list_parameters = [all_parameters[i] for i in order_parameters]\n",
    "no_fitting = {}\n",
    "data = prepare_data_LDH(df_data_LDH)\n",
    "\n",
    "\n",
    "def obj_wrapper_LDH(*parameters):\n",
    "    SSE = objective_LDH(ldh_model.address, data, log_distance, *parameters)\n",
    "    n_samples = len(data[0][0])\n",
    "    return (n_samples/2)*np.log(SSE)\n",
    "\n",
    "\n",
    "ldh_ci = (validation_interval(obj_wrapper_LDH, CRITICAL_VALUE,fitting_parameters, no_fitting, order_parameters, *list_parameters)\n",
    "         )\n",
    "\n",
    "order_parameters = ['k_a', 'k_c', 'k_l', 'k_d', 'alpha', 'beta']\n",
    "fitting_parameters = ['k_a',  'k_d', 'alpha']\n",
    "list_parameters = [all_parameters[i] for i in order_parameters]\n",
    "no_fitting = {key:all_parameters[key] for key in ['k_l', 'beta','k_c']}\n",
    "data = prepare_data_CaL(df_data_sim.query(\"Nb == 0.0\"), all_parameters[\"beta\"],agg=AGG)\n",
    "\n",
    "def obj_wrapper_CaL(*parameters):\n",
    "    SSE = objective_CAL(cal_model.address, data, distance, distance_max,*parameters)\n",
    "    n_samples = np.sum([len(i[1]) for i in prepare_data_CaL(df_data_sim.query(\"Nb == 0.0\"), all_parameters[\"beta\"],agg=AGG)])\n",
    "    return (n_samples/2)*np.log(SSE)\n",
    "\n",
    "\n",
    "cal_ci = (validation_interval(obj_wrapper_CaL, CRITICAL_VALUE, fitting_parameters, no_fitting, order_parameters, *list_parameters)\n",
    "         )\n",
    "\n",
    "\n",
    "order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta']\n",
    "fitting_parameters = ['k_b']\n",
    "list_parameters = [all_parameters[i] for i in order_parameters]\n",
    "no_fitting = {key:all_parameters[key] for key in  ['k_a', 'k_c', 'k_l', 'k_d', 'k_n', 'A','alpha', 'beta']}\n",
    "data = prepare_data_Nb(df_data_sim.query(\"Nb != 0.0\"),df_data_pre.query(\"Nb != 0.0\"),df_data_pre_co.query(\"Nb != 0.0\"),df_data_post.query(\"Nb != 0.0\"), all_parameters[\"beta\"],agg=AGG)\n",
    "\n",
    "\n",
    "def obj_wrapper_Nb(*parameters):\n",
    "    SSE = objective_Nb(cal_nb_model.address, data, distance, distance_max,*parameters)\n",
    "    n_samples = np.sum([len(i[-3]) for i in data[\"sim\"]]+[len(i[-3]) for i in data[\"pre\"]]+[len(i[-3]) for i in data[\"pre_co\"]]+[len(i[-3]) for i in data[\"post\"]])\n",
    "    return (n_samples/2)*np.log(SSE)\n",
    "\n",
    "\n",
    "nb_ci = (validation_interval(obj_wrapper_Nb, CRITICAL_VALUE, fitting_parameters, no_fitting,order_parameters, *list_parameters)\n",
    "        )\n",
    "all_parameters = generate_df(RESULT_FOLDER).query(\"A == 8\").iloc[0].to_dict()\n",
    "order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta','k_s']\n",
    "fitting_parameters = ['k_s']\n",
    "list_parameters = [all_parameters[i] for i in order_parameters]\n",
    "no_fitting = {key:all_parameters[key] for key in list(set(order_parameters)-set(fitting_parameters))}\n",
    "data = prepare_data_Candida(df_data_candida_sim,df_data_candida_pre,df_data_candida_post,all_parameters[\"beta\"])\n",
    "\n",
    "\n",
    "def obj_wrapper_candida(*parameters):\n",
    "    SSE = objective_Candida(candida_model.address, data, distance, distance_max,*parameters)\n",
    "    n_samples = np.sum([len(i[-3]) for i in data[\"sim\"]]+[len(i[-3]) for i in data[\"pre\"]]+[len(i[-3]) for i in data[\"post\"]])\n",
    "    return (n_samples/2)*np.log(SSE)\n",
    "\n",
    "\n",
    "candida_ci = (validation_interval(obj_wrapper_candida, CRITICAL_VALUE, fitting_parameters, no_fitting,order_parameters, *list_parameters)\n",
    "              \n",
    "             )\n",
    "ci_s =(pd.concat([ldh_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  cal_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  nb_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  candida_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"}))])\n",
    "  .reset_index().drop(columns=[\"level_1\"])\n",
    " .pipe(lambda df_: df_.assign(**{key:df_[key].where(df_.Parameter.apply(lambda x:  x not in time_dependent_variables),df_[key]/3600*time_scaling) for key in [\"upper\",\"lower\",\"Mean\"]}))\n",
    "       .assign(Parameter= lambda df_:\"$\"+df_.Parameter.replace(\"alpha\",\"alpha\").replace(\"beta\",\"beta\")+\"$\")\n",
    " [[\"Parameter\",\"Mean\",\"lower\",\"upper\"]]\n",
    "      .query(\"Parameter != 'k_c' and Parameter != 'k_n'\")\n",
    "      .set_index(\"Parameter\")\n",
    "     )\n",
    "fig = px.scatter(ci_s.reset_index(), y=\"Parameter\", x=\"Mean\",log_x=True,\n",
    "                 error_x=ci_s[\"upper\"] - ci_s[\"Mean\"],\n",
    "                 error_x_minus=ci_s[\"Mean\"] - ci_s[\"lower\"],width=1000,height=400)\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = f\"<b>Parameter</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Mean Estimate (CI 95%)</b>\"\n",
    "fig.update_yaxes(showgrid=True,gridcolor=\"darkgray\",gridwidth=1.8)\n",
    "fig.update_xaxes(range=[-15,3],tickmode='array',gridcolor=\"darkgray\",showgrid=True,gridwidth=1.8,\n",
    "    tickvals=10**np.arange(-15,3,1).astype(float),)\n",
    "fig.update_traces(marker=dict(size=12),error_x=dict(\n",
    "    thickness=5,  \n",
    "    width=5      \n",
    "))\n",
    "\n",
    "\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_CI.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_CI.png\"), format=\"png\")\n",
    "display(ci_s)\n",
    "\n",
    "fig = px.scatter((pd.concat([ldh_ci,cal_ci,nb_ci,candida_ci])\n",
    "                  .assign(TestValue = lambda df_:df_[\"Test Value\"].where(df_.Parameter.apply(lambda x:  x not in time_dependent_variables),df_[\"Test Value\"]/3600*time_scaling),\n",
    "                          Parameter = lambda df_:\"$\"+df_.Parameter.map(str)+\"$\",\n",
    "                          CI = lambda df_:100*chi2.cdf(df_[\"Test Statistic\"], 1))),\n",
    "                 x=\"TestValue\",y=\"Test Statistic\",color=\"Parameter\",log_x=True,width=1000,height=400)\n",
    "fig.update_xaxes(range=[-15,3],tickmode='array',tickvals=10**np.arange(-15,3,1).astype(float),)\n",
    "fig.update_yaxes(range=[0,8])\n",
    "fig.update_layout(legend=dict(orientation=\"v\",x=1.0,xanchor=\"left\",yanchor=\"top\",y=1,title=\"<b>Parameter:</b>\"))\n",
    "\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = f\"<b>Test Statistic</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Parameter Value</b>\"\n",
    "\n",
    "fig.add_hline(y=CRITICAL_VALUE,opacity=1.,visible=True,line=dict(width=3,dash=\"dash\",color=\"rgba(255,0,0)\"))\n",
    "fig.add_annotation(\n",
    "    x=-10,  # X-coordinate of the annotation\n",
    "    y=CRITICAL_VALUE + 0.4,  # Adjust the Y-coordinate to position the text\n",
    "    text=\"<b>95% Confidence</b>\",  # The text you want to display\n",
    "    showarrow=False,  # Hide the arrow\n",
    ")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_CI_INTERVALS.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_CI_INTERVALS.png\"), format=\"png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec40c5da-83a8-4597-8e69-d88fd5bfd93a",
   "metadata": {},
   "source": [
    "## Sensitivity analysis using Sobol variance based methode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff16f20c-2f50-497f-8cda-533eec040c63",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fitting_parameters = ['k_a', 'k_l','k_b', 'k_d', 'alpha', 'beta']\n",
    "order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta']\n",
    "no_fitting_dict = {key:all_parameters[key] for key in ['A','k_n','k_c']}\n",
    "\n",
    "boundaries = (pd.concat([ldh_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  cal_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  nb_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  candida_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"}))])\n",
    "              .reset_index().drop(columns=[\"level_1\"])\n",
    "              [[\"Parameter\",\"Mean\",\"lower\",\"upper\"]]\n",
    "              .assign(Bound= lambda df_: df_.apply(lambda x: (x[\"lower\"],x[\"upper\"]),axis=1))\n",
    "             )\n",
    "\n",
    "ordered_boundaries = np.array([boundaries.query(f\"Parameter == '{key}'\").Bound.iloc[0] for key in fitting_parameters], dtype=\"float64\")\n",
    "problem = {\n",
    "  'num_vars': len(fitting_parameters),\n",
    "  'names': fitting_parameters,\n",
    "  'bounds': ordered_boundaries\n",
    "}\n",
    "\n",
    "def f(x,times,cal,nb,col =-1):\n",
    "    CaL_initial = cal/8*system_size*1e-3\n",
    "    Nb_initial = nb*system_size*1e-3\n",
    "    time_points = np.linspace(0,np.max(times),np.max(times)+1,dtype=\"float64\")\n",
    "    parameters = np.array(get_parameters(x, fitting_parameters, no_fitting_dict,order_parameters),dtype=\"float64\")\n",
    "    initials = np.array([CaL_initial,Nb_initial,0.0,1.0,0.0,0.0])\n",
    "    time_points = np.linspace(0,np.max(times),np.max(times)+1)\n",
    "    simulated_data = simulate(parameters, initials, time_points, cal_nb_model.address)[col,times]\n",
    "    return simulated_data\n",
    "\n",
    "SAMPLES = 2**12\n",
    "df_list = list()\n",
    "for CaL in [70]:\n",
    "    for Nb in [0,16]:\n",
    "        times = np.array([3,24,48])* time_scaling\n",
    "        param_values = saltelli.sample(problem, SAMPLES)\n",
    "        Y = np.vstack([f(params,times,CaL,Nb) for params in param_values])\n",
    "        for i in range(len(times)):\n",
    "            Si = sobol.analyze(problem, Y[:,i], print_to_console=False, parallel=False)\n",
    "            df_list.append(pd.DataFrame(fitting_parameters,columns=[\"Parameter\"]).assign(**{\"<b>First Order</b>\" : Si[\"S1\"].flatten(), \n",
    "                                                                  \"<b>Total Order</b>\" : Si[\"ST\"].flatten()},\n",
    "                                                                 CaL = CaL,\n",
    "                                                                 Nb = Nb,\n",
    "                                                                 Time = times[i]))\n",
    "            \n",
    "df = pd.concat(df_list).melt(id_vars=[\"CaL\",\"Nb\",\"Time\",\"Parameter\"])\n",
    "fig = px.bar(df.query(f\"CaL == 70 and Parameter not in {['beta','k_l']}\").assign(Parameter= lambda df_:\"$\"+df_.Parameter.replace({\"beta\":\"β\",\"alpha\":\"α\"})+\"$\"),color=\"variable\",y=\"Parameter\",x=\"value\", barmode='group',facet_col=\"Time\",facet_row=\"Nb\",                      width=1200,height=800, \n",
    "                      facet_col_wrap=2,facet_row_spacing=.15,facet_col_spacing=.12,\n",
    "             labels={\"variable\":\"<b>Sobol Sensitivity:</b>\",\"value\":\"Index\"})\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = anno.text.replace(\"Nb=\",\"\") + f\" µM {synthetic_nb}; 70 µM CaL\" if anno.text.split(\"=\")[0] == \"Nb\" else anno.text.replace(\"Time=\",\"Simultaneous addition; Time point: \") + \" h\"))\n",
    "df.reset_index(drop=True).to_feather(os.path.join(PLOT_FIT_FOLDER,f\"sobol_indicies.feather\"))\n",
    "\n",
    "fig.update_xaxes(\n",
    "    range=[0,1.1],showticklabels=True,\n",
    ")\n",
    "\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = \"<b>Parameter</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Sobol Index</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+\"</b>\"))\n",
    "\n",
    "fig.update_yaxes(showticklabels=True,range=[-1,4])\n",
    "fig.show()\n",
    "\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SI_SYNTHETIC.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SI_SYNTHETIC.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b162f6-426d-41cf-bcd2-74e63e9b3a10",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fitting_parameters = ['k_a', 'k_l','k_b', 'k_d', 'alpha', 'beta','k_s']\n",
    "order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta','k_s']\n",
    "no_fitting_dict = {key:all_parameters[key] for key in ['A','k_n','k_c']}\n",
    "\n",
    "boundaries = (pd.concat([ldh_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  cal_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  nb_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"})),\n",
    "                  candida_ci.query(\"Rejection == False\").groupby(\"Parameter\").apply(lambda grp_:pd.DataFrame((grp_[\"Test Value\"].max(), grp_[\"Test Value\"].min(), grp_[\"Mean\"].iloc[0])).T.rename(columns = {0:\"upper\",1:\"lower\",2:\"Mean\"}))])\n",
    "              .reset_index().drop(columns=[\"level_1\"])\n",
    "              [[\"Parameter\",\"Mean\",\"lower\",\"upper\"]]\n",
    "              .assign(Bound= lambda df_: df_.apply(lambda x: (x[\"lower\"],x[\"upper\"]),axis=1))\n",
    "             )\n",
    "\n",
    "ordered_boundaries = np.array([boundaries.query(f\"Parameter == '{key}'\").Bound.iloc[0] for key in fitting_parameters], dtype=\"float64\")\n",
    "problem = {\n",
    "  'num_vars': len(fitting_parameters),\n",
    "  'names': fitting_parameters,\n",
    "  'bounds': ordered_boundaries\n",
    "}\n",
    "\n",
    "def f(x,times,nb,col =-4):\n",
    "    candida_coverage = df_data_candida_sim.query(\"Time==72\").Readout.mean()/all_parameters[\"beta\"]\n",
    "    Nb_initial = nb*system_size*1e-3\n",
    "    time_points = np.linspace(0,np.max(times),np.max(times)+1,dtype=\"float64\")\n",
    "    parameters = np.array(get_parameters(x, fitting_parameters, no_fitting_dict,order_parameters),dtype=\"float64\")\n",
    "    initials = np.array([0.0,Nb_initial,0.0,candida_coverage,0.0,0.0,2e4,0.0,0.0])\n",
    "    time_points = np.linspace(0,np.max(times),np.max(times)+1)\n",
    "    simulated_data = simulate(parameters, initials, time_points, candida_model.address)[col,times]\n",
    "    return simulated_data\n",
    "\n",
    "SAMPLES = 2**12\n",
    "df_list = list()\n",
    "for Nb in [0,16]:\n",
    "    times = np.array([3,24,48])* time_scaling\n",
    "    param_values = saltelli.sample(problem, SAMPLES)\n",
    "    Y = np.vstack([f(params,times,Nb) for params in param_values])\n",
    "    for i in range(len(times)):\n",
    "        Si = sobol.analyze(problem, Y[:,i], print_to_console=False, parallel=False)\n",
    "        df_list.append(pd.DataFrame(fitting_parameters,columns=[\"Parameter\"]).assign(**{\"<b>First Order</b>\" : Si[\"S1\"].flatten(), \n",
    "                                                                  \"<b>Total Order</b>\" : Si[\"ST\"].flatten()},\n",
    "                                                             Nb = Nb,\n",
    "                                                             Time = times[i]))\n",
    "\n",
    "df = pd.concat(df_list).melt(id_vars=[\"Nb\",\"Time\",\"Parameter\"])\n",
    "fig = px.bar(df.query(f\"Parameter not in {['beta','k_l']}\").assign(Parameter= lambda df_:\"$\"+df_.Parameter.replace({\"beta\":\"β\",\"alpha\":\"α\"})+\"$\"),color=\"variable\",y=\"Parameter\",x=\"value\", barmode='group',facet_col=\"Time\",facet_row=\"Nb\",                      width=1200,height=800, \n",
    "                      facet_col_wrap=2,facet_row_spacing=.15,facet_col_spacing=.12,\n",
    "             labels={\"variable\":\"<b>Sobol Sensitivity:</b>\",\"value\":\"Index\"})\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = anno.text.replace(\"Nb=\",\"\") + f\" µM {native_nb}\" if anno.text.split(\"=\")[0] == \"Nb\" else anno.text.replace(\"Time=\",\"Simultaneous addition; Time point: \") + \" h\"))\n",
    "df.reset_index(drop=True).to_feather(os.path.join(PLOT_FIT_FOLDER,f\"sobol_indicies.feather\"))\n",
    "\n",
    "fig.update_xaxes(\n",
    "    range=[0,1.1],showticklabels=True,\n",
    ")\n",
    "\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = \"<b>Parameter</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Sobol Index</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+\"</b>\"))\n",
    "\n",
    "fig.update_yaxes(showticklabels=True,range=[-1,5])\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SI_NATIVE.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SI_NATIVE.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d16014e9-0483-4fea-9cb0-2eaa0d07df6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_data(funcptr,initials, time_points,**kwargs):\n",
    "    '''Simulate model with a certain parmaater configuration'''\n",
    "    order_parameters = ['k_a', 'k_c', 'k_n','k_l','k_b', 'k_d','A', 'alpha', 'beta']\n",
    "    parameters = np.array([kwargs[op] for op in order_parameters])\n",
    "    data = simulate(parameters, initials, time_points,funcptr)\n",
    "    return data\n",
    "CAL_CONC = 70\n",
    "df_list = list()\n",
    "column_names = [\"CaL\",\"Nb\",\"CaLA\",\"E_Alive\",\"E_Dead\",\"LDH\", \"NCAL\",\"NCALA\"]\n",
    "for funcptr,df in [(cal_nb_with_binding_model.address,generate_df(RESULT_FOLDER))]:\n",
    "    for row in df.to_dict(orient='records'):\n",
    "        for delay in np.arange(0,24):\n",
    "            row.update(delay=delay)\n",
    "            candida_coverage= df_data_candida_sim.query(\"Time==72\").Readout.mean()/df.beta.mean()\n",
    "            time_delay = delay*time_scaling\n",
    "            initials = np.array([CAL_CONC/row[\"A\"]*1e-3*system_size,0,0.0,1.0,0.0,0.0,0.0,0.0])\n",
    "\n",
    "            for Nb in [2,4,8,16]:\n",
    "                if delay != 0.0:\n",
    "                    time_points = np.linspace(0,time_delay,time_delay*10+1,dtype=\"float64\")\n",
    "                    post_initials = np.array([CAL_CONC/row[\"A\"]*1e-3*system_size,0,0.0,1.0,0.0,0.0,0.0,0.0])\n",
    "                    initials = simulate_data(funcptr,post_initials, time_points, **row)[:,-1]\n",
    "                    df_list.append((pd.DataFrame(simulate_data(funcptr,post_initials, time_points, **row).T,columns = column_names)\n",
    "                    .assign(Delay = delay, \n",
    "                      Time = time_points,  \n",
    "                      initNb = Nb, \n",
    "                      Agg= row[\"A\"])))\n",
    "                setup_initials = np.copy(initials)\n",
    "                time_max = 72*time_scaling\n",
    "                time_points = np.linspace(0,time_max,time_max*10+1,dtype=\"float64\")\n",
    "                setup_initials[1] = Nb*system_size*1e-3\n",
    "                df_list.append(pd.DataFrame(simulate_data(funcptr,setup_initials, time_points, **row).T,columns = column_names)\n",
    "                          .assign(Delay = delay, \n",
    "                                  Time = time_points+delay,  \n",
    "                                  initNb = Nb, \n",
    "                                  Agg= row[\"A\"])\n",
    "                 )\n",
    "times = np.arange(0,24)\n",
    "delay_array= [0,3,9,12]\n",
    "\n",
    "data = pd.concat(df_list).assign(Nb = lambda df_:df_.Nb *1/(system_size*1e-3)).query(\"Time <= 48\")\n",
    "fig = px.line(data_frame=data.query(f\"Delay == 0\").assign(Agg = lambda df_:\"<b>\"+df_.Agg.map(int).map(str)+\"</b>\"), x=\"Time\",y=\"Nb\",color=\"Agg\", \n",
    "                      width=1200,height=800, facet_col=\"initNb\",\n",
    "                      facet_col_wrap=2,facet_row_spacing=.2,facet_col_spacing=.12,\n",
    "           labels={\"Time\":\"Time (h)\",\"initNb\":f\"<b>{synthetic_nb}</b>\",\"Agg\":\"<b>Aggregate Size</b>\"})\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \" \".join(anno.text.__add__(\" µM\").split(\"=\"))))\n",
    "\n",
    "fig.update_xaxes(\n",
    "    range=[0,24.3],showticklabels=True,\n",
    "    tickmode='array',\n",
    "    tickvals=np.arange(0,48,2)  # Use all time points as tick values\n",
    ")\n",
    "\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = f\"<b>{synthetic_nb} (µM)</b>\"\n",
    "\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Time (h)</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+f\"; CaL {CAL_CONC} µM</b>\"))\n",
    "\n",
    "fig.update_yaxes(showticklabels=True, matches=None,range=[0,16.2])\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SYNTHETIC_NB_AGG.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_SYNTHETIC_NB_AGG.png\"), format=\"png\")\n",
    "\n",
    "data = pd.concat(df_list).assign(NCALA = lambda df_:(df_.NCALA+df_.NCAL)/CAL_CONC*100 *1/(system_size*1e-3)*df_.Agg).query(\"Time <= 48\")\n",
    "fig = px.line(data_frame=data.query(f\"Delay == 0\").assign(Agg = lambda df_:\"<b>\"+df_.Agg.map(int).map(str)+\"</b>\"), x=\"Time\",y=\"NCALA\",color=\"Agg\", \n",
    "                      width=1200,height=800, facet_col=\"initNb\",\n",
    "                      facet_col_wrap=2,facet_row_spacing=.2,facet_col_spacing=.12,\n",
    "           labels={\"Time\":\"Time (h)\",\"initNb\":f\"<b>{synthetic_nb}</b>\",\"Agg\":\"<b>Aggregate Size</b>\"})\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \" \".join(anno.text.__add__(\" µM\").split(\"=\"))))\n",
    "\n",
    "fig.update_xaxes(\n",
    "    range=[0,24.3],showticklabels=True,\n",
    "    tickmode='array',\n",
    "    tickvals=np.arange(0,48,2)  # Use all time points as tick values\n",
    ")\n",
    "\n",
    "# Update y-axis titles for all subplots\n",
    "for yaxis in fig.layout:\n",
    "    if \"yaxis\" in yaxis:\n",
    "        fig.layout[yaxis].title.text = f\"<b>Neutralized CaL (%)</b>\"\n",
    "\n",
    "# Update x-axis titles for all subplots\n",
    "for xaxis in fig.layout:\n",
    "    if \"xaxis\" in xaxis:\n",
    "        fig.layout[xaxis].title.text = \"<b>Time (h)</b>\"\n",
    "fig.for_each_annotation(lambda anno: anno.update(text = \"<b>\"+anno.text+f\"; CaL {CAL_CONC} µM</b>\"))\n",
    "\n",
    "fig.update_yaxes(showticklabels=True, matches=None,range=[0,100.1])\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_AGG.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_AGG.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdc047e5-ba23-439f-843d-3d39d527538f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "data = (pd.concat(df_list).assign(Nb = lambda df_: df_.initNb-df_.Nb/(system_size*1e-3))\n",
    " .query(\"Time == 24 and Delay == 0\")\n",
    " .assign(Binding = lambda df_:(df_.NCALA/(df_.NCALA+df_.NCAL))*100)\n",
    " [[\"Binding\",\"Agg\",\"initNb\"]]\n",
    "# .pivot(index='Agg', columns='initNb')\n",
    ")\n",
    "\n",
    "fig = px.line(data.assign(Agg=lambda df_:\"<b>\"+df_.Agg.map(int).map(str)+\"</b>\"), x=\"initNb\",y=\"Binding\",color=\"Agg\", labels={\"initNb\":f\"<b>Initial {synthetic_nb} (µM)</b>\",\"Binding\":\"<b>Neutralized CaL in Aggregate (%)</b>\",\"Agg\":\"<b>Aggregate Size</b>\"})\n",
    "fig.update_yaxes(range=[50,100])\n",
    "fig.update_layout(legend=dict(orientation=\"v\",x=1.02,xanchor=\"left\",yanchor=\"top\",y=1))\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_%.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_%.png\"), format=\"png\")\n",
    "\n",
    "data = (pd.concat(df_list).assign(Nb = lambda df_: df_.initNb-df_.Nb/(system_size*1e-3))\n",
    " .query(f\"Time in {[24]} and Delay == 0\")\n",
    " .assign(Binding = lambda df_:(df_.NCALA/(df_.NCALA+df_.NCAL))*100)\n",
    " .assign(Binding= lambda df_:(1-df_.Binding/100)*df_.Agg+df_.Binding/100)\n",
    " [[\"Binding\",\"Agg\",\"initNb\",\"Time\"]]\n",
    "# .pivot(index='Agg', columns='initNb')\n",
    ")\n",
    "\n",
    "fig = px.line(data.assign(Agg=lambda df_:\"<b>\"+df_.Agg.map(int).map(str)+\"</b>\"), x=\"initNb\",y=\"Binding\",color=\"Agg\", labels={\"initNb\":f\"<b>Initial {synthetic_nb} (µM)</b>\",\"Binding\":f\"<b>Binding X Cal: 1 {synthetic_nb}</b>\",\"Agg\":\"<b>Aggregate Size</b>\"})\n",
    "fig.update_yaxes(range=[0,8])\n",
    "fig.update_layout(legend=dict(orientation=\"v\",x=1.02,xanchor=\"left\",yanchor=\"top\",y=1))\n",
    "fig.show()\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_BINDING.svg\"), format=\"svg\")\n",
    "fig.write_image(os.path.join(PLOT_FIT_FOLDER,f\"EXTRA_NCAL_BINDING.png\"), format=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2485dd4f-9527-44f6-a7e7-80c86dd9e967",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:numba]",
   "language": "python",
   "name": "conda-env-numba-py"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}