nvQSP v0.2.0

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.

Distribution

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)

AI Agent Skill

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.)

Requirements

• 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.

Quick Install

Python (from PyPI):

pip install nvqsp

C/C++ (Debian/Ubuntu):

Download nvqsp_0.2.0_amd64.deb from the GitHub release, then:

sudo dpkg -i nvqsp_0.2.0_amd64.deb

See INSTALL.md for full details.

Quick Start

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

# Two-compartment model: dy/dt = A0 + A1*y + A2*(y x y)
neq = 2
A0 = 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 models
A2_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=0
result = 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 patients

See API_REFERENCE.md for the complete Python and C API.

Gradient Exposure

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:

import numpy as np
from nvqsp import CoefficientDerivatives, 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.0

gradient_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.

Dense TSIT5 Solver (with gradients)

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.

import numpy as np
from nvqsp import tsit5, 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.

Model Form

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.

Documentation

• 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

License

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.