Transit Chains Under Time-Varying Drug Exposure

Why Time-Varying Inputs Need a Response Calculation

A drug concentration may rise during absorption, fall during elimination, or vary repeatedly after successive doses. In a transit model, the signal does not pass through the cascade instantaneously: each stage retains and delays part of the incoming signal. Consequently, the output at a given time depends not only on the current concentration but also on the history of the cascade input. This tutorial treats the cascade after its input has been specified, and shows how to calculate that history-dependent response.

The methods below apply to a linear, time-invariant transit chain: the stage rates are constant, each stage transfers signal linearly, and the input enters the first stage. The input itself may vary with time. For example, it may be a concentration-derived signal $u(t)=f(C(t))$. A nonlinear concentration-to-input relationship can be applied before the cascade, but it does not make the cascade output linear in concentration.

Definition: For a chain of $m$ transit stages, let $x_j(t)$ denote the signal in stage $j$, $k_j>0$ its transit rate, and $u(t)$ the input signal. The linear transit-chain equations are $$ \frac{dx_1}{dt}=k_1\bigl(u(t)-x_1(t)\bigr), \qquad \frac{dx_j}{dt}=k_j\bigl(x_{j-1}(t)-x_j(t)\bigr),\quad j=2,\ldots,m. $$ Each $k_j$ has units of inverse time. If $u$ and the $x_j$ are expressed in the same signal units, the equations preserve those units.

A changing input is handled by evaluating $u(t)$ over time, rather than substituting a single concentration into a steady-state expression. If the input is derived from concentration, the concentration profile and the function $f$ must be specified. The transit equations then propagate that input through the stages.

From the ODEs to Convolution

For clarity, first assume the chain has zero initial states at time 0. The first-stage equation can be solved by an integrating factor: $$ x_1(t)=\int_0^t k_1 e^{-k_1(t-s)}u(s)\,ds. $$ Thus stage 1 is a weighted accumulation of past input. For a later stage, the same solution method gives $$ x_j(t)=\int_0^t k_j e^{-k_j(t-s)}x_{j-1}(s)\,ds. $$ Substituting the preceding stage’s integral shows that each stage adds another filtering operation. The total response can be written as one convolution:

x_m(t)=\int_0^t h_m(t-s)u(s)\,ds,

where $h_m(t)$ is the output at time $t$ following a unit impulse into the first stage. This is the convolution representation of the linear chain. It makes explicit that a signal at time $s$ contributes to the output at $t$ according to its age $t-s$.

Result: For distinct transit rates, the impulse response is $$ h_m(t)=\left(\prod_{j=1}^m k_j\right)\sum_{i=1}^m \frac{e^{-k_i t}}{\displaystyle\prod_{\substack{j=1\\j\ne i}}^m(k_j-k_i)},\quad t\geq 0. $$ If all rates are equal to $k$, the limiting expression is $$ h_m(t)=\frac{k^m t^{m-1}e^{-kt}}{(m-1)!},\quad t\geq 0. $$

The distinct-rate formula follows by taking the Laplace transform of the chain with zero initial states. Each stage contributes a factor $k_j/(s+k_j)$, so the input-to-output transfer function is $\prod_j k_j/(s+k_j)$. Partial fractions give the displayed sum. When rates coincide, the partial-fraction terms merge into a repeated pole; taking the equal-rate limit gives the stated polynomial-times-exponential form. For two equal-rate stages, for example, $h_2(t)=k^2te^{-kt}$.

If the initial states are not zero, the solution also contains the response generated by those states in the absence of further input. Convolution of the input alone is therefore sufficient only when the initial states are zero, or when their separate contribution is included. Initial states are especially important when a simulation starts between doses or after a prior period of treatment.

Worked Example: A One-Unit Pulse Followed by Withdrawal

Consider two stages with $k_1=k_2=1$ per minute and $x_1(0)=x_2(0)=0$. Let $u(t)=1$ signal unit from $t=0$ to $t=1$ minute, then $u(t)=0$. Find the output $x_2(2)$.

