I modeled GLP-1 pharmacokinetics in TypeScript (and open-sourced it)

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·Δ)

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·Δ))

• 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 + 1 e - 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 );
}

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

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 ;
}

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 * 864 e5 ,
  pk ,
  200
);
// curve → [{ t, mg }, …] ready to plot

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.