# Design evaluation and optimization in discrete space

# Overview

Data for this example result from a PKPD population study on a non-steroidal molecule (Flores-Murrieta et al. 1998). 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/packagePFIM

# Design evaluation

### Define the model equations of the PKPD model

    modelEquations = 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

    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.

    errorModelRespPK = Combined1( output = "RespPK", sigmaInter = 0, sigmaSlope = 0.21 )

    errorModelRespPD = Constant( output = "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

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

    design1 = Design( name = "design1", arms = list( arm1, arm2, arm3, arm4, arm5, arm6 ) )

### Evaluate the population FIM

    evaluationPop = Evaluation( name = " ",
                                modelEquations = modelEquations,
                                modelParameters = modelParameters,
                                modelError = modelError,
                                outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                                designs = list( design1 ),
                                fimType = "population",
                                odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

    evaluationPop = run( evaluationPop )

### Display the results of the design evaluation

    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

    plotOptions = list( unitTime = c("hour"), unitOutcomes = c("mcg/mL","DI%")  )
    plotEvaluation = plotEvaluation( evaluationPop, plotOptions )
    plotSensitivityIndices = plotSensitivityIndices( evaluationPop, plotOptions )

    SE = plotSE( evaluationPop )
    RSE = plotRSE( evaluationPop )

    # examples
    plotOutcomesEvaluationRespPK = plotEvaluation[["design1"]][["20mg Arm"]][["RespPK"]]
    plotOutcomesEvaluationRespPD = plotEvaluation[["design1"]][["20mg Arm"]][["RespPD"]]
    plotSensitivityIndice_RespPK_Cl = plotSensitivityIndices[["design1"]][["20mg Arm"]][["RespPK"]][["Cl"]]

### Report for the design evaluation

    outputPath = getwd()
    outputFile = "Example01_EvaluationPopFIM.html"
    plotOptions = list( unitTime=c("hour"), unitOutcomes = c("mcg/mL","DI%") )
    Report( evaluationPop, outputPath, outputFile, plotOptions )

### Evaluate the individual and Bayesian FIMs

    evaluationInd = Evaluation( name = " ",
                                modelEquations = modelEquations,
                                modelParameters = modelParameters,
                                modelError = modelError,
                                outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                                designs = list( design1 ),
                                fimType = "individual",
                                odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

    evaluationInd = run( evaluationInd )

    evaluationBay = Evaluation( name = " ",
                                modelEquations = modelEquations,
                                modelParameters = modelParameters,
                                modelError = modelError,
                                outputs = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                                designs = list( design1 ),
                                fimType = "Bayesian",
                                odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

    evaluationBay = run( evaluationBay )

### Display the results of the design evaluation

    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

    outputPath = getwd()

    plotOptions = list( unitTime=c("hour"), unitOutcomes= 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

    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.

    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

    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

    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

    administrationConstraintsRespK = AdministrationConstraints( outcome = "RespPK",
                                                                doses = list( 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

    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 ), numberOfArms = 30 )

### Set the vector of initial proportions or numbers of subjects for each elementary design

    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

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

                                          fimType = "population",

                                          outputs = list( "RespPK" = "Cc","RespPD" = "E" ),

                                          odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

    optimizationFWPopFIM = run( optimizationFWPopFIM )

    saveRDS( optimizationFWPopFIM, "optimizationFWPopFIM.RDS" )

### Display and plot the results of the design optimization

    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

    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

    optimizationMultPopFIM = Optimization( name = "PKPD_ODE_multi_doses_populationFIM",
                                           modelEquations = modelEquations,
                                           modelParameters = modelParameters,
                                           modelError = modelError,

                                           optimizer = "MultiplicativeAlgorithm",

                                           optimizerParameters = list( lambda = 0.99,
                                                                       numberOfIterations = 1000,
                                                                       weightThreshold = 0.01,
                                                                       delta = 1e-04, showProcess = T ),

                                           designs = list( designConstraint ),

                                           fimType = "population",

                                           outputs = 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

    optimizationMultPopFIM = run( optimizationMultPopFIM )
    saveRDS( optimizationMultPopFIM, "optimizationMultPopFIM.RDS" )

### Display and plot the results of the design optimization

    show( optimizationMultPopFIM )
    fisherMatrix = getFisherMatrix( optimizationMultPopFIM )
    getCorrelationMatrix( optimizationMultPopFIM )
    getSE( optimizationMultPopFIM )
    getRSE( optimizationMultPopFIM )
    getShrinkage( optimizationMultPopFIM )
    getDeterminant( optimizationMultPopFIM )
    getDcriterion( optimizationMultPopFIM )

    # plot the weight
    plotWeights = plotWeights( optimizationMultPopFIM )
    plotWeights

### Create and save the report for the design optimization

    outputPath = getwd()
    plotOptions = list( unitTime = c( "hour" ), unitOutcomes= c( "mcg/mL" , "DI%" ) )
    outputFile = "Example01_OptimizationMultPopFIM.html"
    Report( optimizationMultPopFIM, outputPath, outputFile, plotOptions )

# References

Flores-Murrieta, Francisco, Holly Kimko, Dora Flores-Acevedo, Francisco López-Muñoz, William Jusko, Mark Sale, and Gilberto Castañeda-Hernández. 1998. “Pharmacokinetic–Pharmacodynamic Modeling of Tolmetin Antinociceptive Effect in the Rat Using an Indirect Response Model: A Population Approach.” Journal of Pharmacokinetics and Biopharmaceutics 26 (November): 547–57.