Introduction to NCA

Authors

Vijay Ivaturi

Jose Storopoli

Load the necessary libraries

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

Random.seed!(123)

1 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

1.0.1 Load Data

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

# You may need to navigate to the correct path for NCA_tutorial.csv
data = CSV.read("NCA_tutorial.csv", DataFrame, missingstring = ["NA", ".", ""])
first(data, 5)

5×6 DataFrame

Row	ID	TIME	CONC	AMT	ROUTE	II
	Int64	Float64	Float64	Int64	String3	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

1.0.2 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 = @rsubset data :TIME <= 24

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 a specific style of header for a column then you will not need to map that column separately. Otherwise, columns are 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

1.0.3 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).

observations_vs_time(pop[1])

you can switch the axis to log scale by passing axis = (; yscale = log10) as a keyword argument (or use log instead of log10 for natural log).

observations_vs_time(pop[1]; axis = (; yscale = log10))

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.

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"),
)

1.0.4 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]

1.0.4.1 Create Parameters Summaries

params = [:cmax, :aucinf_obs]

2-element Vector{Symbol}:
 :cmax
 :aucinf_obs

param_summary_sd = summarize(report_sd; parameters = params)

2×9 DataFrame

Row	parameters	numsamples	minimum	maximum	mean	std	geomean	geostd	geomeanCV
	String	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	cmax	107	0.34	0.77	0.515888	0.0956897	0.507222	1.20349	18.6827
2	aucinf_obs	107	1.61	4.17	2.61065	0.461849	2.56954	1.19816	18.2274

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

parameters_dist(report_sd; parameter = :aucinf_obs)

1.0.4.2 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.

Volume of Distribution/F, in this case since the drug is given orally

vz = NCA.vz(pop)
first(vz, 5)

5×2 DataFrame

Row	id	vz
	String	Float64
1	2	116.994
2	3	90.694
3	10	86.099
4	14	95.2837
5	15	102.629

Clearance/F, in this case since the drug is given orally

cl = NCA.cl(pop)
first(cl, 5)

5×2 DataFrame

Row	id	cl
	String	Float64
1	2	18.5931
2	3	15.7161
3	10	17.6125
4	14	12.9591
5	15	15.4735

Terminal Elimination Rate Constant, threshold=3 specifies the max no. of time point used for calculation

lambdaz = NCA.lambdaz(pop; threshold = 3)
first(lambdaz, 5)

5×2 DataFrame

Row	id	lambdaz
	String	Float64
1	2	0.137327
2	3	0.173287
3	10	0.22397
4	14	0.156595
5	15	0.20118

slopetimes in this case specifies the exact time point you want for the calculation

lambdaz_1 = NCA.lambdaz(pop; slopetimes = [8, 12, 16])
first(lambdaz_1, 5)

