NONMEM Users Network Archive

Hosted by Cognigen

RE: Code to avoid flip-flop kinetics

From: Jurgen Bulitta <jbulitta>
Date: Mon, 10 Aug 2009 19:54:42 -0400

Dear Nick,
Dear Paolo,
Dear Leonid,
Dear All,

Depending on the application, I think there are the following reasons for c=
onstraining ka and ke:

1) Standard-two-stage (STS) approach: As discussed previously, calculating =
average and SD from STS estimates is likely to yield biased averages and SD=
s and this bias may have a great impact when these descriptive statistics a=
re used as means and variances in a Monte Carlo simulation. Thus, manually =
switching of individual ka & ke (if needed) makes perfect sense for STS.

2) Nonparametric population PK with nonparametric simulation: As the interi=
ndividual distribution of PK parameters can take any shape, it will make no=
 difference for a non-parametric simulation based on the list of support po=
ints whether ka and ke are constrained or not. Even if a support point (wit=
h e.g. 10% probability and ka>ke) splits into two equivalent support points=
 with 5% probability each (5% for ka>ke and 5% for ke>ka) this will not aff=
ect the results of a nonparametric simulation. Thus for a nonparametric sim=
ulation, constraining is not needed.

3) Nonparametric population PK with parametric simulation: Once one compute=
s averages and SDs (etc.) from the list of support points, the same issue o=
f biased averages and SDs arises as for STS, if ka and ke are not constrain=
ed. As averages and SDs are commonly calculated for nonparametric pop PK mo=
dels, constraining at the individual parameter estimate level seems most re=
asonable, as proposed by Vladimir Piotrovsky (option 2 in Paolo's email). I=
 think this was the approach taken at the 2005 PAGE software comparison.
I am not aware of a method of constraining a nonparametric model at the pop=
ulation level.

4) Parametric EM-algorithms: Assume true mean ka and true mean ke are withi=
n 10% of each other and BSV has CV of 30% each. Clearly, there must be subs=
tantial overlap in the true individual estimates of ka and ke.

4.1) As one EM algorithm, the MC-PEM algorithm assumes that the interindivi=
dual distribution of individual model parameters is normal (or lognormal, o=
r can be transformed to be normal). Otherwise, the population mean and popu=
lation variance-covariance matrix obtained from the conditional means and c=
ondition variance matrices of each subject (see eq. 23 and 27 in [1] and re=
f. [2]) does not necessarily minimize the overall log-likelihood. Truncatio=
n of the interindividual distribution of PK parameters, for example by use =
of the EXIT command (option 3 in Paolo's email), seems to violate the norma=
l distribution assumption.

4.2) The conditional means calculated in the expectation-step can be biased=
 due to e.g. a bi-modal distribution caused by the flip-flop (see also Serg=
e's comment at the end of the discussion in 2003). This might be a more sev=
ere concern compared to 4.1).
If such bias in the conditional means and in the conditional variance matri=
ces occurs, the computation of the population means and population variance=
-covariance matrix via eq. (21) and (22) (see Bauer et al. [1]) does no lon=
ger guarantee to minimize the objective function, I think. Constraining ka =
and ke at the population level (option 1 in Paolo's email) is unlikely to s=
olve this problem, if mean ka and mean ke are close.

If one does not constrain ka and ke at the individual parameter level, bias=
 of conditional means is likely to occur, if:
a) the mean ka and mean ke are close given their interindividual variabilit=
b) the data are sparse and one therefore uses a wide envelope function to a=
pproximate an individual's conditional density (this problem is more likely=
 to occur if one uses low values of gefficiency [please see S-ADAPT manual]=
 or if one uses e.g. a t-distribution with small degrees of freedom as enve=
lope function)

Therefore, constraining at the individual parameter level via estimation of=
 CL, V, and (ka-kel) (i.e. the Vladimir Piotrovsky method) and use of a ful=
l var-cov matrix (as proposed by Leonid) seems the best choice for the MC-P=
EM algorithm. I think this also applies to the SAEM and MCMC algorithm. The=
se algorithms usually have no problem with estimating a full var-cov matrix=
 and do not require excessively longer runtimes to estimate a full var-cov =
matrix, as opposed to FOCE.

1) There is no need to constrain ka and ke for a nonparametric estimation w=
ith nonparametric simulation. However, as constraining does not hurt, I wou=
ld still constrain ka and ke, since interpretation of mean parameters is ea=
2) For nonparametric estimation, ka and ke should be constrained at the ind=
ividual level, if means and SDs (etc.) are to be computed from the list of =
support points (as commonly done).
3) Constraining ka and ke at the individual level and estimation of a full =
var-cov matrix seems most appropriate for parametric EM algorithms and for =
nonparametric estimation with subsequent parametric Monte Carlo simulations=

Hope these comments are helpful despite their length.
Nick, I would agree that in many cases, fortunately the issues mentioned ab=
ove only have a minor practical impact.

Best wishes

[1] Bauer RJ et al. AAPS Journal, 2007, 9(1): E60-E83. See bottom left and =
right column on page E64.

