Uncertainty Quantification in Dynamic Models of Biological Systems Using Conformal Prediction

This project is a Rust/C port of the CUQDyn1_Plus MATLAB toolbox, which combines parameter estimation and uncertainty quantification (UQ) in ordinary differential equation (ODE) models where only a subset of state variables are experimentally observed. Methodology described in Portela & Banga (2026) Hybrid Framework for Uncertainty Quantification in Partially Observed Dynamic Biological Systems.

Partial observability

The algorithm builds prediction bands two different ways depending on whether a state is measured:

State	Method	Basis
Observed	Conformal prediction over the leave-one-out ensemble	Distribution-free under exchangeability
Unobserved	Delta method through the parameter covariance	Linearised error propagation

The unobserved states can take their covariance two ways, selected with <uq_method>:

`uq_method`	Covariance	Upstream equivalent
`fim` (default)	Rank-aware FIM / Gauss-Newton	`CUQDyn1_Plus`
`hybridcov`	`D_FIM . R_LOO . D_FIM` — FIM marginal scale, correlation from the LOO ensemble	`CUQDyn1_Plus_HybridCov`

Observability is inferred from the data file: a state is unmeasured when its column is NaN after t = 0. A state is either measured at every point after t = 0 or at none of them — the format cannot express a gap, and a column that mixes NaN with values is rejected rather than silently treated as unobserved.

The same code path handles both cases; nothing branches on observability. With every state observed the bands are exactly the original CUQDyn1, since conformal prediction covers every state and the delta method has nothing left to write. The parameter covariance is still computed and reported (CovP, StdY), so existing .txt inputs keep their bands and gain the FIM diagnostics.

Trajectory sensitivities come from CVODES forward sensitivity analysis (CVodeSensInit), which replaces the complex-step differentiation used by the Matlab toolbox — CVODES cannot integrate in complex arithmetic, but it computes dy/dtheta natively in a single augmented integration.

Four ready-to-run partially observed examples ship in example-files/: linear-cascade, linear-cascade3, lv2-partobs (Lotka-Volterra with the prey measured and the predator hidden) and sir. See Using the CLI.

Note on the optimiser budget. The conformal bands come from the LOO ensemble and tolerate a mediocre global fit; the delta-method bands are anchored at parameters_init and inherit its error directly. If the bands for a hidden state look implausibly wide, raise <maxevaluation> in the eSS config before suspecting anything else.

Project Structure

The project is structured as follows:

deps/: Contains the dependencies of the project.
modules/: Contains the different modules of the project.
- cli/: CLI that enables the use of the cuqdyn-c library.
- cuqdyn-c/: The main library that implements the functionality of the paper project.
- sacess/: External proyect adapted to enable the use of eSS in C.
- cuqdyn-rs/: Rust crate. Parses the model from the XML file to evaluate it.
tests/: Constains the tests of the cuqdyn-c library.
CUQDyn/: Contains the original Matlab project.

Dependencies

xml2-2.9.1
cvodes-7.3.0 (Built by CMake)
hdf5-1.8.12 (Built by CMake)
gsl-2.7 (Built by CMake if not present)

ENV Vars

CUQDYN_DEF_YDOT

When running this software, there is a chance of generating an inf or NaN value. Take for example the following ODE:

<ode_expr y_count="1" p_count="2">
    p0 * y0 * (1 - y0 / p1)
</ode_expr>

If the value of p1 generated by the eSS solver is 0, the model eval produces an inf value. To avoid this, when the evaluation of the ODE generates a non-finite value, the software will set it to 0. If you want to change this behavior, you can set the value you want the variable to take when an invalid value is generated by setting the environment variable CUQDYN_DEF_YDOT to the desired value.

SACESS_SEED

Set this variable if you want to set the seed used by the eSS solver instead of generating a random one.

CUQDYN_CVODES_MIN_STEP

Set this variable if you want to set the minimum step used by the CVODES solver instead of the default value.

CUQDYN_CVODES_MAX_NUM_STEPS

Set this variable if you want to set the maximum number of steps used by the CVODES solver instead of the default value.

Building the project

GCC and GFortran v14 aren't supported by the project, so you need to use GCC v13 or lower.

The project has a scripts/build.sh script. build-[variant]/ directories will be created, each one representing a variant off the project builded. If you call the script without arguments, it will build all the variants, but you can also specify the variant you want, passing it as the first argument. The available variants are:

serial: Builds the project to only execute serial methods.
mpi: Builds the project to execute the sequential solver in parallel using MPI.
mpi2: Builds the project to only include the MPI and OpenMP defined in sacess-library.

