# Comparing to NONMEM (bioav & lags)

**URL:** <https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494>\
**Category:** Simulation\
**Created:** [August 3, 2021, 5:59am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494 "2021-08-03T05:59:57Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 3, 2021, 5:59am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/1 "2021-08-03T05:59:57Z")

</div>

Hello there! I am new to Pumas and have trying to validate some of my previous results from NONMEM with Pumas. I have read through the docs nearly 10 times, and I got some of my other models to work, but for some reason, one of my drugs is producing results that are consistently too low. I cannot seem to identify the structural difference between the models that could be causing this, and I have tested and ensured that the error model is not the cause.

Here is the NONMEM model:

```auto
$PROBLEM Simulation Levofloxacin PK 
$INPUT ID TIME EVID AMT DV TAD DROP SEX AGE WT HT WAZW WHZ HAZW BAZW COUNTRY ; Dropped TAD HIV ... Country    
$DATA ../pop_sample1000_test.csv IGNORE=@
; AMT in mg
; DV in ug/mL or mg/L
; TIME in hours

$SUBROUTINE ADVAN6 TOL=5

$MODEL NCOMP=4
            COMP=(DEPOT DEFDOS) 
            COMP=(CENTRAL DEFOBS)
            COMP=(PERIPH) 
            COMP=(AUC)

$PK
BSIZE1 = WT/12 ; This needed to be changed to 12 kg based on the study. see footnote a in table 2. 
EXPO = 0.75

CLHIV = 1 + 0*THETA(8)

; AGE EFFECT 
GAM = THETA(11)
PMAGE50 = THETA(10)
MAT = AGE**GAM/(PMAGE50**GAM + AGE**GAM) ; AGE in months 

; START PK
TVCL = THETA(1)*CLHIV*MAT*BSIZE1**EXPO ; added maturation on CL, removed HAZ.
CL = TVCL*EXP(ETA(1)) ; CL, L/h [systemic clearance]

TVV = THETA(2)
V = TVV*BSIZE1*EXP(ETA(2)) ; V, L [central volume]

TVKA = THETA(3) ; KA, 1/h [absorption rate]
KA = TVKA*EXP(ETA(3))
 
TVVP = THETA(6)
VP = TVVP*BSIZE1*EXP(ETA(4))

TVQ = THETA(7)
Q = TVQ*BSIZE1**EXPO*EXP(ETA(5))

TVF1 = 1
F1 = TVF1*EXP(ETA(6))

K20 = CL/V
K23 = Q/V
K32 = Q/VP

ALAG1 = THETA(9)*EXP(ETA(7))

$DES

DADT(1) = -KA*A(1)
DADT(2) = KA*A(1)-K20*A(2)-K23*A(2)+K32*A(3) 
DADT(3) = K23*A(2)-K32*A(3)
DADT(4) = A(2)/V

$ERROR 

REP = IREP
AUC = A(4)

; Proportional Residual Error model
CP = A(2)/V
IPRED = CP 
PROP = THETA(4)
ADD = THETA(5)

W = SQRT(PROP **2*IPRED** 2+ADD**2)
IRES = DV - IPRED
IWRES = IRES/W
Y = IPRED + W*EPS(1)

$THETA ; you need to change CL, V, Prop, Add, VP, Q ...
 4.7 ; 1 CL,L/h [systemic clearance]
 19.2 ; 2 Vc,L [central volume]
 1.61 ; 3 Ka
 0.116 ; 4 PROP 
 0.0160 ; 5 ADD 
 2.40 ; 6 VP
 0.796 ; 7 Q
 -0.159 ; 8 HIV on CL, study found 15.9% decrease in CL
 0.242 ; 9 ALAG
 10.6 ; 10 PMAGE50
 3.39 ; 11 GAMMA

$OMEGA  
 0.023 ; CL IIV these need to be converted. 15% CV = .15^2 = 
 0 FIX ; V IIV 
 0.420 ; KA IOV ... we might want to consider inputting the BOV but for now it is fine.
 0 FIX ; VP
 0 FIX ; Q
 0.048 ; F1
 1.69 ; ALAG IOV

$SIGMA 1 FIX ; RV (use 20% CV as init est)
$SIM ONLYSIM NSUB=100 (12345)
$TABLE NOPRINT ONEHEADER FILE=simtab ID TIME EVID AMT TAD DV IPRED CP WT AUC REP WT AGE 

```

And here is the Pumas model;

```auto
model = @model begin
    @param begin
        tvcl ∈ RealDomain(lower = 0) ; 1
        tvv ∈ RealDomain(lower = 0) ; 2
        tvka ∈ RealDomain(lower = 0) ; 3
        σ_prop ∈ RealDomain(lower = 0) ; 4
        σ_add ∈ RealDomain(lower = 0) ; 5
        tvvp ∈ RealDomain(lower = 0) ; 6
        tvq ∈ RealDomain(lower = 0) ; 7
        clhiv ∈ RealDomain(lower = 0) ; 8
        alag ∈ RealDomain(lower = 0) ; 9
        pmage50 ∈ RealDomain(lower = 0) ; 10
        γ ∈ RealDomain(lower = 0) ; 11
        Ω ∈ PDiagDomain(7)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates WT AGE

    @pre begin
        BSIZE1 = WT / 12
        EXPO = 0.75
        
        CLHIV = 1 + 0 * clhiv
        
        VP = tvvp * BSIZE1 * exp(η[4])
        Q = tvq * (BSIZE1^0.75) * exp(η[5])
        MAT = (AGE^γ) / ((pmage50^γ) + (AGE^γ))
        
        CL = tvcl * CLHIV * MAT * (BSIZE1^0.75) * exp(η[1])
        V = tvv * BSIZE1 * exp(η[2])
        KA = tvka * exp(η[3])
        
        K20 = CL / V
        K23 = Q / V
        K32 = Q / VP

        tvf1 = 1
        
        bioav = (Depot = tvf1 * exp(η[6]),)
        lags = (Depot = alag * exp(η[7]),)
    end
    
    @dynamics begin
        Depot' = -KA*Depot
        Central' = KA*Depot - Central*K20 - Central*K23 + K23*Peripheral
        Peripheral' = K23*Central - K32*Peripheral
        AUC' = Central / V
    end

    @vars begin
        conc = Central / V
        auc = AUC
    end

    @derived begin
        cp = @. Central / V
        dv ~ @. Normal(cp, sqrt(σ_add^2 + (cp*σ_prop)^2))
    end
end

fixeffs = (
    tvcl = 4.7,
    tvv = 19.2,
    tvka = 1.61,
    σ_prop = 0.116,
    σ_add = 0.0160,
    tvvp = 2.4,
    tvq = 0.796,
    clhiv = -0.159,
    alag = 0.242,
    pmage50 = 10.6,
    γ = 3.39,
    Ω = Diagonal([0.023, 0, 0.420, 0, 0, 0.048, 1.69]),
)

obs = simobs(model, population, fixeffs, obstimes = 0:4:24) |> DataFrame

```

I apologize in advance if this post is too lengthy or otherwise out of the scope of what everyone can help with, but I appreciate any advice in advance, thank you!

Here is a histogram comparing the distributions of AUC between a single run of Pumas and a single run of NONMEM

 ![Screen Shot 2021-08-02 at 10.56.23 PM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/dbbdaab170ffb2669d52035de1e21c998cd5c8f9.png)

And here is a histogram comparing the median AUC of a single Pumas run to the distribution of median AUC of 100 simulation runs in NONMEM

 ![Screen Shot 2021-08-02 at 10.52.03 PM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/3b587a32c6cb20d1ecac56b47cbd9b7d14e2ef10.png)

Thank you again!

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [August 3, 2021, 9:22am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/2 "2021-08-03T09:22:30Z")

</div>

hi Dhruv,

Welcome to the community and thank you for the detailed message.

Could you please clarify the time vector of the simulation you set up for your NONMEM dataset.

Best,  
Vijay

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 5, 2021, 4:34am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/4 "2021-08-05T04:34:54Z")

