A Quadruply-Asymmetric Sigmoid to Describe the Insulin-Glucose Relationship during Intravenous Insulin Infusion

For hospitalized patients requiring intravenous insulin therapy, an objective is to quantify the intravenous insulin infusion rate (IR) across the domain of blood glucose (BG) values at a single timepoint. The algorithm parameters include low BG (70 mg/dL), critical high BG, target range BG limits, and maintenance rate (MR) of insulin infusion, which, after initialization, depends on rate of change of blood glucose, previous IR, and other inputs. The restraining rate (RR) is a function of fractional completeness of ascent of BG (FCABG) from BG 70 mg/dL to target. The correction rate (CR) is a function of fractional elevation of BG (FEBG), in comparison to elevation of a critical high BG, above target. IR = RR + CR. The proposed mathematical model describing a sigmoidal relationship between IR and BG may offer a safety advantage over the linear relationship currently employed in some intravenous glucose management systems.


INTRODUCTION
Hyperglycemia, defined as an elevation of concentration of blood glucose (BG) above a defined normal range, is a common condition among hospitalized patients, and is an important factor predictive of outcomes [1].In at least some settings and patient populations, treatment of hyperglycemia with intravenous insulin infusion may improve outcomes [2].
Enthusiasm for tight glycemic control in the hospital was sparked by a single center study from a surgical intensive care unit showing reduction of mortality by a protocol with glycemic target range of 80-110 mg/dL [3].However, a subsequent large multicenter randomized trial of intensive vs. standard glucose control in mixed medical-surgical intensive care units showed that attempted attainment of target glucose range of 80-110 mg/dL was complicated by an increased rate of severe hypoglycemia and increased mortality [4].At present, recommended target ranges for optimal control of BG in the intensive care unit have undergone revision and may vary according to patient population [5][6][7].For patients having hyperglycemic emergencies, specifically hyperglycemic hyperosmolar state (HHS) or diabetes ketoacidosis (DKA), glycemic targets generally are higher than those for other critically ill patients, partly because of concern that rapid correction of hyperglycemia could contribute to severe and potentially life-threatening complications, e.g., cerebral edema [8].On general wards as well, although there is less evidence from randomized trials than in the intensive care unit, based on available evidence, specific guidelines for target range control have been suggested [9,10].
The principal risk of insulin treatment is hypoglycemia, generally defined as glucose < 70 mg/dL [11,12].Hypoglycemia is associated with a variety of co-morbidities and may be independently associated with increased hospital mortality [13][14][15][16][17]. Hypoglycemia is a specific and well recognized risk of intravenous insulin infusion, reported to be associated with lower glycemic targets, and potentially attributed to presence of renal failure, use of higher amounts of insulin, design of protocols specifying continuation of insulin despite low glucose readings < 100 mg/dL, reduced frequency of testing, and staff non-adherence to protocol design [18][19][20][21][22]. Since hypoglycemia predicts adverse outcomes and may produce immediate neurologic symptomatology, insulin therapy aiming at specific glycemic targets should be designed to minimize risk for hypoglycemia.Furthermore, in the treatment of patients having hyperglycemic crises and for study and treatment of various other hyperglycemic conditions during hospitalization, a strategy should exist for attainment of specific glycemic targets.
The long term goal of this work is to develop a programmable algorithm for treatment and prevention of hyperglycemia among hospitalized patients, to be used together with glucose monitoring at 1-2 hr intervals, as a decision-support tool that will assign a future intravenous insulin infusion rate (IR) after each iteration of the algorithm.The algorithm must safely and effectively achieve different target ranges of blood glucose appropriate to various patient populations, while minimizing risk for hypoglycemia.A sigmoidal relationship between IR and BG may help achieve these goals.A prototype for the present work, developed as a set of three tabular algorithms, each having a different target range but created from a single mathematical design, has been implemented in clinical practice [23][24][25].The specific objective of this paper is to address an essential step necessary to computerization.For a patient at a particular time during a treatment course whose hourly maintenance insulin requirement (or the maintenance rate, MR) is assumed to be known or determinable and whose IR for the next interval of time during treatment needs to be assigned, we will describe an equation, having as its domain the spectrum of possible BG values, from which a single point-of-care measurement is applicable each time the equation is used, giving as output a value for IR.For the sake of brevity, in this paper, "BG" will be used to refer to a

MR
maintenance rate, insulin infusion rate needed to maintain stable value of blood glucose under ambient conditions of medical stress and carbohydrate exposure, comprised of basal insulin requirements and, for those receiving continuous carbohydrate exposure, also providing nutritional insulin coverage; in order to achieve target range control, IR equals MR when BG value equals a mid-target range blood glucose (BG target ), and IR is close to MR between BG target_minus and BG target_plus .
RR restraining rate of insulin infusion, a positive-valued ascending sigmoidal function of BG that approaches MR asymptotically, having values between zero and MR; when expressed as a function of BG, RR exhibits greatest increase with respect to BG during ascent of BG to BG target_minus and negligible further increase as it asymptotically approaches MR for BG values above BG target_minus .RR is so-named because it restrains the ascent of BG to occur at less than the maximum rate of ascent that otherwise would be seen under ambient conditions of medical stress and carbohydrate exposure.For BG on the target range, the RR provides maintenance insulin coverage (RR ≅ MR).
CR correction rate of insulin infusion, a positive-valued ascending sigmoidal function of BG that approaches the maximum correction rate CR max asymptotically, having values between zero and CR max ; when expressed as a function of BG, CR exhibits negligible increase with respect to BG until BG reaches BG target_plus and greatest increase with respect to BG during correction of hyperglycemia such that BG is greater than BG target_plus .CR is so-named because, when added to RR, it results in a decline of BG from an elevated value toward target.
IR insulin infusion rate, a function of blood glucose (BG) and MR; IR should equal MR at mid-target blood glucose.
the value of MR is identified.Then the MR is used in the equations yielding IR (eqn.1,Table 1).

