Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › Resolving Hysteresis in a Bronchodilator Case

Worked Case Studies · Tutorial 574 of 598

Resolving Hysteresis in a Bronchodilator Case

Learn how to diagnose a delayed bronchodilator response from a time-course plot, calculate effect-site concentrations, and assess whether an estimated equilibration rate supports a delayed-effect interpretation.

Advanced 11 min read

What You'll Learn

  • Diagnose counter-clockwise hysteresis from the ordering of concentration and response over time
  • Specify the effect-compartment model and its key assumptions
  • Calculate effect-site concentrations under a piecewise-constant exposure approximation
  • Fit and compare candidate equilibration rates using residual sums of squares
  • Interpret the equilibration rate and its implied half-time without overstating the evidence
  • Recognize when sampling, model assumptions, or sparse data limit interpretation

Why a Bronchodilator Response Can Lag Concentration

After a bronchodilator dose, plasma concentration may rise and then fall, while bronchodilation develops more slowly and persists after concentration has begun to decline. A plot of response against concentration can then trace a loop rather than a single exposure-response curve. The order in which the points occur matters: the same concentration can be associated with different responses on the rising and falling portions of the time course.

The case below uses hypothetical concentration measurements and baseline-corrected bronchodilation responses. The question is not merely whether concentration and response are associated. It is whether the timing pattern is consistent with a delayed response, and whether a simple effect-compartment model can describe that delay.

Definition: Counter-clockwise hysteresis is a time-ordered concentration-response loop in which, at a given concentration, response is generally greater on the declining-concentration portion of the time course than on the rising portion. It is consistent with response lagging concentration, but by itself does not prove the physiological mechanism or exclude other time-varying processes.

Hypothetical Case Data and the Loop

Suppose concentration is measured in hypothetical concentration units and bronchodilation is recorded as a baseline-corrected fraction of the maximum possible response. The six measurements are:

Time (h)ConcentrationBronchodilation response
000.000
14.00.226
28.00.405
45.00.431
63.00.380
81.50.290

Concentration peaks at 2 h, whereas the largest recorded response is at 4 h. The time ordering therefore suggests a lag. To see the loop direction more directly, linearly interpolate between the 4 h and 6 h measurements: at concentration 4.0, the interpolated response on the declining limb is about \((0.431+0.380)/2=0.4055\). At the measured concentration of 4.0 on the rising limb at 1 h, response is 0.226. The larger response at the same approximate concentration on the declining limb is consistent with counter-clockwise hysteresis. This interpolation is illustrative; it does not add an observed data point or account for measurement uncertainty.

Worked Example: Diagnosing the Hysteresis Direction

At concentration 4.0, the rising-limb response is 0.226 at 1 h. On the declining limb, concentration falls from 5.0 at 4 h to 3.0 at 6 h. Estimate the response at concentration 4.0 by linear interpolation and compare the two responses.

Concentration 4.0 lies halfway between 5.0 and 3.0. The linearly interpolated response is therefore \(0.431+(0.380-0.431)(4.0-5.0)/(3.0-5.0)=0.4055\).

The declining-limb estimate, 0.4055, is greater than the rising-limb observation, 0.226, by \(0.4055-0.226=0.1795\) response units. This ordering supports counter-clockwise hysteresis. It is evidence of a time-dependent concentration-response relationship, not proof that delayed equilibration is the only explanation.

Effect-Compartment Model

A standard way to represent a lag between measured concentration and response is to posit an effect-site concentration, \(C_e(t)\), that approaches the measured concentration, \(C(t)\), at a finite rate. The effect-site concentration is a model quantity; it is not necessarily a directly measured concentration in a physical compartment.

Model: The effect-compartment model is \(dC_e(t)/dt=k_{e0}\{C(t)-C_e(t)\}\), where \(k_{e0}>0\) is the equilibration rate in inverse hours. For the response in this case, use the baseline-corrected Emax relationship \(E(t)=E_{\max}C_e(t)/\{EC_{50}+C_e(t)\}\). Here \(E_{\max}\) is the maximum response fraction and \(EC_{50}\) is the effect-site concentration producing half of \(E_{\max}\).

