https://doi.org/10.5281/zenodo.11373234
all_plots.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "208aa526-161b-4cc8-883d-1c2b801df258",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Package imports\n",
"import pickle\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.patches import Rectangle, FancyArrowPatch\n",
"import matplotlib as mpl\n",
"\n",
"mpl.rcParams['xtick.direction'] = 'in'\n",
"mpl.rcParams['ytick.direction'] = 'in'\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "53b04b38-a64f-4e1e-94ea-66425d121628",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"#SET WHERE YOU WANT THE FIGURES TO BE SAVED\n",
"output_dir = \"../../figures/\"\n",
"\n",
"\n",
"#Global plot settings\n",
"figure_width = 5.6 #inches\n",
"axis_label_size = 12 #pt\n",
"panel_label_size = 12 #pt\n",
"default_left_margin = 64 / 72\n",
"default_bottom_margin = 48 / 72\n",
"default_right_margin = plt.rcParams['axes.linewidth'] /72\n",
"default_top_margin = 4 / 72\n",
"default_inner_margin = 8 / 72 #used when subplots lack axis labels\n",
"\n",
"red = '#da291c'\n",
"blue = '#0033a0'\n",
"purple = '#702963'\n",
"gray = 'gray'\n",
"light_gray = '#D3D3D3'\n",
"\n",
"plt.rcParams['xtick.direction'] = 'in'\n",
"plt.rcParams['ytick.direction'] = 'in'\n",
"plt.rcParams['hatch.linewidth'] = 0.05\n",
"plt.rcParams[\"axes.formatter.use_mathtext\"] = True\n",
"font = {'family' : 'serif',\n",
" 'size' : 12,\n",
" 'serif': 'cmr10'\n",
" }\n",
"\n",
"mpl.rc('font', **font)\n",
"plt.rcParams['text.usetex'] = True\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "64da7eef-fda9-4a9f-b8e0-0809d994ab9b",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 1 (Schematic)\n",
"\"\"\"\n",
"\n",
"fast_colour = \"#D81B60\"\n",
"intermediate_colour = \"#1E88E5\"\n",
"slow_colour = \"#004D40\"\n",
"\n",
"\n",
"left_margin = default_right_margin\n",
"right_margin = default_right_margin\n",
"bottom_margin = default_right_margin\n",
"top_margin = default_right_margin\n",
"inner_margin = default_inner_margin\n",
"\n",
"reactant_fontsize=24\n",
"function_fontsize=18\n",
"\n",
"\n",
"h = figure_width / 3\n",
"reactant_fontsize=16\n",
"function_fontsize=12\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"panel_width = (figure_width - left_margin - inner_margin -right_margin) /3\n",
"panel_height = (h - bottom_margin - top_margin ) \n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin+panel_width+inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax3 = fig.add_axes([(left_margin+2*panel_width+2*inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"for ax in ax1, ax2, ax3:\n",
" ax.set_xlim([0, 1])\n",
" ax.set_ylim([0, 1])\n",
" ax.set_xticks([])\n",
" ax.set_yticks([])\n",
" ax.spines[['bottom','left','right', 'top']].set_visible(False)\n",
" \n",
"x1 = 1/4\n",
"x2 = 3/4\n",
"\n",
"\n",
"ax1.text(x1, 1/2, r\"$G$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"ax1.text(x2, 1/2, r\"$I$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"\n",
"style = \"Simple, tail_width=0.5, head_width=4, head_length=8\"\n",
"style = \"-|>,head_width=3,head_length=6\"\n",
"kw = dict(arrowstyle=style)\n",
"kwB = dict(arrowstyle=\"-[,widthB=4\")\n",
"kwV = dict(arrowstyle=\"Simple, tail_width=0.5\")\n",
"\n",
"r = 0.8*(reactant_fontsize/72)*(1/h)\n",
"\n",
"aGI = FancyArrowPatch((x1+r*np.cos(np.pi/6), 1/2-r*np.sin(np.pi/6)),(x2-r*np.cos(np.pi/6), 1/2-r*np.sin(np.pi/6)),\n",
" connectionstyle=\"arc3,rad=0.5\",color=fast_colour, **kw)\n",
"ax1.add_patch(aGI)\n",
"aIG = FancyArrowPatch((x2-r*np.cos(np.pi/6), 1/2+r*np.sin(np.pi/6)),(x1+r*np.cos(np.pi/6), 1/2+r*np.sin(np.pi/6)),\n",
" connectionstyle=\"arc3,rad=0.5\",color=fast_colour, **kwB)\n",
"ax1.add_patch(aIG)\n",
"\n",
"ax1.text(1/2,1/3,r\"$c \\beta f(G)$\", horizontalalignment=\"center\",verticalalignment=\"top\", fontsize=function_fontsize)\n",
"ax1.text(1/2,1-1/3,r\"$S_I$\", horizontalalignment=\"center\",verticalalignment=\"bottom\", fontsize=function_fontsize)\n",
"\n",
"\n",
"x1=1/4\n",
"x2=3/4\n",
"y1 = 1/4\n",
"y2= 5/6\n",
"\n",
"\n",
"ax2.text(x1, y1, r\"$\\beta$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"ax2.text(x2, y1, r\"$\\beta_{\\mathrm{IN}}$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"ax2.text(1/2, y2, r\"$\\beta_{\\mathrm{D}}$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"\n",
"aBN = FancyArrowPatch((x2-1.3*r*np.cos(np.pi/6), y1+1.3*r*np.sin(np.pi/6)),(x1+r*np.cos(np.pi/6), y1+r*np.sin(np.pi/6)),\n",
" connectionstyle=\"arc3,rad=0.3\", color=intermediate_colour, **kw)\n",
"ax2.add_patch(aBN)\n",
"aNB = FancyArrowPatch((x1+r*np.cos(np.pi/6), y1-r*np.sin(np.pi/6)),(x2-1.3*r*np.cos(np.pi/6), y1-1.3*r*np.sin(np.pi/6)),\n",
" connectionstyle=\"arc3,rad=0.3\", color=intermediate_colour, **kw)\n",
"ax2.add_patch(aNB)\n",
"\n",
"aBD = FancyArrowPatch((x1+r*np.cos(np.pi/3), y1+r*np.sin(np.pi/3)),(1/2-1.3*r*np.cos(np.pi/4), y2-1.3*r*np.sin(np.pi/4)),\n",
" connectionstyle=\"arc3,rad=0\", color=slow_colour, **kw)\n",
"ax2.add_patch(aBD)\n",
"\n",
"aND = FancyArrowPatch((x2-1.3*r*np.cos(np.pi/3), y1+1.3*r*np.sin(np.pi/3)),(1/2+1.3*r*np.cos(np.pi/4), y2-1.3*r*np.sin(np.pi/4)),\n",
" connectionstyle=\"arc3,rad=0\", color=slow_colour, **kw)\n",
"ax2.add_patch(aND)\n",
"\n",
"\n",
"ax2.text(1/2,y1-0.12,r\"$k_{\\mathrm{IN}} g(G)$\", horizontalalignment=\"center\",verticalalignment=\"top\", fontsize=function_fontsize)\n",
"ax2.text(1/2,y1+0.12,r\"$k_{\\mathrm{RE}}$\", horizontalalignment=\"center\",verticalalignment=\"bottom\", fontsize=function_fontsize)\n",
"\n",
"ax2.text((x1+1/2)/2-0.05,(y1+y2)/2-0.05,r\"$k_{\\mathrm{D}} h(G)$\", horizontalalignment=\"right\",verticalalignment=\"bottom\", fontsize=function_fontsize)\n",
"ax2.text((x2+1/2)/2+0.05,(y1+y2)/2-0.05,r\"$k_{\\mathrm{D}} h(G)$\", horizontalalignment=\"left\",verticalalignment=\"bottom\", fontsize=function_fontsize)\n",
"\n",
"\n",
"\n",
"ax3.text(x1, 1/2, r\"$G$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"ax3.text(x2, 1/2, r\"$c$\", verticalalignment=\"center\", horizontalalignment=\"center\", fontsize=reactant_fontsize)\n",
"ax3.text((x2+x1)/2,1/2+0.05,r\"$a(G,c)$\", horizontalalignment=\"center\",verticalalignment=\"bottom\", fontsize=function_fontsize)\n",
"\n",
"\n",
"adaptation_colour = slow_colour\n",
"aGc = FancyArrowPatch((x1+r, 1/2),(x2-r, 1/2),\n",
" connectionstyle=\"arc3,rad=0.0\", color=adaptation_colour, **kw)\n",
"ax3.add_patch(aGc)\n",
"\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.05\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax1.text(panel_label_x,panel_label_y,\"(A)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='left')\n",
"ax2.text(panel_label_x,panel_label_y,\"(B)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='left')\n",
"ax3.text(panel_label_x,panel_label_y,\"(C)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='left')\n",
"\n",
"\n",
"plt.savefig(output_dir+\"schematic_smaller.pdf\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bdc4d454-d454-480f-bc5b-fb72c3111e6d",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 2\n",
"\"\"\"\n",
"with open(\"../data/fig_2.pkl\", \"rb\") as f:\n",
" simulation_data = pickle.load(f)\n",
" \n",
" \n",
"h = 2*figure_width/3\n",
"left_margin = 0.91*default_left_margin\n",
"right_margin = 0.1\n",
"bottom_margin = 0.65*default_bottom_margin\n",
"top_margin = default_right_margin\n",
"inner_margin = default_inner_margin\n",
"\n",
"thin_inner_margin = inner_margin / 2\n",
"thin_panel_height_ratio = 0.2\n",
"\n",
"\n",
"\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin) / 2\n",
"panel_height = (h - bottom_margin - top_margin - 3*thin_inner_margin) / (2+2*thin_panel_height_ratio)\n",
"thin_panel_height = panel_height*thin_panel_height_ratio\n",
"\n",
"#debug_margins(fig)\n",
"\n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin+panel_height+3*thin_inner_margin+2*thin_panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax1a = fig.add_axes([left_margin / figure_width, (bottom_margin+panel_height+2*thin_inner_margin+thin_panel_height) / h, panel_width / figure_width, thin_panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+panel_height+3*thin_inner_margin+2*thin_panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax2a = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+panel_height+2*thin_inner_margin+thin_panel_height) / h, panel_width / figure_width, thin_panel_height / h])\n",
"\n",
"ax3 = fig.add_axes([left_margin / figure_width, (bottom_margin+thin_inner_margin+thin_panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax3a = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, thin_panel_height / h])\n",
"ax4 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+thin_inner_margin+thin_panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax4a = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, thin_panel_height / h])\n",
"\n",
"axes = dict(A=ax1, B=ax2, C=ax3, D=ax4)\n",
"beta_axes = dict(A=ax1a, B=ax2a, C=ax3a, D=ax4a)\n",
"\n",
"t = np.arange(len(simulation_data[\"A\"][\"G\"])) - simulation_data[\"A\"][\"sp\"][\"pre_sugar_days\"]\n",
"\n",
"ax1.xaxis.set_ticklabels([])\n",
"ax1a.xaxis.set_ticklabels([])\n",
"ax1a.yaxis.set_ticks([])\n",
"ax2a.xaxis.set_ticklabels([])\n",
"ax2a.yaxis.set_ticks([])\n",
"#ax3a.xaxis.set_ticklabels([])\n",
"ax3a.yaxis.set_ticks([])\n",
"#ax4a.xaxis.set_ticklabels([])\n",
"ax4a.yaxis.set_ticks([])\n",
"\n",
"ax2.xaxis.set_ticklabels([])\n",
"ax2.yaxis.set_ticklabels([])\n",
"ax4.yaxis.set_ticklabels([])\n",
"ax3.xaxis.set_ticklabels([])\n",
"ax4.xaxis.set_ticklabels([])\n",
"\n",
"\n",
"for key in axes:\n",
" data = simulation_data[key]\n",
" \n",
" axes[key].axvline(x=0, color=red, linestyle=\"--\")\n",
" axes[key].axvline(x=data[\"sp\"][\"sugar_days\"], color=red, linestyle=\"--\")\n",
" \n",
" if data[\"sp\"][\"insulin_days\"] > 0:\n",
" axes[key].axvline(x=data[\"sp\"][\"sugar_days\"]+data[\"sp\"][\"post_sugar_days\"],color=blue,linestyle='--')\n",
" axes[key].axvline(x=data[\"sp\"][\"sugar_days\"]+data[\"sp\"][\"post_sugar_days\"]+data[\"sp\"][\"insulin_days\"],color=blue,linestyle='--')\n",
"\n",
"\n",
" axes[key].plot(t, data[\"G\"], 'k')\n",
" axes[key].plot(t, 80*np.ones(len(data[\"G\"])),linestyle='dotted',color='black')\n",
" axes[key].plot(t, 130*np.ones(len(data[\"G\"])),linestyle='dotted',color='gray')\n",
" \n",
" axes[key].set_xlim([-20,data[\"sp\"][\"sugar_days\"]+data[\"sp\"][\"post_sugar_days\"]+data[\"sp\"][\"insulin_days\"]+data[\"sp\"][\"post_insulin_days\"]])\n",
" axes[key].set_ylim([-20,1200])\n",
"\n",
" if key in beta_axes:\n",
" beta_axes[key].plot(t, np.array(data['β'])/data['β'][0], color=purple)\n",
" beta_axes[key].set_ylim([-0.1,1.25])\n",
" beta_axes[key].set_xlim([-20,data[\"sp\"][\"sugar_days\"]+data[\"sp\"][\"post_sugar_days\"]+data[\"sp\"][\"insulin_days\"]+data[\"sp\"][\"post_insulin_days\"]])\n",
"\n",
"\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.95\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"ax1.text(panel_label_x,panel_label_y,\"(A)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax2.text(panel_label_x,panel_label_y,\"(B)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax3.text(panel_label_x,panel_label_y,\"(C)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax4.text(panel_label_x,panel_label_y,\"(D)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"#print(simulation_data[\"G0\"])\n",
"\n",
"ax1.set_ylabel(\"Average Plasma\\nglucose (mg/dL)\",fontsize=axis_label_size)\n",
"ax3.set_ylabel(\"Average Plasma\\nglucose (mg/dL)\",fontsize=axis_label_size)\n",
"\n",
"ax1a.set_ylabel(r\"$\\beta$\",fontsize=axis_label_size)\n",
"ax3a.set_ylabel(r\"$\\beta$\",fontsize=axis_label_size)\n",
"\n",
"\n",
"ax3a.set_xlabel(\"Time (days)\",fontsize=axis_label_size)\n",
"ax4a.set_xlabel(\"Time (days)\",fontsize=axis_label_size)\n",
"\n",
"\n",
"plt.savefig(output_dir+\"simulations.pdf\")\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "93b44f51-de49-461a-9323-94ef56be7c99",
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 3 of the revised paper.\n",
"\"\"\"\n",
"with open(\"../data/fig_M_bif.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
" \n",
"h = figure_width*0.9\n",
"left_margin = 0.95*default_left_margin\n",
"right_margin = 0.18*default_left_margin\n",
"top_margin = default_top_margin\n",
"inner_margin = default_inner_margin\n",
"bottom_margin = 0.7*default_bottom_margin\n",
"bottom_margin_2 = 0.65*default_bottom_margin\n",
"\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"#debug_margins(fig)\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin) / 2\n",
"panel_height = (h - bottom_margin - bottom_margin_2 - top_margin - inner_margin) / 3\n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin_2+bottom_margin+inner_margin+2*panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+bottom_margin_2+inner_margin+2*panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax3 = fig.add_axes([left_margin / figure_width, (bottom_margin_2 + bottom_margin+panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax4 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+bottom_margin_2+panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"ax5 = fig.add_axes([left_margin / figure_width, (bottom_margin_2) / h, panel_width / figure_width, panel_height / h])\n",
"ax6 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin_2) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"\n",
"ax1.xaxis.set_ticklabels([])\n",
"ax2.xaxis.set_ticklabels([])\n",
"\n",
"ax2.yaxis.set_ticklabels([])\n",
"ax4.yaxis.set_ticklabels([])\n",
"ax6.yaxis.set_ticklabels([])\n",
"\n",
"\n",
"\n",
"ax1.plot(data[\"A\"][\"M_list_l\"], data[\"A\"][\"G_l\"],color=blue,zorder=3)\n",
"ax2.plot(data[\"B\"][\"M_list_l\"], data[\"B\"][\"G_l\"],color=blue,zorder=3)\n",
"ax2.plot(data[\"B\"][\"M_list_h\"], data[\"B\"][\"G_h\"],color=red,zorder=3)\n",
"ax2.plot(data[\"B\"][\"M_list_u\"], data[\"B\"][\"G_u\"],color=purple, linestyle='--', zorder=2)\n",
"\n",
"\n",
"ax3.plot(data[\"A\"][\"M_list_l\"], data[\"A\"][\"beta_l\"],color=blue,zorder=3)\n",
"ax4.plot(data[\"B\"][\"M_list_l\"], data[\"B\"][\"beta_l\"],color=blue,zorder=3)\n",
"ax4.plot(data[\"B\"][\"M_list_h\"], data[\"B\"][\"beta_h\"],color=red,zorder=3)\n",
"ax4.plot(data[\"B\"][\"M_list_u\"], data[\"B\"][\"beta_u\"],color=purple, linestyle='--', zorder=2)\n",
"\n",
"\n",
"#vertical lines\n",
"M_l = data[\"B\"][\"M_list_u\"][0]\n",
"M_h = data[\"B\"][\"M_list_u\"][-1]\n",
"i = np.searchsorted(data[\"B\"][\"M_list_l\"], M_l)\n",
"j = np.searchsorted(data[\"B\"][\"M_list_h\"], M_h)\n",
"\n",
"ax2.vlines(M_l, data[\"B\"][\"G_l\"][i], data[\"B\"][\"G_u\"][0],\n",
" color='k', linestyle='--', zorder=1)\n",
"ax2.vlines(M_h, data[\"B\"][\"G_h\"][j], data[\"B\"][\"G_u\"][-1],\n",
" color='k', linestyle='--', zorder=1)\n",
"\n",
"\n",
"ax4.vlines(M_l, data[\"B\"][\"beta_l\"][i], data[\"B\"][\"beta_u\"][0],\n",
" color='k', linestyle='--', zorder=1)\n",
"ax4.vlines(M_h, data[\"B\"][\"beta_h\"][j], data[\"B\"][\"beta_u\"][-1],\n",
" color='k', linestyle='--', zorder=1)\n",
"\n",
"\n",
"ax1.set_ylim([0, 1500])\n",
"ax2.set_ylim([0, 1500])\n",
"ax3.set_ylim([-0.1, 1.25])\n",
"ax4.set_ylim([-0.1, 1.25])\n",
"ax1.set_xlim([-300, 1000])\n",
"ax2.set_xlim([-300, 1000])\n",
"ax3.set_xlim([-300, 1000])\n",
"ax4.set_xlim([-300, 1000])\n",
"ax5.set_ylim([-0.1, 1.25])\n",
"ax6.set_ylim([-0.1, 1.25])\n",
"ax5.set_xlim([-1.5, 1])\n",
"ax6.set_xlim([-1.5, 1])\n",
" \n",
"\n",
"#arrows\n",
"x_scale = (ax2.get_xlim()[1] - ax2.get_xlim()[0])/panel_width\n",
"G_y_scale =(ax2.get_ylim()[1] - ax2.get_ylim()[0]) / panel_height \n",
"beta_y_scale =(ax4.get_ylim()[1] - ax4.get_ylim()[0]) / panel_height \n",
"l = 0.075\n",
"\n",
"trans_x = -80\n",
"trans_y = 0\n",
"\n",
"arg = np.where(data[\"B\"][\"M_list_h\"] < 0 )\n",
"ax2.plot(data[\"B\"][\"M_list_h\"][arg]+trans_x,\n",
" data[\"B\"][\"G_h\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"M_list_h\"][arg][0]+trans_x\n",
"head_y = data[\"B\"][\"G_h\"][arg][0]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y-(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t1 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax2.add_patch(t1)\n",
"\n",
"trans_x = 80\n",
"trans_y = -50\n",
"arg = np.where(data[\"B\"][\"M_list_l\"] > 100 )\n",
"ax2.plot(data[\"B\"][\"M_list_l\"][arg]+trans_x,\n",
" data[\"B\"][\"G_l\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"M_list_l\"][arg][-1]+trans_x\n",
"head_y = data[\"B\"][\"G_l\"][arg][-1]+trans_y+20\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y+(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t2 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax2.add_patch(t2)\n",
"\n",
"trans_x = 50\n",
"trans_y = 10*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"M_list_l\"] > 100 )\n",
"ax4.plot(data[\"B\"][\"M_list_l\"][arg]+trans_x,\n",
" data[\"B\"][\"beta_l\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"M_list_l\"][arg][-1]+trans_x\n",
"head_y = data[\"B\"][\"beta_l\"][arg][-1]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y-(np.sqrt(3)/2)*l*beta_y_scale]])\n",
"t3 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax4.add_patch(t3)\n",
"\n",
"trans_x = -50\n",
"trans_y = -20*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"M_list_l\"] < 0 )\n",
"ax4.plot(data[\"B\"][\"M_list_h\"][arg]+trans_x,\n",
" data[\"B\"][\"beta_h\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"M_list_h\"][arg][0]+trans_x\n",
"head_y = data[\"B\"][\"beta_h\"][arg][0]+trans_y+10*beta_y_scale/x_scale\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y+(np.sqrt(3)/2)*l*beta_y_scale]])\n",
"t4 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax4.add_patch(t4)\n",
"\n",
"\n",
"#ax1.set_xticks([0, 0.5, 1.0])\n",
"#ax2.set_xticks([0, 0.5, 1.0])\n",
"\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.95\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax1.text(panel_label_x,panel_label_y,\"(A1)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax2.text(panel_label_x,panel_label_y,\"(B1)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"panel_label_x = ax3.get_xlim()[0] + (ax3.get_xlim()[1] - ax3.get_xlim()[0])*0.95\n",
"panel_label_y = ax3.get_ylim()[0] + (ax3.get_ylim()[1] - ax3.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax3.text(panel_label_x,panel_label_y,\"(A2)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax4.text(panel_label_x,panel_label_y,\"(B2)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"panel_label_x = ax5.get_xlim()[0] + (ax5.get_xlim()[1] - ax5.get_xlim()[0])*0.95\n",
"panel_label_y = ax5.get_ylim()[0] + (ax5.get_ylim()[1] - ax5.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax5.text(panel_label_x,panel_label_y,\"(A3)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax6.text(panel_label_x,panel_label_y,\"(B3)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"\n",
"\n",
"x = np.arange(-1000, 0.5)\n",
"ax1.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"ax2.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"ax3.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"ax4.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"ax5.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"ax6.fill_between(x, x+9999, x-9999, \n",
" color=light_gray,zorder=0)\n",
"\n",
"with open(\"../data/fig_F_bif.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
"\n",
"\n",
"\n",
"ax5.plot(data[\"A\"][\"F_list_l\"], data[\"A\"][\"beta_l\"],color=blue,zorder=3)\n",
"ax6.plot(data[\"B\"][\"F_list_l\"], data[\"B\"][\"beta_l\"],color=blue,zorder=3)\n",
"ax6.plot(data[\"B\"][\"F_list_h\"], data[\"B\"][\"beta_h\"],color=red,zorder=3)\n",
"ax6.plot(data[\"B\"][\"F_list_u\"], data[\"B\"][\"beta_u\"],color=purple, linestyle='--', zorder=2)\n",
"\n",
"\n",
"#vertical lines\n",
"B_l = data[\"B\"][\"F_list_u\"][-1]\n",
"B_h = data[\"B\"][\"F_list_u\"][0]\n",
"i = np.searchsorted(data[\"B\"][\"F_list_l\"], B_l)\n",
"j = np.searchsorted(data[\"B\"][\"F_list_h\"], B_h)\n",
"\n",
"\n",
"\n",
"ax6.vlines(B_l, data[\"B\"][\"beta_l\"][i], data[\"B\"][\"beta_u\"][-1],\n",
" color='k', linestyle='--', zorder=1)\n",
"ax6.vlines(B_h, data[\"B\"][\"beta_h\"][j], data[\"B\"][\"beta_u\"][0],\n",
" color='k', linestyle='--', zorder=1)\n",
"\n",
"\n",
"\n",
"x_scale = (ax6.get_xlim()[1] - ax6.get_xlim()[0])/panel_width\n",
"beta_y_scale =(ax6.get_ylim()[1] - ax6.get_ylim()[0]) / panel_height \n",
"l = 0.075\n",
"\n",
"\n",
"trans_x = -0.1\n",
"trans_y = 0.05*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"F_list_l\"] < -1)\n",
"ax6.plot(data[\"B\"][\"F_list_l\"][arg]+trans_x,\n",
" data[\"B\"][\"beta_l\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"F_list_l\"][arg][0]+trans_x\n",
"head_y = data[\"B\"][\"beta_l\"][arg][0]+trans_y\n",
"\n",
"\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y-(np.sqrt(3)/2)*l*beta_y_scale]])\n",
"t5 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax6.add_patch(t5)\n",
"\n",
"\n",
"trans_x = 0.1\n",
"trans_y = -0.05*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"F_list_h\"] > 0)\n",
"ax6.plot(data[\"B\"][\"F_list_h\"][arg]+trans_x,\n",
" data[\"B\"][\"beta_h\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"F_list_h\"][arg][-1]+trans_x\n",
"head_y = data[\"B\"][\"beta_h\"][arg][-1]+trans_y+0.01*beta_y_scale/x_scale\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y+(np.sqrt(3)/2)*l*beta_y_scale]])\n",
"t6 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax6.add_patch(t6)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"ax1.set_ylabel(\"Plasma glucose\\n(mg/dL)\",fontsize=axis_label_size)\n",
"ax3.set_ylabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"ax5.set_ylabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"\n",
"\n",
"ax1.get_yaxis().set_label_coords(-0.2,0.5)\n",
"ax3.get_yaxis().set_label_coords(-0.2,0.5)\n",
"ax5.get_yaxis().set_label_coords(-0.2,0.5)\n",
"\n",
"ax3.set_xlabel(\"Exogenous glucose intake (mg/dL$\\cdot$day)\", fontsize=axis_label_size)\n",
"ax3.get_xaxis().set_label_coords(1.05,-0.175)\n",
"\n",
"#ax4.set_xlabel(\"Exogenous glucose intake (mg / dL day)\", fontsize=axis_label_size)\n",
"\n",
"#ax1.set_xlabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"#ax2.set_xlabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"\n",
"ax5.set_xlabel(\"Exogenous insulin intake\", fontsize=axis_label_size)\n",
"ax5.get_xaxis().set_label_coords(1.05,-0.2)\n",
"\n",
"\n",
"plt.savefig(output_dir+\"sugar_insulin_bifurcations.pdf\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "537aded1-6303-40d2-883e-f8c3d112420b",
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 4 of the revised paper.\n",
"\"\"\"\n",
"\n",
"with open(\"../data/fig_gk_bif.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
"\n",
"h = figure_width/2.5\n",
"left_margin = 0.92*default_left_margin\n",
"right_margin = 0.1*left_margin\n",
"bottom_margin = 0.68*default_bottom_margin\n",
"top_margin = 0.8*default_top_margin\n",
"inner_margin = 0.75*left_margin\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"#debug_margins(fig)\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin) / 2\n",
"panel_height = (h - bottom_margin - top_margin) \n",
"#ax0 = fig.add_axes([default_right_margin/figure_width, default_right_margin/h, (figure_width - 2*default_right_margin)/figure_width, (h - 2*default_right_margin)/h]) #debugging margins\n",
"#ax0.set_xticks([])\n",
"#ax0.set_yticks([])\n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"#ax2.yaxis.set_ticklabels([])\n",
"\n",
"ax1.plot(data[\"B\"][\"gk_list_l\"], data['B']['G_l'], color=blue, zorder=3)\n",
"ax1.plot(data[\"B\"][\"gk_list_u\"], data['B']['G_u'], color=purple, linestyle='--', zorder=2)\n",
"ax1.plot(data[\"B\"][\"gk_list_h\"], data['B']['G_h'], color=red, zorder=3)\n",
"\n",
"ax2.plot(data[\"B\"][\"gk_list_l\"], data['B']['beta_l'], color=blue, zorder=3)\n",
"ax2.plot(data[\"B\"][\"gk_list_u\"], data['B']['beta_u'], color=purple, linestyle='--', zorder=2)\n",
"ax2.plot(data[\"B\"][\"gk_list_h\"], data['B']['beta_h'], color=red, zorder=3)\n",
"\n",
"\n",
"ax1.axvline(1700, linestyle='--',color=gray, zorder=0)\n",
"ax2.axvline(1700, linestyle='--',color=gray, zorder=0)\n",
"\n",
"\n",
"ax1.set_ylim([0,1050])\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.95\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax1.text(panel_label_x,panel_label_y,\"(A)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"panel_label_x = ax2.get_xlim()[0] + (ax2.get_xlim()[1] - ax2.get_xlim()[0])*0.95\n",
"panel_label_y = ax2.get_ylim()[0] + (ax2.get_ylim()[1] - ax2.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"\n",
"ax2.text(panel_label_x,panel_label_y,\"(B)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"\n",
"ax1.set_ylabel(\"Fasting plasma\\nglucose (mg/dL)\",fontsize=axis_label_size)\n",
"ax2.set_ylabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"\n",
"ax1.set_xlabel(r\"$k_{\\mathrm{IN}} /k_{\\mathrm{RE}}$\",fontsize=axis_label_size)\n",
"ax2.set_xlabel(r\"$k_{\\mathrm{IN}} /k_{\\mathrm{RE}}$\",fontsize=axis_label_size)\n",
"\n",
"\n",
"plt.savefig(output_dir+\"gk_bifurcation.pdf\")\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e8616d1f-5e87-47bd-8de5-b9f74e5cd5ee",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 5 of the revised paper\n",
"\"\"\"\n",
"\n",
"with open(\"../data/fig_B_bif.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
" \n",
"h = figure_width*2/3\n",
"left_margin = 0.95*default_left_margin\n",
"right_margin = 2*default_right_margin\n",
"bottom_margin = 0.65*default_bottom_margin\n",
"top_margin = 0.8*default_top_margin\n",
"inner_margin = default_inner_margin\n",
"bottom_margin_2 = 1.2*bottom_margin\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"#debug_margins(fig)\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin) / 2\n",
"panel_height = (h - bottom_margin - bottom_margin_2 - inner_margin) / 2\n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin+bottom_margin_2+panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin+bottom_margin_2+panel_height) / h, panel_width / figure_width, panel_height / h])\n",
"ax3 = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax4 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"#ax1.xaxis.set_ticklabels([])\n",
"#ax2.xaxis.set_ticklabels([])\n",
"\n",
"ax2.yaxis.set_ticklabels([])\n",
"ax4.yaxis.set_ticklabels([])\n",
"\n",
"#print(data[\"B\"][\"SI_list_h\"])\n",
"\n",
"ax1.plot(data[\"A\"][\"B_list_l\"], data[\"A\"][\"G_l\"],color=blue,zorder=3)\n",
"ax2.plot(data[\"B\"][\"B_list_l\"], data[\"B\"][\"G_l\"],color=blue,zorder=3)\n",
"ax2.plot(data[\"B\"][\"B_list_h\"], data[\"B\"][\"G_h\"],color=red,zorder=3)\n",
"ax2.plot(data[\"B\"][\"B_list_u\"], data[\"B\"][\"G_u\"],color=purple, linestyle='--', zorder=2)\n",
"\n",
"\n",
"#vertical lines\n",
"B_l = data[\"B\"][\"B_list_u\"][-1]\n",
"B_h = data[\"B\"][\"B_list_u\"][0]\n",
"i = np.searchsorted(data[\"B\"][\"B_list_l\"], B_l)\n",
"j = np.searchsorted(data[\"B\"][\"B_list_h\"], B_h)\n",
"\n",
"ax2.vlines(B_l, data[\"B\"][\"G_l\"][i], data[\"B\"][\"G_u\"][-1],\n",
" color='k', linestyle='--', zorder=1)\n",
"ax2.vlines(B_h, data[\"B\"][\"G_h\"][j], data[\"B\"][\"G_u\"][0],\n",
" color='k', linestyle='--', zorder=1)\n",
"\n",
"\n",
"\n",
"ax1.set_ylim([0, 1500])\n",
"ax2.set_ylim([0, 1500])\n",
"ax1.set_xlim([0, 1.6])\n",
"ax2.set_xlim([0, 1.6])\n",
"\n",
"\n",
"#arrows\n",
"x_scale = (ax2.get_xlim()[1] - ax2.get_xlim()[0])/panel_width\n",
"G_y_scale =(ax2.get_ylim()[1] - ax2.get_ylim()[0]) / panel_height \n",
"beta_y_scale =(ax4.get_ylim()[1] - ax4.get_ylim()[0]) / panel_height \n",
"l = 0.075\n",
"\n",
"trans_x = 0.05\n",
"trans_y = 0.1*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"B_list_h\"] > 1 )\n",
"ax2.plot(data[\"B\"][\"B_list_h\"][arg]+trans_x,\n",
" data[\"B\"][\"G_h\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"B_list_h\"][arg][-1]+trans_x\n",
"head_y = data[\"B\"][\"G_h\"][arg][-1]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y-(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t1 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax2.add_patch(t1)\n",
"\n",
"trans_x = -0.05\n",
"trans_y = -0.02*G_y_scale/x_scale\n",
"arg = np.where(data[\"B\"][\"B_list_l\"] <1 )\n",
"ax2.plot(data[\"B\"][\"B_list_l\"][arg]+trans_x,\n",
" data[\"B\"][\"G_l\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"B_list_l\"][arg][0]+trans_x\n",
"head_y = data[\"B\"][\"G_l\"][arg][0]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y+(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t2 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax2.add_patch(t2)\n",
"\n",
"\n",
"\n",
"\n",
"#ax1.set_xticks([0, 0.5, 1.0])\n",
"#ax2.set_xticks([0, 0.5, 1.0])\n",
"\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.95\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"ax1.text(panel_label_x,panel_label_y,\"(A1)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax2.text(panel_label_x,panel_label_y,\"(B1)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"panel_label_x = ax3.get_xlim()[0] + (ax3.get_xlim()[1] - ax3.get_xlim()[0])*0.95\n",
"panel_label_y = ax3.get_ylim()[0] + (ax3.get_ylim()[1] - ax3.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"high_panel_label_y = ax3.get_ylim()[0] + (ax3.get_ylim()[1] - ax3.get_ylim()[0])*(0.05*panel_width/panel_height)\n",
"\n",
"\n",
"x = np.arange(-1000, 0.5)\n",
"\n",
"\n",
"\n",
"ax1.set_ylabel(\"Fasting plastma\\nglucose (mg/dL)\",fontsize=axis_label_size)\n",
"\n",
"ax3.set_ylabel(\"Fasting plastma\\nglucose (mg/dL)\",fontsize=axis_label_size)\n",
"ax1.get_yaxis().set_label_coords(-0.2,0.5)\n",
"ax3.get_yaxis().set_label_coords(-0.2,0.5)\n",
"\n",
"ax1.set_xlabel(r\"$\\beta_{\\mathrm{TOT}} / \\beta_{\\mathrm{TOT},0}$\", fontsize=axis_label_size)\n",
"ax2.set_xlabel(r\"$\\beta_{\\mathrm{TOT}} / \\beta_{\\mathrm{TOT},0}$\", fontsize=axis_label_size)\n",
"\n",
"#ax3.get_xaxis().set_label_coords(1.05,-0.175)\n",
"\n",
"#ax4.set_xlabel(\"Exogenous glucose intake (mg / dL day)\", fontsize=axis_label_size)\n",
"\n",
"#ax1.set_xlabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"#ax2.set_xlabel(r\"$\\beta / \\beta_{\\mathrm{TOT}}$\",fontsize=axis_label_size)\n",
"\n",
"\"\"\"\n",
"S_I bifurcation diagram\n",
"\"\"\"\n",
"with open(\"../data/fig_SI_bif.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
"\n",
"\n",
"#print(data[\"B\"][\"SI_list_h\"])\n",
"\n",
"ax3.plot(data[\"A\"][\"SI_list_l\"], data[\"A\"][\"G_l\"],color=blue,zorder=3)\n",
"ax4.plot(data[\"B\"][\"SI_list_l\"], data[\"B\"][\"G_l\"],color=blue,zorder=3)\n",
"ax4.plot(data[\"B\"][\"SI_list_h\"], data[\"B\"][\"G_h\"],color=red,zorder=3)\n",
"ax4.plot(data[\"B\"][\"SI_list_u\"], data[\"B\"][\"G_u\"],color=purple, linestyle='--', zorder=2)\n",
"\n",
"\n",
"#vertical lines\n",
"SI_l = data[\"B\"][\"SI_list_u\"][-1]\n",
"SI_h = data[\"B\"][\"SI_list_u\"][0]\n",
"i = np.searchsorted(data[\"B\"][\"SI_list_l\"], SI_l)\n",
"j = np.searchsorted(data[\"B\"][\"SI_list_h\"], SI_h)\n",
"\n",
"ax4.vlines(SI_l, data[\"B\"][\"G_l\"][i], data[\"B\"][\"G_u\"][-1],\n",
" color='k', linestyle='--', zorder=1)\n",
"ax4.vlines(SI_h, data[\"B\"][\"G_h\"][j], data[\"B\"][\"G_u\"][0],\n",
" color='k', linestyle='--', zorder=1)\n",
"\n",
"\n",
"\n",
"ax3.set_ylim([0, 1500])\n",
"ax4.set_ylim([0, 1500])\n",
"ax3.set_xlim([0, 1.6])\n",
"ax4.set_xlim([0, 1.6])\n",
"\n",
"\n",
"#arrows\n",
"x_scale = (ax4.get_xlim()[1] - ax4.get_xlim()[0])/panel_width\n",
"G_y_scale =(ax4.get_ylim()[1] - ax4.get_ylim()[0]) / panel_height \n",
"l = 0.075\n",
"\n",
"trans_x = 0.05\n",
"trans_y = 0.1*beta_y_scale / x_scale\n",
"arg = np.where(data[\"B\"][\"SI_list_h\"] > 0.8 )\n",
"ax4.plot(data[\"B\"][\"SI_list_h\"][arg]+trans_x,\n",
" data[\"B\"][\"G_h\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"SI_list_h\"][arg][-1]+trans_x\n",
"head_y = data[\"B\"][\"G_h\"][arg][-1]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y-(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t3 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax4.add_patch(t3)\n",
"\n",
"trans_x = -0.05\n",
"trans_y = -0.02*G_y_scale/x_scale\n",
"arg = np.where(data[\"B\"][\"SI_list_l\"] <0.7 )\n",
"ax4.plot(data[\"B\"][\"SI_list_l\"][arg]+trans_x,\n",
" data[\"B\"][\"G_l\"][arg]+trans_y,\n",
" color='k',zorder=3)\n",
"head_x = data[\"B\"][\"SI_list_l\"][arg][0]+trans_x\n",
"head_y = data[\"B\"][\"G_l\"][arg][0]+trans_y\n",
"\n",
"tri_r = np.array([[head_x-l*x_scale/2, head_y], [head_x+l*x_scale/2, head_y], [head_x, head_y+(np.sqrt(3)/2)*l*G_y_scale]])\n",
"t4 = plt.Polygon(tri_r, color='k', zorder=3)\n",
"ax4.add_patch(t4)\n",
"\n",
"\n",
"panel_label_x = ax3.get_xlim()[0] + (ax3.get_xlim()[1] - ax3.get_xlim()[0])*0.95\n",
"panel_label_y = ax3.get_ylim()[0] + (ax3.get_ylim()[1] - ax3.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"high_panel_label_y = ax3.get_ylim()[0] + (ax3.get_ylim()[1] - ax3.get_ylim()[0])*(0.05*panel_width/panel_height)\n",
"\n",
"ax3.text(panel_label_x,panel_label_y,\"(A2)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"ax4.text(panel_label_x,panel_label_y,\"(B2)\", fontsize=panel_label_size, verticalalignment='top', horizontalalignment='right')\n",
"\n",
"\n",
"x = np.arange(-1000, 0.5)\n",
"\n",
"\n",
"\n",
"\n",
"ax1.get_yaxis().set_label_coords(-0.2,0.5)\n",
"ax3.get_yaxis().set_label_coords(-0.2,0.5)\n",
"\n",
"ax3.set_xlabel(\"Insulin sensitivity\", fontsize=axis_label_size)\n",
"ax4.set_xlabel(\"Insulin sensitivity\", fontsize=axis_label_size)\n",
"\n",
"plt.savefig(output_dir+\"B_SI_bifurcation.pdf\")\n",
"\n",
"#print(data[\"B\"][\"SI_list_h\"])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b3de85c8-0806-4691-b350-93977e29ef9e",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 5 of the original paper, 6 of the revised paper\n",
"\"\"\"\n",
"plt.rcParams['xtick.direction'] = 'out'\n",
"plt.rcParams['ytick.direction'] = 'out'\n",
"\n",
"h = figure_width/2.25\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"#debug_margins(fig)\n",
"with open(\"../data/fig_5.pkl\", \"rb\") as f:\n",
" data = pickle.load(f)\n",
"\n",
"cbar_width = 1.5*default_inner_margin\n",
"top_margin = 6*default_top_margin\n",
"bottom_margin = 0.72*default_bottom_margin \n",
"right_margin = 0.765*default_left_margin\n",
"left_margin = 0.64*default_bottom_margin\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin - cbar_width) / 2\n",
"panel_height = (h - bottom_margin - top_margin) \n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"cax = fig.add_axes([(left_margin + 2*panel_width + 2*inner_margin) / figure_width, bottom_margin/h, cbar_width/figure_width, panel_height /h ]) \n",
"\n",
"\n",
"\n",
"high=130\n",
"low=80\n",
"mid = (low+high)/2\n",
"\n",
"axes = [ax1, ax2, cax]\n",
"for i, key in enumerate(['A', 'B', 'C']):\n",
" ax = axes[i]\n",
"\n",
" try:\n",
" SI, beta, G, C = data[key]\n",
" beta = 3*beta/np.max(beta) #all conditions have beta ranging from 0 to 3*the default value\n",
" SI = SI/0.7\n",
" except Exception as e:\n",
" C = 0*C\n",
" beta = np.arange(75,135,0.1)\n",
" G = beta\n",
" G = np.tile(G, (len(SI),1))\n",
"\n",
"\n",
"\n",
" def hex_to_rgb(value):\n",
" value = value.lstrip('#')\n",
" lv = len(value)\n",
" return tuple(int(value[i:i + lv // 3], 16) for i in range(0, lv, lv // 3))\n",
"\n",
"\n",
"\n",
" def cmap(G, red, blue):\n",
" red_vec = np.array(hex_to_rgb(red)) / 256\n",
" blue_vec = np.array(hex_to_rgb(blue)) / 256\n",
" white_vec = np.array((1,1,1))\n",
"\n",
" if G <= low:\n",
" return blue_vec\n",
" elif G >= high:\n",
" return red_vec\n",
" elif G >= mid:\n",
" return ((G - mid)/(high - mid))*red_vec + ((high - G)/(high - mid))*white_vec\n",
" else:\n",
" return ((G - low)/(mid -low))*white_vec + ((mid - G)/(mid - low))*blue_vec\n",
"\n",
" G_0 = 80\n",
" ax.fill_between(SI, 0, np.max(beta),color=cmap(G_0, red, blue))\n",
"\n",
" for G_0 in range(81, 131):\n",
" if G_0 > 80 and G_0 < 130:\n",
" arg = np.where(G.astype(int) == G_0)\n",
" elif G_0 == 130:\n",
" arg = np.where(G.astype(int) >= G_0)\n",
" elif G_0 == 80:\n",
" arg = np.where(G.astype(int) <= G_0)\n",
"\n",
" s = []\n",
" b = []\n",
" for a in np.unique(arg[0]):\n",
" A = np.where(arg[0] == a)\n",
" B = np.argmax(arg[1][A])\n",
"\n",
" s.append(SI[arg[0][A][B]])\n",
" b.append(beta[arg[1][A][B]])\n",
"\n",
"\n",
"\n",
" #ax.plot(s,b)\n",
" if i < 2:\n",
" ax.fill_between(s, 0*np.array(b), b,color=cmap(G_0, red, blue))\n",
" else:\n",
" ax.fill_between(s, b, np.max(beta), color=cmap(G_0, red, blue))\n",
"\n",
" arg = np.where(C == 1)\n",
" s = []\n",
" b = []\n",
" for a in np.unique(arg[0]):\n",
" A = np.where(arg[0] == a)\n",
" B = np.argmax(arg[1][A])\n",
"\n",
" s.append(SI[arg[0][A][B]])\n",
" b.append(beta[arg[1][A][B]])\n",
" ax.fill_between(s, 0*np.array(b), b,alpha=0.0,hatch='////')\n",
"\n",
" arg = np.where(C == 2)\n",
" s = []\n",
" b = []\n",
" for a in np.unique(arg[0]):\n",
" A = np.where(arg[0] == a)\n",
" B = np.argmax(arg[1][A])\n",
"\n",
" s.append(SI[arg[0][A][B]])\n",
" b.append(beta[arg[1][A][B]])\n",
" ax.fill_between(s, 0*np.array(b), b,alpha=0.0,hatch='xxxx')\n",
"\n",
" ax.set_xlim([SI[0], SI[-1]])\n",
" ax.set_ylim([beta[0], beta[-1]])\n",
"\n",
"\n",
"\n",
"cax.set_xticks([])\n",
"#cax.set_ylim([low, high])\n",
"cax.set_yticks([low, mid, high])\n",
"cax.yaxis.tick_right()\n",
"cax.yaxis.set_label_position(\"right\")\n",
"\n",
"cax.set_ylabel(\"Fasting $G$ (mg/dL)\", rotation=270,verticalalignment='bottom',fontsize=axis_label_size)\n",
"\n",
"ax1.set_yticks([1, 2, 3])\n",
"ax1.set_xticks([1, 2, 3, 4])\n",
"ax1.set_ylabel(r\"$\\beta_\\mathrm{TOT} / \\beta_{\\mathrm{TOT},0}$\",fontsize=axis_label_size)\n",
"ax1.set_xlabel(r\"$S_\\mathrm{I} / S_{\\mathrm{I},0}$\",fontsize=axis_label_size)\n",
"\n",
"\n",
"ax2.set_yticks([1, 2, 3])\n",
"ax2.set_xticks([1, 2, 3, 4])\n",
"ax2.set_xlabel(r\"$S_\\mathrm{I} / S_{\\mathrm{I},0}$\",fontsize=axis_label_size)\n",
"ax2.yaxis.set_ticklabels([])\n",
"\n",
"panel_label_x = ax1.get_xlim()[0] + (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.95\n",
"panel_label_y = ax1.get_ylim()[0] + (ax1.get_ylim()[1] - ax1.get_ylim()[0])*(1-0.05*panel_width/panel_height)\n",
"\n",
"\n",
"r_w = (ax1.get_xlim()[1] - ax1.get_xlim()[0])*0.2\n",
"r_h = (ax1.get_ylim()[1] - ax1.get_ylim()[0])*0.2*panel_width/panel_height\n",
"r_x = panel_label_x - r_w/2\n",
"r_y = panel_label_y - r_h/2\n",
"\n",
"R1 = Rectangle([r_x-r_w/2, r_y-r_h/2], r_w, r_h, edgecolor='k',facecolor='white',alpha=0.75)\n",
"R2 = Rectangle([r_x-r_w/2, r_y-r_h/2], r_w, r_h, edgecolor='k',facecolor='white',alpha=0.75)\n",
"\n",
"ax1.add_patch(R1)\n",
"ax2.add_patch(R2)\n",
"\n",
"\n",
"ax1.text(r_x,r_y,\"(A)\", fontsize=panel_label_size, verticalalignment='center', horizontalalignment='center')\n",
"ax2.text(r_x,r_y,\"(B)\", fontsize=panel_label_size, verticalalignment='center', horizontalalignment='center')\n",
"\n",
"\n",
"#debug_margins(fig)\n",
"fig.savefig(output_dir+\"phases.pdf\")\n",
"\n",
"\n",
"\n",
"plt.rcParams['xtick.direction'] = 'in'\n",
"plt.rcParams['ytick.direction'] = 'in'\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f7937daf-7ec8-4f9e-976d-0cfca9bb78c6",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"Figure 7 of the revised paper.\n",
"\"\"\"\n",
"\n",
"h = figure_width/2.5\n",
"left_margin = 0.72*default_left_margin\n",
"right_margin = 0.1*left_margin\n",
"bottom_margin = 0.71*default_bottom_margin\n",
"top_margin = default_top_margin\n",
"inner_margin = 0.75*left_margin\n",
"\n",
"fig = plt.figure(figsize=(figure_width, h))\n",
"#debug_margins(fig)\n",
"\n",
"panel_width = (figure_width - left_margin - inner_margin - right_margin) / 2\n",
"panel_height = (h - bottom_margin - top_margin) \n",
"#ax0 = fig.add_axes([default_right_margin/figure_width, default_right_margin/h, (figure_width - 2*default_right_margin)/figure_width, (h - 2*default_right_margin)/h]) #debugging margins\n",
"#ax0.set_xticks([])\n",
"#ax0.set_yticks([])\n",
"\n",
"ax1 = fig.add_axes([left_margin / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"ax2 = fig.add_axes([(left_margin + panel_width + inner_margin) / figure_width, (bottom_margin) / h, panel_width / figure_width, panel_height / h])\n",
"\n",
"\n",
"import models \n",
"params_B = {'m_0': 8640, 'M': 0, 'S_E': 1.44, 'S_I': 0.72, 'I_0': 5, 'γ': 432, 'f_K': 141.4213562373095, 'f_n': 2, 'k_r': 0.041666666666666664, 'g_K': 4000, 'g_n': 2, 'g_k': 1700, 'beta_tot': 18463.965374557993, 'k_d': 0.001, 'β': 42.74066058925461}\n",
"red_vec = np.array(hex_to_rgb(red)) / 256\n",
"blue_vec = np.array(hex_to_rgb(blue)) / 256\n",
"\n",
"for G in [85, 90, 95, 100]:\n",
" color = red_vec*((G-85)/(100-85)) + blue_vec*((100 - G)/ (100-85))\n",
" t, F = models.fasting_approx_insulin_schedule(params_B, G)\n",
"\n",
" ax2.plot(t,F/params_B['γ'], color=color)\n",
"G, t = models.fasting_max_insulin_timescale(params_B)\n",
"ax1.plot(G,t, color='k')\n",
"\n",
"ax1.set_xlabel(r\"$G_{\\textrm{min}}$\")\n",
"ax2.set_xlabel(r\"Time (days)\")\n",
"\n",
"ax1.set_ylabel(\"Time to achieve\\nremission (days)\")\n",
"ax2.set_ylabel(\"Insulin dose\")\n",
"\n",
"plt.savefig(output_dir+\"protocol.pdf\")\n"
]
}
],
"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.11.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
