Design evaluation and optimization in continuous space

Overview

In this example, we simulate an 1-compartment model with linear elimination for IV infusion over 1 hour (inspired by (Sukeishi et al. 2022)). One hundred and fifty (150) subjects receive a 400mg loading dose on the first day, followed by 4 daily doses of 200mg. Blood samples are taken at the end of the \(1^{st}\) infusion (H1), H20, H44, H66 and H120. By evaluating this design, we will then select 4 sampling times on intervals (0,48) and (72,120) for an optimal design using PSO (Particle Swarm Optimization) algorithm. The same optimization problem will also be run with PGBO (Population Based Genetic Optimization) and Simplex algorithm.

Reports of the design evaluation and optimization are available at https://github.com/packagePFIM

Design evaluation

Define the PK model Linear1InfusionSingleDose_ClV from the library of models


modelFromLibrary = list("PKModel" = "Linear1InfusionSingleDose_ClV")

Define mu and omega for each parameter

 
modelParameters = list(
  ModelParameter( name = "V",    distribution = LogNormal( mu = 50, omega = sqrt( .26 ) ) ),
  ModelParameter( name = "Cl",   distribution = LogNormal( mu = 5, omega = sqrt( .34 ) ) ) )

Define the error model to the response PK RespPK


errorModelRespPK = Combined1( output = "RespPK", sigmaInter = 0.5, sigmaSlope = sqrt( 0.15 ) )
modelError = list( errorModelRespPK )

Define the administration parameters of the response PK


administrationRespPK = Administration( outcome = "RespPK",
                                       Tinf = rep( 1, 5 ),
                                       timeDose = seq( 0, 96, 24 ),
                                       dose = c( 400, rep( 200, 4 ) ) )

Define the sampling times for the response PK RespPK

samplingTimesRespPK = SamplingTimes( outcome = "RespPK", 
                                     samplings = c( 1,12,24,44,72,120 ) )

Define an arm called arm1 of size 150 with administration administrationRespPK and samplings samplingTimesRespPK


arm1 = Arm( name = "arm1",
            size = 150,
            administrations = list( administrationRespPK ) ,
            samplingTimes   = list( samplingTimesRespPK ) )

Add the arm arm1 to the design design1


design1 = Design( name = "design1", arms = list( arm1 ) )

Evaluate the population, individual and Bayesian FIMs


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

evaluationFIMPop = run( evaluationPop )

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

evaluationFIMInd = run( evaluationInd )

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

evaluationFIMBay = run( evaluationBay )

Display the results of the design evaluations


show( evaluationFIMPop )
show( evaluationFIMInd )
show( evaluationFIMBay )

fisherMatrix = getFisherMatrix( evaluationFIMPop )
getCorrelationMatrix( evaluationFIMPop )
getSE( evaluationFIMPop )
getRSE( evaluationFIMPop )
getShrinkage( evaluationFIMPop )
getDeterminant( evaluationFIMPop )
getDcriterion( evaluationFIMPop )

plotOptions = list( unitTime = c("hour"), unitOutcomes = c("mcg/mL")  )

plotEvaluation = plotEvaluation( evaluationFIMPop, plotOptions )

plotSensitivityIndices = plotSensitivityIndices( evaluationFIMPop, plotOptions )

SE = plotSE( evaluationFIMPop )
RSE = plotRSE( evaluationFIMPop )

# examples
plotOutcomesEvaluationRespPK = plotEvaluation[["design1"]][["arm1"]][["RespPK"]]
plotSensitivityIndice_RespPK_Cl =  plotSensitivityIndices[["design1"]][["arm1"]][["RespPK"]][["V"]]

Create and save the report for the design evaluation



outputPath = getwd()

plotOptions = list( unitTime=c("hour"), unitOutcomes = c("mcg/mL") )

outputFile = "Example02_EvaluationPopFIM.html"
Report( evaluationFIMPop, outputPath, outputFile, plotOptions )

outputFile = "Example02_EvaluationIndFIM.html"
Report( evaluationFIMInd, outputPath, outputFile, plotOptions )

outputFile = "Example02_EvaluationBayFIM.html"
Report( evaluationFIMBay, outputPath, outputFile, plotOptions )

Design optimization

Create an sampling times

We create sampling times that will be used in the initial design for comparison during the optimization process. As PSO does not optimize the dose regimens, we keep the same administration.

samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 1, 48, 72, 120  ) )

Define design constraints

We define sampling times constraints that aim to select 4 sampling times: 2 within the interval [1,48], and 2 within [72, 120], with at least a delay of 5 between two points.


samplingConstraintsRespPK  = SamplingTimeConstraints( outcome = "RespPK",
                                                      initialSamplings = c( 1, 48, 72, 120 ),
                                                      numberOfTimesByWindows = c( 2, 2 ),
                                                      samplingsWindows = list( c( 1, 48 ), c( 72, 120 ) ),
                                                      minSampling = 5 )

Create the constraint arm and the associated design


