Solving a single cell in C

In this demo we generate C code for a cardiac cell model, and solve a single cell with a small, standalone C program. This is how you would use the generated code from a solver written in C. If you instead want to call the generated C functions from Python, see Compiling a C-extension.

We use the ToR-ORd model of a human ventricular cardiomyocyte. You can download the model in .ode format here.

First we import what we need

import re
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

Next we load the model and generate C code for it, including the generalized Rush-Larsen scheme that we will use to step the cell forward in time

ode = gotranx.load_ode(here / "ToRORd_dyn_chloride.ode")
code = gotranx.cli.gotran2c.get_code(
    ode, scheme=[gotranx.schemes.Scheme.generalized_rush_larsen]
)
(build / "model.h").write_text(code)

2026-10-21 05:31:39 [info     ] Load ode /home/runner/work/gotranx/gotranx/examples/c-single-cell/ToRORd_dyn_chloride.ode

2026-10-21 05:31:40 [warning  ] Dropping a guard whose replacement failed its numeric check value=0 var=ki

2026-10-21 05:31:40 [warning  ] Dropping a guard whose replacement failed its numeric check value=0 var=nai

2026-10-21 05:31:40 [warning  ] Dropping a guard whose replacement failed its numeric check value=0 var=cli

2026-10-21 05:31:40 [warning  ] Dropping a guard whose replacement failed its numeric check value=0 var=clss

2026-10-21 05:31:49 [info     ] Num states 52                 

2026-10-21 05:31:49 [info     ] Num parameters 144            

259533

This is the same as running the following command on the command line

gotranx ode2c ToRORd_dyn_chloride.ode --scheme generalized_rush_larsen -o build/model.h

The header contains the number of states, parameters and monitored values, functions that set the initial values, functions that look up the index of a state, parameter or monitored value by name, and the functions rhs, monitor_values and generalized_rush_larsen. For example, the signature of the scheme is

signature = re.search(r"void generalized_rush_larsen\([^)]*\)", code).group()
print(" ".join(signature.split()))

void generalized_rush_larsen(const double *__restrict states, const double t, const double dt, const double *__restrict parameters, double *values)

The C program

The program main.c (shown in full at the bottom of this page) includes the generated header, sets the initial states and parameters, and steps the cell for two beats (2000 ms) with a time step of 0.01 ms. Every 0.1 ms it writes the time, the membrane potential v, the intracellular calcium concentration cai and the fast sodium current INa (a monitored value) to a CSV file. It looks these up by name with state_index and monitor_index, so it does not depend on the order of the states.

We compile the program with the C compiler on the system

compiler = shutil.which("cc") or shutil.which("gcc") or shutil.which("clang")
if compiler is None:
    raise RuntimeError("This demo needs a C compiler (cc, gcc or clang)")
executable = build / "main"
subprocess.run(
    [compiler, "-O2", "-I", str(build), str(here / "main.c"), "-o", str(executable), "-lm"],
    check=True,
)

CompletedProcess(args=['/usr/bin/cc', '-O2', '-I', '/home/runner/work/gotranx/gotranx/examples/c-single-cell/build', '/home/runner/work/gotranx/gotranx/examples/c-single-cell/main.c', '-o', '/home/runner/work/gotranx/gotranx/examples/c-single-cell/build/main', '-lm'], 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/c-single-cell/build/main', '/home/runner/work/gotranx/gotranx/examples/c-single-cell/build/results.csv'], returncode=0)

The results

Finally we load the CSV file and plot the results

data = np.loadtxt(results, delimiter=",", skiprows=1)
t, v, cai, INa = data.T

fig, ax = plt.subplots(3, 1, sharex=True, figsize=(8, 8))
ax[0].plot(t, v)
ax[0].set_ylabel("v (mV)")
ax[1].plot(t, cai * 1000)
ax[1].set_ylabel("cai (µM)")
ax[2].plot(t, INa)
ax[2].set_ylabel("INa (µA/µF)")
ax[2].set_xlabel("Time (ms)")
fig.tight_layout()
plt.show()

Source code

main.c

// Solve a single cell with the C code generated by `gotranx ode2c`.
//
// Usage: main <output.csv>
//
// The generated header (model.h) provides everything used below: the number of states and
// monitored values, the initial values, index lookups by name, monitor_values and the
// numerical scheme (here generalized_rush_larsen).
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "model.h"

static int lookup(int index, const char *name)
{
    if (index < 0) {
        fprintf(stderr, "No variable named '%s' in the model\n", name);
        exit(1);
    }
    return index;
}

int main(int argc, char *argv[])
{
    if (argc != 2) {
        fprintf(stderr, "Usage: %s <output.csv>\n", argv[0]);
        return 1;
    }

    const double dt = 0.01;     // ms
    const double t_end = 2000.0;  // ms, two beats with the default 1000 ms stimulus period
    const int save_every = 10;  // save every 0.1 ms

    double *parameters = malloc(NUM_PARAMS * sizeof(double));
    double *states = malloc(NUM_STATES * sizeof(double));
    double *values = malloc(NUM_STATES * sizeof(double));
    double *monitors = malloc(NUM_MONITORED * sizeof(double));
    init_parameter_values(parameters);
    init_state_values(states);

    // Look the variables up by name, so the program does not depend on their order
    const int v = lookup(state_index("v"), "v");
    const int cai = lookup(state_index("cai"), "cai");
    const int INa = lookup(monitor_index("INa"), "INa");

    FILE *out = fopen(argv[1], "w");
    if (out == NULL) {
        fprintf(stderr, "Could not open %s\n", argv[1]);
        return 1;
    }
    fprintf(out, "t,v,cai,INa\n");

    const int num_steps = (int)(t_end / dt);
    for (int step = 0; step <= num_steps; step++) {
        const double t = step * dt;
        if (step % save_every == 0) {
            monitor_values(t, states, parameters, monitors);
            fprintf(out, "%g,%.10g,%.10g,%.10g\n", t, states[v], states[cai], monitors[INa]);
        }
        // The scheme writes the states at t + dt into `values`
        generalized_rush_larsen(states, t, dt, parameters, values);
        memcpy(states, values, NUM_STATES * sizeof(double));
    }

    fclose(out);
    free(parameters);
    free(states);
    free(values);
    free(monitors);
    return 0;
}