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

	ID	TIME	CONC	AMT	ROUTE	II
	Int64	Float64	Float64	Int64	String	Int64
1	2	0.0	0.0	50	ev	0
2	2	0.25	0.35	0	ev	0
3	2	0.5	0.33	0	ev	0
4	2	0.75	0.36	0	ev	0
5	2	1.0	0.37	0	ev	0
6	2	1.5	0.33	0	ev	0

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:

observations_vs_time
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

	parameters	extrema	geomean	geomeanCV	geostd	mean	numsamples	std
	String	Tuple…	Float64	Float64	Float64	Float64	Int64	Float64
1	cmax	(0.34, 0.77)	0.507222	237.271	1.20349	0.515888	107	0.0956897
2	aucinf_obs	(1.61, 4.17)	2.56954	46.6293	1.19816	2.61065	107	0.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

	id	vz	cl	lambdaz	lambdaz_specf	cmax	mrt	aumc
	Int64	Float64	Float64	Float64	Float64?	Float64	Float64	Float64
1	2	116.994	18.5931	0.137327	0.173287	0.0074	6.30874	16.9642
2	3	90.694	15.7161	0.173287	0.22397	0.0146	4.42943	13.9884
3	10	86.099	17.6125	0.22397	0.20118	0.013	4.99676	14.147
4	14	95.2837	12.9591	0.156595	0.103335	0.01	7.93671	31.0819
5	15	102.629	15.4735	0.20118	0.128702	0.0096	6.80482	22.0626
6	16	99.963	19.6887	0.20118	0.18801	0.01	5.14831	13.0689
7	18	112.03	16.6257	0.20118	0.119439	0.0082	7.23706	21.875
8	19	115.472	28.3778	0.22397	0.22397	0.0082	4.17173	7.35114
9	21	98.8231	18.1609	0.0866434	0.233975	0.0096	5.37424	14.6821
10	23	124.628	21.5963	0.173287	0.173287	0.0094	5.11981	11.8743
11	27	106.415	18.4403	0.173287	0.156595	0.0074	6.82738	18.5657
12	28	101.231	23.5235	0.243239	0.243239	0.01	4.41727	9.31446
13	32	61.2703	16.6357	0.287823	0.287823	0.0146	3.8292	11.4152
14	33	89.8264	20.1961	0.137327	0.299737	0.01	4.48313	10.9189
15	34	112.285	23.4704	0.137327	0.274653	0.01	4.69174	9.88074
16	35	93.1374	24.4123	0.22397	0.22397	0.0096	3.94318	8.0806
17	36	129.094	15.3457	0.137327	0.114536	0.008	8.54074	27.9421
18	38	61.8756	24.7525	0.400036	missing	0.0112	3.09497	6.16167
19	40	57.7691	18.0322	0.312142	missing	0.0154	3.18305	8.78804
20	41	91.6955	17.9874	0.20118	0.173287	0.0102	5.21775	14.4968
21	42	100.35	17.3439	0.137327	0.20118	0.0098	5.72516	16.3471
22	43	96.5568	24.7181	0.22397	0.22397	0.0114	3.85494	7.80594
23	44	132.552	31.0771	0.22397	0.22397	0.0076	4.36673	7.01793
24	46	98.903	17.349	0.137327	0.173287	0.0116	5.63225	16.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)

	id	doseamt	tlag	tmax	cmax	tlast	clast	clast_pred	auclast
	Int64	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	2	50	0.0	1.0	0.37	24.0	0.01	0.00965	2.63
