# 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
end
end

``````

``````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:
 throw_exp_domainerror(x::Float64)
@ Base.Math ./math.jl:37
 ^
@ ./math.jl:901 [inlined]
 ^
@ /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)
σ_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
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.

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)

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)
end

@random begin
η         ~ MvNormal(Ω)
end

@pre begin
CL        = tvcl * exp(η)
Vc        = tvvc * exp(η)
end

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

@derived begin
cp_ipred   = @. (Central/Vc)
end
end

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

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

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

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

@random begin
η        ~ MvNormal(Ω)
end

@covariates iCL iVc dose

@pre begin
CL       = iCL
Vc       = iVc

bsl      = tvbsl * exp(η)
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)
end
end

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

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