https://doi.org/10.5281/zenodo.11373234
all_simulations.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "37fe2bce-e5d9-4d70-afef-18901ab2d6ac",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"\"\"\"\n",
"This notebook generates all data found in the paper.\n",
"On a newish laptop, it should run in ~30 min.\n",
"\"\"\"\n",
"\n",
"%load_ext autoreload\n",
"%autoreload 2\n",
"from jax import config\n",
"config.update(\"jax_enable_x64\", True)\n",
"\n",
"import jax.numpy as jnp\n",
"import numpy as np\n",
"import pickle\n",
"\n",
"\n",
"import models\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"#This ignores a FutureWarning due to diffrax doing something that will lead to an error in future JAX releases\n",
"#You may want to comment it out to check that nothing else is going wrong.\n",
"import warnings\n",
"warnings.filterwarnings(\"ignore\")"
]
},
{
"cell_type": "raw",
"id": "aa521a7b-5c41-4dab-a509-b2bfc66f9837",
"metadata": {},
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "e18ad9e0-2d19-4632-925a-8cf025e47015",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def fn(n):\n",
" return \"../data/fig_\"+str(n)+\".pkl\"\n",
"\n",
"\n",
"g = 432\n",
"I0 = 5\n",
"S_E = 1.44\n",
"S_I = 0.72\n",
"m0 = 864*10\n",
"\n",
"common_params = dict(m_0 = m0, M=0, S_E=S_E, S_I=S_I, I_0=I0, γ=g, f_K=np.sqrt(20000), f_n=2)\n",
"uncommon_params = {}\n",
"uncommon_params[\"A\"] = dict(k_r=1/24, g_K=4000, g_n=2, g_k=116)\n",
"uncommon_params[\"B\"] = dict(k_r=1/24, g_K=4000, g_n=2, g_k=1700)\n",
"uncommon_params[\"C\"] = dict(k_r=1/36, g_K=4000, g_n=2, g_k=1700*1.5)\n",
"\n",
"\n",
"#set beta_tot to give fasting G of 110 before the emergency\n",
"for key in [\"A\", \"B\", \"C\"]:\n",
" params = dict(G_0=110.0, **common_params, **uncommon_params[key] )\n",
" B0 = g*models.adaptive_beta(params)\n",
" uncommon_params[key][\"beta_tot\"] = B0\n",
" \n",
"\n",
"small_meal_rate = 300.0\n",
"large_meal_rate = 8500.0\n",
"insulin_rate = 70.0*g/100 \n",
"\n",
"def fig_2_data(params_A, params_B, params_C, params_D, sp_A, sp_B, sp_C, sp_D):\n",
" g_A, _, beta_A, _, _, _ = models.full_simulation(params_A,**sp_A)\n",
" g_B, _, beta_B, _, _, _ = models.full_simulation(params_B, **sp_B)\n",
" g_C, _, beta_C, _, _, _ = models.full_simulation(params_C, **sp_C)\n",
" g_D, _, beta_D, _, _, t = models.full_simulation(params_D, **sp_D)\n",
" \n",
" with open(fn(2), \"wb\") as f:\n",
" pickle.dump(dict(t=t, A=dict(G=g_A, β=beta_A, sp=sp_A), B=dict(G=g_B, β=beta_B, sp=sp_B), C=dict(G=g_C, β=beta_C, sp=sp_C), D=dict(G=g_D, β=beta_D, sp=sp_D)), f)\n",
" \n",
"def fig_5_data(params_A, params_B):\n",
" B0_B = params_B['beta_tot']\n",
"\n",
" data = {}\n",
" data['A'] = models.SI_beta_scan(params_A, SI_range=[0.1, 3.0], beta_range=[0, 3*B0_B], SI_samples=100,beta_samples=100, G_0=80,a=0.0)\n",
" data['B'] = models.SI_beta_scan(params_B, SI_range=[0.1, 3.0], beta_range=[0, 3*B0_B], SI_samples=100,beta_samples=100, G_0=80,a=0.0,table_beta_samples=6000)\n",
"\n",
" with open(fn(5), \"wb\") as f:\n",
" pickle.dump(data, f)\n",
"\n",
"\n",
"def fig_6_data(params_B, sp_B, G_min):\n",
" g_0, I_0, beta_0, _, F_0, t = models.full_simulation(params_B, **sp_B)\n",
"\n",
"\n",
" sp_prime = sp_B | {\"insulin\": models.fasting_maximum_insulin(G_min)}\n",
" g_1, I_1, beta_1, _, F_1, _ = models.full_simulation(params_B, **sp_prime)\n",
"\n",
" with open(fn(6), \"wb\") as f:\n",
" pickle.dump(dict(t=t, A=dict(G=g_0, I = I_0, β=beta_0, F=F_0, sp=sp_B), B=dict(G=g_1, I = I_1, β=beta_1, F=F_1, sp=sp_B)), f)\n",
"\n",
"\n",
"def compute_M(params):\n",
" def computer(M):\n",
" G, f = models.all_eq(params | {\"M\":M })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
" return computer\n",
"\n",
"def compute_SI(params):\n",
" def computer(SI):\n",
" G, f = models.all_eq(params | {\"S_I\":SI })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
" return computer\n",
"\n",
"def compute_B(params):\n",
" def computer(B):\n",
" G, f = models.all_eq(params | {\"β\":params['β']*B })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
" return computer\n",
"\n",
"def compute_F(params):\n",
" def computer(F):\n",
" G, f = models.all_eq(params | {\"F_I\": F})\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
" return computer\n",
"\n",
"def compute_gk(params):\n",
" def computer(k):\n",
" G, f = models.all_eq(params | {\"g_k\": k})\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
" return computer\n",
"\n",
"def fixed_step_bisect(compute_func, bisect_func, x0, max_x, delta):\n",
" X = []\n",
" Y = []\n",
" x = x0\n",
" prev_y = compute_func(x)\n",
" X.append(x)\n",
" Y.append(prev_y)\n",
" x = x0+delta\n",
" y = compute_func(x)\n",
" while (not bisect_func(y, prev_y)):\n",
" X.append(x)\n",
" Y.append(y)\n",
" prev_y = y\n",
" x = x+delta\n",
" y = compute_func(x)\n",
" if x >= max_x:\n",
" return X, Y\n",
"\n",
" return X, Y\n",
"\n",
"def ramp_up_step(compute_func, x0, delta, depth):\n",
" X = []\n",
" Y = []\n",
" x = x0\n",
" prev_y = compute_func(x)\n",
" X.append(x)\n",
" Y.append(prev_y)\n",
" for i in range(depth):\n",
" x = x0+ delta*(2**i)\n",
" y = compute_func(x)\n",
" X.append(x)\n",
" Y.append(y)\n",
" return X,Y\n",
" \n",
"\n",
"def bisecting_sweep(compute_func, bisect_func, x0, delta, max_x, max_depth, snap=True, symmetric=True):\n",
"\n",
" assert snap\n",
" X = []\n",
" Y = []\n",
" x = x0\n",
" depth = 0\n",
" prev_y = compute_func(x)\n",
" X.append(x)\n",
" Y.append(prev_y)\n",
" x = x0+delta\n",
" \n",
" while True:\n",
" temp_X, temp_Y = fixed_step_bisect(compute_func, bisect_func, x, max_x, delta / 2**depth)\n",
"\n",
" X = X + temp_X[1:]\n",
" Y = Y + temp_Y[1:]\n",
"\n",
" if X[-1] >= max_x:\n",
" return X, Y\n",
"\n",
" elif depth < max_depth:\n",
" depth += 1\n",
" x = X[-1]\n",
"\n",
" else: #snap to grid over the bisection point \n",
" if symmetric:\n",
" temp_X, temp_Y = ramp_up_step(compute_func, x, delta / 2**depth, depth)\n",
" X = X + temp_X[1:]\n",
" Y = Y + temp_Y[1:]\n",
" depth = 0\n",
" x = delta*(X[-1]//delta) + delta\n",
"\n",
" return X, Y\n",
"\n",
"def compute_B_M(M):\n",
" G, f = models.all_eq(params_B | {\"M\":M })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
"def compute_B_SI(S_I):\n",
" G, f = models.all_eq(params_B | {\"S_I\":S_I })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
"def compute_B_B(B):\n",
" G, f = models.all_eq(params_B | {\"β\":params_B['β']*B })\n",
" if len(G) > 1:\n",
" return np.array([G[0], G[1], G[2], f[0], f[1], f[2]])\n",
" else:\n",
" return np.array([G[0], -116, -116, f[0], -116, -116])\n",
"\n",
"\n",
"\n",
"def bisect(prev_y, y):\n",
" return not np.sum(y > 0) == np.sum(prev_y > 0)\n",
"\n",
"\n",
"def total_flux(params, G, alt_mode=False):\n",
" f = models.eq_beta_frac(params, G)\n",
" I = f*models.evaluate(models.get_f(params), G)\n",
"\n",
" if alt_mode:\n",
" return (params[\"S_E\"] + params[\"S_I\"]*I)*G\n",
" else:\n",
" return params[\"M\"] + params[\"m_0\"]/(I+params[\"I_0\"])\n",
"\n",
"def process(params, M_list, G_list, β_list, flux=False):\n",
" arg = np.where(G_list > 0)[0]\n",
" m = np.zeros(len(arg))\n",
"\n",
" for i, j in enumerate(arg):\n",
" if flux:\n",
" m[i] = total_flux(params | {\"M\": M_list[j]}, G_list[j])\n",
" else:\n",
" m[i] = M_list[j]\n",
" return m, G_list[arg], β_list[arg]\n",
"\n",
"\n",
"def bifurcation_data(params, compute_func, min_param, param_spacing, max_param, param_name, key):\n",
" input_list, blah = bisecting_sweep(compute_func(params), bisect, min_param, param_spacing, max_param, 32, snap=True)\n",
" input_list = np.array(input_list)\n",
" G_l = np.transpose(blah)[0]\n",
" G_u = np.transpose(blah)[1]\n",
" G_h = np.transpose(blah)[2]\n",
" beta_l = np.transpose(blah)[3]\n",
" beta_u = np.transpose(blah)[4]\n",
" beta_h = np.transpose(blah)[5]\n",
" if not key == \"A\":\n",
" arg = np.where(G_l>np.max(G_h))\n",
" G_h[arg] = G_l[arg]\n",
" G_l[arg] = -116\n",
" beta_h[arg] = beta_l[arg]\n",
" beta_l[arg] = -116\n",
"\n",
" \n",
" \n",
" cut_l = np.where(G_l >= 80)\n",
" cut_u = np.where(G_u >= 80)\n",
" cut_h = np.where(G_h >= 80)\n",
"\n",
" return { param_name+\"_list_l\": input_list[cut_l], \"G_l\": G_l[cut_l], \"beta_l\": beta_l[cut_l],\n",
" param_name+\"_list_u\": input_list[cut_u], \"G_u\": G_u[cut_u], \"beta_u\": beta_u[cut_u], \n",
" param_name+\"_list_h\": input_list[cut_h], \"G_h\": G_h[cut_h], \"beta_h\": beta_h[cut_h] }\n",
" \n",
"\n",
"def fig_M_bif_data(params_A, params_B):\n",
"\n",
" output = {}\n",
" for key, params in zip([\"A\", \"B\"], [params_A, params_B]):\n",
" output[key] = bifurcation_data(params, compute_M, -1000, 1, 1000, \"M\", key) \n",
"\n",
" with open(fn(\"M_bif\"), \"wb\") as f:\n",
" pickle.dump(output, f)\n",
"\n",
"\n",
"def fig_SI_bif_data(params_A, params_B):\n",
"\n",
" output = {}\n",
" for key, params in zip([\"A\", \"B\"], [params_A, params_B]):\n",
" output[key] = bifurcation_data(params, compute_SI, 0.01, 0.01, 2.4, \"SI\", key) \n",
"\n",
" with open(fn(\"SI_bif\"), \"wb\") as f:\n",
" pickle.dump(output, f)\n",
" \n",
"def fig_B_bif_data(params_A, params_B):\n",
"\n",
" output = {}\n",
" for key, params in zip([\"A\", \"B\"], [params_A, params_B]):\n",
" output[key] = bifurcation_data(params, compute_B, 0.1, 0.01, 2.0, \"B\", key) \n",
" \n",
"\n",
" with open(fn(\"B_bif\"), \"wb\") as f:\n",
" pickle.dump(output, f)\n",
" \n",
"def fig_F_bif_data(params_A, params_B):\n",
" output = {}\n",
" for key, params in zip([\"A\", \"B\"], [params_A, params_B]):\n",
" output[key] = bifurcation_data(params, compute_F, -1.5, 0.01, 1.5, \"F\", key) \n",
"\n",
" with open(fn(\"F_bif\"), \"wb\") as f:\n",
" pickle.dump(output, f)\n",
" \n",
" \n",
"def fig_gk_bif_data(params_B):\n",
" output = {}\n",
" for key, params in zip([\"B\"], [params_B]):\n",
" output[key] = bifurcation_data(params, compute_gk, 0.1, 1, 3000, \"gk\", key) \n",
" \n",
"\n",
" with open(fn(\"gk_bif\"), \"wb\") as f:\n",
" pickle.dump(output, f)\n",
" \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b5d109ea-dfcc-4cc5-b692-2e6a102923a6",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"%%time \n",
"params_A = dict(**common_params, **uncommon_params[\"A\"], k_d=1e-3) \n",
"params_A['β'] = params_A['beta_tot'] / params_A['γ']\n",
"params_B = dict(**common_params, **uncommon_params[\"B\"], k_d=1e-3)\n",
"params_B['β'] = params_B['beta_tot'] / params_B['γ']\n",
"params_C = dict(**common_params, **uncommon_params[\"C\"], k_d=1e-3)\n",
"params_C['β'] = params_C['beta_tot'] / params_C['γ']\n",
"params_D = dict(**common_params, **uncommon_params[\"C\"], k_d=7*1e-2)\n",
"params_D['β'] = params_D['beta_tot'] / params_D['γ']\n",
"\n",
"sp_default = dict(pre_sugar_days=14, sugar_days=28, post_sugar_days=7, insulin_days=28, post_insulin_days=140,low_sugar=small_meal_rate, high_sugar=large_meal_rate, insulin=insulin_rate)\n",
"sim_params = {}\n",
"for key, extra_insulin_days in zip([\"A\", \"B\",\"C\",\"D\"], [-28, 0, 28, 56]):\n",
" sim_params[key] = dict(sp_default,insulin_days=sp_default[\"insulin_days\"]+extra_insulin_days, post_insulin_days=sp_default[\"post_insulin_days\"]-extra_insulin_days)\n",
"\n",
"fig_2_data(params_A, params_B, params_C, params_D, sim_params[\"A\"], sim_params[\"B\"], sim_params[\"C\"], sim_params[\"D\"])\n",
"\n",
"output = fig_3_data(params_A, params_B)\n",
"fig_4_data(params_A, params_B, insulin_rate)\n",
"blah = fig_5_data(params_A, params_B)\n",
"fig_6_data(params_B, sim_params[\"B\"], 90)\n",
"fig_M_bif_data(params_A, params_B)\n",
"fig_SI_bif_data(params_A, params_B)\n",
"fig_B_bif_data(params_A, params_B)\n",
"fig_F_bif_data(params_A, params_B)\n",
"fig_gk_bif_data(params_B)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35b4070d-f987-4976-80ad-4800491a31ef",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"#check rate of β-cell decline at 150 mg/dL glucose, as mentioned in discussion\n",
"models.evaluate(models.hill(params_B[\"g_K\"], params_B[\"g_n\"], params_B[\"g_k\"]*params_B[\"k_r\"]), 150.0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3f0e4d6c-19be-4d5e-b2c6-c4b0636dc3fa",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"#estimate grams of sugar produced per day if rate = 2 mg/(kg min), then estimate corresponding rate in units of mg/(dL day)\n",
"#for discussion in methods section\n",
"rate = 2\n",
"body_mass = 70\n",
"minutes_in_day = 24*60\n",
"grams_per_mg = 0.001\n",
"print(\"rough amount of sugar produced per day (grams):\",rate*body_mass*minutes_in_day*grams_per_mg)\n",
"\n",
"body_plasma = 35 #3.5 L, in dL\n",
"print(\"rough HGP rate\", rate*body_mass*minutes_in_day/body_plasma)\n",
"\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
}