Direct ODE solution. During the first minute, the input is constant at 1. Solving the first equation gives $x_1(t)=1-e^{-t}$. Substitution into the second equation gives $x_2(t)=1-(1+t)e^{-t}$. Thus, at the end of the pulse, $$ x_1(1)=1-e^{-1},\qquad x_2(1)=1-2e^{-1}. $$ After withdrawal, $x_1$ decays exponentially. Solving the second-stage equation over the next minute gives $$ x_2(2)=e^{-1}\bigl[x_2(1)+x_1(1)\bigr] =e^{-1}\bigl[(1-2e^{-1})+(1-e^{-1})\bigr] =2e^{-1}-3e^{-2}=0.3298. $$ The final value is rounded only after evaluating the exact expression.

Convolution check. The two-stage impulse response is $h_2(t)=te^{-t}$. The pulse contributes at $t=2$ only from input times $s$ between 0 and 1: $$ x_2(2)=\int_0^1 (2-s)e^{-(2-s)}\,ds =\int_1^2 a e^{-a}\,da =2e^{-1}-3e^{-2}=0.3298. $$ Both calculations agree. Here $x_2$ is in signal units, and time is measured in minutes.

Repeated Doses and Superposition

For a linear chain, the response to a sum of inputs is the sum of the responses to each input. This follows directly from the convolution integral: if $u(t)=u_1(t)+u_2(t)$, then $$ \int_0^t h_m(t-s)\bigl[u_1(s)+u_2(s)\bigr]\,ds =\int_0^t h_m(t-s)u_1(s)\,ds+\int_0^t h_m(t-s)u_2(s)\,ds. $$ Repeated-dose input profiles can therefore be handled by summing the contributions from individual dose-related inputs. This is a statement about the linear cascade receiving those inputs. It does not imply that concentration profiles themselves add linearly after nonlinear pharmacokinetics or a nonlinear concentration-to-signal transformation.

Worked Example: Two Repeated Input Pulses

Use the same two-stage chain with $k=1$ per minute and zero initial states. Let $u(t)=1$ during the intervals from 0 to 1 minute and from 2 to 3 minutes, and zero otherwise. Calculate $x_2(3)$.

The first pulse reaches the output at time 3 with ages from 2 to 3 minutes. Its contribution is $$ \int_0^1 (3-s)e^{-(3-s)}\,ds =\int_2^3 a e^{-a}\,da =3e^{-2}-4e^{-3}=0.2069. $$ The second pulse contributes at ages from 0 to 1 minute: $$ \int_2^3 (3-s)e^{-(3-s)}\,ds \quad\text{with the pulse’s input interval represented by }s\in[2,3], $$ or, using age $a=3-s$, $$ \int_0^1 a e^{-a}\,da=1-2e^{-1}=0.2642. $$ Adding the two contributions gives $$ x_2(3)=3e^{-2}-4e^{-3}+1-2e^{-1}=0.4711. $$ The output is in signal units. The first pulse still contributes at the time of the second pulse’s end, demonstrating that the current output reflects the full input history.

The same calculation works when each dose produces a different input shape: integrate the impulse response against that dose’s time-varying contribution. For many doses, direct summation is often convenient when the single-dose input profile is known. If doses are given at regular intervals and the system has reached a repeating pattern, the contribution from earlier doses may also be included through a repeated-dose sum. A finite simulation should not silently omit doses before its chosen start time if they leave a non-negligible signal in the chain.

Numerically Integrating a Changing Input

Direct ODE integration is useful when the input is available as a time series, the rates differ across stages, or the model includes features that prevent a simple convolution calculation. A numerical solver advances the states using the right-hand sides of the ODEs, evaluating $u(t)$ at the solver’s required times. A convolution calculation instead evaluates an integral using $h_m$ and the full input history. For a linear chain, the two are mathematically equivalent when they use the same input and initial conditions.

Worked Example: Exponentially Changing Input

Let two equal-rate stages have $k=1$ per minute, zero initial states, and $u(t)=e^{-t}$ signal units. Calculate $x_2(2)$ using convolution, then consider numerical integration.

The convolution integrand is $$ h_2(2-s)u(s)=(2-s)e^{-(2-s)}e^{-s}=(2-s)e^{-2}. $$ Therefore, $$ x_2(2)=\int_0^2(2-s)e^{-2}\,ds =e^{-2}\left[2s-\frac{s^2}{2}\right]_0^2 =2e^{-2}=0.2707 $$ signal units. The direct ODE solution provides a check: $x_1(t)=te^{-t}$, and then $x_2(t)=\tfrac{1}{2}t^2e^{-t}$, which also gives $x_2(2)=2e^{-2}$.

