Nested task error while fitting

Hi,
I was trying to fit a PD model but I am getting Domain Error.
Note: Please find the code and error.

@dynamics begin
        Central'     = -((CL/Vc)*Central)  - ((Q/Vc)*Central) - ((Q2/Vc)*Central) + ((Q/Vp)*Peripheral1) + ((Q2/Vp2)*Peripheral2)
        Peripheral1' =  ((Q/Vc)*Central)   - ((Q/Vp)*Peripheral1)
        Peripheral2' =  ((Q2/Vc)*Central)  - ((Q2/Vp2)*Peripheral2) 
        eff_central'  =  Ke0*(Central/Vc) - Ke0*eff_central + Ke21*eff_peripheral - Ke12*eff_central 
        eff_peripheral' = - Ke21*eff_peripheral + Ke12*eff_central
        res' = (E0-((Emax*((Central/Vc)^Sigma))/(((EC50)^Sigma)+((Central/Vc)^Sigma))))
    end

    @derived begin
        mean_conc := @. Central/Vc
        DV ~ @. Normal(mean_conc, sqrt(σ_add^2 + (mean_conc*σ_prop)^2))
        eff := @. res
        BIS ~ @. Normal(eff, σ_add1)
    end
  end

ERROR: TaskFailedException

nested task error: DomainError with -0.14869878833270983:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
  [1] throw_exp_domainerror(x::Float64)
    @ Base.Math ./math.jl:37
  [2] ^
    @ ./math.jl:901 [inlined]
  [3] ^
    @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/ForwardDiff/QOqCN/src/dual.jl:471 [inlined]

I have encountered such errors in the past which were due to some algebraic transformations in the @pre block. Could you share your entire model ?

Hi Rahul,
Please find the code for the entire model.

Jaya Shree

model_PD2 = @model begin
    @param begin
        tvcl   ∈ RealDomain(lower=0)
        tvvc   ∈ RealDomain(lower=0)
        tvvp   ∈ RealDomain(lower=0)
        tvvp2  ∈ RealDomain(lower=0)
        tvq    ∈ RealDomain(lower=0)
        tvq2   ∈ RealDomain(lower=0)
        θ1     ∈ RealDomain(lower=0)
        θ2     ∈ RealDomain(lower=0)
        tve0   ∈ RealDomain(lower=0)
        tvemax ∈ RealDomain(lower=0)
        tvec50 ∈ RealDomain(lower=0)
        tvsigma∈ RealDomain(lower=0)
        tvke0  ∈ RealDomain(lower=0)
        tvke12 ∈ RealDomain(lower=0)
        tvke21 ∈ RealDomain(lower=0)
        ωCL    ∈ RealDomain(lower=0.0001)
        ωVc    ∈ RealDomain(lower=0.0001)
        ωVp    ∈ RealDomain(lower=0.0001)
        ωVp2   ∈ RealDomain(lower=0.0001)
        ωQ     ∈ RealDomain(lower=0.0001)
        ωQ2    ∈ RealDomain(lower=0.0001)
        ωEC50  ∈ RealDomain(lower=0.0001)
        ωSigma ∈ RealDomain(lower=0.0001)
        ωKe0   ∈ RealDomain(lower=0.0001)
        ωKe12  ∈ RealDomain(lower=0.0001)
        ωKe21  ∈ RealDomain(lower=0.0001)
        σ_add  ∈ RealDomain(lower=0.0001)
        σ_add1  ∈ RealDomain(lower=0.0001)
        σ_prop ∈ RealDomain(lower=0.0001)
    end
  
    @random begin
        ηCL        ~ Normal(0.0, ωCL)
        ηVc        ~ Normal(0.0, ωVc)
        ηVp        ~ Normal(0.0, ωVp)
        ηVp2       ~ Normal(0.0, ωVp2)
        ηQ         ~ Normal(0.0, ωQ)
        ηQ2        ~ Normal(0.0, ωQ2)
        ηEC50      ~ Normal(0.0, ωEC50)
        ηSigma     ~ Normal(0.0, ωSigma)
        ηKe0       ~ Normal(0.0, ωKe0)
        ηKe12      ~ Normal(0.0, ωKe12)
        ηKe21      ~ Normal(0.0, ωKe21)
    end
  
    @covariates RF WT
  
    @pre begin
        CL  = (tvcl+θ1*WT) * exp(ηCL)
        Vc  = (tvvc  + θ2*(RF))  * exp(ηVc)
        Vp  = tvvp  * exp(ηVp)
        Vp2 = tvvp2 * exp(ηVp2)
        Q   = tvq   * exp(ηQ)
        Q2  = tvq2  * exp(ηQ2)   
        E0  = tve0
        Emax= tvemax
        EC50= tvec50 * exp(ηEC50)
        Sigma=tvsigma * exp(ηSigma)
        Ke0 = tvke0 * exp(ηKe0)
        Ke12= tvke12 * exp(ηKe12)
        Ke21= tvke21 * exp(ηKe21)
    end

    @dynamics begin
        Central'     = -((CL/Vc)*Central)  - ((Q/Vc)*Central) - ((Q2/Vc)*Central) + ((Q/Vp)*Peripheral1) + ((Q2/Vp2)*Peripheral2)
        Peripheral1' =  ((Q/Vc)*Central)   - ((Q/Vp)*Peripheral1)
        Peripheral2' =  ((Q2/Vc)*Central)  - ((Q2/Vp2)*Peripheral2) 
        eff_central'  =  Ke0*(Central/Vc) - Ke0*eff_central + Ke21*eff_peripheral - Ke12*eff_central 
        eff_peripheral' = - Ke21*eff_peripheral + Ke12*eff_central
        res' = (E0-((Emax*((Central/Vc)^Sigma))/(((EC50)^Sigma)+((Central/Vc)^Sigma))))
    end

    @derived begin
        mean_conc := @. Central/Vc
        DV ~ @. Normal(mean_conc, sqrt(σ_add^2 + (mean_conc*σ_prop)^2))
        eff := @. res
        BIS ~ @. Normal(eff, σ_add1)
    end
  end
