Identifiability of Binding, Internalization, and Target Turnover

Why Identifiability Matters in TMDD

A fitted TMDD model can reproduce observed drug concentrations and target measurements while leaving some underlying rates poorly determined. This is more than a numerical inconvenience. If binding, complex internalization, or target turnover cannot be distinguished by the available observations, a precise-looking estimate may depend mainly on model assumptions, parameter constraints, or optimization settings rather than information in the data.

The question here is not simply whether an optimizer converges. It is whether the observed concentration and target data contain enough information to distinguish the parameters of interest. The full TMDD framework and the roles of target synthesis, degradation, complex internalization, and binding are assumed known from the earlier material on target turnover and complex internalization.

Definition: A model is structurally identifiable for a parameter if ideal, noiseless observations over the specified observation conditions uniquely determine that parameter. A parameter is practically identifiable when the actual finite, noisy data constrain its value sufficiently for useful estimation. A structurally identifiable parameter can still be practically unidentifiable when the data are sparse, noisy, or weakly informative.

Identifiability is conditional on the model, the observation variables, and the observation conditions. For example, observing free drug and total target is not the same information as observing free drug, free target, and complex separately. A conclusion about identifiability therefore needs to say which outputs and observation model were used.

Where the Parameters Enter the Model

For illustration, let \(D\), \(R\), and \(C\) denote free drug, free target, and drug-target complex concentrations. A simple full TMDD system with drug input \(u(t)\), volume \(V\), nonspecific free-drug clearance \(CL\), target synthesis \(k_{\mathrm{syn}}\), free-target degradation \(k_{\mathrm{deg}}\), association rate \(k_{\mathrm{on}}\), dissociation rate \(k_{\mathrm{off}}\), and complex internalization rate \(k_{\mathrm{int}}\) is:

\begin{aligned} \frac{dD}{dt} &= \frac{u(t)}{V}-\frac{CL}{V}D-k_{\mathrm{on}}DR+k_{\mathrm{off}}C,\\ \frac{dR}{dt} &= k_{\mathrm{syn}}-k_{\mathrm{deg}}R-k_{\mathrm{on}}DR+k_{\mathrm{off}}C,\\ \frac{dC}{dt} &= k_{\mathrm{on}}DR-(k_{\mathrm{off}}+k_{\mathrm{int}})C. \end{aligned}

This formulation treats internalization as removal of complex, and therefore removes both bound drug and target through that complex pathway. In the absence of drug, the target steady state is \(R_0=k_{\mathrm{syn}}/k_{\mathrm{deg}}\). That ratio is immediately important: a baseline target measurement alone does not separately reveal its numerator and denominator.

A practical assessment should also specify how model states become observations. For example, measured drug concentrations may correspond to free or total drug, while a target assay may report free target or total target. In a likelihood-based analysis, the observation model specifies the distribution of the measurements around the model predictions, including the relevant measurement error. Identifiability diagnostics are only as meaningful as this connection between states and data.

Two Common Sources of Confounding

The first confounding pattern concerns target turnover. At baseline without drug, the equation for target is \(dR/dt=k_{\mathrm{syn}}-k_{\mathrm{deg}}R\). At steady state, setting the derivative to zero yields \(k_{\mathrm{syn}}=k_{\mathrm{deg}}R_0\), or \(R_0=k_{\mathrm{syn}}/k_{\mathrm{deg}}\). Thus, any pair of positive rates with the same ratio gives the same baseline prediction. A baseline measurement cannot resolve the two rates by itself.

The second pattern involves binding and internalization. If the complex changes quickly relative to the other states, a quasi-steady approximation sets \(dC/dt\) near zero. It follows that:

C \approx \frac{k_{\mathrm{on}}DR}{k_{\mathrm{off}}+k_{\mathrm{int}}}.

