While trying to recreate the nonmem code, the following error occured in the below given code.
ERROR: gradient is exactly zero in tvftrt. This indicates that tvftrt isn’t identified.
any help to resolve is highly appreciated.
thanks.
pop1 = read_pumas(PK_data1,
id = :subject,
time = :act_time,
observations = [:conc],
amt = :dose,
evid = :evid,
covariates = [:wt, :trt, :per,:seq],
route = :route,
cmt = :cmt)
model1 = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvf ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
tvftrt ∈ RealDomain(; lower=0)
tvkatrt ∈ RealDomain(; lower = 0)
tvfseq ∈ RealDomain(; lower=0)
tvkaseq ∈ RealDomain(; lower = 0)
tvfper ∈ RealDomain(; lower=0)
tvkaper ∈ RealDomain(; lower = 0)
Ω_BSV ∈ PSDDomain(3)
Ω_IOV ∈ PSDDomain(2)
σ_prop ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω_BSV)
K ~ MvNormal(Ω_IOV)
end
@covariates wt per trt
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
IOVF = per == "1" ? exp(K[1]) : exp(K[2])
Katrt = trt == "R" ? 1.0 : tvkatrt
Ftrt = trt == "R" ? 1.0 : tvftrt
Kaseq = trt == "R" ? 1.0 : tvkaseq
Fseq = trt == "R" ? 1.0 : tvfseq
Kaper = trt == "R" ? 1.0 : tvkaper
Fper = trt == "R" ? 1.0 : tvfper
F = tvf * exp(IOVF)* Ftrt *Fseq * Fper
Ka = tvka * exp(η[3])* Katrt *Kaseq * Kaper
end
@dynamics Depots1Central1
@derived begin
cp = @. 1000* @. Central/Vc
conc ~ @. Normal(cp, sqrt((cp*σ_prop)^2))
end
end
param1 =
(; tvcl = 6.83, tvvc = 32.4, tvka = 1.2, tvf = 1, tvkatrt = 1.1, tvftrt = 1.1, tvkaseq = 1.1, tvfseq = 1.1, tvkaper = 1.1, tvfper = 1.1,
Ω_BSV = Diagonal([0.03, 0.02, 0.03]), Ω_IOV = Diagonal([0.025, 0.025]),
σ_prop = 0.03)
fit_res = fit(model1,
pop1,
param1,
Pumas.FOCEI(), constantcoef=(tvf=1,))
Hi @yashokrajv,
I assume you are trying to replicate a relative bioavailability model. There cannot be three different bioavailability for the treatment, sequence and period.
Here is an example for that.
In the time vs. concentration plot we can see that test formulation has 50% of exposure compared to reference.
model = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
tvf_ref ∈ RealDomain(; lower = 0)
tvf_test ∈ RealDomain(; lower = 0)
Ω_BSV ∈ PSDDomain(4)
Ω_IOV ∈ RealDomain(lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω_BSV)
κ ~ MvNormal(Diagonal(fill(Ω_IOV, 2)))
end
@covariates begin
Formulation
Period
end
@pre begin
Ka = Formulation == "R" ? tvka * exp(η[1]) : tvka * exp(η[2])
CL = tvcl * exp(η[3])
Vc = tvvc * exp(η[4])
end
@dosecontrol begin
bioav = (; Depot=Formulation == "R" ? tvf_ref * exp(κ[Period]) : tvf_test * exp(κ[Period]))
end
@dynamics Depots1Central1
@derived begin
cp = @. 1000* @. Central/Vc
conc ~ @. Normal(cp, sqrt((cp*σ_prop)^2))
end
end
param =
(; tvcl = 6.83,
tvvc = 32.4,
tvka = 1.2,
tvf_ref = 1,
tvf_test = 0.5,
Ω_BSV = Diagonal([0.04, 0.04, 0.04, 0.04]),
Ω_IOV = 0.04,
σ_prop = 0.03)
sub1 = Subject(; id="1", events=DosageRegimen(5, time=0), covariates=(; Formulation="R", Period=1))
sub2 = Subject(; id="2", events=DosageRegimen(5, time=0), covariates=(; Formulation="T", Period=1))
Pop = [sub1, sub2]
Random.seed!(123)
sim = simobs(model, Pop, param, obstimes=0:1:24)
data(DataFrame(sim))*
mapping(:time => "Time", :cp => "Concentration", color=:Formulation)*
visual(ScatterLines)|>draw
Hi @Jayashree
Thanks for your effort to work on this query and respond to it! We appreciate it!!
I am giving the nonmem code here which i was trying to recreate in PUMAS. It may provide you more understanding on the question of need of additional F and Ka parameters for treatment, sequence and period. Link of the source article is also given here.
Chen X, Nyberg HB, Donnelly M, et al. Development and comparison of
model-integrated evidence approaches for bioequivalence studies with pharmacokinetic end points. CPT Pharmacometrics Syst Pharmacol.
2024;13:1734-1747. doi:10.1002/psp4.13216
Thanks again.
Supplementary Material 1
Example NONMEM code used in the model-integrated bioequivalence method
- Example NONMEM code for estimation step
$PROBLEM PK estimation
$INPUT REP=DROP ID TIME AMT EVID DV TRT SEQ PER
$DATA be_data.csv IGNORE=@
$SUBROUTINE ADVAN2 TRANS2
$PK
FTRT = 1
IF(TRT.EQ.2) FTRT = THETA(5)
KATRT = 1
IF(TRT.EQ.2) KATRT = THETA(6)
FSEQ = 1
IF(SEQ.EQ.2) FSEQ = THETA(7)
KASEQ = 1
IF(SEQ.EQ.2) KASEQ = THETA(8)
FPER = 1
IF(PER.EQ.2) FPER = THETA(9)
KAPER = 1
IF(PER.EQ.2) KAPER = THETA(10)
IF (PER.EQ.1) IOVF = ETA(4)
IF (PER.EQ.2) IOVF = ETA(5)
CL = THETA(1)* EXP(ETA(1))
V = THETA(2)*EXP(ETA(2))
F1 = THETA(3) * EXP(IOVF) * FTRT * FSEQ * FPER
KA = THETA(4) * EXP(ETA(3)) * KATRT * KASEQ * KAPER
$ERROR
IPRED = A(2)/V
IRES = DV-IPRED
PROP = SQRT(SIGMA(1,1))IPRED
W = PROP
IF (W.EQ.0) W=1
IWRES = IRES/W
Y = IPRED + IPREDEPS(1)
$THETA
(0,46) ; 1.CL
(0,400) ; 2.V
1 FIX. ; 3.F
(0, 1.2) ; 4.KA
(0,1.25) ; 5.FTRT
(0, 1.1) ; 6. KATRT
(0, 1.1) ; 7.FSEQ
(0, 1.1) ; 8. KASEQ
(0, 1.1) ; 9.FPER
(0, 1.1) ; 10. KAPER
$OMEGA
0.03. ; 1. CL
0.02 ; 2. V
0.03 ; 3. KA
$OMEGA BLOCK(1) 0.0025 ;IOVF
$OMEGA BLOCK(1) SAME ;IOVF
$SIGMA
0.01 ;PROP
;;; estimation step
$ESTIMATION PRINT=1 MAXEVAL=9999 METHOD=1 INTER NOABORT SADDLE_RESET=1
;;; uncertainty: covariance matrix
$COVARIANCE UNCONDITIONAL
;;; uncertainty: SIR method
$COV SIRSAMPLE=2000 SIRNITER=6 SIRDF=30 IACCEPT=0.4 FILE=output_sir.ext
;;; uncertainty: bootstrap method
$SIML (101) BOOTSTRAP=-1 SUBP=1000 STRAT=SEQ
