Fitting a buprenorphine 3 compartment model with FOCEI

Estimation
1

Hello,
I am trying to estimate a 3 compartment model using FOCEI, however Julia just seems to run for hours and doesn’t seem to give any results (Usually get several warnings of “larger maxiters needed”). I’ve also tried to do this with nlmixr, and it does the same. I’ve gotten mildly good first initial estimates (attached at the bottom), I’m wondering if there’s something wrong with the model or if there is a better approach to shorten the estimation time. Any advice would be appreciated!

using Pumas, DataFrames, Plots, CSV

# 1.) Build the 3 cmt buprenorphine Model
model = @model begin
  @param begin                                      # good initial estimates:
    tfs ∈ RealDomain(lower=0.0,init=.3, upper=1.0)  # unitless  fs = 0.3
    tka ∈ RealDomain(lower=0.0,init=.8)             # /hour     ka = 0.8
    tcl ∈ RealDomain(lower=0.0,init=200)            # L/hour    cl = 200
    tv1 ∈ RealDomain(lower=0.0,init=300)            # L         v1 = 300
    tv2 ∈ RealDomain(lower=0.0,init=100)            # L         v2 = 100
    tv3 ∈ RealDomain(lower=0.0,init=400)            # L         v3 = 400
    tq2 ∈ RealDomain(lower=0.0,init=100)            # L/hour    q2 = 100
    tq3 ∈ RealDomain(lower=0.0,init=50)             # L/hour    q3 = 50
  # Error
    Ω ∈ PDiagDomain(init=[0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])
    σ_prop ∈ RealDomain(init=0.01)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    bioav = tfs*exp(η[1])
    Ka = tka*exp(η[2])
    CL = tcl*exp(η[3])
    V1 = tv1*exp(η[4])
    V2 = tv2*exp(η[5])
    V3 = tv3*exp(η[6])
    Q2 = tq2*exp(η[7])
    Q3 = tq3*exp(η[8])
  end

  @dynamics begin
    Depot' = -Ka*Depot
    Central' = Ka*Depot - (CL/V1)*Central - (Q2/V1)*Central + (Q2/V2)*Peripheral1 - (Q3/V1)*Central + (Q3/V3)*Peripheral2
    Peripheral1' = (Q2/V1)*Central - (Q2/V2)*Peripheral1
    Peripheral2' = (Q3/V1)*Central - (Q3/V3)*Peripheral2
  end

  @derived begin
    conc = @. Central/V1   # ng/L
    dv ~ @. Normal(conc, 0.05)
  end

end

# 2.) Simulate for initial parameters
sdata = CSV.read("Subutex Clinical Pharmacology Review Data Pumas.csv")
sdose = DosageRegimen(8000000, time=0) # 8 mg to nanograms
spop = Population(map(i-> Subject(id=i, evs=sdose),1:24)) # arbitrary # of pts
sparam = init_param(model)
sim = simobs(model, spop, sparam, obstimes=0:1:50)
plot(sim);
plot!(sdata[:time], sdata[:dv], seriestype=:scatter)
plot!(title="Subutex 8 mg single dose Simulation")
plot!(xlabel="Time (hr)",ylabel="Concentration (ng/mL)")

# 3.) Get the fitted parameters for Subutex
ssubject = read_pumas(sdata, dvs=[:dv]) # dv is in ng/L
sfit=fit(model, ssubject, sparam, Pumas.FOCEI())
fittedsim = simobs(model, ssubject, sfit.param, obstimes=0:1:50)# S_res.param
plot(fittedsim)
plot!(sdata[:time],sdata[:dv],seriestype=:scatter)

3cmt initial estimates
3cmt initial estimates728×244 34.5 KB

2

To debug this, it would be helpful if you could try to fit a version of the model without a bioavailability parameter. Please also add optimize_fn=Pumas.DefaultOptimizeFN(show_trace=true) to fit command to see the progress of the optimization.

3

Thanks for the fast reply. I commented out the bioavailability, but i still accounted for it in the dynamics block by doing Central’ = .3KaDepot - …

Looking at the optimize_fn function,the function value started off at 1.812012e+09, so it has made a small amount of progress. I’m not fully sure how to interpret g(x), ~inv(H), and x.

At step 108, im getting a “Warning: First function call produced NaNs. Exiting.”

