Identifiability

When designing a model, we must take care to ensure that there is a unique set of parameters that is better (i.e. has a higher likelihood or lower cost) than all other nearby solutions. If this is the case, the model parameters are said to be identifiable.

If two different sets of parameters can produce an identical fit to the data, and hence an identical likelihood, the model is unidentifiable. (Give a daft example)

Flip-flop is actually something completely different here from the usual meaning of “flip flop” and doesn’t technically agree with what we have written. Specifically, we identify two distinct parameter sets that give the same answer when, in fact, there is an infinite continuum of solutions that give the same answer. Flip-flop kinetics, in contrast, is where the two curves look similar (but not identical) but where the limiting factor in the drug’s reduction is KA or KE, depending on which flip-flop state you are considering.

Generally, avoid this flip-flop example until we’ve sorted out the confusion.

Flip Flop

Wikipedia has a page about Flip Flop Kinetics here:-

https://en.wikipedia.org/wiki/Flip-flop_kinetics

The flipping refers to the difficulty in measuring absorption and elimination rates in PK compartment models.

(We can this tie in with uncertainty and standard errors that have been introduced - you should get a valley in the (x1,x2) error space, resulting in a singular f[X] covariance matrix.)

Identifiability in a one compartment model

Consider one compartmental model with first order absorption.

For dose $AMT$ , absorption rate constant $K_a$ , elimination rate constant $K_e$ and volume $V$ , the resulting plasma concentration profile $C(t)$ is given by:-

() $C(t) = \frac{K_a \times \text{AMT}}{(K_a-K_e)V} \left[ e^{-K_e t} - e^{-K_a t} \right]$

Equation () has three independent variables ( $K_a$ , $K_e$ and $V$ ) which determine the concentration profile $C(t)$ .

Note, however, the equation () rearranges to:-

() $C(t) = \frac{K_e \times \text{AMT}}{ (K_a-K_e) \frac{ V K_e}{K_a} } \left[ e^{-K_e t} - e^{-K_a t} \right]$

Equation () is identical to equation () if:-

() $K_a \to K_e \\ K_e \to K_a \\ V \to \frac{V K_e}{K_a}$

Therefore, the parameters ( $K_a$ , $K_e$ and $V$ ) do not uniquely determine the concentration-time profile as a different set of parameters ( $\alpha$ , $\epsilon$ and $\nu$ ) can be constructed which would give exactly the same profile too:-

Solution	$K_a$	$K_e$	$V$
1	$\alpha$	$\epsilon$	$\nu$
2	$\epsilon$	$\alpha$	$\frac{\nu \epsilon}{\alpha}$

Flip Flop can create a problem when fitting a model as there are two equally valid solutions to the same problem. One solution will have a ka which is higher than ke and the other a ka which is lower than ke. The values of ke and ka are interchangeable and the solution reached depends on the initial parameters used in the model.

Flip Flop Good Fit

See Flip Flop tutorial with low initial estimates of KA, V and CL for example using PoPy.

The ground truth generating parameters put into the tutorial script under the GEN_EFFECTS->POP section are:

EFFECTS:
    POP: |
        c[ID] = 1
        c[AMT] = 100.0
        t[RESET] = 0.0
        t[DOSE] = 1.0
        t[OBS] ~ unif(1.0, 40.0; 100)
        f[KA] = 0.15
        f[KE] = 0.3
        f[ANOISE_STD] = 1.0

The initial parameters used in the tutorial script under the FIT_EFFECTS->POP section are:

EFFECTS:
    POP: |
        f[KA] ~ unif(0.001, 1) 0.5
        f[KE] ~ unif(0.001, 10) 0.5
        f[ANOISE_STD] ~ unif(0.001, 10) 5.0

Note in the section above the initial values of KE is much less than the initial values of KA, as is the case with the ground truth data.

The final fitted parameters are:

f[KA] = 0.1526
f[KE] = 0.3014
f[ANOISE_STD] = 0.9207

The final objective function value was

83.4768484306088

The results are summarised in the table below.

Name	Initial	Fitted	True	Abs. Error	Prop. Error
f[KA]	0.5	0.153	0.15	2.57e-03	1.72%
f[KE]	0.5	0.301	0.3	1.41e-03	0.47%
f[ANOISE_STD]	5	0.921	1	7.93e-02	7.93%

You can see that the final fitted values (third row) are close to the original ground truth (first row).

Flip Flop Bad Fit

DERIVATIVES: |
    conc = s[CENTRAL]/c[V]
    d[CENTRAL] = @bolus{amt:c[AMT], lag:0} - c[CL]*conc
    d[DIS] = exp(m[ALPHA]) - m[BETA]*conc

