Linear polynomial representations of Insulin on Board for Apidra, Novolog, and Humalog

[Monica C.]

Team:

I have created linear polynomial bisquare approximations (9th order--best fit) of insulin on board (IOB) for Apidra so far, for durations of insulin action (DIA) ranging from 3-8 hours, in 15 minute time intervals. It is based off of Clinical Pharmacokinetics and Pharmacodynamics of Insulin Glulisine. Per Matthias' suggestion, I used papers modeling euglycemic clamp methods, which is truly the only way to reliably model insulin pharmacokinetics, versus dosing off of just weight, as detailed here:

Principle of the Euglycaemic Clamp Method

Characterization of the concentration-time and action-time profiles of any insulin product is fundamental in understanding the clinical implication for therapy. An estabilished method of simultaneously assessing pharmacokinetic charachteristics and quantifying pharmacodynamic characteristics of insulin products is the euglycemic clamp technique. In this experimental paradigm, the decrease in the blood glucose level in response to the injected insulin is prevented by an intravenous glucose infusion, which can manually or automatically adjusted (by a Biostator), while conventional pharmacokinetic variables are derived from concentration profiles obtained from samples taken at defined times after injection. The rate and duration of the glucose infusion required to maintain the blood glucose concentration at or slightly below the fasting level reflects the metabolic activity of the injected insulin.

Of course, there is an equation for each 15 minute time interval, and you have to pick which DIA is an appropriate "fit". I did basic calculations in Excel, and did more advanced calculations in MATLAB. I copied and pasted the equations from MATLAB (and verified they were correct) into Excel, along with the statistical information and a graph of the curve fit.

More will be coming soon, and I anticipate it will be done by morning.

[Monica C.]

I forgot to say that this Apidra was modeled at 0.15 U/kg roughly. Humalog and Apidra can be modeled at ranges approximately from 0.1-0.4 U/kg. Novolog is modeled at only 0.2 U/kg, which is a relatively moderate dose. I might add more to these documents if time persists.

[Monica C.]

One more thing: Time is the input, and with the time elapsed you get a percentage value representing insulin on board as output for the equation you use. To get the fractional percentage value, you divide by the output of the equation by 100. To get the IOB remaining in units, you multiply the fractional percentage value, representing remaining IOB, times the initial insulin dosage. For "stacked boluses" or aggregate IOB (basal + bolus) you would ideally use loops, typically a for or a while loop to deal with that.

On Sunday, February 22, 2015 at 2:43:33 AM UTC-6, Monica C. wrote:

[Monica C.]

Same thing, just with Novolog.

The data comes from a euglycemic clamp study: Insulin aspart has a shorter duration of action than human insulin over a wide dose-range

I am going to take a break for an hour. This modeling is actually really boring and not engaging for me. I anticipate being done with all of it, aka the final Humalog part, at 10 am Central Time.

[Monica C.]

Actually, Novolog and Humalog are completely interchangable and are virtually identical, according to this document of a euglycemic clamp trial (see the first figure on page 4). If anyone has objections to me not modeling Humalog, please let me know.

It would honestly, in my opinion, be a waste of time to model both Humalog and Novolog. Apidra, however, is different.

On Sunday, February 22, 2015 at 2:43:33 AM UTC-6, Monica C. wrote:

[Monica C.]

I didn't sleep all night and I forgot to do the last 4 excel sheets for Humalog and Novolog. I haven't checked my work, but theoretically it should all work. I am going to rest a little bit and work on it some more in a couple of hours,to confirm it's all correct. For the record, the corrected version Humalog/Novolog can be found here.

[Bert Berla]

Hi Monica,

Is there a reason for fitting a 9th order polynomial to these curves as opposed to modeling with a simple decay? How much better of a fit are you getting?

It seems like this would require individuals to derive a huge number of parameters to personalize the equations, and it couldn't be done easily. Or am I misunderstanding?

[Matthias Granberry]

I’ll put a tl;dr at top. A 9th-order polynomial does’t require more input from a user, but it might not be appropriate when used directly in a system.

A 9th order fit doesn’t require extra parameters because it is in terms of higher powers of a single variable, but it might be susceptible to significant overfitting. A much smoother and likely more appropriate fit for most data would be a second-order polynomial or possibly a higher-order fit paired with a regularization term to minimize the “wiggle” that extreme values in high order terms will give to the function.

I suspect that minimizing the “jerk” of the IOB graph will result in a more stable system than having really accurate predictions for the IOB estimate once it is coupled with a nonlinear controller. That balance often gets struck with quadratic or possibly second-order fittings. Another concern is that the IOB from a single bolus should be monotonically decreasing, which becomes more tedious to guarantee with a high-order fitting.

A B-spline that ensures continuity along the function allows some extra helpful properties (continuous, differentiable even at knots) and they can be constructed to guarantee monotonicity and still obtain a very accurate fit. The graphs I have seen look like a 2- or 3-segment quadratic B-spline would probably work just fine. If your control system doesn’t require the curve to be continuous and differentiable, a quadratic spline with a knot at the point of transition from increasing to decreasing activity will provide a much tamer function than a higher-order polynomial.

Having said this, it should provide a good match for the data as long as it isn’t overfitted. A simple graph of the function over time will reveal this. If it looks like a low-frequency pulse added to a ringing bell on an oscilloscope, then it is overfitted. If it works well, it should be simple enough to fit a new curve with necessary properties to the high-order curves.

> --
> You received this message because you are subscribed to the Google Groups "OpenAPS Dev" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to openaps-dev...@googlegroups.com.
> To post to this group, send email to opena...@googlegroups.com.
> Visit this group at http://groups.google.com/group/openaps-dev.
> To view this discussion on the web visit https://groups.google.com/d/msgid/openaps-dev/138d7dfc-5604-453a-a3bb-8260d8388779%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Monica C.]

Team:

Disregard my above work. Insulin on board is a really funky function that cannot be approximated well off of graphs, and all of the work above has to be disregarded/discarded.

In lieu of the past technique, I did a 5 hour experiment using 4 insulin pumps (Tandem t:slim, Animas 2020, Medtronic 530G, Medtronic 522) where I recorded the IOB data every 5 minutes and extrapolated the data that way. Using the 61 data points from each respective pump, I "stretched" or "contracted" the time intervals, at equal lengths, to fit in time intervals from 3-8 hours, in 30 minute increments, while keeping the 61 data points for IOB and the whole number percentage IOB remaining.

So far, I have the equation for the whole number percentage IOB remaining for the Tandem t:slim, which I attached. Due to file size, I will have to make multiple posts to attach the remainder of my work, so I apologize for the annoyance.

Note that the equations are exactly the same for each of the DIAs. This is because the value of x is normalized as the variable, that is multiplied by the constant p*, for whatever corresponding coefficient in the 5th degree linear autoregressive polynomial approximation. I had to use the "center and fit" feature, in order for the 5th degree linear autoregressive polynomial approximation to fit and converge, which required normalization of the variable time (time elapsed from initial bolus--or series of boluses or basal deliveries--if using programming loops).

So, to get the nomalized variable:

x = (time in minutes - mean)/std

where x is to be substituted in the 5th degree linear autoregressive polynomial approximation. time is the time elapsed from an initial bolus (or series of boluses or basal deliveries--if using programming loops) where the whole number percentage value of IOB remaining is to be found, mean is the mean value provided for the given DIA used for the 5th degree linear autoregressive polynomial approximation, and std is the standard deviation provided for the given DIA for the 5th degree linear autoregressive polynomial approximation.

The equation