Under this approximation, complex concentration at a given \(D\) and \(R\) depends on the combination \((k_{\mathrm{off}}+k_{\mathrm{int}})/k_{\mathrm{on}}\). A measurement of complex at one condition may constrain that combination without separately determining all three rates. This does not mean the full dynamic model always has the same structural ambiguity: time-varying outputs can contain additional information. It does mean that data dominated by the quasi-steady relationship can support a long, shallow trade-off among parameters.

Worked Example: Baseline Target Does Not Separate Turnover Rates

Suppose a hypothetical study measures a baseline target concentration of \(100\ \mathrm{pM}\), with no drug present. Consider two possible parameter pairs: \(k_{\mathrm{syn}}=10\ \mathrm{pM/day}\) and \(k_{\mathrm{deg}}=0.10\ \mathrm{day}^{-1}\); or \(k_{\mathrm{syn}}=20\ \mathrm{pM/day}\) and \(k_{\mathrm{deg}}=0.20\ \mathrm{day}^{-1}\).

For the first pair, \(R_0=k_{\mathrm{syn}}/k_{\mathrm{deg}}=10/0.10=100\ \mathrm{pM}\). For the second, \(R_0=20/0.20=100\ \mathrm{pM}\). Both pairs predict the observed baseline exactly. With only this baseline observation, the data determine the ratio, not either rate separately.

Now suppose a hypothetical target-depletion observation starts at \(R(0)=0\), and the measured concentration is \(63.2\ \mathrm{pM}\) after 5 days. Under drug-free turnover, \(R(t)=R_0+(R(0)-R_0)e^{-k_{\mathrm{deg}}t}\). Substitution gives \(63.2=100(1-e^{-5k_{\mathrm{deg}}})\), so \(e^{-5k_{\mathrm{deg}}}=0.368\) and \(k_{\mathrm{deg}}=-\ln(0.368)/5\approx0.200\ \mathrm{day}^{-1}\). Combining this rate with the baseline gives \(k_{\mathrm{syn}}=k_{\mathrm{deg}}R_0=0.200(100)=20.0\ \mathrm{pM/day}\). The recovery observation adds information about the rate of approach to baseline; the baseline alone did not.

Worked Example: A Quasi-Steady Complex Measurement Constrains a Combination

Suppose hypothetical measurements at one condition give \(D=10\ \mathrm{pM}\), \(R=20\ \mathrm{pM}\), and \(C=40\ \mathrm{pM}\). Under the quasi-steady approximation, \(C=k_{\mathrm{on}}DR/(k_{\mathrm{off}}+k_{\mathrm{int}})\). One possible parameter set is \(k_{\mathrm{on}}=0.20\ \mathrm{pM}^{-1}\mathrm{day}^{-1}\), \(k_{\mathrm{off}}=0.40\ \mathrm{day}^{-1}\), and \(k_{\mathrm{int}}=0.60\ \mathrm{day}^{-1}\).

Substituting gives \(C=0.20(10)(20)/(0.40+0.60)=40/1.00=40\ \mathrm{pM}\). A different set, \(k_{\mathrm{on}}=0.10\ \mathrm{pM}^{-1}\mathrm{day}^{-1}\), \(k_{\mathrm{off}}=0.20\ \mathrm{day}^{-1}\), and \(k_{\mathrm{int}}=0.30\ \mathrm{day}^{-1}\), gives \(C=0.10(10)(20)/(0.20+0.30)=20/0.50=40\ \mathrm{pM}\). Both sets give the same predicted complex concentration under this relation.

The example demonstrates a limitation of this single quasi-steady relationship, not proof that these parameter sets predict every output of the full dynamic model identically. A practical identifiability assessment must use the actual observed outputs and their time courses.

Profile Likelihood: A Direct Practical Diagnostic

For a parameter of interest, the profile likelihood fixes that parameter at a sequence of values and re-estimates the other parameters at each value. Let \(\theta\) denote the full parameter vector, and let \(\ell(\theta)\) be the log likelihood under the specified observation model. For parameter \(\theta_j\), define:

\ell_{\mathrm{prof}}(\theta_j) = \max_{\theta_{-j}}\ell(\theta_j,\theta_{-j}), \qquad \Delta(\theta_j) = 2\left[\ell(\hat{\theta})-\ell_{\mathrm{prof}}(\theta_j)\right].

Here, \(\theta_{-j}\) denotes all parameters other than \(\theta_j\), and \(\hat{\theta}\) is the maximum-likelihood estimate. Re-optimizing the other parameters is essential: it shows how well the model can fit when the parameter of interest is constrained, while allowing the remaining parameters to compensate.

Under regularity conditions, a likelihood-ratio interval for one parameter includes values for which \(\Delta(\theta_j)\leq 3.84\), the approximate 95th percentile of a chi-square distribution with one degree of freedom. The approximation may be unreliable near parameter boundaries, under strong nonregularity, or when the model is weakly identified; in those cases, calibration by simulation may be needed. A profile that rises clearly on both sides of its minimum supports a bounded estimate. A flat profile, a one-sided profile, or a long shallow valley indicates weak information or a boundary-limited estimate.

Worked Example: Interpreting a Profile for Internalization

In a hypothetical fit to drug and target observations, the estimated \(k_{\mathrm{int}}\) is \(0.50\ \mathrm{day}^{-1}\). At each fixed value below, all other parameters are re-optimized. The resulting profile deviances are: \(k_{\mathrm{int}}=0.10\), \(\Delta=4.8\); \(0.25\), \(\Delta=1.6\); \(0.50\), \(\Delta=0\); \(1.00\), \(\Delta=1.1\); and \(2.00\), \(\Delta=4.2\), with rates in \(\mathrm{day}^{-1}\).

Using the approximate 95% cutoff of \(3.84\), the values \(0.25\), \(0.50\), and \(1.00\ \mathrm{day}^{-1}\) are inside the profile interval, while \(0.10\) and \(2.00\ \mathrm{day}^{-1}\) are outside it. Linear interpolation between the two nearest values on each side gives approximate crossings. On the lower side, the fraction of the way from \(\Delta=4.8\) at \(0.10\) to \(\Delta=1.6\) at \(0.25\) is \((4.8-3.84)/(4.8-1.6)=0.300\); the crossing is \(0.10+0.300(0.15)=0.145\ \mathrm{day}^{-1}\). On the upper side, the fraction from \(\Delta=1.1\) at \(1.00\) to \(\Delta=4.2\) at \(2.00\) is \((3.84-1.1)/(4.2-1.1)=0.884\); the crossing is approximately \(1.88\ \mathrm{day}^{-1}\).

The approximate interval is therefore \(0.145\) to \(1.88\ \mathrm{day}^{-1}\), conditional on this model, observation model, and likelihood-ratio approximation. It is broad relative to the point estimate: the data support a range of internalization rates, not a finely determined rate. The interval is not a 95% probability that the fixed parameter lies within these particular endpoints.

Simulation-Based Parameter Recovery

A complementary diagnostic is to simulate data from the fitted model under the actual observation schedule and measurement-error assumptions, then fit each simulated data set using the intended estimation procedure. This asks whether the analysis can recover known parameter values under conditions resembling the data. It can reveal strong parameter correlations, frequent boundary estimates, failed fits, or large variation in estimates even when a single fit appears well behaved.

Simulation is not a substitute for defining a plausible observation model or checking the fit to observed data. It evaluates recovery under the assumed model and simulation conditions; it cannot establish that those assumptions are true. It is most informative when the simulated parameter values and error distributions are stated, and when recovery is assessed for each parameter rather than only by whether the optimizer converges.

Specify the fitted model and observation assumptions.
State which drug and target outputs are measured, the residual-error distributions, and the parameter bounds used in fitting.

Generate replicate data sets.
Simulate model predictions at the observed conditions and add measurement error according to the stated observation model.