The Maintenance Rate (MR)
Over the treatment course of a single patient, BG will change and IR will be reassigned repeatedly.After initialization, every iteration of the algorithm requires discovery of a value for an intermediary variable, MR.MR is expected to change over time, in relation to carbohydrate exposure and conditions of illness or diurnal variation that might affect insulin resistance.MR is identified from antecedent interactions by the rate of change of blood glucose, previous IR, and potentially other inputs.MR, once assigned its value for the coming iteration, is used for that iteration of the algorithm as a fixed parameter in the equation that will give IR as a function of BG.Therefore, the IR here may be referred to as "MR-Dependent IR." The rate of insulin infusion necessary to prevent unchecked glucose production and release to the bloodstream by the liver is the basal rate (BR).The rate of insulin infusion necessary to promote cellular uptake and retention of administered carbohydrate is the nutritional rate (NR).For patients receiving continuous carbohydrate exposure, MR is the sum of BR and NR.MR differs among patients and patient populations, and changes in the course of treatment of a given patient.Determinants of MR include insulin resistance and any continuous carbohydrate exposure.An algorithm utilizing a computational method of estimating MR has been evaluated "in silico" [27].

Insulin Infusion Rate (IR) under Progenitor Algorithms
For a given MR, in a published method of assigning IR which we will call the "three segment method," the IR function utilized a different rule on each of three BG domains [23][24][25][26].Additionally, there were rules for response to hypoglycemia with BG < 70 mg/dL.The three segments of the IR function were demarcated by two possible values for the variable BG, at which the function was not differentiable: the BG values marking the lower and upper limits of the target range.In the present work, the terms BG target_minus and BG target_plus replace terminology of earlier publications, in which the limits of the target range were designated as BG true_target and BG upper_target [25,26].

•
For BG < BG target_minus , the fractional completeness of ascent of BG (FCABG) between BG = 70 mg/dL and BG target_minus was defined as follows: (2) For BG < BG target_minus , the IR was an ascending exponential function of FCABG with the following characteristics: (a) lower asymptote equal to zero, (b) a negligible rate of insulin delivery (0.1 units/hr for adults) at BG = 70 mg/dL, and (c) value approaching MR at BG target_minus .The curve was characterized by upward concavity.

•
For BG target_plus ≤ BG, the fractional elevation of BG (FEBG) above BG target_plus compared to the elevation of BG critical_high above BG was defined as follows: For BG target_plus ≤ BG, IR was an ascending bounded decaying exponential function of FEBG with the following characteristics: (a) value equal to MR at BG target_plus , (b) value equal to MR plus a correction rate (CR) of insulin infusion, designed to induce a specified rate of descent (ROD) of blood glucose during the next iteration, at a BG value called BG critical_high , and (c) value approaching an upper asymptote that equaled MR plus an MR-dependent algorithm parameter, defined as the maximum correction rate, CR max.The curve was characterized by upward convexity.

•
For BG target_minus BG < BG target_plus , IR had a constant value equal to MR; FCABG equaled 1, and FEBG equaled zero.The curve was a straight horizontal line.
Additionally, assumptions were made about the rate of ascent and ROD of blood glucose, which might be testable in actual clinical trials.The spontaneous rate of ascent of blood glucose, in the absence of hypoglycemia or a counter-regulatory response to hypoglycemia, was assumed to have a maximum value at BG ഡ 70 mg/dL, occurring during negligible insulin infusion at the boundary between euglycemia and hypoglycemia.Spontaneous ascent was assumed to be completely restrained if IR = MR.ROD of BG was assumed to be predicted by the amount by which the hourly IR exceeded MR, called the "correction rate" or CR [26].The values of CR under this earlier definition were determined by analogy to an ambulatory rule, variably termed the rule of 1500, the rule of 1800, or the rule of 1700, predicting an inverse relationship between the total daily dose of insulin and the decline of BG that might result from administration of 1 unit of insulin [28].Using a parameter, G-per-DIEM (glucose per daily dose of insulin exogenously mediated), the rule states that: (decline of glucose)/(1 unit insulin) equals (G-per-DIEM)/(total daily dose of insulin).
Assuming CR was added to the true MR, then under both the old definition of CR and the present definition of the CR function, the targeted rate of descent during hyperglycemia for each coming iteration (ROD ideal,next ) would bear a direct relationship to CR, given by: (4) The mathematical design, intended for eventual computerization, first was used to construct tabular paper algorithms having three distinct target ranges appropriate to different medical conditions: HHS, DKA, and other critical care conditions treated with a standard default algorithm.Within the paper tabular algorithms, each column was identified with a value for MR.The caregiver had rules for switching columns based on rate of change of BG.Each row represented a range of BG values.A BG value within the range represented by each row was used to compute an IR for each cell of the table.Within each column, going from rows showing the lowest to highest BGs, the IR could be represented as a step function having roughly sigmoidal shape, giving IR as a function of BG, with inflection point roughly at the target range of BG [23][24][25].

Insulin Infusion Rate as a Quadruply-Asymmetric Sigmoid
Algorithm parameters and intermediary variables, previously described for use in a method of assignment of IR having discontinuities of the first derivative and piloted as paper algorithms, were adapted for use in the present algorithm [25,26].The preconceived shape of a graph of IR as a function of BG ideally would be a strictly ascending asymmetric double sigmoid, with a middle that was relatively flat.The relatively flat middle of the curve would correspond to a brief segment of the domain of BG regarded by clinicians as a target range of BG, at which a maintenance rate of insulin should be infused.
To achieve the desired shape, the function IR would be constructed as the sum of two asymmetrical sigmoidal curves, named the restraining rate (RR) and the correction rate (CR), with their inflection points slightly offset from each other.Adaptations of the Richards differential equation would be used from which appropriate solutions would be derived, which would be the equations for RR and CR [29].In the differential equations, the derivatives of the RR or CR with respect to BG were equated with a function of RR or CR, respectively (eqns.18-19, Table 5).
Each differential equation contains 2 undefined constants and would give rise to a 3 rd undefined constant of integration.Reasonable values for the constants were suggested by clinical experience in using intravenous insulin infusion to treat hyperglycemic crisis, or hyperglycemia which may accompany critical illness, and by favorable experience with progenitor paper algorithms [23][24][25].For each equation, after rearrangement the unknown constants finally will drop out, to yield three defined constants that depend upon parameters of the algorithm.With the use of the three constants, in the final simplified equations for RR and CR, respectively, the IR as a function of FCABG and FEBG will appear as a parametric function of BG, such that IR equals the sum of RR and CR (eqn. 1, Table 1).For purposes of graphing RR, CR, and IR, Mathematica ® (version 9.0) was used [30].

