Tutorials › Pharmacokinetic-Pharmacodynamic (PK/PD) Modeling › Integrating-Factor Solution for Time-Varying Drug Effects

Indirect Response Models · Tutorial 129 of 598

Integrating-Factor Solution for Time-Varying Drug Effects

Learn to derive and apply one exact solution for indirect response when production, loss, or both vary continuously over time.

Advanced 10 min read

What You'll Learn

  • State the assumptions for a linear indirect-response model with time-varying production and loss.
  • Derive the exact response trajectory from an integrating factor.
  • Calculate responses when production varies, loss varies, or both vary.
  • Check the units, signs, and integrated loss terms in a solution.
  • Recognize when constant-rate or numerical approaches are more suitable.

Why Use an Integrating Factor?

In an indirect response model, a drug changes the rate at which a response is produced, the rate at which it is lost, or both. When those rates change continuously, a single constant-rate formula does not describe the full trajectory. The integrating-factor method gives an exact solution for the linear model while allowing either rate to vary with time.

The method is useful when a time course of drug exposure implies continuously changing production or loss, and the rates can be specified as functions of time. It also makes the roles of the initial response, production history, and accumulated loss explicit. The earlier material on piecewise-constant exposure and recursive solutions handles intervals of constant rates; here the rates need not be constant within an interval.

Definition: Let \(R(t)\) be the response at time \(t\), \(q(t)\) the production rate, and \(k(t)\) the first-order loss-rate constant. The time-varying linear indirect-response model is $$ \frac{dR(t)}{dt}=q(t)-k(t)R(t), \qquad R(t_0)=R_0. $$ Here \(q(t)\) has units of response per unit time, \(k(t)\) has units of inverse time, and \(R_0\) is the response at the starting time \(t_0\).

Drug effects may be represented through either rate function. For example, inhibition or stimulation of production changes \(q(t)\); inhibition or stimulation of loss changes \(k(t)\). The equation assumes that the response is produced at the specified rate and that loss is proportional to the current response. It does not add a separate delay mechanism beyond the turnover of \(R(t)\).

For the derivation and examples, assume \(q(t)\) and \(k(t)\) are integrable over the time interval of interest. Continuous functions satisfy this condition, as do piecewise-continuous functions. For ordinary drug-induced production and loss, \(q(t)\geq 0\) and \(k(t)\geq 0\). The derivation itself requires integrability, not nonnegativity.

Deriving the General Solution

The equation is linear in \(R(t)\). Its integrating factor is chosen to cancel the term involving the product of the loss rate and the response. Define

$$ M(t)=\exp\left(\int_{t_0}^{t} k(u)\,du\right). $$

Since \(dM(t)/dt=k(t)M(t)\), the product rule gives

$$ \frac{d}{dt}\{M(t)R(t)\} =M(t)\frac{dR(t)}{dt}+k(t)M(t)R(t) =M(t)\{q(t)-k(t)R(t)\}+k(t)M(t)R(t) =M(t)q(t). $$

Integrating both sides from \(t_0\) to \(t\), and using \(M(t_0)=1\), yields

$$ M(t)R(t)-R_0=\int_{t_0}^{t}M(s)q(s)\,ds. $$

Divide by \(M(t)\). Because \(M(s)/M(t)=\exp\{-\int_s^t k(u)\,du\}\), the resulting solution is

$$ R(t)=R_0\exp\left(-\int_{t_0}^{t}k(u)\,du\right) +\int_{t_0}^{t}q(s)\exp\left(-\int_s^t k(u)\,du\right)\,ds. $$

The first term is the part of the initial response remaining after loss has accumulated from \(t_0\) to \(t\). The second term adds production at each earlier time \(s\), with that contribution reduced by loss accumulated between \(s\) and \(t\). The inner integral is always over the interval from the time a contribution enters to the time the response is evaluated; reversing those limits would give the wrong sign.

Key takeaway: The same solution covers time-varying production, time-varying loss, and simultaneous variation in both. Evaluate the accumulated loss over the correct time interval, then account for production across the interval.

Worked Applications

Worked Example: Continuously Increasing Production

