# 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)
ΟΒ²_prop     β RealDomain(lower=0)
end

@random begin
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
``````

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.

It is not clear to me how you got the data frame?

It appears to me, that you are maybe outputting the DataFrame of an icoef? However, that seems to be only evaluated at 0 because your subject only has constant covariates, and icoef returns the output a the covariate times. Please see the docstring that explains this

`````` icoef(fpm::AbstractFittedPumasModel; obstimes)::Vector{ConstantInterpolationStructArray}

Return the individual coefficients from fpm. The individual coefficients are the variables defined in the @pre block in the @model macro. The function returns a vector of covariate interpolation objects defined on obstimes. If no obstimes are
given, the time points used to define the covariates will be used. Each of them can be evaluated at any time point, but the interpolant is piecewise constant between the obstimes. Each of the covariate interpolation objects can be converted to a
DataFrame by calling the DataFrame constructor on each element, or a DataFrame with all subjects can be obtained by calling the DataFrame constructor on the vector output by this function.
``````

You can try to increase the time resolution by supplying an obstimes vector and Iβm sure youβll see the that it switches after st_emp and again after st_emp+dur