Is there coding error?

Pumas issues
1

Hello,
I have a question - once the dug is absorbed into circulation, the fate should be same - same CL & Volume of distribution despite the dose route, normal PK or flip-flop PK. Is this true? If IM displays flip-flop PK and PO displays normal/conventional PK, is it good idea to simultaneously fit PO/IM PK to avoid switch of kel and ka? Here is the code I used to simulate IM PK profile using fitted parameters. Depot is for PO route here (Set oral bioavailability Fpo=1, solve FIM for fitting). IM PK (IR and SR) is described as parallel absorption - slow process with transit compartment delay. The code works for simulation but the simulation output doesn’t look right (pkdata_im = DataFrame(obs), some parameters show NA in the output csv file) . Could you help check the code? Many Thanks!

two_parallel_foabs = @model begin
  @param begin
    tvcl  ∈ RealDomain(lower=0)
    tvvc  ∈ RealDomain(lower=0)
    tvka ∈ RealDomain(lower=0)
    tvka1 ∈ RealDomain(lower=0)
    tvka2 ∈ RealDomain(lower=0)
    tvlag ∈ RealDomain(lower=0)
    tvFracIM ∈ RealDomain(lower=0)
    tvvp ∈ RealDomain(lower=0)
    tvq  ∈ RealDomain(lower=0)
    tvvp2 ∈ RealDomain(lower=0)
    tvq2  ∈ RealDomain(lower=0)
    tvktr ∈ RealDomain(lower=0)
    tvfIM ∈ RealDomain(lower=0)
    Ω    ∈ PDiagDomain(8)
    σ_prop    ∈ RealDomain(lower=0)
  end

  @random begin
    η ~ MvNormal(Ω)
  end


  @pre begin
    CL  = tvcl*exp(η[1])
    Vc  = tvvc
    Vp = tvvp 
    Q  = tvq 
    Vp2 = tvvp2 
    Q2  = tvq2*exp(η[2])
    Ka =  tvka*exp(η[3])
    Ka1 = tvka1*exp(η[4])
    Ka2 = tvka2*exp(η[5])
    Ktr = tvktr*exp(η[7])
  end

  @dosecontrol begin
    lags    = (Depot = tvlag,)
    bioav = ( Depot=1, IR = tvfIM*exp(η[8])*(1-tvFracIM)*exp(η[6]), SR = tvfIM*exp(η[8])*tvFracIM*exp(η[6]))
  end

  @dynamics begin
    Depot'   = -Ka*Depot
    IR'      = -Ka1*IR
    SR'    = -Ktr*SR
    Transit1' = Ktr*SR   - Ktr*Transit1
    Transit2' = Ktr*Transit1 - Ktr*Transit2
    Transit3' = Ktr*Transit2 - Ktr*Transit3
    Transit4' = Ktr*Transit3 - Ktr*Transit4
    Transit5' = Ktr*Transit4 - Ktr*Transit5
    Transit6' = Ktr*Transit5 - Ktr*Transit6
    Transit7' = Ktr*Transit6 - Ka2*Transit7
    Central' = Ka*Depot +  Ka1*IR + Ka2*Transit7 - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral1  - (Q2/Vc)*Central + (Q2/Vp2)*Peripheral2
    Peripheral1' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral1
    Peripheral2' =  (Q2/Vc)*Central  - (Q2/Vp2)*Peripheral2
  end

  @derived begin
    cp = @. (Central/Vc)
    dv ~ @. Normal(abs(cp), abs(cp)*σ_prop)
end

end


param = ( tvcl = 14, tvvc = 43, tvq = 1,
          tvvp = 8.7,
          tvq2 = 0.7,
          tvvp2 = 206,
          tvka = 1.28, tvka1 = 0.0024, tvka2 = 0.001, tvlag = 0.18, tvFracIM = 0.49, tvktr=0.0166, tvfIM=3, 
           Ω = Diagonal([0, 0, 0, 0, 0, 0,0,0]),σ_prop =0.00001)

## simulate a profile with this type of absorption:

dose_fo1 = DosageRegimen(112500, cmt = 2, time = 0)
dose_fo2 = DosageRegimen(112500, cmt = 3, time = 0)
dose     = DosageRegimen(dose_fo1, dose_fo2)  # Actual dose is made up of 2 virtual doses
s3_IM = Subject(id=3, events=dose,observations = (conc = nothing,))

subj_with_covariates = map(1:10) do i
    Subject(id = i,
            events = dose,
            covariates = choose_covariates(),
            observations = (conc = nothing,))
end

obs = simobs(two_parallel_foabs, s3_IM, param, obstimes = 0:1:4800)

## plot the results

sim_plot(two_parallel_foabs, obs, observations = :cp)

pkdata_im = DataFrame(obs)

CSV.write("../data/pk_im.csv", pkdata_im)
2

Hi Jenny,

The code looks fine to me, but overall since you have 0 on your etas you don’t end up getting values on the dataframe. If you slightly fix the etas to a very low value you should see them.

param = ( tvcl = 14, tvvc = 43, tvq = 1,
          tvvp = 8.7,
          tvq2 = 0.7,
          tvvp2 = 206,
          tvka = 1.28, tvka1 = 0.0024, tvka2 = 0.001, tvlag = 0.18, tvFracIM = 0.49, tvktr=0.0166, tvfIM=3, 
           Ω = Diagonal([0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001,0.0001,0.0001]),σ_prop =0.00001)
3

Yes, exactly. I made the same changes and it worked. Thanks so much!

4

I’m not sure that is true. In flip-flop kinetic, we could end up with different parameter estimate if we keep changing the starting point so you could end up with different set of solution (different PK parameters) each time so your CL and V could be different in terms of estimation process (flop vs flip-flop). Things might be the same in theoretical point of view in the elimination. If you want to do simultaneous fit, I encourage you to think about the way you parameterize your IM parameters that will avoid identifiability and end up with different solution. Take a look at this fine commentary: The influence of flip‐flop in population pharmacokinetic analyses - PMC

5

Thanks so much for sharing your thoughts and paper! Much appreciated!