Suppose a hypothetical response starts at \(R_0=10\) response units at \(t_0=0\). Production rises linearly, \(q(t)=(4+2t/(1\,\mathrm{h}))\) response units per hour, while loss is constant at \(k(t)=0.5\,\mathrm{h}^{-1}\). Find \(R(2\,\mathrm{h})\).

The production function is nonnegative over the interval, the loss rate is positive, and both functions are continuous. Thus the model assumptions hold. For the initial-response term, the integrated loss is \(0.5\times2=1\):

$$ R_0\exp\left(-\int_0^{2\,\mathrm{h}}0.5\,du\right) =10e^{-1} =3.6788\ \text{response units}. $$

For production, let \(x=s/(1\,\mathrm{h})\), so \(x\) runs from 0 to 2. The production contribution is

$$ \int_0^{2\,\mathrm{h}}q(s)e^{-0.5(2\,\mathrm{h}-s)/(1\,\mathrm{h})}\,ds =\int_0^2(4+2x)e^{-0.5(2-x)}\,dx. $$

With \(u=2-x\), this becomes \(\int_0^2(8-2u)e^{-0.5u}\,du\). Evaluating the two terms,

$$ 8\int_0^2e^{-0.5u}\,du-2\int_0^2ue^{-0.5u}\,du =8\{2(1-e^{-1})\}-2\{4-8e^{-1}\} =8.000\ \text{response units}. $$

Therefore \(R(2\,\mathrm{h})=3.6788+8.000=11.6788\), or approximately 11.68 response units. The response can exceed its initial value because production over the interval outweighs the loss.

Worked Example: Continuously Varying Loss

Now let production be constant at \(q(t)=6\) response units per hour, with \(R_0=5\) response units. Let the loss rate decline according to \(k(t)=\{1/[1+t/(1\,\mathrm{h})]\}\,\mathrm{h}^{-1}\). Find the response at \(t=2\,\mathrm{h}\).

Both rates are continuous on the interval, production is positive, and the loss rate remains positive. For the initial response, the accumulated loss is

$$ \int_0^{2\,\mathrm{h}} k(u)\,du =\int_0^2\frac{dx}{1+x} =\ln(3), $$

where \(x=u/(1\,\mathrm{h})\). The remaining initial contribution is therefore \(5e^{-\ln(3)}=5/3=1.6667\) response units. For a production contribution entering at time \(s\), the attenuation to \(2\,\mathrm{h}\) is

$$ \exp\left(-\int_s^{2\,\mathrm{h}}k(u)\,du\right) =\frac{1+s/(1\,\mathrm{h})}{3}. $$

Thus the production contribution is

$$ 6\frac{\text{response units}}{\mathrm{h}} \int_0^{2\,\mathrm{h}}\frac{1+s/(1\,\mathrm{h})}{3}\,ds =6\cdot\frac{1}{3}\cdot4\,\mathrm{h} =8\ \text{response units}. $$

Adding the two parts gives \(R(2\,\mathrm{h})=1.6667+8=9.6667\), or approximately 9.667 response units. In this example, the loss rate declines over time, so later production is reduced less before \(t=2\,\mathrm{h}\) than earlier production.

Worked Example: Both Production and Loss Vary

Consider a hypothetical model with \(R_0=4\) response units, \(k(t)=\{1/[1+t/(1\,\mathrm{h})]\}\,\mathrm{h}^{-1}\), and \(q(t)=2\{1+t/(1\,\mathrm{h})\}\) response units per hour. Calculate \(R(2\,\mathrm{h})\).

The rates are continuous and positive throughout the interval, so they satisfy the stated assumptions. The initial response is attenuated by \(e^{-\ln(3)}=1/3\), giving \(4/3=1.3333\) response units. As in the previous calculation, production added at \(s\) is reduced by the factor \(\{1+s/(1\,\mathrm{h})\}/3\). With \(x=s/(1\,\mathrm{h})\), the production contribution is

$$ \int_0^2 2(1+x)\frac{1+x}{3}\,dx =\frac{2}{3}\int_0^2(1+x)^2\,dx =\frac{2}{9}\{3^3-1^3\} =\frac{52}{9} =5.7778\ \text{response units}. $$

Therefore