RESULTS
The intended characteristics of the functions for MR, RR, CR, and IR are described in Table 1 (a Glossary of Terms), and the parameters and intermediary variables of the equations are defined in Table 2.The asymptotes and defining points for the logistic equations are given in Table 3.In Table 4, solutions to the differential equations are summarized, incorporating the parameters that will yield the desired functions, with the derivation given in Table 5.For the parameters described in Table 2, population-specific values are provided in Table 6.
The two columns of Table 5 show points of parallelism between the derivation of the equation for RR and for CR, each having a known inflection point and another known point that was assigned by design.
The asymmetry factors "n" and "m" for the functions RR and CR, respectively, were derived by setting the second derivatives with respect to BG equal to zero, at those values of RR and CR that we judged clinically appropriate to be inflection points, relative to the lower asymptote, zero in each case, and the upper asymptotes MR and CR max , for RR and CR, respectively.It was judged clinically appropriate that RR should terminate its ascent sharply and begin its approach toward its upper asymptote (MR) for BG > BG target_minus , and that CR should terminate its ascent most sharply at BG critical_high .Therefore, the inflection points of the functions were chosen to occur at  BG critical_high -C R @BGcritical_high Artificial BG and upper asymptote.

Table 2. Algorithm parameters and intermediary variables (Continued)
CR @BGcritical_high Value assigned to CR at its inflection point. ( CR max (9) IR max *Maximum permitted IR for algorithm, overriding higher rates that might be computed.* A patient population parameter, see Table 6.
BG target_minus for RR and at BG critical_high for CR, where it had been determined clinically that the values of RR and CR, respectively, should be MR minus and CR @BGcritical_high .When the second derivatives of the differential equations were set to zero and when MR minus and CR @BGcritical_high were entered into the equations as the desired values for RR and CR at these inflection points, equalities were obtained giving MR/MR minus and CR max /CR @BGcritical_high as functions of the asymmetry factors "n" and "m", respectively.Determination of the values of "n" and "m" as the inverse functions, respectively of MR/MR minus and CR max /CR @BGcritical_high , is not given by any explicit expression and requires computational methods.
After transposition of the left-hand denominator dBG and the right-hand side of the differential equations for dRR/dBG and for dCR/dBG, and after integration of the two sides of each transposed equation, values for the constants of integration were provided as known points, assigned through clinical judgment.For RR, it was judged that BG ≤ 70 mg/dL should yield an IR so low as to be negligible, specifically RR @BG70 .For CR, it was judged that BG target_plus should be the highest BG value at which CR would be negligible, close in value to RR @BG70.CR should rise from its lower asymptote (zero)

RR CR
Asymmetry constants "n" and "m" (Inverse function to be given by a computational method) (10) (11) Sum of logistic equations yielding strictly increasing asymmetric double sigmoid and begin to assume values indicating clinically effective rates of insulin infusion for BG ≥ BG target_plus .The points (70 mg/dL, RR @BG70 ) and (BG target_plus, ∆ MR ) therefore were assigned as known points on the curves RR and CR respectively.
There was a remaining unknown in the equations, specifically the rate constant.By providing another known point for each equation, namely (BG target_minus , MR minus ) for RR, now incorporating BG target_minus as the BG value corresponding to RR at its previously identified inflection point, and by identifying (BG critical_high , CR @BGcritical_high ) for CR, now incorporating BG critical_high as the BG value corresponding to CR at its previously identified inflection point, it was possible to solve the equations giving a value for the rate constants.
RR and CR are strictly increasing functions having positive values at all BG values.The value of RR at BG target_minus is (MR -∆ MR ).The value of CR at BG target_plus is ∆ MR .Based on eqn. 1, Table 1, it is evident that there exists a value of BG at which IR equals MR, that may be called "BG target ", such that: BG target_minus < BG target < BG target_plus A Quadruply-Asymmetric Sigmoid to Describe the Insulin-Glucose Relationship during Intravenous Insulin Infusion

Second derivatives and designer-defined inflection points
By design, RR = MR minus at inflection point; and By design, CR = CR @BGcritical_high at inflection point; and must equal 0 at inflection point.must equal 0 at inflection point.
Asymmetry constants "n" and "m", functions defined for values of real argument > 1

"n" is a function dependent upon "m" is a function dependent upon
There is no explicit formula for the asymmetry constant, even with domain restricted to values > 1.
The values of "n" and "m" at given values of the argument may be estimated using Newton's method.Shorthand within RR; let S be defined as: Shorthand within CR; let T be defined as:

Integration of eqn. 18 for RR and eqn. 19 for CR
The rate constants are defined by using the second designer-defined point of RR and CR Solve eqn.22 for "k" and use designer Solve eqn.23 for "j" and use designer -defined point to assign its value -defined point to assign its value Replacing the variables in eqn.26 with Replacing the variables in eqn.27 with the known designer-defined the known designer-defined point = (BG target_minus , MR minus ), which point = (BG critical_high , CR @BGcritical_high ), which is also an inflection point, we learn is also an inflection point, we learn the value of "j": the value of "k": Shorthand within RR; let P be defined as: Shorthand within CR; let Q be defined as:   See also Table 2 for definitions of parameters.
An individual treatment course under the adult algorithms will be discontinued if re-estimation of MR during therapy repeatedly yields MR< 0.5 unit/hr.The algorithm will assign the IR to be the lesser of • IR max • IR as computed by the equation for MR-Dependent IR.
* For patients judged a priori to have insulin sensitivity or insulin resistance, MR initial for DKA might be 2 or 4 units/hr respectively, rather than the default 3 units/hr.No upward revision of MR to a computed value greater than MR initial will be made until rehydration time has elapsed.
(The population value that has been designated BG lower_target in previous publications is renamed here as BG lowest_acceptable , the value designated as BG true_target in the previous publications is renamed here as BG target_minus, and the value designated BG upper_target is renamed here as BG target_plus [24][25][26]).
It is hypothesized without proof that there exists no explicit formula for BG target , which generally will have an irrational value and cannot be regarded as a realistically exact therapeutic true target, despite the fact that it is identified mathematically with MR as the point (BG target , MR).
In Figure 1, the curvature of one limb of RR is too subtle to be discerned in the graph, but actually RR and CR both are sigmoidal.The tight curvature about the inflection point of RR is shown in a magnified view (Figure 2).The curve for IR has two segments that are upwardly convex and two that are downwardly convex and thus has 3 inflection points, one each at the boundaries of the target range, and one between those boundaries.Rotation by 180 degrees about each of the 3 inflection points does not uncover symmetry between any pair of segments having opposite convexity, when the four possible comparisons are made.Therefore, the function for IR will be designated as quadruply asymmetric.
For any given set of parameters defining glycemic targets, different possible values of MR will generate a family of functions for IR which, when graphed, will appear as a family of strictly increasing asymmetric sigmoidal curves, with each member of the family identified by its MR (Figure 3).

Parameters for Pilot Studies
The G-per-DIEM parameter was chosen in our case to be 1800 mg/dL [26].The RR @BG70 (0.1 units/hr under the adult algorithms), considered to be negligible for most adults, may be at the lower limits achievable by some adult infusion pumps.Use of lower rates may require use of neonatal or pediatric equipment.Tentatively, in order to ensure that MR is at least 5-fold greater than RR @BG70 , it is suggested that use of the algorithm for adults should be discontinued if MR is estimated < 0.5 unit/hr.A long term strategy for some patients having low MR and low body weight could be to convert to use of the weight-based pediatric algorithm, which would assign potentially lower values for RR @BG70 and would recommend smaller incremental adjustments of MR and IR during therapy.
Under the adult algorithms, if MR = 0.5, 1.0, or 5.0 units/hr, then at BG target_plus under eqn.4, the hypothetically predicted ROD of BG attributable to ∆ MR = CR at BG target_plus would be 12.00, 6.75, or 1.47 mg/dl per hr, respectively.As shown in Table 6, specific values for ROD @BGcritical_high and for ROD ideal,next,max were assigned by clinical judgment as algorithm parameters giving targeted ROD for higher blood glucose.

Constraints upon Parameters and Intermediary Variables
Users should be able to establish preferences for algorithm parameters that might differ from those chosen at the pilot sites.Within clinically reasonable guidelines, for institutional programming, flexibility should be assumed to be achievable.A reassignment of target range would require not only new values for BG target_minus and Restraining rate, correction rate, and infusion rate as functions of blood glucose.Vertical dashed lines demarcate the boundaries of the target range using parameters of the standard default algorithm (Table 6), shown here for MR = 3 units/hr.Restraining rate as a function of blood glucose (BG) in magnified scale.The tight curvature is shown at BG values just above the inflection point for RR, which is at BG 120 mg/dL for the standard default algorithm, shown here for MR = 3 units/hr.
BG target_plus , but also potentially new values for BG critical_high , for ROD @BGcritical_high and for ROD ideal,next,max , in order to ensure that after 1-2 hr use of the CR during hyperglycemia, the targeted rate of descent (ROD ideal,next ) would not be predicted to deliver a blood glucose value < BG target_plus .For patients having chronic kidney disease, the parameter ROD ideal,next,max might be decreased without necessarily making any adjustment to target range BG values.Pediatric values for parameters would best be established by practitioners familiar with the treatment of infants and children.Additionally, the algorithm mathematically requires constraints upon parameters.Most obviously, it is necessary that 70 mg/dL < BG lowest_acceptable < BG target_minus < BG target_plus < BG critical_high .As revisions of MR are made during treatment, it is required that the intended relationship should be preserved for all values of MR, that ∆ MR < RR @BG70 << MR minus < MR.The algorithm cannot be used if computation of MR leads to contradictions of the relationships above.
The value of ∆ MR , equaling the difference between MR and MR minus , was intended to have a value close to RR @BG70 .However, it could be found that some insulinsensitive patients may have an MR of the same order of magnitude as RR @BG70 .The formula for ∆ MR (eqn.6, Table 2) has been constructed to ensure that MR minus (i.e., MR -∆ MR ) will exceed RR @BG70 (Appendix 1).
We further require that MR minus should be substantially greater than RR @BG70 .As revisions of MR are made during treatment, it also is required that ∆ MR (the value of CR at BG target_plus ) << CR @BGcritical_high < CR max .By design, for BG < BG target_minus , the value of the CR function should be negligible.Additionally, for BG target_plus < BG, the CR function giving the correction rate above RR should closely approximate the correction rate above MR, given by solving eqn. 4 for CR, which is equivalent to requiring that Relationship during Intravenous Insulin Infusion  MR-RR be negligible for BG target_plus < BG.These goals are achieved by (a) requiring that RR @BG70 << MR, (b) defining ∆ MR to have a value nearly equal to RR @BG70 , and (c) using ∆ MR in the definitions of RR at BG target_minus and CR at BG target_plus .
The designer or user does not directly assign values for ∆ MR , CR @BGcritical_high or CR max .Nevertheless, a user seeking to revise parameters must preserve the relationships above.Any permissible MR must substantially exceed RR @BG70 .The predicted ROD ideal,next attributable to CR at BG target_plus must be exceeded by values selected for ROD @BGcritcal_high and ROD ideal,next,max , which will determine the values for the CR @BGcritical_high and CR max .
To prevent appearance of undefined quantities in the equations for RR and CR when fractional roots are taken of the expressions [1 + S × (S × P) -FCABG ] and [1 + T × (T × Q) FEBG ], a sufficient condition is that P and Q should have positive values for all BG and all choices of parameters.The requirement will be met if "n" and "m" both have positive values, a condition which is assured if the parameter constraints above are observed and if the value of MR/MR minus and CR max /CR @BGcritical_high each is less than the base e of the natural logarithm.

DISCUSSION
Efforts at strict glycemic control using conventional insulin infusion algorithms have been thwarted by oscillations of blood glucose, such that retrenchment from a philosophy of tight control for hospitalized patients has been recommended by professional societies.Nevertheless, reasonable control still is recommended, with targets differing according to patient population, and with safeguards against hypoglycemia [5][6][7][8][9][10].
A modular design presently is in trials in the ambulatory setting, such that the basal requirements of an artificial pancreas might be separately managed from bolus therapy or interception of threatened hypoglycemia [31,32].A successful model for assigning IR in the hospital for patients whose maintenance requirements have been estimated might be viewed as a model for assigning corrective adjustments in the ambulatory setting for patients whose basal rates have been established.The linear rule, used in ambulatory medicine for determination of the correction dose component of bolus therapy, perhaps could be improved during continuous subcutaneous insulin infusion if it were replaced by a sigmoidal rule for correction dose assignment.

Algorithm Description and Comparison to 3-Segment Method
The present method refines the previously published concept of MR-dependent IR [23][24][25][26] so as to yield IR as a strictly increasing, continuous and differentiable function on the entire domain of BG, consisting of the sum of the values of RR and CR at any value of BG (Eqn.1).The FCABG (eqn.2) and FEBG (eqn.3) are intermediary variables depending on BG that appear in the expressions for RR and CR, respectively.There will be a BG value (BG target ), generally an irrational number, not required to be known or used in any of the equations, such that BG target_minus < BG target < BG target_plus , and at which IR = MR.The BG target is defined as the BG that would be yielded by the inverse of the function for IR, at the value IR = MR.A graph of IR as a function of BG will show a strictly increasing but nearly flattened mid-region, for BG target_minus ≤ BG < BG target_plus , displaying a reduction of the rate of increase of IR as a function of BG, when BG is in the neighborhood of BG target .A graph of IR will display no symmetries about its three inflection points and will have roughly the form of an asymmetric doubly sigmoidal curve.The graphs of the functions RR and IR appear sharply bent about the inflection point at BG target_minus , but in fact they are curved (Figure 2).Depending on the ratio of CR max to CR @BGcritical_high , the curve for CR and IR also may be fairly sharply bent about BG critical_high , as seen in the algorithms for HHS and DKA (Figure 3).
The CR function in the present work, equaling the value of (IR minus RR), is not identical to CR in the 3-segment method but closely approximates the value of (IR minus MR).The value of (IR minus MR) is the hourly insulin rate used for correction, hypothesized to predict ROD ideal,next according to eqn. 4 under either method.
In contrast to the previously described three-segment method, under the method of assigning MR-dependent IR described here, yielding a smooth function for IR, the equations for FCABG, FEBG, and IR do not differ according to whether BG is below, at, or above a value demarcating the boundaries of the target range.In case values for BG target_minus , BG target_plus , and BG critical_high are redefined, in order to apply the equations that are unique to each segment under the 3-segment method, there is a need to demarcate segments of the BG domain with correctly-matched choices of inequality symbols ( "<" vs. "≤" ) for the FCABG and FEBG, and for the definition of the dependent variable itself, the IR.Under the asymmetric sigmoid method, one set of equations applies to the entire domain of BG values, so that replacement of BG parameters according to institutional choice of target ranges is less likely to introduce programming error.

BG as the Independent Variable, Itself not a Function of Time
In the biological sciences, for a variable undergoing growth, the Richards differential equation commonly shows a derivative of that variable with respect to time, equated to a function of the variable itself.As a growth equation, the integrated function is designed to show the effect that at low values of the variable there is rapid growth, and at high values the growth rate declines.The solutions to the Richards differential equation showing growth over time are expressed with time as the independent variable.Here, in the adaptations of the Richards equation that we are using, the variable undergoing "growth" is the IR.The independent variable in the solutions is BG, which is not a function of time.To prevent confusion, it has to be mentioned that the patient's blood glucose is continuously changing over time and could itself be represented in other contexts as a function of time.The equations for IR here are reused repeatedly at each timepoint when a change of IR is considered.Nevertheless, at each utilization of the equations for IR, the meaning of BG as an independent variable in the solutions to the Richards equations will refer to the spectrum of possible values that could be detected by measurement of the current BG at that single timepoint.The second and first derivatives of IR will be identified with respect to BG, not time.The value of IR as a function will be defined by the current BG.

Meaning of the Equations
Many clinicians acknowledge the desirability of recognizing physiologically and therapeutically that there should be not a single target, but a target range of BG.Stimulation of counter-regulation marks what might be considered the lower end of the

42
A Quadruply-Asymmetric Sigmoid to Describe the Insulin-Glucose Relationship during Intravenous Insulin Infusion physiologic range in normal physiology (BG of 70-72 mg/dL) [11,12].In clinical therapeutics, for safety, a target range may be established that is broader and higher than the normal physiologic range.A doubly sigmoidal curve was favored to describe the insulin-glucose relation under our algorithm, having capability of programming a relatively flat mid-region or target range between BG target_minus and BG target_plus , over which the rate of insulin infusion does not sharply increase with BG.
The method of assigning MR-dependent IR described here requires a straightforward summation of two equations, RR and CR (Figure 1, bottom panel).Each equation, designed to replicate a physiologic pattern of time-averaged insulin release by the beta cell over the short term, contributes to a model of therapy for patients who are not eating.Simply stated, beta cells provide some insulin coverage when glucose levels are below target.The rate of ascent of blood glucose to target depends upon hepatic glucogenesis, the rate of infusion of carbohydrate (if any), and the peripheral and hepatic disposal of glucose, which are influenced in turn by several biological regulators.The rate of ascent of BG responds to the restraining action of insulin.During ascent of glucose from levels below target, the model for insulin response is reflected in the curve for RR, the insulin infusion rate used to restrain the ascent of glucose to occur at less than the maximum rate of ascent that otherwise would be seen, and which approaches the full maintenance rate MR as the blood glucose approaches target (Figure 1, top panel).When glucose levels are above target, beta cells provide maintenance plus correction insulin.The ROD of BG, from an elevated level to the target range, depends upon hepatic and peripheral glucose disposal in response to correction insulin, added to the maintenance insulin requirement.The model for insulin response to hyperglycemia is reflected in the curve for CR, the insulin infusion rate added to the RR, in nonlinear proportion to the elevationof BG above the target (Figure 1, middle panel).When the glucose is close to target, the incremental contribution of CR is negligible, such that IR approximates MR, the infusion rate that preserves BG close to its current value.

Paper Predecessor
A predecessor of a sigmoidal design has been studied using a step function having roughly a sigmoidal shape, fashioned according to the "three-step" method [25].Preliminary clinical results have been promising, supporting future development of the method proposed in this work.The "three-segment method" for computing IR was used to construct the tabular paper protocols of the Saint Francis Hospital in Evanston, IL [23][24][25].A set of paper tabular column-based algorithms, each having a different target range, has been studied in actual practice for 58 patients [25].Each of the 3 paper algorithms had a distinct target range: 200-299 mg/dL for HHS, 150-199 for DKA, and 100-149 for general ICU care.The results of use during a predefined data collection interval showed that after the first attainment of the target range BG, the desired targets were maintained without severe hypoglycemia and with clinically acceptable glycemic variability.By modifying the shape of the curve for CR in comparison to the bounded decaying rule of the predecessor algorithm, now having upward concavity between BG target_plus and BG critcal_high , we propose in the present work to improve upon an already-satisfactory design for assignment of IR by reducing risk for excessive dosing as target is approached from BG above target.

Limitations and Areas for Future Technical Work
Parameters for use in initial piloting of the present algorithm have been identified for the standard default adult, DKA, and HHS algorithms (Table 6).Appropriateness of these parameters is supported by clinical use of paper algorithms that were structured according to the closely-related 3-segment method of IR assignment.For the pediatric parameters shown in Table 6, and for a possible future option permitting real-time parameter reassignment by users, in silico testing will be required in advance of human trials to prove the computability of the asymmetry parameters "n" and "m" and variables that are dependent upon them, under the range of expected clinical scenarios.
In the context of finding an asymmetry constant to describe an asymmetrical curve that would fit observational botanical data, the claim has been made that use of the asymmetry constant yields results that are "unstable" [33,34].If observational data are used to infer the asymmetry constants, sampling and imprecision of measurement could lead to inflation of errors.Ideally, in the equations, the same source of blood should be used for every BG [35].In the present paper, it is proposed to move from predefined clinical parameters, to infer the value of the asymmetry parameter.The "n" and "m" asymmetry factors require only one input each: the ratios MR/MR minus or CR max /CR @BGcritical_high , respectively.The values assigned to MR minus , CR max , and CR @BGcritical_high are dependent upon reassignment of MR after each iteration of the algorithm.Measurement errors could affect the value assigned to MR.We will have a predefined value for RR @BG70 , clinical limits on increments that can be made to MR, and a permissible upper limit of value that IR or MR can assume, thus limiting the value that MR and the MR-dependent CR may assume.The greatest clinical precaution is that the program should choose the more conservative among possible estimates for any upward readjustment of MR.
The ideal method of implementation might be that a hand-held device would be dedicated to each patient during the treatment course, used by the nurse.Such devices could be in communication with a calculator capable of performing the computations.
Calculation of the asymmetry constants in real time requires software that will correctly handle a floating point decimal.Alternatively, for each target range, calculation could be performed in advance of implementation, and the function could be expressed as a set of ordered pairs on a countable domain of uniquely-identified duples of MR and BG values.Advance computation of a table for each target range is the approach we will favor for initial piloting.For example, for use of the standard default algorithm, the set of duples could consist of all combinations of 236 MR values between 0.5 to 24.0 units per hour (given in increments of 0.1 unit per hour, such that RR @BG70 < MR ≤ IR max ) and 1930 BG values between 70 to 1999 mg/dL (given as integers).Each duple would be matched to a value of IR, thus defining the function for all practical purposes with a set of 455,480 ordered pairs.
It is recognized that some algorithms for insulin infusion capture the pulsatile nature of insulin secretion in their design [36,37].The present algorithm does not replicate physiologic pulsatility [38,39].Rather, the algorithm attempts to capture changes in average delivery by making adjustments of IR at 1−2 hr intervals, in effect creating an ultradian pattern that is slowly oscillating about the moving value of MR.
It is recognized that complex algorithms may announce the optimal timing for future blood glucose sampling, using a strategy that predicts future disturbances and considers the potential error of predictions [40][41][42].The algorithm described in the present work would be enhanced by evidence about rate of ascent or ROD of BG in relation to IR, so as to optimize predictions about timing of future sampling.Our rules for upward revision of MR, presently constrained conservatively, potentially could be modified based on future evidence about error of predictions.

Rationale in Relation to Physiology of the Glucoregulatory System
In biologic systems that are characterized by feedback inhibition, an approximately sigmoidal relationship may be seen between levels of a regulator and a measure of its biologic effect [43].The insulin-blood glucose regulatory system is characterized by negative feedback between the regulator (insulin) and its effect upon the regulated substance (reduction of blood glucose), as discussed under 4.3.Peripheral insulin receptors exhibit saturation behavior.Since the action of insulin may persist after its administration has ceased, excessive dosing should be avoided.Treatment achieving a sigmoidal relationship between BG and IR would help prevent excessive insulin dosing at both lower and higher BG values, thus theoretically reducing the risk for hypoglycemia.
Commercialized medical devices in the USA have been developed using a multiplier design, and similar algorithms have been programmed by institutions for local implementation.Multiplier algorithms seek to discover the insulin sensitivity of the patient and make the IR dependent upon the sensitivity, but differ from the algorithm proposed here in the method of attaining specific targets.Under multiplier algorithms, a common design is to assign a multiplier such that the IR = (multiplier) × (BG -60 mg/dL).Reassignment of the value of the multiplier is governed by response to therapy [44].Under a multiplier design, at different levels of insulin sensitivity, the algorithm thus assigns IR as a linear function of BG [45].Attainment of a specific glycemic target may require interruption of the infusion for readings below target, requiring a waiting time interval before resumption of insulin infusion.In contrast, when managing glycemic excursions below target, the sigmoidal design proposed here may obviate the necessity for interruption of insulin infusion.For patients having type 1 diabetes, the consequence of prolonged interruption could be development or exacerbation of ketosis.

Inputs Determining and Modifying Assignment of MR
After initialization, the value of MR is inferred by rules hypothesizing a relationship between MR, the observed rate of change of blood glucose [referred to as (BG current -BG previous )/(∆time) in previously published descriptions], the previous IR (formerly referred to as IR previous ), and the parameter G-per-DIEM [24][25][26]28].Although the computation or assignment of value of MR requires knowledge of rate of change of blood glucose and therefore incorporates its current value (i.e., BG, as used in this paper), nevertheless, it is the rate of change of blood glucose with respect to time, and the current BG status in relation to target (above, at, or below target), that determine the value of MR, not the specific value of BG.Under the method, the MR is reassigned repeatedly in the course of treatment, usually every 1-2 hr.
If the previous blood glucose has exceeded the upper target, the MR may be computed tentatively based on previous ROD of BG and previous IR.However, under several conditions, rules of the algorithm would prevent computation or would force a conservative reduction of the computed value of MR.These situations might include: first iteration, rehydration time not elapsed, BG between lowest acceptable value (BG lowest_acceptable ) and BG target_minus , BG below BG lowest_acceptable , hypoglycemia, status post treatment for hypoglycemia, interim interruption of insulin infusion, and other nursing inputs.(The population value that has been designated BG lower_target in previous publications is renamed here as BG lowest_acceptable , being the level of BG below which a forced reduction of MR is indicated [24][25][26].)Nursing inputs predicting risk for hypoglycemia might include planned reduction or interruption of intravenous dextrose infusions, interruption of enteral feedings, interruption of glucocorticoid therapy, increase of insulin as additive to total parenteral nutrition, onset of renal failure, or other perceived risk factors for hypoglycemia, as permitted by institutional programming design.An alert from nursing staff would result in reduction of MR from the computed value.
The estimated value of MR must exceed RR @BG70 by at least 5-fold or else treatment should be revised (see Section 3.1).A maximum value permitted for MR is determined by the maximum allowed IR, the algorithm population parameter IR max (Table 6).
It is recognized that estimates of MR are limited by the accuracy of measurements of point-of-care glucose that are used to determine the rate of change of BG and the correct recording of test times, and that the true rate of change may not be captured by single measurements made at 1-2 hr intervals.Therefore, in order to reduce the risk of overestimation of MR, we plan to use a damping rulethat may retard upward revisions of MR.We will aim to allow increasing values of MR to be approached over several iterations, if the need for a higher MR value is computed repeatedly.

