LoopKit / Loop

An automated insulin delivery app for iOS, built on LoopKit
https://loopdocs.org
Other
1.51k stars 1.3k forks source link

Dosing algorithm #533

Closed dm61 closed 6 years ago

dm61 commented 7 years ago

This is inspired in part by the discussion in #388 where @ps2 and others noted that longer-tail exponential insulin absorption curves may expose some shortcomings of the dosing calculation, which is currently based on BG predicted at t=DIA in the future. Currently, the amount of bolus insulin suggested after a meal is entered is found as follows: image Such DIA-based dosing works well in many cases. However, it ignores the effect of the insulin delivered on BG between now (t=0) and DIA, which may result in too high or too low amount suggested. To address this issue, one may consider an iterative dosing algorithm that would search for an insulin dose so that BG predicted after the dose is delivered does not go below a target. Such iteration can in fact be captured by a relatively simple modification of the above dosing calculation: image The argument of the min function has the meaning of the amount of insulin required to drive BG to the target value at time t. The 1-IOB(t) term in the denominator represents the fraction of to-be-dosed insulin expended from 0 to time t. By choosing the minimum of the doses suggested over the time interval from t=0 to t=DIA we guarantee that the predicted BG, after the dose is delivered, will just touch the target value BGtarget at some point in time between now and DIA (not necessarily at DIA).

To illustrate how the suggested dynamic dosing calculation would work, consider a simple scenario: ISF=50, CR=10, DIA=6 hours, and meal of 50g that absorbs slowly at constant rate over a period of 5 hours. Since the entire meal is absorbed within DIA, the current dosing calculation would suggest 50/10 = 5 U bolus. A simple simulation shown below indicates that this dose is too large for the given meal, with BG dropping from initial 100 mg/dL to as low as 56 mg/dL (neglecting any low temping, which is not simulated here, but which would be insufficient to prevent the low) too_large_dosing_at_dia In the same example, the dynamic dosing suggests 3.7 U, and the simulation shows that BG just touches the target (100 mg/dL) back at around 200min, which is t=140min after the dose is delivered. dynamic_dosing Since this initial dose is insufficient to cover all carbs, BG eventually drifts up to a value above the target, but this would be corrected by Loop high-temping (which is not included in this simple calculation).

I think the suggested modification of the dosing algorithm could apply to both boluses and temps, and I think everything necessary to perform the calculation is available in DoseMath.

A concern is that the dynamic dosing described above may be too conservative, especially in view of the Loop dynamic carb absorption algorithm, which includes a 50% extension to a user-entered carb absorption time. The dynamic dosing formula could be modified to be more aggressive in at least two ways: by replacing BGtarget with BGmin (or perhaps even BGguard just for boluses), or by taking into account the effect future zero temping could have on predicted BG, thus enabling sort of dynamic super-bolusing as a part of the algorithm. The formula below includes both of these modifications: image

elnjensen commented 7 years ago

Interesting idea. I think I understand what you're suggesting here, but I have a question, which may be mostly about notation. Your IOB(t) term doesn't actually have units of insulin, right? Rather, it's the fractional area under the insulin action curve at time t. So it's the ratio of the integral of that curve from 0 to t to that of the integral from 0 to DIA. Is that correct?

Edit: caught up on reading all the comments in #388 and I see that the usage there is indeed that IOB(t) is unitless, descending from 1 to 0. So this all makes sense to me, given that definition.

elnjensen commented 7 years ago

Another thought on this: if you start out at or below target, then this formula will recommend a zero (or negative! πŸ˜‰) bolus, since the minimum of the series will occur at the first time step. If you're sitting on a nice flat BG line just slightly below target, then zero bolus for a meal probably isn't what you want. So maybe it would make more sense to start that search for the minimum at some time > 0, though I'm not sure where, exactly. Maybe starting at/around the peak of the action curve (esp. if we parametrize the curves as suggested in #388, so this peak value is easily available)?

ps2 commented 7 years ago

@elnjensen No, that's not correct. Bolus dosing allows full bolus unless min bg is below guard.

dm61 commented 7 years ago

@elnjensen the formulas in my top post are just suggestions on how the core dosing calculation could be modified to a more general form beyond the current DIA-only based calculation; it would be necessary to add provisions to these formulas to account for various conditions, including pre-bolusing, starting form an initial bg below target, etc, just as is currently done in Loop. Many options could be considered, but would be good to stay away from ad-hoc conditions as much as possible. For example, BGtarget in the formula above could actually be BGtarget(t) moving from BGguard at t=0 to BGtarget at t = DIA, which would allow the entire Loop logic to stay as is, including bolus recommendation allowed as long as bg is above guard, as noted by @ps2

