Likelihood Estimation for Indirect Response Models

From a Turnover Model to an Observation Likelihood

An indirect response model predicts how a response changes over time because drug action alters its production or loss. Estimation connects that deterministic prediction to noisy measurements. This requires more than choosing an optimization algorithm: the analyst must specify what is random, how measurement errors behave, and which parameter values are admissible.

Consider the canonical model in which drug concentration inhibits response production while loss remains first order. Let $R(t)$ be the response, $C(t)$ the drug concentration, $k_{\mathrm{in}}$ the baseline production rate, $k_{\mathrm{out}}$ the first-order loss-rate constant, $I_{\max}$ the maximum fractional inhibition, and $IC_{50}$ the concentration giving half of that maximum inhibition. The model is

\frac{dR(t)}{dt} = k_{\mathrm{in}}\left(1-\frac{I_{\max}C(t)}{IC_{50}+C(t)}\right) -k_{\mathrm{out}}R(t).

When $C(t)=0$ and the system is at baseline steady state, production equals loss: $k_{\mathrm{in}}=k_{\mathrm{out}}R_0$, where $R_0$ is the baseline response. Thus one may parameterize the model using $R_0$ and $k_{\mathrm{out}}$, setting $k_{\mathrm{in}}=k_{\mathrm{out}}R_0$. This parameterization makes the baseline condition explicit; it does not assert that a measured predose value is exactly $R_0$.

Definition: An observation likelihood is the joint probability density of the measured responses, viewed as a function of the model parameters with the observed data held fixed. For independent continuous observations $Y_i$ with model-predicted means $R(t_i;\theta)$ and Gaussian measurement error with standard deviation $\sigma$, the likelihood is $$ L(\theta,\sigma)=\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma} \exp\left[-\frac{\{y_i-R(t_i;\theta)\}^2}{2\sigma^2}\right], $$ where $y_i$ is the observed response, $n$ is the number of observations, and $\theta$ denotes the model parameters.

This likelihood assumes that, conditional on the model and its parameters, observations are independent and normally distributed with a common variance $\sigma^2$, and that the model prediction is their conditional mean. The standard deviation $\sigma$ has the same units as the response. If measurement variability changes with response magnitude, or measurements are correlated, this likelihood is not appropriate without modification.

Constraints and Initial Conditions

For the inhibitory model, a common biologically constrained parameter space is $R_0>0$, $k_{\mathrm{out}}>0$, $IC_{50}>0$, $0\le I_{\max}\le1$, and $\sigma>0$. The upper bound on $I_{\max}$ represents inhibition no greater than complete suppression of production. If the scientific mechanism permits a different range, that range must be justified rather than imposed automatically. Concentration and $IC_{50}$ must use the same units; $k_{\mathrm{out}}$ has units of inverse time.

The initial condition is part of the prediction, not a harmless software default. If observations begin after a drug-free period long enough to establish baseline, set $R(0)=R_0$. If the response at the start is not known to be at steady state, use an initial response $R_{\mathrm{init}}$ and set $R(0)=R_{\mathrm{init}}$. Treating $R_{\mathrm{init}}$ as known when it is not can force the model to explain an initial offset by changing turnover or drug-effect parameters. Conversely, estimating both $R_{\mathrm{init}}$ and $R_0$ can weaken information about baseline if the data contain little predose or recovery information.

For fixed parameters and piecewise constant concentration over an interval of length $\Delta$, the equation has a convenient exact propagation. Define the interval-specific production-adjusted equilibrium as $R_{\mathrm{ss}}=R_0[1-I_{\max}C/(IC_{50}+C)]$. Then

R(t+\Delta)=R_{\mathrm{ss}}+ \{R(t)-R_{\mathrm{ss}}\}\exp(-k_{\mathrm{out}}\Delta).

For continuously varying concentration, the same likelihood is used, but the predicted mean is obtained by solving the differential equation numerically or by a valid analytic solution. In either case, the prediction at each observation time must use the dosing and concentration history and the specified initial condition.

Why Maximizing the Likelihood Gives Least Squares Here

Taking logarithms turns the product of densities into a sum. With $\mu_i=R(t_i;\theta)$, the Gaussian log-likelihood is

\ell(\theta,\sigma) =-n\log(\sqrt{2\pi}\sigma) -\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i-\mu_i)^2.

For fixed $\sigma$, maximizing this expression over $\theta$ is equivalent to minimizing the sum of squared residuals, provided the same independent, constant-variance Gaussian assumptions hold. This equivalence explains why nonlinear least squares can estimate the mean-model parameters in this setting. It does not make least squares universally equivalent to likelihood estimation: unequal variances, non-Gaussian errors, or correlation change the likelihood and generally change the objective function.

For a fixed mean model, maximizing over $\sigma$ gives $\hat{\sigma}^2=\sum_i(y_i-\mu_i)^2/n$. This is the maximum-likelihood estimate, not the degrees-of-freedom-adjusted variance estimator sometimes used for unbiased estimation in a linear model.

Worked Example: A Likelihood Contribution at One Time

Suppose a hypothetical response measurement is $y=84$ response units and the model predicts $R(t;\theta)=85$ response units. Under the stated Gaussian observation model, suppose $\sigma=4$ response units. The single-observation log-likelihood contribution is

\ell_i=-\log(\sqrt{2\pi}\sigma) -\frac{\{y-R(t;\theta)\}^2}{2\sigma^2} =-\log(4\sqrt{2\pi})-\frac{(84-85)^2}{2(4^2)} \approx -2.3365.

The squared standardized discrepancy is $1^2/16=0.0625$, so its contribution after multiplication by one half is $0.03125$. The log density is therefore approximately $-2.3052-0.03125=-2.3365$. A density, and hence its logarithm, is defined relative to the response measurement scale; the likelihood is not itself a probability that the parameter is true.

Worked Estimation for a Canonical Response Shape

A useful conditional calculation arises when the drug-effect and turnover parameters are fixed, so the model-predicted response at each time is proportional to the baseline $R_0$. Write the predictions as $\mu_i=R_0 f_i$, where $f_i$ are known response factors from the model solution and $f_1=1$ for a baseline observation. With independent Gaussian errors of common variance, maximizing likelihood over $R_0$ is equivalent to minimizing $\sum_i(y_i-R_0f_i)^2$. Differentiating that sum with respect to $R_0$ and setting the derivative to zero yields the unconstrained minimizer. Under the strict constraint $R_0>0$, this formula is the constrained estimate only when $\sum_i f_i y_i>0$.

\frac{d}{dR_0}\sum_i(y_i-R_0f_i)^2 =-2\sum_i f_i(y_i-R_0f_i)=0, \qquad \hat R_0=\frac{\sum_i f_i y_i}{\sum_i f_i^2}.

When the numerator is nonpositive, the strict constraint $R_0>0$ means the likelihood has no attained maximum; its supremum is approached as $R_0\to0^+$. This conditional estimate is not a replacement for fitting all parameters jointly. It isolates the baseline-scale calculation and shows exactly how the likelihood uses measurements across time.

Worked Example: Conditional Maximum-Likelihood Estimate of Baseline

Suppose a hypothetical data set has responses $100,\ 84,\ 72,\ 67$ response units at four times. With other model parameters fixed, the corresponding predicted response factors are $1,\ 0.835160,\ 0.724664,\ 0.650597$. Assume independent Gaussian errors with common variance and use $R_0$ as the only unknown mean-model parameter.

The least-squares expression, and hence the conditional Gaussian maximum-likelihood estimate, is

\hat R_0= \frac{(1)(100)+(0.835160)(84)+(0.724664)(72)+(0.650597)(67)} {1^2+0.835160^2+0.724664^2+0.650597^2} = \frac{265.919247}{2.6459066} \approx 100.5021\ \text{response units}.

The fitted means are approximately $100.5021,\ 83.9353,\ 72.8303,\ 65.3864$ response units. Subtracting these from the observations gives residuals of approximately $-0.5021,\ 0.0647,\ -0.8303,\ 1.6136$ response units. Their squared sum is approximately $3.5494$ response-units squared. If $\sigma$ is also estimated by maximum likelihood, then

\hat{\sigma}=\sqrt{\frac{3.5494}{4}} \approx 0.9420\ \text{response units}.

The divisor is $n=4$ because this is the Gaussian maximum-likelihood estimate. The calculation estimates $R_0$ conditional on the specified factors; uncertainty or misspecification in the fixed factors is not represented by this conditional calculation.

Worked Example: Prediction and Likelihood Under Constant Concentration

Consider a hypothetical interval beginning at $R(0)=100$ response units. Let $k_{\mathrm{out}}=0.20\ \text{h}^{-1}$, $I_{\max}=0.60$, $IC_{50}=2\ \text{mg/L}$, and let the concentration remain at $3\ \text{mg/L}$ for two hours. Assume $R_0=100$ response units, so $k_{\mathrm{in}}=k_{\mathrm{out}}R_0=20$ response units per hour.

The fractional inhibition of production is $0.60(3)/(2+3)=0.36$, leaving $1-0.36=0.64$ of baseline production. Thus $R_{\mathrm{ss}}=100(0.64)=64$ response units. The predicted response after two hours is

R(2)=64+(100-64)\exp[-(0.20)(2)] =64+36\exp(-0.4) \approx 88.13\ \text{response units}.

If the hypothetical observed response is $90$ response units and $\sigma=3$ response units, its residual is $90-88.13=1.87$ response units. The log-likelihood contribution is approximately

-\log(3\sqrt{2\pi})-\frac{(1.87)^2}{2(3^2)} \approx -2.212.

This contribution is calculated after checking the model conditions: concentration is constant on the interval, the starting response is specified, and the chosen parameters satisfy the stated positivity and inhibition constraints.

Practical Estimation, Diagnostics, and Alternatives

In a full fit, the optimizer evaluates the model at the observation times for a candidate parameter vector, calculates the log-likelihood, and searches the constrained parameter space for a maximum. Positive parameters are often represented on a log scale during optimization; a fractional parameter such as $I_{\max}$ can be represented through a logistic transformation when strict interior values are desired. If a parameter can genuinely equal a boundary value, the fitting procedure must allow that boundary rather than exclude it inadvertently. Starting values should be scientifically plausible, and multiple starting points help reveal sensitivity to initialization or local optima.

Specify the mean model.
State the turnover equation, concentration input, observation times, and whether the initial response is fixed at baseline or estimated.

Specify the observation model.
State the distribution, independence or covariance assumptions, and the units and role of the residual scale.

Constrain and optimize.
Use scientifically justified parameter bounds, evaluate predictions at every observation time, and maximize the log-likelihood.

Check the fit.
Confirm numerical convergence, inspect whether estimates sit at bounds, and assess whether the assumed error scale and mean trajectory are compatible with the observed data.

Warning signs include a parameter repeatedly reaching a constraint, markedly different estimates from different starting values, a nearly flat likelihood over a wide parameter range, or strong dependence of estimates on the assumed initial condition. These can indicate weak information, practical non-identifiability, an unsuitable model, or numerical problems; a successful optimizer termination alone does not establish that the parameters are well determined. The earlier material on profile likelihood provides a useful way to examine uncertainty when a simple local approximation is doubtful.

Under independent, equal-variance Gaussian errors, nonlinear least squares and Gaussian maximum likelihood give the same estimates for the mean-model parameters. Likelihood estimation is more general because its form can be adapted to heteroscedastic Gaussian errors, correlated errors, or another scientifically justified observation distribution. Such changes require a new likelihood and new assumptions; they are not merely different optimizer settings. This tutorial’s Gaussian likelihood also treats the supplied concentration history as known. Uncertainty in that history or additional sources of variability require an expanded model.

Common Pitfall and Limitations

Common mistake: “The baseline response is measured before dosing, so I can fix $R_0$ to that observed value without affecting estimation.” A baseline observation is noisy under the stated likelihood. Fixing $R_0$ to it treats that one measurement as error-free and can distort the estimated turnover and drug-effect parameters. Instead, include the baseline observation in the likelihood and estimate $R_0$, unless an independent design or a justified model constraint makes it known.

A related error is to fit an unconstrained optimizer and interpret a negative $k_{\mathrm{out}}$, negative $IC_{50}$, or biologically impossible inhibition as a meaningful estimate. Such values violate the model’s parameter definitions. Constraints make the model coherent, but they do not create information: if observations do not distinguish parameter combinations, estimates may remain highly uncertain even within the admissible region. Finally, a Gaussian likelihood can be sensitive to extreme observations, and a single common $\sigma$ may be inadequate when variability changes with response. The observation model should be assessed rather than accepted solely because it is convenient.

The likelihood provides a disciplined link between an indirect-response trajectory and measured data. Its value depends on three choices being explicit and defensible: the response mechanism used to calculate the mean, the initial condition that starts that trajectory, and the probability model for deviations of observations from it.

Check Your Understanding

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

Under independent Gaussian errors with common $\sigma$, what quantity is minimized when maximizing the likelihood over the mean-model parameters with $\sigma$ fixed?
For the inhibitory model, what baseline relationship connects $k_{\mathrm{in}}$, $k_{\mathrm{out}}$, and $R_0$ at steady state?
Why does the maximum-likelihood estimate of $\sigma^2$ divide the residual sum of squares by $n$, not by $n-p$?
If the response is not at baseline steady state when observations begin, what should be reconsidered in the model specification?
Does a converged optimizer establish that every parameter is precisely estimated? Explain.

Solutions. (1) It minimizes $\sum_i\{y_i-R(t_i;\theta)\}^2$, because the other terms in the Gaussian log-likelihood are constant with respect to $\theta$. (2) At drug-free steady state, $k_{\mathrm{in}}=k_{\mathrm{out}}R_0$. (3) Direct maximization of the Gaussian likelihood over $\sigma^2$ gives the residual sum of squares divided by the number of observations $n$; $n-p$ is an adjustment used in a different variance-estimation context. (4) Use or estimate an initial response $R_{\mathrm{init}}$, rather than automatically imposing $R(0)=R_0$. (5) No. Convergence only indicates that the numerical procedure met its stopping criterion; weak information, flat likelihoods, boundary estimates, and sensitivity to starting values still require investigation.