L2 Correlations

parssh.mehta · August 5, 2023, 7:41pm

Hi I am trying to use L2 correlation, since the pk and pd concentrations are taken from the same sample. The simulations work and it also starts the estimation but then fails with an error. Can you please help me fix this?

Attached model code and failed estimation error

true_model_pk_delay_pd_sim = @model begin
    @param begin
        tvcl      ∈ RealDomain(lower=0.0)
        tvvc      ∈ RealDomain(lower=0.0)
        tvq       ∈ RealDomain(lower=0.0)
        tvvp      ∈ RealDomain(lower=0.0)
        tvka      ∈ RealDomain(lower=0.0)
        tvkin     ∈ RealDomain(lower=0.0)
        tvkout    ∈ RealDomain(lower=0.0)
        tvImax    ∈ RealDomain(lower=0.0)
        tvIC50    ∈ RealDomain(lower=0.0)
        Ω         ∈ PDiagDomain(8)
        Σ         ∈ PSDDomain(2)
    end

    @random begin
        η      ~ MvNormal(Ω)
    end

    @pre begin
        CL     = tvcl * exp(η[1])
        Vc     = tvvc * exp(η[2])
        Q      = tvq * exp(η[3])
        Vp     = tvvp * exp(η[4])
        Ka     = tvka * exp(η[5])
        Kin    = tvkin * exp(η[6])
        Kout   = tvkout * exp(η[7])
        Imax   = tvImax
        IC50   = tvIC50 * exp(η[8])

        σ_prop_pk = Σ[1]
        σ_prop_pd = Σ[2]
    end

    @init begin
      Response = Kin/Kout
    end

    @vars begin
      Central_conc = (Central/Vc)*1000
    end

    @dynamics begin
      Depot'      = -Ka*Depot
      Central'    =  Ka*Depot - (CL/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral 
      Peripheral' =                               (Q/Vc)*Central - (Q/Vp)*Peripheral
      Response'   =  Kin*(1-((Imax * Central_conc)/(IC50 + Central_conc))) - Kout*Response
    end

    @derived begin
      Cp  = @. Central_conc
      DV  ~ @. Normal(Cp, Cp*σ_prop_pk)

      Resp    = @. Response
      DV_Resp ~ @. Normal(Resp, Resp*σ_prop_pd)
    end
end

##Simulation
final_params = (tvcl      = 11.3,
                tvvc      = 37.6,
                tvq       = 31.4,
                tvvp      = 47.8,
                tvka      = 0.8,
                tvkin     = 20,
                tvkout    = 0.8,
                tvImax    = 1,
                tvIC50    = 2000,
                Ω         = Diagonal([0.0625,0.09,0.0625,0.09,0.09,0.16,0.09,0.09]),
                Σ         = [0.1  0.02
                              0.02 0.2])

N=30
ev0  = DosageRegimen(500, time=0, cmt=1)
protocol_sampling_time = [0.20,1.33,1.66,1.99,2.95,6.96,7.30,7.63,19.22]
pop0 = map(i -> Subject(id=i, events=ev0, time=protocol_sampling_time), 1:N)
pop0_df = DataFrame(pop0)

Random.seed!(123)
sim_base  = simobs(true_model_pk_delay_pd_sim, 
                    pop0,   
                    final_params)
df_base_sim   = DataFrame(sim_base)


##Estimation
initial_params = (  tvcl      = 10.3,
                    tvvc      = 39.6,
                    tvq       = 35.4,
                    tvvp      = 50.8,
                    tvka      = 0.8,
                    tvkin     = 20,
                    tvkout    = 0.7,
                    tvImax    = 0.8,
                    tvIC50    = 2100,
                    Ω         = Diagonal([0.0625,0.09,0.0625,0.09,0.09,0.16,0.09,0.09]),
                    Σ          = [0.1  0.02
                                  0.02 0.2])

read_fit_df = read_pumas(df_base_sim,
                        observations = [:DV, :DV_Resp])

fit_foce = fit(true_model_pk_delay_pd_sim, read_fit_df, initial_params, 
                        FOCE(), ensemblealg=EnsembleThreads())

Error

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     2.810998e+03     2.300546e+04
 * time: 2.5987625122070312e-5
ERROR: MethodError: no method matching iterate(::Nothing)
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at range.jl:872
  iterate(::Union{LinRange, StepRangeLen}, ::Integer) at range.jl:872
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at dict.jl:712
  ...
Stacktrace:
  [1] indexed_iterate(I::Nothing, i::Int64)
    @ Base ./tuple.jl:91
  [2] _unidentified_exception(A::PDMats.PDMat{Float64, Matrix{Float64}}, i::Int64, j::Int64, d::Int64, k::Symbol)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3510
  [3] _check_zero_gradient(vparam::Vector{Float64}, g::Vector{Float64}, fixedtrf::TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3475
  [4] (::Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}})(state::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:4047
  [5] #588
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3000 [inlined]
  [6] update!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, iteration::Int64, f_x::Float64, grnorm::Float64, dt::Dict{Any, Any}, store_trace::Bool, show_trace::Bool, show_every::Int64, callback::Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}, trace_simplex::Bool)
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/utilities/update.jl:23
  [7] update!
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/utilities/update.jl:11 [inlined]
  [8] trace!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, iteration::Int64, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, curr_time::Float64)
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/solvers/first_order/bfgs.jl:191
  [9] optimize(d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}})
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/optimize/optimize.jl:50
 [10] optimize
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/optimize/optimize.jl:36 [inlined]
 [11] DefaultOptimizeFN
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:2993 [inlined]
 [12] _fit(m::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, tvka = 1, tvkin = 1, tvkout = 1, tvImax = 1, tvIC50 = 1, Ω = 8, Σ = 3), 8, (:Depot, :Central, :Peripheral, :Response), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, population::Vector{Subject{NamedTuple{(:DV, :DV_Resp), Tuple{Vector{Union{Missing, Float64}}, Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, param::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, approx::FOCE, ensemblealg::EnsembleThreads, optimize_fn::Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, fixedparamset::ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, fixedparam::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_vrandeffsorth::Vector{Vector{Float64}}, verbose::Bool)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:4062
 [13] __fit
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3804 [inlined]
 [14] fit(m::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, tvka = 1, tvkin = 1, tvkout = 1, tvImax = 1, tvIC50 = 1, Ω = 8, Σ = 3), 8, (:Depot, :Central, :Peripheral, :Response), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, _population::Vector{Subject{NamedTuple{(:DV, :DV_Resp), Tuple{Vector{Union{Missing, Float64}}, Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}}, Vector{Float64}}}, param::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, approx::FOCE; optim_alg::Nothing, optim_options::Nothing, optimize_fn::Nothing, constantcoef::NamedTuple{(), Tuple{}}, omegas::Tuple{}, ensemblealg::EnsembleThreads, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_randeffs::Nothing, init_vrandeffsorth::Nothing, verbose::Bool)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3768
 [15] top-level scope
    @ ~/data/code/Methodology-Project/model/L2 Correlation/l2_correlation.jmd:161

storopoli · August 7, 2023, 11:10am

This is probably an indexing issue, since Σ is a matrix. Try:

σ_prop_pk = Σ[1, 1]
σ_prop_pd = Σ[2, 2]

benjaminrich · August 7, 2023, 1:03pm

I don’t think this is the right approach. If you want the errors for DV and DV_Resp to be correlated, you would need to use a multivariate normal distribution for the error vector, something like this:

ϵ ~ @. MvNormal(Σ)

Cp  = @. Central_conc
DV  = @. Cp*(1 + ϵ[1])

Resp    = @. Response
DV_Resp = @. Resp*(1 + ϵ[2])

I’m not sure if this is the exact way to write it in Pumas, maybe @storopoli can help?

storopoli · August 7, 2023, 7:28pm

@andreasnoack and I tried to do a minimum working example of such correlated error models. It would look something like this:

mod = @model begin
  @param begin
    σ_PK in RealDomain(; lower=0)
    σ_PD in RealDomain(; lower=0)
    ρ in RealDomain()
  end

  @derived begin
    DV ~ @. MvNormal([σ_PK^2 σ_PK*σ_PD*ρ; σ_PK*σ_PD*ρ σ_PD^2])
  end
end

Of course you would need to supply your DV as a vector of length-2.

However, this does not currently works in Pumas. We’ll support this type of model in a future Pumas release.

parssh.mehta · August 8, 2023, 1:04am

Thanks @storopoli , that works.

jogagobburu · August 8, 2023, 1:10pm

@parssh.mehta Hope you are doing well. Glad to see the model specification part is taken care of.

I am unable to appreciate the scientific basis for estimating a correlation between the random effects of PK and PD, albeit being bioanalyzed from the same sample. Even if one were to assume PD is say another concentration of some endogenous biomarker- can you please explain the scientific rationale for attempting a correlation of residual errors of PK and PD?
J

parssh.mehta · August 9, 2023, 7:36pm

Hi @jogagobburu I’m doing great, hope you are doing well too.

The current project that I am working on is impact of assessing the influence of systematic errors on Pop PK/PD model. The scientific basis for imputing this correlation is that historically the assay variability contributed the highest to the residual error, but that is not the case these days with improvement in analytical techniques. But still the residual error reported in literature is greater than the assay error and since we are modeling clinical trial perturbations (dose variability, dosing time, sampling time, and bioanalytical errors) assuming both scenarios of PK and PD samples collected at the same time and also collected at different times. If PK and PD are collected at the same time, in the same sample, we would assume that if any errors in dose variability, dosing time, sampling time, and bioanalytical errors can impact both PK and PD and not one. This can contribute that these errors are not independent from PK or PD if they are collected at the same time. Moreover, If we stay at the same sample, any misspecification in PK can propagate to PD and this shows the dependence between these dynamics. Thus we wanted to see if if there is a correlation which can bias our estimation of the PK and PD parameters from the various perturbations.

Happy to hear your thoughts on this any feedback is appreciated.

Parssh

mjaber · August 9, 2023, 7:53pm

Hi @storopoli,

Thanks for investigating this. In simulation exercise, things are promising with this implementation, not yet validated.

model = @model begin
    @param begin
        θ in VectorDomain(6)
        Ω in PDiagDomain(6)
        Σ in PSDDomain(2)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        CL = θ[1] * exp(η[1])
        V = θ[2] * exp(η[2])
        ka = θ[3] * exp(η[3])

        kin = θ[4] * exp(η[4])
        kout = θ[5] * exp(η[5])
      
        IC = θ[6] * exp(η[6])
    end
    @init begin
        y3 =kin/kout
    end

    @dynamics begin
        y1' = - ka * y1 
        y2' = ka * y1 - (CL/V) * y2 
        y3'  = kin * (1 - y2/(y2 + IC)) - kout * y3  #'
    end

    @derived begin
        cp = @. log(y2/V + 0.001) 
        pd = @. log(y3 + 0.01)  

        vec := @. MvLogNormal([cp, pd], Σ)
        dv1 ~ @. vec[1]
        dv2 ~ @. vec[2]
    end
end

dose = 500 
param = (θ = [8.0, 25.0, 0.8, 20.0, 0.2, 1.3], 
        Ω = Diagonal(ones(6) * 0.1), 
        Σ = [0.2 0.1; 0.1 0.3], )

ev = DosageRegimen(dose)       
sub = [Subject(;id=i, events=ev) for i in 1:50]
sim = simobs(model, sub, param, obstimes=0:24) |> DataFrame

However, as you mentioned read_pumas need to be modified to handle this type of data. I also want to ask about your first post, in pre chunk. What does pumas do under the hood? How it can handle the covariance or calculate the correlation since it sounds like you are using standard deviation rather than variance term?

This implementation can be helpful not only in PK/PD systems, but with parent/metabolite and plasma/urine when we do not assume independence (or at least independence is not valid assumption).

storopoli · August 10, 2023, 12:52pm

We don’t know how this would be implemented.
But there are some issues with the code that you’ve presented.

Every ~ in the derived block is some sort of “likelihood statement”. Hence you would need to use ~ with any distribution. Hence it is analogous to a probabilistic statement (since you are familiar to PPLs).
= assignments are only deterministic statements, thus, cannot have any distribution on the right hand side.

so this would probably become something like this:

DVs ~ @ MvLogNormal([cp, pd], Σ)
dv1 = @. vec[1]
dv2 = @. vec[2]

And in the read_pumas you would need to pass a column of length-2 vector elements to the observations kwarg:

pop = read_pumas(
  df;
  observations=[:DVs]
  ...
)

mjaber · August 11, 2023, 3:44am

Thanks Jose! Does it matter when we conduct simulations?

jogagobburu · August 11, 2023, 10:53am

@parssh.mehta Thank you. I would be curious to see what your research teaches us. Still wrapping my head around what practical problem such a correlation of random effects will solve (ignoring the pain of simultaneous PKPD modeling in the real world). One potential application would be for performing more realistic clinical trial simulations (CTS) - however handling residual error in CTS is by itself a challenging topic. Once you have some findings, we can discuss if the juice is worth the squeeze!
On a personal note, thank you for using Pumas for such innovative research - it is built for power users such as yourself!
J

Topic		Replies	Views
Non linear kinetics flow dependant Concepts	11	1232	April 20, 2020
Target Mediated Drug Disposition (TMDD)	4	839	April 18, 2020
Error while running working example	8	674	November 12, 2019
PK model with transit compartments for absorption	4	301	January 31, 2022
Fix data error,	3	479	November 12, 2019

L2 Correlations

Related Topics