Response to Hypoglycemia
We believe both algorithm design and clinical precepts may prevent or reduce the occurrence of hypoglycemia (BG < 70 mg/dL).The algorithm should require that nursing staff contact the provider in case hypoglycemia does occur.In any case of hypoglycemia, there is a need for evaluation of whether an unexpected change of carbohydrate exposure has occurred, whether insulin has been delivered by another route (such as an additive to total parenteral nutrition), or whether the medical condition of the patient has changed.Realistically, the value of RR @ BG70 is at the lower limits of the capabilities of infusion pumps used in clinical practice and is physiologically negligible.After hypoglycemia, the MR is adjusted according to the rules of the algorithm.The hypoglycemia is treated with concentrated dextrose.For BG < 70 mg/dL, the following hypoglycemia rule replaces computation of IR: regardless of immediate post-treatment BG measurements that might be higher than 70 mg/dL, these immediate post-treatment values are not used for reassignment of IR; rather, the IR for the next iteration of the algorithm is reassigned as the physiologically negligible rate, RR @BG70 .This conservative rule overrides other rules for IR assignment that might have delivered a higher rate of insulin infusion in response to a higher BG measured shortly after treatment of hypoglycemia.For insulin-dependent patients, recurrences of hypoglycemia could necessitate use of an algorithm having a higher target range.

