PyBioNetFit: Biological Model Parameterization

• PyBioNetFit is an open-source software that parameterizes biological models using structured text workflows and supports both BNGL and SBML formats.
• It integrates quantitative simulation data with qualitative constraints through BPSL to enhance reproducibility and extend classical residual-based fitting.
• The tool employs parallelized metaheuristic optimization and adaptive MCMC methods to tackle large-scale, non-gradient, and mixed-data inference problems.

PyBioNetFit, often abbreviated PyBNF, is an open-source software package for parameterization of mechanistic biological models that also supports uncertainty quantification, model checking, and design. It was developed for systems biology workflows in which model evaluation depends on repeated simulation and in which the available evidence may include both conventional quantitative measurements and formally encoded qualitative biological knowledge. A defining feature of the software is that it works directly with standard model formats—especially BioNetGen Language (BNGL) and Systems Biology Markup Language (SBML)—and uses structured text inputs to make calibration and inference workflows reproducible. Its most distinctive methodological contribution is the Biological Property Specification Language (BPSL), which allows qualitative system properties to be incorporated into parameter estimation, model checking, and design alongside standard numerical data (Mitra et al., 2019). Later work extended this foundation with likelihood-based Bayesian treatment of qualitative observations (Mitra et al., 2019), an adaptive Markov chain Monte Carlo method in PyBNF version 1.1.9 (Neumann et al., 2021), and a detailed mixed-data case study demonstrating parameterization and uncertainty quantification in a MEK/ERK signaling model (Miller et al., 26 Aug 2025).

1. Origins, scope, and historical position

PyBioNetFit was introduced as a general-purpose tool for systems biology model analysis with four principal task classes: parameterization, uncertainty quantification, model checking, and design. The motivating problem was not ordinary least-squares fitting in small, smooth ODE models alone, but a broader class of calibration settings in which gradients are unavailable, expensive, or scientifically awkward. The development rationale emphasized three problem types that prior workflows served poorly: very large ODE systems, especially those implied by rule-based models; stochastic simulation models, including network-free rule-based simulations; and use of unconventional data, especially qualitative or non-numeric observations. In such settings, manual parameter tuning or bespoke scripts were common, but both reduced reproducibility and scaled poorly (Mitra et al., 2019).

Historically, PyBioNetFit followed BioNetFit, an earlier tool for parameter identification in BNGL-encoded rule-based models. BioNetFit used a file-based workflow organized around BNGL models, EXP data files, CONF configuration files, and optional bootstrap analysis. It implemented a metaheuristic population-based global optimization procedure through an evolutionary algorithm and supported deterministic and stochastic simulators within the BioNetGen framework (Hlavacek et al., 2018). PyBioNetFit retained the emphasis on structured, plain-text workflows but broadened the scope in two ways. First, it natively supported both BNGL and SBML rather than BNGL alone. Second, it introduced BPSL as a domain-specific language for qualitative biological properties, thereby extending the range of usable evidence beyond purely numerical datasets (Mitra et al., 2019).

The software was positioned as complementary rather than universally superior to existing parameter-estimation environments. Data2Dynamics and COPASI were described as supporting fitting for conventional ODE models and including some metaheuristics, but not as being designed around parallelized metaheuristic optimization for difficult biological models, the same breadth of stochastic or rule-based workflows, or qualitative constraints in the PyBioNetFit sense. This positioning matters because PyBioNetFit’s role is not primarily to replace gradient-based tools in all contexts, but to make hard simulation-based inference problems tractable without problem-specific code (Mitra et al., 2019).

2. Software architecture and reproducible workflow

PyBioNetFit operates through a structured text-based workflow. In the forms described across its core papers, the principal inputs are a configuration file with extension .conf, one or more model-definition files in .bngl or .xml format, one or more quantitative data files with extension .exp, and, when qualitative knowledge is used, one or more .prop files containing BPSL statements. The command-line interface is given as pybnf -c <filename>.conf, after which the software reads the configuration and model or data files, proposes parameter sets, dispatches simulations, computes objective values, and updates the optimizer (Mitra et al., 2019).

