Likelihoods for ODE-Based QSP Observations

From a Mechanistic Trajectory to a Likelihood

A QSP model may describe how latent concentrations, cell populations or pathway activities change continuously over time. Experimental data, by contrast, are measured at particular times, often with error and sometimes on a different scale from the model state. Calibration requires a link between these two objects: an observation model that assigns a probability density to the measured values, conditional on the ODE predictions and the model parameters.

The distinction matters. An ODE state is not automatically a measurement. A model may predict an unobserved tissue concentration while the assay measures a related plasma biomarker; an observation function must connect them. Likewise, the model state exists along a continuous trajectory, but the likelihood is constructed for the values actually recorded at specified observation times.

Definition: Let \(x(t;\theta)\) denote the latent state vector obtained by solving an ODE system with parameters \(\theta\), specified initial conditions and known inputs. At observation time \(t_i\), let \(h_i(x(t_i;\theta),\theta)\) be the model-predicted mean on the measurement scale. The observed value \(y_i\) is linked to that prediction by an observation model. The likelihood \(L(\theta;y)\) is the joint probability density of the observed data \(y\), viewed as a function of \(\theta\) with the data held fixed.

For example, an ODE may describe a latent biomarker state \(B(t)\), while an assay measures \(B(t)\) plus instrument error. In another model, an observation function may map several latent states to one measured biomarker. The time \(t_i\) tells us when to evaluate the trajectory; it is not itself the measured biomarker value. If the assay occurs at \(t_i\), the prediction is formed from the ODE solution at that time, not from a nearby time chosen for convenience.

Independent Gaussian Observation Likelihood

A useful starting model assumes that each measured value differs from its predicted mean by an independent, normally distributed error with known standard deviation \(\sigma_i\). The standard deviation is expressed in the same units as the measurement. Under this model, the density for one observation is a normal density centred on the ODE-based prediction.

y_i \mid \theta \sim N\!\left(h_i(x(t_i;\theta),\theta),\,\sigma_i^2\right)

Conditional independence means that, given the model predictions and parameters, the measurement errors for different observations are independent. With this assumption, the joint density is the product of the individual densities. For a single biomarker with \(m\) observations, this gives:

L(\theta;y) = \prod_{i=1}^{m} \frac{1}{\sqrt{2\pi}\sigma_i} \exp\!\left[ -\frac{\left(y_i-h_i(x(t_i;\theta),\theta)\right)^2}{2\sigma_i^2} \right]

The likelihood is a density, not generally a probability of the exact observed data. Its units are the inverse of the measurement units multiplied across observations. For numerical work, it is often more convenient to calculate the log likelihood, because products become sums:

\log L(\theta;y) = -\sum_{i=1}^{m} \left[ \log(\sqrt{2\pi}\sigma_i) + \frac{\left(y_i-h_i(x(t_i;\theta),\theta)\right)^2}{2\sigma_i^2} \right]

This expression makes the assumptions visible: the ODE and observation function determine the means, and the specified \(\sigma_i\) values determine the error scale. If a common standard deviation \(\sigma\) is used, the likelihood depends on the sum of squared residuals. That relationship explains why, with independent Gaussian errors and fixed common \(\sigma\), parameter values that produce smaller squared residuals have larger likelihoods.

Worked Example: Likelihood for Three Measurements Along a Decay Trajectory

Suppose a hypothetical biomarker follows \(dB/dt=-kB\), with initial state \(B(0)=10\) pg/mL and \(k=\log(2)/2\) per hour. The model solution is \(B(t)=10e^{-kt}\). Measurements are taken at 1, 2 and 3 hours. The corresponding model predictions are approximately 7.0711, 5.0000 and 3.5355 pg/mL; the measured values are 7.2711, 4.7000 and 3.6355 pg/mL. Assume independent normal errors with known \(\sigma=0.5\) pg/mL for each measurement.

The residuals, defined as observed minus predicted values, are \(0.2000\), \(-0.3000\) and \(0.1000\) pg/mL. Their squared sum is \(0.2000^2+(-0.3000)^2+0.1000^2=0.1400\;(\text{pg/mL})^2\). Substituting into the joint likelihood gives:

L = \left(\frac{1}{\sqrt{2\pi}(0.5)}\right)^3 \exp\!\left(-\frac{0.1400}{2(0.5)^2}\right) \approx 0.3839\;(\text{pg/mL})^{-3}