2	3	50	0.0	1.0	0.73	20.0	0.01	0.01	3.12
3	10	50	0.0	0.5	0.65	20.0	0.01	0.00983	2.79
4	14	50	0.0	0.75	0.5	24.0	0.02	0.0212	3.71
5	15	50	0.0	1.0	0.48	24.0	0.01	0.0131	3.16
6	16	50	0.0	0.75	0.5	20.0	0.01	0.00971	2.49
7	18	50	0.0	0.75	0.41	24.0	0.01	0.0133	2.94
8	19	50	0.0	0.75	0.41	16.0	0.01	0.00873	1.72
9	21	50	0.0	0.5	0.48	24.0	0.01	0.00626	2.7
10	23	50	0.0	0.5	0.47	20.0	0.01	0.01	2.26
11	27	50	0.0	0.5	0.37	24.0	0.01	0.01	2.65
12	28	50	0.0	0.5	0.5	16.0	0.01	0.0121	2.08
13	32	50	0.0	0.5	0.73	16.0	0.01	0.0111	2.97
14	33	50	0.0	0.5	0.5	20.0	0.01	0.00634	2.43
15	34	50	0.0	0.25	0.5	20.0	0.01	0.00672	2.08
16	35	50	0.0	0.5	0.48	16.0	0.01	0.00835	2.01
17	36	50	0.0	1.0	0.4	24.0	0.02	0.0227	3.09
18	38	50	0.0	0.75	0.56	12.0	0.01	0.01	2.0
19	40	50	0.0	0.25	0.77	12.0	0.02	0.02	2.71
20	41	50	0.0	0.5	0.51	20.0	0.01	0.0111	2.73
21	42	50	0.0	1.5	0.49	24.0	0.01	0.00743	2.82
22	43	50	0.0	0.5	0.57	16.0	0.01	0.00844	1.98
23	44	50	0.0	0.25	0.38	16.0	0.01	0.00905	1.57
24	46	50	0.0	1.0	0.58	24.0	0.01	0.0073	2.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

	id	auc0_12	auc12_24
	Int64	Float64	Float64
1	2	2.2645	0.329483
2	3	2.89846	0.193123
3	10	2.54614	0.223346
4	14	3.00976	0.666289
5	15	2.6432	0.484746
6	16	2.26543	0.208671
7	18	2.40497	0.501636
8	19	1.63556	0.0777078
9	21	2.42495	0.228671
10	23	2.05621	0.193123
11	27	2.23295	0.387556
12	28	1.96172	0.106562
13	32	2.82771	0.106562
14	33	2.26213	0.132819
15	34	1.93152	0.132819
16	35	1.91536	0.0777078
17	36	2.42019	0.642679
18	38	1.95222	0.02
19	40	2.66782	0.04
20	41	2.45431	0.249428
21	42	2.49977	0.269428
22	43	1.89294	0.0777078
23	44	1.47949	0.0777078
24	46	2.49727	0.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

	ID	TIME	CONC	AMT	ROUTE	II
	Int64	Float64	Float64	Int64	String	Int64
1	2	96.0	0.01	50	ev	24
2	2	96.25	0.43	0	ev	0
3	2	96.5	0.46	0	ev	0
4	2	96.75	0.4	0	ev	0
5	2	97.0	0.37	0	ev	0
6	2	97.5	0.33	0	ev	0
7	2	98.0	0.31	0	ev	0
8	2	98.5	0.27	0	ev	0
9	2	99.0	0.23	0	ev	0
10	2	100.0	0.23	0	ev	0
11	2	102.0	0.16	0	ev	0
12	2	104.0	0.13	0	ev	0
13	2	108.0	0.06	0	ev	0
14	2	112.0	0.04	0	ev	0
15	2	116.0	0.02	0	ev	0
16	2	120.0	0.01	0	ev	0
17	3	96.0	0.0	50	ev	24
18	3	96.25	0.66	0	ev	0
19	3	96.5	0.67	0	ev	0
20	3	96.75	0.62	0	ev	0
21	3	97.0	0.61	0	ev	0
22	3	97.5	0.61	0	ev	0
23	3	98.0	0.52	0	ev	0
24	3	98.5	0.39	0	ev	0
⋮	⋮	⋮	⋮	⋮	⋮	⋮

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)

	parameters	extrema	geomean	geomeanCV	geostd	mean	numsamples
	String	Tuple…	Float64	Float64	Float64	Float64	Int64
1	half_life	(1.91, 5.94)	3.43684	36.2521	1.24592	3.5172	107
2	tmax	(0.25, 1.5)	0.523904	291.54	1.52739	0.570093	107
3	cmax	(0.32, 0.84)	0.51204	235.745	1.20711	0.521215	107
4	auclast	(1.47, 4.11)	2.54092	47.9994	1.21963	2.58972	107
5	auc_tau_obs	(1.49, 4.11)	2.55509	47.5629	1.21527	2.60243	107
6	vz_f_obs	(57.7, 150.0)	95.5629	1.27251	1.21604	97.3551	107
7	cl_f_obs	(12.0, 33.0)	19.2703	6.32373	1.2186	19.6561	107
8	aumcinf_obs	(5.0, 29.7)	12.9381	11.1715	1.44538	13.8078	107
9	tau	(24.0, 24.0)	24.0	4.16667	1.0	24.0	107
10	cavgss	(0.0631, 0.174)	0.108151	1126.78	1.21862	0.110213	107

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