figure:: /case-studies/indiv_examples/disease_prog/dp_exponential/images/fit_sim_grph_outputs/indOBS_vs_TIME/000001.* :width: 75% :align: center

See Flip Flop tutorial with inaccurate starting values for example using PoPy.

The ground truth generating parameters are the same as in section Flip Flop Good Fit. Note that we used:

METHOD_OPTIONS: {py_module: tut, rand_seed: 314159}

In both tutorial scripts, this ensures that the seed used to generate random data is the same, so both tutorials have exactly the same ground truth data.

The initial parameters used in the bad fit tutorial script under the FIT_EFFECTS->POP section are now:

EFFECTS:
    POP: |
        f[KA] ~ unif(0.001, 1) 0.1
        f[KE] ~ unif(0.001, 10) 1.0
        f[ANOISE_STD] ~ unif(0.001, 10) 5.0

Note in the sections above the initial values of KE is greater than the initial values of KA. Where as KE is less than the ground truth synthetic data.

The final fitted parameters this time are:

f[KA] = 0.1526
f[KE] = 0.3014
f[ANOISE_STD] = 0.9207

The final objective function value was

83.47684854744028

The results are summarised in the table below.

Name	Initial	Fitted	True	Abs. Error	Prop. Error
f[KA]	0.1	0.153	0.15	2.58e-03	1.72%
f[KE]	1	0.301	0.3	1.41e-03	0.47%
f[ANOISE_STD]	5	0.921	1	7.93e-02	7.93%

The final fitted values (third tow) have flipped with the value of ka being close to the initial value (first row) of ke used to generate the data and vice versa.

$K_a = 0.6 \approx 0.59833 (\text{Fitted}\enspace K_e) \\ K_e = 0.05 \approx 0.048954 (\text{Fitted}\enspace K_a)$

Also note that the flipped Volume of distribution is computed (from the ground truth parameters) as follows.

$V = \frac{V \times K_e}{K_a} = 40*0.05/0.6 = 3.33 \approx 3.23 (\text{Fitted}\enspace V)$

Note that we are comparing the f[X] population variables here, which will agree only approximately with the theoretical flipped results based on individual compartment models.

Conclusion

[THIS IS OUT OF DATE AND NEEDS REVISING]

The Flip Flop Bad Fit example is a worse result having an objective function of 123.162477043, compared to -281.748419579 with the Flip Flop Good Fit example.

The fitting in PoPy (and other PKPD software) is a local optimisation process. Therefore as with all local optimisation there is a risk of reaching false minima, if the initial starting values are incorrect.

Flip flop models are prone to such false minima on the likelihood surface, however the initial values have to be very wrong to end up in the flipped minima.

Overparameterization

When fitting a model, it is important to ensure all the parameters are determinable.

The compartment diagram is:

The parameters in this model are derived as follows:

$\text{ka} &= \text{Absorption rate from depot compartment to central compartment} \\ \text{k} &= \frac{\text{CL}}{\text{V1}} \\ \text{k12} &= \frac{\text{Q}}{\text{V1}} \\ \text{k21} &= \frac{\text{Q}}{\text{V2}} \\ \text{CL} &= \text{Clearance from Central compartment} \\ \text{V1} &= \text{Volume of distribution in Central compartment} \\ \text{Q} &= \text{Inter-compartmental Clearance} \\ \text{V2} &= \text{Volume of distribution in Peripheral compartment} \\ \text{dur} &= \text{Duration of infusion}$

The parameters used in the gen script are:

$\text{ka} &= 0.3 \\ \text{CL} &= 2 \\ \text{V2} &= 10 \\ \text{Q} &= 1.5 \\ \text{V3} &= 50 \\ \text{dur} &= 20$

The estimated parameters from the fit script are:

$\text{ka} &= 0.17 \\ \text{CL} &= 1.84 \\ \text{V2} &= 17.19 \\ \text{Q} &= 1.18 \\ \text{V3} &= 42.65 \\ \text{dur} &= 8.59 \\ \text{objective function} &= -2056$

The PoPy fit script does not return the values of the parameters used in the generation script.

This is because ka and dur both affect the rate of the drug entry into the central compartment. If we look at the compartment diagram for this model, we can see that infusion duration affects the rate into the DEPOT compartment and ka affects the rate between the DEPOT and CENTRAL compartments. However, all our concentration measurements are taken in the CENTRAL compartment so we have no information on the different effects of these two parameters. Therefore, it is not possible for the fit script to separate these two effects, so the two parameters are not determinable.

