Non Compartmental Analysis

2021-08-21

Load the necessary libraries

using Random
using DataFrames
using CSV
using NCA
using NCAUtilities
using NCA.Unitful
using Dates

Introduction

This tutorial will introduce the workflow of performing Non-Compartmental Analysis (NCA) in a command-line based workflow. A graphical user interface for NCA is also available for a more interactive experience. Please contact sales@pumas.ai for more information.

For the purpose of this tutorial, we will be using the following dataset:

Blood samples for pharmacokinetics were collected in 107 subjects who received 50 mg of DrugX q24 hours for 5 days, on Day 1 (first dose) and Day 5 (last dose at steady-state). We will compute Single Dose and Multiple Dose NCA parameters

In the given dataset

  • Time = hrs

  • Conc = mg/L

  • Amt = mg

Load Data

Load the NCA dataset using the CSV.read function and convert it to a DataFrame.

data = CSV.read("NCA_tutorial.csv", DataFrame, missingstrings=["NA", ".", ""])
first(data, 6)

6 rows × 6 columns

IDTIMECONCAMTROUTEII
Int64Float64Float64Int64StringInt64
120.00.050ev0
220.250.350ev0
320.50.330ev0
420.750.360ev0
521.00.370ev0
621.50.330ev0

Single-dose NCA

Samples on Day 1 will be used for single dose NCA. So, we filter data for time <= 24 hrs

single_dose_data  = filter(x -> x.TIME <= 24 , data)

We will specify the units of time, concentration and amount which will be used in the analysis.

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg";

Next, we map the columns in the dataset to the required arguments in the read_nca function

  • id

  • observations

  • time

  • ii - if you have multiple doses

  • amt

  • route - Infusion "inf", IV "iv", Oral "ev"

  • duration - Only for Infusion

  • group - any other grouping variables that you require in the output

Note:

  • If your dataset has the specific style of headers for the column you will not need to map them separately. They will be

automatically determined by Pumas. eg: id, time, observations, amt, route, duration(if needed), group(if needed)

  • Do not specify the Inter Dose Interval in this case because we are computing Single-Dose NCA parameters

pop        = read_nca(  single_dose_data,
                        id    = :ID,
                        time  = :TIME,
                        observations  = :CONC,
                        amt   = :AMT,
                        route = :ROUTE,
                        llq   = 0.001
                      )
NCAPopulation (107 subjects):
Number of missing observations: 0
Number of blq observations: 100

Exploratory Analysis

We start by different ways of visualizing the concentration time profiles. Depending on how you are performing your analysis, plots can be either interactive or static. A detailed documentation on NCA related plotting is available here. There are two main plots that can be used:

  1. observations_vs_time

  2. summary_observations_vs_time

If you are working in VSCode, you could enjoy the interactive experience by calling interactive. Let's look at the observations vs time plots for one subject (you can toggle through the subject index).

obsvstimes = observations_vs_time(pop[1])

you can switch the axis to log scale by passing axis = (yscale = log,) as a keyword argument.

obsvstimes = observations_vs_time(pop[1], axis = (yscale = log,))

However, that can become cumbersome for many subjects. Here we show how to view many subjects at once in the non-interactive mode. Let's plot the concentration vs time for a few subject picked randomly.

Random.seed!(123)
rand_subjs = rand(1:length(pop), 9)
ctplots = observations_vs_time(pop[rand_subjs], 
                                axis = (xlabel = "Time hr", 
                                        ylabel = "DrugX Concentration mg/L"),
                                paginate = true, 
                                columns = 3, rows = 3,
                                facet = (combinelabels = true,))
ctplots[1]

You can also visualize the mean concentration vs time plot by using

summary_observations_vs_time(pop,
                        axis = (xlabel = "Time  hr", 
                                ylabel = "DrugX Concentration mg/L"))

NCA analysis

The workhorse function that performs the analysis is run_nca (previously called NCAReport). While it is sufficient to use run_nca(pop), there are a lot of arguments that the function accepts. These are mostly metadata information that is used to generate an automated report that will be covered later.

report_sd    = run_nca(pop, sigdigits=3,
                studyid="STUDY-001",
                studytitle="Phase 1 SAD of DrugX", # required
                author = [("Vijay", "Pumas-AI")], # required
                sponsor = "PumasAI",
                date=Dates.now(),
                conclabel="DrugX Concentration (mg/L)",
                timelabel="Time hr",
                versionnumber=v"0.1",)

You can view the individual subject fits in an interactive way, but for the tutorial we will showcase the non-interactive way

individual_fits = subject_fits(report_sd,
             axis = (xlabel = "Time hr", 
                     ylabel = "DrugX Concentration mg/L",
                     yscale = log10),
             separate = true, paginate = true,
             limit = 16, columns = 4, rows = 4, 
             facet = (combinelabels = true,))

