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, 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), 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 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 and make regular insulin deliveries through an insulin pump 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. 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.

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)

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.

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 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.

Now we'll use the same technique to also compute carbohydrate absorption. For this, we'll use the classic Walsh curves. We use a set of parameters for each glycaemic category: these are added together to compute the aggregate carbohydrate absorption, ca (g/5m).

We can now use Tensorflow to define and fit a model to these data. While Tensorflow is overkill here, we'll use some of its more advanced features soon.

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 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 loss function.

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.)

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.

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).

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.

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.

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.

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.

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.

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, 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.