If we think about why there would be differences in the duration of infusion biologically, then the reason for this becomes clearer. Differences in duration could arise because the drug is infused into the body at different rates, but this would be set by the clinician and it should not be necessary to model this. However, the infusion may not be straight into the central compartmental and the different times taken to reach the central compartment in different individuals are what we are modelling. This could be regarded as analagous to a rate of absorption, which also varies in individuals due to the time the drug takes to reach the central compartment.

The solution to this problem would be to model either ka or dur, but not both.

Bioavailability and Lag

As an example, we consider a one compartment model with first order absorption (see One Compartment Model with Absorption) that has two additional parameters: bioavailability and a lag.

Bioavailability

If a drug is administrated extravascularly, it is not guaranteed that the whole amount will reach the plasma. Bioavailability is the fraction of the drug that does reach the Central compartment.

If a drug is taken orally, for example, not all of the drug may be absorbed due to the first pass effect, where the drug concentration is reduced in either the digestive tract or liver, before reaching blood circulation. This will reduce the amount of drug entering the plasma and therefore the bioavailability.

Bioavailability can vary both between individuals and in the same individual over time, and it is important to include it in a PK/PD model when this is the case. Because it is a fraction, bioavailability can take a value ranging from 0 (no drug reaches the Central compartment) and 1 (all drug reaches the Central compartment).

In PoPy, bioavailability can be introduced as a model parameter in the MODEL_PARAMS section (e.g. m[BIO]) and multiplied by the dose in the DERIVATIVES section:

d[DEPOT] = @bolus{ amt: c[AMT]*m[BIO] } + ...

Visually, a change in bioavailability scales the Amount Administered vs Time curve vertically (table_bioavailability).

Two curves with different bioavailabilities: (left) BIO=1.0 and (right) BIO=0.583. (Note the difference in scale on the y-axis.) All observations sampled at identical time points.

Lag Time

Lag time is the time taken for a drug to appear in the system following administration. This can vary between individuals or in the same individual over time (or both).

In PoPy, lag time can be introduced as a model parameter in the MODEL_PARAMS section (e.g. m[LAG]) and then used as a dose parameter in the DERIVATIVES section:

d[DEPOT] = @bolus{ amt: c[AMT], lag: m[LAG] } + ...

Visually, a change in the lag shifts the Amount Administered vs Time curve horizontally (table_lag).

Two curves with different lags: (left) LAG=0.0 and (right) LAG=5.0. All observations sampled at identical time points.

This example - a one compartment model with IV administration - makes no sense. BIO would be 1.0 (by definition) and LAG would be 0.0 (by definition).

Bioavailability and Lag Time

Bioavailability and Lag is an example of a one compartment model with absorption containing both bioavailability and lag time.

Assuming a bolus dose, the compartment diagram is:

../_images/compartment_diagram25.svg — Compartment diagram for a model with bioavailability and a time lag

The MODEL_PARAMS section is:

MODEL_PARAMS: |
    m[CL] = f[CL]
    m[V] = f[V]
    m[ANOISE_STD] = f[ANOISE_STD]
    m[BIO] = f[BIO]
    m[LAG] = f[LAG]

And the DERIVATIVES section:

DERIVATIVES: |
    d[CENTRAL] = @bolus{amt:c[AMT]*m[BIO], lag: m[LAG]}  # dose in
    d[CENTRAL] -= m[CL]/m[V]*s[CENTRAL]  # elimination out

Confounded Bioavailability and Lag Time

This tutorial outlines some problems that can occur when bioavailability and lag time are both included in a model.

This can happen when a bolus dose feeds directly into the measurement compartment of a model containing both bioavailability and a lag time.

A diagram of the most simple, one compartment model, is shown below:

The amount of dose reaching the measurement compartment is affected by both BIO and LAG.

Bioavailability (BIO) is the fraction of the dose reaching the compartment. For example, if a dose was 100mg and bioavailability was 0.8, then only 80mg would reach the compartment.

A lag time (LAG) is a delay in the dose reaching the compartment. For example, if the LAG was 1 hour, then the dose would not arrive in the compartment until one hour after the dose was administered.

In this case, the amount of the drug in the body at time t, has a closed form solution as follows:

$A(t) = \Biggl \lbrace { F \times D \times \exp (-k(t-L)) ,\text{ if } { t \geq L } \atop 0.0, t < L }$

Where:

Variable	Description
A(t)	Amount at time t
F	Bioavailability
D	Dose
k	Rate constant
t	Time
L	Time lag