$$ R(2\,\mathrm{h})=\frac{4}{3}+\frac{52}{9} =\frac{64}{9} =7.111\ \text{response units}. $$

This example checks both time-varying functions in one calculation. The production contribution has units of response units because the production rate is multiplied by a time interval; the attenuation factor is dimensionless.

Assumptions, Checks, and Alternatives

The solution is exact for the stated differential equation, not for every indirect-response model. It relies on a single response state, production that enters additively as \(q(t)\), and first-order loss \(k(t)R(t)\). If loss is nonlinear in the response, if the rate depends on the response in another way, or if unmodelled delay processes matter, this integrating-factor result does not apply directly.

Several checks help catch errors before interpreting a trajectory:

  • Check units. The exponent must be dimensionless: \(k(u)\,du\) has units \(\mathrm{time}^{-1}\times\mathrm{time}\). The production integral must have response units.
  • Check the limits. A contribution produced at \(s\) is affected by loss from \(s\) to the evaluation time \(t\), not from \(t\) back to \(s\).
  • Check signs and bounds. If \(k(t)\geq0\), each exponential attenuation factor lies between 0 and 1 for \(t\geq s\). A factor greater than 1 signals a likely sign or limit error.
  • Check the initial condition. Setting \(t=t_0\) must return \(R(t_0)=R_0\); the production integral then has zero length.
  • Check numerical implementation. Compare the formula with a numerical solution of the differential equation at selected times. Large disagreement can indicate a time-unit mismatch, inaccurate integration of a rapidly changing rate, or an incorrectly coded rate function.

A changing rate does not imply that the response must move monotonically. Its instantaneous slope is \(q(t)-k(t)R(t)\): the response rises when production exceeds loss at that moment and falls when loss exceeds production. Nor does the integrating-factor expression remove uncertainty in \(q(t)\), \(k(t)\), or the measured response; it gives the trajectory conditional on those inputs and the model assumptions.

For constant rates, the familiar constant-rate solution is simpler to evaluate and is a useful check. For piecewise-constant rates, recursive solutions can be convenient when each interval has a known rate. When rates vary continuously but the integrals cannot be evaluated analytically, numerical quadrature of the exact expression or a numerical differential-equation solver can be used. These numerical methods approximate the solution; they do not change the assumptions of the model. A numerical result should be checked for adequate time-step or integration accuracy, especially when rates change rapidly.

Common mistake: Treating the current loss rate \(k(t)\) as though it had applied throughout the whole interval. With time-varying loss, the relevant quantity is the accumulated rate \(\int_s^t k(u)\,du\). Replacing it with \(k(t)(t-s)\) is valid only if loss is constant over that interval.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

  1. In the model \(dR/dt=q(t)-k(t)R(t)\), what units must \(q(t)\) and \(k(t)\) have if \(R(t)\) is measured in response units and time in hours?
  2. For \(R_0=2\) response units, constant \(q=3\) response units per hour, and constant \(k=0.5\,\mathrm{h}^{-1}\), calculate \(R(1\,\mathrm{h})\).
  3. If \(q(t)=0\), \(R_0=12\) response units, and \(k(t)=\{2/[1+t/(1\,\mathrm{h})]\}\,\mathrm{h}^{-1}\), calculate \(R(3\,\mathrm{h})\).
  4. If \(R(t_0)=6\) response units, \(q(t_0)=4\) response units per hour, and \(k(t_0)=0.5\,\mathrm{h}^{-1}\), what is the instantaneous slope at \(t_0\)? What does its sign imply at that instant?

Solutions. (1) \(q(t)\) must have units of response units per hour, and \(k(t)\) must have units of inverse hours, so both terms in the differential equation have units of response units per hour. (2) \(R(1)=2e^{-0.5}+3(1-e^{-0.5})/0.5=1.2131+2.3608=3.574\) response units. (3) The accumulated loss is \(\int_0^{3\,\mathrm{h}}k(u)\,du=2\ln(4)\), so \(R(3)=12e^{-2\ln(4)}=12/16=0.7500\) response units. (4) \(dR/dt=4-(0.5)(6)=1\) response unit per hour, so the response is increasing at \(t_0\).

Clinical Biostats