Designing a general function that describes pharmacodynamic data

Question

I am interested in coming up with a function describing the dose-dependent action of a drug over time. In this specific case, the image attached below shows the glucose infusion rate (GIR) of the insulin "Lyumjev", where the GIR represents the amount of glucose that needs to be infused into a patient to keep that patient's blood glucose level constant. It's therefore a sort of surrogate marker for insulin action.

If you look at the image, you can see that that there's an extremely steep, linear climb during the first hour, followed by a broad peak or plateau before the required GIR slowly drops back to zero. Depicted are the GIRs for three doses (7, 15 or 30 units of insulin), and if you look closely you can see that the plateaus become higher and broader, and the duration of insulin action becomes longer in dependence of the dose.

I could use simple software like Excel to describe each of the three curves using, e.g., simple polynomial fits, but this won't let me predict the likely curve for, say, 5 U or 10 U of this insulin. To really solve this issue, I will have to come up with a general function that takes the insulin dose as an argument, but I have no idea how to do this and where to start.

What I know

1. The dose (I would like to plot curves for whichever dose I enter)
2. Area under the curve (put simply, the AUC more or less scales linearly with the dose such that 2 U of insulin will have double the AUC of 1 U, but the distribution will be slightly shifted towards later hours as shown in the image)

I am not sure if my question is naive, but is there a way to "design" an equation like this? I did check pharmacokinetic equations that I found on the web, but the issue is that these normally describe serum drug levels over time, and those I have tried did not alter the distribution of the drug over time, but instead only altered the height of the peak. I appreciate any hint you may have on where to start, or where to look... thanks a lot in advance!

Glucose infusion rate over time, plotted for three insulin doses:

The picture is taken from the EMA's assessment report for the drug Lyumjev. (https://www.ema.europa.eu/en/documents/assessment-report/liumjev-epar-public-assessment-report_en.pdf page 69)

Answer 1

While this does seem like the kind of phenomenon that has been modeled already, just from looking at the graphs I might suggest functions of the form $xe^{-x}$ with various scaling constants inserted. For example, here are the graphs of $xe^{-x/\sqrt{7}}$ (lowest curve), $xe^{-x/\sqrt{15}}$, and $xe^{-x/\sqrt{30}}$ (highest curve), which seem like reasonable facsimiles of your plot. The peaks of the functions occur at times $\sqrt7$, $\sqrt{15}$, and $\sqrt{30}$, respectively.

Notice that for small values of time ($x$-coordinate), all these functions look identical; this would be appropriate if the initial GIR is independent of the dose (and differences only appear over time), but that might not be accurate. The function $Cxe^{-x/\sqrt B}$ has a graph under which the total area equals $BC$, so you can "split up" the doses $D$ into the product of two positive factors $D=BC$ in different ways; the slope of the graph at $x=0$ will equal $C$, while the peak of the graph will be at $x=\sqrt B$. (The above examples all have $C=1$ and $B=D$.)

Answer 2

Caveat: as this is a medical issue, mathematics cannot predict what will really happen, and all reccomendations should be taken with a grain of salt (or, since we are dealing with Insulin, a grain of sugar).

You might try a "scaling" approach: assume that all the plots are related to each other though a scale. That is, assume there is a "base response" function $GIR_{base}=F(t)$ and then relate the various dosage plots to the base response by scaling time and total response $$GIR(t,{\text dose})= b({\text dose}) F({t \over a({\text dose})})$$ For example, pick the 30U response as the base function (since it is the largest, and probably has the smallest relative error, and then try to scale it until it can be superimposed onto the other doses, 15U, 7U. You can estimate $a$ independently by adjusting until the peak time positions match, and then scale the heights using $b$ until the peak GIR match (alternatively, until the AUC matches). The found values of $a,b$ can then be interpolated to estimate the graph for other dosages. This is similar with the proposal of @Greg Martin, but with a non parametric base response model.

EDIT

To clarify what I mean by setting the base function to the 30U response, I mean using the values aa a lookup table. that is, treat the 30U as a table [time, height]. If you want $a$ to be say $2$, make a new table, then read the 30U response at time $t$ , and place it at time $t\over2$ in the new, scaled table. After you have went over all times, multiply the height values in the new table by $b$.