Non linear kinetics flow dependant

monika · April 19, 2020, 12:19pm

Hello, can anyone help me with code this is a multi compartment model and the clearance is concentration dependant and also they have given a proportionality constant (a) for the changes in plasma concentration, which i have given in dynamics but still not getting the results.

PK24= @model begin
  @param begin
    tvvc    ∈ RealDomain(lower=0)
    tvclo   ∈ RealDomain(lower=0)
    tvcl    ∈ RealDomain(lower=0)
    tvq1    ∈ RealDomain(lower=0)
    tvq2    ∈ RealDomain(lower=0)
    tvvp1   ∈ RealDomain(lower=0)
    tvvp2   ∈ RealDomain(lower=0)
    tva     ∈ RealDomain(lower=0)
    Ω       ∈ PDiagDomain(8)
    σ_prop  ∈ RealDomain(lower=0)
  end
  @random begin
    η ~ MvNormal(Ω)
  end
  @pre begin
    Vc    = tvvc*exp(η[1])
    CL    = tvcl*exp(η[2])
    CLo   = tvclo*exp(η[3])
    Q1    = tvq1*exp(η[4])
    Q2    = tvq2*exp(η[5])
    Vp1   = tvvp1*exp(η[6])
    Vp2   = tvvp2*exp(η[7])
    A     = tva*exp(η[8])

  end

  @dynamics begin
    Central' =  - (CLo -(A*(Central/Vc))) * (Central/Vc) - Q1*(Central/Vc) + Q1*(Shallow/Vp1) - Q2*(Central/Vc) + Q2*(Deep/Vp2)
    Shallow' =    Q1*(Central/Vc) - Q1*(Shallow/Vp1)
    Deep'    =    Q2*(Central/Vc) - Q2*(Deep/Vp2)
  end
  
  @derived begin
    cp = @. 1000*Central/Vc
    dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
  end
end

jogagobburu · April 19, 2020, 4:49pm

What are the initial conditions? I dont see them…J

monika · April 19, 2020, 5:49pm

Actually sir, I tried adding the initial conditions but might be i was doing wrong it was giving the concentration in negative.

monika · April 19, 2020, 5:59pm

Cl = Cl0 - a* plasma concentration
this was the equation which was given in the book.

ramya · April 19, 2020, 5:59pm

Dear Monika,
Your dynamics is balanced using Proportionality constant for Clearance in central compartment. Im able to produce results using same model. Kindly check the other areas of simulation.

PK24= @model begin
@param begin
tvvc ∈ RealDomain(lower=0)
tvclo ∈ RealDomain(lower=0)
tvq1 ∈ RealDomain(lower=0)
tvq2 ∈ RealDomain(lower=0)
tvvp1 ∈ RealDomain(lower=0)
tvvp2 ∈ RealDomain(lower=0)
tva ∈ RealDomain(lower=0)
Ω ∈ PDiagDomain(7)
σ_prop ∈ RealDomain(lower=0)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
Vc = tvvcexp(η[1])
CLo = tvcloexp(η[2])
Q1 = tvq1exp(η[3])
Q2 = tvq2exp(η[4])
Vp1 = tvvp1exp(η[5])
Vp2 = tvvp2exp(η[6])
A = tva*exp(η[7])

end

@dynamics begin
Central’ = - (CLo -(A*(Central/Vc))) * (Central/Vc) - Q1*(Central/Vc) + Q1*(Shallow/Vp1) - Q2*(Central/Vc) + Q2*(Deep/Vp2)
Shallow’ = Q1*(Central/Vc) - Q1*(Shallow/Vp1)
Deep’ = Q2*(Central/Vc) - Q2*(Deep/Vp2)
end

@derived begin
cp = @. (Central/Vc)
dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
end
end
param= (tvvc=0.68,
tvclo=6.61,
tvq1=5.94,
tvq2=0.93,
tvvp1=1.77,
tvvp2=3.18,
tva=0.0025,
Ω= Diagonal([0.001,0.001,0.001,0.001,0.001,0.001,0.001]),
σ_prop=0.01)

ev1=DosageRegimen(10000,time=0,cmt=1,duration=2)
s1=Subject(id=1,evs=ev1)
obs=simobs(PK24, s1, param, obstimes=[0.25, 0.5, 0.75, 1.05, 1.25, 1.49, 1.75, 1.99, 2.16, 2.35, 2.4, 2.65, 2.81, 2.95, 3.11, 3.56, 4.15, 6, 7])
plot(obs)

