Numerical Stability at Extreme Drug Effects

Why Extreme Drug Effects Challenge Computation

Indirect-response models describe a response that changes because drug alters the rate of production, the rate of loss, or both. When production is almost completely inhibited, the production rate can be very small. When loss is strongly stimulated, the response may change on a very short time scale. These are valid model behaviors, but they can expose weaknesses in both the algebra used to calculate drug effects and the numerical method used to solve the differential equation.

A computed curve can look smooth and still be wrong. It may have lost a small but important residual production rate through subtraction, become negative because an integration step was too large, or appear to reach equilibrium only because the solver stopped resolving a rapid transient. Numerical stability is therefore not just a programming concern: it affects whether an estimated model can be interpreted biologically.

Definition: Consider the indirect-response model $$ \frac{dR}{dt}=k_{\mathrm{in}}q(C(t))-k_{\mathrm{out}}g(C(t))R(t), $$ where $R(t)$ is the response, $C(t)$ is drug concentration, $k_{\mathrm{in}}>0$ is the baseline production rate, $k_{\mathrm{out}}>0$ is the baseline loss-rate constant, $q(C)\geq 0$ is the production multiplier, and $g(C)>0$ is the loss multiplier. With no drug effect, $q(0)=g(0)=1$, and the baseline response is $R_0=k_{\mathrm{in}}/k_{\mathrm{out}}$.

This tutorial focuses on computation near extreme values of $q(C)$ or $g(C)$, not on deciding which mechanism is biologically appropriate. The earlier material on indirect-response mechanisms and parameter identifiability supplies that modeling context. The numerical checks below assume that concentration and parameters are finite, the response begins nonnegative, and the rate functions are nonnegative over the concentrations being evaluated.

Stable Calculation of Nearly Complete Inhibition

For an inhibitory production effect with a Hill coefficient $\gamma>0$, a common form is $$ q(C)=1-I_{\max}\frac{C^\gamma}{IC_{50}^{\gamma}+C^\gamma}, $$ where $0\leq I_{\max}\leq 1$. When $I_{\max}$ is close to one and $C$ is large relative to $IC_{50}$, the two terms in the subtraction are nearly equal. In finite-precision arithmetic, subtracting nearly equal numbers can discard significant digits. The result may be an inaccurate small positive value, or even zero when the exact value is positive.

An algebraic rearrangement exposes the residual production directly. Define $\varepsilon=1-I_{\max}$, the fraction of baseline production remaining at arbitrarily high concentration when the effect is saturable. With $x=(C/IC_{50})^\gamma$, $$ q(C)=\varepsilon+\frac{1-\varepsilon}{1+x}. $$ This form adds nonnegative quantities rather than subtracting two nearly equal ones. If complete inhibition is allowed, $\varepsilon=0$; at any finite $x$, the residual $1/(1+x)$ remains positive.

The ratio $x$ itself can overflow when $C/IC_{50}$ is very large or $\gamma$ is high. One stable alternative is to calculate $\ell=\gamma[\log(C)-\log(IC_{50})]$, so $x=\exp(\ell)$, and evaluate the fraction $1/(1+\exp(\ell))$ using a stable logistic-function implementation. For large positive $\ell$, compute it as $\exp(-\ell)/(1+\exp(-\ell))$, avoiding an unnecessarily large exponential. At $C=0$, use the model's limiting value $q(0)=1$ rather than taking $\log(0)$.

Worked Example: Preserving a Small Residual Production Rate

Suppose a hypothetical inhibitory model has $I_{\max}=1-10^{-12}$, Hill coefficient $\gamma=1$, and $C/IC_{50}=10^8$. Calculate the production multiplier in a form that retains its small residual.

Here $\varepsilon=10^{-12}$ and $x=(10^8)^1=10^8$. The rearranged expression gives $$ q(C)=10^{-12}+\frac{1-10^{-12}}{1+10^8} \approx 10^{-12}+9.9999999\times 10^{-9} \approx 1.0001\times 10^{-8}. $$

The multiplier is small, but it is not zero. If a calculation instead rounded the inhibitory fraction to exactly one before subtracting it from one, the residual contribution from the concentration-dependent term would be lost. The stable expression makes clear that production remains positive at this finite concentration, as required by the stated model.

Exact Constant-Exposure Solution and Its Numerical Value

When concentration is constant over an interval, $q(C)$ and $g(C)$ are constant there as well. Write $q=q(C)$, $g=g(C)$, and $\lambda=k_{\mathrm{out}}g$. The differential equation becomes $dR/dt=k_{\mathrm{in}}q-\lambda R$. Its equilibrium under that exposure is $R_{\mathrm{eq}}=k_{\mathrm{in}}q/\lambda$, and the exact solution over an interval of length $t$ is $$ R(t)=R_{\mathrm{eq}}+\{R(0)-R_{\mathrm{eq}}\}\exp(-\lambda t). $$

A useful equivalent form is $$ R(t)=R(0)+\{R_{\mathrm{eq}}-R(0)\}\left[1-\exp(-\lambda t)\right]. $$ When $\lambda t$ is very small, calculating $1-\exp(-\lambda t)$ directly can lose precision because the exponential is close to one. Use $-\operatorname{expm1}(-\lambda t)$, where $\operatorname{expm1}(z)$ computes $\exp(z)-1$ accurately near zero. This is a numerical implementation of the same exact solution, not a different model.

Justification: Let $W(t)=R(t)-R_{\mathrm{eq}}$. Since $k_{\mathrm{in}}q=\lambda R_{\mathrm{eq}}$, the differential equation implies $dW/dt=-\lambda W$. Separating variables gives $W(t)=W(0)\exp(-\lambda t)$, which yields the stated solution. If $R(0)\geq0$, $R_{\mathrm{eq}}\geq0$, and $\lambda>0$, the solution is a weighted average of $R(0)$ and $R_{\mathrm{eq}}$ with weights $\exp(-\lambda t)$ and $1-\exp(-\lambda t)$, both between zero and one. Thus the exact solution remains nonnegative.

Worked Example: Near-Complete Inhibition and Very Large Loss Stimulation

Consider hypothetical parameters $k_{\mathrm{in}}=10$ response units per hour, $k_{\mathrm{out}}=0.1$ per hour, and $R(0)=100$ response units. During a constant-concentration interval, let the production multiplier be $q=0.001$ and the loss multiplier be $g=1001$. Find the response after $0.05$ hour.

The baseline is consistent with the specified rates: $k_{\mathrm{in}}/k_{\mathrm{out}}=10/0.1=100$ response units. Under exposure, the loss constant and equilibrium are $$ \lambda=k_{\mathrm{out}}g=0.1(1001)=100.1\ \text{h}^{-1}, \qquad R_{\mathrm{eq}}=\frac{10(0.001)}{0.1(1001)} =\frac{0.01}{100.1} \approx 0.00009990\ \text{response units}. $$

Substituting into the exact solution gives $$ R(0.05)=0.00009990+(100-0.00009990)\exp[-100.1(0.05)]. $$ The exponent is $-5.005$, and $\exp(-5.005)\approx0.00670434$. Therefore, $$ R(0.05)\approx0.00009990+99.99990010(0.00670434) \approx0.6705\ \text{response units}. $$

The response remains well above its very low exposure-specific equilibrium because the initial response has not yet fully decayed. The calculation also illustrates why the exponential and all dependent arithmetic should be checked: a small error in a reported exponential can produce a materially wrong response.

Large Loss Stimulation, Stiffness, and Integration

For time-varying concentration, a numerical ODE solver may be needed. Strong loss stimulation makes the effective rate $\lambda(t)=k_{\mathrm{out}}g(C(t))$ large. The response then changes rapidly over a time scale of about $1/\lambda(t)$, potentially alongside slower changes in concentration or production. A system is commonly called stiff when widely separated time scales cause fast modes to impose restrictive stability limits on explicit numerical methods, even if those modes have little effect on the solution of interest. Stiffness does not mean the model is invalid; it means that some numerical methods must take very small steps to remain stable.

For the simple equation $dR/dt=-\lambda R$, the explicit Euler method with step size $h$ multiplies the current response by $1-\lambda h$ at each step. Stability requires $|1-\lambda h|\leq1$, or $0\leq\lambda h\leq2$. Even within that stability range, $\lambda h>1$ makes the multiplier negative, which can produce nonphysical negative responses. Accuracy usually requires a smaller step than the stability limit. This condition is a diagnostic for this simple decay equation, not a universal step-size rule for every nonlinear PK/PD system.

Worked Example: Detecting an Unstable Euler Step

Use the parameters from the preceding constant-exposure example, so $\lambda=100.1\ \text{h}^{-1}$, $R_{\mathrm{eq}}\approx0.00009990$, and $R(0)=100$. Suppose explicit Euler uses $h=0.03$ hour for the decay toward that equilibrium.

The Euler update can be written $R_{n+1}-R_{\mathrm{eq}}=(1-\lambda h)(R_n-R_{\mathrm{eq}})$. Here, $$ 1-\lambda h=1-100.1(0.03)=1-3.003=-2.003. $$ Its magnitude exceeds one, so deviations from equilibrium grow in magnitude and change sign at each step. The numerical response can oscillate and become negative even though the exact constant-exposure solution remains positive and decays monotonically toward equilibrium. The step fails the explicit Euler stability condition because $\lambda h=3.003>2$; reducing the step or using a suitable stiff solver is necessary.

Adaptive solvers adjust step size to meet local error controls, while implicit methods such as backward differentiation formulas are often more effective for stiff equations. Neither choice removes the need to check results. Absolute tolerance matters when $R$ approaches zero: a tolerance appropriate for responses near 100 units may be too loose to resolve a response near $10^{-4}$ units. Relative tolerance alone is also insufficient near zero, because relative error is ill-defined at zero and can demand extreme precision as the response becomes tiny. Set tolerances with the response scale and scientific precision in mind.

Diagnostics, Alternatives, and Limitations

A reliable computation should be checked in ways that target the failure mode, rather than judged only by whether the software reports success. Useful diagnostics include:

Refinement: rerun with smaller maximum steps or tighter error tolerances. Important predictions and parameter estimates should not shift materially relative to the required precision.
Physical constraints: inspect whether the response becomes negative, jumps discontinuously without a modeled event, or exceeds plausible bounds. A violation can indicate inadequate resolution, but should not be silently corrected by truncating values to zero.
Rate scales: compare step sizes with the fastest effective rate, including $k_{\mathrm{out}}g(C)$. A step that spans many rapid decay time scales may miss a transient.
Independent calculation: for a constant-concentration interval, compare the ODE solver with the exact solution. Agreement tests the implementation and solver settings under a case with a known answer.
Arithmetic range: check for overflow, underflow, or non-finite values when evaluating high powers, exponentials, and concentration ratios.

There are several complementary approaches, each with limitations. Algebraically stable effect functions reduce cancellation and overflow before integration begins. An exact solution is preferable when its assumptions hold—for example, constant concentration and constant effect multipliers over the interval—but it does not directly solve a model with continuously changing concentration. Adaptive explicit solvers can be efficient for nonstiff intervals, whereas implicit stiff solvers can take larger stable steps in rapidly decaying systems; they may require more computation per step and still need sensible tolerances. Smaller fixed steps can improve an explicit solution, but may be inefficient and do not guarantee accuracy unless convergence is checked.

Numerical stability cannot resolve structural or practical identifiability. A stable calculation may still produce uncertain parameters, and a visually plausible fit does not establish that the mechanism is correct. Conversely, numerical artifacts can distort residuals and make competing mechanisms appear more or less similar than they are. Separate numerical verification from model evaluation and parameter inference.

Common Mistake

Common mistake: “The response became negative, so I will set all negative solver values to zero; that is the model’s prediction.” Correction: For nonnegative production and positive loss multiplier, the exact solution starting from a nonnegative response cannot become negative. A negative value may reveal an overly large step, inadequate error control, or a coding error. Diagnose and correct the numerical calculation rather than altering its output after the fact.

Worked Example: Avoiding Cancellation for a Tiny Time Step

Suppose the response starts at $100$ units and, for a constant exposure, its equilibrium is zero. Let $\lambda=10^{-4}\ \text{h}^{-1}$ and consider $t=10^{-4}$ hour, so $z=\lambda t=10^{-8}$. The exact response is $R(t)=100\exp(-z)$.

Directly forming $1-\exp(-z)$ subtracts two values that are nearly equal; in finite precision, $\exp(-z)$ may round very close to one, erasing digits in the small change. Instead calculate $w=-\operatorname{expm1}(-z)$. Since $-\operatorname{expm1}(-10^{-8})\approx9.99999995\times10^{-9}$, $$ R(t)=100[1-w]\approx100(1-9.99999995\times10^{-9}) \approx99.9999990\ \text{response units}. $$ The small decrease is retained. This matters when evaluating early-time changes or comparing two curves whose difference is small.

Stabilize the rate functions.
Rewrite near-cancelling effects in terms of residual fractions, and evaluate extreme Hill functions without overflowing intermediate powers.

Choose computation to match the time scales.
Use the exact constant-exposure solution when applicable; otherwise consider an adaptive or stiff solver when loss stimulation creates rapid dynamics.

Verify the result.
Check positivity and finiteness, compare with analytic cases, and repeat with tighter tolerances or refined steps to assess convergence.

Key takeaway: Near complete inhibition, preserve the small residual production rate algebraically. Under very large loss stimulation, resolve the resulting fast time scale with an appropriate solver. In both cases, compare computations across tolerances or against an exact solution when available.

Check Your Understanding

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

Why can direct calculation of $1-I_{\max}C^\gamma/(IC_{50}^\gamma+C^\gamma)$ be inaccurate at high concentration, and what rearrangement helps?
For constant exposure, name the equilibrium response and the effective loss constant in the exact solution.
For explicit Euler applied to $dR/dt=-\lambda R$, what stability condition applies, and why might a stable step still produce a negative response?
Give two checks that can reveal an inaccurate numerical solution near an extreme drug effect.

Solutions. (1) The calculation subtracts a term close to one from one, which can lose significant digits. Set $\varepsilon=1-I_{\max}$ and $x=(C/IC_{50})^\gamma$, then use $q=\varepsilon+(1-\varepsilon)/(1+x)$. (2) The effective loss constant is $\lambda=k_{\mathrm{out}}g(C)$, and the equilibrium is $R_{\mathrm{eq}}=k_{\mathrm{in}}q(C)/\lambda$. (3) Stability requires $|1-\lambda h|\leq1$, or $0\leq\lambda h\leq2$. If $1<\lambda h\leq2$, the multiplier is negative, so the Euler response can change sign despite satisfying this stability bound. (4) Examples include tightening tolerances or refining steps and checking whether the result converges, checking nonnegativity and finite values, and comparing the solver with the exact constant-exposure solution.