so I insert the specific values the paper mentions (I'm not 100% on what "fractional slow component" is see answer figured it out).
Blue is expected, pinkred is what I get froma interpolated version of expected. Then the diamond lines are various convolutions with either or both of the two half lives active in the convolution, and obviously they don't match.
$$ y=72\%e^{-5.0t} + 28\%e^{-2.8t} $$$$ y=72\%e^{-log(2)t/5.0} + 28\%e^{-log(2)t/2.8} $$
$$ y=0.72e^{-5.0t} + 0.28e^{-2.8t} $$
$$ y=0.72e^{-log(2)t/5.0} + 0.28e^{-log(2)t/2.8} $$ Putting it together
I get the following three graphs, first is a picture of the original. Second is the ISR data extracted from the graph, and finally I have the original insulin concentration in blue and the convolution calculation in pink... and things don't look right.
final graph has,
Pink is calculated using Python (NumPy)
- original insulin content in blue
- interpolated version in red to 110 data points
- three convolutions with either both biexponential slopes or just one of them in red
- and a convolution of the original data in royal blue.
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
ax[2].plot(x, np.convolve(y, biExp(y, 0.72, 5.0, 0.28, 2.8), mode = "same"), color = "pink")
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
import scipy as sp
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
new_length = 110
new_x = np.linspace(data_ins.iloc[:,0].min(), data_ins.iloc[:,0].max(), new_length)
new_y = sp.interpolate.interp1d(data_ins.iloc[:,0], data_ins.iloc[:,1], kind='cubic')(new_x)
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(-b * xx*b)) + (c * np.exp(-d * xx*d))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = -np.log(2)/5.80 # 2.8 1/minutes
c = 0.28 # fractional slow component of 0.28?
d = 5-np.0log(2)/2.8 # 5.0 1/minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
# only graph the slow part (second)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0 , 5.0, 0.28, 2.8),
marker = '_', colorx_window = 'orange')
# only graph the fast part (first)100
ax.plot (np.linspace(-10,x_values 10),= biExp(np.linspace(-10, 10), 0.72 ,5.0, 0, 2.8),
marker = '_'x_window, color = 'green'x_window)
# graph both
ax.scatterplot(np.linspacex_values, biExp(-10x_values, 10a,b,c,d), biExp
marker = '_', color = 'blue')
# only graph the fast part (npfirst)
ax.linspaceplot(-10x_values, 10)biExp(x_values, a,b,0.72,d),
marker = 'x',5 color = 'green')
# only graph the slow part (second)
ax.0plot(x_values, 0.28biExp(x_values, 2.80,b,c,d),
marker = 'x', color = 'blue''orange')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,3424))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# turn off axis since they really don't add much
ax[0].set_axis_off()
# read the image file in a numpy array
aimg = plt.imread(f, format='jpg')
ax[0].imshow(aimg, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(new_x, new_y, color = "red")
ax[2].plot(data_ins.iloc[:, 0], data_ins.iloc[:, 1], color = "blue")
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], 0.72a,b,c,d), mode = "same")/5, # bi exponential function
color = "royalblue", marker = "d")
ax[2].0plot(new_x, np.convolve(new_y,
biExp(new_y, a,b,c,d), mode = "same") / 10, # bi exponential function
color = "maroon", marker = "d")
# break them down to components
ax[2].plot(new_x, np.convolve(new_y,
biExp(new_y, 0,b,1,d), mode = "same") / 10, # bi exponential function
color = "firebrick", marker = "d")
ax[2].28plot(new_x, 2np.8convolve(new_y,
biExp(new_y, 1,b,0,d), mode = "same") / 10, # bi exponential function
color = "pink""tomato", marker = "d")
if __name__ == "__main__":
main()
Blue is expected, pink is what I get from the convolution, and obviously they don't match.
$$ y=72\%e^{-5.0t} + 28\%e^{-2.8t} $$
$$ y=0.72e^{-5.0t} + 0.28e^{-2.8t} $$
Putting it together
I get the following three graphs, first is a picture of the original. Second is the ISR data extracted from the graph, and finally I have the original insulin concentration in blue and the convolution calculation in pink... and things don't look right.
Pink is calculated using Python (NumPy)
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
ax[2].plot(x, np.convolve(y, biExp(y, 0.72, 5.0, 0.28, 2.8), mode = "same"), color = "pink")
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = 2.8 # 2.8 minutes
c = 0.28 # fractional slow component of 0.28?
d = 5.0 # 5.0 minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
# only graph the slow part (second)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0 , 5.0, 0.28, 2.8),
marker = '_', color = 'orange')
# only graph the fast part (first)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0, 2.8),
marker = '_', color = 'green')
# graph both
ax.scatter(np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0.28, 2.8),
marker = 'x', color = 'blue')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,34))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# read the image file in a numpy array
a = plt.imread(f, format='jpg')
ax[0].imshow(a, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(data_ins.iloc[:, 0], data_ins.iloc[:, 1])
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], 0.72, 5.0, 0.28, 2.8), mode = "same"), # bi exponential function
color = "pink")
if __name__ == "__main__":
main()
Blue is expected, red is a interpolated version of expected. Then the diamond lines are various convolutions with either or both of the two half lives active in the convolution.
$$ y=72\%e^{-log(2)t/5.0} + 28\%e^{-log(2)t/2.8} $$
$$ y=0.72e^{-log(2)t/5.0} + 0.28e^{-log(2)t/2.8} $$ Putting it together
I get the following three graphs, first is a picture of the original. Second is the ISR data extracted from the graph, and the final graph has,
- original insulin content in blue
- interpolated version in red to 110 data points
- three convolutions with either both biexponential slopes or just one of them in red
- and a convolution of the original data in royal blue.
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
import scipy as sp
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
new_length = 110
new_x = np.linspace(data_ins.iloc[:,0].min(), data_ins.iloc[:,0].max(), new_length)
new_y = sp.interpolate.interp1d(data_ins.iloc[:,0], data_ins.iloc[:,1], kind='cubic')(new_x)
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(x*b)) + (c * np.exp(x*d))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = -np.log(2)/5.0 # 1/minutes
c = 0.28 # fractional slow component of 0.28?
d = -np.log(2)/2.8 # 1/minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
x_window = 100
x_values = np.linspace(-x_window, x_window)
# graph both
ax.plot(x_values, biExp(x_values, a,b,c,d),
marker = '_', color = 'blue')
# only graph the fast part (first)
ax.plot(x_values, biExp(x_values, a,b,0,d),
marker = 'x', color = 'green')
# only graph the slow part (second)
ax.plot(x_values, biExp(x_values, 0,b,c,d),
marker = 'x', color = 'orange')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,24))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# turn off axis since they really don't add much
ax[0].set_axis_off()
# read the image file in a numpy array
img = plt.imread(f, format='jpg')
ax[0].imshow(img, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(new_x, new_y, color = "red")
ax[2].plot(data_ins.iloc[:,0], data_ins.iloc[:,1], color = "blue")
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], a,b,c,d), mode = "same")/5, # bi exponential function
color = "royalblue", marker = "d")
ax[2].plot(new_x, np.convolve(new_y,
biExp(new_y, a,b,c,d), mode = "same") / 10, # bi exponential function
color = "maroon", marker = "d")
# break them down to components
ax[2].plot(new_x, np.convolve(new_y,
biExp(new_y, 0,b,1,d), mode = "same") / 10, # bi exponential function
color = "firebrick", marker = "d")
ax[2].plot(new_x, np.convolve(new_y,
biExp(new_y, 1,b,0,d), mode = "same") / 10, # bi exponential function
color = "tomato", marker = "d")
if __name__ == "__main__":
main()
and on the description of figure 2
ISR was estimated by deconvolution analysis from plasma insulin concentrations.
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = 2.8 # 2.8 minutes
c = 0.28 # fractional slow component of 0.28?
d = 5.0 # 5.0 minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
# only graph the slow part (second)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0 , 5.0, 0.28, 2.8),
marker = '_', color = 'orange')
# only graph the fast part (first)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0, 2.8),
marker = '_', color = 'green')
# graph both
ax.scatter(np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0.28, 2.8),
marker = 'x', color = 'blue')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,34))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# read the image file in a numpy array
a = plt.imread(f, format='jpg')
ax[0].imshow(a, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(data_ins.iloc[:, 0], data_ins.iloc[:, 1])
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], 0.72, 5.0, 0.28, 2.8), mode = "same"), # bi exponential function
color = "pink")
if __name__ == "__main__":
main()
CSV Export of the data
Time (M),Insulin (mU/ml),Time (M),ISR (mU/ml/min)
4.143,12.821,1.746,1.680
13.954,12.919,3.710,2.651
23.984,3.649,5.675,2.533
34.014,1.381,7.639,2.840
44.044,1.381,9.821,2.959
54.073,4.635,11.786,2.036
64.103,24.951,13.968,1.491
74.133,37.081,15.714,2.580
83.944,29.586,17.897,0.757
93.974,27.219,19.861,0.828
104.004,15.878,21.825,0.379
113.816,8.481,23.790,1.964
123.845,18.639,25.754,1.254
133.875,24.260,27.718,0.331
144.123,12.032,29.683,1.112
153.935,12.426,31.865,0.970
163.964,15.582,33.829,1.302
173.994,23.570,35.794,0.663
184.024,16.765,37.758,1.680
194.054,9.172,39.722,0.734
203.865,5.030,41.687,0.237
213.895,8.383,43.651,1.420
,,45.833,2.746
,,47.798,1.349
,,49.980,0.355
,,51.726,1.420
,,53.909,1.893
,,55.873,1.302
,,57.837,2.331
,,59.802,2.036
,,61.766,2.864
,,63.730,2.012
,,65.913,11.811
,,67.877,4.166
,,69.841,2.438
,,71.806,2.935
,,73.770,5.870
,,75.734,5.278
,,77.698,5.515
,,79.772,2.970
,,81.627,2.130
,,83.810,1.538
,,85.774,1.325
,,87.738,5.751
,,89.921,5.728
,,91.885,2.462
,,93.849,3.243
,,95.813,1.657
,,97.778,3.598
,,99.742,1.112
,,101.706,1.491
,,103.671,1.112
,,105.853,2.651
,,107.817,3.243
,,109.673,0.734
,,111.746,0.521
,,113.710,0.781
,,115.893,3.479
,,117.639,0.805
,,119.603,0.805
,,121.786,5.065
,,123.750,5.254
,,125.714,1.325
,,127.897,3.148
,,129.861,1.728
,,131.935,5.479
,,133.899,1.882
,,135.972,4.923
,,137.718,1.586
,,139.683,3.053
,,141.647,0.592
,,143.611,0.450
,,145.903,2.391
,,147.758,1.018
,,149.940,1.207
,,151.905,2.331
,,153.651,2.982
,,155.615,1.373
,,157.798,3.574
,,159.762,1.538
,,161.726,2.225
,,163.909,3.740
,,165.873,1.799
,,167.837,2.840
,,169.802,1.633
,,171.984,7.669
,,173.948,2.036
,,175.913,1.562
,,177.877,1.728
,,179.841,4.639
,,181.806,1.041
,,183.988,1.716
,,185.952,2.769
,,187.917,0.852
,,189.772,3.574
,,191.845,0.734
,,193.810,0.568
,,195.992,2.746
,,197.956,0.663
,,199.921,0.805
,,201.885,1.444
,,204.067,0.710
,,206.032,2.225
,,207.887,0.876
,,209.960,2.201
,,212.143,2.296
,,214.107,1.988
,,215.853,0.615
,,218.036,2.367
,,220.000,2.947
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = 2.8 # 2.8 minutes
c = 0.28 # fractional slow component of 0.28?
d = 5.0 # 5.0 minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
# only graph the slow part (second)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0 , 5.0, 0.28, 2.8),
marker = '_', color = 'orange')
# only graph the fast part (first)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0, 2.8),
marker = '_', color = 'green')
# graph both
ax.scatter(np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0.28, 2.8),
marker = 'x', color = 'blue')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,34))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# read the image file in a numpy array
a = plt.imread(f, format='jpg')
ax[0].imshow(a, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(data_ins.iloc[:, 0], data_ins.iloc[:, 1])
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], 0.72, 5.0, 0.28, 2.8), mode = "same"), # bi exponential function
color = "pink")
if __name__ == "__main__":
main()
and on the description of figure 2
ISR was estimated by deconvolution analysis from plasma insulin concentrations.
import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
'weight' : 'bold',
'size' : 22}
mpl.rc('font', **font)
# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only
#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')
#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')
insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)
def main():
# function of a biexponential decay
# https://swharden.com/blog/2020-09-24-python-exponential-fit/
# https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
# https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
# http://websites.umich.edu/~elements/07chap/html/07prof5.htm
#
# "This type of behavior is observed, for example, in the radioactive decay
# of a mixture of two nuclides with different half-lives, one short
# lived and the other relatively longer-lived."
# y=ae^(-bt) +ce^(-dt)
# https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
#
# A short introduction to pharmacokinetics
# R. URSO, P. BLARDI, G. GIORGI
# https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
# https://www.certara.com/knowledge-base/simplifying-deconvolution/
def biExp(x, a, b, c, d):
return (a * np.exp(-b * x)) + (c * np.exp(-d * x))
################################################################
# 2.8 and 5.0 minutes and 0.28 percent from #
# Ultradian Oscillations of Insulin Secretion in Humans #
# -------------------------------------------------------------#
# "A deconvolution method based on a biexponential #
# disappearance rate of insulin, assuming half-lives #
# for insulin of 2.8 and 5 min with a #
# fractional slow component of 28% (19)" #
# ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/ #
# https://doi.org/10.2337/diabetes.51.2007.S258 #
# -------------------------------------------------------------#
# Direct measurement of pulsatile insulin secretion from #
# the portal vein in human subjects #
# -------------------------------------------------------------#
# "a biexponential insulin disappearance model in the #
# systemic circulation, consisting of earlier directly #
# estimated half-lives of 2.8 and 5.0 min and a #
# fractional slow component of 0.28 in healthy fasting humans" #
# https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above) #
# -------------------------------------------------------------#
# In humans at least 75% of insulin secretion arises from #
# punctuated insulin secretory bursts #
# -------------------------------------------------------------#
# "insulin kinetics of 2.8 min (first half-life), #
# 5.0 min (second half-life), and #
# a fractional slow component of 0.28" #
# https://pubmed.ncbi.nlm.nih.gov/9374676/ #
# -------------------------------------------------------------#
a = 1 - 0.28 # I assume from graphpad link that speaks of percent
b = 2.8 # 2.8 minutes
c = 0.28 # fractional slow component of 0.28?
d = 5.0 # 5.0 minutes
# See if the exponential function looks correct compared to
# Concepts in Clinical Pharmacokinetics
# By Joseph T. DiPiro
# Page 78 Figure 6-9
fig, ax = plt.subplots(1, figsize=(14,10))
ax.semilogy()
ax.title.set_text("Both components in blue, and lines showing each constituent slope")
ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
ax.xaxis.set_label_text("Time (minutes)")
# only graph the slow part (second)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0 , 5.0, 0.28, 2.8),
marker = '_', color = 'orange')
# only graph the fast part (first)
ax.plot (np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0, 2.8),
marker = '_', color = 'green')
# graph both
ax.scatter(np.linspace(-10, 10), biExp(np.linspace(-10, 10), 0.72 ,5.0, 0.28, 2.8),
marker = 'x', color = 'blue')
# now calculate the actual convolution
fig, ax = plt.subplots(3, figsize=(14,34))
# create a file-like object from the url
f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")
# read the image file in a numpy array
a = plt.imread(f, format='jpg')
ax[0].imshow(a, cmap='gray')
# display ISR which was calculated orginally from the deconvolution
# of measured Insulin Concentration
ax[1].title.set_text("Insulin Secretion Rate")
ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
ax[1].xaxis.set_label_text("Time (minutes)")
ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])
# display non convolution version
# (needs some smoothing)
ax[2].title.set_text("Insulin Content")
ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
ax[2].xaxis.set_label_text("Time (minutes)")
ax[2].plot(data_ins.iloc[:, 0], data_ins.iloc[:, 1])
# convolution of ISR
ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
biExp(data_isr.iloc[:,1], 0.72, 5.0, 0.28, 2.8), mode = "same"), # bi exponential function
color = "pink")
if __name__ == "__main__":
main()
CSV Export of the data
Time (M),Insulin (mU/ml),Time (M),ISR (mU/ml/min)
4.143,12.821,1.746,1.680
13.954,12.919,3.710,2.651
23.984,3.649,5.675,2.533
34.014,1.381,7.639,2.840
44.044,1.381,9.821,2.959
54.073,4.635,11.786,2.036
64.103,24.951,13.968,1.491
74.133,37.081,15.714,2.580
83.944,29.586,17.897,0.757
93.974,27.219,19.861,0.828
104.004,15.878,21.825,0.379
113.816,8.481,23.790,1.964
123.845,18.639,25.754,1.254
133.875,24.260,27.718,0.331
144.123,12.032,29.683,1.112
153.935,12.426,31.865,0.970
163.964,15.582,33.829,1.302
173.994,23.570,35.794,0.663
184.024,16.765,37.758,1.680
194.054,9.172,39.722,0.734
203.865,5.030,41.687,0.237
213.895,8.383,43.651,1.420
,,45.833,2.746
,,47.798,1.349
,,49.980,0.355
,,51.726,1.420
,,53.909,1.893
,,55.873,1.302
,,57.837,2.331
,,59.802,2.036
,,61.766,2.864
,,63.730,2.012
,,65.913,11.811
,,67.877,4.166
,,69.841,2.438
,,71.806,2.935
,,73.770,5.870
,,75.734,5.278
,,77.698,5.515
,,79.772,2.970
,,81.627,2.130
,,83.810,1.538
,,85.774,1.325
,,87.738,5.751
,,89.921,5.728
,,91.885,2.462
,,93.849,3.243
,,95.813,1.657
,,97.778,3.598
,,99.742,1.112
,,101.706,1.491
,,103.671,1.112
,,105.853,2.651
,,107.817,3.243
,,109.673,0.734
,,111.746,0.521
,,113.710,0.781
,,115.893,3.479
,,117.639,0.805
,,119.603,0.805
,,121.786,5.065
,,123.750,5.254
,,125.714,1.325
,,127.897,3.148
,,129.861,1.728
,,131.935,5.479
,,133.899,1.882
,,135.972,4.923
,,137.718,1.586
,,139.683,3.053
,,141.647,0.592
,,143.611,0.450
,,145.903,2.391
,,147.758,1.018
,,149.940,1.207
,,151.905,2.331
,,153.651,2.982
,,155.615,1.373
,,157.798,3.574
,,159.762,1.538
,,161.726,2.225
,,163.909,3.740
,,165.873,1.799
,,167.837,2.840
,,169.802,1.633
,,171.984,7.669
,,173.948,2.036
,,175.913,1.562
,,177.877,1.728
,,179.841,4.639
,,181.806,1.041
,,183.988,1.716
,,185.952,2.769
,,187.917,0.852
,,189.772,3.574
,,191.845,0.734
,,193.810,0.568
,,195.992,2.746
,,197.956,0.663
,,199.921,0.805
,,201.885,1.444
,,204.067,0.710
,,206.032,2.225
,,207.887,0.876
,,209.960,2.201
,,212.143,2.296
,,214.107,1.988
,,215.853,0.615
,,218.036,2.367
,,220.000,2.947