[ Info: ID 38 errored: lambdaz calculation failed because of non-positive concentration
[ Info: ID 40 errored: lambdaz calculation failed because of non-positive concentration
[ Info: ID 76 errored: lambdaz calculation failed because of non-positive concentration
[ Info: ID 107 errored: lambdaz calculation failed because of non-positive concentration
[ Info: ID 116 errored: lambdaz calculation failed because of non-positive concentration
[ Info: ID 192 errored: lambdaz calculation failed because of non-positive concentration

5×2 DataFrame

Row	id	lambdaz
	String	Float64?
1	2	0.173287
2	3	0.22397
3	10	0.20118
4	14	0.103335
5	15	0.128702

Half-life calculation for 4th individual

thalf = NCA.thalf(pop[4])
first(thalf, 5)

1-element Vector{Float64}:
 5.096481600264849

Dose Normalized Cmax

cmax_d = NCA.cmax(pop, normalize = true)
first(cmax_d, 5)

5×2 DataFrame

Row	id	cmax
	String	Float64
1	2	0.0074
2	3	0.0146
3	10	0.013
4	14	0.01
5	15	0.0096

Mean residence time

mrt = NCA.mrt(pop)
first(mrt, 5)

5×2 DataFrame

Row	id	mrt
	String	Float64
1	2	6.36645
2	3	4.53536
3	10	4.98956
4	14	7.53794
5	15	6.69361

AUMC calculation, using :linlog method

aumc = NCA.aumc(pop; method = :linlog)
first(aumc, 5)

5×2 DataFrame

Row	id	aumc
	String	Float64
1	2	17.1523
2	3	14.3832
3	10	14.1224
4	14	29.1743
5	15	21.6259

Note

# since we have two lambdaz calculations rename one column to merge in a dataframe
rename!(lambdaz_1, Dict(:lambdaz => "lambdaz_specf"))

individual_params =
    innerjoin(vz, cl, lambdaz, lambdaz_1, cmax_d, mrt, aumc; on = [:id], makeunique = true)
first(individual_params, 5)

5×8 DataFrame

Row	id	vz	cl	lambdaz	lambdaz_specf	cmax	mrt	aumc
	String	Float64	Float64	Float64	Float64?	Float64	Float64	Float64
1	2	116.994	18.5931	0.137327	0.173287	0.0074	6.36645	17.1523
2	3	90.694	15.7161	0.173287	0.22397	0.0146	4.53536	14.3832
3	10	86.099	17.6125	0.22397	0.20118	0.013	4.98956	14.1224
4	14	95.2837	12.9591	0.156595	0.103335	0.01	7.53794	29.1743
5	15	102.629	15.4735	0.20118	0.128702	0.0096	6.69361	21.6259

1.0.4.3 Calculate AUC at specific intervals

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

# various other methods are :linear, :linlog
auc0_12 = NCA.auc(pop; interval = (0, 12), method = :linuplogdown)
first(auc0_12, 5)

5×2 DataFrame

Row	id	auc0_12
	String	Float64
1	2	2.2645
2	3	2.89846
3	10	2.54614
4	14	3.00976
5	15	2.6432

auc12_24 = NCA.auc(pop; interval = (12, 24), method = :linuplogdown)
first(auc12_24, 5)

5×2 DataFrame

Row	id	auc12_24
	String	Float64
1	2	0.329483
2	3	0.193123
3	10	0.223346
4	14	0.666289
5	15	0.484746

final = innerjoin(report_sd.reportdf, auc0_12, auc12_24; on = [:id], makeunique = true)
first(final, 5)

5×40 DataFrame

Row	id	dose	tlag	tmax	cmax	tlast	clast	clast_pred	auclast	kel	half_life	aucinf_obs	aucinf_pred	vz_f_obs	cl_f_obs	vz_f_pred	cl_f_pred	n_samples	cmax_dn	auclast_dn	aucinf_dn_obs	auc_extrap_obs	aucinf_dn_pred	auc_extrap_pred	aumclast	aumcinf_obs	aumc_extrap_obs	aumcinf_pred	aumc_extrap_pred	n_samples_kel	rsq_kel	rsq_adj_kel	corr_kel	intercept_kel	kel_t_low	kel_t_high	span	route	auc0_12	auc12_24
	String	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	String	Float64	Float64
1	2	50	0.0	1.0	0.37	24.0	0.01	0.00965	2.63	0.159	4.36	2.69	2.69	117.0	18.6	117.0	18.6	16	0.0074	0.0525	0.0538	2.34	0.0537	2.26	15.2	17.1	11.1	17.1	10.8	11	0.997	0.997	0.998	-0.826	1.5	24.0	5.16	EV	2.2645	0.329483
2	3	50	0.0	1.0	0.73	20.0	0.01	0.01	3.12	0.173	4.0	3.18	3.18	90.7	15.7	90.7	15.7	16	0.0146	0.0625	0.0636	1.81	0.0636	1.81	12.9	14.4	10.3	14.4	10.3	3	1.0	1.0	1.0	-1.14	12.0	20.0	2.0	EV	2.89846	0.193123
3	10	50	0.0	0.5	0.65	20.0	0.01	0.00983	2.79	0.205	3.39	2.84	2.84	86.1	17.6	86.1	17.6	16	0.013	0.0558	0.0568	1.72	0.0568	1.69	12.9	14.2	8.59	14.1	8.46	8	0.994	0.993	0.997	-0.531	2.5	20.0	5.16	EV	2.54614	0.223346
4	14	50	0.0	0.75	0.5	24.0	0.02	0.0212	3.71	0.136	5.1	3.86	3.87	95.3	13.0	95.1	12.9	16	0.01	0.0742	0.0772	3.81	0.0773	4.03	24.5	29.1	15.9	29.4	16.6	10	0.993	0.993	0.997	-0.59	2.0	24.0	4.32	EV	3.00976	0.666289
5	15	50	0.0	1.0	0.48	24.0	0.01	0.0131	3.16	0.151	4.6	3.23	3.25	103.0	15.5	102.0	15.4	16	0.0096	0.0633	0.0646	2.05	0.065	2.68	19.6	21.6	9.39	22.3	12.0	11	0.988	0.986	0.994	-0.714	1.5	24.0	4.89	EV	2.6432	0.484746

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)
first(partial_aucs, 5)

5×3 DataFrame

Row	id	auc0_12	auc12_24
	String	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

1.0.5 Multiple-dose NCA

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

multiple_dose_data = @rsubset data :TIME > 24
first(multiple_dose_data, 5)

5×6 DataFrame

Row	ID	TIME	CONC	AMT	ROUTE	II
	Int64	Float64	Float64	Int64	String3	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

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

NCAPopulation (107 subjects):
  Number of missing observations: 0
  Number of blq observations: 168

report_md = run_nca(
    pop_md;
    sigdigits = 3,
    studyid = "STUDY-002",
    studytitle = "Phase 1 MAD 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",
)

Similar to the Single-NCA case 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;
    parameters = [
        :half_life,
        :tmax,
        :cmax,
        :auclast,
        :auc_tau_obs,
        :vz_f_obs,
        :cl_f_obs,
        :aumcinf_obs,
        :tau,
        :cavgss,
    ],
)

10×9 DataFrame

Row	parameters	numsamples	minimum	maximum	mean	std	geomean	geostd	geomeanCV
	String	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	half_life	107	1.91	5.94	3.5172	0.741081	3.43684	1.24592	22.2562
2	tmax	107	0.25	1.5	0.570093	0.232385	0.523904	1.52739	44.3287
3	cmax	107	0.32	0.84	0.521215	0.100599	0.51204	1.20711	18.9908
4	auclast	107	1.47	4.11	2.58972	0.498504	2.54092	1.21963	20.0518
5	auc_tau_obs	107	1.49	4.11	2.60243	0.492456	2.55509	1.21527	19.6837
6	vz_f_obs	107	57.7	150.0	97.3551	18.6422	95.5629	1.21604	19.7489
7	cl_f_obs	107	12.0	33.0	19.6561	4.05638	19.2703	1.2186	19.9654
8	aumcinf_obs	107	5.0	29.7	13.8078	4.9646	12.9381	1.44538	38.1231
9	tau	107	24.0	24.0	24.0	0.0	24.0	1.0	4.46179e-14
10	cavgss	107	0.0631	0.174	0.110213	0.021162	0.108151	1.21862	19.9671

1.0.6 Documentation

Please refer to the NCA Documentation for any further help.

1.0.7 Save NCA Report

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

# Single-NCA final report
report(report_sd, param_summary_sd)

[ Info: Generating Tables
[ Info: Generating Plots
[ Info: Gathering System Information
[ Info: Ordering Report Sections

"output/Phase 1 SAD of DrugX.pdf"

# Multiple-NCA final report
report(report_md, param_summary_md)

[ Info: Generating Tables
[ Info: Generating Plots
[ Info: Gathering System Information
[ Info: Ordering Report Sections

"output/Phase 1 MAD of DrugX.pdf"

Reuse

CC BY-SA 4.0