Enterohepatic recirculation

[MathangiCTM]

Hello,
I am trying to model a drug undergoing enterohepatic recirculation like shown below. The emptying of drug from bile to gut only happens during a time interval, ie., ETime to ETime+Dur_empty. Either before ETime or after ETime+Dur_empty, the empty from bile compartment is not allowed.

In NONEME, the MTIME function does the turning off and on and I found a related discussion here.:
https://www.mail-archive.com/nmusers@globomaxnm.com/msg01349.html

Any help to implement the EHC model in Pumas is appreciated.
Thanks.

[parssh.mehta]

@MathangiCTM

EHC model in Pumas is appreciated.

Hi Mathangi,

I might have a similar example from the book by Drs.Gabrielsson and Weiner. (PK40). Hope this helps.

pk_40            = @model  begin
   @param begin
     tvcl        ∈ RealDomain(lower=0)
     tvvc        ∈ RealDomain(lower=0)
     tvvp        ∈ RealDomain(lower=0)
     tvQ         ∈ RealDomain(lower=0)
     tvka        ∈ RealDomain(lower=0)
     tvklg       ∈ RealDomain(lower=0)
     tvτ         ∈ RealDomain(lower=0)
     Ω           ∈ PDiagDomain(7)
     σ²_prop     ∈ RealDomain(lower=0)
   end

   @random begin
     η           ~ MvNormal(Ω)
   end

   @pre begin
     Cl          = tvcl * exp(η[1])
     Vc          = tvvc * exp(η[2])
     Vp          = tvvp * exp(η[3])
     Q           = tvQ * exp(η[4])
     Ka          = tvka * exp(η[5])
     Klg         = tvklg * exp(η[6])
     τ           = tvτ * exp(η[7])
     Kempt       = (t>10 && t<(10+τ))*(1/τ)
   end

   @dynamics begin
     Central'    = Ka*Depot - (Cl/Vc)*Central + (Q/Vp)*Peripheral - (Q/Vc)*Central - Klg*Central
     Peripheral' = (Q/Vc)*Central - (Q/Vp)*Peripheral
     Bile'       = Klg*Central - Bile*Kempt
     Depot'      = Bile*Kempt - Ka*Depot
   end

   @derived begin
     cp          = @. Central/Vc
     dv          ~ @. Normal(cp, sqrt(cp^2*σ²_prop))
   end
end

[storopoli]

Currently, there is no way to tell the underlying ODE solver about τ.

The model that @parssh.mehta kindly contributed might not finish the fit or have problems with any successful fit.

We will sure implement a way to run this model in Pumas.

[parssh.mehta]

Thanks for the info @storopoli , just curious to know what you mean by there is no way of telling the underlying ODE solver about τ.

Like in many cases where we have an infusion, we estimate the duration of the infusion, which still does not go directly into the ODE solver. Just wondering how both these cases are different?

Thanks,

Parssh

[patrick]

Usually, the emptying of the gall bladder is modeled as a sudden release - in nonmem often controlled by “mtime”. mtime allows the user to programmatically add events to the subjects at time points controlled by a fixed effect (and maybe a random effect to add offsets that account for different meal intake etc).

In Pumas we have no way of adding these events. We can control a switch of parameters by comparing t to some number in the pre block, but if we do that, we should tell the ODE solver about this point in time, because it causes a non-differentiability in the system ( a discontinuity in the right hand side of X’ = … equations in the dynamics block). We can add tstops in simobs and fit, but the problem is that they would have to be the same for all subjects the way Pumas works right now. This is not good, because this means that we cannot vary it throughout the fit for example, and if we do not tell the solver about these discontinuities in the time derivatives, the solution can be quite wrong and you may end up with warnings and errors.

