Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dosing algorithm #533

Closed
dm61 opened this issue Jul 27, 2017 · 34 comments
Closed

Dosing algorithm #533

dm61 opened this issue Jul 27, 2017 · 34 comments

Comments

@dm61
Copy link
Contributor

dm61 commented Jul 27, 2017

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
Copy link
Contributor

elnjensen commented Jul 30, 2017

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
Copy link
Contributor

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
Copy link
Collaborator

ps2 commented Jul 31, 2017

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

@dm61
Copy link
Contributor Author

dm61 commented Jul 31, 2017

@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
Copy link
Contributor

@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
Copy link
Collaborator

ps2 commented Aug 1, 2017

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
Copy link
Contributor Author

dm61 commented Aug 1, 2017

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

@ps2
Copy link
Collaborator

ps2 commented Aug 2, 2017

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
Copy link
Contributor Author

dm61 commented Aug 2, 2017

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
Copy link
Collaborator

ps2 commented Aug 2, 2017

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
Copy link
Contributor Author

dm61 commented Aug 2, 2017

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
Copy link

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
Copy link
Contributor Author

dm61 commented Aug 6, 2017

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
Copy link

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
Copy link

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
Copy link

@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
Copy link
Collaborator

ps2 commented Aug 8, 2017

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
Copy link

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
Copy link
Contributor Author

dm61 commented Aug 10, 2017

@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
Copy link

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
Copy link

trixing commented Aug 10, 2017 via email

@jeremybarnum
Copy link

@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
Copy link
Collaborator

ps2 commented Sep 2, 2017

@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
Copy link
Contributor Author

dm61 commented Sep 2, 2017

@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
Copy link

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
Copy link

trixing commented Sep 5, 2017 via email

@dm61
Copy link
Contributor Author

dm61 commented Sep 5, 2017

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
Copy link

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
Copy link
Collaborator

ps2 commented Sep 6, 2017

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

@trixing
Copy link

trixing commented Sep 6, 2017 via email

@scottleibrand
Copy link

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
Copy link
Contributor Author

dm61 commented Sep 6, 2017

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
Copy link

rsilvers129 commented Sep 6, 2017

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
Copy link

trixing commented Sep 6, 2017 via email

@ps2 ps2 mentioned this issue Sep 20, 2017
@dm61 dm61 closed this as completed Dec 18, 2017
ps2 added a commit that referenced this issue Nov 28, 2022
* Prevent marking onboarding as finished when pausing

* Update comment for clarity
ps2 added a commit that referenced this issue Nov 28, 2022
* [LOOP-4205] always display time changed alert when detected. retract when needed (#523)

* Adjust method to match updated protocol (#524)

* [LOOP-4289-4348] Critical alert modal and banner update (#525)

* initial pass at adding go to settings button in the critical alert modal

* remove foreground content from the riskMitigationAlert, since it is handle by the AlertManager

* added placeholder text for the alert premissions disabled banner

* add dismiss button to alert permissions disabled warning

* acknowledge alert with close option

* response to PR comments

* clean up

* more cleanup

* updated notification permissions banner

* updated style of alert permission disabled alert

* copy update

* [LOOP-4348] adjusted layout of notification permissions disabled banner (#527)

* [LOOP-4289] Also issue iOS notification (#528)

* issue notification permission alert when in background and retract in-app alert as needed

* remove unused code

* response to PR comments

* [LOOP-4348] notification permission banner for narrow display (#529)

* Fix crash in ZipArchiveTests (#530)

* [LOOP-4400] made status highlight container for CGM and pump pill the same size (#531)

* [LOOP-4349] Temporary mute alerts (#526)

* added placeholder UI to configure mute alerts duration

* checkpoint

* checkpoint

* initial commit at the alert muter

* minor clean-up

* removed unused file from project

* corrected predicate syntax

* all alerts at least vibrate

* checkpoint

* insteadOf = https://github.com/
updating unit tests

* checkpoint

* tighter coupling between alert manager and in-app/user notification alert issuers

* corrected rescheduling loop not running notifications and use timers for UI

* make the mute end period timer internval to the alert muter

* corrected typo

* updated unit tests

* added duration to temp mute alerts setting

* minor clean-up

* using source of truth for last loop date

* renamed Issuer -> Scheduler for in-app and user notification

* [LOOP-4349] corrected loop not looping rescheduling (#532)

* corrected identifier for loop not looping notifications

* updated test to run for both DIY and Tidepool

* COASTAL-1076 Support PumpManagerUI pausing onboarding (#533)

* Prevent marking onboarding as finished when pausing

* Update comment for clarity

* Use new iOS 16 method to notify VC of supported orientation change update (#534)

* [COASTAL-1124] Time alert copy update (#535)

* updated copy of time alert

* using device model

* Disable mute alerts feature (wip) and fix tests

Co-authored-by: Nathaniel Hamming <nate@hmtconsulting.org>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

7 participants