46
A Quadruply-Asymmetric Sigmoid to Describe the Insulin-Glucose Relationship during Intravenous Insulin Infusion An alternative for responding to hypoglycemia would be to interrupt the infusion of insulin.However, timing of resumption would have to be studied.The argument against turning off the insulin infusion is that in a busy intensive care unit, outside of a research context, there is risk of failure by personnel to remember to resume the infusion.We argue that interruption of insulin infusion therapy should require a new order by the provider specifying discontinuation of the insulin infusion protocol.

The Example of Diabetic Ketoacidosis
As an example, under the DKA algorithm, a patient treated for DKA might be assigned to have the default algorithm parameters including an assumed initial MR (MR initial ) of 3 units/hr.There would be corresponding highest computable IR ≅ 9 units/hr and, for an initial critical high BG = 350 mg/dL, an initial IR = IR @BGcritical_high ≅ 8.0 units/hr.For an assumed MR initial of 4 units/hr, as might be chosen for a patient judged a priori to have more severe insulin resistance, the corresponding rates for highest computable IR and IR @critical_high would be ≅ 12.0 and 10.7 units/hr, respectively.For an assumed MR initial of 2 units/hr, as might be chosen for a patient judged a priori to be especially insulin sensitive, the corresponding rates for highest computable IR and IR @critical_high would be ≅ 6.0 and 5.3 units/hr, respectively.
As the blood glucose declines, there also is a decline of the targeted ideal rate of descent of blood glucose for the next iteration of the algorithm (ROD ideal,next ), given by eqn.4, until the value for CR approaches a negligible value, between zero and ∆ MR , upon attainment of target range control, at which point IR ≅ MR.Therefore, even if MR is unchanged, the IR declines as BG approaches target.
It is intended that the MR be revised during therapy.According to response, if required under algorithm rules, at any time during treatment, the MR may be adjusted downward, and after 4 hr of rehydration, MR may be adjusted upward.Upon attainment of target range control, the MR will cover both basal needs under conditions of medical stress, and also the needs imposed by delivery of dextrose containing fluids.As a result, at a given BG, according to revisions of MR, the IR varies among patients.
By institutional preference, instead of using a default MR upon initiation of treatment, the assignment of MR initial could be made among adults according to body weight, as is recommended among children.For a 70 kg man, an assumption of MR initial = 3 units/hr would translate to MR = 0.043 units/kg per hr, as might be preferred as a default assumption.An MR initial of 4 units/hr for a 70 kg man would translate to 0.057 units/kg per hr, which might be preferred for patients suspected a priori to be insulin resistant.An MR initial of 2 units/hr for a 70 kg man would translate to 0.029 units/kg per hr, which might be preferred for patients suspected a priori to be insulin sensitive.
Limitations of the method include the fact that the severity of DKA or requirement for insulin is not determined solely by the presenting blood glucose or body weight.For example, under the proposed algorithm, very insulin-sensitive patients could receive initially unnecessarily higher IRs than required.On the other hand, under the algorithm, unless a modification is introduced for management of relatively euglycemic diabetic ketoacidosis, patients having DKA with relatively euglycemic values of presenting glucose and requiring insulin could nevertheless receive lower than optimal initial IRs.Although "euglycemic " may be a misnomer, glucose values < 200 or < 250 mg/dL together with appropriate history might be used to define a case of DKA as relatively euglycemic, in comparison to the more typical presentation of patients having greater severity of hyperglycemia at the time of presentation with DKA [8,46].It has been observed that in relatively euglycemic DKA, increased glucose administration may be required to facilitate the concomitant administration of relatively large amounts of insulin that are required [47].Also in relatively euglycemic DKA, it would require study to learn whether the relatively low initial IRs recommended by the present algorithm would be sufficient.There could be a need for either a priming bolus of insulin or selection of higher MR initial than for other patients of similar body weight.As is true in any case of DKA, intravenous insulin infusion should not be abandoned until ketosis has largely cleared, as indicated by closing of the anion gap.

