--- title: "Design evaluation and optimization in discrete space" classoption: openany output: rmarkdown::html_vignette: toc: true bibliography: references.bib biblio-style: apalike link-citations: yes linkcolor: blue urlcolor: green vignette: > %\VignetteIndexEntry{Design evaluation and optimization in discrete space} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r global_options, echo = FALSE, include = FALSE} backup_options <- options() options(width = 3000) knitr::opts_chunk$set(collapse = TRUE, comment = "#>",echo = FALSE, warning = FALSE, message = FALSE, cache = FALSE, tidy = FALSE, size = "small") library(PFIM) ``` # Overview Data for this example result from a PKPD population study on a non-steroidal molecule [@FloresMurrieta1998]. Single doses of 1, 3.2, 10, 31.6, 56.2, or 100mg/kg per os were given respectively to 6 groups, each of which contained at least 6 rats (weighing 200g on average), following a parallel design. Blood sampling and drug response (DI score) evaluation were conducted at 0, 15, 30, and 45 min and at 1, 1.25, 1.5, 2, 3 and 4 hours after administration. Based on the evaluation results, we choose to optimize a design for 30 rats receiving a single dose of either 100mg/kg or 320 mg/kg, selecting only 3 from previously defined time points for blood sampling and response evaluation by using Fedorov-Wynn method and Multiplicative Algorithm. Reports of the design evaluation and optimization are available at https://github.com/iame-researchCenter/PFIM # Design evaluation ### Define the model equations of the PKPD model ```{r, echo = TRUE, comment=''} modelEquations = list( outcomes = list( "RespPK" = "Cc", "RespPD" = "E" ), equations = list( "Deriv_Cc" = "dose_RespPK/V*ka*exp(-ka*t) - Cl/V*Cc", "Deriv_E" = "Rin*(1-Imax*(Cc**gamma)/(Cc**gamma + IC50**gamma))-kout*E" ) ) ``` ### Set mu and omega for each parameter ```{r, echo = TRUE, comment=''} modelParameters = list( ModelParameter( name = "V", distribution = LogNormal( mu = 0.74, omega = 0.316 ) ), ModelParameter( name = "Cl", distribution = LogNormal( mu = 0.28, omega = 0.456 ) ), ModelParameter( name = "ka", distribution = LogNormal( mu = 10, omega = sqrt( 0 ) ), fixedMu = TRUE ), ModelParameter( name = "kout", distribution = LogNormal( mu = 6.14, omega = 0.947 ) ), ModelParameter( name = "Rin", distribution = LogNormal( mu = 614, omega = sqrt( 0 ) ), fixedMu = TRUE ), ModelParameter( name = "Imax", distribution = LogNormal( mu = 0.76, omega = 0.439 ) ), ModelParameter( name = "IC50", distribution = LogNormal( mu = 9.22, omega = 0.452 ) ), ModelParameter( name = "gamma", distribution = LogNormal( mu = 2.77, omega = 1.761 ) ) ) ``` ### Create the error model to the responses PK and PD. ```{r, echo = TRUE, comment=''} errorModelRespPK = Combined1( outcome = "RespPK", sigmaInter = 0, sigmaSlope = 0.21 ) errorModelRespPD = Constant( outcome = "RespPD", sigmaInter = 9.6 ) modelError = list( errorModelRespPK, errorModelRespPD ) ``` ### Define the arms and for each arm + create the administration parameters for the response PK + create the sampling times for the responses PK and PD ```{r, echo = TRUE, comment=''} # sampling times samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4 ) ) samplingTimesRespPD = SamplingTimes( outcome = "RespPD", samplings = c( 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4 ) ) # Define the arms and the administrations administrationRespPK1 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 0.2 ) ) arm1 = Arm( name = "0.2mg Arm", size = 6, administrations = list( administrationRespPK1 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) administrationRespPK2 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 0.64 ) ) arm2 = Arm( name = "0.64mg Arm", size = 6, administrations = list( administrationRespPK2 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) administrationRespPK3 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 2 ) ) arm3 = Arm( name = "2mg Arm", size = 6, administrations = list( administrationRespPK3 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) administrationRespPK4 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 6.24 ) ) arm4 = Arm( name = "6.24mg Arm", size = 6, administrations = list( administrationRespPK4 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) administrationRespPK5 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 11.24 ) ) arm5 = Arm( name = "11.24mg Arm", size = 6, administrations = list( administrationRespPK5 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) administrationRespPK6 = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 20 ) ) arm6 = Arm( name = "20mg Arm", size = 6, administrations = list( administrationRespPK6 ) , samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), initialCondition = list( "Cc" = 0, "E" = 100 ) ) ``` ### Add the arms to the design `design1` ```{r, echo = TRUE, comment=''} design1 = Design( name = "design1", arms = list( arm1, arm2, arm3, arm4, arm5, arm6 ) ) ``` ### Evaluate the population FIM ```{r, echo = TRUE, eval = FALSE, comment=''} evaluationPop = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" = "Cc", "RespPD" = "E" ), designs = list( design1 ), fim = "population", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationPop = run( evaluationPop ) ``` ### Display the results of the design evaluation ```{r, echo = TRUE, eval = FALSE, comment=''} show( evaluationPop ) fisherMatrix = getFisherMatrix( evaluationPop ) getCorrelationMatrix( evaluationPop ) getSE( evaluationPop ) getRSE( evaluationPop ) getShrinkage( evaluationPop ) getDeterminant( evaluationPop ) getDcriterion( evaluationPop ) ``` ### Graphs of the results of the design evaluation ```{r, echo = TRUE, eval = FALSE, comment=''} plotOptions = list( unitTime = c("hour"), unitOutcomes = c("mcg/mL","DI%") ) plotEvaluation = plotEvaluation( evaluationPop, plotOptions ) plotSensitivityIndice = plotSensitivityIndice( evaluationPop, plotOptions ) SE = plotSE( evaluationPop ) RSE = plotRSE( evaluationPop ) # examples plotOutcomesEvaluationRespPK = plotEvaluation[["design1"]][["20mg Arm"]][["RespPK"]] plotOutcomesEvaluationRespPD = plotEvaluation[["design1"]][["20mg Arm"]][["RespPD"]] plotSensitivityIndice_RespPK_Cl = plotSensitivityIndice[["design1"]][["20mg Arm"]][["RespPK"]][["Cl"]] SE = SE[["design1"]] RSE = RSE[["design1"]] ``` ### Report for the design evaluation ```{r, echo = TRUE, eval = FALSE, comment=''} outputPath = "C:/...." outputFile = "Example01_EvaluationPopFIM.html" plotOptions = list( unitTime=c("hour"), unitResponses= c("mcg/mL","DI%") ) Report( evaluationPop, outputPath, outputFile, plotOptions ) ``` ### Evaluate the individual and Bayesian FIMs ```{r, echo = TRUE, eval = FALSE, comment=''} evaluationInd = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" = "Cc", "RespPD" = "E" ), designs = list( design1 ), fim = "individual", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationInd = run( evaluationInd ) evaluationBay = Evaluation( name = "", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outcomes = list( "RespPK" = "Cc", "RespPD" = "E" ), designs = list( design1 ), fim = "Bayesian", odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) evaluationBay = run( evaluationBay ) ``` ### Display the results of the design evaluation ```{r, echo = TRUE, eval = FALSE, comment=''} show( evaluationInd ) show( evaluationBay ) fisherMatrix = getFisherMatrix( evaluationInd ) getCorrelationMatrix( evaluationInd ) getSE( evaluationInd ) getRSE( evaluationInd ) getShrinkage( evaluationInd ) getDeterminant( evaluationInd ) getDcriterion( evaluationInd ) fisherMatrix = getFisherMatrix( evaluationBay ) getCorrelationMatrix( evaluationBay ) getSE( evaluationBay ) getRSE( evaluationBay ) getShrinkage( evaluationBay ) getDeterminant( evaluationBay ) getDcriterion( evaluationBay ) ``` ### Reports of the results of the design evaluation ```{r, echo = TRUE, eval = FALSE, comment=''} outputPath = "C:/...." plotOptions = list( unitTime=c("hour"), unitResponses= c("mcg/mL","DI%") ) outputFile = "Example01_EvaluationIndFIM.html" Report( evaluationInd, outputPath, outputFile, plotOptions ) outputFile = "Example01_EvaluationBayFIM.html" Report( evaluationBay, outputPath, outputFile, plotOptions ) ``` # Design optimization ### Create an administration for response PK and samplings times ```{r, echo = TRUE, comment=''} administrationRespPK = Administration( outcome = "RespPK", timeDose = c(0), dose = c( 6.24 ) ) ``` We create sampling times that will be used in the initial design for comparison during the optimization process. ```{r, echo = TRUE, comment=''} samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 0.25, 0.75, 1, 1.5, 2, 4, 6 ) ) samplingTimesRespPD = SamplingTimes( outcome = "RespPD", samplings = c( 0.25, 0.75, 1.5, 2, 3, 6, 8, 12 ) ) ``` ### Create a sampling constraint for the responses PK and PD + Set the vector of allowed sampling times for the responses PK and PD + Fix certain time values + Set the number of optimizable sampling times ```{r, echo = TRUE, comment=''} samplingConstraintsRespPK = SamplingTimeConstraints( outcome = "RespPK", initialSamplings = c( 0.25, 0.75, 1, 1.5, 2, 4, 6 ), fixedTimes = c( 0.25, 4 ), numberOfsamplingsOptimisable = 4 ) samplingConstraintsRespPD = SamplingTimeConstraints( outcome = "RespPD", initialSamplings = c( 0.25, 0.75, 1.5, 2, 3, 6, 8, 12 ), fixedTimes = c( 2, 6 ), numberOfsamplingsOptimisable = 4 ) ``` ### Set the initial vectors for the FedorovWynn algorithm ```{r, echo = TRUE, comment=''} initialElementaryProtocols = list( c( 0.25, 0.75, 1, 4 ), c( 1.5, 2, 6, 12 ) ) ``` ### Create an administration constraint with the vector of allowed amount of doses for the response PK ```{r, echo = TRUE, comment=''} administrationConstraintsRespK = AdministrationConstraints( outcome = "RespPK", doses = c( 0.2, 0.64, 2, 6.24, 11.24, 20 ) ) ``` ### Create the constraint design to the project + total number of individuals to be considered + administration constraint + sampling constraint ```{r, echo = TRUE, comment=''} armConstraint = Arm( name = "armConstraint", size = 30, administrations = list( administrationRespPK ), samplingTimes = list( samplingTimesRespPK, samplingTimesRespPD ), administrationsConstraints = list( administrationConstraintsRespK ), samplingTimesConstraints = list( samplingConstraintsRespPK, samplingConstraintsRespPD ), initialCondition = list( "Cc" = 0, "E" = "Rin/kout" ) ) designConstraint = Design( name = "designConstraint", arms = list( armConstraint ) ) ``` ### Set the vector of initial proportions or numbers of subjects for each elementary design ```{r, echo = TRUE, comment=''} numberOfSubjects = c(30) proportionsOfSubjects = c(30)/30 ``` ### Set the the parameters of the FedorovWynn algorithm and run the algorithm for the design optimization with a population FIM ```{r, echo = TRUE, eval = FALSE, comment=''} optimizationFWPopFIM = Optimization( name = "PKPD_ODE_multi_doses_populationFIM", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "FedorovWynnAlgorithm", optimizerParameters = list( elementaryProtocols = initialElementaryProtocols, numberOfSubjects = numberOfSubjects, proportionsOfSubjects = proportionsOfSubjects, showProcess = T ), designs = list( designConstraint ), fim = "population", outcomes = list( "RespPK" = "Cc","RespPD" = "E" ), odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) optimizationFWPopFIM = run( optimizationFWPopFIM ) ``` ### Display and plot the results of the design optimization ```{r, echo = TRUE, eval = FALSE, comment=''} show( optimizationFWPopFIM ) fisherMatrix = getFisherMatrix( optimizationFWPopFIM ) getCorrelationMatrix( optimizationFWPopFIM ) getSE( optimizationFWPopFIM ) getRSE( optimizationFWPopFIM ) getShrinkage( optimizationFWPopFIM ) getDeterminant( optimizationFWPopFIM ) getDcriterion( optimizationFWPopFIM ) # plot the frequencies plotFrequencies = plotFrequencies( optimizationFWPopFIM ) plotFrequencies ``` ### Reports for the design optimization ```{r, echo = TRUE, eval = FALSE, comment=''} outputFile = "Example01_OptimizationFWPopFIM.html" Report( optimizationFWPopFIM, outputPath, outputFile, plotOptions ) ``` ### Set the the parameters of the Multiplicative Algorithm algorithm and run the algorithm for the design optimization with a population FIM ```{r, echo = TRUE, eval = FALSE, comment=''} optimizationMultPopFIM = Optimization( name = "PKPD_ODE_multi_doses_populationFIM", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "MultiplicativeAlgorithm", optimizerParameters = list( lambda = 0.99, numberOfIterations = 1000, delta = 1e-04, showProcess = T ), designs = list( designConstraint ), fim = "population", outcomes = list( "RespPK" = "Cc","RespPD" = "E" ), odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) ``` ### Run the Multiplicative Algorithm algorithm for the optimization with a population FIM ```{r, echo = TRUE, eval = FALSE, comment=''} optimizationMultPopFIM = run( optimizationMultPopFIM ) ``` ### Display and plot the results of the design optimization ```{r, echo = TRUE, eval = FALSE, comment=''} show( optimizationMultPopFIM ) fisherMatrix = getFisherMatrix( optimizationMultPopFIM ) getCorrelationMatrix( optimizationMultPopFIM ) getSE( optimizationMultPopFIM ) getRSE( optimizationMultPopFIM ) getShrinkage( optimizationMultPopFIM ) getDeterminant( optimizationMultPopFIM ) getDcriterion( optimizationMultPopFIM ) # plot the weight plotWeights = plotWeights( optimizationMultPopFIM, threshold = 0.001 ) plotWeights ``` ### Create and save the report for the design optimization ```{r, echo = TRUE, eval = FALSE, comment=''} outputPath = "C:/...." plotOptions = list( unitTime = c( "hour" ), unitResponses= c( "mcg/mL" , "DI%" ), threshold = 0.01 ) outputFile = "Example01_OptimizationMultPopFIM.html" Report( optimizationMultPopFIM, outputPath, outputFile, plotOptions ) ``` # References ```{r global_options_end, echo = FALSE, include = FALSE} options(backup_options) ```