Tutorials › Stochastic Differential Equations in Pharmacokinetics › Worked Construction Of A Two-Compartment SDE

Writing a Pharmacokinetic SDE · Tutorial 34 of 100

Worked Construction Of A Two-Compartment SDE

Build and check a two-compartment amount SDE by assigning each drift and noise term to its physical process and tracing its contribution to total drug amount.

Intermediate 10 min read

What You'll Learn

  • Write the central and peripheral drift terms from their transfer and elimination rates.
  • Represent a fluctuating intercompartmental transfer with one shared Brownian increment and opposite signs.
  • Distinguish a transfer-noise channel from a central-elimination-noise channel.
  • Derive the total-amount SDE and identify which mechanisms cancel from it.
  • Check noise units, covariation, mass balance, and physical limitations of the model.

From Compartment Diagram to Stochastic Model

A two-compartment pharmacokinetic model has two amount states: \(A_c(t)\), the amount in the central compartment, and \(A_p(t)\), the amount in the peripheral compartment. Deterministic intercompartmental transfer moves drug between these states, while elimination removes drug from the body. A stochastic model must also say which physical process each random increment represents. This matters for mass balance: random variation in transfer redistributes amount, whereas random variation in elimination changes the total amount.

We construct a model between dose events. Let \(k_{cp}>0\) be the first-order transfer rate from central to peripheral, \(k_{pc}>0\) the rate from peripheral to central, and \(k_e>0\) the central elimination rate. Let \(u(t)\) be a deterministic input rate, such as a bounded infusion rate, into the central compartment. A bolus dose is instead represented as a jump in amount; it is not included in these continuous-time equations.

Definition: Let \(W_T(t)\) and \(W_E(t)\) be independent standard Wiener processes. The transfer-noise channel and elimination-noise channel are distinct, independent sources of continuous process noise. For constants \(\sigma_T,\sigma_E\geq 0\), define the amount model $$ \begin{aligned} dA_c(t)={}&\big[u(t)-(k_{cp}+k_e)A_c(t)+k_{pc}A_p(t)\big]\,dt\\ &-\sigma_T A_c(t)\,dW_T(t)-\sigma_E A_c(t)\,dW_E(t),\\ dA_p(t)={}&\big[k_{cp}A_c(t)-k_{pc}A_p(t)\big]\,dt +\sigma_T A_c(t)\,dW_T(t). \end{aligned} $$ Here \(A_c,A_p\) are measured in amount, the rate constants in inverse time, \(u\) in amount per time, and \(\sigma_T,\sigma_E\) in inverse square-root time. The initial state is assumed square-integrable and independent of the future Brownian increments. The input \(u\) is deterministic, measurable, and bounded on each finite time interval.

The equations use one \(W_T\) increment with opposite signs. It represents a fluctuating central-to-peripheral transfer flux: the same random amount lost from the central compartment is gained by the peripheral compartment. Its amplitude \(\sigma_T A_c\) is a modeling choice, not a universal law. The elimination channel acts only on the central amount, with negative sign, because elimination leaves the modeled body system. Since the coefficients are linear in the states, they are globally Lipschitz; under the stated input and initial-state assumptions, the standard existence-and-uniqueness theorem for Itô SDEs gives a unique global strong solution.

Assembling and Checking Each Term

The deterministic central drift contains three contributions: central-to-peripheral transfer removes \(k_{cp}A_c\), elimination removes \(k_eA_c\), and return transfer adds \(k_{pc}A_p\). The peripheral drift receives \(k_{cp}A_c\) and loses \(k_{pc}A_p\). The input \(u(t)\) enters only the central compartment. This accounting gives the drift terms before any stochastic terms are added.

For the noise, first identify the mechanism and then its destinations. Transfer noise belongs to both compartments, with opposite signs; elimination noise belongs only to the central compartment. Distinct drivers are used because the two random mechanisms are modeled as independent. The coefficients must also have the right units: multiplying an amount by \(\sigma\,dW\) gives an amount, since \(dW\) has units of square-root time.

Worked Example: Writing the Instantaneous Model

Suppose \(k_{cp}=0.30\ \mathrm{h}^{-1}\), \(k_{pc}=0.20\ \mathrm{h}^{-1}\), \(k_e=0.10\ \mathrm{h}^{-1}\), \(\sigma_T=0.15\ \mathrm{h}^{-1/2}\), and \(\sigma_E=0.08\ \mathrm{h}^{-1/2}\). At a particular time, let \(A_c=100\) mg, \(A_p=40\) mg, and \(u=0\). The central drift is \( -(0.30+0.10)(100)+(0.20)(40)=-40+8=-32\) mg/h. The peripheral drift is \((0.30)(100)-(0.20)(40)=30-8=22\) mg/h.

At these state values, the central diffusion coefficients are \(-15\) mg/\(\sqrt{\mathrm h}\) for \(W_T\) and \(-8\) mg/\(\sqrt{\mathrm h}\) for \(W_E\). The peripheral coefficient for \(W_T\) is \(+15\) mg/\(\sqrt{\mathrm h}\). Thus the model at this instant is $$ \begin{aligned} dA_c&=(-32\ \mathrm{mg/h})\,dt -(15\ \mathrm{mg}/\sqrt{\mathrm h})\,dW_T -(8\ \mathrm{mg}/\sqrt{\mathrm h})\,dW_E,\\ dA_p&=(22\ \mathrm{mg/h})\,dt +(15\ \mathrm{mg}/\sqrt{\mathrm h})\,dW_T. \end{aligned} $$ The transfer drift sums to \(-30+30=0\) mg/h, and the transfer-noise coefficients sum to \(-15+15=0\) mg/\(\sqrt{\mathrm h}\). The remaining total drift is \(-10\) mg/h, exactly the central elimination rate at this state.

Mass Balance: A Pathwise Result

Define total amount \(A_{\mathrm{tot}}(t)=A_c(t)+A_p(t)\). The balance equation follows by adding the two SDEs. Addition is valid directly for Itô differentials: there is no nonlinear transformation here and hence no second-derivative correction.

Theorem: For the model above, on every finite time interval and almost surely, $$ dA_{\mathrm{tot}}(t)=\big[u(t)-k_eA_c(t)\big]\,dt-\sigma_EA_c(t)\,dW_E(t). $$ In particular, both deterministic intercompartmental exchange and transfer noise cancel from the total-amount equation.
1
Add the drift terms.
The sum of the two drifts is \(\,u-(k_{cp}+k_e)A_c+k_{pc}A_p+k_{cp}A_c-k_{pc}A_p=u-k_eA_c.\) Each transfer term has an equal and opposite contribution.
2
Add the Brownian terms by driver.
The coefficient of \(dW_T\) is \(-\sigma_TA_c+\sigma_TA_c=0\). The coefficient of \(dW_E\) is \(-\sigma_EA_c\), because elimination noise appears only in the central equation.
3
Integrate over time.
For any \(0\leq s\leq t\), linearity of the time and Itô integrals gives $$ A_{\mathrm{tot}}(t)-A_{\mathrm{tot}}(s) =\int_s^t [u(r)-k_eA_c(r)]\,dr -\int_s^t \sigma_EA_c(r)\,dW_E(r). $$ This is the asserted balance, including its remaining elimination fluctuation. \(\square\)

The theorem is a pathwise accounting identity, not a claim that the amount decreases monotonically on every sample path. The Brownian term can have either sign over a time interval. It expresses uncertainty in the continuous elimination process around its specified drift; it is not a literal guarantee that every infinitesimal stochastic increment removes a nonnegative amount.

Worked Example: Checking an Infusion-Step Balance

Use the parameters and state values from the first example, but now set the infusion rate to \(u=5\) mg/h. The central drift is \(-32+5=-27\) mg/h and the peripheral drift remains \(22\) mg/h, so the total drift is \(-5\) mg/h. For a left-endpoint Euler step of \(\Delta t=0.1\) h, take \(\Delta W_T=0.4\sqrt{\mathrm h}\) and \(\Delta W_E=-0.25\sqrt{\mathrm h}\). The central increment is $$ \Delta A_c\approx(-27)(0.1)-(15)(0.4)-(8)(-0.25) =-2.7-6+2=-6.7\ \mathrm{mg}. $$ The peripheral increment is $$ \Delta A_p\approx(22)(0.1)+(15)(0.4)=2.2+6=8.2\ \mathrm{mg}. $$ Their sum is \(1.5\) mg. The balance equation gives the same result: the deterministic total increment is \((-5)(0.1)=-0.5\) mg, while the elimination-noise increment is \(-8(-0.25)=2\) mg. Thus \(-0.5+2=1.5\) mg. The transfer increment cancels exactly in the total.

Why the Shared Transfer Driver Matters

Using opposite signs is not enough if the two compartments are assigned independent Brownian drivers. The same random transfer event must be represented by the same increment in both equations. This can also be checked through the instantaneous covariance contributed by the transfer channel.

Result: At a fixed state, the transfer-noise contribution to the covariance matrix per unit time for \((A_c,A_p)\) is $$ \sigma_T^2 A_c^2 \begin{pmatrix} 1&-1\\ -1&1 \end{pmatrix}. $$ Its row sums are zero, reflecting that this noise redistributes amount without changing the total.

