Indirect Response Models

Shamir Kalaria
July 18, 2019

Introduction

This is an introduction to how to excute individual and population level simulations for indirect response models driven by continuous PK. For the following simulation, we have chosen use an oral drug formulation that is intended to

  • inhibit the production of a potential biomarker - irm1

  • inhibit of degredation of a potential biomarker - irm2

  • stimulate the production of a potential biomarker - irm3

  • stimulate the degredation of a potential biomarker - irm4

Getting Started

using Pumas, LinearAlgebra, Plots

Model Code

The following provides specifics to the model code:

  • One compartment pharmacokinetic model was used to drive pharmacodynamics

  • Cp is defined in "derived" and represents the plasma concentration in the central compartment

  • IMAX/EMAX parameter was fixed to 1 for maximal inhibition/simulation

  • Random effects on Ka1, CL, Vc, Kin, Kout, and IC50, EC50 were assumed to follow a log-normal distribution

irm1

irm1 = @model begin
    @param begin
        θ  VectorDomain(6)
        Ω   PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        IC50    = θ[6]*exp(η[5])
        IMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        inhibition = 1-(IMAX*cp/(IC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin*inhibition - Kout*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end

Specified Parameters

The initial estimates of the parameters to simulate from. The fixed effects are provided in the θ vector and the between-subject variability parameteres are provided in the Ω vector as variances. A variance of 0.04 suggests a 20% coefficient of variation.

param = (θ = [
          0.5, # Ka  Absorption rate constant 1 (1/time)
          5, # CL   Clearance (volume/time)
          30, # Vc   Central volume (volume)
          10, # Kin  Response in rate constant (1/time)
          0.02, # Kout Response out rate constant (1/time)
          10, # IC50 Concentration for 50% of max inhibition (mass/volume)
          1, # IMAX Maximum inhibition
          ],
    Ω = Diagonal([0.04,0.04,0.04,0.04,0.04]));

Simulation

For the purpose of this tutorial, a simulation of a single subject receiving a single oral dose is shown below:

  • regimen1 provides the single oral dose at time=0 into the gut compartment (cmt=1)

  • subject1 provides a single subject receiving regimen1

regimen1 = DosageRegimen(1500, time=0,cmt=1)
subject1 = Subject(id=1,evs=regimen1)
Subject
  ID: 1
  Events: 1

The simobs function will call in the model object irm1, subject characteristics, parameters. For this simulation, a 100 hour simulation legnth was selected, with observations at every 0.1 hours. The plot() function can be used to visualize the simulation output.

sim = simobs(irm1,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

You can now try out the other indirect response models by changing the @dynamics.

irm2

irm2 = @model begin
    @param begin
        θ  VectorDomain(6)
        Ω   PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        IC50    = θ[6]*exp(η[5])
        IMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        inhibition = 1-(IMAX*cp/(IC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin - Kout*inhibition*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end
PumasModel
  Parameters: θ, Ω
  Random effects: η
  Covariates: 
  Dynamical variables: Gut, Cent, Resp
  Derived: cp, inhibition, resp
  Observed: cp, inhibition, resp
sim = simobs(irm2,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

irm3

Change IMAX and IC50 to EMAX and EC50 in the @param block and the subsequent equations

irm3 = @model begin
    @param begin
        θ  VectorDomain(6)
        Ω   PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        EC50    = θ[6]*exp(η[5])
        EMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        stimulation = 1 + (EMAX*cp/(EC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin*stimulation - Kout*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end
PumasModel
  Parameters: θ, Ω
  Random effects: η
  Covariates: 
  Dynamical variables: Gut, Cent, Resp
  Derived: cp, stimulation, resp
  Observed: cp, stimulation, resp
sim = simobs(irm3,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

irm3

irm4

irm4 = @model begin
    @param begin
        θ  VectorDomain(6)
        Ω   PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        EC50    = θ[6]*exp(η[5])
        EMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        stimulation = 1 + (EMAX*cp/(EC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin - Kout*stimulation*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end
PumasModel
  Parameters: θ, Ω
  Random effects: η
  Covariates: 
  Dynamical variables: Gut, Cent, Resp
  Derived: cp, stimulation, resp
  Observed: cp, stimulation, resp
sim = simobs(irm4,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])
using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

Appendix

These tutorials are part of the PumasTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using PumasTutorials
PumasTutorials.weave_file("pkpd","indirect_response_models.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\accou\AppData\Local\atom\app-1.38.2\atom.exe"  -a
  JULIA_NUM_THREADS = 4

Package Information:

Status `C:\Users\accou\.julia\environments\v1.1\Project.toml`
[621f4979-c628-5d54-868e-fcf4e3e8185c] AbstractFFTs 0.4.1
[c52e3926-4ff0-5f6e-af25-54175e0327b1] Atom 0.8.8
[f0abef60-9ec0-11e9-27de-db6506a91768] AutoOffload 0.1.0
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[4ece37e6-a012-11e8-38cd-91247efc2c34] Bioequivalence 0.1.0
[336ed68f-0bac-5ca0-87d4-7b16caf5d00b] CSV 0.5.9
[c5f51814-7f29-56b8-a69c-e4d8f6be1fde] CUDAdrv 3.0.1
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.1
[49dc2e85-a5d0-5ad3-a950-438e2897f1b9] Calculus 0.5.0
[7057c7e9-c182-5462-911a-8362d720325c] Cassette 0.2.5
[34da2185-b29b-5c13-b0c7-acf172513d20] Compat 2.1.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.1.0
[667455a9-e2ce-5579-9412-b964f529a492] Cubature 1.4.0
[a93c6f00-e57d-5684-b7b6-d8193f3e46c0] DataFrames 0.18.4
[82cc6244-b520-54b8-b5a6-8a565e85f1d0] DataInterpolations 0.2.0
[31a5f54b-26ea-5ae9-a837-f05ce5417438] Debugger 0.5.0
[bcd4f6db-9728-5f36-b5f7-82caef46ccdb] DelayDiffEq 5.9.1
[2b5f629d-d688-5b77-993f-72d75c75574e] DiffEqBase 5.16.3
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.2.0
[31c91b34-3c75-11e9-0341-95557aab0344] DiffEqBenchmarks 0.1.0
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2+
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.13.0
[01453d9d-ee7c-5054-8395-0335cb756afa] DiffEqDiffTools 0.14.0
[aae7a2af-3d4f-5e19-a356-7da93b79d9d0] DiffEqFlux 0.5.2
[071ae1c0-96b5-11e9-1965-c90190d839ea] DiffEqGPU 0.1.0
[c894b116-72e5-5b58-be3c-e6d8d4ac2b12] DiffEqJump 6.1.1+
[8f2b45d5-b17b-5532-9e92-98ae0077e2e3] DiffEqMachineLearning 0.1.0
[78ddff82-25fc-5f2b-89aa-309469cbf16f] DiffEqMonteCarlo 0.15.1
[77a26b50-5914-5dd7-bc55-306e6241c503] DiffEqNoiseProcess 3.3.1
[9fdde737-9c7f-55bf-ade8-46b3f136cc48] DiffEqOperators 3.5.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.2.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.3.0
[41bf760c-e81c-5289-8e54-58b1f1f8abe2] DiffEqSensitivity 3.3.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.6.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[e30172f5-a6a5-5a46-863b-614d45cd2de4] Documenter 0.23.0
[587475ba-b771-5e3f-ad9e-33799f191a9c] Flux 0.8.3
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3+
[ba82f77b-6841-5d2e-bd9f-4daf811aec27] GPUifyLoops 0.2.5
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.1.0
[bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26] GitHub 5.1.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[42fd0dbc-a981-5370-80f2-aaf504508153] IterativeSolvers 0.8.1
[033835bb-8acc-5ee8-8aae-3f567f8a3819] JLD2 0.1.2
[e5e0dc1b-0480-54bc-9374-aad01c23163d] Juno 0.7.0
[2d691ee1-e668-5016-a719-b2531b85e0f5] LIBLINEAR 0.5.1
[7f56f5a3-f504-529b-bc02-0b1fe5e64312] LSODA 0.4.0
[6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f] Libtask 0.3.0
[c7f686f2-ff18-58e9-bc7b-31028e88f75d] MCMCChains 0.3.10
[33e6dc65-8f57-5167-99aa-e5a354878fb2] MKL 0.0.0
[cc2ba9b6-d476-5e6d-8eaf-a92d5412d41d] MLDataUtils 0.5.0
[eb30cadb-4394-5ae3-aed4-317e484a6458] MLDatasets 0.3.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.5.0
[4886b29c-78c9-11e9-0a6e-41e1f4161f7b] MonteCarloIntegration 0.0.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[872c559c-99b0-510c-b3b7-b6c96a88d5cd] NNlib 0.6.0
[8faf48c0-8b73-11e9-0e63-2155955bfa4d] NeuralNetDiffEq 0.1.0
[09606e27-ecf5-54fc-bb29-004bd9f985bf] ODEInterfaceDiffEq 3.3.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.12.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.2.0
[14b8a8f1-9102-5b29-a752-f990bacb7fe1] PkgTemplates 0.6.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.3
[92933f4c-e287-5a05-a399-4b506db050ca] ProgressMeter 1.0.0
[d7b8c89e-ad89-52e0-b9fd-d0ed321fa021] Pumas 0.1.0
[b7b41870-aa11-11e9-048a-09266ec4a62f] PumasTutorials 0.0.1
[438e738f-606a-5dbb-bf0a-cddfbfd45ab0] PyCall 1.91.2
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[1fd47b50-473d-5c70-9696-f719f8f3bcdc] QuadGK 2.1.0
[612083be-0b0f-5412-89c1-4e7c75506a58] Queryverse 0.3.1
[6f49c342-dc21-5d91-9882-a32aef131414] RCall 0.13.3
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[37e2e3b7-166d-5795-8a7a-e32c996b4267] ReverseDiff 0.3.1
[295af30f-e4ad-537b-8983-00126c2a3abe] Revise 2.1.6
[2b6d1eac-7baa-5078-8adc-e6a3e659f14f] SingleFloats 0.1.3
[47a9eef4-7e08-11e9-0b38-333d64bd3804] SparseDiffTools 0.4.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[4c63d2b9-4356-54db-8cca-17b64c39e42c] StatsFuns 0.8.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[9672c7b4-1e72-59bd-8a11-6ac3964bc41f] SteadyStateDiffEq 1.5.0
[789caeaf-c7a9-5a7d-9973-96adeb23e2a0] StochasticDiffEq 6.6.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[fd094767-a336-5f1f-9728-57cf17d0bbfb] Suppressor 0.1.1
[6fc51010-71bc-11e9-0e15-a3fcc6593c49] Surrogates 0.1.0
[9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c] Tracker 0.2.2
[fce5fe82-541a-59a6-adf8-730c64b5f9a0] Turing 0.6.18
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.16.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.1
[e88e6eb3-aa80-5325-afca-941959d7151f] Zygote 0.3.2