Nvqsp

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.

Which solver

nvQSP ships two solvers. Pick before writing code.

Default to sparse for QSP/PBPK. These models are stiff, and it is the path that works straight out of pip install. Reach for tsit5 only when the model does not fit the polynomial form, or when it is genuinely non-stiff.

The one constraint that governs everything

nvQSP solves only systems of this exact form:

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

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.

Gradients (sparse RODAS4)

nvqsp.gradients.solve() returns the primal trajectories and derivatives with respect to A0, the sparse A1/A2 values, and the initial conditions. These are continuous forward sensitivities computed analytically

• not finite differences - so they are accurate to solver tolerance and suitable for parameter fitting and model training.

You describe what to differentiate by supplying CoefficientDerivatives: the derivative of each coefficient array with respect to your parameters. Omitted fields are treated as zero.

import numpy as np
from nvqsp import gradients
from nvqsp.gradients import CoefficientDerivatives

# One parameter, CL, which enters through a single A1 value.
dA1 = np.zeros((A1_nnz, 1))          # (A1_nnz, P)
dA1[cl_index, 0] = -1.0              # d(A1_val[cl_index]) / d(CL)

result = gradients.solve(
    A0, A1_csr, A2_csr, y0, times,   # first five as in sparse.solve()
    CoefficientDerivatives(parameter_names=["CL"], dA1=dA1),
    doses=doses,
)
result.y          # primal,   (batch, n_times, neq)
result.dy_dtheta  # Jacobian, (batch, n_times, neq, P)  <- parameter-last
result.names      # parameter axis labels

derivatives is the sixth positional argument and is required. Shapes are (neq, P) for dA0/dy0, (A1_nnz, P) for dA1, (A2_nnz, P) for dA2, or the same with a leading batch dimension for per-patient derivatives.

Do not confuse the two gradient entry points. CoefficientDerivatives drives the sparse analytic path here. GradientRequest belongs to the TSIT5 finite-difference path below. They are not interchangeable.

Cost scales linearly in the number of differentiated parameters. The augmented state has neq * (P + 1) entries, so differentiating 50 parameters is roughly 50x the work of a plain solve. Keep P to the parameters you actually need to fit.

For PyTorch training loops use the optional bridge: nvqsp.torch exposes an autograd Function over coefficient tensors, and SensitivityTargets limits which slots are differentiated. Ordinary inference skips the augmented solve entirely. Install with pip install nvqsp[torch]; importing nvqsp alone never imports PyTorch.

Dense TSIT5

TSIT5 is specialized per model: you generate and compile a CUDA library once, then solve from it many times.

from nvqsp import tsit5
from nvqsp.gradients import GradientRequest, GradientTarget

build = tsit5.build_model(model, output_dir="./build")   # needs nvcc
sol = tsit5.solve(                       # all args after library_path are keyword-only
    build.library_path, y0=y0, theta=theta, times=times,
)
sol.y                                     # (batch, n_times, neq)

sol_g = tsit5.solve_with_gradients(
    build.library_path, y0=y0, theta=theta, times=times,
    request=GradientRequest(target=GradientTarget.THETA),   # or GradientTarget.Y0
)

solve_with_gradients() returns a (batch, time, state, selected_input) Jacobian with respect to selected theta or y0, by central finite differences - so unlike the sparse path, accuracy depends on the step size. Cross-check with tsit5.reference_solve_model_with_gradients() and tsit5.validate_gradients(), which compare against a tight SciPy CPU reference (pip install nvqsp[reference]). tsit5.solve_torch() gives a differentiable PyTorch op.

build_model() needs nvcc, and nothing else

The wheel bundles the dense TSIT5 CUDA sources (nvqsp/src/dense/, nvqsp/src/substrate/), so build_model() works from a plain pip install nvqsp provided a CUDA Toolkit is present. It compiles a model-specialized .so once; solve() and solve_with_gradients() then consume that library and need no toolkit at all.

If the user has no nvcc, say so plainly rather than emitting a build call that cannot work — but do not tell them they need a source checkout. (Earlier 0.2.0 builds did require one; that was fixed before release.) build_model() also accepts repo_root= to point at a checkout explicitly, which is only needed for development against modified CUDA sources.

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

pip install nvqsp (v0.2.0; the wheel bundles the compiled sparse library and needs only NumPy). Extras: nvqsp[torch] for the autograd bridge, nvqsp[reference] for the SciPy CPU cross-check used by the TSIT5 validators.

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 the sparse solver; nvcc is needed only for tsit5.build_model(), whose CUDA sources ship inside the wheel. 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.