Estimating lag time in my model

I am working on creating a model that uses an equation to estimate transdermal flux over time which is then fed into the dynamics block to get plasma concentrations. I want to estimate the parameters that go into the flux equation. I need the absorption to be zero up to some time point (lag) then follow a simple equation. My model runs, but it is not finding a good lag time for some subjects.
When I look at inspect(fit), I can see that the flux being used in @dynamics is negative. Why is my equation in @pre not making flux zero before the lag time? Is there a better way to do this?
See here how the lag time is estimated well for all subjects except subject 2-NOOCC.


Here is my current model:

model_fabs_noocc = @model begin
    @param begin
        tv_lag   ∈  RealDomain(lower=0.1, upper=6) 
        tv_R0    ∈  RealDomain(lower=30,upper=100) 
        tv_kout  ∈  RealDomain(lower=0.01,upper=10) 

        tv_vc    ∈  RealDomain(lower=0) 
        tv_cl    ∈  RealDomain(lower=0) 
        tv_vp1   ∈  RealDomain(lower=0) 
        tv_qp1   ∈  RealDomain(lower=0) 
        tv_vp2   ∈  RealDomain(lower=0) 
        tv_qp2   ∈  RealDomain(lower=0) 

        Ω²_lag   ∈  RealDomain(lower=0.001, upper=0.9) 
        Ω²_R0    ∈  RealDomain(lower=0.001) 
        Ω²_kout  ∈  RealDomain(lower=0.001) 

        σ²_add   ∈  RealDomain(lower=0.001) 
        σ²_prop  ∈  RealDomain(lower=0.001) 
    end

    @random begin
        η_R0    ~ Normal(0,sqrt(Ω²_R0))
        η_kout  ~ Normal(0,sqrt(Ω²_kout))
        η_lag   ~ Normal(0,sqrt(Ω²_lag))
    end
 
    @covariates CATEGORY 

    @pre begin  
        lag = tv_lag  * exp(η_lag[1])
        R   = tv_R0   * exp(η_R0[1]) 
        k   = tv_kout * exp(η_kout[1]) 
        
        flux =  R - R*exp(-k*(t-lag))
        flux .= ifelse.(flux .<= 0, 0, flux) #limits time > 0

        Vc  = tv_vc  
        CL  = tv_cl
        Vp1 = tv_vp1  
        Qp1 = tv_qp1  
        Vp2 = tv_vp2 
        Qp2 = tv_qp2 
    end

    @dynamics begin
        Central'     = flux - (CL/Vc)*Central - (Qp1/Vc)*Central + (Qp1/Vp1)*Peripheral1 - (Qp2/Vc)*Central + (Qp2/Vp2)*Peripheral2 
        Peripheral1' =                          (Qp1/Vc)*Central - (Qp1/Vp1)*Peripheral1
        Peripheral2' =                                                                     (Qp2/Vc)*Central - (Qp2/Vp2)*Peripheral2
    end

    @derived begin
        Cp      =  @. Central/Vc 
        CONC ~ @. Normal(Cp, sqrt(((abs(Cp)^2)*σ²_prop) + σ²_add))

        flux_dv = @. flux
    end
end

@KVTobin333 have you considered using the more robust delay functionality?

Also, the way you coded flux almost seems like this is an infusion. How does your dataset look for dosing?

I will look into using that macro, but it is best if I can use this exact equation for flux since the parameters in this equation give us information about the delivery. I am running with model with no dosing events.

Yes, this is for a transdermal product. So it is similar to infusion, but with a changing rate over time described by the flux equation. It would be similar to an infusion starting at a slow rate as the rate is slowly increased to a constant rate. But I don’t know the exact so need to deconvolute them from the plasma concentration data.
All patients were given the same dose:1 mg applied to skin.

Over-writing flux like this is not allowed. Try replacing the 2 flux lines with:

flux =  max(0, R - R*exp(-k*(t-lag)))

Also using . in this case is not necessary because the code in the @pre block is evaluated at every time point and the numbers here are all scalars.

That makes sense. I will also try this, but I suspect it won’t be able to calculate the bsv on lag since it will be zero for the beginning of the curve and won’t know when to stop it from being zero. When I had my flux in one line like this but using ifelse, it always gave an error saying omega_lag was not defined.