nvqsp

name

nvqsp

description

Build and run GPU-accelerated Quantitative Systems Pharmacology (QSP) and PBPK population ODE simulations with the nvQSP library (RODAS4 stiff solver). Use this whenever the user wants to simulate compartmental PK/PBPK models, run virtual-population or batch patient simulations, do population PK variability studies, solve stiff pharmacology ODEs on an NVIDIA GPU, or mentions nvQSP, `nvqsp`, `sparse.solve`, RODAS4, or translating a pharmacology model into A0/A1/A2 coefficient form. Trigger this even when the user just describes a compartmental model ("two-compartment IV infusion across 1000 patients") without naming the library, since encoding the model into nvQSP's sparse-tensor form and getting the dosing/CSR conventions right is the hard part.

nvQSP

nvQSP is a GPU-accelerated RODAS4 stiff ODE solver for QSP and PBPK population studies. It solves a whole batch (virtual population) of patients in parallel on one NVIDIA GPU. The work this skill exists to get right is translating a pharmacology model into nvQSP's restricted polynomial form and building the sparse coefficient triples correctly — that is where mistakes happen.

The one constraint that governs everything

nvQSP solves only systems of this exact form:

dy/dt = A0 + A1·y + A2·(y⊗y)

Term	Shape	Captures
A0	(neq,) or (batch, neq)	Zeroth-order: constant synthesis, zero-order infusion
A1	(neq, neq) sparse CSR	First-order: linear elimination, inter-compartment rates
A2	(neq, neq, neq) sparse	Second-order: bilinear / mass-action (binding) terms

Before writing any code, check the model fits. nvQSP covers all linear PBPK models, first-order absorption, IV bolus/infusion, and bimolecular mass-action kinetics (drug–receptor binding, second-order TMDD approximation).

It does NOT cover, and you must stop and tell the user, if the model needs: Michaelis–Menten elimination, Hill-function PD, TMDD with quasi-steady-state, indirect-response (turnover) models, or any DAE system. These are not representable as a degree-2 polynomial and silently forcing them in will give wrong answers. If the user needs one of these, say so plainly rather than approximating without flagging it.

Encoding a model — the core procedure

For each state equation dy_i/dt, sort every term by polynomial degree:

1. Constant term (no y) → set A0[i]. Example: zero-order infusion rate.
2. Linear term k·y_j → add CSR entry at row i, column j, value k. Elimination -k·y_i is a diagonal entry A1[i,i] = -k; transfer into i from j is A1[i,j] = +k_ji.
3. Quadratic term k·y_a·y_b → add an A2 entry (row=i, col1=a, col2=b, val=k). The solver evaluates each entry literally as val * y[col1] * y[col2] with no implicit symmetrization, so add the term exactly once with its rate constant as written. A squared term k·y_a² is (i, a, a, k).

Sign convention: a species that is consumed gets a negative value; one that is produced gets a positive value. Mass-action binding L + R ⇌ C with on-rate kon, off-rate koff becomes: -kon on L and R (A2), +kon on C (A2), +koff on L and R (A1), -koff on C (A1).

Always sanity-check the encoding against a case with a known analytic answer (e.g. one-compartment decay y = y0·exp(-k·t)) before trusting a complex model. The verification snippet in references/api.md does exactly this.

Two mandatory gotchas

These cause the most failures:

• A2 must have ≥ 1 entry even for a purely linear model. If your model has no quadratic terms, insert a single negligible placeholder (A2_val = [1e-30] pointing at (i=0, col1=0, col2=0)). Omitting it raises ValueError: A2_nnz == 0.
• Doses always add to state index 0. The dosing compartment must be state 0. If the drug enters a different compartment, reorder your state vector so the dosing compartment is index 0. Doses are (time, amount) pairs; per-patient PK variation goes through A0/A1/A2 values, not through the dose schedule (which is shared across the batch).

CSR format details