</div>

Hi Vijay! Thank you for the prompt response.

The NONMEM dataset has dose times [0, 24] and obs times [0, 4, 8, 12, 16, 20, 24]

Here is the population constructor I used in Julia for Pumas (reads from a DataFrame of individuals containing demographic data in corresponding columns)

```auto
function dose(WT, AGE)
    function amt(WT, AGE)
        return WT * 20
    end
    return DosageRegimen(amt(WT, AGE); time = 0:24:24)
end

function population_constructor(subjects)
    population = Subject[]
    for row in eachrow(subjects)
        subject = Subject(;
            id = row.ID, 
            events = dose(row.WT, row.AGE), 
            covariates = (WT = row.WT, AGE = row.AGE)
        )
        push!(population, subject)
    end
    return population
end

population = population_constructor(subjects[1:1000, :])

```

---

<div class="post-metadata">

**Author:** ![benjaminrich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/benjaminrich/32/156_2.png) [@benjaminrich](https://discourse.pumas.ai/u/benjaminrich)\
**Post date:** [August 6, 2021, 12:17pm UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/5 "2021-08-06T12:17:32Z")

</div>

The first thing I would do is simulate a single subject with each model, setting all random effects to 0, and see if the simulated concentrations match up. Have you tried that? On inspection, the models do look the same to me (I don’t see any obvious errors).

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 7, 2021, 5:48am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/6 "2021-08-07T05:48:09Z")

