Workflows¶

import numpy as np
import pandas as pd

from pkpdutils import NCAOptions, Route, Timecourses, nca, summary_table
from pkpdutils.plot import plot_mean_timecourse

# the study as it arrives: one row per event (a dose or a sample), three dose
# groups of four subjects, the dose group as a column of its own
rng = np.random.default_rng(1)
time = np.array([0.25, 0.5, 1, 2, 3, 4, 6, 8, 12, 24])
ke = rng.uniform(0.15, 0.3, size=4)
ka = rng.uniform(1.0, 3.0, size=4)
records = []
for amount in (50.0, 100.0, 200.0):
    for j in range(4):
        curve = (
            amount
            / 40
            * ka[j]
            / (ka[j] - ke[j])
            * (np.exp(-ke[j] * time) - np.exp(-ka[j] * time))
            * rng.lognormal(0, 0.05, size=time.size)
        )
        group, subject = f"{amount:.0f}", f"{amount:.0f}-s{j + 1}"
        records.append(
            {
                "ID": subject,
                "TIME": 0.0,
                "DV": np.nan,
                "AMT": amount,
                "EVID": 1,
                "dose": group,
            }
        )
        records += [
            {"ID": subject, "TIME": t, "DV": c, "AMT": 0.0, "EVID": 0, "dose": group}
            for t, c in zip(time, curve, strict=True)
        ]
events = pd.DataFrame(records)
events.to_csv("study.csv", index=False)  # the table of the Quickstart
print(events.head())

# the table becomes a batch of twelve curves, each with its dosing protocol;
# a column which is constant within a subject becomes a coordinate
batch = Timecourses.from_events(
    events,
    time_unit="hr",
    unit="mg/l",
    dose_unit="mg",
    route=Route.ORAL,
    substance="drug",
    covariates=["dose"],
)

# every curve is analysed in one vectorized call
result = nca(batch, options=NCAOptions())

# the parameter table of the report: geometric mean and CV per dose group
print(
    summary_table(
        result,
        "individual",
        by="dose",
        parameters=["auc_inf_obs", "cmax", "thalf", "cl_f"],
        stats=("n", "geomean", "geocv"),
    ).to_string(index=False)
)

result.to_dataframe().to_csv("study_parameters.csv", index=False)
plot_mean_timecourse(batch, by="dose").savefig("study_curves.png", dpi=120)

import numpy as np

from pkpdutils import Route, Timecourses, bioequivalence, nca
from pkpdutils.plot import plot_ratio
from pkpdutils.stats import ratio_table

# twelve subjects: sequence RT takes the reference in period 1 and the test in
# period 2, sequence TR the other way round; the test formulation has a lower
# bioavailability (0.93) and a slower absorption
time = np.array([0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 24])
subjects = [f"s{i:02d}" for i in range(12)]
sequence = np.array(["RT"] * 6 + ["TR"] * 6)
period_test = np.where(sequence == "RT", 2, 1)
rng = np.random.default_rng(12)
subject_scale = rng.lognormal(0, 0.25, 12)  # between-subject variability


def period_batch(bioavailability: float, ka: float, period: np.ndarray) -> Timecourses:
    ke = 0.15
    scale = subject_scale * np.where(period == 2, 1.05, 1.0)  # period 2 runs higher
    values = np.stack(
        [
            s
            * bioavailability
            * 100
            * ka
            / (ka - ke)
            * (np.exp(-ke * time) - np.exp(-ka * time))
            / 30
            * rng.lognormal(0, 0.06, time.size)
            for s in scale
        ]
    )
    # the period and the sequence of every subject travel to the result as
    # coordinates along the individual dimension and make the design a 2x2
    return Timecourses.from_arrays(
        time,
        values,
        time_unit="hr",
        unit="mg/l",
        dims=("individual",),
        coords={
            "individual": subjects,
            "period": ("individual", period),
            "sequence": ("individual", sequence),
        },
        dose={"amount": np.full(12, 100.0), "unit": "mg"},
        route=Route.ORAL,
        substance="drug",
    )


reference = nca(period_batch(1.0, 1.5, 3 - period_test))
test = nca(period_batch(0.93, 0.9, period_test))

be = bioequivalence(test, reference, parameters=["auc_inf_obs", "auc_last", "cmax"])
print(ratio_table(be).to_string(index=False))
print("design:", be["cmax"].design, "| bioequivalent:", be.bioequivalent)

plot_ratio(be).savefig("bioequivalence.png", dpi=120)

import numpy as np

from pkpdutils import Route, Timecourses, ddi_classification, nca, ratio
from pkpdutils.plot import plot_ratio
from pkpdutils.stats import DDIThresholds, ddi_table

# the substrate alone and with the perpetrator, two parallel groups of ten
# subjects; the inhibitor lowers the elimination of the substrate to 35 %
time = np.array([0.5, 1, 2, 4, 6, 8, 12, 24, 36, 48])
rng = np.random.default_rng(8)


def arm(clearance_factor: float, label: str) -> Timecourses:
    ke = 0.12 * clearance_factor
    values = np.stack(
        [
            rng.lognormal(np.log(8), 0.2)
            * np.exp(-ke * time)
            * rng.lognormal(0, 0.05, time.size)
            for _ in range(10)
        ]
    )
    return Timecourses.from_arrays(
        time,
        values,
        time_unit="hr",
        unit="mg/l",
        dims=("individual",),
        coords={"individual": [f"{label}{i}" for i in range(10)]},
        dose={"amount": np.full(10, 100.0), "unit": "mg"},
        route=Route.IV_BOLUS,
        substance="substrate",
    )