We are aware of a request from users to have a feature like mtime. In the mean time, if the project allows it, maybe consider a more model based approach with a gall bladder emptying controlled by a continuous onset function ? Benjamin Guiastrennec’s dissertation summary called “Mechanism-based modeling of biological processes involved in oral absorption” describes one approach (found here https://www.diva-portal.org/smash/get/diva2:1178146/FULLTEXT01.pdf) (see eq 8 for the onset function ) but I’m sure there are others.

[MathangiCTM]

Thank you @parssh.mehta @patrick @storopoli.
Will check on Patrick’s suggestion.

[patrick]

The paper is here https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5192972/pdf/PSP4-5-692.pdf

Again, I’m not sure a similar approach will be applicable in your situation, but maybe it’s worth investigating.

[parssh.mehta]

Thanks @patrick for the detailed explanation. That’s very helpful.

[MathangiCTM]

Hi,
I am still in the process of exploring the enterohepatic recirculation (EHC) model depicted in my original message. I used a conditional time statement to turn on and off the bile compartment and was able to simulate concentration - time profile emulated the EHC. The simulation code for is below:

ehcmodel = @model begin
    @param begin
        tvcl            ∈   RealDomain(lower=0.0001)
        tvvc            ∈   RealDomain(lower=0.0001)
        tvka            ∈   RealDomain(lower=0.0001)
        tvkpb           ∈   RealDomain(lower=0.0001)
        tvkgp           ∈   RealDomain(lower=0.0001)
        tvk0bile        ∈   RealDomain(lower=0.0001)
        st_emp          ∈   RealDomain(lower=0.0001)
        dur             ∈   RealDomain(lower=0.0001)
        # RUV
         σ_prop          ∈   RealDomain(lower=0.0001)
    end
    
    @pre begin
    Cl          = tvcl 
    Vc          = tvvc 
    Ka          = tvka 
    Kpb         = tvkpb 
    Kgp         = tvkgp 
    k0bile      = (t < st_emp || t > st_emp+dur) ? 0 : tvk0bile
    end

    @dynamics begin
    Depot'      = - Ka*Depot
    Central'    = Ka*Depot - (Cl/Vc)*Central - Kpb*Central + Kgp * Gut
    Bile'       = Kpb*Central - k0bile
    Gut'        = k0bile - Kgp * Gut
    end

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

params = (;
         tvcl = 4.6,
         tvvc = 15,
         tvka = 4.0,
         tvkpb = 0.5,
         tvkgp = 0.5,
         tvk0bile = 20,
         st_emp = 2,
         dur = 3, 
         σ_prop = 0.2,
)

## Dosing regimen
ev = DosageRegimen(75; time=0)

## subjects
s1 = Subject(;
	id=1,
	events=ev,
	#covariates=(; Wt=70),
	observations=(; dv=nothing),
)
s2 = Subject(;
	id=2,
	events=ev,
	#covariates=(; Wt=70),
	observations=(; dv=nothing),
)
twosubjs = [s1,s2]

## SImulation

sims1 = simobs(ehcmodel, twosubjs, params; obstimes=[0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 5,5.5, 6, 8, 10, 12, 14, 16, 18, 24])

I was also able to reestimate the parameters using a NaivePooled approach.
The estimation model is below:

## Estimation of EHC model parameters
check_ehc = DataFrame(sims1)
check_ehc[!, :cmt] .= ifelse.(check_ehc.evid .== 1, 1, 2)
@rsubset! check_ehc !(ismissing(:amt) && :time == 0)

checkehc_fitdata        =   read_pumas(check_ehc,
id                      =   :id,
time                    =   :time,
amt                     =   :amt,
cmt                     =   :cmt,
observations   =   [:dv],
evid                    =   :evid)

ehcmodel_est = @model begin
    @param begin
        tvcl            ∈   RealDomain(lower=0.0001)
        tvvc            ∈   RealDomain(lower=0.0001)
        tvka            ∈   RealDomain(lower=0.0001)
        tvkpb           ∈   RealDomain(lower=0.0001)
        tvkgp           ∈   RealDomain(lower=0.0001)
        tvk0bile         ∈   RealDomain(lower=0.0001)
        st_emp          ∈   RealDomain(lower=0.0001)
        dur             ∈   RealDomain(lower=0.0001)
         # RUV
         σ_prop          ∈   RealDomain(lower=0.0001)
    end

    @pre begin
    Cl          = tvcl 
    Vc          = tvvc 
    Ka          = tvka 
    Kpb         = tvkpb 
    Kgp         = tvkgp 
    k0bile       = (t < st_emp || t > st_emp+dur) ? 0 : tvk0bile
    end

    @dynamics begin
    Depot'      = - Ka*Depot
    Central'    = Ka*Depot - (Cl/Vc)*Central -  Kpb*Central + Kgp*Gut
    Bile'       = Kpb*Central - k0bile
    Gut'        = k0bile - Kgp * Gut
    end

    @derived begin
    cp          = @. Central/Vc
    dv          ~ @. Normal(cp, sqrt((cp*σ_prop)^2))
    end
end
params_ehc = (;
         tvcl = 2.5,
         tvvc = 15,
         tvka = 2.5,
         tvkpb = 0.5,
         tvkgp = 0.5,
         tvk0bile = 5,
         st_emp = 2,
         dur = 3, 
         σ_prop = 0.2
)

pkfit_1cmt_ehc_np = fit(ehcmodel_est,
                 checkehc_fitdata,
                 params_ehc,
                 Pumas.NaivePooled(),
                 constantcoef = (st_emp = 2.0, dur = 3.0,))

I get reasonable estimates and the fit seems adequate.

However, when I do a two stage analysis, k0bile is not estimated rather it takes the value 0 as mentioned in the conditional statement.

## Zero order bile - 2 stage      
pkfit_2cmt_ehc = [fit(ehcmodel_est,
            checkehc_fitdata[i],
            params_ehc,
            constantcoef = (st_emp = 2.0, dur = 3.0,)) for i in 1:length(checkehc_fitdata)]

2×8 DataFrame
 Row │ id      time     Cl       Vc       Ka       Kpb       Kgp       k0bile  
     │ String  Float64  Float64  Float64  Float64  Float64   Float64   Float64 
─────┼─────────────────────────────────────────────────────────────────────────
   1 │ 1           0.0  3.92074  17.4281  5.71104  0.435616  0.521555      0.0
   2 │ 2           0.0  4.34398  13.8105  2.983    0.526967  0.508399      0.0

Any reason why that could be the case with the 2-stage approach?

Thanks.