After this, running scripts/test.sh is a good way to know if the cuqdyn library works as expected.

There is also a Dockerfile and a Docker Compose file to build and run the project in a container.

docker compose up

Using the CLI

Every example lives in its own folder under example-files/

example-files/<example>/
    data.txt               measurements, one column per state
    cuqdyn-fim.xml         algorithm config
    cuqdyn-hybridcov.xml   hybrid-covariance variant, where it applies
    sacess-serial.xml      optimiser config, serial
    sacess-mpi.xml         optimiser config, MPI

`<example>`	states	hidden	params	MPI	hybrid
`logistic`	1	0	2	yes
`lotka-volterra`	2	0	4	yes
`alpha-pinene`	5	0	5	yes
`linear-cascade`	2	1	2
`linear-cascade3`	3	2	3
`lv2-partobs`	2	1	4		yes
`sir`	3	2	2		yes
`nfkb`	15	5	29	yes

Examples with hidden states are partially observed: their data.txt carries NaN in the columns of the states that are never measured, so they exercise both band types. The rest are fully observed and reduce to the original CUQDyn1.

After building with scripts/build.sh, run any example by setting EXAMPLE:

EXAMPLE=lv2-partobs
VARIANT=serial
BUILD_TYPE=debug
mkdir -p output/${EXAMPLE}
./build/${BUILD_TYPE}-${VARIANT}/modules/cli/cli solve \
    -c example-files/${EXAMPLE}/cuqdyn-fim.xml \
    -s example-files/${EXAMPLE}/sacess-serial.xml \
    -d example-files/${EXAMPLE}/data.txt \
    -o output/${EXAMPLE}/

The MPI build takes the same shape, with sacess-mpi.xml and mpirun:

EXAMPLE=lotka-volterra
BUILD_TYPE=debug
mkdir -p output/${EXAMPLE}
mpirun -np 6 --use-hwthread-cpus ./build/${BUILD_TYPE}-mpi/modules/cli/cli solve \
    -c example-files/${EXAMPLE}/cuqdyn-fim.xml \
    -s example-files/${EXAMPLE}/sacess-mpi.xml \
    -d example-files/${EXAMPLE}/data.txt \
    -o output/${EXAMPLE}/

Swap cuqdyn-fim.xml for cuqdyn-hybridcov.xml on the examples that offer it to get the hybrid covariance instead: the hidden states keep the FIM marginal scale but take their correlation structure from the leave-one-out ensemble. That run also emits the plain FIM bands, so the two can be drawn together.

NF-kB is the largest case, and its FIM is close to singular (cond(J) ~ 6e8), so its hidden-state bands should be read together with the reported rank and condition number rather than taken at face value.

Every run prints a one-line FIM summary:

FIM: rank 4/4, condition number 7.050e+01, ridge 9.803e-08, sigma2 1

cuqdyn-hybridcov.xml adds the marginal standard deviations of both covariances, so the two can be compared at a glance:

HybridCov: FIM marginal scale with LOO correlation
   marginal std devs  FIM: 0.02264 0.004309 0.00155 0.03492
   marginal std devs  hybrid: 0.02264 0.004309 0.00155 0.03492

Everything above also runs inside the container:

docker compose up -d
docker compose exec cuqdyn_c ./build/debug-serial/modules/cli/cli solve \
    -c example-files/lv2-partobs/cuqdyn-fim.xml \
    -s example-files/lv2-partobs/sacess-serial.xml \
    -d example-files/lv2-partobs/data.txt \
    -o output/lv2-partobs/

After this, the file output/cuqdyn-results.txt contains the results of the algorithm but reading it as a plain text is not very useful. To fix this, you can run: (Needs python3 and matplotlib installed)

python3 plot.py output/cuqdyn-results.txt

Note: Be carefull when executing with mpirun, the number of precesses must be divisor of m - 1, where m is the number of rows in the input data matrix.

This will save a graphic representation for each y(t) in different png files inside the directory where the results are (output folder in this example). Each panel is labelled with the band type it carries, conformal or delta/FIM, and coloured accordingly.

To get information about all the options the cli supports, you can run the following command:

./build/{type}-{variant}/modules/cli/cli help

You can also run

./build/{type}-{variant}/modules/cli/cli version

to get the version of the cuqdyn-c lib, sacess lib and cli you are using.

Input files

There are three types of input files that must be provided:

eSS Solver config xml: This file contains the configuration of the eSS solver used in the sacess library. The specifications of this xml and how to build it are in this link.
cuqdyn config xml: This file contains the configuration of the cuqdyn solver used in the cuqdyn-c library.
- tolerances: rtol and atol used by the cvodes library.
- ode_expr: ODE model expression or identifier.
- time_scaling: Scaling factor for time. Using lower than 1 helps cvodes speed. You may have to increase the number of sacess maxevals. The sacess lb, up, and point constraints also get scaled. The same for the params median written in the output file. Params that are dividing shouldn't be scaled as the testing shows but they get scaled as well. This helps the solver converge faster but hasn't been proved to be valid for all cases. This has been made with testing purposes, not recommended to use it.
- y0: Initial conditions of the ODE.
- observed states: inferred from the data file rather than declared here. A state is unmeasured when its column is NaN after t = 0; see Partial observability.
There is an option to accelerate the process of evaluating the ODE by defining it inside the cuqdyn-rs module and compiling it. We will talk about this later.

CUQDyn1_Plus options (all optional, shown with their defaults):
- alp: predictive region level. Bands are nominally 1 - 2*alp, so 0.025 gives 95%. Defaults to 0.025.
- uq_method: fim (default) or hybridcov. Only affects the states that are never measured; see Partial observability.
- cost: how residuals are weighted before the fit and the FIM.
  - residual_model: none (raw residuals, the default), known_sigma (divide by the measurement standard deviations) or state_weights (multiply by explicit weights).
  - sigma: one value per observed state, or a single value applied to all.
  - sigma_is_known: when true the residuals are already standardized and the FIM variance scale is fixed at 1 instead of being estimated. Defaults to true.
  - observed_state_weights: one weight per observed state, for state_weights.
- fim: parameter covariance construction, used only when some state is hidden.
  - parameterization: log (default) or natural. The log form is more stable for positive rate constants; the covariance is always returned in natural units.
  - relative_ridge: regularizer relative to the largest squared singular value. Defaults to 1e-12.
  - rank_tol_factor: singular-value cutoff factor for the rank test. Defaults to 100.
```
<alp>0.025</alp>
<uq_method>fim</uq_method>
<cost>
    <residual_model>known_sigma</residual_model>
    <sigma>2.4531439541663325</sigma>
    <sigma_is_known>true</sigma_is_known>
</cost>
<fim>
    <parameterization>log</parameterization>
    <relative_ridge>1e-12</relative_ridge>
</fim>
```
```
<?xml version="1.0" encoding="UTF-8" ?>

<cuqdyn-config>
    <tolerances>
        <rtol>1e-8</rtol>
        <atol>1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8</atol>
    </tolerances>
    <ode_expr y_count="15" p_count="29">
        p19 - p20 * y0 - p16 * y0
        p16 * y0 - p18 * y1 - p17 * y1 * y7 - p20 * y1 - p1 * y1 * y9 + p2 * y3 - p3 * y1 * y12 + p4 * y4
        p18 * y1 + p17 * y1 * y7 - p20 * y2
        p1 * y1 * y9 - p2 * y3
        p3 * y1 * y12 - p4 * y4
        p10 * y12 - p0 * y5 * y9 + p4 * y4 - p22 * y5
        p22 * p21 * y5 - p0 * y10 * y6
        p14 * y8 - p15 * y7
        p12 + p11 * y6 - p13 * y8
        -p1 * y1 * y9 - p0 * y9 * y5 + p8 * y11 - p9 * y9 - p24 * y9 + p25 * y10
        -p0 * y10 * y6 + p24 * p21 * y9 - p25 * p21 * y10
        p6 + p5 * y6 - p7 * y11
        p0 * y9 * y5 - p10 * y12 - p3 * y1 * y12 + p23 * y13
        p0 * y10 * y6 - p23 * p21 * y13
        p27 + p26 * y6 - p28 * y14
    </ode_expr>
    <time_scaling>0.001</time_scaling> 
    <y0> 
        0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783, 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003
    </y0>
</cuqdyn-config>
```
Data file: The data file containing a mtrix of observed data and the initial value needed to solve the ODE. The data file should be a txt file written with the following format:
```
31 3            # matrix dimensions, so the parsing is easier
0 10 5          # row 1 is the initial condition for every state
1 16.27 -0.03196
.
.
30 9.117 4.548
```
Column 1 is time and every column after it is a model state, in order: y1, y2, ... yn. Row 1 (t = 0) must be finite for every state, hidden ones included, because it is the full initial condition used to integrate the ODE. If y0 is set in the cuqdyn config it takes precedence.

A state that is never measured carries NaN from t > 0 onwards:
```
41 3
0 10 0
0.5 NaN 1.69132433390405
1 NaN 2.66665964152082
```
That is the only thing that marks a state as unobserved. A file with no NaN is fully observed, and the bands reduce to the original CUQDyn1.

