R/ModelODEDoseNotInEquations.R

In PFIM: Population Fisher Information Matrix

#' Class "ModelODEDoseNotInEquations"
#' @description ...
#' @name ModelODEDoseNotInEquations-class
#' @aliases ModelODEDoseNotInEquations
#' @docType class
#' @include ModelODE.R
#' @export

ModelODEDoseNotInEquations = setClass("ModelODEDoseNotInEquations", contains = "ModelODE")

# ======================================================================================================
# EvaluateModel
# ======================================================================================================

#' @rdname EvaluateModel
#' @export

setMethod(f = "EvaluateModel",
          signature = "ModelODEDoseNotInEquations",
          definition = function( object, arm )
          {
            # ===============================================
            # model parameters
            # ===============================================

            inputsModel = list()

            parameters = getParameters( object )
            modelParametersNames = lapply( parameters, function(x) getName(x) )
            numberOfParameters = getNumberOfParameters( object )

            # assign parameter values
            for ( parameter in parameters )
            {
              parameterMu = getMu( parameter )
              parameterName = getName( parameter )
              assign( parameterName, parameterMu )
            }

            # ===============================================
            # outcomes and variables and initial conditions
            # ===============================================

            modelOutcomes = getOutcomesForEvaluation( object )
            outcomes = names( modelOutcomes )

            administrations = getAdministrations( arm )
            samplingTimes = getSamplingTimes(  arm  )

            initialConditions = getInitialConditions( arm )
            variables = names( initialConditions )

            # variables with administration
            variablesWithAdministration = unlist( lapply( administrations, function(x) getOutcome(x) ) )

            # variables without administration
            variablesWithNoAdministration = variables[!(variables %in% variablesWithAdministration)]

            # variable with sampling times
            variablesWithSamplingTimes = unlist( lapply( samplingTimes, function(x) getOutcome(x) ) )

            numberOfOutcomes = length( modelOutcomes )

            # ===============================================
            # convert model equations string to expression
            # ===============================================

            modelEquations = list()

            modelEquations = getEquations( object )

            equationNames = names( modelEquations )

            for ( equationName in equationNames )
            {
              modelEquations[[equationName]] = parse( text = modelEquations[[equationName]] )
            }

            # ===============================================
            # sampling times variables With No Administration
            # ===============================================

            # model sampling times = sampling times of all responses
            samplingTimesVariable = list()

            for ( variable in variablesWithSamplingTimes )
            {
              samplingTimesVariable[[variable]] = getSamplings( getSamplingTime( arm, variable ) )
            }

            samplingTimesModel = sort( unique( c( 0, unlist( samplingTimesVariable ) ) ) )
            colnames( samplingTimesModel ) = NULL

            # ===============================================
            # time matrix for ODE
            # ===============================================

            # parameters matrix for ode evaluation
            timeDose = list()
            indexTime = list()

            for ( variable in variablesWithAdministration )
            {
              # time dose, dose and tau
              administration = getAdministration( arm, variable )
              timeDose[[variable]] = getTimeDose( administration )
              dose = getDose( administration )
              tau = getTau( administration )

              # for repeated doses
              if ( tau !=0 )
              {
                maxSamplingTimeVariable = max( unlist( samplingTimesVariable ) )
                n = maxSamplingTimeVariable%/%tau
                inputsModel$dose[[variable]] = rep( dose, n+1 )
                timeDose[[variable]] = sort( unique( seq( 0, maxSamplingTimeVariable, tau ) ) )
              }
              else
              {
                # for multi doses
                inputsModel$dose[[variable]] = dose
                timeDose[[variable]] = sort( unique( c( timeDose[[variable]] ) ) )
              }
            }

            # ===============================================
            # Initial conditions
            # ===============================================

            for ( variable in variablesWithAdministration )
            {
              administration = getAdministration( arm, variable )
              dose = getDose( administration )

              initialConditions[[variable]] = dose[1]
            }

            for ( variable in variablesWithNoAdministration )
            {
              assign( "dose",  inputsModel$dose[[variable]][1])
              initialConditions[[variable]] = eval( parse( text = initialConditions[[variable]]  ) )
            }

            # order the names of the variables
            orderNameVariable = gsub("Deriv_","",names( modelEquations ) )
            orderNameVariable = orderNameVariable[order(factor(orderNameVariable))]
            initialConditions = unlist( initialConditions )
            initialConditions = initialConditions[orderNameVariable]

            # number of variables
            numberOfVariables = length( initialConditions )

            # ===============================================
            # Dose event
            # ===============================================

            doseEvent = list()

            for ( variable in variablesWithAdministration )
            {
              doseEvent[[variable]] = data.frame( var = variable,
                                                  time = timeDose[[variable]],
                                                  value = inputsModel$dose[[variable]],
                                                  method = c("add") )
            }

            doseEvent = do.call( rbind, doseEvent )
            doseEvent = doseEvent[ order( doseEvent$time ), ]

            # ensure that event time in times to avoid warning
            uniqueTimes = cleanEventTimes( doseEvent$time, samplingTimesModel )
            samplingTimesModel = sort( unique( c( samplingTimesModel, uniqueTimes ) ) )

            # ===============================================
            # create and solve ode model
            # ===============================================

            modelEvaluation = list()

            modelODE = function( samplingTimesModel, initialConditions, parameters )
            {
              with(as.list(c(initialConditions)),{

                for ( i in 1:numberOfVariables )
                {
                  modelEvaluation[[i]] = eval( modelEquations[[i]] )
                }

                return( list( c( modelEvaluation ) ) )
              })
            }

            # ===============================================
            # evaluate the model
            # ===============================================

            # parameters atol and rtol for ode solver
            odeSolverParameters = getOdeSolverParameters( object )

            atol = odeSolverParameters$atol
            rtol = odeSolverParameters$rtol

            out = ode( initialConditions,
                       samplingTimesModel,
                       modelODE,
                       atol = atol, rtol = rtol,
                       method = "lsodar",
                       hmax = 0,
                       events = list( data = doseEvent ) )

            # ===================================================
            # scale evaluation model ODE with outcomes evaluation
            # take the sampling times
            # ===================================================

            evaluationOutcomes = list()

            for ( variable in variables )
            {
              assign( variable, out[,variable] )
            }

            iterVariableName = 1

            for ( outcome in outcomes )
            {
              variable = variablesWithSamplingTimes[iterVariableName]

              evaluationOutcomes[[outcome]] = eval( parse( text = modelOutcomes[[outcome]]) )

              indexSamplingTimes = which( samplingTimesModel %in%  samplingTimesVariable[[variable]] )

              evaluationOutcomes[[outcome]] = as.data.frame( cbind( samplingTimesModel[indexSamplingTimes],
                                                                    evaluationOutcomes[[outcome]][indexSamplingTimes] ) )

              colnames( evaluationOutcomes[[outcome]] ) = c( "time", outcomes[iterVariableName] )

              iterVariableName = iterVariableName+1
            }

            # =================================================
            # substitute for outcomes evaluation with scaling
            # =================================================

            subsituteTmp = list()

            modelEquationsTmp = getEquations( object )
            modelEquationsNames = names( modelEquations )
            modelEquationsTmpNames = str_remove(names(modelEquationsTmp), "Deriv_")
            names( modelEquationsTmp ) = modelEquationsTmpNames

            for( outcome in outcomes )
            {
              variableOutcome = modelOutcomes[[outcome]]

              for ( variable in variablesWithSamplingTimes )
              {
                modelEquationsTmp[[variable]] = paste0("(", modelEquationsTmp[[variable]],")")
                variableOutcome = gsub( variable, modelEquationsTmp[[variable]], variableOutcome )
              }
              subsituteTmp[[outcome]] = parse( text = variableOutcome )
            }

            # rename equations
            names( subsituteTmp ) = paste0( "Deriv_", c( variablesWithSamplingTimes ) )

            for ( name in  names( subsituteTmp ) )
            {
              modelEquations[[name]] = subsituteTmp[[name]]
            }

            # ===============================================
            # compute gradient
            # ===============================================

            parameters = getParameters( object )

            parametersGradient = parametersForComputingGradient( object )

            shiftedParameters = parametersGradient$shifted
            Xcols = parametersGradient$Xcols
            frac = parametersGradient$frac

            outcomesGradient = list()

            for( variable in variablesWithSamplingTimes )
            {
              resultsGrad = list()

              for ( iterShifted in 1:dim( shiftedParameters)[2] )
              {
                valuesParameters = shiftedParameters[1:numberOfParameters,iterShifted]

                # assign parameter values
                for( iterParameter in 1:numberOfParameters )
                {
                  parameterMu = valuesParameters[iterParameter]
                  parameterName = getName( parameters[[iterParameter]] )
                  assign( parameterName, parameterMu )
                }

                out = ode( initialConditions,
                           samplingTimesModel,
                           modelODE,
                           atol = atol, rtol = rtol,
                           method = "lsoda",
                           events = list( data = doseEvent ) )

                resultsGrad[[iterShifted]] = out[,variable]
              }

              resultsGrad = do.call( cbind, resultsGrad )

              coefs = list()

              for ( i in 1 :dim( resultsGrad )[1] )
              {
                coefs[[i]] = solve( do.call("cbind", Xcols), resultsGrad[i,])/frac
                coefs[[i]] = coefs[[i]][1 + seq_along( parameters )]
              }

              outcomesGradient[[variable]] = do.call( rbind, coefs )

              indexSamplingTimes = which( samplingTimesModel %in%  samplingTimesVariable[[variable]] )
              outcomesGradient[[variable]] = as.data.frame( cbind( samplingTimesModel[indexSamplingTimes], outcomesGradient[[variable]][indexSamplingTimes,] ) )
              colnames(  outcomesGradient[[variable]] ) = c( "time", modelParametersNames )
            }

            names( outcomesGradient ) = outcomes

            # -----------------------------------------------
            # outcomesAllGradient
            # select with model error
            # -----------------------------------------------

            outcomesAllGradient = list()

            modelError = getModelError( object )

            for( outcome in outcomes )
            {
              index = which( sapply( modelError, function (x) getOutcome(x) == outcome ) )

              if ( length( index ) != 0 )
              {
                outcomesAllGradient[[outcome]] = outcomesGradient[[outcome]]
              }
            }

            outcomesAllGradient = as.data.frame( do.call( rbind, outcomesAllGradient ) )
            rownames( outcomesAllGradient ) = NULL

            return( list( evaluationOutcomes = evaluationOutcomes,
                          outcomesGradient = outcomesGradient,
                          outcomesAllGradient = outcomesAllGradient ) )
          })