A larger \(k_{e0}\) makes \(C_e(t)\) track \(C(t)\) more quickly. As \(k_{e0}\) becomes very large, the model approaches an immediate concentration-response relationship, \(C_e(t)\approx C(t)\). A smaller positive rate produces a greater lag. This model assumes a constant equilibration rate, a stable response relationship, and no additional time-varying influence such as tolerance or a changing baseline. The case calculation also treats the measured concentrations as known and the response residuals as suitable for least-squares fitting.

Calculating Effect-Site Concentrations

Concentrations are observed only at discrete times. For a transparent worked fit, assume concentration is constant within each interval and equal to its value at the interval endpoint. If concentration is \(C_i\) for an interval of length \(\Delta t_i\), solving the differential equation gives:

$$ C_{e,i}=C_i+(C_{e,i-1}-C_i)\exp(-k_{e0}\Delta t_i). $$

This is an exact solution under the stated piecewise-constant input assumption. It is not an exact solution for a continuously changing concentration profile. In practice, denser PK samples or a specified interpolation of the concentration-time curve may be preferable.

Worked Example: Propagating Effect-Site Concentration

Use \(k_{e0}=0.50\ \mathrm{h}^{-1}\), \(C_e(0)=0\), and the piecewise-constant interval convention. Calculate \(C_e\) at 1 h and 2 h.

For the first interval, concentration is 4.0 for 1 h. Thus \(C_e(1)=4.0+(0-4.0)\exp(-0.50\times1)=4.0(1-0.6065)=1.574\) concentration units.

For the second interval, concentration is 8.0 for 1 h. Therefore \(C_e(2)=8.0+(1.574-8.0)\exp(-0.50\times1)\approx 4.102\) concentration units.

Effect-site concentration is below measured concentration at both times, but it is rising toward it. That lag is the feature the effect compartment is intended to represent.

Fitting the Equilibration Rate

For this illustration, set \(E_{\max}=0.80\) and \(EC_{50}=4.0\) concentration units in advance, and estimate only \(k_{e0}\). This isolates the rate calculation; it does not establish that the fixed response parameters are known in a real analysis. For each candidate rate, recursively calculate \(C_e\), predict response, and calculate the residual sum of squares:

$$ \widehat{E}_i=\frac{0.80C_{e,i}}{4.0+C_{e,i}}, \qquad RSS(k_{e0})=\sum_i\{E_i-\widehat{E}_i(k_{e0})\}^2. $$

The rate is selected to minimize \(RSS\) over positive values. The comparison uses the same observations and fixed response parameters for every candidate. Least squares is most directly interpreted when response errors are independent with a common variance; repeated measurements may instead have correlated errors, which this simple calculation does not model.

Worked Example: Comparing Candidate Equilibration Rates

For \(k_{e0}=0.25\ \mathrm{h}^{-1}\), applying the recursion at times 1, 2, 4, 6, and 8 h gives predicted responses of approximately 0.145, 0.305, 0.371, 0.361, and 0.314. The corresponding residuals, observed minus predicted, are about 0.081, 0.100, 0.060, 0.01967, and \(-0.024\). Thus \(RSS\approx0.081^2+0.100^2+0.060^2+0.01967^2+(-0.024)^2\approx0.02112\), with unrounded predictions giving \(RSS=0.02121\).

For \(k_{e0}=1.00\ \mathrm{h}^{-1}\), the predictions are approximately 0.310, 0.480, 0.450, 0.361, and 0.243. Using the unrounded predictions in the same residual calculation gives \(RSS=0.01554\). The latter candidate fits better than \(0.25\ \mathrm{h}^{-1}\), but this comparison alone does not show it is optimal.

Minimizing the same one-dimensional objective over positive rates gives a rate close to \(0.50\ \mathrm{h}^{-1}\) for these hypothetical data and fixed parameters. At that rate, the predictions are close to the five observed responses. The decrease in RSS relative to the two displayed candidates is consistent with the time lag seen in the data. Because responses are rounded and there are only five post-baseline observations, the estimate should not be interpreted as precisely determined.

What the Estimated Rate Does and Does Not Mean

The estimated \(k_{e0}\) describes the speed at which modeled effect-site concentration approaches plasma concentration under this model. Its reciprocal is a characteristic time scale, while the time for a first-order effect compartment to close half the gap between \(C_e\) and a fixed \(C\) is \(\ln(2)/k_{e0}\). At \(k_{e0}\approx0.50\ \mathrm{h}^{-1}\), that half-time is approximately \(0.693/0.50=1.39\) h. This is a model-based equilibration half-time, not the half-life of the drug in plasma or a direct measurement of receptor binding.

