Three-Compartment Tissue Distribution With Stochastic Exchange

Why Model Tissue Exchange Stochastically?

A three-compartment model distinguishes drug in a central compartment from drug in shallow and deep tissue. In a deterministic model, pairwise exchange is driven by fixed rate constants. An SDE can represent additional random variation in the exchange itself. The key modeling decision is to represent each exchange by a shared stochastic flux: the amount randomly transferred out of one compartment is the same amount transferred into its partner.

This is different from adding independent noise to each compartment. Independent noise can create or destroy total drug even when the system is meant to describe redistribution alone. The model below isolates stochastic redistribution; it has no dose input, elimination, or measurement error. Those processes would require additional terms and separate mass-balance accounting.

Pairwise Exchange SDE

Let $A_c(t)$, $A_s(t)$, and $A_d(t)$ denote amounts of drug in the central, shallow-tissue, and deep-tissue compartments. For each connected pair $i,j$, define $k_{ij}$ as the first-order rate constant for transfer from $i$ to $j$. Thus, $k_{cs}A_c$ is the central-to-shallow transfer rate, while $k_{sc}A_s$ is the reverse rate. All rate constants are nonnegative and have units of inverse time.

To make the stochastic flux vanish when either compartment has no drug, choose a scale $K>0$, measured in amount units, and define $h(x,y)=xy/(K+x+y)$ for nonnegative amounts $x,y$. For each pair, use a Brownian motion $W_{ij}$ and a nonnegative noise scale $\sigma_{ij}$, with units of inverse square-root time. The same $W_{ij}$ must appear with opposite signs in the two compartments connected by that pair.

Definition: A closed three-compartment stochastic exchange model is $$ \begin{aligned} dA_c={}&[-(k_{cs}A_c-k_{sc}A_s)-(k_{cd}A_c-k_{dc}A_d)]\,dt\\ &-\sigma_{cs}h(A_c,A_s)\,dW_{cs} -\sigma_{cd}h(A_c,A_d)\,dW_{cd},\\ dA_s={}&[(k_{cs}A_c-k_{sc}A_s)-(k_{sd}A_s-k_{ds}A_d)]\,dt\\ &+\sigma_{cs}h(A_c,A_s)\,dW_{cs} -\sigma_{sd}h(A_s,A_d)\,dW_{sd},\\ dA_d={}&[(k_{cd}A_c-k_{dc}A_d)+(k_{sd}A_s-k_{ds}A_d)]\,dt\\ &+\sigma_{cd}h(A_c,A_d)\,dW_{cd} +\sigma_{sd}h(A_s,A_d)\,dW_{sd}. \end{aligned} $$ Here $W_{cs},W_{cd},W_{sd}$ are independent standard Brownian motions. The initial amounts are nonnegative. The function $h$ has amount units, so each stochastic term has amount units, as required.

The drift for a pair is its net deterministic flux: for example, $k_{cs}A_c-k_{sc}A_s$ is positive when the net drift is from central to shallow tissue. The stochastic term for that same pair is a random flux whose sign is not fixed in advance. Its direction on any small time interval depends on the Brownian increment.

Worked Example: Drift and Instantaneous Variability

Take $K=50$ mg and current amounts $A_c=100$ mg, $A_s=20$ mg, and $A_d=10$ mg. Let the rates, in $\text{h}^{-1}$, be $k_{cs}=0.20$, $k_{sc}=0.10$, $k_{cd}=0.05$, $k_{dc}=0.02$, $k_{sd}=0.08$, and $k_{ds}=0.04$.

The three net deterministic fluxes are $q_{cs}=0.20(100)-0.10(20)=18$ mg/h, $q_{cd}=0.05(100)-0.02(10)=4.8$ mg/h, and $q_{sd}=0.08(20)-0.04(10)=1.2$ mg/h. Thus the drift rates are $-q_{cs}-q_{cd}=-22.8$ mg/h centrally, $q_{cs}-q_{sd}=16.8$ mg/h in shallow tissue, and $q_{cd}+q_{sd}=6.0$ mg/h in deep tissue. Their sum is zero.

