Bayesian Updating in a Mechanistic Dose-Response Model

From Mechanistic Predictions to Posterior Probabilities

A mechanistic dose-response model links a biological parameter to measurable outcomes. Bayesian updating combines that link with a prior distribution and observed data: parameter values that predict the observations well gain posterior probability, while values that predict them poorly lose probability. The calculation is explicit for a small set of candidate parameter values, making it useful for seeing exactly how prior information and data contribute.

The example below uses a one-compartment concentration model followed by a saturable exposure-response function. Pharmacokinetic parameters are treated as known, while the maximal response is uncertain. This separation keeps the calculation focused: the measured response data update uncertainty about the response parameter, conditional on the specified exposure model.

Definition: For parameter \(\theta\), data \(y\), prior density or probability \(p(\theta)\), and likelihood \(p(y\mid\theta)\), Bayes’ rule gives the posterior \(p(\theta\mid y)=p(y\mid\theta)p(\theta)/p(y)\). The likelihood measures how compatible the observed data are with each parameter value; the denominator \(p(y)\) is the marginal probability of the data and normalizes the posterior so its probabilities sum or integrate to one.

Model and Assumptions

For a single intravenous bolus in a one-compartment model, concentration at time \(t\) is \(C(t)=(D/V)e^{-kt}\), where \(D\) is dose, \(V\) is apparent distribution volume, and \(k\) is the elimination rate constant. The dose and pharmacokinetic parameters are fixed here. The response mean at concentration \(C\) follows a one-site saturable relationship:

\mu(C;E_{\max})=E_{\max}\frac{C}{EC_{50}+C}.

In this tutorial, \(EC_{50}=1\) concentration unit is fixed, so the response is half of \(E_{\max}\) at \(C=1\). The observations \(Y_i\) are modeled as conditionally independent normal variables with this mean and known standard deviation \(\sigma=8\) response units:

Y_i\mid E_{\max}\sim N\left(\mu(C_i;E_{\max}),8^2\right), \qquad p(\mathbf y\mid E_{\max}) =\prod_{i=1}^{n}\frac{1}{8\sqrt{2\pi}} \exp\left[-\frac{\{y_i-\mu(C_i;E_{\max})\}^2}{2(8^2)}\right].

This observation model assumes the mean function is appropriate, residual variability is normal with constant known variance, and errors at the two sampling times are independent conditional on the parameter. It also treats the concentrations as known rather than uncertain. These are substantive modeling assumptions, not consequences of Bayes’ rule. The calculations are appropriate only if the response scale is continuous and the error description is defensible for the hypothetical data.

To make the update hand-calculable, suppose the prior assigns probabilities to three candidate values rather than spreading probability continuously:

P(E_{\max}=85.334)=0.40,\qquad P(E_{\max}=40)=0.35,\qquad P(E_{\max}=22.222)=0.25.

These values form a deliberately small illustrative grid, not a claim that biological knowledge usually supports only three possibilities. For a discrete prior, Bayes’ rule is applied to each candidate and the resulting unnormalized weights are divided by their sum.

Worked Example: Translating Sampling Times Into Exposures

Suppose the hypothetical bolus dose and volume give \(D/V=8\) concentration units, and \(k=0.6931\ \text{h}^{-1}\). Compute concentrations at \(t_1=1\) h and \(t_2=3\) h.

Using \(C(t)=(D/V)e^{-kt}\), the first concentration is \(8e^{-0.6931(1)}\), approximately \(8(0.5000)=4.000\) concentration units. The second is \(8e^{-0.6931(3)}\), approximately \(8(0.1250)=1.000\) concentration unit. Thus the model uses \(C_1=4\) and \(C_2=1\). This calculation assumes the stated bolus model and fixed pharmacokinetics; it does not estimate \(k\) or \(D/V\) from response data.

At \(C=4\), the response fraction is \(4/(1+4)=0.8000\). At \(C=1\), it is \(1/(1+1)=0.5000\). For example, candidate \(E_{\max}=40\) therefore predicts mean responses of \(32\) and \(20\) response units at the two times.