108     1.812009e+09     2.570311e+01
 * Current step size: 0.2890559494397361
 * time: 4688.763999938965
 * g(x): [10.822483569868893, -5.146782670283571, 25.70311488194548, -0.213140858742625, 0.30042403405228413, -0.1365250052313456, -3.920713837784877, -0.02723066681084063, 0.00547909458785506, -0.19483871013087392, -1.252096516904058, -0.006029698849097986, 0.9323236364607219, 0.0851854160972148, 0.0]
 * ~inv(H): [1.702410943500223e-5 -0.0005358811882408853 1.7924984971802362e-5 -0.01043850392239277 0.00033923719858368497 -0.009748594354678342 0.00019371269940528712 -0.007183006671199019 0.00028403968623342456 -0.01388304563803872 0.05408608758241748 -0.0007093134433752714 0.05263859092406791 -0.00564656851916979 0.0; -0.0005358811882408853 0.022423060210639643 -0.0007414585344283164 0.43087105179035745 -0.014027740407651987 0.401664581124126 -0.007978948470543718 0.29637715954294 -0.011662135187656797 0.5722904817732553 -2.230768247790538 0.029129242432340247 -2.1710643873346975 0.23272093414977824 0.0; 1.7924984971802796e-5 -0.0007414585344283164 2.7818997980747202e-5 -0.014413739825238416 0.0004637352261464159 -0.013432443382045127 0.0002677256935980234 -0.009950022392568439 0.00038587759856026756 -0.01913755574428883 0.07476437701577626 -0.0009813047917807127 0.07276143291534191 -0.007788372645593027 0.0; -0.01043850392239277 0.43087105179035745 -0.014413739825238416 8.349458831291223 -0.26942960365680224 7.762255093534141 -0.15369789788040109 5.71692823224385 -0.22815411018750456 11.088536586144512 -43.114986229365854 0.5726069155592541 -41.9621755372456 4.5106849958141435 0.0; 0.0003392371985836848 -0.014027740407651987 0.0004637352261464162 -0.26942960365680224 0.009641620275980005 -0.25153058743284595 0.004974782277003171 -0.18558282210771046 0.007346066149057995 -0.3592856501621107 1.3988465390028468 -0.018929643254172172 1.3613943415503633 -0.14643583914060299 0.0; -0.009748594354678342 0.40166458112412606 -0.013432443382045127 7.762255093534141 -0.25153058743284595 7.2486418225862534 -0.14318364307513512 5.3273593798525845 -0.21245376383169798 10.330615755227484 -40.173881852297235 0.5320668403779819 -39.09976253220818 4.204938392333567 0.0; 0.0001937126994052871 -0.007978948470543716 0.00026772569359802685 -0.15369789788040109 0.004974782277003169 -0.14318364307513512 0.0029208184067513313 -0.10518688126976904 0.004235134120574836 -0.20461656179727697 0.7942663073128018 -0.010643339639102217 0.7730477573777738 -0.08330143039304279 0.0; -0.007183006671199019 0.29637715954294 -0.009950022392568437 5.716928232243849 -0.18558282210771043 5.327359379852584 -0.10518688126976902 3.9901146790186854 -0.14826262439203652 7.547380866464317 -29.595996494753003 0.36763995186051845 -28.80420297640562 3.0656943285684677 0.0; 0.0002840396862334245 -0.011662135187656797 0.00038587759856026756 -0.2281541101875045 0.007346066149057993 -0.21245376383169798 0.0042351341205748355 -0.14826262439203652 0.03632682009029203 -0.31700242105349674 1.1778059994352597 -0.020913288405818384 1.1463109967167762 -0.1291084722415903 0.0; -0.01388304563803872 0.5722904817732553 -0.01913755574428883 11.088536586144512 -0.3592856501621107 10.330615755227484 -0.20461656179727697 7.547380866464317 -0.31700242105349674 14.886448340499516 -57.35553717564917 0.8001974156518219 -55.820871216065676 6.039498861534043 0.0; 0.05408608758241748 
-2.230768247790538 0.07476437701577626 -43.11498622936585 1.3988465390028468 -40.17388185229723 0.7942663073128018 -29.595996494753003 1.1778059994352599 -57.35553717564917 223.10903020064987 -2.9480377453542537 217.14108851755543 -23.304438967232237 0.0; -0.0007093134433752713 0.029129242432340247 -0.0009813047917807125 0.5726069155592541 -0.018929643254172172 0.5320668403779819 -0.010643339639102216 0.36763995186051845 -0.020913288405818384 0.8001974156518219 -2.9480377453542537 0.08227944955091257 -2.8687329299327526 0.3259943672468034 0.0; 0.0526385909240679 -2.1710643873346975 0.07276143291534193 -41.96217553724559 1.3613943415503633 -39.09976253220818 0.7730477573777736 -28.804202976405616 1.1463109967167762 -55.82087121606567 217.14108851755543 -2.8687329299327526 211.33298628501475 -22.680948413465018 0.0; -0.005646568519169789 0.23272093414977824 -0.007788372645593027 4.510684995814142 -0.14643583914060296 4.204938392333567 -0.08330143039304279 3.0656943285684677 -0.12910847224159028 6.039498861534042 -23.304438967232237 0.3259943672468034 -22.680948413465018 2.4771907353036613 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.028295160737636375]
 * x: [-0.12939770546327528, 4.505033697127493, 5.8749869335810985, 4.452057705278781, 8.002407036963088, 1.748516778903211, 5.295639220599646, -10.363896690290215, -5.9481110744428385, -10.756214198628768, -1.6931710247152316, 
-4.785547875190446, -1.9898913688283661, -7.250191457639581, 0.01]
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\initdt.jl:135
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\solve.jl:455
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\DiffEqBase\k3AhB\src\integrator_interface.jl:323
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\initdt.jl:135
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\solve.jl:455
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\DiffEqBase\k3AhB\src\integrator_interface.jl:323
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\initdt.jl:135
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\OrdinaryDiffEq\kFkXd\src\solve.jl:455
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\asybi\.juliapro\JuliaPro_v1.4.1-1\packages\DiffEqBase\k3AhB\src\integrator_interface.jl:323
4