To verify the result, the transfer diffusion vector is \(g_T=(-\sigma_TA_c,\ \sigma_TA_c)^\mathsf{T}\). A single scalar Brownian driver contributes the outer product \(g_Tg_T^\mathsf{T}\) per unit time. Multiplication gives the displayed matrix. In particular, the transfer-induced instantaneous variances of the two states are both \(\sigma_T^2A_c^2\), and their covariance is \(-\sigma_T^2A_c^2\). The negative covariance is the mathematical signature of opposite changes from one shared transfer fluctuation.

Worked Example: Detecting Independent-Noise Misplacement

Suppose the transfer terms are incorrectly written as \(-15\,dW_1\) in the central equation and \(+15\,dW_2\) in the peripheral equation, where \(W_1\) and \(W_2\) are independent. The resulting transfer contribution to total amount is \(-15\,dW_1+15\,dW_2\), not zero. Over a time interval of length \(\Delta t\), its variance is $$ 225\Delta t+225\Delta t=450\Delta t\ \mathrm{mg}^2, $$ because the increments are independent and each has variance \(\Delta t\). This invents random changes in total amount from an internal transfer.

The correction is to use one driver, \(-15\,dW_T\) centrally and \(+15\,dW_T\) peripherally. Their sum is identically zero, not merely zero on average. If instead the goal is to model independent compartment-specific environmental disturbances, those can be included as additional channels, but they should not be labelled a mass-conserving transfer.

Assumptions, Diagnostics, and Limitations

The construction depends on a clear definition of the system boundary: elimination leaves the two-compartment body system, while intercompartmental transfer stays inside it. If a metabolite or another tissue compartment is added to the system, the accounting boundary and equations must be reconsidered. Likewise, an infusion contributes positively to total amount, whereas an external dose bolus produces a jump in total amount at its dosing time.

Several checks help identify a model that has been assembled incorrectly:

  • Sum the drift coefficients. Exchange should cancel from total amount; central elimination should remain as \(-k_eA_c\).
  • Group diffusion terms by Brownian driver. A transfer channel should have coefficients that sum to zero across the participating compartments.
  • Check dimensions. Drift coefficients must yield amount per time, and diffusion coefficients must yield amount per square-root time.
  • Inspect simulated paths. A transfer-only simulation should not change total amount except for numerical integration error. A persistent total change suggests incorrect signs, mismatched transfer amplitudes, or mistakenly independent drivers.
  • Check state plausibility. This linear-noise model does not guarantee nonnegative amounts. Negative simulated values may indicate that the chosen diffusion structure or numerical step is unsuitable for the intended application.

The last limitation is important. A Brownian fluctuation proportional to \(A_c\) can move an amount in either direction, and an SDE with this structure does not automatically preserve the nonnegative state space. Clipping negative values after simulation changes the model and can distort its distribution. If strict nonnegativity is essential, a positivity-preserving model or method must be chosen and justified; that is a different modeling decision, not a consequence of the mass-balance calculation.

The main alternatives answer different questions. A deterministic two-compartment ODE has no continuous process noise and gives a deterministic mass balance. Independent noise in each state may represent separate disturbances, but does not represent a conserved internal transfer. A shared transfer channel, as used here, represents uncertain redistribution while preserving its exact contribution to total amount. The elimination channel adds a separate source of uncertainty to total amount. Neither channel should be included just to make simulated trajectories look variable: its interpretation and placement are part of the model specification.

Key takeaway: Write a transfer fluctuation once as a shared Brownian increment with equal-and-opposite compartment coefficients. When the compartment equations are added, that channel cancels; only external input and processes that cross the system boundary remain in the total-amount balance.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. With no input, write the total-amount drift for the model and explain why neither exchange rate appears.
  2. If the transfer diffusion coefficient is \(q(A_c,A_p)\) rather than \(\sigma_TA_c\), what signs and shared-driver structure preserve mass balance?
  3. For a transfer-only interval, what is the transfer contribution to the variance of \(A_c+A_p\) when the same Brownian increment is used with opposite coefficients?
  4. What changes in the total balance if a deterministic infusion \(u(t)\) is added to the peripheral compartment instead of the central compartment?

Solutions. (1) The drift is \(-k_eA_c\): the central-to-peripheral loss and gain are \(-k_{cp}A_c\) and \(+k_{cp}A_c\), and the reverse-transfer terms \(+k_{pc}A_p\) and \(-k_{pc}A_p\) also cancel. (2) Use \(-q(A_c,A_p)dW_T\) centrally and \(+q(A_c,A_p)dW_T\) peripherally, with the same driver; their sum is zero. (3) It is zero, since the total transfer increment is \(-q\,dW_T+q\,dW_T=0\) for every increment. (4) The total balance is unchanged except that the same input \(u(t)\) is still added to total amount; the input’s compartment location changes the separate state equations, not the total contribution.

Clinical Biostats