Turnover PK model simulation

Dawoon · March 13, 2022, 8:26pm

Hello,
I am simulating a turnover PK model using an example shown in tutorial notebooks in pumas.ai. The model I am trying is Exercise PK49. I modified PK parameters to replicate my data. If I use 1 cmt model as shown in the example, the PK curve is plotted as I expected. However, if I modified the model to 2cmt model, the PK curves become weird. The only difference between IDs is dose but ID1 has lower concentrations than baseline after a peak and rebounds to the baseline. (I think ID2 looks good) The baseline value is about 90. Also, I don’t see distribution and elimination phase in the semi-log scale. Does anyone see something wrong in my code? I think the issue might have come from @initial block.
Here’s my simulation codes.

data = DataFrame(id = [1,2],
                                      time = [0,0],
                                     amt = [160, 320],
                                     evid = [1,1],
                                     cmt = [1,1],
                                     rate = [640,1280],
                                     weight_kg = [3.2,3.2])

data[!, :at3_plasma] .= missing

pkdata = read_pumas(data,
                          observations  = [:at3_plasma,],
                          id            = :id,
                          time          = :time,
                          amt           = :amt,
                          rate          = :rate,
                          evid          = :evid,
                          cmt           = :cmt,
                          covariates    = [:weight_kg])

rhAT3_model = @model begin
    @param begin
        tvcl              ∈ RealDomain(lower=0.001)
        tvvc             ∈ RealDomain(lower=0.001)
        tvvp             ∈ RealDomain(lower=0.001)
        tvq               ∈ RealDomain(lower=0.001)
        tvkin           ∈ RealDomain(lower=0.001)
        Ω                 ∈ PDiagDomain(2)
        σ_prop      ∈ RealDomain(lower=0.001)
    end

    @random begin
        η   ~   MvNormal(Ω)
    end

    @covariates weight_kg

    @pre begin
      CL        = tvcl * weight_kg
      Vc        = tvvc * weight_kg
      Vp        = tvvp * weight_kg *exp(η[1])
      Q         = tvq * weight_kg * exp(η[2])
      Kin       = tvkin * weight_kg      
    end

    @init begin
      Central   =  Kin/(CL/Vc)
    end

    @dynamics begin
      # 2cmt
      Central'    =  Kin - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
      Peripheral' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral
    end

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


