The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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


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


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 )

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


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

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


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

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

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 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, threshold = 0.001 )
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,
                                                                   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 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, threshold = 0.001 )
plotWeights

Create and save the report for the design optimization


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

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.