The above code generates a vector of 16 plots (4 x 4). To view each set of plots, you can access via the vector index of 1:16 e.g.

individual_fits[1]
Create Parameters Summaries
parms = [:cmax, :aucinf_obs]
param_summary_sd = summarize(report_sd.reportdf, parameters=parms)

2 rows × 8 columns

parametersextremageomeangeomeanCVgeostdmeannumsamplesstd
StringTuple…Float64Float64Float64Float64Int64Float64
1cmax(0.34, 0.77)0.507222237.2711.203490.5158881070.0956897
2aucinf_obs(1.61, 4.17)2.5695446.62931.198162.610651070.461849

One can use some of the plotting functions to visualize the distribution of the parameters.

parameters_dist(report_sd, parameter = :aucinf_obs)
Compute Specific Parameters

We will now compute the specific parameters that we require and merge all of them to a dataframe. Note how you have to qualify each function with NCA. to avoid variables named as those parameters.

vz        = NCA.vz(pop, sigdigits=3)  # Volume of Distribution/F, in this case since the drug is given orally
cl        = NCA.cl(pop, sigdigits=3)  # Clearance/F, in this case since the drug is given orally
lambdaz   = NCA.lambdaz(pop, threshold=3, sigdigits=3)  # Terminal Elimination Rate Constant, threshold=3 specifies the max no. of time point used for calculation
lambdaz_1 = NCA.lambdaz(pop, slopetimes=[8,12,16], sigdigits=3) # slopetimes in this case specifies the exact time point you want for the calculation
thalf     = NCA.thalf(pop[4], sigdigits=3) # Half-life calculation for 4th individual
cmax_d    = NCA.cmax(pop, normalize=true, sigdigits=3) # Dose Normalized Cmax
mrt       = NCA.mrt(pop, sigdigits=3) # Mean residence time
aumc      = NCA.aumc(pop, method=:linlog, sigdigits=3) # AUMC calculation, using :linlog method
rename!(lambdaz_1, Dict(:lambdaz => "lambdaz_specf")) #since we have two lambdaz calculation rename one column to merge in a dataframe
individual_params      = innerjoin(vz,cl,lambdaz, lambdaz_1,cmax_d,mrt,aumc, on=[:id], makeunique=true)

107 rows × 8 columns

idvzcllambdazlambdaz_specfcmaxmrtaumc
Int64Float64Float64Float64Float64?Float64Float64Float64
12116.99418.59310.1373270.1732870.00746.3087416.9642
2390.69415.71610.1732870.223970.01464.4294313.9884
31086.09917.61250.223970.201180.0134.9967614.147
41495.283712.95910.1565950.1033350.017.9367131.0819
515102.62915.47350.201180.1287020.00966.8048222.0626
61699.96319.68870.201180.188010.015.1483113.0689
718112.0316.62570.201180.1194390.00827.2370621.875
819115.47228.37780.223970.223970.00824.171737.35114
92198.823118.16090.08664340.2339750.00965.3742414.6821
1023124.62821.59630.1732870.1732870.00945.1198111.8743
1127106.41518.44030.1732870.1565950.00746.8273818.5657
1228101.23123.52350.2432390.2432390.014.417279.31446
133261.270316.63570.2878230.2878230.01463.829211.4152
143389.826420.19610.1373270.2997370.014.4831310.9189
1534112.28523.47040.1373270.2746530.014.691749.88074
163593.137424.41230.223970.223970.00963.943188.0806
1736129.09415.34570.1373270.1145360.0088.5407427.9421
183861.875624.75250.400036missing0.01123.094976.16167
194057.769118.03220.312142missing0.01543.183058.78804
204191.695517.98740.201180.1732870.01025.2177514.4968
2142100.3517.34390.1373270.201180.00985.7251616.3471
224396.556824.71810.223970.223970.01143.854947.80594
2344132.55231.07710.223970.223970.00764.366737.01793
244698.90317.3490.1373270.1732870.01165.6322516.2119
Calculate AUC at specific intervals

Calculation of AUC at specific time intervals and merge to the final report dataframe

auc0_12   = NCA.auc(pop, interval=(0,12), method=:linuplogdown, sigdigits=3) #various other methods are :linear, :linlog
auc12_24  = NCA.auc(pop, interval=(12,24), method=:linuplogdown, sigdigits=3)
final     = innerjoin(report_sd.reportdf, auc0_12, auc12_24, on = [:id], makeunique=true)

107 rows × 41 columns (omitted printing of 32 columns)

