Introduction

Allen, R.J., Rieger, T.R., & Musante, C.J. (2016). Efficient Generation and Selection of Virtual Populations in Quantitative Systems Pharmacology Models. CPT: pharmacometrics & systems pharmacology, 5 3, 140-6.

Definitions

Definitions of plausible and virtual patients and populations

Workflow

Overview of MAPK Virtual Population

• Try to re-create oberved (simulated) distribution of tumor sizes for GDC treatment
  • start with mapk model
  • only data is tumor size at 56 days

Plausible Patients

• Parameters which are to be varied in the model are identified and physiologic limits on the bounds of these parameters are established based:
  • intiution
  • data
  • experience
• In the mapk model the following 27 parameters are considered:

##  [1] "BRAFt" "CRAFt" "dmax"  "G13"   "Gspry" "Gdusp" "ke2"   "ke3"  
##  [9] "ke4"   "MEKb"  "MEKt"  "PI3Kb" "PI3Kt" "RASb"  "RASt"  "RTK1t"
## [17] "RTK2t" "wOR"   "wRAS"  "ka2"   "V2"    "ka3"   "q2"    "V3"   
## [25] "V3b"   "ka4"   "V4"

Plausible Patients

• Parameter limits set at 50% smaller than the minimum and 50% greater than the maximum values from the reported virtual population
  • Often have to do this based on physiology, prior knowledge
• Initially, no limits on upper bound of tumor size, lower limit = 0.
  • Can set limits at multiple time point (i.e constrained trajectories)

Plausible Patients

params <- readRDS("../mapk/s10vpop_pk.RDS")
pLower <- params %>% summarise_all(.funs = function(x){min(x)*0.5})
pUpper <- params %>% summarise_all(.funs = function(x){max(x)*1.5})
pLower %<>% gather(key="Name",value="Lower")
pUpper %<>% gather(key="Name",value="Upper")
paramLims <- left_join(pLower,pUpper,by="Name")
paramLims %<>% filter(Name %in% p_names)
kable(paramLims)

Plausible Patients

stateLims <- data.frame(Name = 'TUMOR','Lower'=0,'Upper'=Inf,Time=56)
kable(stateLims)

Plausible Patient Generation

A set of initial random parameters is used as the initial guess for a simulated annealing (SA) optimization algorithm

• Typically: \[f(p) = \sum_{i=1...m}(M_i(p) - d_i)^2\]

• Want to converge to plausible region, \(L\), as defined by the output constraints: \[L = X_1 \times ... \times X_m\] where \(X_i=[l_i,u_i]\) are the lower and upper bounds, respectively for ouput \(i\).

• new objective function, \(g(p)\), flat-bottomed quadratic: \[g(p) = \sum_{i=1...m}max((M_i(p) - c_i)^2-(u_i-c_i)^2,0)\] where \(c_i=\frac{l_i+u_i}{2}\). The result of this is \(g(p) = 0\) if \(M(p) \in L\).

• SA algorithm stops when OFV = 0 (or some threshold)
  • parameter set that meets all parameter and state constraints has been found
  • repeated many times with different initial parameter guesses
    • stochastic nature of SA algorithm leads to unique and independent plausible patient population

Simulation Function

model <- mread("mapk", '../mapk/', soloc = '../mapk')
sim_fcn <- function(parameters=NULL,model,pnames,dosing,ICs,simulate=0){
  loadso(model) # Load model
  # Case with parameters defined by plausible patient algorithm
  if(!is.null(parameters)){
    param_in <- data.frame(Names = pnames,Values=parameters)
    param_in <- spread(param_in,key=Names,value=Values)
  # Create case for testing with default parameters
  }else{
    param_in <- data.frame(ID=1)
  }
  # Bind parameters with initial conditions
  param_in %<>% cbind(ICs)
  # Simulate and extract tumor size at dat 56
  output <- model%>%idata_set(param_in) %>%Req(TUMOR)%>%
    obsonly%>%mrgsim(delta=56,end=56,events=as.ev(dosing))%>%
    filter(time==56)%>%as.data.frame()
  # For VP generation, return a list of steady state and non-steady state outputs
  if(simulate==0){
    return(list(NSS=output))
  # Otherwise return the dataset
  }else{
    return(output)
  }
}