![pk24|576x384](upload://pzb2weYXesfCEoEZQRlDo4R8qsH.png)

monika · April 19, 2020, 6:01pm

Yes Ramya, the plot is coming but the Cmax if you see can you tell me at which concentration your Cmax is coming?

ramya · April 19, 2020, 6:03pm

Cmax is 1000ug.

pk24

monika · April 19, 2020, 6:20pm

I am not getting plot and also in the diagonal you have used 0.001 which is not as mentioned in book and in sigma also

ramya · April 19, 2020, 6:23pm

Hello all,
In the above model 10mg of dose is given. If I scale it in derived block my cmax is approx. 600ug. And if I did not scale it in derived and give it as 10000ug in dosage regimen then cmax is 1000. what may be the reason?. Can anyone help me with this?

@derived begin
cp=@. 1000* (Central/Vc)
end

ev1=DosageRegimen(10000,time=0,cmt=1,duration=2)
s1=Subject(id=1,evs=ev1)

jogagobburu · April 19, 2020, 8:41pm

Dear Monica - the following modified code seems to produce some results, not sure if that is what you were looking for. My plots are similar to those shown above. Few pointers for consideration:

I think there was an issue with the character you used to denote derivative. When I copied & pasted your code, I had to change the '. You may want to check if its the right character.
Please include notes block upfront, provide details about what the model is supposed to do, definition of the parameters, units etc. This is very important for documentation, and future reference.
Technically, it is not clear to me what you were trying to achieve with - (CLo -(A*(Central/Vc))). CLo seems to be a zero-order elimination; A seems to be first-order. If that is the case, then what the is the rationale to call the first-order clearance A? I suggest you use mechanistically identifiable nomenclature. Scientifically, what is not clear is if the zero-order and first-order processes are separate physiologic processes or not? Typically concentration-dependent clearance (this is a better term than ‘nonlinear’ kinetics), is well dependent on the concentration range. Lower conc - first order, higher conc - zero-order. Not sure how both processes are at play simultaneously. Again, notes block upfront will be helpful esp for readers like me.

using Pumas
using Plots
######################################
##### Definitions and Units ###############
#####
######################################
PK24= @model begin
    @param begin
        tvvc ∈ RealDomain(lower=0)
        tvclo ∈ RealDomain(lower=0)
        tvq1 ∈ RealDomain(lower=0)
        tvq2 ∈ RealDomain(lower=0)
        tvvp1 ∈ RealDomain(lower=0)
        tvvp2 ∈ RealDomain(lower=0)
        tva ∈ RealDomain(lower=0)
        Ω ∈ PDiagDomain(7)
        σ_prop ∈ RealDomain(lower=0)
    end
    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Vc = tvvc*exp(η[1])
        CLo = tvclo*exp(η[2])
        Q1 = tvq1*exp(η[3])
        Q2 = tvq2*exp(η[4])
        Vp1 = tvvp1*exp(η[5])
        Vp2 = tvvp2*exp(η[6])
        A = tva*exp(η[7])
    end

#    @init begin
#       Central    = 0
#        Shallow    = 0
#        Peripheral = 0
#    end

    @dynamics begin
        Central' = - (CLo -(A*(Central/Vc))) * (Central/Vc) - Q1*(Central/Vc) +
                        Q1*(Shallow/Vp1) - Q2*(Central/Vc) + Q2*(Deep/Vp2)
        Shallow' = Q1*(Central/Vc) - Q1*(Shallow/Vp1)
        Deep'    = Q2*(Central/Vc) - Q2*(Deep/Vp2)
    end

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

param= (tvvc=0.68,
        tvclo=6.61,
        tvq1=5.94,
        tvq2=0.93,
        tvvp1=1.77,
        tvvp2=3.18,
        tva=0.0025,
        Ω= Diagonal([0.001,0.001,0.001,0.001,0.001,0.001,0.001]),
        σ_prop=0.01)

ev1=DosageRegimen(10000,time=0,cmt=1,duration=2)
s1=Subject(id=1,evs=ev1)
obs=simobs(PK24, s1, param, obstimes=[0.25, 0.5, 0.75, 1.05, 1.25, 1.49, 1.75,
                            1.99, 2.16, 2.35, 2.4, 2.65, 2.81, 2.95, 3.11, 3.56, 4.15, 6, 7])
plot(obs)

monika · April 20, 2020, 5:25am

Okay sir, I will remember this for future to make note.
So, I will try to explain what the model is all about, the aim is to model non linear kinetics with respect to flow. An IV infusion dose of 10mg/kg is given to a subject. The drug is displaying muti-compartment disposition showing perfusion limited behaviour and also known to reduce cardiac output and hepatic blood flow when plasma concentration increases.
In the exercise they have not included any parameter related to blood flow but instead they have given the above equation i.e. CL = CL0 -A * Cp, where A is proportionality constant between CL and Cp and it accounts for the change in the plasma concentration.
I have included this equation in the dynamics to correlate the change is plasma concentration.
When I tried to run the model the plot was coming but the Cmax was not reaching to 1000 and I was looking for this only that where I am doing wrong or am I missing something.

parssh.mehta · April 20, 2020, 7:20am

Dear Ramya,

In this example you are using A as a proportionality constant between Cl and cp. The dose that you are giving is in mg and after all the necessary derivations are done you are scaling it by factor of 1000. So only your final output of cp is scaled by the necessary factor.

If you want the similar output you will need to scale you proportional constant (A) accordingly in the parameters to account for the necessary concentration. Hope this helps.

Topic		Replies	Views
Target Mediated Drug Disposition (TMDD)	4	839	April 18, 2020
PK model with transit compartments for absorption	4	301	January 31, 2022
Non Linear Model - Phenytoin Simulation	3	601	February 18, 2020
L2 Correlations	10	224	August 11, 2023
Parent-metabolite + dual absorption: fitting error	5	86	November 28, 2023

Non linear kinetics flow dependant

Related Topics