# ======================================================================================================
# definePKModel
# ======================================================================================================

#' @rdname definePKModel
#' @export
#'
setMethod("definePKModel",
          signature("ModelODE"),
          function( object, outcomes )
          {
            # -------------------------------------------
            # change names: responses, variables, doses
            # -------------------------------------------

            # original and new outcomes
            newOutcomes = outcomes
            originalOutcomes = getOutcomes( object )

            if ( length( outcomes ) != 0 )
            {
              # variable names
              variablesNames = unlist( originalOutcomes )
              variablesNewNames = unlist( newOutcomes )

              # change equation names
              equations = getEquations( object )
              names( equations ) = paste0( "Deriv_", variablesNewNames )

              # response names old and new
              responsesNames = names( originalOutcomes )
              responsesNewNames = names( newOutcomes )

              for ( iterEquation in 1:length( equations ) )
              {
                # change response names
                for( iterResponseName in 1:length( responsesNames ) )
                {
                  equations[[iterEquation]] = gsub( responsesNames[iterResponseName],
                                                    responsesNewNames[iterResponseName], equations[[iterEquation]] )
                }

                # change variable names
                for( iterVariableName in 1:length( variablesNewNames ) )
                {
                  equations[[iterEquation]] = gsub( variablesNames[iterVariableName],
                                                    variablesNewNames[iterVariableName], equations[[iterEquation]] )
                }

              }
              object = setEquations( object, equations )
              object = setOutcomes( object, newOutcomes )

            }else{

              # change only dose name
              equations = getEquations( object )
              responseNames = names( originalOutcomes )

              # set equation and outcome
              object = setOutcomes( object, originalOutcomes )
              object = setEquations( object, equations )
            }

            return( object )
          })

