continuous markov process simulation

IQR Tools supports and increases the efficiency, quality, and compliance of model-based analyses in pharmacometrics, systems parmacology, and systems biology by incorporating and extending the capabilities of existing tools. The user-friendliness of the package considerably lowers the threshold for the conduct of modeling & simulation based analyses, which also makes it useful both for complex analyses and for educational purposes.

Here you can ask questions and get answers ...

Moderators: daniel.kaschek, intiquan

Post Reply
poppharm@gmail.com
Posts: 2
Joined: Thu Jun 20, 2019 10:20 am

continuous markov process simulation

Post by poppharm@gmail.com » Thu Jun 20, 2019 10:46 am

Again me as I am learning how to simulate wioth IQRtools complex models.

This is a continuous markov process I want to simulate. Easy to do with Phoenix NLME but here I have some issues.

First of all I need to reset after each simulated outcome that is a categorical some compartments to 0 and then reset values in each of these compartments, depending on the simulated outcome.

How can I reset some sate variables to 0?
How can I redefine the initial conditions of each state variable depending on the simulated outcome?

Is there a way to do it? If so, what is the algorithm?

Thanks in advance
Serge Guzy

The Phoenix code to simulate is shown below

test(){



# we put the PK model
cfMicro(A1, Cl / V, first = (Aa = Ka))
dosepoint(Aa)
C = A1 / V

stparm(Ka = tvKa * exp(nKa))
stparm(V = tvV * exp(nV))
stparm(Cl = tvCl * exp(nCl))

fixef(tvKa = c(, 0.379, ))
fixef(tvV = c(, 6.28, ))
fixef(tvCl = c(, 0.141, ))
ranef(diag(nV, nCl, nKa) = c(0.15, 0.1, 0.23))


fcovariate(timecov,dosageintervalmonth,dosing)
# Ka=(timecov<0?0:kacov)
# V=(timecov<0?0:vcov)
# Cl=(timecov<0?0:clcov)
C = (V<=0?0:A1 / V)


# we put dose as covariate
fcovariate(dose)
# 0 - nodisease, 1 - mild , 2 - moderate , 3 - severe , D - Drug effect
# deriv(nodisease = -(nodisease*K01-mild*K10*D)-(nodisease*K02-moderate*K20*D)-(nodisease*K03-severe*K30*D)
# deriv(mild = (nodisease*K01-mild*K10*D)+(moderate*K21*D-mild*K12)-(mild*K13-severe*K31*D)
# deriv(moderate = (nodisease*K02-moderate*K20*D)-(moderate*K21*D-mild*K12)-(moderate*K23-severe*K32*D)
# deriv(severe = (nodisease*K03-severeK30*D)+(moderate*K23-severe*K32*D)+(mild*K13-severe*K31*D)

# we put emax model for dose effect
# we define the drug effect 0 for Placebo group else the drug effect
# C/(C+ed50)=1/(1+ed50/C)
# We prevent to have C=0
term=(C<=0?0:1/(1+ed50/C))
de=(term==0?0:dmax*term)
# we define the placebo effect during the run in and different than after randomizaiton
#k10runin*(1+imaxrunin*t)
disease10=(k10*(1+imax*t))
disease20=(k20*(1+imax*t))
disease30=(k30*(1+imax*t))
disease21=(k21*(1+imax*t))
disease31=(k31*(1+imax*t))
disease32=(k32*(1+imax*t))
deriv(nodisease = -(nodisease * k01 - mild * disease10*(1+de))- (nodisease * k02 - moderate * disease20*(1+de)) - (nodisease *k03 -severe * disease30*(1+de)))
deriv(mild = (nodisease * k01 - mild * disease10*(1+de)) + (moderate * disease21*(1+de) - mild * k12) - (mild * k13 - severe * disease31*(1+de)))
deriv(moderate = (nodisease * k02 - moderate * disease20*(1+de)) - (moderate * disease21*(1+de) - mild * k12) - (moderate * k23 - severe * disease32*(1+de)))
deriv(severe = (nodisease * k03 - severe * disease30*(1+de)) + (moderate * k23 - severe * disease32*(1+de)) + (mild * k13 - severe * disease31*(1+de)))

fixef(tvdmax = c(, 2.09167, ))
stparm(dmax=tvdmax*exp(ndmax))
ranef(diag(ndmax)=c(1.4072371))

fixef(tvk02 = c(, 0.537678, ))
stparm(k02=tvk02*exp(nk02))
ranef(diag(nk02)=c(2.1479864))

fixef(tvk03 = c(, 0.0725198, ))
stparm(k03=tvk03*exp(nk03))
ranef(diag(nk03)=c(2.5249955))

fixef(tvk12 = c(, 0.691888, ))
stparm(k12=tvk12*exp(nk12))
ranef(diag(nk12)=c(0.34061999))

fixef(tvk13 = c(, 0.016028, ))
stparm(k13=tvk13*exp(nk13))
ranef(diag(nk13)=c(1.3463482))

fixef(tvk23 = c(, 0.654812, ))
stparm(k23=tvk23*exp(nk23))
ranef(diag(nk23)=c(1.0203568))

fixef(tvk20 = c(, 0.0505547, ))
stparm(k20=tvk20*exp(nk20))
ranef(diag(nk20)=c(1.0523402))

fixef(tvk30 = c(, 0.282396, ))
stparm(k30=tvk30*exp(nk30))
ranef(diag(nk30)=c(2.2056248))

fixef(tvk21 = c(, 0.824128, ))
stparm(k21=tvk21*exp(nk21))
ranef(diag(nk21)=c(0.3442199))

fixef(tvk31 = c(, 0.0113712, ))
stparm(k31=tvk31*exp(nk31))
ranef(diag(nk31)=c(2.0856592))

fixef(tvk32 = c(, 0.991065, ))
stparm(k32=tvk32*exp(nk32))
ranef(diag(nk32)=c(0.42273243))

fixef(tved50 = c(, 313159, ))
stparm(ed50=tved50*exp(ned50))
ranef(diag(ned50)=c(2.1327861))

fixef(tvimax = c(, 0.00826507, ))
stparm(imax=tvimax*exp(nimax))
ranef(diag(nimax)=c(1.6724364))


fixef(tvimaxrunin = c(, 0.000929115, ))
stparm(imaxrunin=tvimaxrunin*exp(nimaxrunin))
ranef(diag(nimaxrunin)=c(2.3217324))

stparm(k10 = tvk10 * exp(nk10))
stparm(k10runin = tvk10runin * exp(nk10runin))
stparm(k01 = tvk01 * exp(nk01))
# covariate(dvbase)
fcovariate(dvcovariate)
fixef(tvk10 = c(, 1.16088, ))
fixef(tvk10runin = c(, 0.706373, ))
fixef(tvk01 = c(, 0.614314, ))
ranef(diag(nk10, nk01) = c(1.5476961,0.40932173))
ranef(diag(nk10runin)=c(0.21708326))
double(u,tempheadache,tempnoheadache)

double(tempnodisease,tempmild,tempmoderate,tempsevere,prob0,prob01,prob012,dvtemp)
# define the likelihiid and initial conditions
LL(dv,log(dv==0?max(nodisease,0.00001):
dv==1?max(mild,0.00001)
:dv==2?max(moderate,0.00001)
:max(severe,0.00001))

,simulate={
# simulate(0-1) uniform
u=unif()

# # headache value (0 or 1) and nomheadache
# tempheadache=(u<=headache?1:0)
# tempnoheadache=(u>noheadache?1:0)
# dv=(tempheadache==1?1:0)

# t=0 we purt dvcovariate
# dv=(t==0?dvcovariate:dv)

# 0 - nodisease, 1 - mild , 2 - moderate , 3 - severe
# sum of all prob=1 then we define prob0 prob01 and prob012
prob0=nodisease
prob01=prob0+mild
prob012=prob01+moderate

# now we sim uniform01 and if less than prob- then dv
dv=(u<prob0?0:u<prob01?1:u<prob012?2:3)
dvtemp=dv
}


,doafter=
{nodisease=0
mild=0
moderate=0
severe=0

nodisease=(dvtemp==0?1:0)
mild=(dvtemp==1?1:0)
moderate=(dvtemp==2?1:0)
severe=(dvtemp==3?1:0)


}


)
# initial conditions(time3=0)
# we will use the cumulative probabilities seen in the data as initial conditions to sim
double(ran1,level)

sequence{

# we simulate 0-1 uniform
# 0.13333333 level0
# 0.26666667 level1
# 0.44952381 level2
# 0.15047619 level3

ran1=unif()
# we first initiate all to 0
nodisease=0
mild=0
moderate=0
severe=0
# if ran1 < 0.13 then nodisease=1, other=0
# else ran1<0.4 ?mild=1, else ran1<0.85 ? moderate=1 else severe=1

level=(ran1<0.13?0:ran1<0.40?1:ran1<0.85?2:3)
# level defined, now we reset the state variables
nodisease=(level==0?1:0)
mild=(level==1?1:0)
moderate=(level==2?1:0)
severe=(level==3?1:0)


}
}
intiquan
Site Admin
Posts: 8
Joined: Mon Sep 03, 2018 11:50 am
Location: Basel, Switzerland
Contact:

Re: continuous markov process simulation

Post by intiquan » Thu Jul 04, 2019 11:00 am

Hi Serge,

again a very good and important question!

To be honest, I am not familiar with Phoenix. Personally I do not like GUIs to conduct modeling and simulation work - especially in this controlled environment where we need to ensure that we can reproduce all our analysis (from data preparation to the final Word report) by executing a single "RUNALL" script from within R.
How can I reset some sate variables to 0?
IQRmodels do have the section "********** MODEL EVENTS". Here you are free to define conditions that will trigger a reset or a new setting of both model state variables and parameters.
Example:

********** MODEL EVENTS
EVENT1 = gt(time,30),X,0

This code above will reset a state X to 0 when time gets larger than 30 time units. The name of the event "EVENT1" is arbitrary but should follow sum standard rules (no space, not special signs, starting with A-Za-z, etc. Think of it to follow the same rules as a variable in R ... without ".".

********** MODEL EVENTS
EVENT1 = gt(A+2,B),X,X+5, Y,X+Y

The event above will assume a state A is present in the model and if (A+2) becomes larger than B (which can be a state or a parameter or a variable) two things are done:
1) the state X is increased by 5
2) the state Y is set to a new value of X+Y

More complex conditions can be used.

How can I redefine the initial conditions of each state variable depending on the simulated outcome?
You can on one hand do as I show above ... but you could also put a small R script around the simulations and suddenly there is nothing that is impossible.

Hope that helps little,

Henning
Post Reply