elnjensen commented 7 years ago

@dm61 Right, that makes sense. @ps2 Sorry if I wasn't clear - I wasn't trying to describe current Loop behavior, but rather noting what would happen if the formula that @dm61 gives above were used, unmodified, to set bolus behavior. As @dm61 just noted in his most recent comment, there would have to be a more complicated set of conditions considered in order to use a dosing algorithm like this, one of which is the case I was noting of starting below target.

ps2 commented 7 years ago

Thanks @dm61! This has been a potential issue for a long time (bolus for eventualBG results in near term low). I had been keeping it in the back of my mind, but I haven't taken issue, since we don't see it very often and when we do, it's very short term, and not much under target. With a faster peak insulin curve, it will be more pronounced, and I think there are already users who prebolus that probably see a deeper dip for some meals. (We're not great at pre-bolus).

My comment in that last ticket was partially a thinly veiled attempt to get you thinking about this. ;) The suggestion for an iterative approach was because that's all I could think of. I am wondering if there is a more direct way of computing this. Perhaps not. I'm guessing that any direct approach would be quite dependent on the actual carb absorption and insulin effect functions, and would break if we wanted to switch those out. So maybe iterative is still the best approach.

dm61 commented 7 years ago

@ps2 the calculation suggested above is direct - it does not require iteration

ps2 commented 7 years ago

Oh, I obviously misunderstood it. Now I see. And we can use BGpredicted sans insulin, and the IOB curve directly. Efficient and clever! Still convincing myself there aren't edge cases, but it sounds great. Also still considering the modifications for replacing target with the min of the range.

dm61 commented 7 years ago

How about sliding the target value in that equation linearly from BGguard at t=0 to BGtarget at t=DIA? That would keep the existing Loop logic as is: it would allow bolus suggestion as long as the initial point is greater than BGguard, and it would ultimately drive BG to BGtarget just as Loop does now.

ps2 commented 7 years ago

I've been playing with this, and so far it's looking good. A couple notes. At t = 0, we have a divide by zero; I can deal with that; it doesn't really make sense to look for a scale at that time. Second, I get an updated value of 3.0 when I run the equation, but I'm guessing you might have a slightly different insulin curve in operation.

https://github.com/ps2/LoopIOB/blob/master/Dosing%20algorithm%20%23533.ipynb

dm61 commented 7 years ago

Yeah, calculation at zero does not make sense, start from the next element of the array, I'll correct those equations. I think my insulin curve was the same (td = DIA = 6h, tp=75min). Oh, I've just realized I had a typo in my meal example - it actually extended over 5 hours not 6; will edit that typo now. That might explain the difference 3.0 in your simulation versus my 3.7. edit just checked: I get 3.0 for the 6-hour meal.

jeremybarnum commented 7 years ago

It's pretty clear to me this is the way forward. I think the next piece is a framework for safely taking into account "dry powder" from the ability to low temp in the back of the DIA, which allows for more aggressiveness in the front of the DIA. It should be something like "I am comfortable relying on being able to cut the basal in half, but I want to keep the rest for errors". Sort of a souped up version of low BG guard, integrated with the above.

dm61 commented 7 years ago

I've played some more with the suggested dosing approach. This is a version that seems to work well: image where target BGt(t) starts from BGguard at time of dosing (t=0) and ends at BGtarget at t=DIA, while aggressiveness parameter alpha blends in effects of potential future zero temping, thus allowing for a more aggressive bolus suggestion up front. As mentioned by @jeremybarnum, alpha = 0.5 might be a default choice. Or, aggressiveness alpha could be another user-adjustable parameter.

A Matlab script used to perform simulation and test various options can be found here: https://github.com/dm61/BasicLoopSimulator.

Below are couple of examples to illustrate typical differences between current DIA-based dosing and the above dynamic dosing.

Scenario 1: 50g meal absorbed over 5 hours, ISF=50, CR=10, DIA=6, BGtarget=100, BGguard=70, 20min pre-bolus. Dynamic dosing (right) delivers a smaller bolus, and BG never drops below BGguard 5h_meal_08062017

Scenario 2: 50g meal absorbed over 2 hours, ISF=50, CR=10, DIA=6, BGtarget=100, BGguard=70, alpha =0.5, 20min pre-bolus. In this case, dynamic dosing safely delivers a larger bolus, which reduces the spike. 2h_meal_08062017

jeremybarnum commented 7 years ago

This looks fantastic @dm61. I'm looking forward to playing with it tomorrow and developing some intuition. Thanks for your efforts and clear explanations.