iddoseamttlagtmaxcmaxtlastclastclast_predauclast
Int64Int64Float64Float64Float64Float64Float64Float64Float64
12500.01.00.3724.00.010.009652.63
23500.01.00.7320.00.010.013.12
310500.00.50.6520.00.010.009832.79
414500.00.750.524.00.020.02123.71
515500.01.00.4824.00.010.01313.16
616500.00.750.520.00.010.009712.49
718500.00.750.4124.00.010.01332.94
819500.00.750.4116.00.010.008731.72
921500.00.50.4824.00.010.006262.7
1023500.00.50.4720.00.010.012.26
1127500.00.50.3724.00.010.012.65
1228500.00.50.516.00.010.01212.08
1332500.00.50.7316.00.010.01112.97
1433500.00.50.520.00.010.006342.43
1534500.00.250.520.00.010.006722.08
1635500.00.50.4816.00.010.008352.01
1736500.01.00.424.00.020.02273.09
1838500.00.750.5612.00.010.012.0
1940500.00.250.7712.00.020.022.71
2041500.00.50.5120.00.010.01112.73
2142500.01.50.4924.00.010.007432.82
2243500.00.50.5716.00.010.008441.98
2344500.00.250.3816.00.010.009051.57
2446500.01.00.5824.00.010.00732.82

One can also pass in a vector of partial AUC's directly as seen below

partial_aucs = NCA.auc(pop, interval = [(0,12), (12,24)],method=:linuplogdown, sigdigits=3)

107 rows × 3 columns

idauc0_12auc12_24
Int64Float64Float64
122.26450.329483
232.898460.193123
3102.546140.223346
4143.009760.666289
5152.64320.484746
6162.265430.208671
7182.404970.501636
8191.635560.0777078
9212.424950.228671
10232.056210.193123
11272.232950.387556
12281.961720.106562
13322.827710.106562
14332.262130.132819
15341.931520.132819
16351.915360.0777078
17362.420190.642679
18381.952220.02
19402.667820.04
20412.454310.249428
21422.499770.269428
22431.892940.0777078
23441.479490.0777078
24462.497270.285943

Multiple-dose NCA

Filter data for time >= 24 hrs to perform the Multiple-dose NCA

multiple_dose_data           = filter(x -> x.TIME > 24 , data)

1,712 rows × 6 columns

IDTIMECONCAMTROUTEII
Int64Float64Float64Int64StringInt64
1296.00.0150ev24
2296.250.430ev0
3296.50.460ev0
4296.750.40ev0
5297.00.370ev0
6297.50.330ev0
7298.00.310ev0
8298.50.270ev0
9299.00.230ev0
102100.00.230ev0
112102.00.160ev0
122104.00.130ev0
132108.00.060ev0
142112.00.040ev0
152116.00.020ev0
162120.00.010ev0
17396.00.050ev24
18396.250.660ev0
19396.50.670ev0
20396.750.620ev0
21397.00.610ev0
22397.50.610ev0
23398.00.520ev0
24398.50.390ev0

In the case of multiple-dose NCA the extra parameters which are calculated based on the dosing interval (τ) will be included in the output. You can select the NCA parameters that you wish to see in the output with the parameters keyword argument to run_nca

pop_md       = read_nca( multiple_dose_data,
                       id    = :ID,
                       time  = :TIME,
                       observations  = :CONC,
                       amt   = :AMT,
                       route = :ROUTE,
                       ii    = :II,        # please specify II for Multiple-dose NCA
                       llq   = 0.001)
report_md    = run_nca(pop_md, sigdigits=3)

Similar to the Single-NCA you can compute individual parameters and merge them to the final report

Compute the summary statistics of the selected variables that you require

param_summary_md  = summarize(report_md.reportdf, 
                            parameters = [:half_life, 
                                          :tmax, :cmax, 
                                          :auclast, :auc_tau_obs, 
                                          :vz_f_obs, :cl_f_obs, 
                                          :aumcinf_obs, :tau, :cavgss])

10 rows × 8 columns (omitted printing of 1 columns)

parametersextremageomeangeomeanCVgeostdmeannumsamples
StringTuple…Float64Float64Float64Float64Int64
1half_life(1.91, 5.94)3.4368436.25211.245923.5172107
2tmax(0.25, 1.5)0.523904291.541.527390.570093107
3cmax(0.32, 0.84)0.51204235.7451.207110.521215107
4auclast(1.47, 4.11)2.5409247.99941.219632.58972107
5auc_tau_obs(1.49, 4.11)2.5550947.56291.215272.60243107
6vz_f_obs(57.7, 150.0)95.56291.272511.2160497.3551107
7cl_f_obs(12.0, 33.0)19.27036.323731.218619.6561107
8aumcinf_obs(5.0, 29.7)12.938111.17151.4453813.8078107
9tau(24.0, 24.0)24.04.166671.024.0107
10cavgss(0.0631, 0.174)0.1081511126.781.218620.110213107

Documentation

Please refer to the NCA Documentation for any further help.

Save NCA Report

The NCA report can be saved generated with report function on the output of run_nca, this generates a pdf report as discussed here.

report(report_sd, param_summary_sd)      # Single-NCA final report
report(report_md, param_summary_md)  # Multiple-NCA final report