Design evaluation and optimization 01

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/iame-researchCenter/PFIM

Design evaluation

Define the model equations of the PKPD model

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

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

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

## arms

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


# Administrations
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 ) )

# Administrations
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 ) )

# Administrations
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 ) )

# Administrations
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 ) )

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

evaluationFIM = Evaluation( name = "",
                            modelEquations = modelEquations,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outcomes = list( "RespPK"  = "Cc", "RespPD"  = "E" ),
                            designs = list( design1 ),
                            fim = "population",
                            odeSolverParameters = list( atol = 1e-4, rtol = 1e-4 ) )

evaluationFIM = run( evaluationFIM )

Display the results of the design evaluation

show( evaluationFIM )

Create and save the report for the design evaluation


outputFile = "Example01_EvaluationPopFIM.html"
plotOptions = list( unitTime=c("hour"), unitResponses= c("mcg/mL","DI%") )

Report( evaluationFIM, outputPath, outputFile, plotOptions )

Evaluate the individual and Bayesian FIMs

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

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-4, rtol = 1e-4 ) )

evaluationBay = run( evaluationBay )

Display the results of the design evaluation

show( evaluationInd )
show( evaluationBay )

Reports for the design evaluation


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

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

samplingConstraintsRespPK  = SamplingTimeConstraints( outcome = "RespPK",
                                                     initialSamplings = c( 0.25, 0.75, 1, 1.5, 2, 4, 6 ),
                                                     fixedTimes = c( 0.25, 4 ),
                                                     numberOfsamplingsOptimisable = 3 
                                                     )

samplingConstraintsRespPD  = SamplingTimeConstraints( outcome = "RespPD",
                                                     initialSamplings = c( 0.25, 0.75, 1.5, 2, 3, 6, 8, 12 ),
                                                     fixedTimes = c( 2, 6 ),
                                                     numberOfsamplingsOptimisable = 3 )

Set the initial vectors for the FedorovWynn algorithm

initialElementaryProtocols = list( c( 0.25, 1, 4 ), c( 1.5, 2, 6 ) )

Create an administration constraint with the vector of allowed amount of doses for the response PK

administrationConstraintsResp = AdministrationConstraints( outcome = "RespPK",
                                                             doses = c( 0.2, 0.64, 2, 6.24, 11.24, 20 )
                                                           )

Create the constraint design to the project


armConstraint = Arm( name = "armConstraint",
            size = 30,
            administrations = list( administrationRespPK1 ),
            samplingTimes   = list( samplingTimesRespPK, samplingTimesRespPD ),
            administrationsConstraints = list( administrationConstraintsResp ),
            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 ),
                             fim = "population",
                             outcomes = list( "RespPK" = "Cc","RespPD" = "E" ),
                             odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

optimizationFWPopFIM = run( optimizationFWPopFIM )

Display the results of the design optimization

show( optimizationFWPopFIM )

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

optimizationMultPopFIM = run( optimizationMultPopFIM )

Display the results of the design optimization

show( optimizationMultPopFIM )

Create and save the report for the design optimization


plotOptions = list( unitTime=c("hour"), unitResponses= c("mcg/mL","DI%"), threshold = 0.01 )

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.