[2] Schumitzky A . EM algorithms and two stage methods in pharmacokinetics =
population analysis. In: D'Argenio DZ , ed. Advanced Methods of Pharmacokin=
etic and Pharmacodynamic Systems Analysis. vol. 2. Boston, MA : Kluwer Acad=
emic Publishers ; 1995 :145- 160.

Jurgen Bulitta, Ph.D., Senior Scientist
Ordway Research Institute, Albany, NY, USA

-----Original Message-----
From: owner-nmusers
 Behalf Of Nick Holford
Sent: Friday, August 07, 2009 5:03 PM
To: nmusers
Subject: Re: [NMusers] Code to avoid flip-flop kinetics


I agree with Leonid -- especially this:
> If population K and KA are close, then it is unclear why individual K
> and KA should be consistently ordered as K < KA and not vice versa for
> some subjects. If so, why not to allow the model to decide whether to
> have flip flop or not?
I do not understand what the concern is with flip-flop kinetics and why
people want to force K<KA. If you have a good reason then various
methods are available (as reviewed again by Paulo).

I wonder if someone would like to describe why it is of interest to
force K<KA?


Leonid Gibiansky wrote:
> Paolo,
> Option 1 is very helpful. Option 2 is not attractive for the reason
> that you stated, especially, extra correlation introduced by this
> trick. If used, it should be used with the full OMEGA block (CL-V-KA).
> I used it once but only as the last resort when nothing else worked. I
> have not tried option 3 but this option artificially restricts
> distributions, so I am not sure whether it is good or not even when it
> is working.
> On the other hand, if population K and KA are sufficiently far apart,
> you are unlikely to get individual K and KA flip-flopped.
> If population K and KA are close, then it is unclear why individual K
> and KA should be consistently ordered as K < KA and not vice versa for
> some subjects. If so, why not to allow the model to decide whether to
> have flip flop or not?
> I would try to use option (1), and if you like the model, diagnostic
> plots, etc, I would not worry about individual K and KA relation. One
> of the diagnostics could be the fraction of patients with the
> flip-flop. If it is small, this would justify the approach. Another
> possible diagnostic is ETA_KA - ETA_K differences, If it is close to
> normal you are OK, if it has two mirror-symmetric (relative to the
> y-axis) peaks, then flip flop is interfering with the model.
> Thanks
> Leonid
> --------------------------------------
> Leonid Gibiansky, Ph.D.
> President, QuantPharm LLC
> web:
> e-mail: LGibiansky at
> tel: (301) 767 5566
> Paolo Denti wrote:
>> Dear NMUsers,
>> having to deal with the flip-flop kinetics phenomenon, I had a look
>> at previous posts on the NMusers list.
>> I found this post particularly enlightening:
>> Some code was proposed to avoid the flip-flop at population and
>> individual level. Here's a not-so-brief summary.
>> This parameterization solves the issue at population level:
>> CL=THETA(1)*EXP(ETA(1))
>> V=THETA(2)*EXP(ETA(2))
>> However, it does not prevent the phenomenon occurring at individual
>> level. Vlamidir Piotrovsky proposed the code below, which does solve
>> the problem at individual level, but makes the interpretation of the
>> results a bit awkward and introduces correlation among the model
>> parameters. In particular, variance of ETA3 was greatly increased.
>> CL=THETA(1)*EXP(ETA(1))
>> V=THETA(2)*EXP(ETA(2))
>> KE=CL/V
>> Another approach, suggested by Nick Holford, implements error
>> recovery using EXIT 1. The code is reported below:
>> CL=THETA(cl)*EXP(ETA(cl))
>> V=THETA(v)*EXP(ETA(v))
>> KA=THETA(ka)*EXP(ETA(ka))
>> K=CL/V
>> IF (KA.LE.K) EXIT 1 101 ; try again (PREDERR message error code 101)
>> As far as I understand, this interrupts the computation whenever the
>> flip-flop occurs and lets NONMEM restart. However, if such an error
>> arises at initialization, NONMEM does not recover and the run goes no
>> further. Nick probably experienced something similar, but apparently
>> received no answer
>> Does anyone know of a way around this drawback? Or have other code to
>> deal with flip-flop kinetics?
>> Thank you in advance,
>> Paolo
>> --
>> ------------------------------------------------
>> Paolo Denti, Post-Doctoral Fellow
>> Division of Clinical Pharmacology
>> Department of Medicine
>> University of Cape Town
>> K45 Old Main Building
>> Groote Schuur Hospital
>> Observatory, Cape Town
>> 7925 South Africa
>> phone: +27 21 404 7719
>> fax: +27 21 448 1989
>> email: paolo.denti
>> ------------------------------------------------

Nick Holford, Professor Clinical Pharmacology
Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealan=
mobile: +64 21 46 23 53
Received on Mon Aug 10 2009 - 19:54:42 EDT

The NONMEM Users Network is maintained by ICON plc. Requests to subscribe to the network should be sent to:

Once subscribed, you may contribute to the discussion by emailing: