# Nvqsp > Build and run GPU-accelerated Quantitative Systems Pharmacology (QSP) and PBPK population ODE simulations with the nvQSP library, and differentiate them. nvQSP provides a sparse polynomial RODAS4 stiff solver (prebuilt) and a model-specialized dense adaptive TSIT5 solver, both with first-class gradients. Use this whenever the user wants to simulate compartmental PK/PBPK models, run virtual-population or batch patient simulations, do population PK variability studies, fit or train model parameters against data, solve stiff pharmacology ODEs on an NVIDIA GPU, or mentions nvQSP, `nvqsp`, `sparse.solve`, `gradients.solve`, `tsit5`, RODAS4, TSIT5, parameter sensitivities, 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. - Skill: `nvidia-bionemo/nvqsp` (Agent Skill, multi-file: 4 files) - Install (CLI): `npx skillmds@latest add nvidia-bionemo/nvqsp` - Raw SKILL.md: https://api.skillmd.com/api/skills/nvidia-bionemo/nvqsp/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: AI & ML - Author: nvidia-bionemo (https://skillmd.com/u/nvidia-bionemo) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/nvidia-bionemo/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. | | `nvqsp.sparse` (RODAS4) | `nvqsp.tsit5` (dense TSIT5) | | --- | --- | --- | | Method | 4th-order Rosenbrock, implicit | 5th-order explicit Runge-Kutta (Tsitouras) | | Use for | **stiff** systems | non-stiff / mildly stiff | | Model form | restricted polynomial `A0 + A1·y + A2·(y⊗y)` | general model via codegen | | Ready to run? | **yes** - prebuilt library ships in the wheel | **no** - must compile a model-specialized library first | | Needs CUDA Toolkit? | no | yes, `nvcc`, for the build step only | | Gradients | analytic forward sensitivities | central finite differences | **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) ``` | 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 ```python 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): ```python 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. ```python 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. ```python 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.