 # Non linear kinetics flow dependant

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(η)
CL    = tvcl*exp(η)
CLo   = tvclo*exp(η)
Q1    = tvq1*exp(η)
Q2    = tvq2*exp(η)
Vp1   = tvvp1*exp(η)
Vp2   = tvvp2*exp(η)
A     = tva*exp(η)

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
``````

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

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

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

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(η)
CLo = tvclo
exp(η)
Q1 = tvq1exp(η)
Q2 = tvq2
exp(η)
Vp1 = tvvp1exp(η)
Vp2 = tvvp2
exp(η)
A = tva*exp(η)

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)``

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

Cmax is 1000ug. 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

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) 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:

1. 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.
2. 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.
3. 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(η)
CLo = tvclo*exp(η)
Q1 = tvq1*exp(η)
Q2 = tvq2*exp(η)
Vp1 = tvvp1*exp(η)
Vp2 = tvvp2*exp(η)
A = tva*exp(η)
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)
``````
2 Likes

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.

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.

1 Like