The exponent is \(-0.1400/0.5000=-0.2800\); the density prefactor is approximately \(0.5079\;(\text{pg/mL})^{-3}\), and \(0.5079e^{-0.2800}\approx0.3839\;(\text{pg/mL})^{-3}\). This value is the joint density at these three measurements under the stated parameter and error assumptions. It is not the probability that the model is correct.

Multiple Biomarkers and Observation Schedules

For multiple biomarkers, each observation needs its own identity, time, measured value, model prediction and error specification. Write \(y_{ij}\) for the measurement of biomarker \(j\) at its observation time \(t_{ij}\), and let \(h_j(x(t_{ij};\theta),\theta)\) be its predicted mean. Under independent normal errors with standard deviation \(\sigma_j\), the likelihood is the product over all recorded biomarker-time pairs. Different biomarkers may have different measurement units and error scales.

L(\theta;y) = \prod_{j=1}^{J}\prod_{i=1}^{m_j} \frac{1}{\sqrt{2\pi}\sigma_j} \exp\!\left[ -\frac{\left(y_{ij}-h_j(x(t_{ij};\theta),\theta)\right)^2}{2\sigma_j^2} \right]

Worked Example: Combining Two Biomarkers With Different Units

A hypothetical ODE model predicts biomarker A at two times as 4.0 and 3.0 ng/mL, and biomarker B as 20.0 and 18.0 pmol/L. The observations are 4.1 and 2.8 ng/mL for A, and 20.5 and 17.5 pmol/L for B. Assume independent normal errors, with \(\sigma_A=0.2\) ng/mL and \(\sigma_B=1.0\) pmol/L.

For A, the residuals are \(0.1\) and \(-0.2\) ng/mL. Its likelihood contribution is:

L_A = \left(\frac{1}{\sqrt{2\pi}(0.2)}\right)^2 \exp\!\left[-\frac{0.1^2+(-0.2)^2}{2(0.2)^2}\right] \approx 2.130\;(\text{ng/mL})^{-2}

Here the exponent is \(-0.0500/0.0800=-0.6250\). For B, the residuals are \(0.5\) and \(-0.5\) pmol/L, giving:

L_B = \left(\frac{1}{\sqrt{2\pi}(1.0)}\right)^2 \exp\!\left[-\frac{0.5^2+(-0.5)^2}{2(1.0)^2}\right] \approx 0.1240\;(\text{pmol/L})^{-2}

The joint likelihood is \(L_A L_B\approx0.2640\;(\text{ng/mL})^{-2}(\text{pmol/L})^{-2}\). Multiplication is justified by the assumed conditional independence. If observations were missing at some scheduled times, those absent measurements would not contribute factors: the product is over recorded observations, not every time point at which the ODE could have been evaluated.

Why the Likelihood Has This Form

The likelihood construction follows directly from the observation model. For each observation, subtract its predicted mean from the measured value to obtain a residual. Under the normal model, the density of that value is the normal density evaluated at its residual from the mean. Conditional independence then permits multiplication of those densities. Taking logarithms converts the product to a sum without changing which parameter value maximises it, because the logarithm is strictly increasing.

Solve the state equations.
For a proposed \(\theta\), solve the ODE using the specified initial conditions and inputs to obtain \(x(t;\theta)\).

Evaluate the observation function.
At each actual observation time, compute the predicted measurement \(h_i(x(t_i;\theta),\theta)\).

Form residual densities.
For every measured value, evaluate its density using the stated error model and scale.

Combine the contributions.
Multiply the densities under conditional independence, or add their logarithms to obtain the log likelihood.

For example, suppose the residuals are \(r_i=y_i-h_i(x(t_i;\theta),\theta)\) and the common error standard deviation \(\sigma\) is fixed. The parameter-dependent part of the log likelihood is \(-\sum_i r_i^2/(2\sigma^2)\); the remaining normalising term is constant in \(\theta\). Thus, for two parameter vectors \(\theta_A\) and \(\theta_B\), their log-likelihood difference is:

\log L(\theta_B;y)-\log L(\theta_A;y) = -\frac{\sum_i r_i(\theta_B)^2-\sum_i r_i(\theta_A)^2}{2\sigma^2}

This identity is a direct consequence of subtracting the two log-likelihood expressions: their identical normalising terms cancel. It applies only when the same fixed \(\sigma\) and the same observations are used for both parameter vectors.