When L = 0.0, F = 1.0, K = 0.5 and D = 100mg, the time course of the drug in the body is shown in Figure (a). Introducing a bioavailability of 0.8 reduces the height of the curve by 20%, see Figure(b). Introducing a lag time of 1.0 would move the curve to the right, see Figure (c). Both a lag of 1.0 and bioavailability of 0.8 are added in Figure(d).

Now consider a Figure (e) with a lag time of 0.83 and a bioavailability of 0.9. This includes a measurement after 0.83 time units. If the first measturement is at 1.0, then you have Figure (f).

In the Figures (a-d), the first measurement is after 1.0 time unit. The plotted points after time 1.0 in Figure (f) are identical to the points plotted in Figure (d). Figure(d) has bioavilability 0.8 and lag 1.0, but Figure (f) has bioavailability 0.9 and the lag time 0.83.

If we were trying to solve this problem and only had measurements for time points after 1.0 an answer of bioavailability of 0.8 and a lag of 1 or an answer of bioavailability of 0.9 and a lag of 0.83 are equally valid.

Other estimates of both these parameters would also be possible. For this example, bioavailability could range between 0.8 and 1 and lag time between 0.678 and 1. As the estimate of bioavailability increases, the estimate of time lag decreases. Mathematically, bioavailability be anything over 0.8 and lag time anything under 1, but biologically this would be nonsense and bioavailablity is limited to a maximum of 1 (when lag time = 0.678).

The problem occurs because the bolus dose suddenly enters the central compartment between the time and units of zero and one, but we have no information on the amount of drug in the body between the time units of zero and one.

This makes it impossible to determine accurate estimates of bioavailability and lag time.

PoPy Example of a Confounded Model

The DERIVATIVES section in PoPy is shown below:

DERIVATIVES: |
    d[CENTRAL] = @bolus{amt: c[DOSE]*m[BIO], lag: m[LAG]} - m[K]*s[CENTRAL]

If we examine the tutorial script we can see that the ground truth values for BIO and LAG were 0.8 and 1.

The parameter estimates from the fit script are BIO = 0.623 and LAG = 1.510.

Clearly PoPy was unable to recover the original parameters. The fit script also pushes up the additional noise to the maximum we specified (1.0). This should alert us to the fact that something went wrong with the model.

It should be noted that this is not due to a bug in PoPy. If the data generated by this tutorial is fitted using NONMEM, the ground truth values of BIO and LAG are not recovered either.

This plot, though very illustrative, represents the relationship between BIO and LAG for a one compartment model with IV administration. However, with this model BIO=1 and LAG=0 by definition. Therefore these parameters would never actually vary for this structure.

Confounding between BIO and LAG

The Bioavailability and Lag tutorial demonstrates, using a one compartment model with first-order absorption into the Central compartment, how bioavailability and lag interact.

Because not all of the drug will be absorbed by the Depot compartment, and some will be metabolized due to the first pass effect, the amount of drug will be a fraction, m[BIO], of the amount administered.

Furthermore, because of the time taken to absorb the drug its effect will be seen some time, m[LAG], after it was administered.

../_images/compartment_diagram5.svg — Compartment diagram for a one compartment model with first order absorption, plus variable bioavailability and lag.

These new parameters are included in the dose line of the DERIVATIVES section in the PoPy script:

DERIVATIVES: |
    d[DEPOT] = @bolus{amt:c[AMT]*m[BIO], lag: m[LAG]}  # dose in
    d[CENTRAL] = 0.0
    d[DEPOT->CENTRAL] += m[KA]*s[DEPOT]
    d[CENTRAL] -= m[CL]/m[V]*s[CENTRAL]  # elimination out

Looking at the cost surface for these two parameters, m[LAG] and m[BIO], we see that there is a strong relationship between the two.

Specifically, we see a ‘valley’ running through the cost surface where many (LAG,BIO) pairings give an almost equally low cost. This is because for the simulated data we have, where the absorption phase of the curve is not observed, scaling (via BIO) and simultaneously shifting (via LAG) the curve can generate similar observations (table_biolag).

Two curves with different parameters: (left) LAG=0.0, BIO=1.0 and (right) LAG=5.0, BIO=0.583. (Note the difference in scale on the y-axis.) All observations sampled at identical time points, and the observed concentrations are very similar.

As a result, for these data the bioavailability and lag parameters are hard to estimate with any certainty.

Weibull Dosing Parameters

New tutorial that shows that Weibull shape parameters can easily be confounded with the LAG parameter. Good demonstration of the Weibull dosing function (that NONMEM does not have).