jeremybarnum commented 7 years ago

Love the simulator! This one is awfully promising looking: Fiasp with no prebolus but aggressive setting:

screen shot 2017-08-08 at 1 39 21 am

And I wonder whether in reality the "true" carb absorption curve has some delay in it, which might not be that different from the Fiasp delay. Interesting to note this - Fiasp with 20 minute pre-bolus and aggressive alpha sends you uncomfortably low - even with fast burning carbs. When you buy it, they make a special call to tell you not to pre-bolus...

screen shot 2017-08-08 at 1 44 24 am
jeremybarnum commented 7 years ago

@dm61 can you validate that this is an accurate plain english explanation of this code:

% dynamic dosing algorithm: zero-temping BG impact functions
low_temp_scale = -(min_temp_rate/60)*ISF; % asymptote of the bg rate of change due to zero temping
% bg impact of zero temping as a function of time
BGTempImpact = @(t) low_temp_scale*(S/td)*(tau/td)*tau*...
    (exp(-t/tau).*(-6*tau^2+2*tau*(td-2*t)+t.*(td-t))+...
    6*tau^2-2*tau*td-2*tau*t+td*t);

Is it essentially saying: the zero temp can be seen as a mini-bolus = basal rate/60 for every minute minute it's in place with a total eventual BG impact of = basal rate/60*ISF.

Then taking that quantity and turning it into a vector of BGIs using the insulin activity function, which are then proliferated through the simulation?

ps2 commented 7 years ago

This is nice work. Thanks @dm61!

I understand the desire to have a dosing algo that includes all of the above features, but it does come at some expense of ease of understanding. I first want to get a better idea of the current propensity for overbolusing or predicting an overbolus with the exponential insulin curves. We have been running the exponential insulin curves for a week or so now, but surprisingly haven't hit the large-meal-bolus-causes-forecasted-low situation with any severity.

The next step would be to implement a bg guard (or low glucose suspend as it will likely be renamed) back-off approach as discussed at the beginning of this ticket. I think the aggressive bolus would have to live on an experimental branch for validation before being considered.

jeremybarnum commented 7 years ago

Despite my enthusiasm, I do agree that the more aggressive bolus should probably not be a mainstream feature. It might be nice if the necessary modification to the code were relatively straightforward so that those who wanted to use it could access it without worrying too much about breaking other things. Ultimately I think there should be some relatively simple approximations for how much to scale up the bolus as a function of CA. People who really wanted to could simply override the bolus manually and then Loop would naturally low temp as a result.

dm61 commented 7 years ago

@ps2, @jeremybarnum thanks for the feedback - I agree aggressive bolus should be considered highly experimental.

@jeremybarnum question:

Is it essentially saying: the zero temp can be seen as a mini-bolus = basal rate/60 for every minute minute it's in place with a total eventual BG impact of = basal rate/60*ISF.

That's how I first approached it, and that's how I think Loop (and OpenAPS) are dealing with temp basals, numerically, as an array of tiny 5min doses. For the purposes of the above dosing calculations, turns out the effect of zero temping on BG can be found analytically by integration of the 1-IOB(t) function. The end result is a bit messy (that's the BGTempImpact function in the simulator), but it can be computed directly, so it's not too bad. A simple piecewise linear approximation is also an option, as illustrated below.

Starting from a steady-state BG, with no other disturbances, the graphs below show DeltaBG in response to zero temping, for the case when basal rate is 0.5 U/h, ISF = 50, rapid insulin tp = 75, td = 5 hours. zero_temping The simulation curve shows the response for zero-temping froom t=0 to t=DIA, while the analytical curves assume persistent zero temping starting at t=0. The slope of BG ultimately approaches (rate*ISF) per hour, as one would intuitively expect. Fun stuff. [edit: corrected a typo in slope statement, thanks @jeremybarnum]

jeremybarnum commented 7 years ago

Very cool - thanks for the answer and explanation. Your point about the curve approaching rate*ISF per hour (I think rate/ISF is a typo?) is intuitive in hindsight, but I had missed it! Today I spent way too much time trying to implement this in Excel the ugly way just to build some intuition. I need better skills! πŸ˜„ Tomorrow's task will be to understand your "messy" analytical implementation of the low temp effect. Fun stuff indeed.

trixing commented 7 years ago

Has any you do some simulations of the effect of wrongly inputting carbs? Especially with kids this seems potentially dangerous to take the carbs by face value and remove any chance of correcting afterwards by low temping?

Also the effect of time delay (entering 15 minute before, during, after lunch, or not at all) would be interesting especially with respect to Fiasp.

The model looks perfect though for perfect data :)