Set $\sigma_{cs}=0.30$, $\sigma_{cd}=0.20$, and $\sigma_{sd}=0.25$, in $\text{h}^{-1/2}$. The exchange functions are $h(100,20)=2000/170=11.7647$ mg, $h(100,10)=1000/160=6.25$ mg, and $h(20,10)=200/80=2.5$ mg. The corresponding stochastic amplitudes are $3.5294$, $1.25$, and $0.625$ mg/$\sqrt{\text{h}}$. Since the two Brownian motions affecting $A_c$ are independent, its instantaneous variance rate is $3.5294^2+1.25^2=12.4567+1.5625=14.0192$ mg$^2$/h.

Mass Balance and Mean Behavior

The pairwise construction gives an exact conservation result, not merely an approximate one. This requires that every pairwise drift flux and Brownian flux occur with opposite signs in its two compartments. The result holds for all nonnegative rate and noise parameters and for any solution of the displayed system.

Theorem: If $T(t)=A_c(t)+A_s(t)+A_d(t)$, then $T(t)=T(0)$ almost surely for every $t\geq0$. If the initial amounts have finite second moments, the model has finite second moments on finite time intervals, and the mean amounts satisfy the deterministic exchange equations obtained by removing the Brownian terms.

Proof. Add the three SDEs. The central–shallow drift terms sum to $-q_{cs}+q_{cs}=0$, the central–deep terms to $-q_{cd}+q_{cd}=0$, and the shallow–deep terms to $-q_{sd}+q_{sd}=0$. Each Brownian flux cancels in the same way; for example, the $W_{cs}$ terms are $-\sigma_{cs}h(A_c,A_s)dW_{cs}$ and $+\sigma_{cs}h(A_c,A_s)dW_{cs}$. Therefore $dT=0$, which implies $T(t)=T(0)$ almost surely.

For the mean statement, $h(x,y)$ is locally Lipschitz on the nonnegative orthant and has at most linear growth there. The drift is linear and the diffusion coefficients also have at most linear growth. The standard existence and moment bounds for SDEs with these conditions give finite second moments on finite time intervals. Consequently, each Brownian integral has mean zero: its integrand is square-integrable over such an interval. Taking expectations in the integral form of the SDE is therefore valid. Linearity of the drift gives, for example, $\frac{d}{dt}E[A_c]=-(k_{cs}+k_{cd})E[A_c]+k_{sc}E[A_s]+k_{dc}E[A_d]$; the other two mean equations follow in the same way. The stochastic terms affect variability and joint dependence, but not these mean equations. $\square$

For nonnegative initial amounts, the continuous-time model also has an inward boundary structure. If, for example, $A_c=0$, both stochastic coefficients involving $A_c$ vanish because $h(0,y)=0$; the central drift is $k_{sc}A_s+k_{dc}A_d$, which is nonnegative. The analogous property holds at the other two boundaries. With the stated locally Lipschitz coefficients and nonnegative rates, the standard boundary criterion for SDEs therefore keeps the solution in the nonnegative orthant. This property of the continuous-time model does not guarantee that every numerical time step will stay nonnegative.

Worked Example: A Short-Time Mean Prediction

Use the amounts and rates from the first example, and consider a step of $\Delta t=0.10$ h. The initial mean drift vector is $(-22.8,16.8,6.0)$ mg/h. A forward Euler approximation to the mean over this short step is

\begin{aligned} E[A_c(\Delta t)]&\approx100+0.10(-22.8)=97.72\ \text{mg},\\ E[A_s(\Delta t)]&\approx20+0.10(16.8)=21.68\ \text{mg},\\ E[A_d(\Delta t)]&\approx10+0.10(6.0)=10.60\ \text{mg}. \end{aligned}

The approximate total is $97.72+21.68+10.60=130.00$ mg, matching the initial total. This is a first-order Euler prediction, not an exact solution of the mean ODE over 0.10 h: the drift changes as the expected amounts change. The stochastic exchange changes the distribution around the mean, even though it does not add a term to these mean equations.

Reading a Stochastic Redistribution Step

For each pair, a positive Brownian increment in the chosen orientation produces a transfer in that orientation: from central to shallow for $W_{cs}$, central to deep for $W_{cd}$, and shallow to deep for $W_{sd}$. A negative increment reverses the direction. In either case, the paired compartments receive equal and opposite stochastic increments.

Worked Example: One Euler–Maruyama Redistribution Step

Start again from $(100,20,10)$ mg with the first example’s parameters. Take $\Delta t=0.01$ h and illustrative Brownian increments $\Delta W_{cs}=0.04$, $\Delta W_{cd}=-0.02$, and $\Delta W_{sd}=0.03$, each in $\sqrt{\text{h}}$. Freeze the coefficients at the start of the step, as in Euler–Maruyama. The three pairwise stochastic flux increments are $\xi_{cs}=3.5294(0.04)=0.141176$ mg, $\xi_{cd}=1.25(-0.02)=-0.025$ mg, and $\xi_{sd}=0.625(0.03)=0.01875$ mg.

The drift increments are $(-22.8,16.8,6.0)(0.01)=(-0.228,0.168,0.060)$ mg. Combining drift and paired fluxes gives

\begin{aligned} \Delta A_c&=-0.228-\xi_{cs}-\xi_{cd}=-0.344176\ \text{mg},\\ \Delta A_s&=0.168+\xi_{cs}-\xi_{sd}=0.290426\ \text{mg},\\ \Delta A_d&=0.060+\xi_{cd}+\xi_{sd}=0.053750\ \text{mg}. \end{aligned}

The resulting Euler–Maruyama amounts are approximately $(99.655824,20.290426,10.053750)$ mg. The increments sum to $-0.344176+0.290426+0.053750=0$, and the amounts sum to 130 mg. The numerical step illustrates the pairwise accounting; because it is an approximation, it does not by itself establish that a simulated step will always remain nonnegative.

Diagnostics, Alternatives, and Limitations

A useful first diagnostic is to calculate total amount at every simulated time point. In a closed exchange-only model, any change signals an implementation error, such as using separate Brownian increments for the two sides of a pair or assigning the same sign to both sides. Also check that the mean trajectories agree with the deterministic mean equations and that simulated variability increases or decreases sensibly as the $\sigma_{ij}$ values change.

A deterministic compartment ODE is the simpler alternative when exchange is adequately described by fixed fluxes and trajectory-level randomness is not needed. Adding independent noise to each compartment is not an equivalent alternative: it changes total amount and describes a different process. Correlated multivariate noise can be used instead, but its covariance structure must encode the same equal-and-opposite pairwise fluxes if redistribution alone is intended.

The chosen function $h$ is a modeling choice, not a universal biological law. It makes exchange variability vanish when either connected compartment is empty, but its scale $K$ and noise strengths need scientific justification and data support. The closed model also excludes input and elimination; including either requires adding its flux to the relevant compartment and accounting for its effect on total amount. Finally, Euler–Maruyama updates can cross a boundary even when the continuous-time model has an inward boundary structure. Monitor negative simulated amounts and step-size sensitivity; simply clipping negative values can introduce bias and conceal a numerical or model problem.

Common Pitfall

Correction: “I added independent noise to each amount, so the randomness averages to zero and total drug is still conserved.” Zero-mean noise does not conserve total amount pathwise. If $dA_c$ and $dA_s$ have independent Brownian terms, their sum generally contains a nonzero random increment. For a pairwise exchange, use the same Brownian motion with opposite signs; then the increments cancel exactly in the total.

Check Your Understanding

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

For the central–shallow pair, write the drift contributions to $dA_c$ and $dA_s$, and explain why their sum is zero.
If the shallow–deep Brownian increment is negative, which direction does its stochastic flux move drug under the stated orientation?
Suppose $A_c=0$ while $A_s,A_d\geq0$. What are the central stochastic coefficients and the central drift?
Why do the stochastic terms not appear in the mean equations under the stated moment conditions?
A simulation’s total amount rises on one step even though the model has no input or elimination. Give one likely implementation error.

Solutions. (1) The contributions are $-(k_{cs}A_c-k_{sc}A_s)dt$ centrally and $+(k_{cs}A_c-k_{sc}A_s)dt$ in shallow tissue; they cancel. (2) The negative increment reverses the chosen shallow-to-deep orientation, so the stochastic flux moves drug from deep to shallow. (3) Both central diffusion coefficients vanish because $h(0,y)=0$; the central drift is $k_{sc}A_s+k_{dc}A_d$, which is nonnegative. (4) The Brownian integrals have mean zero because their integrands are square-integrable on finite intervals. (5) A likely error is using different Brownian increments or the same sign for the two sides of a pairwise flux.