arm2 = Arm( name = "arm2",
            size = 150,
            administrations = list( administrationRespPK ) ,
            samplingTimes   = list( samplingTimesRespPK ),
            samplingTimesConstraints = list( samplingConstraintsRespPK ) )

design2 = Design( name = "design2", arms = list( arm2 ), numberOfArms = 150 )

Set the the parameters of the PSO algorithm and run the algorithm for the design optimization with a population FIM


optimization = Optimization( name = "",
                             
                             modelFromLibrary = modelFromLibrary,
                             modelParameters = modelParameters,
                             modelError = modelError,
                             
                             optimizer = "PSOAlgorithm",
                             
                             optimizerParameters = list(
                               maxIteration = 100,
                               populationSize = 10,
                               personalLearningCoefficient = 2.05,
                               globalLearningCoefficient = 2.05,
                               seed = 42,
                               showProcess = T  ),
                             
                             designs = list( design2 ),
                             
                             fimType = "population",
                             
                             outputs = list( "RespPK") )

Run the PSO algorithm for the optimization with a population FIM


optimizationPSO = run( optimization )
saveRDS( optimizationPSO, "optimizationPSO.RDS" )

Display the results of the design optimization


show( optimizationPSO )

fisherMatrix = getFisherMatrix( optimizationPSO )
getCorrelationMatrix( optimizationPSO )
getSE( optimizationPSO )
getRSE( optimizationPSO )
getShrinkage( optimizationPSO )
getDeterminant( optimizationPSO )
getDcriterion( optimizationPSO )

Create and save the report for the design optimization



outputPath = "C:/..."

outputFile = "Example02_OptimizationPSOPopFIM.html"

Report( optimizationPSO, outputPath, outputFile, plotOptions )

Set the the parameters of the PGBO algorithm and run the algorithm for the design optimization with a population FIM


optimizationPGBO = Optimization( name = "",
                                 modelFromLibrary = modelFromLibrary,
                                 modelParameters = modelParameters,
                                 modelError = modelError,
                                 
                                 optimizer = "PGBOAlgorithm",
                                 
                                 optimizerParameters = list(
                                   N = 30,
                                   muteEffect = 0.65,
                                   maxIteration = 1000,
                                   purgeIteration = 200,
                                   seed = 42,
                                   showProcess = TRUE  ),
                                 
                                 designs = list( design2 ),
                                 
                                 fimType = "population",
                                 
                                 outputs = list( "RespPK") )

Run the PGBO algorithm for the optimization with a population FIM


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

fisherMatrix = getFisherMatrix( optimizationPGBO )
getCorrelationMatrix( optimizationPGBO )
getSE( optimizationPGBO )
getRSE( optimizationPGBO )
getShrinkage( optimizationPGBO )
getDeterminant( optimizationPGBO )
getDcriterion( optimizationPGBO )

Display the results of the design optimization


show( optimizationPGBO )

Create and save the report for the design optimization



outputPath = "C:/..."

outputFile = "Example02_OptimizationPGBOPopFIM.html"

Report( optimizationPGBO, outputPath, outputFile, plotOptions )

Set the the parameters of the Simplex algorithm and run the algorithm for the design optimization with a population FIM


optimizationSimplex = Optimization( name = "",
                                   modelFromLibrary = modelFromLibrary,
                                   modelParameters = modelParameters,
                                   modelError = modelError,
                                   
                                   optimizer = "SimplexAlgorithm",
                                   
                                   optimizerParameters = list( pctInitialSimplexBuilding = 10,
                                                               maxIteration = 1000,
                                                               tolerance = 1e-10,
                                                               showProcess = TRUE  ),
                                   
                                   designs = list( design2 ),
                                   
                                   fimType = "population",
                                   
                                   outputs = list( "RespPK") )

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

Display the results of the design optimization


show( optimizationSimplex )

fisherMatrix = getFisherMatrix( optimizationSimplex )
getCorrelationMatrix( optimizationSimplex )
getSE( optimizationSimplex )
getRSE( optimizationSimplex )
getShrinkage( optimizationSimplex )
getDeterminant( optimizationSimplex )
getDcriterion( optimizationSimplex )

Create and save the report for the design optimization



outputPath = "C:/..."

outputFile = "Example02_OptimizationSimplexPopFIM.html"

Report( optimizationSimplex, outputPath, outputFile, plotOptions )

saveRDS( optimizationSimplex, file = paste0( outputPath, "Example02_OptimizationSimplexPopFIM.rds" ) )

References

Sukeishi, Asami, Kotaro Itohara, Atsushi Yonezawa, Yuki Sato, Katsuyuki Matsumura, Yoshiki Katada, Takayuki Nakagawa, et al. 2022. “Population Pharmacokinetic Modeling of GS-441524, the Active Metabolite of Remdesivir, in Japanese COVID-19 Patients with Renal Dysfunction.” CPT: Pharmacometrics & Systems Pharmacology 11 (1): 94–103.