{ "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 }