Generate Plausible Patients

control=list(runParallel="parallel",nCores = 4,
                       parallel_libs="mrgsolve") # Setup parallel run
plausiblePatients <- generatePPs(model_fn = sim_fcn, NP=1e3, paramLims=paramLims,stateLims = stateLims,
                                 method='SA',
                                 model_args = model_args,
                                 scoreThreshold = 0)

Virtual Population

• Unlikely that plausible population distribution(s) will match data
  • need to decide which plausible patients to include in virtual population -> resampling

• Consider the selection function \(S:V_{pl} \rightarrow \{0,1\}\) as: \[ S(p_j) = \begin{cases} 1, \text{ if } p_j \in V_{pop}\\ 0, \text{ if otherwise} \end{cases} \]
• From Bayes’ thereom we have: \[ P\left(S(p_j) = 1 | M(p_j)\in R\right) = \frac{\left(M(p_j)\in R|S(p_j) = 1\right)P\left( S(p_j)=1\right)}{P\left( M(p_j)\in R\right)} \] this gives the probability of a plausible patient \(p_j\), given that \(M(p_j)\) (the simulated outcomes) is in the region \(R\) and is also a virtual patient.

• If \(p_j\) is in the virtual population then the probability of \(M(p_j)\) being in the same region \(R\) should be approximately equal to the probability of observing the clinical or empirical data in the region \(R\).
  • if \(d \in D\), Bayes’ rule from above can be re-written as: \[ P\left(S(p_j) = 1|M(p_j)\in R\right) \approx \frac{P(d\in R)P\left(S(p_j)=1\right)}{P\left(M(p_j)\in R\right)} \]

Virtual Population

• For a fixed plausible population \(P\left(S(p_j)=1\right)\) should be a constant

• Reformulate Bayes’ Rule again: \[ P\left(S(p_j)=1|M(p_j)\in R\right) \approx \alpha\frac{P(d\in R)}{P\left(M(p_j)\in R\right)} \]

• Need to compute \(P(d\in R)\) and \(P\left(M(p_j) \in R\right)\)
  • Use kernel density estimate (KDE)
    • univariate for mapk problem, can be multivariate (warning: computationally expensive)

• \(\alpha\) ~ size of the virtual population (or the inclusion rate into the virtual population)
  • Here \(\alpha\) will be optimized using particle swarm optimization (PSO) for speed and robustness
  • a vector of random numbers \(r_1,r_2...r_i...r_N\) is generated from the uniform distribution and if \(r_i < \alpha\frac{P(d\in R)}{P\left(M(p_j)\in R\right)}\) \(p_i\) is selected as a possible member of the virtual population
  • “Score” or OFV: two-sample Kolmogorov-Smirnov Test statistics between \(p_{select}\) and the empirical distribution for each measured outcome in the data (not strictly applicable for multivariate problems)
  • find \(\alpha\) to minimize OFV

• Note: for robustness, rejection-sampling is repeated between 10 and 50 times for each value of \(\alpha\) and the final objective function value is taken as the average of the summed KS scores over all trials.

• Final VP generated by resampling with optimal \(\alpha\)

Virtual Population (Results)

Virtual Population Efficiency

Yield: \(N_{VP}/N_{PP}\) = 21.8%

• How can we increase VPop efficiency? * More closely approximate final distribution with plausible patients -> better bounds

stateLims <- data.frame(Name = 'TUMOR','Lower'=0,'Upper'=2.5,Time=56)
kable(stateLims)

Virtual Population Efficiency

Virtual Population Effiency

Yield = 25.5%

Further Improvements

See Rieger TR, Allen RJ, Bystricky L, Chen Y, Colopy GW, Cui Y, Gonzalez A, Liu Y, White RD, Everett RA, Banks HT, Musante CJ: Improving the generation and selection of virtual populations in quantitative systems pharmacology models. Prog Biophys Mol Biol 139:15–22, 2018

• Alternative optimization algorithms for plausible patient generation
  • nested simulated annealing: augments the SA method by targeting PP generation using the probability density function of the target distribution
  • genetic algorithm: iteratively builds a plausible population according to a fitness function defined by the desired distribution
  • Metropolis-Hastings: importance samping

Further Improvements and Performance