```

I do not see an issue with the model. I would recommend a data quality check of all columns in the dataset that are directly/indirectly involved in the calculation of an exponent in the model.
For e.g. If the RF column is negative someplace (which we might not expect) then Vc can turn out to be -ve leading to the error which will stem from performing the operation: (Central/Vc)^Sigma

Another suggestion for writing a relatively cleaner model is to use the @vars block :

    @vars begin
        conc := Central/Vc
        e_drug := Emax*conc^Sigma/(EC50^Sigma + conc^Sigma)
    end

    @dynamics begin
        Central'     = -(CL/Vc)*Central  - (Q/Vc)*Central - (Q2/Vc)*Central + (Q/Vp)*Peripheral1 + (Q2/Vp2)*Peripheral2
        Peripheral1' =  (Q/Vc)*Central   - (Q/Vp)*Peripheral1
        Peripheral2' =  (Q2/Vc)*Central  - (Q2/Vp2)*Peripheral2
        eff_central'  =  Ke0*(Central/Vc) - Ke0*eff_central + Ke21*eff_peripheral - Ke12*eff_central 
        eff_peripheral' =                                   - Ke21*eff_peripheral + Ke12*eff_central
        res' = E0 - e_drug
    end

I have fixed all the PK parameters. Only PD parameters need to be estimated in this model.

I tried with @vars block also but I have the same error and if I remove some etas, the estimation runs but gives the same parameter given in the param block.

@jayashreed1412
Are you trying to fit a direct effect (Emax) model to the concentration-effect data? If so, I suggest the following: You do not need a differential equation for the res. Use the expression you have e_drug in the @derived block. Let me know how it goes.
J

Hi,
I tried running the fit by using expression e_drug in @derived block but the fit runs and had given the same parameters which was initially given to estimate.

Jaya Shree

@jayashreed1412

  1. If you can share your data with me, I can setup the model and share with you. If that is not possible;
  2. I am unsure why you are trying to estimate PK and PD. Only under some very unique circumstances, one might need to estimate both pk and pd parameters together (called ‘simultaneous pk and pd estimation’). You should start with sequential pk and pd estimation. That is, first estimate PK and then fix the PK parameters to drive the PD model. I can help you set it up, once I understand what you are trying to achieve.
    J

Hi J,
I’m estimating the PK and PD in a sequential manner. I had already estimated the PK parameters, and now I have fixed all PK parameters in the PD estimation. As a result, I encounter this error when estimating only the PD parameter.

So, I’d appreciate it if you could assist me with this. Can you provide me with your email address to send you my data and code?

Jaya Shree

Hi Jay Shree,

I am sharing an example. I hope this will be helpful. In this example the PKPD estimation is divided into two parts - In the first part, the PK parameters are estimated. Individual parameters from the pk fit are used for PD estimation. Individual clearance and volume of distribution are used to calculate the concentration for the PD analysis. In the second part, the PD parameters are estimated using the PD model. Please let me know if this approach works for your model.

using Pumas
using CSV
using StatsBase
using Random
using DataFramesMeta
using StatsPlots
using PumasPlots
using Statistics

pkdata          = CSV.File("iv_inf_pd_direct_add.csv"; missingstring="") |> DataFrame

pkdata.conc_obs = ifelse.(pkdata.time .== 0.0, missing, pkdata.conc_obs)

data4pkest   =  read_pumas(pkdata,
                           id           = :id,
                           time         = :time,
                           amt          = :amt,
                           rate         = :rate,
                           cmt          = :cmt,
                           evid         = :evid,
                           observations = [:conc_obs])
pk1cmt = @model begin
    @param   begin
        tvcl      ∈ RealDomain(lower=0.1)
        tvvc      ∈ RealDomain(lower=0.1)
        Ω         ∈ PDiagDomain(2)
        σ_add_pk  ∈ RealDomain(lower=0)
    end

    @random begin
        η         ~ MvNormal(Ω)
      end
  
    @pre begin
        CL        = tvcl * exp(η[1])
        Vc        = tvvc * exp(η[2])
    end

     @dynamics begin
         Central'  =  - (CL/Vc)*Central
     end

    @derived begin
        cp_ipred   = @. (Central/Vc)
        conc_obs   ~ @. Normal(cp_ipred, sqrt(σ_add_pk))
    end
end

pkparams     = (tvcl     = 1.5,
                tvvc     = 25.0,
                Ω        = Diagonal([0.09, 0.09]),
                σ_add_pk = 2.0)

pk_1cmt_results      = fit(pk1cmt, data4pkest, pkparams, Pumas.FOCE())

pk_1cmt_inspect = DataFrame(inspect(pk_1cmt_results))
pk_1cmt_icoef   = DataFrame(icoef(pk_1cmt_results))
select!(pk_1cmt_icoef, ([:id,:CL,:Vc]))
rename!(pk_1cmt_icoef,:CL => :iCL, :Vc => :iVc)
pk_1cmt_icoef.id        = map(x->parse(Int,x), pk_1cmt_icoef.id)

pkpddata_pkindvparams   = innerjoin(pkpddata, pk_1cmt_icoef, on=:id)

# PD Modeling using Individual PK estimates

data4pdest   =  read_pumas(pkpddata_pkindvparams,
                           id           = :id,
                           time         = :time,
                           amt          = :amt,
                           rate         = :rate,
                           cmt          = :cmt,
                           evid         = :evid,
                           observations = [:resp_obs],
                           covariates   = [:iCL,:iVc,:dose]
                           )

pd_add = @model begin
    @param   begin
        tvbsl    ∈ RealDomain(lower=0)
        tvec50   ∈ RealDomain(lower=0)
        tvemax   ∈ RealDomain(lower=0)
        Ω        ∈ PDiagDomain(1)
        σ_add_pd ∈ RealDomain(lower=0)
    end

    @random begin
        η        ~ MvNormal(Ω)
      end

    @covariates iCL iVc dose

    @pre begin
        CL       = iCL
        Vc       = iVc

        bsl      = tvbsl * exp(η[1])
        ec50     = tvec50
        emax     = tvemax
    end

     @dynamics begin
         Central' =  - (CL/Vc)*Central
     end

    @derived begin
        cp_ipred    = @. (Central/Vc)
        resp_ipred  = @. bsl + emax*cp_ipred/(ec50 + cp_ipred)
        resp_obs    ~ @. Normal(resp_ipred, sqrt(σ_add_pd))
    end
end

pd_params =    (tvbsl     = 10,
                tvec50    = 25,
                tvemax    = 10,
                Ω         = Diagonal([0.09]),
                σ_add_pd  = 2)

pd_fit        = fit(pd_add, data4pdest, pd_params, Pumas.FOCE())

pd_fit_infer  = infer(pd_fit)

I just checking to be sure, but are you sure that all values of WT and RF are positive?

There are a few places where you’re using a power function and they include only Central, Vc, Sigma, and EC50. So if something is negative either WT, tvcl, θ1, exp(ηCL), tvvc, RF, θ2, exp(ηVc), tvec50, exp(ηEC50), tvsigma, or exp(ηSigma) is negative. All parameters are restricted to be positive, and all random effects are in exp(), so that means that the error is in WT, RF, or Central if I’m not mistaken. If Central is negative you would have to look at your three compartment system and see if something is wrong (it looks fine), or if the solver makes it ever so slightly negative and when you divide by Vc it becomes -0.1486… if Vc < 1. In that case, you can wrap your Central variables in abs such that you res' calculation becomes

    res' = (E0-((Emax*((abs(Central)/Vc)^Sigma))/(((EC50)^Sigma)+((abs(Central)/Vc)^Sigma))))

but I would make sure that it’s not RF that’s negative first.

If you can provide more of the error message, we can see exactly where the calculation went wrong.

Hi Nikita,
Thank you for providing the example.
My credit limit in Julia is currently exceeded, so once the credit limit is increased and I gain access to Julia. I’ll give that a try and update you on the status.
Thanks,
Jaya Shree

2 Likes

Hi Patrick,
I am unable to log in to Julia due to a technical issue. Once the problem has been resolved. I will keep you updated on the same.
Thanks,
Jaya Shree

1 Like

Hello, Nikita and Patrick.
I tried using the example shared by Nikita. Now I can run the PD model. Thank you for the help.

Regards,
Jaya Shree

2 Likes