Refit and summarize recovery.
Record convergence, bias, spread, boundary estimates, and correlations among parameter estimates. Persistent trade-offs indicate that the data-generation and fitting setup do not reliably separate those parameters.

For example, if many simulated fits recover \(k_{\mathrm{on}}\) and \(k_{\mathrm{off}}\) with a strong negative correlation while predictions remain similar, those rates are trading off in the fitted model. The correlation is evidence of a practical problem, but it does not by itself prove structural nonidentifiability. A profile likelihood can show whether the trade-off leaves a broad range of likelihood-supported values.

Reading Diagnostics and Avoiding Overinterpretation

Use more than one signal when deciding whether a parameter is well determined. A flat or asymmetric profile, a confidence interval extending over a large range, estimates repeatedly reaching parameter bounds, or large variability across simulated refits all point to weak practical identifiability. A narrow local standard error from an optimizer is not conclusive: it may rely on a local quadratic approximation that fails to describe a long curved ridge in the likelihood.

A useful distinction is between a poor fit and weak identifiability. Poor fit means the model predictions do not adequately describe the observations under the chosen error model. Weak identifiability means multiple parameter values give comparably good fits. The issues can occur together, but they are not interchangeable; a good fit does not guarantee uniquely estimated mechanisms.

A common mistake is to treat a converged fit with plausible parameter values as evidence that every mechanism has been identified. Convergence only indicates that the numerical procedure met its stopping criteria. If \(k_{\mathrm{int}}\) can vary substantially while the other rates re-optimize and the likelihood barely changes, the data do not support a narrow mechanistic conclusion about internalization, even if the optimizer returns several decimal places.

The main alternatives provide different views rather than interchangeable guarantees. The local covariance matrix or standard errors are computationally convenient but describe uncertainty near the optimum and can miss nonlinear ridges. Profile likelihood directly examines the fit over a parameter range, at greater computational cost. Simulation-based recovery tests the behavior of the whole fitting procedure under specified conditions, but conclusions remain conditional on the simulation model. No diagnostic can compensate for a misspecified model or an observation model that does not represent what the assay measures.

Key takeaway: Interpret a binding, internalization, or turnover estimate as mechanistically resolved only when the model and observations constrain it. Check profile shape and parameter-recovery behavior, and report broad or unbounded uncertainty as limited information rather than as a precise rate.

Check Your Understanding

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

A drug-free baseline target concentration is \(80\ \mathrm{pM}\). Give two positive synthesis and degradation rate pairs that predict this baseline, and state what the baseline identifies.
Under the quasi-steady complex approximation, what parameter combination appears in the denominator of the expression for \(C\)? Why can one complex measurement fail to separate the rates?
For a one-parameter profile likelihood, what does a broad region with \(\Delta\leq3.84\) indicate, under the usual regularity approximation?
A model fit converges, but simulation refits show large variation in \(k_{\mathrm{int}}\) and a strong correlation between \(k_{\mathrm{on}}\) and \(k_{\mathrm{off}}\). Does convergence establish that all three rates are well identified?

Solutions. (1) For example, \(k_{\mathrm{syn}}=8\ \mathrm{pM/day}\) and \(k_{\mathrm{deg}}=0.10\ \mathrm{day}^{-1}\) give \(8/0.10=80\ \mathrm{pM}\); \(16\ \mathrm{pM/day}\) and \(0.20\ \mathrm{day}^{-1}\) give the same value. Baseline identifies the ratio \(k_{\mathrm{syn}}/k_{\mathrm{deg}}\), not the rates separately. (2) The denominator is \(k_{\mathrm{off}}+k_{\mathrm{int}}\), and the observed complex also depends on \(k_{\mathrm{on}}\); different rate combinations can preserve the predicted \(C\). (3) It indicates a broad range of parameter values compatible with the data under the approximate likelihood-ratio criterion, so practical precision is limited. (4) No. Convergence is a numerical outcome, not evidence of identifiability; the simulation variability and rate correlation indicate that the fitted data may not distinguish the processes.