control = nca(arm(1.0, "control"))
inhibited = nca(arm(0.35, "inhibitor"))

# the exposure ratio with over without the perpetrator, and its class
auc_ratio = ratio(
    inhibited.sample("auc_inf_obs", "individual"),
    control.sample("auc_inf_obs", "individual"),
)
cmax_ratio = ratio(
    inhibited.sample("cmax", "individual"), control.sample("cmax", "individual")
)
ddi = ddi_classification(auc_ratio, cmax_ratio=cmax_ratio)
print(f"{ddi.strength} {ddi.kind}, uncertain: {ddi.uncertain}")

# the same over several parameters at once, one row each
print(
    ddi_table(inhibited, control, ["auc_inf_obs", "cmax"], dim="individual").to_string(
        index=False
    )
)

plot_ratio(
    {"auc_inf_obs": auc_ratio, "cmax": cmax_ratio},
    limits=None,
    thresholds=DDIThresholds.fda(),
).savefig("ddi.png", dpi=120)

import numpy as np

from pkpdutils import AUCMethod, Dose, Dosing, NCAOptions, Route, Timecourses, nca
from pkpdutils.nca import accumulation_ratio
from pkpdutils.plot import plot_intervals, plot_troughs

# four subjects on a twice daily oral regimen over two days, sampled in every
# dosing interval; the protocol is what makes this a multiple dose analysis
dose = Dose(amount=100, unit="mg", time=0, route=Route.ORAL)
protocol = Dosing.regimen(dose, interval=12, n_doses=4)
subjects = ["s1", "s2", "s3", "s4"]
offsets = np.array([0.5, 1, 2, 4, 8, 12])
time = np.concatenate([dose_time + offsets for dose_time in protocol.times])

rng = np.random.default_rng(3)
ke = rng.uniform(0.12, 0.18, size=4)
ka = rng.uniform(0.8, 1.2, size=4)
values = np.empty((4, time.size))
for j in range(4):
    elapsed = time[None, :] - protocol.times[:, None]
    single = (
        100.0
        / 20.0
        * ka[j]
        / (ka[j] - ke[j])
        * (np.exp(-ke[j] * elapsed) - np.exp(-ka[j] * elapsed))
    )
    values[j] = np.where(elapsed >= 0, single, 0.0).sum(axis=0) * rng.lognormal(
        0, 0.03, size=time.size
    )


def batch_of(times: np.ndarray, data: np.ndarray, dosing: Dose | Dosing) -> Timecourses:
    return Timecourses.from_arrays(
        times,
        data,
        time_unit="hr",
        unit="mg/l",
        dims=("individual",),
        coords={"individual": subjects},
        dose=dosing,
        route=Route.ORAL,
        substance="drug",
    )


options = NCAOptions(auc_method=AUCMethod.LOG)
result = nca(batch_of(time, values, protocol), options=options)

# one row per subject and dosing interval
print(
    result.intervals()[["individual", "interval", "interval_auc", "interval_ctrough"]]
    .head(8)
    .to_string(index=False)
)

# the steady state parameters describe the last complete interval
print(
    result.summary_table(
        "individual",
        parameters=["tau", "auc_tau", "cavg", "ctrough", "fluctuation"],
        stats=("n", "mean", "sd", "cv"),
    ).to_string(index=False)
)

# the accumulation, observed within the protocol and against the first interval
# analysed on its own as a single dose curve
first = time <= 12
single = nca(
    batch_of(time[first], values[:, first], dose),
    options=NCAOptions(auc_method=AUCMethod.LOG, tau=12),
)
print(result["accumulation_ratio_obs"].values.round(3))
ratios = accumulation_ratio(result, single)
print(ratios["accumulation_ratio"].values.round(3))
print(ratios["stationarity_ratio"].values.round(3))

plot_troughs(result, x="interval").savefig("troughs.png", dpi=120)
plot_intervals(result, "interval_ctrough").savefig("intervals.png", dpi=120)

import numpy as np

from pkpdutils import Power, Route, Timecourses, fit_table, nca, proportionality_test
from pkpdutils.fit import proportionality_table
from pkpdutils.plot import plot_dose_proportionality

# a dose escalation whose exposure grows slightly faster than the dose
time = np.array([0.5, 1, 2, 4, 6, 8, 12, 24])
doses = np.array([25.0, 50.0, 100.0, 200.0, 400.0])
rng = np.random.default_rng(4)
values = np.stack(
    [
        d**1.15 / 10 * np.exp(-0.25 * time) * rng.lognormal(0, 0.04, time.size)
        for d in doses
    ]
)
batch = Timecourses.from_arrays(
    time,
    values,
    time_unit="hr",
    unit="mg/l",
    dims=("dose",),
    coords={"dose": doses},
    dose={"amount": doses, "unit": "mg"},
    route=Route.IV_BOLUS,
    substance="drug",
)

result = nca(batch)
# the dose coordinate of the result carries no unit, the fit needs one
ds = result.ds.assign_coords(dose=("dose", doses, {"units": "mg"}))
power = fit_table(Power(), ds, "dose", "auc_inf_obs", dim="dose")
test = proportionality_test(power, dose_range=(25.0, 400.0))
print(proportionality_table(test).to_string(index=False))

plot_dose_proportionality(power, test=test).savefig("dose_proportionality.png", dpi=120)