Justification: If \(C\) is held constant, the differential equation has solution \(C_e(t)=C+\{C_e(0)-C\}\exp(-k_{e0}t)\). The remaining gap is \(C-C_e(t)=\{C-C_e(0)\}\exp(-k_{e0}t)\). Setting this gap to half its initial value gives \(\exp(-k_{e0}t)=1/2\), hence \(t=\ln(2)/k_{e0}\).

A finite estimated rate and a better fit than the direct-effect limit support a delayed-effect description of these data. They do not identify the biological location of the delay. The same pattern could arise from an unmeasured active metabolite, measurement timing, a changing baseline, or another time-dependent process. The rate may also depend on assumptions about concentration between sample times and on the fixed values of \(E_{\max}\) and \(EC_{50}\). A real analysis should examine parameter uncertainty, residual patterns, sensitivity to PK interpolation, and whether allowing response parameters to be estimated changes the rate.

Diagnostics, Alternatives, and Common Pitfalls

Useful checks begin with the time-course plot and the residuals, not just the minimized objective. Plot response against time and concentration against time, then connect concentration-response points in chronological order. For a fitted model, inspect observed-minus-predicted residuals over time and concentration. A persistent signed residual pattern, a missed peak, or large deviations at the late samples can indicate a poor input approximation, a misspecified response curve, or a process beyond simple equilibration. Sparse sampling can also make the loop direction or rate difficult to establish.

A direct concentration-response model is the main simpler alternative: it assumes response follows concentration without a lag. It is parsimonious, but cannot represent different responses at the same concentration on rising and falling limbs. A more flexible time-varying response curve can describe a loop, but may not provide a compact, interpretable equilibration rate and can overfit sparse data. The effect-compartment model offers a mechanistic-looking summary with few parameters, but its effect-site state is latent and its rate is conditional on the model assumptions.

  • Pitfall: treating a loop as proof of an effect compartment. The wrong reasoning is that counter-clockwise hysteresis uniquely establishes delayed equilibration. The correction is that the loop supports a time-dependent response; alternative biological and measurement explanations must still be considered.
  • Pitfall: reading the rate as a plasma elimination rate. \(k_{e0}\) governs modeled equilibration between measured concentration and effect-site concentration, not plasma elimination.
  • Pitfall: reporting a precise rate from a small case data set. An optimizer returns a value even with sparse, rounded observations. Report the conditional estimate, assess uncertainty and sensitivity, and avoid interpreting numerical precision as evidential certainty.
Key takeaway: In this hypothetical bronchodilator case, the response peaks after measured concentration and is higher on an interpolated declining limb than on the rising limb at the same concentration. An effect-compartment fit with a finite equilibration rate is consistent with that delay, but the interpretation remains conditional on the exposure interpolation, response model, and limited data.

Check Your Understanding

Answer each question using the definitions and assumptions above. Short solutions follow the questions.

  1. At a matched concentration, response is 0.20 on the rising limb and 0.35 on the declining limb. What hysteresis direction does this support, and what does it not prove?
  2. For \(k_{e0}=0.50\ \mathrm{h}^{-1}\), \(C_e=2.0\), and a constant concentration of 6.0 over the next hour, calculate the new effect-site concentration.
  3. If \(k_{e0}=0.70\ \mathrm{h}^{-1}\), calculate the equilibration half-time in hours.
  4. Why must candidate rates be compared using the same observations and the same fixed response parameters?
  5. Name one diagnostic sign that the effect-compartment model or its input approximation may be inadequate.

Solutions. (1) The ordering supports counter-clockwise hysteresis because response is greater on the declining limb; it does not prove delayed equilibration is the unique cause. (2) \(C_{e,\mathrm{new}}=6.0+(2.0-6.0)\exp(-0.50)=6.0-4.0(0.6065)=3.574\) concentration units. (3) \(t_{1/2}=\ln(2)/0.70=0.990\) h. (4) Otherwise differences in RSS could reflect changed data or response assumptions rather than the candidate equilibration rate. (5) Examples include systematic residuals over time, a missed response peak, or sensitivity of the estimate to how concentration is interpolated between samples.

Clinical Biostats