GPU-accelerated stiff ODE solvers for Quantitative Systems Pharmacology (QSP) and PBPK population studies, with first-class gradient support. This release provides a sparse polynomial RODAS4 solver (prebuilt; no CUDA Toolkit needed at runtime) and a model-specialized dense adaptive TSIT5 solver, each with gradient exposure.
| Channel | Install |
|---|---|
| PyPI | pip install nvqsp |
| GitHub Release | nvqsp_0.2.0_amd64.deb (C/C++ headers + lib) |
| GitHub Release | libsparse_rodas4.so (standalone shared library) |
This repository includes source for a skill that helps AI agents translate
QSP/PBPK compartmental models into nvQSP's A0, A1, and A2 coefficient
form.
The reviewable source for the skill lives in skills/nvqsp/. Maintainers
should edit that source tree directly.
The skill does not change nvQSP runtime requirements; users still need the Python package or C/C++ library installed as described below.
All binaries are fat binaries with native code for:
- sm_80 — Ampere (A100, A10)
- sm_89 — Ada Lovelace (L4, L40, RTX 4090)
- sm_90 — Hopper (H100, H200)
- compute_90 PTX — forward compatibility for future architectures (Blackwell, etc.)
- Linux x86_64
- NVIDIA GPU: Ampere (sm_80), Ada Lovelace (sm_89), or Hopper (sm_90)
- NVIDIA driver 525+ (CUDA runtime 12.0+)
- Python 3.8+ with NumPy (for the Python API)
- PyTorch 2.0+ only when using the optional autograd bridge
- CUDA Toolkit (
nvcc) only when building a dense TSIT5 model library
The sparse RODAS4 solver needs no CUDA Toolkit at runtime — it ships as a
prebuilt library. Only nvqsp.tsit5.build_model() (dense TSIT5 model
specialization) and building from source require nvcc.
Python (from PyPI):
pip install nvqspC/C++ (Debian/Ubuntu):
Download nvqsp_0.2.0_amd64.deb from the GitHub release, then:
sudo dpkg -i nvqsp_0.2.0_amd64.debSee INSTALL.md for full details.
importnumpyasnpfromscipy.sparseimportcsr_matrixfromnvqspimportsparsefromnvqsp.optionsimportSparseOptions# Two-compartment model: dy/dt = A0 + A1*y + A2*(y x y)neq=2A0=np.array([0.0, 0.0])
A1=csr_matrix([[-0.3, 0.1], [0.3, -0.1]])
A1_rowptr=A1.indptr.astype(np.int32)
A1_col=A1.indices.astype(np.int32)
A1_val=A1.data.astype(np.float64)
# A2 must have >= 1 entry; use epsilon for purely linear modelsA2_rowptr=np.array([0, 1, 1], dtype=np.int32)
A2_col1=np.array([0], dtype=np.int32)
A2_col2=np.array([0], dtype=np.int32)
A2_val=np.array([1e-30])
# 100 patients, 48 time points, dose of 100 mg at t=0result=sparse.solve(
A0=A0,
A1_csr=(A1_rowptr, A1_col, A1_val),
A2_csr=(A2_rowptr, A2_col1, A2_col2, A2_val),
y0=np.tile([10.0, 0.0], (100, 1)),
times=np.linspace(1.0, 24.0, 48),
doses=[(0.0, 100.0)],
opts=SparseOptions(rtol=1e-6, atol=1e-9),
)
print(result.y.shape) # (100, 48, 2)print(result.steps) # total ODE steps across all patientsSee API_REFERENCE.md for the complete Python and C API.
nvqsp.gradients computes continuous forward sensitivities of the same
polynomial model in a single augmented GPU solve. Supply the derivatives of the
direct coefficients and the initial state with respect to each user parameter:
importnumpyasnpfromnvqspimportCoefficientDerivatives, gradients# The parameter axis P is always last. This example differentiates a single# clearance parameter carried through the sparse A1 values.dA1=np.zeros((A1_val.size, 1)) # (A1_nnz, P)dA1[0, 0] =-1.0gradient_result=gradients.solve(
A0=A0,
A1_csr=(A1_rowptr, A1_col, A1_val),
A2_csr=(A2_rowptr, A2_col1, A2_col2, A2_val),
y0=np.array([10.0, 0.0]),
times=np.linspace(1.0, 24.0, 48),
derivatives=CoefficientDerivatives(
parameter_names=("clearance",),
dA1=dA1,
),
)
print(gradient_result.y.shape) # (1, 48, 2)print(gradient_result.dy_dtheta.shape) # (1, 48, 2, 1)The augmented system holds the original neq states plus one neq-state
sensitivity block per parameter, so memory and work grow linearly with the
number of differentiated parameters P. Forward mode is intended for a modest
number of parameters. These are derivatives of the continuous ODE solution;
adaptive step-size and controller decisions are not differentiation targets.
For training, install the optional dependency (pip install "nvqsp[torch]") and
use nvqsp.torch.solve. It accepts direct coefficient tensors and y0, returns
an ordinary PyTorch tensor, and uses the same continuous sensitivities during
backpropagation. SensitivityTargets restricts differentiation to selected
coefficient slots to bound augmented-system size. Importing nvqsp never
imports PyTorch, so inference-only deployments keep the minimal dependency set.
The bridge supports first-order VJPs and rejects higher-order autograd.
Fixed dose schedules may be used during a gradient solve, but dose times and dose amounts are not differentiation targets in this release.
For general (non-stiff and mildly stiff) systems, nvqsp.tsit5 runs an adaptive
5th-order TSIT5 integrator as a model-specialized CUDA library and exposes
central-finite-difference gradients with respect to parameters (theta) or
initial conditions (y0).
Unlike the sparse RODAS4 solver — which ships as one prebuilt library — dense
TSIT5 is specialized per model: nvqsp.tsit5.build_model() generates and
compiles CUDA for your model, so it requires a CUDA Toolkit (nvcc). Solving and
differentiating an already-built library does not.
importnumpyasnpfromnvqspimporttsit5, GradientRequest, GradientTarget# 1) Build a model-specialized library once (requires nvcc).build=tsit5.build_model(model, "artifacts/", cuda_arch="sm_80")
# 2) Solve a batch of trajectories on the GPU.solve=tsit5.solve(
build.library_path,
y0=y0, # (neq,) or (batch, neq)theta=theta, # (P,) or (batch, P)times=np.linspace(0.0, 10.0, 64),
)
print(solve.y.shape) # (batch, n_times, neq)# 3) Gradients w.r.t. selected parameters.grad=tsit5.solve_with_gradients(
build.library_path,
y0=y0,
theta=theta,
times=np.linspace(0.0, 10.0, 64),
request=GradientRequest(target=GradientTarget.THETA, indices=None),
)
print(grad.gradients.shape) # (batch, time, state, n_selected)nvqsp.tsit5.reference_solve_model_with_gradients() and validate_gradients()
cross-check TSIT5 gradients against a tight SciPy CPU reference, and
nvqsp.tsit5.solve_torch exposes the solver as a differentiable PyTorch
operation. The gradient method is central finite differences; cost scales with
the number of requested coordinates.
The solver handles polynomial ODE systems of the form:
dy/dt = A0 + A1 * y + A2 * (y x y)
| Term | Shape | Meaning |
|---|---|---|
| A0 | (neq,) | Zeroth-order: constant synthesis, zero-order infusion |
| A1 | (neq, neq) sparse CSR | First-order: linear elimination, transfer rates |
| A2 | (neq, neq, neq) sparse CSR | Second-order: bilinear / mass-action terms |
Covers: all linear PBPK models, first-order absorption, IV bolus/infusion, bimolecular mass-action kinetics (drug-receptor binding, target-mediated disposition with second-order approximation).
Does not cover: Michaelis-Menten elimination, Hill-function PD, TMDD with quasi-steady-state, indirect response models, DAE systems.
- INSTALL.md — Installation and verification
- API_REFERENCE.md — Complete Python and C API reference
- CHANGELOG.md — Release notes and known limitations
- skills/nvqsp — AI agent skill source for nvQSP model encoding
This software is licensed under the NVIDIA Software License Agreement and the Product-Specific Terms for AI Products. By downloading, installing, or using this software you agree to the terms of both licenses.