# Enterohepatic recirculation

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.

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(η)
Vc          = tvvc * exp(η)
Vp          = tvvp * exp(η)
Q           = tvQ * exp(η)
Ka          = tvka * exp(η)
Klg         = tvklg * exp(η)
τ           = tvτ * exp(η)
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
``````

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.

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,

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.

1 Like

Will check on Patrick’s suggestion.

1 Like

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

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

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)

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.