Six Weeks Inside a Los Alamos Simulation Engine
Since September 21, 2026, I have been contributing to bngsim, Los Alamos National Laboratory's open-source simulation engine for biochemical reaction networks. In that time I opened 204 pull requests against it. 171 have been merged, 25 are open as I write this, and 8 I closed myself. The merged ones added about 77,000 lines and removed about 6,800.
This post is about what that work actually was. It covers what bngsim is, where my changes landed, and the handful of bugs that taught me the most.
What bngsim Is
Cell biologists describe signaling pathways as networks of reactions. A receptor binds a ligand, a kinase phosphorylates a protein, a complex forms and falls apart. BioNetGen is the long-standing tool for writing those models as rules and expanding them into networks. For years, actually simulating one meant BioNetGen shelling out to a separate program, run_network, through files and a Perl toolchain.
bngsim replaces that driver. It is a C++ simulation kernel with Python bindings, built on SUNDIALS 7 (CVODE and CVODES), that runs models in-process and hands the results back as NumPy arrays. It loads BioNetGen .net and .xml files, Antimony, and SBML. It runs deterministic ODE simulations, stochastic ones (SSA), and two network-free simulators, and it computes forward sensitivities and steady states. It is a Los Alamos open-source release (LANL reference O5098) under the MIT license.
Here is what using it looks like.
Watching It Work
A quickstart only proves so much, so I ran a real model through it. bngsim's benchmark set includes Dolan2015.bngl, taken from the published-models collection of RuleWorld's RuleHub. It is a model of how a cell responds to DNA damage. It couples the p53 signaling network to the non-homologous end joining repair machinery, with Ku, DNA-PK, Ligase IV, PARP and Ligase III working on 50 sites of DNA. It has 2,339 reaction rules and 176 observables. Its authors ran it network-free, with NFsim simulating the molecules themselves instead of a pre-built list of every possible complex. bngsim ships NFsim in-process, so I did the same.
BioNetGen turned the rules into an XML description in 2.5 seconds, and bngsim simulated 40 hours of the unirradiated cell in 55 seconds. As shipped, though, the model is mostly quiet. Its radiation source starts at zero, and the authors' irradiation protocol, a pulse from minute 174 to minute 180, is commented out in the file.
So I turned it back on from the outside. bngsim's NfsimSession keeps a single simulation alive across segments, which means you can add molecules to a running cell and take them away again, without editing the model.
Here is what the 50 DNA sites did, replayed frame by frame from that run.
VIEW DATA
| Peak complete foci | 33 sites |
| Reactive oxygen during pulse | 95,954 molecules |
| Frames | 263, every 4 minutes from minute 150 to 1,200 |
The pulse floods the cell with about 96,000 molecules of reactive oxygen, and within three minutes nearly all of the Ku repair protein is oxidized. DNA breaks appear across the sites at once, each marked by H2AX, and as ATM binds those marks they become complete damage foci, peaking at 35 sites around minute 224. Repair is slow. It is past minute 1,300 before the count falls back to five or fewer.
The signal that matters for the cell's fate comes later. Phosphorylated p53 peaks at 150 molecules around minute 480. That drives up p21, the protein that stops the cell cycle, which crosses the model's threshold of 15 around minute 1,174. By minute 1,704 the senescence switch has flipped. The cell stops dividing for good and scales back its Ku and PARP.
VIEW DATA
| Dose 1,000: peak complete foci | 25 at min 230 |
| Dose 1,000: peak phospho-p53 | 142 at min 360 |
| Dose 1,000: p21 at 40 h | 37 |
| Dose 2,000: peak complete foci | 31 at min 220 |
| Dose 2,000: peak phospho-p53 | 206 at min 360 |
| Dose 2,000: p21 at 40 h | 37 |
| Dose 6,000: peak complete foci | 33 at min 230 |
| Dose 6,000: peak phospho-p53 | 150 at min 480 |
| Dose 6,000: p21 at 40 h | 22 |
The model is stochastic, so a single run is one possible cell, not the answer. Running the same unirradiated cell six times with six different seeds shows how much the chemistry alone varies.
VIEW DATA
| Seed 1 | 55s to simulate, p53 at 40 h: 136 |
| Seed 2 | 121.7s to simulate, p53 at 40 h: 63 |
| Seed 3 | 94.1s to simulate, p53 at 40 h: 149 |
| Seed 4 | 156.9s to simulate, p53 at 40 h: 152 |
| Seed 5 | 87.3s to simulate, p53 at 40 h: 96 |
| Seed 6 | 97.3s to simulate, p53 at 40 h: 168 |
That spread is the point of simulating molecule by molecule. A deterministic solver would draw one smooth line through the middle and hide how differently individual cells can behave.
Bill Hlavacek
bngsim is led by Bill Hlavacek at Los Alamos, who has written the large majority of it, more than 2,000 commits to my 247. It is hard to think of anyone better placed to build it, because he helped build the thing it runs.
His whole career since his Ph.D. has been at Los Alamos. According to his profile at the lab's Center for Nonlinear Studies, he finished his Ph.D. at the University of Michigan in 1996 and came to the lab that year as a postdoctoral research associate in the Theoretical Biology and Biophysics Group. He stayed on as a limited-term technical staff member from 2000, became a technical staff member in the same group in 2004, and served as acting deputy director of the Center for Nonlinear Studies from 2006 to 2007. He sits on the center's executive committee and has held an adjunct appointment in the University of New Mexico's biology department since 2003.
The publication that matters most here is the 2004 paper in Bioinformatics that introduced BioNetGen, "software for rule-based modeling of signal transduction based on the interactions of molecular domains," which he wrote with Blinov, Faeder and Goldstein. The idea is to describe a model as rules about how molecular domains interact, instead of listing every possible complex by hand. A year earlier, with Faeder, Blinov, Perelson and Goldstein, he had written "The complexity of complexes in signal transduction" in Biotechnology and Bioengineering, and his Center for Nonlinear Studies profile lists "Rules for Modeling Signal-Transduction Systems" as a 2006 research highlight. More than two decades later, bngsim is the engine that runs those models.
His work goes well beyond one tool. In 2004 alone he co-authored a review of mathematical and computational models of immune-receptor signaling in Nature Reviews Immunology with Goldstein and Faeder, and "Design of gene circuits: lessons from bacteria" in Nature Reviews Genetics with Wall and Savageau. His research interests cover signal transduction in the immune system, the biophysics of ligand-receptor binding, the design principles of bacterial gene regulation, and metabolomics and metabolic networks in bacteria.
He merged 152 of my 171 merged pull requests himself. Working on a codebase with one of BioNetGen's original authors changes how you write code. Getting a number almost right isn't an option when the person reviewing it knows what the model is supposed to do.
Ely Miller, a systems biology Ph.D. student at Northern Arizona University, also contributes. Most of my own work has been on the parts of a numerical engine that are easy to get almost right, like solver edge cases, model loading, and the testing that proves two answers agree.
Six Weeks in Numbers
The work came in two big bursts, the week I started and the week of October 19, with a quiet stretch in between.
VIEW DATA
| Week of Sep 21 | 64 |
| Week of Sep 28 | 7 |
| Week of Oct 5 | 0 |
| Week of Oct 12 | 20 |
| Week of Oct 19 | 74 |
| Week of Oct 26 (through Oct 28) | 6 |
Most of the line count is not engine code. Of those 77,000 added lines, 43% are tests and 31% are parity harnesses, benchmarks and recorded reference data from other simulators. Only about a fifth is the engine itself. In a project like this, proving the answer is right is most of the work.
VIEW DATA
| Through the week of Sep 21 | 13,001 |
| Through the week of Sep 28 | 15,429 |
| Through the week of Oct 5 | 15,429 |
| Through the week of Oct 12 | 23,364 |
| Through the week of Oct 19 | 58,039 |
| Through Oct 28 | 76,999 |
By area, SBML loading took the most merges. SBML is the standard interchange format for these models, and it allows a lot, including compartments whose sizes change over time, rules that depend on other rules, events that fire mid-step. Each of those is a place where an engine can quietly read a model differently from the way its author meant it.
VIEW DATA
| SBML loading | 23 |
| Events | 12 |
| CVODE solver | 12 |
| Steady states | 10 |
| Sensitivities | 9 |
| JAX translator | 9 |
| Jacobians | 7 |
| SSA | 7 |
| Simulator | 6 |
| Parity checks | 5 |
| Everything else | 71 |
The Bug That Only Happened in Germany
My favorite fix is PR #1238. bngsim's C++ read numbers with strtod, std::stod and atof, and wrote them with snprintf. All of those follow the C library's current locale. Python's own float() does not, which is why nobody had noticed.
So on a machine where a program had set the numeric locale to German or French, where the decimal separator is a comma, bngsim's C++ read 0.5 as 0. The damage depended on where the number was.
- The statistical factor BioNetGen writes for every homodimer reaction,
0.5*kf1, made the model fail to load. - An Antimony model with
k = 0.1loaded withk = 0, so the reaction had no rate at all, and no error was raised. - An inline table's values were read wrong and the run finished with the wrong answer, also without an error.
Silent wrong answers are the worst kind in scientific software. The fix creates one C locale once and routes every number the C++ side reads or writes through it, whatever locale the host process has set.
The pull request also added a test file that runs the affected constructs under a comma-decimal locale, so the bug cannot come back unnoticed.
Two Paths, One Answer
A pattern kept showing up in bngsim's bug reports. Two code paths that should give the same answer, and don't. Running a simulation in one go versus in two halves. Compiled code versus the interpreter. A sensitivity run versus a plain run. Each pair should agree to within solver tolerance, and when one disagrees, that is a bug, even with no external reference to compare against.
PR #747 pinned those invariants down on nine small inline models, including one stochastic check: the mean of many SSA runs against the exact mean from the chemical master equation. Here is the simplest one.
Several of those invariants failed on the day they were written. They went in marked as expected failures, each tied to its own issue, so that fixing a bug flips its test from expected-fail to pass, and the suite says so.
PR #1291, the largest change I've made at about 15,800 added lines, extended the idea to every model in the parity corpus. Each model now runs through the analytic Jacobian and a finite-difference one, the sparse KLU solver and the dense one, compiled code and the interpreter, and a German-locale child process against the C locale. It also checks three things no comparison of answers can catch, such as whether the backend you asked for is the one that actually ran.
Making CI Catch What Review Can't
Some of the most useful changes weren't to the engine at all.
PR #1163 closed a gap where no CI job compared a run against the committed reference trajectories. A change that shifted every trajectory would have left CI green. Every pull request is now verified against those references, with a tolerance that tells platform-level noise apart from real drift.
PR #1157 runs the C++ test suite under AddressSanitizer and UndefinedBehaviorSanitizer, and builds the library with compiler warnings on. Memory errors in a numerical kernel often show up only as slightly wrong numbers, so this is the cheapest insurance there is.
A series of pull requests under issue #805 recorded the outputs of the reference simulators bngsim is checked against (RoadRunner, AMICI and the legacy BioNetGen path) so that a parity sweep can replay them later without installing those tools at all.
Differentials
Almost everything bngsim does comes down to differential equations. A model's reactions say how fast each species changes, dx/dt = f(x, p), where x is the state and p the parameters, and CVODE steps that system forward in time. A sensitivity asks a second question, how much the answer moves when a parameter does. That derivative, ∂x/∂p, obeys a differential equation of its own, built from the Jacobian of f with respect to the state and the partial derivative of f with respect to each parameter. CVODES solves the two together.
The catch is where those partial derivatives come from. bngsim writes them out exactly, with sympy, when it can. When it can't, CVODES estimates them with a difference quotient, nudging each parameter a little and watching how f changes. That works for a smooth rate law. It goes wrong at a corner, where the slope on one side is not the slope on the other, and a narrow pulse is almost nothing but corners.
Issue #1227 was exactly that. A rate law that switched a reaction on with clamp() had no sympy counterpart, so the run fell back to the difference quotient. With a pulse 0.01 wide, the sensitivity it returned was off by as much as five-fold, and tightening the tolerance didn't make it converge. The same pulse written with min() and max() was fine, because those were already differentiated exactly.
PR #1242 gave clamp a derivative of its own. clamp(lo, x, hi) is lo when x is below lo, hi when x is above hi, and x in between. So its derivative in each argument is a step, 1 in the region where that argument is the value and 0 everywhere else, split evenly at the edge. Here is the diff, which is git's own kind of differential.
One case needed a guard. With its bounds crossed, lo above hi, a clamp isn't a bend anymore but a jump, and differentiating a jump as if it were a bend would quietly drop the jump from the answer. So the pull request admits a clamp to the analytic path only when both bounds are numbers in order, and turns everything else away by name, the way bngsim treats anything it can't differentiate honestly.
VIEW DATA
| rtol default, before | -1.174339 (21%) |
| rtol default, after | -0.974258 (0.078%) |
| rtol 1e-6, before | -5.987923 (515%) |
| rtol 1e-6, after | -1.041962 (7.0%) |
| rtol 1e-8, before | -0.980685 (0.74%) |
| rtol 1e-8, after | -0.973615 (0.012%) |
| rtol 1e-9, before | -0.90723 (6.8%) |
| rtol 1e-9, after | -0.973532 (0.0032%) |
After the change, clamp tracks the min and max spelling of the same pulse in every row of the pull request's table. At the tightest tolerance the error went from 6.8% to 0.003%. Bill merged it on October 25.
Steady States, Events and Sensitivities
The rest of the work is spread across the numerical core.
- Every steady state, not just one. PR #1165 enumerates all the steady states of a model by starting Newton's method from many points sampled on the model's conservation laws. Bistable switches, like the classic genetic toggle, have more than one.
- Delayed events in stochastic runs. The SSA engine used to refuse any event with a delay. PR #1159 made it run them with the same semantics as the ODE engine.
- Sensitivities through the awkward parts. Beyond clamp, several fixes taught the sensitivity code to carry a run past a conditional bend where it used to stall. Where a rate law sits exactly on an equality, it now refuses by name instead of quietly returning a derivative of zero.
What's Next?
I'm not done. I'm going to keep working on bngsim and keep contributing to it, and there is plenty queued up already. 25 of my pull requests are open right now, and most of them are differential testing, which means checking bngsim against itself and against other tools in ways nobody has to write by hand. That includes the analytic Jacobian against finite differences of the right-hand side, metamorphic transformations of a model file that should not change its answer, randomly generated models with every finding shrunk to a minimal reproducer, and the two network-free simulators, NFsim and RuleMonkey, compared in distribution.
After those land, I want to keep going where bngsim can still be wrong without anyone noticing, in the solvers, the model loaders, and the checks that prove two answers agree. I'll write about what I find as I go.
All of it is public at github.com/lanl/bngsim, and the documentation is at bngsim.readthedocs.io. If you simulate biochemical networks in Python, it's worth a look.
Sources: github.com/lanl/bngsim (README, contributor list, and pull requests #747, #1157, #1159, #1163, #1165, #1238 and #1291, plus the #805 series), William Hlavacek's profile at the Los Alamos Center for Nonlinear Studies, and Dolan2015.bngl from bngsim's benchmark set (RuleWorld RuleHub, Published), simulated with bngsim 0.16.0, all retrieved or run October 28, 2026. Pull request counts and line totals are as reported by GitHub for pull requests authored by Montana.