The model-format layer is central to PyBioNetFit’s design. For BNGL models, PyBioNetFit uses BioNetGen and supports ODE simulation, stochastic simulation algorithm (SSA), and NFsim network-free stochastic simulation. For SBML models, it uses libRoadRunner, which by default uses CVODE for numerical integration. This architecture allows a single fitting framework to cover conventional biochemical ODE models, large rule-based ODE systems, and stochastic or network-free simulations without rewriting the mechanistic model in custom inference code (Mitra et al., 2019).

The workflow imposes explicit conventions. In BNGL jobs, free parameters to be fit must have names ending in __FREE, model outputs to be compared to data must be declared as observables or functions, and each simulation command must have a unique suffix. In SBML jobs, free parameters must be either SBML parameters or initial concentrations of species; outputs compared to measurements must be SBML species; and simulation protocols must be declared in the CONF file because SBML does not itself define them. Every simulation must have a suffix so that simulation output can be linked to the appropriate EXP file or BPSL reference. PyBioNetFit can also handle jobs with multiple model files, which is useful when different model instances share parameters, such as wild-type and mutant versions of the same signaling system (Mitra et al., 2019).

The file-based organization is treated as a reproducibility mechanism rather than a mere implementation detail. In the MEK/ERK study, PyBioNetFit was used with five BNGL files corresponding to WT, KO, N78G, T292A, and T292D model variants; quantitative time-series data were stored in EXP files; qualitative constraints were stored in PROP files; and a CONF file declared the adjustable parameters, box constraints, mappings between files and model variants, and optimization settings. The study explicitly argued that the model, the data, the fitting algorithm, the parameter bounds, and the mappings between data and model outputs are thereby captured in standardized text files rather than hidden in manual tuning decisions (Miller et al., 26 Aug 2025).

PyBioNetFit was also built to exploit parallel hardware. It interfaces with dask.distributed, is designed to run on multicore workstations and computer clusters, and treats each simulation or objective-function evaluation as a separate worker process. This is consequential because the intended use cases include large ODE systems, SSA runs, and NFsim simulations, all of which can be computationally expensive (Mitra et al., 2019).

3. BPSL and the formalization of qualitative biological knowledge

The Biological Property Specification Language is PyBioNetFit’s defining extension beyond conventional residual-based fitting. BPSL is a domain-specific language for declaring biological system properties in a human-readable but machine-evaluable form. It was introduced so that qualitative biological knowledge could be used in parameterization, model checking, and design, either alone or in combination with quantitative data (Mitra et al., 2019).

A BPSL statement has three parts: an inequality, an enforcement condition, and a weight. The inequality may compare a variable with a constant or two variables using one of <, >, <=, or >=. The enforcement condition determines where in the independent-variable domain the property is enforced. The core enforcement keywords are always, once, at <condition>, and between <condition1>, <condition2>. BPSL can also refer to outputs from different simulations or model variants, enabling case-control comparisons such as wild type versus mutant or untreated versus treated conditions (Mitra et al., 2019).

Computationally, the original PyBioNetFit treatment of qualitative constraints was a static penalty method. If a qualitative constraint is represented as $g(x) < 0$, where $x$ is the trial parameter set, then violation adds the penalty

$C \cdot g(x),$

with $C$ the user-specified weight. This makes qualitative observations neither hard feasibility conditions nor informal narrative comments: they become explicit objective-function terms whose magnitude reflects the extent of violation. The paper gave concrete examples such as

$X>5 \ \text{always weight } 2$

and

$X<1 \ \text{between time=8,Y=5 weight } 3,$

which translate into penalties based on the minimum or maximum simulated value over the relevant interval (Mitra et al., 2019).

