Every "medication level" chart I saw in a GLP-1 tracking app was wrong in the same two ways. So while building Tiro, a tracker for people on semaglutide and tirzepatide, I wrote a proper pharmacokinetic model — and I just pulled the pure math out into a zero-dependency package: glp1-pk.

Here's the interesting part of the problem.

The two bugs in every naive "drug level" curve

If you plot "how much drug is in your body" as a single exponential decay, you get this:

level(Δ) = dose · e^(−k·Δ)

Enter fullscreen mode Exit fullscreen mode

That's wrong twice.

Bug 1 — it ignores absorption. A subcutaneous injection isn't instantly in your bloodstream. The level rises to a peak over hours or days, then falls. A pure decay starts at maximum, which never happens with a depot injection.

Bug 2 — it ignores route. Oral semaglutide has a bioavailability of roughly 0.8%. Injected, it's about 89%. So "14 mg" taken orally puts about two orders of magnitude less drug on board than 14 mg injected. A decay curve keyed only on milligrams gets this exactly backwards.

The fix: the Bateman function

A one-compartment model with first-order absorption and first-order elimination gives you the Bateman function. For a single dose D at elapsed time Δ:

c(Δ) = D · kA/(kA − kE) · (e^(−kE·Δ) − e^(−kA·Δ))

Enter fullscreen mode Exit fullscreen mode

  • kE is the elimination rate — ln(2) / half_life. It's a property of the compound (semaglutide's half-life is ~1 week; liraglutide's is ~13 hours).
  • kA is the absorption rate, solved from the drug's time-to-peak (tMax), which depends on route — a slow subcutaneous depot peaks in days, an oral dose in ~1 hour.

The catch: you can't invert tMax = ln(kA/kE) / (kA − kE) for kA in closed form. It's transcendental. So you solve it numerically with bisection (tMax is strictly decreasing in kA, so it's bulletproof):

export function absorptionRateFromTmax(target, kE) {
  const ceiling = 1 / kE;                 // limit of tMax as kA → kE⁺
  const t = target >= ceiling ? 0.98 * ceiling : target;
  let lo = kE * (1 + 1e-12), hi = kE * 2;
  while (tMaxFromRates(hi, kE) > t) hi *= 2;   // bracket the root
  for (let i = 0; i < 200; i++) {
    const mid = 0.5 * (lo + hi);
    if (tMaxFromRates(mid, kE) - t > 0) lo = mid; else hi = mid;
  }
  return 0.5 * (lo + hi);
}

Enter fullscreen mode Exit fullscreen mode

Two things I care about here:

  • Totality. A real kA > kE only exists when tMax < 1/kE. Rather than return NaN on an out-of-range input, it clamps to just under the ceiling so you always get a finite, sane rate. Health-adjacent code should never surface a NaN to a chart.
  • The kA ≈ kE singularity. When the two rates are nearly equal the formula divides by ~0, so the level function falls back to the L'Hôpital limit D · kE · Δ · e^(−kE·Δ).

From "mg on board" to nmol/L

Milligrams-on-board is fine for a relative curve, but to compare against published exposures you want a concentration:

C[mg/L]   = absorbedMg / Vd            // volume of distribution
C[nmol/L] = C[mg/L] · 1e6 / molarMass  // mg → nmol

Enter fullscreen mode Exit fullscreen mode

The molar mass matters more than you'd guess. Dulaglutide is a ~59.7 kDa antibody-Fc fusion — about 15× heavier than semaglutide (~4.1 kDa). For the same mass on board, its molar concentration is ~15× lower. A model that hard-codes one "nmol per mg" factor is wrong for every compound but one.

Superposition = a real dose history

Because the model is linear, a full history is just the sum of each dose's curve:

export function levelAt(doses, t, pk) {
  let sum = 0;
  for (const d of doses) sum += doseLevelAt(d.amountMg, t - d.takenAt, pk);
  return sum;
}

Enter fullscreen mode Exit fullscreen mode

Doses in the future contribute zero (negative elapsed time → guarded to 0). That's the whole "estimated medication level" line.

Using it 

import { pkFor, levelAt, sampleLevelSeries } from "glp1-pk";

const pk = pkFor("tirzepatide", "injection");
const doses = [
  { amountMg: 2.5, takenAt: Date.parse("2026-06-01T09:00:00Z") },
  { amountMg: 5.0, takenAt: Date.parse("2026-06-08T09:00:00Z") },
];

const mgNow = levelAt(doses, Date.now(), pk);
const curve = sampleLevelSeries(doses, Date.now(), Date.now() + 14 * 864e5, pk, 200);
// curve → [{ t, mg }, …] ready to plot

Enter fullscreen mode Exit fullscreen mode

Zero dependencies, ships types, pure functions, 16 tests on Node's built-in runner (no jest). MIT: github.com/navidmosleminiya/glp1-pk.

The disclaimer that actually matters

The per-compound constants (half-life, tMax, bioavailability, Vd, molar mass) are central population-PK estimates from public literature and labels. Individual pharmacokinetics vary a lot. This is a tracking/visualisation/teaching tool — an estimate from logged doses, not a measured concentration, and not a dosing tool. That framing is load-bearing for anything health-adjacent, and it's in the README and the types.

I build this full-time on Tiro, a GLP-1 companion that unifies your shots, protein, and a private body scan. If you want the model without the app, the package is right there. PRs on the parameter tables welcome — especially if you have better retatrutide numbers.