rhAT3_param =   (tvcl        =  0.0383*60/100,
                                   tvvc        =  41.1/100,
                                   tvvp        =  74.3/100,
                                   tvq         =  0.0763*60/100,
                                   tvkin       =  2.1,
                                   Ω           =  Diagonal([0.1785, 0.0275]),
                                   σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(rhAT3_model, rhAT3_param, pkdata)
sims_zero = simobs(rhAT3_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:0.25:96)

e, f, p = sim_plot(rhAT3_model, [sims_zero][1], 
        observations = :CP, 
        color = :redsblues,
        linewidth = 4,
        axis = (xlabel = "Time (hrs)", 
                ylabel = "Conc",
                xticks = 0:24:96))
axislegend(f) 
e

Dawoon · March 13, 2022, 8:27pm

FYI, I am attaching the result of 1cmt model simulation. 1cmt looks good.

Rahulub3r · March 13, 2022, 10:30pm

Hi Dawoon - If you look from a mass balance perspective, there is a constant rate of input of endogenous AT3 in the form of Kin (in your case if will be 2.1 x 3.2 = 6.72 mg/hr), the rate of output is the second part of your differential equation for dCentral/dt (- (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral). For ID=1, the rate of input is higher than the rate of output at one point of time because of which you see that increase in concentrations.

To get more intuition for this, I ran an example code (attached below). You can see that in the output of the last line of code, the column dCentral_dt turns positive at time 22.5 hrs which means there is a net positive rate in the central compartment.

Let me know if something is not clear.

Code

using Pumas
using CairoMakie
using PumasUtilities

data = DataFrame(id = [1,2,3],
                time = [0,0,0],
                amt = [160,320,320],
                evid = [1,1,1],
                cmt = [1,1,1],
                rate = [640,1280,2560],
                weight_kg = [3.2,3.2,3.2])

data[!, :at3_plasma] .= missing

pkdata = read_pumas(data,
                    observations  = [:at3_plasma,],
                    id            = :id,
                    time          = :time,
                    amt           = :amt,
                    rate          = :rate,
                    evid          = :evid,
                    cmt           = :cmt,
                    covariates    = [:weight_kg])

rhAT3_model = @model begin
    @param begin
        tvcl        ∈ RealDomain(lower=0.001)
        tvvc        ∈ RealDomain(lower=0.001)
        tvvp        ∈ RealDomain(lower=0.001)
        tvq         ∈ RealDomain(lower=0.001)
        tvkin       ∈ RealDomain(lower=0.001)
        Ω           ∈ PDiagDomain(2)
        σ_prop      ∈ RealDomain(lower=0.001)
    end

    @random begin
        η   ~   MvNormal(Ω)
    end

    @covariates weight_kg

    @pre begin
      CL        = tvcl  * weight_kg
      Vc        = tvvc  * weight_kg
      Vp        = tvvp  * weight_kg * exp(η[1])
      Q         = tvq   * weight_kg * exp(η[2])
      Kin       = tvkin * weight_kg      
    end

    @init begin
      Central   =  Kin/(CL/Vc)
    end

    @dynamics begin
      # 2cmt
      Central'    =  Kin - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
      Peripheral' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral
    end

    @derived begin
      CP            =  @. Central/Vc
      at3_plasma    ~  @. Normal(CP, abs(CP) * σ_prop)
      Kin_perhr = @. Kin
      CentralOut_perhr = @. - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
    end
end


rhAT3_param =   (tvcl       =  0.0383*60/100, #0.02298
                tvvc        =  41.1/100, #0.411
                tvvp        =  74.3/100,
                tvq         =  0.0763*60/100,
                tvkin       =  2.1,
                Ω           =  Diagonal([0.1785, 0.0275]),
                σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(rhAT3_model, rhAT3_param, pkdata)
sims_zero = simobs(rhAT3_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:0.25:48)

e, f, p = sim_plot(rhAT3_model, [sims_zero][1], 
                    observations = :CP, 
                    color = :redsblues,
                    linewidth = 4,
                    axis = (xlabel = "Time (hrs)", 
                            ylabel = "Conc", yscale=log10,
                            xticks = 0:24:96))
axislegend(f) 
e

sims_zero_df = DataFrame(sims_zero)

# Find the minimum row for ID = 1
findmin(filter(x -> x.evid == 0 && x.id == "1", sims_zero_df).CP)

# Checkout rows
sims_zero_df[:, :dCentral_dt] .= sims_zero_df[:, :Kin_perhr] .+ sims_zero_df[:, :CentralOut_perhr]
filter(x -> x.id == "1" && x.time > 21 && x.time < 24, sims_zero_df) # min at 22.25

Dawoon · March 14, 2022, 12:05am

Hello Rahul,
Thank you so much for your explanation and codes. Your codes provided me a better idea how to interpret this mathematically. And I apologize for my unclear statement. Actually my concern was that why the concentrations are below the baseline after a peak in ID 1. In other words, I don’t think this adds up physiologically because the AT3 level should be always equal to or higher than the baseline when subject was given exogenous AT3. Also, as we can see in the semi-log scale plot, I couldn’t capture distribution and elimination phase. It looks like 1 cmt model although I built 2cmt model. Do you have any ideas for this phenomenon? Thank you for your advice!

Rahulub3r · March 14, 2022, 1:35am

That’s a good point. I do not know why the values go below baseline and why the two-compartment model shape does not appear.
I tried to figure out the baseline issue using a simple model like the 1-cmp model as in the PK49 tutorial and saw a pulsating trend which I couldn’t understand. Hope we can discuss why we see such a pattern as well along with your two questions (plot and code below)

CODE

pk_49            = @model begin
  @param begin
    tvcl         ∈ RealDomain(lower=0)
    tvvc         ∈ RealDomain(lower=0)
    tvsynthesis  ∈ RealDomain(lower=0)
    Ω            ∈ PDiagDomain(3)
  end
  @random η ~ MvNormal(Ω)
  @pre begin
    CL           = tvcl * exp(η[1])
    Vc           = tvvc * exp(η[2])
    Synthesis    = tvsynthesis * exp(η[3])
  end
  @init Central = Synthesis/(CL/Vc)
  @dynamics Central' = Synthesis - (CL/Vc)*Central
  @derived cp = @. Central/Vc
end

param = ( tvcl        = 0.14204,
          tvvc        = 3.59259,
          tvsynthesis = 16.2272,
          Ω           = Diagonal([0.0,0.0,0.0]))

ev1  = DosageRegimen(1, cmt=1) #0.3166
sub1 = Subject(id=1, events=ev1)

sim_sub1 = simobs(pk_49, sub1, param, obstimes=0:1:2055)

f, a, p = sim_plot(pk_49, [sim_sub1], 
        observations = :cp, 
        color = :redsblues,
        linewidth = 4,
        axis = (xlabel = "Time (hrs)", 
                ylabel = "PK49 Concentrations (μg/L)",
                #yscale = log10,
                #xticks = 0:20:160, 
                ))
axislegend(a) 
f

Dawoon · March 14, 2022, 2:16am

Thank you, Rahul.
It’s very interesting. I can’t comp up with any reasons which cause the fluctuation in terms of modeling. I think the fluctuation is negligible physiologically because the magnitude of fluctuation is about 0.2 ug/L and the baseline is about 114 ug/L. I noticed that a given dose in your codes was 1 and we may not be able to capture the fluctuation if we give a larger dose than 1 (like 100) and adjust y axis to a larger scale. But, I agree with you that the fluctuation per se seems not to make sense in terms of modeling. Like Rahul said, any opinions to our discussion would be appreciated!

shli6161 · March 14, 2022, 2:42am

Interesting discussion! I tried a different approach from @Rahulub3r i tested three scenarios. All of them are without giving At3. Just to see how the turnover model looks like. So the short answer is i think your two compartment code wasn’t really at steady state . Here is how i reached the conclusion.

Firstly i tested the two-compartment without giving AT3 using your code. Below are the code and graph. So essentially we see a bouncing-back all the way till 192 hours. Basically, it is the time for a two-compartment turnover model to reach ss. This is probably why you see the bouncing for ID1.

code for original two cmpt-no drug

using Pumas
using CairoMakie
using PumasUtilities

data = DataFrame(id   = 1,
                 time = 0,
                 amt  = 0,
                 evid = 1,
                 cmt  = 1,
                 rate = 640,
                 weight_kg = 3.2)

data[!, :at3_plasma] .= missing

pkdata = read_pumas(data,
                          observations  = [:at3_plasma,],
                          id            = :id,
                          time          = :time,
                          amt           = :amt,
                          rate          = :rate,
                          evid          = :evid,
                          cmt           = :cmt,
                          covariates    = [:weight_kg])
############## two cmpt #####################
twocmpt_model = @model begin
    @param begin
        tvcl              ∈ RealDomain(lower=0.001)
        tvvc             ∈ RealDomain(lower=0.001)
        tvvp             ∈ RealDomain(lower=0.001)
        tvq               ∈ RealDomain(lower=0.001)
        tvkin           ∈ RealDomain(lower=0.001)
        Ω                 ∈ PDiagDomain(2)
        σ_prop      ∈ RealDomain(lower=0.001)
    end

    @random begin
        η   ~   MvNormal(Ω)
    end

    @covariates weight_kg

    @pre begin
      CL        = tvcl * weight_kg
      Vc        = tvvc * weight_kg
      Vp        = tvvp * weight_kg *exp(η[1])
      Q         = tvq * weight_kg * exp(η[2])
      Kin       = tvkin * weight_kg      
    end

    @init begin
      Central   =  Kin/(CL/Vc)
    end

    @dynamics begin
      # 2cmt
      Central'    =  Kin - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
      Peripheral' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral
    end

    @derived begin
      CP            =  @. Central/Vc
      Cperi         =  @. Peripheral/Vp
      at3_plasma    ~  @. Normal(CP, abs(CP) * σ_prop)
    end
end


rhAT3_param =   (tvcl        =  0.0383*60/100,
                 tvvc        =  41.1/100,
                 tvvp        =  74.3/100,
                tvq         =  0.0763*60/100,
                 tvkin       =  2.1,
                 Ω           =  Diagonal([0.1785, 0.0275]),
                 σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(twocmpt_model, rhAT3_param, pkdata)
sims_zero = simobs(twocmpt_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:1:596)

e, f, p = sim_plot(twocmpt_model, [sims_zero][1], 
        observations = :CP, 
        color = :redsblues,
        linewidth = 4,
        axis = (xlabel = "Time (hrs)", 
                ylabel = "2cmpt_Conc",
                xticks = 0:96:596))
axislegend(f) 
e

Profile for original two cmpt-no drug

And i applied the same logic to the one-cmpt turnover model. We see a straight line, which is the way it supposed to be. Below are the one-cmpt code and graph.

one cmpt code

############## one cmpt #####################
#
data = DataFrame(id   = 1,
                 time = 0,
                 amt  = 0,
                 evid = 1,
                 cmt  = 1,
                 rate = 640,
                 weight_kg = 3.2)

data[!, :at3_plasma] .= missing
#
pkdata = read_pumas(data,
                          observations  = [:at3_plasma,],
                          id            = :id,
                          time          = :time,
                          amt           = :amt,
                          rate          = :rate,
                          evid          = :evid,
                          cmt           = :cmt,
                          covariates    = [:weight_kg])
#
onecmpt_model = @model begin
  @param begin
      tvcl              ∈ RealDomain(lower=0.001)
      tvvc             ∈ RealDomain(lower=0.001)
      #tvvp             ∈ RealDomain(lower=0.001)
      #tvq               ∈ RealDomain(lower=0.001)
      tvkin           ∈ RealDomain(lower=0.001)
      Ω                 ∈ PDiagDomain(2)
      σ_prop      ∈ RealDomain(lower=0.001)
  end

  @random begin
      η   ~   MvNormal(Ω)
  end

  @covariates weight_kg

  @pre begin
    CL        = tvcl * weight_kg
    Vc        = tvvc * weight_kg
    #Vp        = tvvp * weight_kg *exp(η[1])
    #Q         = tvq * weight_kg * exp(η[2])
    Kin       = tvkin * weight_kg      
  end

  @init begin
    Central   =  Kin/(CL/Vc)
  end

  @dynamics begin
    # 1cmt
    Central'    =  Kin - (CL/Vc)*Central 
  end

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


rhAT3_param = (tvcl        =  0.0383*60/100,
               tvvc        =  41.1/100,
               #tvvp        =  74.3/100,
               #tvq         =  0.0763*60/100,
               tvkin       =  2.1,
               Ω           =  Diagonal([0.1785, 0.0275]),
               σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(onecmpt_model, rhAT3_param, pkdata)
sims_zero = simobs(onecmpt_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:0.01:96)

e, f, p = sim_plot(onecmpt_model, [sims_zero][1], 
      observations = :CP, 
      color = :redsblues,
      linewidth = 4,
      axis = (xlabel = "Time (hrs)", 
              ylabel = "1cmpt_Conc",
              xticks = 0:24:96))
axislegend(f) 
e

Therefore i added a command to your original code to make sure that peripheral cmpt is also at ss in the beginning, below are the code and graph. Basically, it is to make the Peripheral’ = 0 => Peripheral = Central * Vp/Vc, then plug in Central = Kin/ (CL/Vc) => Peripheral = Kin * Vp/CL

code for modified two cmpt - no drug

data = DataFrame(id   = 1,
                 time = 0,
                 amt  = 0,
                 evid = 1,
                 cmt  = 1,
                 rate = 640,
                 weight_kg = 3.2)

data[!, :at3_plasma] .= missing

pkdata = read_pumas(data,
                          observations  = [:at3_plasma,],
                          id            = :id,
                          time          = :time,
                          amt           = :amt,
                          rate          = :rate,
                          evid          = :evid,
                          cmt           = :cmt,
                          covariates    = [:weight_kg])
############## two cmpt #####################
twocmpt_model = @model begin
    @param begin
        tvcl              ∈ RealDomain(lower=0.001)
        tvvc             ∈ RealDomain(lower=0.001)
        tvvp             ∈ RealDomain(lower=0.001)
        tvq               ∈ RealDomain(lower=0.001)
        tvkin           ∈ RealDomain(lower=0.001)
        Ω                 ∈ PDiagDomain(2)
        σ_prop      ∈ RealDomain(lower=0.001)
    end

    @random begin
        η   ~   MvNormal(Ω)
    end

    @covariates weight_kg

    @pre begin
      CL        = tvcl * weight_kg
      Vc        = tvvc * weight_kg
      Vp        = tvvp * weight_kg *exp(η[1])
      Q         = tvq * weight_kg * exp(η[2])
      Kin       = tvkin * weight_kg      
    end

    @init begin
      Central    =  Kin/(CL/Vc)
      Peripheral =  Kin * Vp/CL
    end

    @dynamics begin
      # 2cmt
      Central'    =  Kin - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
      Peripheral' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral
    end

    @derived begin
      CP            =  @. Central/Vc
      Cperi         =  @. Peripheral/Vp
      at3_plasma    ~  @. Normal(CP, abs(CP) * σ_prop)
    end
end


rhAT3_param =   (tvcl        =  0.0383*60/100,
                 tvvc        =  41.1/100,
                 tvvp        =  74.3/100,
                tvq         =  0.0763*60/100,
                 tvkin       =  2.1,
                 Ω           =  Diagonal([0.1785, 0.0275]),
                 σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(twocmpt_model, rhAT3_param, pkdata)
sims_zero = simobs(twocmpt_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:1:96)

e, f, p = sim_plot(twocmpt_model, [sims_zero][1], 
        observations = :CP, 
        color = :redsblues,
        linewidth = 4,
        axis = (xlabel = "Time (hrs)", 
                ylabel = "2cmpt_Conc",
                xticks = 0:24:96))
axislegend(f) 
e

Then we see this beautiful straight line in the profile for the two cmpt.

So i further updated the code by giving AT3. Below are the code and the profile for ID1. Doesn’t have any abnormality any more!!

modified two cmpt-with drug

data = DataFrame(id   = 1,
                 time = 0,
                 amt  = 160,
                 evid = 1,
                 cmt  = 1,
                 rate = 640,
                 weight_kg = 3.2)

data[!, :at3_plasma] .= missing

pkdata = read_pumas(data,
                          observations  = [:at3_plasma,],
                          id            = :id,
                          time          = :time,
                          amt           = :amt,
                          rate          = :rate,
                          evid          = :evid,
                          cmt           = :cmt,
                          covariates    = [:weight_kg])
############## two cmpt #####################
twocmpt_model = @model begin
    @param begin
        tvcl              ∈ RealDomain(lower=0.001)
        tvvc             ∈ RealDomain(lower=0.001)
        tvvp             ∈ RealDomain(lower=0.001)
        tvq               ∈ RealDomain(lower=0.001)
        tvkin           ∈ RealDomain(lower=0.001)
        Ω                 ∈ PDiagDomain(2)
        σ_prop      ∈ RealDomain(lower=0.001)
    end

    @random begin
        η   ~   MvNormal(Ω)
    end

    @covariates weight_kg

    @pre begin
      CL        = tvcl * weight_kg
      Vc        = tvvc * weight_kg
      Vp        = tvvp * weight_kg *exp(η[1])
      Q         = tvq * weight_kg * exp(η[2])
      Kin       = tvkin * weight_kg      
    end

    @init begin
      Central    =  Kin/(CL/Vc)
      Peripheral =  Kin * Vp/CL
    end

    @dynamics begin
      # 2cmt
      Central'    =  Kin - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
      Peripheral' =  (Q/Vc)*Central  - (Q/Vp)*Peripheral
    end

    @derived begin
      CP            =  @. Central/Vc
      Cperi         =  @. Peripheral/Vp
      at3_plasma    ~  @. Normal(CP, abs(CP) * σ_prop)
    end
end


rhAT3_param =   (tvcl        =  0.0383*60/100,
                 tvvc        =  41.1/100,
                 tvvp        =  74.3/100,
                tvq         =  0.0763*60/100,
                 tvkin       =  2.1,
                 Ω           =  Diagonal([0.1785, 0.0275]),
                 σ_prop      =  0.114)

zero_rand_eff = zero_randeffs(twocmpt_model, rhAT3_param, pkdata)
sims_zero = simobs(twocmpt_model, pkdata, rhAT3_param, zero_rand_eff, obstimes = 0:1:96)

e, f, p = sim_plot(twocmpt_model, [sims_zero][1], 
        observations = :CP, 
        color = :redsblues,
        linewidth = 4,
        axis = (xlabel = "Time (hrs)", 
                ylabel = "2cmpt_Conc",
                xticks = 0:24:96))
axislegend(f) 
e

profile for ID1

Dawoon · March 14, 2022, 3:11am

Thank you so much, Shuhui!!! I spent a lot of time but couldn’t figure out why. Thanks again for your help, Rahul and Shuhui .

shli6161 · March 14, 2022, 3:22am

No problem, glad we could help

Topic		Replies	Views
Is there coding error? Pumas issues	4	195	July 3, 2023
Comparing to NONMEM (bioav & lags) Simulation	13	538	August 11, 2021
Post-processing steps for full Bayesian model bayesian	6	204	July 19, 2022
Simultaneous PKPD Modeling with Missing PD Data for Some Individuals maximum-likelihood	3	242	April 25, 2022
IOV modeling and simulation How-to	8	192	August 29, 2023

Turnover PK model simulation

Related Topics