Advantages of Computerization
Computerization increasingly has been adopted by many hospitals to provide a decision support system guiding insulin adjustments during intravenous insulin infusion [45].A number of studies have shown performance improvement in comparison to manual methods of determining insulin treatment [41,[48][49][50][51][52][53][54][55].Intravenous insulin infusion necessitates frequent monitoring of BG and adjustment of IR and therefore is laborintensive for nursing staff.Although many hospitals restrict use of intravenous insulin infusion to the intensive care unit setting, some hospitals conduct intravenous insulin therapy on general wards.Correct use of paper tabular algorithms requires detailed training of new providers and nursing personnel.Some institutions assign performance of point-of-care BG and adjustment of IR to different personnel.Users are burdened by complexity of the paper tabular protocols, which require making a correct choice among institutional algorithms, assigning the starting column and IR, assigning test times and remembering to perform the frequent point-of-care BGs, adhering to columnchange rules, recognizing the correct cell of the table, adjusting IR rate at the right time (often with supervision of one nurse by another), recording the time of changes of IR, and, in case the rules of the algorithm require interruption of IR for low BG, recognizing the conditions and timing after which the infusion must be re-started.In reality, a reasonable rule for routine administration of ICU medications may provide 30-minute grace period for performance, so that the actual execution may deviate from the intent of the algorithm even in the absence of error.However, actual errors in fact are a prevalent finding, possibly reduced by computerization.It may be that computerization as a decision support tool would enhance nursing satisfaction and utilization of intravenous insulin algorithms [56][57][58].
A further potential advantage of computerization is that the intended mathematical design of the algorithm can be implemented with greater precision than with use of tabular protocols.When full computerization is achieved, the column change rules will be replaced by computation of MR.The stepwise increments of IR according to BG, selected from a table, will be replaced by computation of a continuous function for IR.

Limitations and Areas for Future Clinical Work
There is relatively little literature about the pharmacodynamic delay of insulin under regimens of continuous infusion [59][60][61].Some algorithms may account for pharmacodynamics, but in general, many present day algorithms use the IR

48
A Quadruply-Asymmetric Sigmoid to Describe the Insulin-Glucose Relationship during Intravenous Insulin Infusion (or a multiplier assigning IR according to the elevation of BG above 60 mg/dL), without evaluating pharmacodynamic delay, to interpret observed change of blood glucose.The previous IR or a multiplier is used to interpret changes of BG that have occurred in the past, and an IR or multiplier is assigned to achieve the desired effects for the future.An additional protection for our algorithm, under development, will be consideration of the rate of delivery of the hypoglycemic effect of intravenous insulin so as to account for delay in action of insulin delivered intravenously.It is anticipated that the response of a patient to carbohydrate exposure could be learned so as to introduce a predictive component to the estimation of MR, in case of change of carbohydrate exposure.
It is proposed that the rate of ascent from BG below target, ROD from BG above target, and value of G-per-DIEM could be studied for various populations and under various conditions of treatment, to permit refinement of assumed values.

CONCLUSION
We provide essential steps necessary for the computerization of a programmable continuous intravenous insulin delivery system with the use of a differentiable equation to estimate the IR.The changing insulin requirement and nutritional support in hospitalized patients with hyperglycemia and diabetes result in large glucose oscillations and increased risk of hypoglycemia.The proposed equation that assumes a sigmoidal relationship between IR and BG concentration may offer a safety advantage over the linear relationship currently employed in some intravenous glucose management systems.
Figure 1.Restraining rate, correction rate, and infusion rate as functions of blood glucose.Vertical dashed lines demarcate the boundaries of the target range using parameters of the standard default algorithm (Table6), shown here for MR = 3 units/hr.
Figure 2.Restraining rate as a function of blood glucose (BG) in magnified scale.The tight curvature is shown at BG values just above the inflection point for RR, which is at BG 120 mg/dL for the standard default algorithm, shown here for MR = 3 units/hr.

Figure 3 :
Figure 3: Relationships between insulin infusion rate (IR) and blood glucose (BG) for three medical conditions treated in the hospital with intravenous insulin infusion, with target ranges demarcated by vertical dashed lines that correspond to hyperglycemic hyperosmolar state (a), diabetic ketoacidosis (b), and standard default critical care use (c).The curves shown in top and middle panel in descending order correspond to maintenance rate (MR) = 5, 4, 3, 2, or 1 unit/hr; and in the lower panel correspond to MR = 8, 6, 4, 3, 2, 1.5, and 1 unit/hr.The initial maintenance rate (MR initial ) is 3 units/hr for (a) and (b), and 1.5 units/hr for (c).

Table 5 . Derivation of the equation for IR next, for given MR and possible values of BG at one point in time (Continued) Solve for "c" using a designer-defined Solve for "d" using a designer-defined point for RR point for CR
Replace the variables in eqn.20 with the known Replace the variables in eqn.21 with the known designer-defined point = (70 mg/dL, RR @BG70 ) designer-defined point = (BG target_plus, ∆ MR ); = (BG TP , ∆ MR ); (BG target_plus = BG TP = BGTP) Solve