Linear model Poly5:
     f(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6

returns the whole number percentage value of IOB remaining. So to determine the actual IOB remaining, you would divide the return value of the equation, f(x) by 100, and multiply f(x)/100 by the original bolus amount. If using series of boluses or basal deliveries, you would do this iteratively using programming loops and subsequently sum up the resulting IOB remaining.

Please note that I have not verified my results and "tested it out", but I have no reason to believe it is wrong. I will check my work when I complete all of the modeling for the remaining 3 pumps.

[Monica C.]

So this is the final product. I confirmed they were all correct. The order of the polynomials range from the 4th-7th degree, which depends on the equation converging in MATLAB. Linear and least squares autoregression was used to fit the polynomials to the percent insulin on board (IOB) remaining equation. To calculate IOB remaining in units for a single bolus given at a specified time, input the time elapsed in minutes from the initial bolus and find the solution percent IOB remaining.  Divide percent IOB remaining by 100 to get the fractional IOB remaining and multiply the solution by the initial bolus dose in units. This will give you IOB remaining in units.

Note that the polynomial equations are exactly the same for each individual pump, irrespective of the duration of insulin action (DIA) chosen. This is because the variable x (which is not t, or time elapsed from initial bolus), is normalized for each DIA time. The normalization values differ for each DIA. I had to normalize the equations in order for the fit to converge in MATLAB. To normalize and find x:

x = (t - mean)/std

Where t is time elapsed (typically from the initial bolus dosage, t = 0), mean is the mean value of t, and std is the standard deviation of t. x is the normalized value of t, which is plugged into a typical polynomial equation (not an actual approximation) such as:

f(x) = p1*(x^2) + p2*(x) + p3

I verified my work was all correct. If you have any questions, please ask.

[Monica C.]

Medtronic 522 pump

[Monica C.]

Animas 2020 pump

[Bert Berla]

On Monday, February 23, 2015 at 6:03:42 PM UTC-5, Monica C. wrote:
> Animas 2020 pump

My other concern is that calculating IOB from clamp studies but making dose decisions based on CGM data is a bad match. If we assume that interstitial glucose lags behind veinous, then this approach will lead to underestimation of the effect of the elapsed part of a dose and insulin stacking. I guess this can be prevented by just requiring a fixed spacing between doses but can you give a more basic explanation of why you're taking the current approach to estimating IOB as opposed to something far simpler like exponential decay plus a lag, for example? Or something that might have more of a mechanistic underpinning?

[Monica C.]

I actually ditched the clamp studies, for similar reasons, and just extrapolated IOB information from 4 pumps by recording the IOB values for 5 hours at 5 minute intervals, providing 61 data points. I approximated the IOB functions that way.

The IOB calculation interfaces directly with Nightscout for OpenAPS users, like Dana, Scott, and Toby. The Nightscout web interface only uses JavaScript and does not include the math.js library, so exponential functions at this point of time are not a logical option. I have no control over what libraries are included in the Nightscout web interface, leaving linear and least squares autoregression with a polynomial approximation the most viable option given the limitations of Nightscout. The RMSE (root mean square estimate) is acceptable on the approximations, and I am not going to worry about it.

According to Clinical Performance of Three Bolus Calculators in Subjects with Type 1 Diabetes Mellitus: A Head-to-Head-to-Head Comparison, Animas and Roche are the most accurate representations of IOB.

[Monica C.]

From Method for Automatic Adjustment of an Insulin Bolus Calculator: In Silico Robustness Evaluation under Intra-Day Variability (http://www.cmpbjournal.com/article/S0169-2607%2815%2900027-9/abstract), IOB can be calculated (curvilinear formula):



where IOB is the IOB remaining in units, B_sub_i is administered bolus at time T_sub_B_sub_i, n is the number of boluses previously administered, and tau is the predefined exponential decay constant (e.g. tau = 5 h). t is the current time. All time quantities are in terms of hours.

[Bert Berla]

If exponential functions aren't part of nightscout, can't you just add one in as a table? You only need 5-minute time resolution, so it wouldn't take much data.

Also, the paper you cited above about the different algorithms didn't suggest that Medtronic's was less accurate, just that it was less successful in avoiding hyperglycemia because it corrected to a higher target. All three curvilinear IOB representations were equally good. It's a shame they didn't compare to a linear representation, although their protocol of correcting at 2 hours post-meal (as opposed to at more random intervals like a real user might) would also probably mask any differences.

Again, I'm not trying to bash your choice, but I do have to say I'm not comfortable using an equation I don't understand that seems to basically have the same shape as another to calculate an insulin dose. I'd much rather use an equation I fully comprehend and that has a basis in chemical kinetics.

[spiph23]

Hi Bert,

I've been trying to read through the back and forth, but I think I'm missing something.  When you wrote, "I'm not comfortable using an equation I don' understand," what equation are you referring too exactly?"  Also, did I understand you correctly when you said you would rather have a table, did you mean a table for duration of insulin activity that has the percentage of insulin remaining listed in 5 minutes increments? 

Sorry for my confusion!

Toby 

[Monica C.]

First, exponential functions are not a part of NS, as the math.js library is not included, although you can fork a copy of Nightscout and include the math.js function and program the exponential function rather easily. You can approximate the exponential function without adding the math.js library, similar to what I did with the curvilinear percent IOB left curves provided with all modern "smart" insulin pumps that stay constant regardless of initial insulin dose/bolus amount. The problem with approximating the exponential function, is you would have to create a matrix for each DIA with IOB remaining, based on each respective insulin dose, ranging from 0-20 units, likely in 0.25 unit increments.

So, the matrix would look something like:
                                       0.25           0.50              1.00             1.25             1.5   .... 10  ...       20
_____________________________________________________________________________________
DIA_5h =    0 min     |    0.25           0.50              1.00             1.25             1.5   .... 10  ...       20
                   ....
                   150 min |    0.125         0.25              0.50             0.625            0.75     0.50         10
                    ...        
                   300  min|      0               0                    0                   0                  0          0             0


And yes, while this is more dynamic than the curvilinear IOB percentage calculator that is provided with modern insulin pumps, IOB in itself is only an approximation, and it cannot be reliably taken from euglycemic insulin clamp studies, even with the respective fast-acting insulins:

1. Because there are limited studies, with different set parameters, making it impossible to get a reliable approximation in between insulins.
2. IOB is a very strange function that can only be approximated (see this: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3045234/figure/fig1/ ) (Original article here: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3045234/ )
3. No IOB equation will have a true basis in chemical kinetics because IOB is a systematic concept of human metabolism and pharmacokinetics and pharmacokinetics and is accounted for and approximated by dosing "healthy" subjects with high doses of insulin, because of the inability to extrapolate the endogenous glucose production when a "basal insulin rate" is set to 0 in healthy subjects. Also, IOB is based in terms of two variables: time and initial insulin dose. It does not account for variability in absorption or diurnal variations, or accumulated IOB, and most importantly, it is a system of several variables that IOB is dependent on.

Either way, using curvilinear percent IOB representation based as a function of time, or exponential decay, as a function of both time and initial bolus dose(s), is merely a rudimentary approximation that does not account for things as a dynamic system. That is why I am programming Hovorka's models (again), which are based on metabolic modeling of at least 5 separate systems working in tandem, in compartments. I am designing my controller based off of Hovorka's models too. I have pretty much passed the part of being overwhelmed by the information for the mathematical models, and I hope to be done with the mathematical modeling in a couple of weeks. It's gotten to a point where it is common sense to me now conceptually. I start school in a couple of weeks, and during those 8 weeks that I will be taking a course, I hope to get my controller finished. A month after I finish the controller, I hope to be done with the interfacing and "go live".

As for the paper, I knew that the results were skewed, and that was a poor example for me to use.

[Bert Berla]

Hi Toby, I'm referring to the set of curvilinear equations for representing IOB. It's just too complicated for my taste, and I don't see the benefit relative to something MUCH simpler like a linear decay, given all the uncertainties about which model is actually accurate. But I guess this is everybody's choice to make for themselves.

Yes, that's the sort of table I'm talking about.

[Scott Leibrand]

Bert,

Our curvilinear IOB decay models are very much based in actual physiology (pharmicookinetics/pharmicodynamics/whatever it is).  A number of studies have empirically determined the insulin activity curves of various insulins, usually for purposes of comparing them to each other.  The key insight for our purposes is to realize that IOB is nothing but the integral of the insulin activity rate (and conversely, insulin activity is the slope, or derivative, of the IOB curve).  So to derive a decent IOB curve, you can start with an estimate of the insulin activity curve (starting at zero, peaking at 60-90 minutes, and decaying to negligible levels over 3-5 hours), and then integrate that to get a corresponding IOB curve. That's exactly what the IOB formula in wip/iob-cob does.  It starts with an insulin activity curve that (for DIA=3) starts at 0 and rises linearly to a peak at 75m, and then falls linearly back to zero at 180m. Taking the integral of those two lines (in my case, by calculating it algorithmically in Excel and curve-fitting to find the two R^2=1.000 quadratic equations) gives you the corresponding IOB curve.

Obviously insulin activity is more of a curve than an inverted V, so you can derive more complicated formulas by starting with a function that more closely matches experimental insulin activity curves, but at the end of the day your IOB curve is going to look a lot closer to the one Nightscout and most insulin pumps use (an S-curve) than a linear IOB decay function.

It's also worth noting that in #DIYPS we started with a linear IOB decay function, and only introduced the S-curve IOB decay once we saw that the predictions were always wrong in the 30-60m interval before insulin activity really picks up.  We tested both linear and S-curve IOB decay, and found a definite improvement in predictions with the curvilinear function.

-Scott

--

You received this message because you are subscribed to the Google Groups "OpenAPS Dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openaps-dev...@googlegroups.com.
To post to this group, send email to opena...@googlegroups.com.
Visit this group at http://groups.google.com/group/openaps-dev.

To view this discussion on the web visit https://groups.google.com/d/msgid/openaps-dev/af9b263c-9cb3-4def-8a76-72d5744300c1%40googlegroups.com.

[Monica C.]

While I do take a lot of interest in how you derived the IOB, Scott, and I used a similar approach for the first set of equations I did for the Apidra and Novolog/Humalog modeling, Walsh in the previous paper said that IOB is based off of glucose infusion rate in the last post of the paper I referenced.

I don't mean to cause any further disagreement, as obviously your technique has worked and is working. But, after reading that paper, I really questioned my work.

[Monica C.]

That is why I extrapolated the data from the insulin pumps. This is a graphical representation of IOB. Check out the Axes, particularly the y one. The glucose infusion rate used in euglycenic clamp studies, at high dosages, 0.3-0.4 U/kg, is used to avoid less accurate results because of the unknown endogenous glucose infusion rates by the human body. Anyways, IOB is not a function of insulin appearance rate in the plasma, although it is dependent on that function. It is a function of the glucose infusion rate of a study at 0.3-0.4 U/kg of fast acting insulin of a euglycenic clamp study, plus the endogenous glucose production/infusion rate, that has to be extrapolated (Please see: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3045234/figure/fig1/ ) . The entire paper is a good read. All in all, IOB is an approximation, and is not strictly based off of pharmacokinetics and pharmacodynamics.

[Bert Berla]

Scott, I think your explanation makes sense, and I'm glad it's working. I'm gonna leave it there.

[Monica C.]

First, I apologize for the 2 deletes. I made 2 minor errors in my equations that I wrote that I copied and pasted from MathType that might cause some confusion.

Bert, I guess what you have been wondering about is the relationship between the IOB function and the curvilinear percent IOB remaining equations, all along but that was unclear to me. Here is how it is derived. I thought about it a little bit when I was driving and all it took was doing basic calculus in my head and relating things graphically. Of course, I checked my work.

I defined the IOB function as IOB, which is essentially either the glucose infusion rate of a euglycemic clamp study when 0.3-0.4 U/kg is given, interpolated to a mathematical function representing what is going on graphically, or the IOB function I provided in a previous post which is an exponential decay equation.

Either way, the curvilinear percent remaining IOB function is based off of pharmacokinetics and pharmacodynamics, but not in a direct or "pure" way. That's one of the reasons why I am doing the Hovorka models. Plus, my medical situation is too complex to trust anything less than a controller based on control theory. But part of that is the engineering student in me.

So, the Percent remaining IOB equation f(x) equals: