Metropolis-Hastings for QSP Parameters

Sampling a QSP Posterior

Bayesian calibration of a quantitative systems pharmacology model produces a posterior distribution over parameters, not just one best-fitting parameter vector. For a nonlinear ODE model, evaluating that distribution at a proposed parameter vector may require solving the ODE and computing the observation likelihood. Metropolis–Hastings (MH) uses these evaluations to construct a Markov chain whose stationary distribution is the posterior, even when the posterior’s normalizing constant is unavailable.

Let $\theta$ denote the parameter vector, $y$ the observed data, $p(\theta)$ the prior density and $L(y\mid\theta)$ the likelihood. Write the posterior target density as $\pi(\theta\mid y)$. The ODE model and observation-error assumptions determine the likelihood; the MH algorithm uses that likelihood together with the prior.

Definition: Given a current state $\theta$, a proposal density $q(\theta'\mid\theta)$ generates a candidate state $\theta'$. The Metropolis–Hastings acceptance probability is $$ \alpha(\theta,\theta')= \min\left\{1,\frac{\pi(\theta'\mid y)\,q(\theta\mid\theta')} {\pi(\theta\mid y)\,q(\theta'\mid\theta)}\right\}. $$ The candidate is accepted with probability $\alpha(\theta,\theta')$; otherwise the chain remains at $\theta$. The proposal-density terms account for proposals that are more likely in one direction than the other.

In Bayesian calibration, the posterior is proportional to likelihood times prior:

\pi(\theta\mid y)\propto L(y\mid\theta)p(\theta).

The unknown normalizing constant of the posterior cancels from the ratio of posterior densities. Thus, the algorithm needs only an unnormalized target, such as $L(y\mid\theta)p(\theta)$, evaluated at the current and proposed states.

Why the Acceptance Rule Works

The key requirement is that probability flow between any two states balances in both directions. For two distinct states $\theta$ and $\theta'$, the probability flow from $\theta$ to $\theta'$, weighted by the target density, is proportional to $\pi(\theta\mid y)q(\theta'\mid\theta)\alpha(\theta,\theta')$. Substituting the acceptance rule gives

\pi(\theta\mid y)q(\theta'\mid\theta)\alpha(\theta,\theta') = \min\left\{ \pi(\theta\mid y)q(\theta'\mid\theta), \pi(\theta'\mid y)q(\theta\mid\theta') \right\}.

The right-hand side is unchanged if $\theta$ and $\theta'$ are interchanged. Therefore,

\pi(\theta\mid y)q(\theta'\mid\theta)\alpha(\theta,\theta') = \pi(\theta'\mid y)q(\theta\mid\theta')\alpha(\theta',\theta).

This detailed-balance equality implies that the target is invariant under the transition: if the chain has the target distribution before a step, it has that distribution after the step. Rejections contribute probability of staying at the current state and are part of the transition rule, not missing data. Detailed balance establishes invariance; convergence from other starting states additionally requires suitable conditions, including that the proposal can reach the relevant target support and that the chain is not trapped in a periodic pattern.

For calculation, it is usually safer to work on the log scale. Define the log acceptance ratio

r= \log L(y\mid\theta')-\log L(y\mid\theta) +\log p(\theta')-\log p(\theta) +\log q(\theta\mid\theta')-\log q(\theta'\mid\theta).

Then $\alpha=\min\{1,\exp(r)\}$. If $r\geq 0$, accept with probability 1; if $r<0$, accept with probability $\exp(r)$. This avoids multiplying many small likelihood and prior values. For conditionally independent observations, their likelihood contributions multiply, so their log-likelihood contributions add.

Worked Example: An ODE Likelihood and a Log-Scale Prior

Consider a hypothetical one-compartment decay model, $\frac{dC}{dt}=-kC$, with initial concentration $C(0)=10$ mg/L. Its solution is $C(t;k)=10\exp(-kt)$. At $t=2$ days, the observed concentration is $y=7.6013$ mg/L. Assume normally distributed observation error with known standard deviation $\sigma=0.5$ mg/L, so $Y\mid k\sim N(C(2;k),\sigma^2)$. Use the prior $\log k\sim N(\log(0.10),0.5^2)$, with $k$ measured in inverse days. The current value is $k=0.10$, and the proposed value is $k'=0.08$.

The proposed move is a symmetric random walk in $x=\log k$, so its proposal-density ratio in $x$ is 1. At the current value the model prediction is $10\exp(-0.20)=8.1873$ mg/L; at the proposed value it is $10\exp(-0.16)=8.5214$ mg/L. With a normal likelihood, the log-likelihood ratio is

\log\frac{L(y\mid k')}{L(y\mid k)} = -\frac{(7.6013-8.5214)^2-(7.6013-8.1873)^2}{2(0.5)^2} \approx -1.0065.

The current log parameter equals the prior mean. The prior log-density ratio is therefore

\log\frac{p(\log k')}{p(\log k)} = -\frac{\{\log(0.08/0.10)\}^2}{2(0.5)^2} = -\frac{(-0.22314)^2}{0.5} \approx -0.0996.

Combining the likelihood and prior contributions, with no proposal correction in $x$, gives

r\approx -1.0065-0.0996=-1.1061, \qquad \alpha=\exp(-1.1061)\approx 0.3308.

The proposed parameter is accepted with probability about 0.3308. For example, if the algorithm draws $u=0.20$ from a uniform distribution on $[0,1]$, it accepts because $u<\alpha$. If it draws $u=0.50$, it rejects and retains $k=0.10$. An accepted or retained state is one step in the chain.

The prior ratio above is a density ratio for $x=\log k$, the coordinate in which the proposal is symmetric and the prior is specified. If instead the chain is described using densities in $k$, the transformed prior density includes the Jacobian factor $1/k$; proposal densities must also be expressed in the same coordinate used to calculate the ratio. Mixing densities from different coordinates gives an incorrect acceptance probability.

Proposal Asymmetry and the Hastings Correction

A symmetric proposal satisfies $q(\theta'\mid\theta)=q(\theta\mid\theta')$, so the proposal terms cancel. This is convenient, but it is not a general property of proposals. Whenever the forward and reverse proposal probabilities differ, omitting their ratio changes the target distribution.

Worked Example: Correcting an Asymmetric Proposal

For a simple two-state demonstration, suppose the target probabilities are $\pi(A)=0.4$ and $\pi(B)=0.6$. The proposal moves from $A$ to $B$ with probability $q(B\mid A)=0.8$, and from $B$ to $A$ with probability $q(A\mid B)=0.2$. The acceptance probability for $A$ to $B$ is

\alpha(A,B)= \min\left\{1,\frac{0.6(0.2)}{0.4(0.8)}\right\} =\min\{1,0.375\}=0.375.

For the reverse move,

\alpha(B,A)= \min\left\{1,\frac{0.4(0.8)}{0.6(0.2)}\right\} =1.

The target-weighted flows match: $0.4(0.8)(0.375)=0.12$ from $A$ to $B$, and $0.6(0.2)(1)=0.12$ from $B$ to $A$. If the proposal ratio were incorrectly omitted, the $A$-to-$B$ move would be accepted with probability 1, and the balancing equality would fail.

Putting the ODE Evaluation Into the Chain

At each iteration, the algorithm evaluates a proposed parameter vector through the same model and observation assumptions used to define the posterior. For an ODE-based likelihood, this commonly means solving the ODE at the proposed parameters, extracting model predictions at observation times, and evaluating the observation-error density. The following sequence makes the calculation explicit.

Initialize in the target support.
Choose a starting parameter vector with a defined prior and a valid likelihood. For positive rate constants, for example, use positive values or a suitable unconstrained transformation.

Generate a candidate.
Draw $\theta'$ from the specified proposal $q(\theta'\mid\theta)$, keeping track of the coordinate and proposal density used.

Evaluate the target contributions.
Solve the ODE at $\theta'$, calculate $\log L(y\mid\theta')$ and $\log p(\theta')$, and calculate the reverse and forward log proposal densities.

Accept or retain.
Calculate $r$, set $\alpha=\min\{1,\exp(r)\}$, and accept if a uniform draw $u$ satisfies $u<\alpha$. Otherwise, set the next state equal to $\theta$.

Worked Example: Combining Two ODE Observations

Use the same decay model $C(t;k)=10\exp(-kt)$, now with hypothetical observations at $t=1$ and $t=2$ days. Let the observed concentrations equal the model predictions at the current value $k=0.20$: $y_1=8.1873$ mg/L and $y_2=6.7032$ mg/L. Assume independent normal observation errors, each with known standard deviation $0.2$ mg/L. Propose $k'=0.25$ using a symmetric random walk in $\log k$. Let the prior be $\log k\sim N(\log(0.20),0.4^2)$.

At $k'=0.25$, the predictions are $10\exp(-0.25)=7.7880$ mg/L and $10\exp(-0.50)=6.0653$ mg/L. Independence means the log-likelihood ratio is the sum of the two normal log-density differences:

\log\frac{L(y\mid k')}{L(y\mid k)} = -\frac{(8.1873-7.7880)^2+(6.7032-6.0653)^2}{2(0.2)^2} \approx -7.079.

The log-prior ratio is

\log\frac{p(\log k')}{p(\log k)} = -\frac{\{\log(0.25/0.20)\}^2}{2(0.4)^2} = -\frac{(0.22314)^2}{0.32} \approx -0.1556.

The proposal ratio is 1 in log coordinates. Thus $r\approx-7.079-0.1556=-7.235$, and $\alpha=\exp(-7.235)\approx0.00072$. A uniform draw $u=0.001$ exceeds this acceptance probability, so the proposal is rejected and the chain remains at $k=0.20$. This is a low-probability move under the specified posterior and proposal; it is not evidence that the parameter is impossible.

Assumptions, Diagnostics, and Limitations

The acceptance rule is valid when the target density is defined and the proposal probabilities used in the ratio are correct. A proposed state outside the prior support has zero target density and must be rejected. The chain’s proposal mechanism must also be able to reach the portions of the target one intends to sample; a proposal that permanently excludes a posterior region cannot explore it. For convergence to the target from a chosen initialization, invariance alone is not enough: appropriate irreducibility and aperiodicity conditions are also needed.

For QSP models, a failed or unreliable ODE solve is a practical failure mode. If a solver returns non-finite states, violates model constraints, or has substantial numerical error, the resulting likelihood may be meaningless. Treat genuinely undefined model predictions as having zero likelihood, and check solver tolerances and output validity rather than silently using corrupted values. Numerical error that varies with parameters can distort the target, especially when likelihood differences are small.

Monitor whether proposed states are accepted or rejected and inspect parameter values for obvious pathologies such as repeated solver failures or a chain that remains at one value for long stretches. These are warning signs, not a complete convergence assessment. A high acceptance fraction by itself does not establish that the chain has explored the posterior, and a low fraction can reflect a poor proposal scale or difficult posterior geometry. Detailed proposal tuning and formal convergence diagnostics require additional methods.

MH is useful when the posterior can be evaluated up to a constant but is not available in a form that permits direct sampling. It differs from optimization: optimization seeks a high-density point, whereas MH produces dependent draws intended to represent the posterior distribution. It also differs from a grid approximation, which can become impractical as the number of parameters grows. The basic method is broadly applicable but may explore slowly when parameters are strongly correlated, the posterior has separated regions, or ODE evaluation is expensive. A chain’s draws are not independent, and the MH rule alone does not guarantee that a finite run has adequately explored the posterior.

Common Mistakes

Dropping the proposal correction. This is valid only for a proposal symmetric in the coordinates used for the ratio. Otherwise include $q(\theta\mid\theta')/q(\theta'\mid\theta)$.
Using a prior in one coordinate and a proposal in another. Specify the target and proposal densities consistently; transformed densities require the appropriate Jacobian.
Adding densities instead of log densities. Independent likelihood contributions multiply, while their log contributions add.
Calling a rejection a failed iteration. A rejection produces a repeat of the current state, which is necessary for the chain to have the intended target distribution.
Interpreting an acceptance probability as a posterior probability. It is the conditional probability of accepting this particular proposal from this particular current state, not the probability that a parameter value is true.

Key takeaway: For each proposed QSP parameter vector, combine the likelihood ratio, prior ratio and reverse-to-forward proposal ratio. Accept with probability equal to the smaller of one and this full ratio; otherwise retain the current state.

Check Your Understanding

Answer each question, then compare with the solutions below.

If a proposal is symmetric in the parameter coordinates used, which term in the log acceptance ratio cancels?
A proposed state has a finite likelihood but zero prior density. What is its acceptance probability?
For an asymmetric proposal, suppose the target ratio is $0.5$ and the reverse-to-forward proposal ratio is $0.4$. Find the acceptance probability.
Why does an MH rejection leave the chain at its current state rather than remove that iteration?

Solutions. (1) The proposal contribution $\log q(\theta\mid\theta')-\log q(\theta'\mid\theta)$ is zero. (2) Its posterior target density is zero, so the acceptance probability is zero. (3) The full ratio is $0.5(0.4)=0.20$, giving acceptance probability $0.20$. (4) The retained state is part of the transition rule and helps preserve the target distribution; deleting rejected iterations would change the chain’s distribution.