For a numerical convolution, midpoint quadrature with four intervals of width $\Delta t=0.5$ evaluates the integrand at $s=0.25,0.75,1.25,1.75$: $$ x_2(2)\approx 0.5e^{-2}(1.75+1.25+0.75+0.25)=2e^{-2}=0.2707. $$ This particular quadrature is exact because the integrand is linear in $s$; that exactness should not be expected for a general input.

A direct ODE solver approximates the evolving states rather than this integral. For instance, explicit Euler with a coarse step of 0.5 minute gives $x_1(0.5)=0.5$, $x_2(0.5)=0$, and continues using the current input and states at each step. That step size is too coarse for an accurate final value here; the result should be checked by reducing the step size or using an adaptive solver with tighter tolerances. Agreement with the convolution result is a useful verification when both methods are available.

Choosing and Checking a Calculation

Approach	Strength	Check or limitation
Direct ODE integration	Handles changing inputs and general cascade equations flexibly.	Check solver tolerances, step-size sensitivity, and input interpolation.
Convolution	Shows how the entire input history contributes; efficient when the impulse response is available.	Requires a linear, time-invariant cascade and correct treatment of initial states.
Closed-form segments	Can be exact for simple inputs such as constant or exponential segments.	Algebra becomes cumbersome for many segments or complex dosing schedules.

A dose may create a discontinuity in $u(t)$, even when the concentration profile itself is continuous. A numerical solver should be informed of known discontinuity times, or the integration should be split at those times. Otherwise, a solver may take a step across an abrupt change and represent the input poorly. For concentration measurements or model predictions supplied on a grid, the interpolation rule matters too: a stepwise-constant input and a linearly interpolated input define different $u(t)$ between observations.

Useful checks include confirming that the input and output have compatible signal units, that rates have inverse-time units, and that a constant input eventually produces the same steady state at every stage. Numerical results should be stable to reasonable refinement of step size and solver tolerances. A negative state under a nonnegative input can indicate numerical error or an inappropriate model formulation; it should not be dismissed without checking the equations and the integration settings.

Assumptions, Limitations, and a Common Pitfall

Convolution and dose-by-dose superposition rely on constant rates and linear transfer within the chain. If transfer between stages is nonlinear, the output from two inputs generally is not the sum of the separate outputs. Direct numerical solution may still be possible, but the linear impulse-response convolution does not represent that nonlinear system. Likewise, if rates change with time, the response depends on when the signal entered the chain, not just on its age, so a single time-invariant impulse response is not sufficient.

Common mistake: “The current output equals the current input because the cascade’s steady-state output equals its steady-state input.” This confuses an equilibrium relationship with a time-dependent response. A cascade can have equal input and output at steady state while still showing substantial delay after a changing input. The correction is to solve the ODEs over time or use the convolution integral, with the input history and initial states specified.

Key takeaway: For a linear transit chain with fixed rates, a time-varying input can be propagated either by integrating the cascade ODEs or by convolving the input history with the chain’s impulse response. Correct handling of input timing, initial states, and numerical resolution is essential for the two calculations to agree.

Check Your Understanding

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

For a two-stage chain with equal rate $k$, what is the impulse response, and what does its time argument represent?
Why can the responses to repeated input pulses be added in a linear cascade?
A solver’s output changes substantially when its maximum step size is halved. What should be checked before interpreting the output?
Does a convolution with one fixed impulse response generally apply if a transit rate changes over time? Explain.
For two equal-rate stages with zero initial states and constant input $u(t)=u_0$, what steady-state value does each stage approach?

Solutions. (1) $h_2(t)=k^2te^{-kt}$; $t$ is the elapsed time since an impulse entered the first stage. (2) The convolution integral is linear in its input, so the integral of a sum is the sum of the separate integrals. (3) Check step-size sensitivity, error tolerances, input interpolation, and whether the solver steps across a known discontinuity; refine or split the integration until results stabilize. (4) Not generally: with time-varying rates, the response depends on the input time as well as the elapsed time, so one fixed time-invariant impulse response is inadequate. (5) Each stage approaches $u_0$, because at steady state every derivative is zero, requiring $x_1=u_0$ and $x_j=x_{j-1}$ for later stages.