The later Bayesian-UQ paper reformulated qualitative observations probabilistically. In PyBNF version 1.1.0, BPSL could specify qualitative likelihood terms through [confidence](https://www.emergentmind.com/topics/confidence) and tolerance for the two-category likelihood, or pmin, pmax, and tolerance for the many-category likelihood. This replaced the earlier heuristic interpretation of “constraint violation penalties” with explicit observation models based on latent Gaussian variables, thresholds, and discrepancy probabilities. The resulting formulation enabled posterior sampling, credible intervals, parameter-correlation analysis, and prediction uncertainty for models informed by qualitative or mixed qualitative and quantitative datasets (Mitra et al., 2019).

The 2025 MEK/ERK case study shows how BPSL functions in practice. Qualitative observations were encoded as propositions such as WT.MEK_pRDS at time = 300 < N78G.MEK_pRDS at time = 300, KO.MEK_pRDS at time = 300 > KO.MEK_pRDS at time = 1800, and N78G.MEK_pRDS at time = 300 > T292D.MEK_pRDS at time = 300. Supplementary Tables 1–5 reportedly enumerate 90 BPSL statements across the five cell-line conditions. In that study, qualitative constraints were combined with quantitative time-series data in a single parameter-estimation workflow, and the paper emphasized that this replaced “tedious trial-and-error manual tuning” with an automated and statistically interpretable procedure (Miller et al., 26 Aug 2025).

4. Optimization, likelihoods, and uncertainty quantification

PyBioNetFit’s optimization layer is centered on parallelized derivative-free metaheuristics. The recommended fitting algorithms in the 2019 introduction were Differential Evolution (DE), Asynchronous Differential Evolution (aDE), Particle Swarm Optimization (PSO), and Scatter Search (SS), with Simulated Annealing (SA) and Island-based Differential Evolution (iDE) also implemented. A parallelized simplex local optimizer is available for refinement. The stated reason for emphasizing metaheuristics is that they do not require gradients, remain applicable when gradients are expensive or unavailable, work with stochastic simulators, and parallelize naturally. The tradeoff is explicit: they do not guarantee convergence to the global optimum (Mitra et al., 2019).

The uncertainty-quantification layer evolved over time. Early PyBioNetFit supported Bayesian MCMC through Metropolis–Hastings and Parallel Tempering, as well as bootstrapping. At that stage, however, the software warned that Bayesian MCMC assumes the objective function represents a log likelihood and that this assumption was valid for the chi-squared objective used in the uncertainty example, but not recommended for arbitrary objectives or objectives including qualitative-data penalty contributions (Mitra et al., 2019). The qualitative-likelihood work addressed this gap by deriving explicit likelihood functions for binary and ordinal qualitative observations and by showing that mixed-data inference can be performed by summing the qualitative and quantitative negative log-likelihood terms (Mitra et al., 2019).

The MEK/ERK study made this mixed-data structure explicit. Its central objective was

$$F(\theta)=F_{\text{qual}}(\theta)+F_{\text{quant}}(\theta),$$

with the fitted parameter estimate defined by

$\hat{\theta}=\arg\min_{\theta\in\Theta}F(\theta),$

where $\Theta$ is the feasible region defined by box constraints. The paper stressed that this was not merely heuristic optimization: $F$ was connected to a likelihood, with the quantitative term arising from Gaussian noise assumptions and the qualitative term approximated by a logistic/Bernoulli observation model for pairwise comparisons (Miller et al., 26 Aug 2025).

A major later addition was the adaptive MCMC method am, introduced in PyBNF version 1.1.9. This sampler uses a multivariate Gaussian random-walk proposal and adapts the proposal covariance during sampling so that the learned covariance approaches the posterior covariance and the acceptance rate approaches the theoretical high-dimensional optimum of $x$0. The algorithm supports warm starts, continuation, user-specified initial locations and covariance matrices, multiple chains in parallel on Linux clusters, and likelihoods including normal and negative binomial forms. The paper recommends choosing the stationary-phase offset to be at least $x$1 iterations and emphasizes routine use of diagnostics such as trace plots, pairs plots, acceptance behavior, and effective sample size (Neumann et al., 2021).

Bootstrapping remains part of the PyBioNetFit workflow. In the earlier description, it is enabled in the CONF file by bootstrap=<number>, and the resulting confidence intervals quantify uncertainty in the location of the best-fit parameter set while reflecting both data variability and optimizer performance. This caveat is nontrivial: if the optimizer is imperfect, bootstrap intervals may widen because of optimization error rather than data uncertainty alone (Mitra et al., 2019).

5. Applications, benchmarks, and empirical behavior

The initial PyBioNetFit paper reported 31 example problems spanning published biological models, synthetic fitting tasks, all four supported simulation modes, and all four high-level application classes: parameterization, uncertainty quantification, model checking, and design. All were said to be solved acceptably with at least one algorithm using default settings. The benchmark suite also examined scaling on clusters using up to 288 cores. Median runtime generally decreased with additional cores, but no single optimizer dominated all problems: at 288 cores, PSO performed best on the stochastic SSA and NFsim benchmarks, whereas aDE and SS outperformed PSO on the other benchmarks. A major mixed-data demonstration was fitting a Tyson-group yeast cell-cycle control model in a 153-dimensional parameter space, using both qualitative mutant-phenotype constraints and quantitative expression data. PyBioNetFit achieved a minimum objective value of 80 versus 70 in the earlier custom-coded study and still produced reasonable consistency with both qualitative constraints and quantitative time-course data (Mitra et al., 2019).

PyBioNetFit was also demonstrated as a model-checking and design platform. In the autophagy-signaling example, ten qualitative behaviors were encoded in BPSL and evaluated against a reference model and eight hypothetical variants. Variants 1–4 remained consistent with all ten properties, whereas variants 5–8 failed one or more properties. In the design setting, drug concentrations $x$2 through $x$3 were treated as decision variables; the desired reduction of autophagic vesicle count below 20 under high stress was encoded as a BPSL inequality; and total drug dose was minimized as a quantitative term. The reported optimum used less total drug than any single-drug treatment and had only two nonzero drugs even though all six were allowed to vary (Mitra et al., 2019).

The adaptive-MCMC paper extended the empirical scope beyond canonical signaling models. It validated am on a synthetic linear-regression problem with known analytical posterior, a nonlinear HIV viral-dynamics model, a mast-cell degranulation benchmark with 16 adjustable parameters, and COVID-19 forecasting problems for 15 metropolitan statistical areas. Relative to the earlier fixed-proposal mh sampler, am was described as substantially more efficient. In the 16-parameter degranulation benchmark, the ratio $x$4 ranged from 3.03 to 27.7, with mean 11.5 and median 10.5. For the HIV and COVID-19 examples, PyBNF’s marginal posteriors and posterior predictive summaries were reported as indistinguishable from or consistent with those obtained using independent, problem-specific Bayesian code (Neumann et al., 2021).

The 2025 MEK/ERK case study is the most detailed account of PyBioNetFit’s mixed-data workflow. Starting from the Kocieniewski and Lipniacki model of the EGFR–SOS1–RAS–RAF–MEK–ERK cascade, the study fit 31 adjustable quantities—28 model parameters plus 3 scaling factors—across five model variants and 90 BPSL statements. Maximum-likelihood fitting used PyBioNetFit’s parallelized DE algorithm on a high-performance cluster with 25 CPUs on one node. The reported average wall-clock time was about 2 minutes with 2,825 objective-function evaluations. Supplementary Table 7 gave an overall objective score of 40.0 for the original manual parameterization and 15.0 for the PyBioNetFit parameterization; Supplementary Table 9 gave global RMSE values of 1.47 and 1.076, respectively, on the WT quantitative data. The resulting trajectories were described as qualitatively very similar to those of the original study despite substantial differences in parameter estimates, highlighting non-uniqueness of fitted parameter sets (Miller et al., 26 Aug 2025).

That same study also carried out both frequentist-style and Bayesian uncertainty quantification. Profile likelihoods for the rate constants $x$5 and $x$6 indicated that both were practically identifiable, with $x$7 more tightly constrained. For Bayesian analysis, the authors sampled a reduced six-dimensional parameter set consisting of $x$8, $x$9, three scaling factors, and a noise-model hyperparameter $C \cdot g(x),$0, using adaptive MCMC with a proper uniform prior defined by box constraints. Each of five independent chains used 25,000 burn-in iterations, 25,000 adaptation iterations, and 250,000 production iterations. The jobs were run on 25 CPUs, took about 5.15 days, and required 1.25 million objective-function evaluations. Reported convergence diagnostics included, for $C \cdot g(x),$1, ESS bulk 6160, ESS tail 8620, and $C \cdot g(x),$2, and, for $C \cdot g(x),$3, ESS bulk 2150, ESS tail 2120, and $C \cdot g(x),$4 (Miller et al., 26 Aug 2025).

6. Limitations, standards, and interoperability

PyBioNetFit’s strengths are accompanied by explicit limitations. The core optimization methods are metaheuristics and therefore do not guarantee global optimality. Performance is problem-dependent, and algorithm choice often requires experimentation. The 2019 paper also stated that gradient-based tools may be better for small, smooth ODE problems; in its comparison with Data2Dynamics, one benchmark was solved by Data2Dynamics in under 5 minutes, much faster than PyBioNetFit. BPSL was deliberately designed to be more readable than temporal logic formalisms such as LTL, but the tradeoff is reduced expressivity. Practical use also requires model-format adaptations such as naming conventions for free parameters, observables or species setup, suffix assignments, and, for SBML, protocol declarations in the CONF file (Mitra et al., 2019).

The statistical treatment of qualitative data likewise carries caveats. In the likelihood-based framework, the user must choose confidence, pmin, pmax, and tolerance values so that the implied sampling model is coherent. The software evaluates the formulas, but responsibility for specifying a valid qualitative observation model remains with the user. The qualitative-likelihood paper also made clear that the current treatment assumes deterministic models $C \cdot g(x),$5; stochastic models would require additional work beyond the paper’s scope (Mitra et al., 2019).

With respect to standards and exchange formats, PyBioNetFit occupies an interesting position. PEtab was introduced as an interoperable specification of parameter-estimation problems in systems biology, using an SBML model together with TSV files for conditions, observables, measurements, parameters, and optional visualization, tied together by a YAML 1.2 problem file. The PEtab paper listed eight supporting software packages—COPASI, AMICI, pyPESTO, pyABC, Data2Dynamics, dMod, parPE, and MEIGO—but not PyBioNetFit. PEtab was presented as a problem-specification layer rather than an optimizer or simulator, and, at the time of writing, it only allowed models defined in SBML format (Schmiester et al., 2020).

In relation to PyBioNetFit, this means that direct native PEtab support was not claimed in the cited literature. A plausible implication is that PEtab could serve as an interoperability target for the SBML side of PyBioNetFit’s model support if parser or exporter functionality were implemented. Because PEtab provides validation tools, COMBINE-archive packaging, example problems, and test cases for semantics such as pre-equilibration, parameter overrides, observable transformations, and likelihood-related metadata, such an interface could function as a common exchange layer between PyBioNetFit and a wider ecosystem of simulation, optimization, and uncertainty-quantification tools. In the absence of such an implementation, however, PEtab remains an external specification standard rather than part of PyBioNetFit proper (Schmiester et al., 2020).

Taken together, the literature presents PyBioNetFit as a systems biology analysis environment specialized for difficult simulation-based inference problems: large or rule-based models, stochastic simulations, mixed qualitative and quantitative evidence, and workflows that demand explicit reproducibility. Its central methodological identity lies in the combination of standard model formats, file-based job specification, derivative-free global optimization, BPSL-based qualitative knowledge representation, and progressively more rigorous uncertainty-quantification methods.