• A1_csr = (rowptr, col, val) — standard CSR. rowptr: int32 (neq+1,), col: int32 (nnz,), val: float64 (nnz,).
• A2_csr = (rowptr, col1, col2, val) — a CSR-like 3-tensor. rowptr: int32 (neq+1,), col1/col2: int32 (nnz,), val: float64 (nnz,).
• Build A1 from a dense or scipy.sparse matrix; use the helpers in scripts/build_model.py to avoid hand-constructing CSR arrays.
• Dtypes matter: index arrays must be int32, value arrays float64.

Minimal working example

import numpy as np
from scipy.sparse import csr_matrix
from nvqsp import sparse
from nvqsp.options import SparseOptions

# Two-compartment model, 100 patients
A1 = csr_matrix([[-0.3, 0.1], [0.3, -0.1]])
A1_csr = (A1.indptr.astype(np.int32), A1.indices.astype(np.int32),
          A1.data.astype(np.float64))

# No quadratic terms -> epsilon placeholder so A2_nnz >= 1
A2_csr = (np.array([0, 1, 1], dtype=np.int32),  # rowptr
          np.array([0], dtype=np.int32),         # col1
          np.array([0], dtype=np.int32),         # col2
          np.array([1e-30]))                      # val

result = sparse.solve(
    A0=np.array([0.0, 0.0]),
    A1_csr=A1_csr,
    A2_csr=A2_csr,
    y0=np.tile([10.0, 0.0], (100, 1)),     # initial conditions per patient
    times=np.linspace(1.0, 24.0, 48),       # strictly increasing, positive
    doses=[(0.0, 100.0)],                    # 100 mg bolus into state 0 at t=0
    opts=SparseOptions(rtol=1e-6, atol=1e-9),
)

print(result.y.shape)   # (100, 48, 2)  -> (batch, n_times, neq)
print(result.steps)     # total ODE steps across all patients

result.y[b, t, i] is state i of patient b at times[t].

Population (virtual patient) variation

The sparsity structure (rowptr, columns) is shared across the batch; only the values differ per patient. Pass (batch, nnz)-shaped value arrays to vary parameters across the population (e.g. lognormal clearance):

batch = 1000
cl = np.random.lognormal(np.log(0.5), 0.2, size=batch)   # per-patient clearance
A1_val_batch = np.tile(A1_val_base, (batch, 1))          # (batch, A1_nnz)
A1_val_batch[:, cl_index] = -cl                          # vary the CL entry
A1_csr = (A1_rowptr, A1_col, A1_val_batch)

A0 accepts (batch, neq), and both A1/A2 val accept (batch, nnz). Results are bit-exact vs running each patient singly. See references/api.md for the full per-patient layout.

Running and tuning

• Defaults (opts=None) work for most QSP/PBPK models.
• Too many steps / too slow → relax rtol toward 1e-4.
• Oscillation or instability → tighten linsol_rtol or raise linsol_max_iters.
• Very stiff systems (stiffness ratio > 1e8) → raise max_steps to ~50000.
• times must be strictly increasing and positive.

Install & environment (only if the user is setting up)

pip install nvqsp (wheel bundles the compiled library; needs NumPy, optional SciPy). Requires Linux x86_64 + an NVIDIA Ampere/Ada/Hopper GPU (sm_80/sm_89/ sm_90), driver 525+, CUDA runtime 12.0+. No CUDA Toolkit needed to run. Volta and Turing are not supported (invalid device function error). C/C++ users get a .deb or standalone .so. Full install matrix and troubleshooting in references/api.md.

When to read the reference

references/api.md has the complete parameter and return tables, the full SparseOptions field list with defaults, the C API signature and memory layout, backend-selection notes, dosing conventions, the library search order, and install/troubleshooting detail. Read it when you need an exact signature, a field default, the C API, or are debugging an install/runtime error.

scripts/build_model.py provides dense_to_a1_csr, build_a2, and epsilon_a2 helpers plus a runnable verification example — prefer these over hand-writing CSR arrays.