The objective function value is pretty large which sometimes indicates a problem with the data. Could you try to run Pumas.findinfluential(model, ssubject, sparam, Pumas.FOCEI(), k=10)? It will give you a list of the subjects that are contributing most to the objective function. Sometimes you’ll see that one or two completely dominates. Please post the output here, if possible.

I also noticed that you include σ_prop in the estimation but that the parameter isn’t used in the model. Notice the zero gradient element in g(x). I don’t think this is the reason why you are having problem fitting this model but it might be a good idea to remove the parameter if you aren’t using it.

5

I used k=8. k=9 & 10 was giving me a bounds error. I’ve also reattached the model I used to account for the changes.
8-element Array{Tuple{String,Float64},1}:
(“2”, 4.6330697500688064e8)
(“4”, 4.569537416167642e8)
(“3”, 4.5610205518180084e8)
(“6”, 4.559816859018432e8)
(“5”, 4.5529599238322175e8)
(“1”, 4.5333662761216176e8)
(“7”, 4.4644052300875634e8)
(“8”, 4.3660682253425854e8)

model = @model begin
  @param begin                                      # good initial estimates:
    #tfs ∈ RealDomain(lower=0.0,init=.3, upper=1.0)  # unitless  fs = 0.3
    tka ∈ RealDomain(lower=0.0,init=.8)             # /hour     ka = 0.8
    tcl ∈ RealDomain(lower=0.0,init=200)            # L/hour    cl = 200
    tv1 ∈ RealDomain(lower=0.0,init=300)            # L         v1 = 300
    tv2 ∈ RealDomain(lower=0.0,init=100)            # L         v2 = 100
    tv3 ∈ RealDomain(lower=0.0,init=400)            # L         v3 = 400
    tq2 ∈ RealDomain(lower=0.0,init=100)            # L/hour    q2 = 100
    tq3 ∈ RealDomain(lower=0.0,init=50)             # L/hour    q3 = 50
  # Error
    Ω ∈ PDiagDomain(init=[0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01])
    #σ_prop ∈ RealDomain(init=0.01)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    #bioav = tfs*exp(η[1])
    Ka = tka*exp(η[1])
    CL = tcl*exp(η[2])
    V1 = tv1*exp(η[3])
    V2 = tv2*exp(η[4])
    V3 = tv3*exp(η[5])
    Q2 = tq2*exp(η[6])
    Q3 = tq3*exp(η[7])
  end

  @dynamics begin
    Depot' = -Ka*Depot
    Central' = .3*Ka*Depot - (CL/V1)*Central - (Q2/V1)*Central + (Q2/V2)*Peripheral1 - (Q3/V1)*Central + (Q3/V3)*Peripheral2
    Peripheral1' = (Q2/V1)*Central - (Q2/V2)*Peripheral1
    Peripheral2' = (Q3/V1)*Central - (Q3/V3)*Peripheral2
  end

  @derived begin
    conc = @. Central/V1   # ng/L
    dv ~ @. Normal(conc, 0.05)
  end

end
6

Without seeing the data it’s a difficult to know for sure but I think the small standard deviation could cause some numerical issues. I.e. when you are measuring the concentration in ng then a standard deviation of 0.05 is very small. I’d consider bumping it to say 100 and also maybe include it in the estimation. It should bring the objective function values down to more reasonable magnitudes.

1 Like
7

That worked! Thank you so much.

image
image724×222 22 KB