On Thu, Aug 10, 2017 at 4:29 AM, jeremybarnum notifications@github.com wrote:

Very cool - thanks for the answer and explanation. Your point about the curve approaching rate*ISF per hour (I think rate/ISF is a typo?) is intuitive in hindsight, but I had missed it! Today I spent way too much time trying to implement this in Excel the ugly way just to build some intuition. I need better skills! πŸ˜„ Tomorrow's task will be to understand your "messy" analytical implementation of the low temp effect. Fun stuff indeed.

β€” You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/LoopKit/Loop/issues/533#issuecomment-321433204, or mute the thread https://github.com/notifications/unsubscribe-auth/AAE2ajPgvKf5ibwq1H-ouDDppTEAkJEhks5sWmsMgaJpZM4OlzHM .

-- Jan Dittmer jdi@l4x.org, http://l4x.org

jeremybarnum commented 7 years ago

@trixing it's definitely potentially dangerous (what you are referring to is equivalent to an alpha setting of 1 - no room for error, relying on zero temping to the maximum extent possible) which is why I think we all agree that this is highly experimental and not a mainstream option. Any modification to the algorithm should be resilient to the real world where all inputs and assumptions are subject to error and this applies equally to kids and adults.

The only reason it's even being discussed in my view is that if we are updating insulin curves to reflect small amounts of activity in a tail between, roughly, DIA = 4 and new exponential curve td = 5 or 6 (60 minutes), and all dosing decisions are based on a prediction horizon = td, it seems like it could be a tiny bit conservative in some cases by failing to rely on the ability to low temp late in the DIA to offset that tail which is of the order of 5-10% of the bolus. And in practice, if the new curves are "true", then to the extent people are using shorter DIAs now they are doing this anyway* because in effect they are ignoring the tail of the insulin - but in a more dangerous way, because the tail will come as a surprise.

Personally I am motivated in part by the desire to avoid people being unhappy with improvements in accuracy/reality because the improvements are removing unintentional aggressiveness which might be, for some people, producing good outcomes, albeit in an un-transparent way.

ps2 commented 7 years ago

@dm61 The difficulty we are having in implementing this is that we want to allow a certain threshold to avoid (suspend threshold), while targeting another separate target for eventual bg. I'm wondering if there is any way to implement a similar algorithm with that in mind.

dm61 commented 7 years ago

@ps2 In the dosing equation of the top post I suggest you replace the fixed BGtarget with a function of time BGt(t) that slides from suspend threshold to eventual target over DIA. This can be done in any number of ways. Two approaches I've tried in simple simulations are illustrated below. The first option is to just linearly move from suspend to eventual. The second option is to stay at suspend threshold until some fraction x of DIA and then linearly move to eventual target. The fraction x is somewhat arbitrary (could be 0.5 for example), or we could select it based on some fraction of IOB expired by that time. I like that second option better because it is less conservative, but still very safe.

bgt

jeremybarnum commented 7 years ago

Do we agree this approach should apply equally to bolusing and temp basal recommendations? I can imagine doing it for temp basals is a lower priority because for a range of reasons it seems less likely to create problematic lows. But in theory the framework should apply equally, so that people who want to run with very high max basal settings and basically operate in "unannounced" mode can benefit from this enhancement too. Again, in practice, it's probably an unimportant edge case, but I thought worth clarifying what everyone thinks the "ideal" state is.

trixing commented 7 years ago

I believe they both should be exactly the same. I'd actually propose to change DoseMath to have the formula only once.

Currently there is a subtle difference I think that Bolus targets the upper target range, whereas TempBasal uses the middle of the range. I think otherwise they try to do the same safety checks.

It would be nice to split up DoseMath to a function which calculates the necessary Insulin amount and then have the Basal / Bolus calculation use that amount to either suggest a single bolus or a rate and duration respectively.

On Tue, Sep 5, 2017 at 4:54 PM, jeremybarnum notifications@github.com wrote:

Do we agree this approach should apply equally to bolusing and temp basal recommendations? I can imagine doing it for temp basals is a lower priority because for a range of reasons it seems less likely to create problematic lows. But in theory the framework should apply equally, so that people who want to run with very high max basal settings and basically operate in "unannounced" mode can benefit from this enhancement too. Again, in practice, it's probably an unimportant edge case, but I thought worth clarifying what everyone thinks the "ideal" state is.

β€” You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/LoopKit/Loop/issues/533#issuecomment-327201502, or mute the thread https://github.com/notifications/unsubscribe-auth/AAE2alMHK9KckeDxwdADPLcqUJTeFTjkks5sfWC6gaJpZM4OlzHM .

-- Jan Dittmer jdi@l4x.org, http://l4x.org

dm61 commented 7 years ago

Applying this approach to temp basals would be fine but the benefits, if any, would likely be relatively small. Loop reevaluates temps every 5 minutes, and so it already has built-in corrective mechanisms, which do not exist for boluses. In basic idealized simulations, applying the same algorithm to temp basals results in smoother looking temp patterns, but the resulting bg responses look very similar to the responses with the existing Loop algorithm. With dynamic carbs, I've raised my temp limit considerably, and have not noticed any issues in actual practice. So, imo, applying this approach to temps would be a low priority.

scottleibrand commented 7 years ago

This sounds like the kind of thing that would (only) be relevant if at some future time you decided to implement an SMB-like approach. At that point, having insulin requirements calculated independently of the delivery method (user-verified meal bolus, temp basals, or microboluses) would be very useful. But I don't think microboluses are even on the roadmap yet, so I'm guessing this will probably be a low priority until after things like automatic sensitivity detection are done. (Speaking only as an interested observer.)

ps2 commented 7 years ago

Hey guys, it's already in the code and working.😜 On dev, both bolus and temp basals use this algo.

trixing commented 7 years ago

How is a temp with a high maxbasal limit different from a bolus in practice?

I believe fundamentally the insulin needs of the body should be the same in both cases. Sure you might want to do different 'aggressiveness ' optimizations in both cases but even their having a higher than necessary temp basal is a trick useful similar to higher bolus.

(Fwiw Microboluses are just a small step to add which I actually had working for a while in Loop now mostly saving me from very high temp rates)

On Wed, Sep 6, 2017 at 1:41 AM Scott Leibrand notifications@github.com wrote:

This sounds like the kind of thing that would (only) be relevant if at some future time you decided to implement an SMB-like approach. At that point, having insulin requirements calculated independently of the delivery method (user-verified meal bolus, temp basals, or microboluses) would be very useful. But I don't think microboluses are even on the roadmap yet, so I'm guessing this will probably be a low priority until after things like automatic sensitivity detection are done. (Speaking only as an interested observer.)

β€” You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/LoopKit/Loop/issues/533#issuecomment-327333306, or mute the thread https://github.com/notifications/unsubscribe-auth/AAE2aiSBVlDl7iMEUvnzz0kKjI2RAtxHks5sfdwTgaJpZM4OlzHM .

-- Jan Dittmer jdi@l4x.org, http://l4x.org

scottleibrand commented 7 years ago

If you could do arbitrarily short temp basals, there wouldn't be a difference. But because the minimum temp basal length on Medtronic pumps is 30m, and since we have to assume that we'll lose connectivity and any high temps will run to completion, we can only set high temps high enough to deliver the required insulin over 30 minutes. With microboluses, you can deliver the a sizeable fraction of the required insulin every 5m with each new BG, such that you can deliver all of it in 15-20m (up to twice as fast as with high temps) while still having the ability to stop if a sudden BG downtick indicates you maybe didn't need all the insulin you first thought.

I'll leave it at that to avoid rabbit-holing any further off-topic, but feel free to jump over to Gitter if you'd like to discuss further all the nuances of how we use SMBs in OpenAPS.

dm61 commented 7 years ago

Oh, wow, @ps2, that's truly amazing 😲, I thought it only applied to bolus calculations in dev. Well, that wraps up this discussion I guess. πŸ‘

rsilvers129 commented 7 years ago

I know OpenAPS only sets temp basals that it calculates are ok to run to completion, but does Loop also do that? Or is Loop willing to set very high ones and then counts on being able to stop it five minutes later? If the latter, then just sending it as a bolus would seem better.

trixing commented 7 years ago

Thanks for the explanation of the difference. I still don't see how that changes the required (not given) insulin, but as you said, not worth discussing this further.

On Wed, Sep 6, 2017 at 4:22 AM, rsilvers129 notifications@github.com wrote:

I know OpenAPS only sets temp basals that it calculates are ok to run to completion, but does Loop also do that? Or is Loop willing to set very high ones and then counts on being able to stop it five minutes later? If the latter, then an SMB feature would be nice.

β€” You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/LoopKit/Loop/issues/533#issuecomment-327355800, or mute the thread https://github.com/notifications/unsubscribe-auth/AAE2ahdoC5g7Gnc4OwDNaPPU2rLAh7_4ks5sfgHngaJpZM4OlzHM .

-- Jan Dittmer jdi@l4x.org, http://l4x.org