# ======================================================================================================
# definePKPDModel
# ======================================================================================================

#' @rdname definePKPDModel
#' @export
#'
setMethod("definePKPDModel",
          signature("ModelODEDoseNotInEquations","ModelODE"),
          function( PKModel, PDModel, outcomes )
          {
            model = ModelODEDoseNotInEquations()

            if ( length( outcomes ) != 0 )
            {
              # original and new outcomes
              newOutcomes = outcomes
              originalOutcomesPKModel = getOutcomes( PKModel )
              originalOutcomesPDModel = getOutcomes( PDModel )

              originalOutcomesPKModel = unlist( originalOutcomesPKModel )
              originalOutcomesPDModel = unlist( originalOutcomesPDModel )

              originalOutcomes = as.list( c( originalOutcomesPKModel, originalOutcomesPDModel ) )

              # variable names
              variablesNames = unlist( originalOutcomes )
              variablesNewNames = unlist( newOutcomes )

              # model equation
              PKModelEquations = getEquations( PKModel )
              PDModelEquations = getEquations( PDModel )
              equations = c( PKModelEquations, PDModelEquations )
              equations = lapply( equations, function(x) parse( text = x ) )
              names( equations ) = paste0( "Deriv_", variablesNewNames )
              numberOfEquations = length( equations )

              # response names old and new
              responsesNames = names( originalOutcomes )
              responsesNewNames = names( newOutcomes )

              # variables substitution
              variablesNewNames = lapply( variablesNewNames, function(x) parse( text = x ) )
              variablesNewNames = lapply( variablesNewNames, function(x) x[[1]] )

              # RespPK change for PD Model with PK ode Michaelis-Menten
              variablesNewNames = append( variablesNewNames, variablesNewNames[[1]] )
              names( variablesNewNames ) = c( variablesNames, "RespPK" )

              for ( iterEquation in 1:numberOfEquations )
              {
                equations[[iterEquation]] = as.expression(do.call( 'substitute',
                                                                   list( equations[[iterEquation]][[1]], variablesNewNames ) ) )
              }

              # convert equations from expression to string
              equations = lapply( equations, function(x) x[[1]] )
              equations = lapply( equations, function(x) paste( deparse( x ), collapse = " " ) )
              equations = lapply( equations, function(x) str_replace_all( x, " ", "" ) )

              # set outcomes and equations
              model = setEquations( model, equations )
              model = setOutcomes( model, newOutcomes )

            }else{

              # outcomes
              newOutcomes = outcomes
              originalOutcomesPKModel = getOutcomes( PKModel )
              originalOutcomesPDModel = getOutcomes( PDModel )

              originalOutcomesPKModel = unlist( originalOutcomesPKModel )
              originalOutcomesPDModel = unlist( originalOutcomesPDModel )

              originalOutcomes = as.list( c( originalOutcomesPKModel, originalOutcomesPDModel ) )

              # response names old and new
              responsesNames = names( originalOutcomes )

              # variable names
              variablesNames = unlist( originalOutcomes )
              variablesNames = lapply( variablesNames, function(x) parse( text = x ) )

              # model equations
              PKModelEquations = getEquations( PKModel )
              PDModelEquations = getEquations( PDModel )
              equations = c( PKModelEquations, PDModelEquations )
              equations = lapply( equations, function(x) parse( text = x ) )
              numberOfEquations = length( equations )

              # RespPK change for PD Model with PK ode Michaelis-Menten
              variableSubstitutedMichaelisMenten = list()
              variableSubstitutedMichaelisMenten[[1]] = variablesNames[[1]][[1]]
              names( variableSubstitutedMichaelisMenten ) = "RespPK"

              for ( iterEquation in 1:numberOfEquations )
              {
                equations[[iterEquation]] = as.expression(do.call( 'substitute', list( equations[[iterEquation]][[1]],
                                                                                       variableSubstitutedMichaelisMenten ) ) )
              }

              # convert equations from expression to string
              equations = lapply( equations, function(x) x[[1]] )
              equations = lapply( equations, function(x) paste( deparse( x ), collapse = " " ) )
              equations = lapply( equations, function(x) str_replace_all( x, " ", "" ) )

              model = setEquations( model, equations )
              model = setOutcomes( model, originalOutcomes )
            }

            return( model )
          })

##########################################################################################################
# END Class ModelODEDoseNotInEquations
##########################################################################################################