Worked Example: Comparing Two Parameter-Dependent Trajectories

Return to the three decay measurements in the first example, with \(\sigma=0.5\) pg/mL. Suppose two parameter vectors produce residuals \((0.2,-0.3,0.1)\) and \((0.4,-0.1,0.2)\) pg/mL, respectively. The first squared-residual sum is \(0.1400\;(\text{pg/mL})^2\); the second is \(0.1600+0.0100+0.0400=0.2100\;(\text{pg/mL})^2\).

The log-likelihood difference, second parameter vector minus first, is:

\log L(\theta_B;y)-\log L(\theta_A;y) = -\frac{0.2100-0.1400}{2(0.5)^2} = -0.1400

Therefore \(L(\theta_B;y)/L(\theta_A;y)=e^{-0.1400}\approx0.8694\). For these data and assumptions, the second vector has about 0.8694 times the likelihood of the first. This comparison is not a posterior probability or a probability that either vector is true; a prior distribution is needed to form a posterior distribution.

Assumptions, Alternatives and Diagnostics

A likelihood is meaningful only relative to its assumptions. For the independent normal construction above, the conditional mean must be adequately represented by the ODE and observation function, the error distribution and scales must be appropriate, and the errors must be independent if the product factorisation is used. The model also requires specified initial conditions and inputs, and a numerical ODE solution accurate enough at the observation times. These assumptions concern different parts of the construction: a poor trajectory prediction is not repaired merely by changing the error distribution, and an incorrect independence assumption is not resolved by fitting the mean more closely.

The independent normal model is not the only possible observation model. For a positive biomarker whose variability is more plausibly additive on the log scale, one may model \(\log(y_i)\) as normal around a predicted log mean. If the likelihood is then expressed as a density for the original \(y_i\), the change of variables contributes a factor \(1/y_i\) for each observation. Alternatively, a covariance model can represent correlated measurement errors; its joint density is not the product of independent marginal densities. These choices alter the likelihood and must be stated rather than silently substituted.

Useful checks during construction include:

Confirm that each measured value is paired with the correct biomarker, units and actual observation time.
Plot residuals against time and predicted values. Persistent trends, changing spread or clusters can indicate that the mean or error assumptions are inadequate.
Check that predicted means and measurements are on the same scale. A latent state and an assay readout may require an observation function between them.
Check solver tolerances and evaluate predictions at the exact observation times. Solver inaccuracies can distort residuals and hence the likelihood, especially when measurement error is small.
Verify that standard deviations use the measurement units and that the likelihood includes all and only the recorded observations.

A common mistake is to treat the ODE state values at every integration step as independent data points. The solver’s internal time steps are numerical calculations, not additional observations. The correction is to evaluate the prediction at each recorded time and include one likelihood contribution for each recorded measurement. Another mistake is to multiply marginal densities when errors are known to be correlated; in that case, the joint error distribution must represent that dependence.

Key takeaway: An ODE-based likelihood combines model-predicted measurement means at the recorded times with an explicit probability model for measurement error. State variables, observation functions, measured values and error assumptions are distinct ingredients; each must be identified before the likelihood is evaluated.

Check Your Understanding

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

An ODE predicts a latent tissue concentration continuously, but an assay records plasma biomarker values at three times. What quantities enter the likelihood, and what role does the observation function play?
Three independent measurements have residuals \(0.1\), \(-0.1\) and \(0.2\) units, with common known \(\sigma=0.2\) units. What is their sum of squared residuals?
If an observation time is missing and no value was measured, should its normal density factor be included in the likelihood?
Under a common fixed \(\sigma\), two parameter vectors have squared-residual sums 0.30 and 0.50. Which has the larger likelihood, and why?
Why is the likelihood value not the probability that a parameter vector is true?

Solutions. (1) Include the measured plasma values, their recorded times, and model predictions at those times; the observation function maps the latent tissue or other model states to the predicted plasma measurement scale. (2) The sum is \(0.1^2+(-0.1)^2+0.2^2=0.0600\) squared units. (3) No. The likelihood product contains contributions for recorded measurements, not unobserved scheduled values. (4) The vector with squared-residual sum 0.30 has the larger likelihood, because with fixed common \(\sigma\), the log likelihood decreases as the sum of squared residuals increases. (5) The likelihood is the data density as a function of the parameter, with the observed data fixed. A probability distribution for parameters requires additional ingredients, such as a prior distribution combined with the likelihood.