Two models in one C++ program#
In this demo we generate C++ code for two different cardiac cell models and solve a single
cell of each, side by side, in one C++ program. This is something the
C code cannot do: every model’s C code defines functions with the
same names (rhs, init_state_values, …) and global constants (NUM_STATES, …), so two
models cannot be linked into one program. The C++ code instead puts each model in its own
namespace, and defines its functions inline and its sizes constexpr, so a header can also
be included from several source files.
We use the ToR-ORd model
(ToRORd_dyn_chloride.ode) and the
ten Tusscher-Panfilov model
(tentusscher_panfilov_2006_M_cell.ode).
First we import what we need
import shutil
import subprocess
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import gotranx
Generated files go in a separate build directory
here = Path.cwd()
build = here / "build"
build.mkdir(exist_ok=True)
Generating the C++ code#
We generate one header per model, each in a namespace we choose, with the generalized Rush-Larsen scheme
models = {
"tor_ord": "ToRORd_dyn_chloride.ode",
"ten_tusscher": "tentusscher_panfilov_2006_M_cell.ode",
}
for namespace, odefile in models.items():
ode = gotranx.load_ode(here / odefile)
code = gotranx.cli.gotran2cpp.get_code(
ode,
scheme=[gotranx.schemes.Scheme.generalized_rush_larsen],
namespace=namespace,
)
(build / f"{namespace}.hpp").write_text(code)
2026-10-21 05:32:21 [info ] Load ode /home/runner/work/gotranx/gotranx/examples/cpp-single-cell/ToRORd_dyn_chloride.ode
2026-10-21 05:32:22 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=clss
2026-10-21 05:32:22 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=nai
2026-10-21 05:32:22 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=ki
2026-10-21 05:32:22 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=cli
2026-10-21 05:32:31 [info ] Num states 52
2026-10-21 05:32:31 [info ] Num parameters 144
2026-10-21 05:32:38 [info ] Load ode /home/runner/work/gotranx/gotranx/examples/cpp-single-cell/tentusscher_panfilov_2006_M_cell.ode
2026-10-21 05:32:39 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=K_i
2026-10-21 05:32:39 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=Na_i
2026-10-21 05:32:39 [warning ] Dropping a guard whose replacement failed its numeric check value=0 var=Ca_i
2026-10-21 05:32:39 [info ] Num states 19
2026-10-21 05:32:39 [info ] Num parameters 53
This is the same as running
gotranx ode2cpp ToRORd_dyn_chloride.ode --scheme generalized_rush_larsen --namespace tor_ord -o build/tor_ord.hpp
gotranx ode2cpp tentusscher_panfilov_2006_M_cell.ode --scheme generalized_rush_larsen --namespace ten_tusscher -o build/ten_tusscher.hpp
Without --namespace, the namespace is the name of the ODE (here ToRORd_dyn_chloride and
tentusscher_panfilov_2006_M_cell). Each header starts like this
print("\n".join((build / "tor_ord.hpp").read_text().splitlines()[:10]))
#pragma once
#include <cmath>
#include <cstring>
namespace tor_ord
{
constexpr int NUM_STATES = 52;
constexpr int NUM_PARAMS = 144;
constexpr int NUM_MONITORED = 357;
// Parameter index
The C++ program#
The program main.cpp (shown in full at the bottom of this page) includes both headers. It
sets up each model through its namespace (tor_ord::init_state_values,
ten_tusscher::NUM_STATES, …), steps both cells for one beat (1000 ms) with a time step of
0.01 ms, and every 0.1 ms writes the membrane potential of each model to a CSV file. The
potential is called v in one model and V in the other, and each model’s own
state_index finds it.
We compile the program with the C++ compiler on the system
compiler = shutil.which("c++") or shutil.which("g++") or shutil.which("clang++")
if compiler is None:
raise RuntimeError("This demo needs a C++ compiler (c++, g++ or clang++)")
executable = build / "main"
subprocess.run(
[compiler, "-std=c++17", "-O2", "-I", str(build), str(here / "main.cpp"), "-o", str(executable)],
check=True,
)
CompletedProcess(args=['/usr/bin/c++', '-std=c++17', '-O2', '-I', '/home/runner/work/gotranx/gotranx/examples/cpp-single-cell/build', '/home/runner/work/gotranx/gotranx/examples/cpp-single-cell/main.cpp', '-o', '/home/runner/work/gotranx/gotranx/examples/cpp-single-cell/build/main'], returncode=0)
and run it
results = build / "results.csv"
subprocess.run([str(executable), str(results)], check=True)
CompletedProcess(args=['/home/runner/work/gotranx/gotranx/examples/cpp-single-cell/build/main', '/home/runner/work/gotranx/gotranx/examples/cpp-single-cell/build/results.csv'], returncode=0)
The results#
Finally we plot the action potential of both models
data = np.loadtxt(results, delimiter=",", skiprows=1)
t, v_tor_ord, v_ten_tusscher = data.T
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(t, v_tor_ord, label="ToR-ORd")
ax.plot(t, v_ten_tusscher, label="ten Tusscher-Panfilov")
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Membrane potential (mV)")
ax.legend()
fig.tight_layout()
plt.show()
Source code#
main.cpp#
// Solve a single cell of two different models in one C++ program, with the code generated
// by `gotranx ode2cpp`.
//
// Usage: main <output.csv>
//
// Each generated header puts its model in a namespace (tor_ord and ten_tusscher), so the
// two models' functions and constants (rhs, NUM_STATES, ...) do not clash.
#include <cstdio>
#include <cstdlib>
#include <vector>
#include "ten_tusscher.hpp"
#include "tor_ord.hpp"
static int lookup(int index, const char *name)
{
if (index < 0) {
std::fprintf(stderr, "No state named '%s'\n", name);
std::exit(1);
}
return index;
}
int main(int argc, char *argv[])
{
if (argc != 2) {
std::fprintf(stderr, "Usage: %s <output.csv>\n", argv[0]);
return 1;
}
const double dt = 0.01; // ms
const double t_end = 1000.0; // ms, one beat
const int save_every = 10; // save every 0.1 ms
// ToR-ORd
std::vector<double> tor_ord_parameters(tor_ord::NUM_PARAMS);
std::vector<double> tor_ord_states(tor_ord::NUM_STATES);
std::vector<double> tor_ord_values(tor_ord::NUM_STATES);
tor_ord::init_parameter_values(tor_ord_parameters.data());
tor_ord::init_state_values(tor_ord_states.data());
const int tor_ord_v = lookup(tor_ord::state_index("v"), "v");
// ten Tusscher-Panfilov
std::vector<double> ten_tusscher_parameters(ten_tusscher::NUM_PARAMS);
std::vector<double> ten_tusscher_states(ten_tusscher::NUM_STATES);
std::vector<double> ten_tusscher_values(ten_tusscher::NUM_STATES);
ten_tusscher::init_parameter_values(ten_tusscher_parameters.data());
ten_tusscher::init_state_values(ten_tusscher_states.data());
const int ten_tusscher_v = lookup(ten_tusscher::state_index("V"), "V");
std::FILE *out = std::fopen(argv[1], "w");
if (out == nullptr) {
std::fprintf(stderr, "Could not open %s\n", argv[1]);
return 1;
}
std::fprintf(out, "t,v_tor_ord,v_ten_tusscher\n");
const int num_steps = static_cast<int>(t_end / dt);
for (int step = 0; step <= num_steps; step++) {
const double t = step * dt;
if (step % save_every == 0) {
std::fprintf(out, "%g,%.10g,%.10g\n", t, tor_ord_states[tor_ord_v],
ten_tusscher_states[ten_tusscher_v]);
}
// Each scheme writes the states at t + dt into its `values`
tor_ord::generalized_rush_larsen(tor_ord_states.data(), t, dt, tor_ord_parameters.data(),
tor_ord_values.data());
tor_ord_states.swap(tor_ord_values);
ten_tusscher::generalized_rush_larsen(ten_tusscher_states.data(), t, dt,
ten_tusscher_parameters.data(),
ten_tusscher_values.data());
ten_tusscher_states.swap(ten_tusscher_values);
}
std::fclose(out);
return 0;
}