A column must be all NaN after t = 0 or have no NaN at all. A single stray one would otherwise demote the state to unobserved and throw away every real measurement it had, so it is reported as an error instead.

Output file

cuqdyn-results.txt holds one section per array. The first five are the original CUQDyn1 output and the rest are added by CUQDyn1_Plus:

Section	Shape	Contents
`Params`	n_params	Median of the leave-one-out parameter estimates
`Data`	m x n_states	Median of the leave-one-out trajectories
`Q_low` / `Q_up`	m x n_states	Prediction bands for every state
`Times`	m	Time points
`ObservedIdx`	n_obs	0-based indices of the measured states
`ParamsInit`	n_params	Best-fit parameters from the full-data fit
`MediaTot`	m x n_states	Trajectory at those best-fit parameters
`CovP`	n_params x n_params	Parameter covariance, only when a state is hidden
`StdY`	m x n_states	Delta-method standard deviations, same condition
`LooParams`	(m-1) x n_params	Per-replicate leave-one-out parameter estimates
`Q_low_fim` / `Q_up_fim`	m x n_states	Plain FIM bands, only with `uq_method=hybridcov`

ObservedIdx is what tells a reader which bands are conformal and which come from the delta method.

Defining a new model

Using strings and evaluate them is slow compared to compiled instructions, so, to make this possible, we designed a way to define the models in Rust and compile them to machine code.

Let's dig into it with an example of the Lotka Volterra model:

First of all, we need to write some Rust code. We will be using the following file modules/cuqdyn-rs/src/models.rs. Inside, we create a new unit struct and implement the Model trait like this:

#[derive(Default)]
struct LotkaVolterra;

impl Model for LotkaVolterra {
    fn eval(&self, _t: f64, y: &[f64], ydot: &mut [f64], p: &[f64]) {
        ydot[0] = y[0] * (p[0] - p[1] * y[1]);
        ydot[1] = - y[1] * (p[2] - p[3] * y[0]);
    }
}

Once the model is defined, we should give it an identifier:

pub fn eval_model_fun(model: &str, ode_expr: &OdeExpr) -> Box<dyn Model> {
    match model {
        "lotka-volterra" => Box::new(LotkaVolterra::default()),
        "alpha-pinene" | "α-pinene" => Box::new(AlphaPinene::default()),
        _ => Box::new(GenericModel::new(ode_expr))
    }
}

After you establish an identifier in the match at the bottom of the file, the model can be used indicating the identifier in the XML config file like this:

<?xml version="1.0" encoding="UTF-8" ?>

<cuqdyn-config>
    <tolerances>
        <rtol>1e-8</rtol>
        <atol>1e-9, 1e-10</atol>
    </tolerances>
    <ode_expr y_count="2" p_count="4">
        lotka-volterra
    </ode_expr>
</cuqdyn-config>

Note that y_count and p_count are still present.

Name		Name	Last commit message	Last commit date
.githooks		.githooks
.github/workflows		.github/workflows
CUQDyn		CUQDyn
CUQDyn1_Plus		CUQDyn1_Plus
benches		benches
deps		deps
example-files		example-files
example-results		example-results
modules		modules
scripts		scripts
tests		tests
toolchains		toolchains
.clang-format		.clang-format
.clang-tidy		.clang-tidy
.gitignore		.gitignore
CMakeLists.txt		CMakeLists.txt
Dockerfile		Dockerfile
LICENSE		LICENSE
README.md		README.md
TODO		TODO
deps.cmake		deps.cmake
docker-compose.yml		docker-compose.yml
matconverter.py		matconverter.py
plot.py		plot.py

Navigation Menu

Repository files navigation

Uncertainty Quantification in Dynamic Models of Biological Systems Using Conformal Prediction

Index

Partial observability

Project Structure

Dependencies

ENV Vars

CUQDYN_DEF_YDOT

SACESS_SEED

CUQDYN_CVODES_MIN_STEP

CUQDYN_CVODES_MAX_NUM_STEPS

Building the project

Using the CLI

Input files

Output file

Defining a new model

About

Releases

Packages

Used by

Contributors

Languages

Latest commit

History

Folders and files

Repository files navigation

Uncertainty Quantification in Dynamic Models of Biological Systems Using Conformal Prediction

Index

Partial observability

Project Structure

Dependencies

ENV Vars

CUQDYN_DEF_YDOT

SACESS_SEED

CUQDYN_CVODES_MIN_STEP

CUQDYN_CVODES_MAX_NUM_STEPS

Building the project

Using the CLI

Input files

Output file

Defining a new model

About

Resources

Stars

Watchers

Forks

Releases

Packages

Used by

Contributors

Languages