#!/usr/bin/env python
# coding: utf-8
# # Model Predictive Control for Closed-Loop Insulin Delivery
#
# **(Or: how I learned to stop worrying and replaced my daughter's pancreas with TensorFlow)**
#
# Marius Eriksen ([@marius](http://twitter.com/marius), [marius@monkey.org](mailto:marius@monkey.org)), March 2018.
#
# People with Type-1 Diabetes do not have a functioning pancreas,
# the organ responsible for regulating insulin levels in your body.
# Insulin is a hormone that regulates the metabolism of glucose (and also fats and proteins);
# it works by signaling, among others, fat and liver cells to absorb glucose.
# Blood glucose needs to remain in a fairly narrow "Goldilocks" zone:
# high levels of blood glucose can have severe long-term health consequences;
# low levels of blood glucose can lead to unconsciousness or even death.
#
# Without a functioning pancreas,
# people living with Type-1 Diabetes must assume its job.
# Insulin must be provided when blood glucose is high
# (it is usually injected subcutaneously with a syringe or an [insulin pump](https://www.medtronicdiabetes.com/treatments/insulin-pump-therapy)),
# but not too much or blood glucose levels will go too low.
# (If glucose levels do go low, fast-acting carbohydrates must be quickly consumed
# to counteract the effects of the excess insulin.)
# Typically, people with Type-1 Diabetes deliver insulin through an insulin pump or through subcutaneous injections.
#
# What's more,
# whereas the pancreas can secrete insulin directly into the blood stream,
# so that cells are signaled with little delay,
# externally administered insulin takes much longer to reach its destination.
# [FiASP](https://www.myfiasp.com/) is the fastest acting subcutaneously delivered insulin on the market today,
# and while it enters the blood stream quite quickly,
# its maximum effect on blood glucose is between 1 and 3 hours after injection.
# This hints at yet another complication: the pharmacokinetics of insulin analogs like used to treat Type-1 Diabetes.
# After a unit of insulin is delivered, it trickles into the blood stream over the course of many hours.
# The following shows how the concentration of FiASP
# evolves in blood for the first two hours after an injection.
#
#
#
# Carbohydrates consumed in food, by the way, has its own kinetics.
# Depending on the amount of carbohydrates and the food's glycaemic index,
# a meal can take anywhere between an hour or six to fully digest:
# the carbohydrates in the meal enter the blood stream over the same course of time.
#
# And these are far from the only factors that affect blood glucose!
# Stress levels, exercise, diurnal rhythms, and illness
# all have a significant effect on how your body metabolizes sugar,
# or how your cells respond to insulin.
# So too do imperfections in insulin manufacturing,
# the site at which insulin is injected,
# the temperature at which insulin is kept,
# and even the mechanical design of an insulin pump.
#
# It's easy to see how Type-1 (T1D) diabetes is a therapy intense disorder.
# People living with T1D must constantly adjust insulin delivery (or drink some juice),
# depending on any of the above factors.
# People with T1D must also anticipate meals
# and estimate the amount of carbohydrates that they will eat.
# Ideally, they inject insulin in good time before starting a meal
# so that the insulin has time to enter the blood stream in sufficient
# quantities.
# (But don't wait too long, or your blood glucose will go too low!)
#
# Closed-loop insulin delivery systems,
# or "artificial pancreases"
# are designed to assume much of this burden.
# Such systems automate insulin delivery
# in order to maintain "Goldilocks-zone" blood glucose levels.
# These systems read glucose levels from a
# [continuous glucose monitor](https://www.dexcom.com/)
# and make regular insulin deliveries through an [insulin pump](https://www.animas.com/diabetes-insulin-pump-and-bloog-glucose-meter/onetouch-ping-blood-glucose-monitor)
# based on current and predicted glucose levels.
#
# We describe a model for predicting future blood glucose values,
# and a technique for using this model to determine how much insulin should be delivered
# to maintain a good blood glucose range.
#
# Our initial model will be quite simple,
# but the framework we set up is flexible:
# the model can be improved or even replaced later.
# A follow-up notebook will explore some additional ideas
# to improve the model.
# # A simple model for predicting blood glucose
#
# Our first model is going to be a linear approximation of
# [Bergman's minimal model](https://www.ncbi.nlm.nih.gov/pubmed/443421).
# Bergman's is a compartmental model, which presents some difficulty in
# parameter estimation ("training").
# Instead, we're going to replace Bergman's insulin and glucose compartments
# (i.e., the parts of the model that deal with the pharmacokinetics of insulin and glucose)
# with commonplace approximations
# (which are themselves models: later revisions of our model will incorporate these directly).
#
# We begin by assuming that the _change_ in blood glucose over a 5-minute period
# is linearly related to
# (1) the amount of insulin absorbed into blood;
# (2) the amount of carbohydates that are absorbed into blood;
# and (3) the amount of endogenous glucose production —
# i.e., glucose that is secreted into the blood by your liver.
#
# We can write this as:
#
# $$
# \delta_t = c_e + c_C*C_t - c_I*I_t
# $$
#
# where:
# - $C_t$: the amount of carbohydates absorbed in period $t$;
# - $I_t$: the amount of insulin absorbed in period $t$; and
# - $c_e$: the amount of endogenous production (independent of $t$).
#
# It turns out that this model is _too simple_.
# In order to compensate for this,
# we observe that blood glucose typically moves consistently:
# that is, $\delta_t$, the change in blood glucose at time $t$
# is highly dependent on the change in blood glucose at time $t-1$, $\delta_{t-1}$.
#
# So, we'll add this as another term in our model, i.e.:
#
# $$
# \delta_t = c_e + c_C*C_t - c_I*I_t + c_{\delta}*\delta_{t-1}
# $$
#
# In order to make it simpler to computing with this model,
# we sample all of our data at 5-minute intervals.
# Every 5 minutes, we have a row of data that includes:
# the amount of insulin delivered in the 5-minute period,
# the blood glucose change in the 5-minute period,
# the amount of carbohydrates delivered in the 5-minute period,
# and so on.
#
# Thus when we talk about time,
# we are talking about a single period, or row, of data.
#
# The coefficients ($c_e$, $c_C$, $c_I$, and $c_{\delta}$)
# are learned (as we'll see) by the model to weigh
# the relative importance of the various inputs.
#
# This model is too simple to be useful in practice,
# but it's a good starting point that we can build on.
# ## Building the simple model
# Let's implement this model with actual code,
# and see how it does.
#
# First, we import some Python packages that can read and clean up raw data
# that comes from our artificial pancreas system (which already samples
# data every 5 minutes).
# We also import some other utilities.
# In[1]:
from IPython.display import HTML, display
import tabulate
import numpy as np
import sys
sys.path.extend(["/Users/marius/src/tinyap.org/model"])
from tapmodel import data
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
import dateutil
tz = dateutil.tz.gettz('US/Pacific')
# We then read some recent data,
# and sample a few rows to show what's there.
#
# The columns are:
#
# - `deltaipid`: pump-delivered insulin (mU/5m)
# - `sgv`: sensor glucose value (mg/dL)
# - `ubi`: user bolus input (externally administered insulin, as we often do with FiASP for meals) (mU/5m)
# - `uciXS`, `uciS`, `uciM`, `uciL`: user carb input (entered carb values and their estimated "t-shirt sized" glycaemic index) (g/5m)
#
# In[2]:
frame = data.read_raw_csv("/Users/marius/Documents/tinyAP/frame.2018-02-10.csv", tz)
frame = frame[['sgv', 'deltaipid', 'ubi', 'uciXS', 'uciS', 'uciM', 'uciL']]
frame[np.isfinite(frame['sgv'])].tail(5)
# It's useful to look at glucose values as _deltas_:
# i.e., the change in blood glucose over a 5-minute period.
# In fact, this is what our model is going to be predicting.
# We compute these and add them as new columns to the frame.
#
# We also compute a rolling mean of glucose changes;
# this is used to smooth the effect of glucose changes.
# In[3]:
# Maxdelta is the maximum allowable glucose delta.
# We clip deltas greater than maxdelta or smaller than -maxdelta.
maxdelta = 15
frame['deltasgv'] = frame['sgv'] - frame['sgv'].shift(1)
frame['deltasgv'] = np.maximum(np.minimum(frame['deltasgv'], maxdelta), -maxdelta)
frame['meandeltasgv'] = frame['deltasgv'].rolling(window=4, min_periods=1).mean()
frame['prevmeandeltasgv'] = frame['meandeltasgv'].shift(1)
frame[np.isfinite(frame['sgv'])].tail(5)
# The model we proposed above is parameterized in terms of _active_ insulin and glucose,
# so before we proceed we have to convert glucose and insulin deliveries
# (e.g., at time $t$) into _active_ glucose and insulin values
# (at times greater than $t$).
#
# This requires us to model the pharmacodynamics of insulin.
# For the time being we'll use the model worked by [Dragan Maksimovic for Loop](https://github.com/ps2/LoopIOB/blob/master/ScalableExp.ipynb)
# with some standard parameters;
# we'll see later how we can learn its parameters from our data.
#
# Since our data contains two types of insulin deliveries
# (pump deliveries are Humalog, externally administerd insulin is FiASP),
# we'll use two different sets of parameters.
#
# Once we're done we have a new column, `ia`, in our frame
# that computes insulin activity (insulin absorption per 5 minutes)
# per our chosen pharmacodynamic model.
# In[4]:
def expia(t, tp, td):
"""Exponential insulin curve, parameterized by peak and duration,
due to Dragan Maksimovic (@dm61).
Worked by Pete (@ps2) in the notebook
https://github.com/ps2/LoopIOB/blob/master/ScalableExp.ipynb
Reworked here to be vectorized by numpy.
"""
tau = tp*(1-tp/td)/(1-2*tp/td)
a = 2*(tau/td)
S = 1/(1-a+(1+a)*np.exp(-td/tau))
return np.maximum(0., (S/np.power(tau,2))*t*(1-t/td)*np.exp(-t/tau))
# W is the window of history that we consider.
# Our periods are 5 minutes, so 12*6 = 6 hours.
W = 12*6
Wtime = np.linspace(0., 1.*W, W, endpoint=False, dtype='float32')
humalog_coeffs = expia(Wtime, 13., 48.) # 1:05h peak, 4h duration
fiasp_coeffs = expia(Wtime, 7., 45.) # 35m peak, 3:45h duration
# We flip them because we're going to be applying these to "trailing" data.
humalog_coeffs = np.flip(humalog_coeffs, 0)
fiasp_coeffs = np.flip(fiasp_coeffs, 0)
# Ia indicates insulin activity. This is computed by adding up the
# contributions of each delivery over a rolling window. Conveniently,
# this is equivalent to taking the dot product of the deliveries in the
# window with the coefficients computed above.
frame['ia'] = frame['deltaipid'].rolling(window=W).apply(lambda pids: np.dot(pids, humalog_coeffs)) + frame['ubi'].rolling(window=W).apply(lambda pids: np.dot(pids, fiasp_coeffs))
frame[np.isfinite(frame['sgv'])].tail(5)
# Now we'll use the same technique to also compute carbohydrate absorption.
# For this, we'll use the classic [Walsh curves](https://www.amazon.com/Pumping-Insulin-Everything-Need-Success/dp/1884804845).
# We use a set of parameters for each glycaemic category:
# these are added together to compute the aggregate carbohydrate absorption, `ca` (g/5m).
# In[5]:
def walshca(t, tdel, tdur):
"""Walsh carb absorption curves with provided delays and duration."""
return ((t >= tdel) & (t <= tdel+tdur/2)) * (4*(t-tdel)/np.square(tdur)) + ((t > tdel+tdur/2) & (t <= tdel+tdur)) * (4/tdur * (1-(t-tdel)/tdur))
# Carbohydrate parameters (delay, duration): XS, S, M, L.
carb_params = {
'uciXS': (2., 6.),
'uciS': (3., 12.),
'uciM': (3., 24.),
'uciL': (3., 36.),
}
carb_coeffs = {which: walshca(Wtime, tdel, tdur) for (which, (tdel, tdur)) in carb_params.items()}
carb_coeffs = {which: np.flip(coeffs, 0) for (which, coeffs) in carb_coeffs.items()}
frame['ca'] = np.zeros(len(frame))
# Add up the contributions from each carbohydrate type.
for column, coeffs in carb_coeffs.items():
frame['ca'] += frame[column].rolling(window=W).apply(lambda ucis: np.dot(ucis, coeffs))
frame[np.isfinite(frame['sgv'])].tail(5)
# We can now use [Tensorflow](https://www.tensorflow.org)
# to define and fit a model to these data.
# While Tensorflow is overkill here,
# we'll use some of its more advanced features soon.
# In[6]:
import tensorflow as tf
import os
# Quelch annoying Tensorflow warnings.
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"
class model:
# The following are the columns we're using in our model.
# The use of 'None' here indicates that the value of the dimension
# is unknown until runtime. In our case, it's the number of samples
# that we run through each iteration.
car = tf.placeholder(tf.float32, [None, 1], name='car')
ia = tf.placeholder(tf.float32, [None, 1], name='ia')
prevdelta = tf.placeholder(tf.float32, [None, 1], name='prevdelta')
# And the following are the (learned) parameters:
egp_effect = tf.Variable(0., name='egp_effect')
carb_coeff = tf.Variable(0., name='carb_coeff')
insulin_coeff = tf.Variable(0., name='insulin_coeff')
prevdelta_coeff = tf.Variable(0., name='delta_coeff')
# Our model, as defined above:
predicted_delta = egp_effect + carb_coeff*car - insulin_coeff*ia + prevdelta_coeff*prevdelta
# Now we're ready to train:
# we need to fit the coefficients (variables)
# to values that best fit our observed data.
#
# We'll do this with one of the optimizers that are bundled with Tensorflow.
# These perform [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) on a _loss function_
# in order to find the values of the variables that minimize the loss.
# We'll begin by using the standard [mean squared error](https://en.wikipedia.org/wiki/Mean_squared_error) loss function.
# In[7]:
actual_delta = tf.placeholder(tf.float32, [None, 1], name='delta')
# m is the batch size. We have to determine this dynamically
# since our input sizes may vary.
m = tf.cast(tf.shape(actual_delta)[0], tf.float32)
loss = tf.reduce_sum(tf.square(actual_delta - model.predicted_delta)) / (2*m)
# Then we train the model using a gradient descent optimizer.
# First we have to create a subset of our dataset that contain only defined values.
# (We'll later apply other filters to the data as well,
# to clean it up in various ways.)
# In[8]:
# Select all the rows that are defined for the columns we use in our model.
# Note that this can make our data frame noncontiguous, but that's ok since
# we've already computed the values (ca and ia) that depend on contiguous data.
rows = np.isfinite(frame['meandeltasgv']) & np.isfinite(frame['prevmeandeltasgv']) & np.isfinite(frame['ca']) & np.isfinite(frame['ia'])
filtered = frame[rows]
filtered.tail(5)
# We now shuffle the data and split it into a training and test set.
# The training set will be used to perform model training (parameter estimation);
# the test set is used to measure the performance of the model.
# We reserve 5% of the data (randomly) sampled for testing.
# In[9]:
shuffled = filtered.sample(frac=1)
n = int(len(shuffled)*0.95)
train, test = shuffled[:n], shuffled[n:]
# In[10]:
def train_feed_dict(batch):
"""Construct a data dictionary from a batch (a pandas frame)."""
return {
model.car: batch['ca'].values.reshape([-1, 1]),
model.ia: batch['ia'].values.reshape([-1, 1]),
model.prevdelta: batch['prevmeandeltasgv'].values.reshape([-1, 1]),
actual_delta: batch['meandeltasgv'].values.reshape([-1, 1]),
}
# The number of epochs to run through. Each epoch sees the entirety
# of the (training) data set.
nepoch = 1000
learning_step = tf.train.AdagradOptimizer(1.0).minimize(loss)
# batch_size is the number of samples we run through each mini-batch
# when training.
batch_size = 5000
nbatch = int(len(train)/batch_size)
if nbatch < 1:
raise ValueError('not enough data')
feed_dicts = []
for i in range(nbatch):
batch = train[i*batch_size:(i+1)*batch_size]
feed_dicts.append(train_feed_dict(batch))
test_feed_dict = train_feed_dict(test)
# Tensorflow needs a session within which to run our
# computation graphs.
sess = tf.InteractiveSession()
sess.run(tf.global_variables_initializer())
for e in range(nepoch):
for i in range(nbatch):
sess.run(learning_step, feed_dict=feed_dicts[i])
# Print the loss periodically so we can monitor progress.
if e%500 == 0:
est_loss = sess.run(loss, test_feed_dict)
print('epoch {} loss {}'.format(e, est_loss))
est_loss = sess.run(loss, test_feed_dict)
print('final loss: {}'.format(est_loss))
est_egp_effect, est_carb_coeff, est_insulin_coeff, est_prevdelta_coeff = sess.run([model.egp_effect, model.carb_coeff, model.insulin_coeff, model.prevdelta_coeff])
tab = [
["EGP", est_egp_effect],
["Carb", est_carb_coeff],
["Insulin", est_insulin_coeff],
["Prevdelta", est_prevdelta_coeff],
]
display(HTML(tabulate.tabulate(tab, tablefmt='html')))
# We can interpret our estimated coeffients as follows:
#
# 1. The EGP and insulin coefficients can give a _basal rate_.
# 2. The carb and insulin coefficients can give a _carb ratio_.
# 3. Prevdelta as the amount which blood glucose tends to move in the direction it's already moving.
# Since this is generally true, you can also think of this as glucose changes that probably should
# be explained by something that isn't modeled explicitly in the model.
#
# The estimated loss of $1.7$ is the mean squared error we expect from our model.
# In our case, this means that the mean square of the error in our predictions is $1.7$ every five minutes.
# ## Predicting future blood glucose values
# Now we can use our model to visualize how it fares in a real-world scenario
# by feeding the model past data (insulin and carb history for the last 6 hours),
# and then asking for its prediction.
#
# It's not terribly meaningful to infer _too_ far into the future,
# since this prediction doesn't account for _future_ carb or insulin deliveries
# (including what the AP controller will deliver; more on this soon).
# In[11]:
def plot_glucose(time, actual, predicted):
"""Plot actual vs. predicted glucose."""
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_ylabel('blood glucose (mg/dL)')
ax.set_xlabel('time')
ax.plot(time, actual, label='actual')
ax.plot(time, predicted, label='predicted')
ax.xaxis.set_major_locator(plt.MaxNLocator(10))
fmt = DateFormatter('%H:%M', tz=tz)
ax.xaxis.set_major_formatter(fmt)
plt.legend()
plt.show()
# ## Predicting a single time step
#
# The following illustrates how our model performs
# when predicting a single time step: the next glucose delta.
#
# As expected, it's generally pretty good at this.
# We can see that it follows the trends quite well,
# though it may miss some of the "spikyness" of the actual data.
# (This isn't a big problem in practice.)
#
# The first plot shows the _deltas_ predicted directly,
# the second shows these deltas accumulated over a 3 hour period
# (and thus errors are also accumulated).
# Again, this isn't necessarily a super meaningful way to visualize
# how the model is doing,
# but it gives you a good idea of cumulative error.
#
# You'll also note that the timeframe chosen was during the night time,
# which is easier for any model:
# it is a more controlled environment as
# there is (usually) little to no carbohydrate activity,
# and also little physical activity.
# In[12]:
testframe = frame[(frame.index >= '2018-02-06 00:00 PST') & (frame.index <= '2018-02-06 03:00 PST')]
predicted_deltas = sess.run(model.predicted_delta, feed_dict=train_feed_dict(testframe))
plot_glucose(testframe.index, testframe.deltasgv, predicted_deltas)
# In[13]:
predicted_sgvs = testframe.iloc[0].sgv + np.cumsum(predicted_deltas)
plot_glucose(testframe.index, testframe.sgv, predicted_sgvs)
# We can also see how the model performs if we ask if to infer 3 hours in the future.
# This is not a meaningful thing to do in this context,
# as it does not account for insulin (or carb) deliveries
# during this time period,
# but it's nevertheless a fun visualization.
# In[14]:
def feed_dict_row_delta(row, prevdelta):
"""Construct a data dictionary for inference
from a single row and a prevdelta."""
return {
model.car: np.array([row['ca']]).reshape([-1, 1]),
model.ia: np.array([row['ia']]).reshape([-1, 1]),
model.prevdelta: np.array([prevdelta]).reshape([-1, 1]),
}
predicted_deltas = np.zeros(len(testframe))
for i in range(len(testframe)):
row = testframe.iloc[i]
if i == 0:
# Start off with the initial mean delta.
prevdelta = row['prevmeandeltasgv']
else:
prevdelta = predicted_deltas[i-1]
predicted_deltas[i] = sess.run(model.predicted_delta, feed_dict=feed_dict_row_delta(row, prevdelta))
predicted_sgvs = testframe.iloc[0].sgv + np.cumsum(predicted_deltas)
plot_glucose(testframe.index, testframe.sgv, predicted_sgvs)
# ## Visualizing the model by synthesizing scenarios
# Another useful "sanity check" on our model is
# to see how it responds to certain common scenarios.
# We have certain ideas for how your body should respond
# to, say, eating carbohydrates without insulin,
# or giving a correction bolus.
#
# The following sets up code to synthesize (future) events,
# and then visualize the model's predictions for
# how the body responds to these events.
# In[15]:
from numpy.lib.stride_tricks import as_strided
def rolldot(events, coeffs, w=W):
"""Helper function to compute a rolling dot-product,
used to compute current insulin and carbohydrate activity
according to our curves."""
events = np.append(np.zeros(w), events)
stride, = events.strides
strides = as_strided(events, shape=[len(events) - w + 1, w], strides=[stride, stride])
strides = strides[1:]
return np.dot(strides, coeffs)
def feed_dict_synthetic(carbs, boluses, prevdelta):
"""Construct a data dictionary with the
provided carb and bolus events. The returned
feed dict contains only the last entry."""
return {
model.car: rolldot(carbs, carb_coeffs['uciM'])[-1].reshape([1, 1]),
model.ia: rolldot(boluses, humalog_coeffs)[-1].reshape([-1, 1]),
model.prevdelta: np.array([prevdelta]).reshape([-1, 1]),
}
def predict_glucose(carbs, boluses, npad=12*3):
"""Infers future glucose values based on the provided
carbohydrate and bolus (insulin) events."""
pad = np.zeros(npad)
carbs = np.append(carbs, pad)
boluses = np.append(boluses, pad)
prevdelta = 0
deltas = np.zeros(len(carbs))
for i in range(len(carbs)):
prevdelta = sess.run(
model.predicted_delta,
feed_dict=feed_dict_synthetic(carbs[:i+1], boluses[:i+1], prevdelta),
)
deltas[i] = prevdelta
return np.cumsum(deltas)
def plot_predict_glucose(title, events):
"""Plot the model's glucose predictions given
a set of events."""
maxtime = 0
for (time, _, _) in events:
if time > maxtime:
maxtime = time
carbs, boluses = np.zeros(maxtime+1), np.zeros(maxtime+1)
for (time, event, v) in events:
if event == 'C':
carbs[time] = v
elif event == 'B':
boluses[time] = v
else:
raise Exception("unknown event type")
glucose = predict_glucose(carbs, boluses)
n = len(glucose)
time = np.linspace(0., 5.*n, n, endpoint=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_ylabel('change in blood glucose (mg/dL)')
ax.set_xlabel('time (m)')
ax.set_title(title)
ax.plot(time, glucose)
plt.show()
plot_predict_glucose('prebolus for food', [(0, 'B', 100), (4, 'B', 2000), (6, 'C', 20)])
plot_predict_glucose('do nothing', [])
plot_predict_glucose('give big bolus', [(0, 'B', 2000)])
plot_predict_glucose('eat carbs then bolus later', [(0, 'C', 10), (4, 'B', 1000)])
# As we can see,
# our simple model follows intuition:
#
# - when we pre-bolus for food, glucose dips a little before rising
# and settling at a slightly higher level;
# - if we do nothing (i.e., no insulin delivery, no carb delivery),
# endogenous glucose production ensures that glucose levels steadily
# rise;
# - if we give a big bolus, glucose levels drop significantly;
# - if we eat carbs before bolusing, glucose levels rise and then level off.
# # Creating an insulin delivery schedule
#
# We are now ready to use our blood glucose model to create an insulin delivery controller.
# The controller should use the predictions given by the glucose model
# to learn a _schedule_ of insulin deliveries that
# maximizes the amount of time where glucose is in range (near the setpoint).
#
# Here's how it's going to work:
# given the last 6 hours of insulin and carb deliveries,
# we use our model to predict the next 6 hours.
# While the previous 6 hours of insulin delivery are fixed immutable,
# we can _train_ the schedule for the next 6 hours.
# We use a loss function that minimizes the (squared) distance from our desired setpoint.
# In[16]:
def roll_indices(n):
m = np.zeros([n, n], dtype=int)
for i in range(0, n):
m[i, :] = np.arange(i, i+n, dtype=int)
return m
# rolled is a W-by-W matrix of indices used to construct
# lookback windows for carb and insulin curves.
# These cannot be computed outside of Tensorflow anymore because they
# are part of the computation that is involved in optimization.
rolled = roll_indices(W)
class schedule:
# future_dpids is the next 6 hours of insulin deliveries, as learned
# by the schedule model.
future_dpids = tf.Variable(tf.zeros([W, 1]), name='future_dpids')
# These are placeholders used to feed the glucose model.
past_dpids = tf.placeholder(tf.float32, [W, 1], name='past_dpids')
past_ucis = tf.placeholder(tf.float32, [W, 1], name='past_ucis')
# Last_delta is used to initialize the pdgsv factor from the glucose model.
last_delta = tf.placeholder(tf.float32, [1], name='last_delta')
# Deviation is the deviation from the desired setpoint: i.e., the amount
# (up or down) glucose should move and settle to, relative to the starting
# point.
deviation = tf.placeholder(tf.float32, [], name='deviation')
dpids = tf.gather(
tf.reshape(tf.concat([past_dpids, future_dpids], 0), [-1]),
rolled,
)
ucis = tf.gather(
tf.reshape(tf.concat([past_ucis, tf.zeros([W, 1])], 0), [-1]),
rolled,
)
ia = tf.matmul(dpids, humalog_coeffs.reshape([-1, 1]))
car = tf.matmul(ucis, carb_coeffs['uciM'].reshape([-1, 1]))
egp_effect = model.egp_effect
carb_effect = car*model.carb_coeff
insulin_effect = ia*model.insulin_coeff
# Delta contains the predicted deltas, sans the prevdelta term.
# Since prevdelta uses the previous SGV delta, we must use the
# predictions yielded by the model itself to compute this term
# for future data points. We do this by scanning over the tensor
# of deltas, accumulating contributions.
delta = egp_effect + carb_effect - insulin_effect
delta = tf.scan(lambda a, x: a*model.prevdelta_coeff+x, delta, last_delta)
cumdelta = tf.scan(lambda a, x: a + x, delta)
# Define a loss function. This is the function we will optimize.
# Negative insulin deliveries are not permitted (they are not physically
# possible—it would be nice if they were!), and so we must penalize
# these heavily so that the model does not find schedules that would
# subtract insulin. (Interestingly, if it could, it could almost always
# correct positive deviations immediately, since it could give a large
# amount of insulin and then cut the tail.)
scaleneg = 100000.
loss = tf.reduce_sum(tf.pow(deviation+cumdelta, 2))/(2*W) - scaleneg*tf.minimum(tf.reduce_min(future_dpids), 0)
optimizer = tf.train.GradientDescentOptimizer(5.0)
grads, v = zip(*optimizer.compute_gradients(loss, [future_dpids]))
grads, _ = tf.clip_by_global_norm(grads, 1.)
#assert len(grads) == 1
#grads = [tf.clip_by_value(grads[0], -1., 1.)]
learning_step = optimizer.apply_gradients(zip(grads, v))
clip = future_dpids.assign(tf.maximum(future_dpids, 0.))
sess.run(future_dpids.initializer)
# With a scheduling model in hand,
# we can go through the same exercise as above:
# synthesize a scenario, and see what insulin delivery schedule is learned.
# The big difference here is that the glucose values reflect how the model
# believes glucose is going to respond not just to recent history but to the
# insulin schedule that's been learned by the model itself.
#
# We'll plot several scenarios:
# each includes both the predicted glucose schedule
# as well as the learned insulin delivery schedule.
# In[17]:
def schedule_insulin(ucis, dpids, deviation, niter=500, display_loss=False):
"""Create a 6-hour insulin schedule based on the past carb and insulin inputs,
correcting for the provided deviation."""
# Reset the future DPIDS back to their initial value.
sess.run(schedule.future_dpids.assign(np.zeros([W, 1])))
# Extend history to go W in the past.
feed_dict = {
schedule.past_ucis: ucis.reshape([-1, 1]),
schedule.past_dpids: dpids.reshape([-1, 1]),
schedule.last_delta: np.array([0.]).reshape([1]),
schedule.deviation: deviation,
}
# Learn the insulin delivery schedule.
for i in range(niter):
sess.run(schedule.learning_step, feed_dict=feed_dict)
if display_loss and i%100 == 0:
loss = sess.run(schedule.loss, feed_dict=feed_dict)
print("iter {} loss {}".format(i, loss))
dpids, cumdelta = sess.run([schedule.future_dpids, schedule.cumdelta], feed_dict=feed_dict)
return dpids, cumdelta
def plot_schedule_insulin(title, deviation, events):
"""Plot the learned insulin schedule given the provided deviation and
past events."""
ucis, dpids = np.zeros(W), np.zeros(W)
for (time, event, v) in events:
if time >= 0:
raise Exception("bad time")
if event == 'C':
ucis[time] = v
elif event == 'B':
dpids[time] = v
else:
raise Exception("unknown event type")
dpids, cumdelta = schedule_insulin(ucis, dpids, deviation)
n = len(dpids)
time = np.linspace(0., 5.*n, n, endpoint=False)
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(time, deviation+cumdelta)
ax1.set_title(title)
ax2.bar(time, dpids, width=5)
ax1.set_xlabel("minutes")
ax1.set_ylabel("deviation (mg/dL) ")
ax1.set_ylim([-40, 40])
ax2.set_ylabel("delivery (mU/5m)")
ax2.set_ylim([0, 500])
ax1.legend()
plt.show()
plot_schedule_insulin('hold steady', 0, [])
plot_schedule_insulin('correct deviation of 30', 30, [])
plot_schedule_insulin('correct deviation of -30', -30, [])
plot_schedule_insulin('hold steady after bolus', 0, [(-1, 'B', 400)])
# We can see here that the scheduling model is able to
# learn a schedule that minimizes deviation from our desired setpoint (0 on the y axis).
# Since the model minimizes the mean-squared-error,
# this means that it also takes into account cumulative effects.
# For example, in the second scenario, the learned insulin schedule
# frontloads most of the insulin (and, presumably, withholds delivery later)
# in order to quickly bring the blood glucose down to the setpoint.
# Thus, the traditional distinction between "basal" and "bolus" insulin deliveries is gone:
# the model simply delivers insulin that maximizes the amount
# of time spent near the desired setpoint.
#
# A closed-loop insulin delivery system computes a schedule at regular intervals
# in order to account for new inputs (insulin delivered, carbs consumed) and for
# unexplained deviations in the model
# (i.e., where the actual deviations differ from the expected ones).
# The system then delivers insulin according to the schedule.
# # Epilogue
#
# In this notebook,
# we have derived a (very) simple model for blood glucose,
# and shown how to make use of this model in an insulin controller,
# such as may be used by a closed-loop insulin delivery system.
#
# At a very high level,
# this is the strategy that is used in [tinyAP](https://github.com/tinyap),
# but reality demands a much more complicated model,
# to account for time-related effects, glucose-related effects,
# and other nonlinearities.
#
# The framework set up here is a flexible one:
# any modelled effect is automatically accounted for when computing insulin deliveries.
# For example, if we make EGP a function of time (as is done in the model used by tinyAP),
# time is automatically accounted for:
# if, in an hour, the estimated EGP increases or decreases,
# the learned schedule will anticipate this in good time and deliver (or withhold) insulin
# so that the desired setpoint is maintained.
# This, of course, is similar to a traditional "basal schedule",
# but with a key difference:
# Basal schedules are _delivery schedules_, that is,
# they specify the rate of insulin to be delivered at specific times,
# in anticipation of increased (or decreased) insulin needs in the future.
# With model predictive control, we instead model EGP and other effects directly;
# the insulin scheduler learns a _delivery schedule_ based instead on the directly modelled effects.
#
# Another advantage of this approach is its composability.
# Because we have a general framework for relating inputs
# (insulin, carbohydrates, time, glucose levels, etc.)
# to predicted blood glucose,
# improving the model (in the sense that it accounts for more of the observed glucose level variance)
# automatically improves the AP's performance.
# As an example, the model used by tinyAP includes parameters for both
# time (to account for different insulin sensitivies and endogenous glucose production rates throughout the day)
# and glucose levels (to account for observed effects of increased insulin resistance at higher glucose levels, presumably due to hepatic activity).
# In order to account for these observed effects,
# only the model itself was modified;
# the rest of the system remained unchanged.
#
# In the next installment,
# we'll build on the simple model here to account for some of the aforementioned effects.