</div>

Hi Benjamin!

This was very insightful advice. I completed the steps to set random effects to zero (set the Omega vector to all zeros in Pumas and fixed the Omega block to all zeros in NONMEM). The results seem to show that the concentration in Pumas decays much faster after the dose event as compared to NONMEM. This seems to be causing the Pumas AUC to be significantly lower than the AUC in NONMEM as referenced in my first post. However, I cannot seem to identify the component of my model in Pumas / NONMEM that could be causing this difference.

Here is the plot for comparison:

 ![Screen Shot 2021-08-06 at 10.42.15 PM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/73becc0fbe64128f739ca370c22a8d0fcd4a60ec.jpeg)

Here is the difference in AUC:

 ![Screen Shot 2021-08-06 at 10.46.57 PM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/c085895c8e0908c436f66f4c9d18e1a8b6d2c747.png)

And here is the difference in DV for reference:

 ![Screen Shot 2021-08-06 at 10.47.01 PM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/ed5df9e29372ccde3ac9ab99be26d486f21093a6.png)

I guess my question is: Is this due to some misspecification in my model or simply due to NONMEM and Pumas converging to different solutions?

---

<div class="post-metadata">

**Author:** ![benjaminrich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/benjaminrich/32/156_2.png) [@benjaminrich](https://discourse.pumas.ai/u/benjaminrich)\
**Post date:** [August 7, 2021, 11:11am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/7 "2021-08-07T11:11:50Z")

</div>

I think I’ve spotted the error. It’s in the Pumas model, in this line:

```auto
        Central' = KA*Depot - Central*K20 - Central*K23 + K23*Peripheral

```

The last term has `K23` where it should be `K32`. In other words, it should read:

```auto
        Central' = KA*Depot - Central*K20 - Central*K23 + K32*Peripheral

```

---

<div class="post-metadata">

**Author:** ![pkofod](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/pkofod/32/59_2.png) [@pkofod](https://discourse.pumas.ai/u/pkofod)\
**Post date:** [August 10, 2021, 7:22am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/8 "2021-08-10T07:22:16Z")

</div>

That does seem right @benjaminrich ! Good catch.

Let me add another way to check or verify this. You could use the `Depots1Central1Periph1` analytical solution and add on the AUC calculations as a numerical calculation

```auto
@dynamics Depots1Central1Periph1 begin
    AUC' = Central/V
end

```

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 10, 2021, 8:07am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/9 "2021-08-10T08:07:04Z")

</div>

Hi All!

Thank you very much for the advice. In fact, @benjaminrich that was the issue. I am now getting values much more closely in the range of what I would expect! Thank you also to @pkofod! I tried running the model with this preset, but after consulting the documentation, I determined that the built in specification for that model is different from how I have defined it in NONMEM (uses (CL + Q)/Vc for K20 whereas I am using CL/Vc

Here are the final results:  
**TIME vs. CP**

 ![Screen Shot 2021-08-10 at 1.03.53 AM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/71caf7c80810b7f4d8eb1c00c535ec0839e707a4.png)

**AUC Distribution**

 ![Screen Shot 2021-08-10 at 1.03.59 AM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/2e3310d39090d7dbc804578b3fb87e76622811b2.png)

**Median AUC’s**

 ![Screen Shot 2021-08-10 at 1.05.32 AM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/416787357e6bbb9511bad308b0438f6f83090f75.png)

These Pumas results appear slightly lower on average than the NONMEM results, but within the margin of error. Thank you all for all your help!

---

<div class="post-metadata">

**Author:** ![benjaminrich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/benjaminrich/32/156_2.png) [@benjaminrich](https://discourse.pumas.ai/u/benjaminrich)\
**Post date:** [August 10, 2021, 12:41pm UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/10 "2021-08-10T12:41:13Z")

</div>

Hi @dvaish. This is not resolved in my opinion. The CP curves should completely overlap (the “margin of error” should be much smaller than what you are seeing). So I would continue to look for an explanation for the difference you are observing.

---

<div class="post-metadata">

**Author:** ![benjaminrich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/benjaminrich/32/156_2.png) [@benjaminrich](https://discourse.pumas.ai/u/benjaminrich)\
**Post date:** [August 10, 2021, 12:43pm UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/11 "2021-08-10T12:43:40Z")

</div>

For the TIME vs. CP figure you generated, what AGE and WT did you use?

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 11, 2021, 6:53am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/12 "2021-08-11T06:53:36Z")

</div>

Hi @benjaminrich  
WT = 8, AGE = 5

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 11, 2021, 7:19am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/14 "2021-08-11T07:19:45Z")

</div>

Actually you are absolutely correct @benjaminrich, it appears that in rerunning the model I misspecified one of the parameters, and in fact, upon rerunning it just now, I have achieved exactly the correct answer.

 ![Screen Shot 2021-08-11 at 12.18.39 AM](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/47320ce67bfd50b88363c90caa4dea53d841315f.png)

The lines are completely indistinguishable! Thank you so much, this was incredibly helpful I look forward to using this tool much more in the future!

---

<div class="post-metadata">

**Author:** ![dvaish](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/dvaish/32/151_2.png) [@dvaish](https://discourse.pumas.ai/u/dvaish)\
**Post date:** [August 11, 2021, 7:27am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/15 "2021-08-11T07:27:15Z")

</div>

As a followup @benjaminrich or @pkofod, part of my modeling workflow is to use the NSUB=100 feature in NONMEM simulation, which allows for 100 separate simulations of the data (where I got the AUC distributions). Can this be replicated by some option in Pumas or should I just be bootstrapping and rerunning the simulation 100 times ?

---

<div class="post-metadata">

**Author:** ![benjaminrich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/benjaminrich/32/156_2.png) [@benjaminrich](https://discourse.pumas.ai/u/benjaminrich)\
**Post date:** [August 11, 2021, 11:09am UTC](https://discourse.pumas.ai/t/comparing-to-nonmem-bioav-lags/494/16 "2021-08-11T11:09:55Z")

</div>

That is excellent, just what I was hoping and expecting.

For the simulation replicates, you can try something like this:

```auto
nsub = 100
sim_replicates = [DataFrame(simobs(model, population, fixeffs, obstimes = 0:4:24)) for i in 1:nsub]
sim_replicates = vcat(sim_replicates...)

```