Deriving the Discrete Update

Let the candidate values be \(E_j\), with prior probabilities \(\pi_j\). For independent observations under the specified model, the likelihood at candidate \(E_j\) is the product of its observation densities. The Gaussian normalizing factor is identical across candidates because \(\sigma\) is fixed, so it cancels when posterior probabilities are normalized. The posterior can therefore be calculated from likelihood kernels:

w_j=\pi_j\prod_{i=1}^{n} \exp\left[-\frac{\{y_i-\mu(C_i;E_j)\}^2}{128}\right], \qquad P(E_{\max}=E_j\mid\mathbf y)=\frac{w_j}{\sum_{\ell}w_\ell}.

The denominator is the sum of all unnormalized weights. This is a finite version of the marginal likelihood calculation: the probability of the data is obtained by averaging their likelihood over the prior candidates.

The expression follows directly from Bayes’ rule. For each candidate, the numerator is its prior probability multiplied by the likelihood. Summing these numerators over all mutually exclusive candidates gives the total probability of the data. Dividing each numerator by that total makes the posterior probabilities sum to one:

\sum_j P(E_j\mid\mathbf y) =\sum_j\frac{p(\mathbf y\mid E_j)P(E_j)} {\sum_\ell p(\mathbf y\mid E_\ell)P(E_\ell)} =1.

Worked Example: Updating With Two Response Measurements

Consider hypothetical responses \(y_1=50\) units at \(C_1=4\), and \(y_2=38\) units at \(C_2=1\). Use the prior and observation model above. First calculate the predicted means. For \(E_{\max}=85.334\), they are \(0.8(85.334)=68.2672\) and \(0.5(85.334)=42.667\). For \(E_{\max}=40\), they are \(32\) and \(20\). For \(E_{\max}=22.222\), they are \(17.7776\) and \(11.111\).

For the first observation, the likelihood kernel is \(\exp[-(y_1-\mu_1)^2/128]\). The candidate \(85.334\) gives \(\exp[-(50-68.2672)^2/128]\), approximately \(0.07376\). The candidate \(40\) gives \(\exp[-(50-32)^2/128]=\exp(-2.53125)\), approximately \(0.07956\). The candidate \(22.222\) gives \(\exp[-(50-17.7776)^2/128]\), approximately \(0.00030004\).

For the second observation, the corresponding kernels are \(\exp[-(38-42.667)^2/128]\), approximately \(0.84353\); \(\exp[-(38-20)^2/128]\), approximately \(0.07956\); and \(\exp[-(38-11.111)^2/128]\), approximately \(0.003523\). In particular, the last value uses the squared residual \((26.889)^2\) in the exponent, giving \(\exp(-5.6485)\), approximately \(0.003523\).

Multiply the two observation kernels and the prior probability for each candidate. The unnormalized weights are approximately:

\(0.40(0.07376)(0.84353)=0.02489\) for \(E_{\max}=85.334\).
\(0.35(0.07956)(0.07956)=0.002215\) for \(E_{\max}=40\).
\(0.25(0.00030004)(0.003523)=0.0000002642\) for \(E_{\max}=22.222\).

Their sum is approximately \(0.02710\). Dividing each weight by that sum gives posterior probabilities of approximately \(0.91825\), \(0.08174\), and \(0.00000975\), respectively. The data substantially favor the largest candidate in this particular grid. The smallest candidate has very little posterior probability because both observed responses are poorly predicted by it. The result is conditional on the prior candidates and all the stated model assumptions; it is not a probability that the model itself is true.

Interpreting Posterior Predictions

A posterior probability summarizes uncertainty among the candidates after seeing the data. A posterior prediction averages model predictions over that updated uncertainty. For a future response at a specified concentration, the posterior predictive mean is the posterior-weighted mean response. It is not the same as selecting the candidate with the largest posterior probability and pretending that all other candidates are impossible.

Worked Example: A Posterior-Averaged Response Prediction

Predict the mean response at a future concentration of \(C=2\) concentration units. The response fraction is \(2/(1+2)=2/3\). The candidate means are therefore approximately \(56.89\) units for \(E_{\max}=85.334\), \(26.67\) units for \(E_{\max}=40\), and \(14.81\) units for \(E_{\max}=22.222\).

Weight these means by the posterior probabilities from the preceding calculation:

E[Y_{\mathrm{new}}\mid \mathbf y,C=2] \approx (0.91825)(56.89)+(0.08174)(26.67) +(0.00000975)(14.81) \approx 54.42\ \text{response units}.

This is the posterior predictive mean under the specified model. A future measured response would also include residual variability with standard deviation 8 units, so the predicted mean alone does not describe the full uncertainty in a future observation.

Diagnostics, Alternatives, and Limitations

A numerical update should be accompanied by checks of the model-data relationship. Plot observed responses against concentration and compare them with the candidate model curves. Check whether residuals \(y_i-\mu(C_i;E_{\max})\) show systematic concentration patterns, whether their spread is compatible with the assumed constant standard deviation, and whether measurements plausibly behave independently. With only two observations, these checks have little power: their main value is revealing obvious mismatches, not certifying the assumptions.

A very small likelihood for every candidate is a warning that the candidate set or the observation model may not represent the data adequately. A posterior concentrated on a grid boundary can mean the data favor values beyond the grid; it does not establish that the boundary value is biologically correct. Sensitivity to prior probabilities should also be considered, especially when the data are sparse. Posterior odds between two candidates equal prior odds multiplied by their likelihood ratio, so a strong prior preference can persist if the data provide only modest discrimination.

A common mistake is to say, “The most probable candidate has posterior probability 0.918, so it is a 91.8% chance that this mechanistic model is true.” The calculation supports a narrower statement: conditional on this candidate set, prior probabilities, and observation model, the posterior probability assigned to \(E_{\max}=85.334\) is about 0.918. It does not account for omitted mechanisms, uncertainty in fixed pharmacokinetic parameters, or alternatives outside the grid.

A plug-in approach might choose a single best-fitting parameter and then report predictions as though that parameter were known. Bayesian averaging instead propagates uncertainty across the specified candidates. A continuous-parameter analysis can represent values between grid points, but requires specifying a continuous prior and evaluating or approximating the resulting posterior. Neither approach repairs a misspecified response function or an incorrect error model. The discrete calculation is transparent and useful for learning the update, but its resolution is limited by the candidates supplied.

Key takeaway: Bayesian updating multiplies prior support by the likelihood of the data, then normalizes. The posterior is meaningful only relative to the parameter values considered and the mechanistic and observation assumptions used to calculate that likelihood.

Check Your Understanding

Answer each question, then compare with the concise solutions below.

With \(D/V=8\), \(k=0.6931\ \text{h}^{-1}\), and \(t=2\) h, what concentration does the one-compartment bolus model predict?
At \(C=3\) and \(EC_{50}=1\), what fraction of \(E_{\max}\) is the predicted mean response?
Why does the Gaussian normalizing constant cancel when calculating posterior probabilities over the three candidates?
If two candidates have equal prior probability but one has four times the likelihood of the other, what are their posterior probabilities when these are the only candidates?
What does a posterior probability near one at the largest candidate in a finite grid fail to establish?

Solutions. (1) \(C=8e^{-0.6931(2)}\approx8(0.2500)=2.000\) concentration units. (2) The fraction is \(3/(1+3)=0.7500\), or 75% of \(E_{\max}\). (3) It is the same factor for every candidate, so it appears in both numerator and denominator and cancels. (4) The normalized likelihood weights are \(4/(4+1)=0.8000\) and \(1/(4+1)=0.2000\). (5) It does not establish that the boundary is biologically correct or that omitted parameter values, mechanisms, or observation models are implausible.