PhysiBoSS 2.0: a sustainable integration of stochastic Boolean and agent-based modelling frameworks

The paper addresses the need to integrate intracellular signalling/regulatory dynamics with multicellular population behaviour in complex microenvironments. Boolean modelling is widely used to simulate signalling and regulatory networks at the intracellular level, while agent-based modelling (ABM) captures multicellular interactions, diffusion, mechanics, and phenotype changes. The authors present PhysiBoSS 2.0 as a redesign of the original PhysiBoSS to overcome maintainability and extensibility issues stemming from tightly coupled code, predefined phenotypes, and a lack of clear interfaces. PhysiBoSS 2.0 is implemented as a decoupled add-on that bridges PhysiCell (multiscale ABM) and MaBoSS (stochastic Boolean simulator), enabling model-agnostic integration of signalling/regulatory logic with environmental cues such as drugs and substrate transport. The purpose is to provide a sustainable, flexible framework for studying cancer and other systems by linking microenvironmental conditions to intracellular decision-making and population dynamics, supporting drug effect and synergy studies validated with experimental data.

The work builds on prior use of Boolean models for signalling/regulatory networks and precision oncology, including patient-specific logical models to predict treatment responses. At the multicellular scale, ABMs (e.g., PhysiCell) have been used to simulate diffusion, mechanics, and cellular processes in tissues and tumours. The original PhysiBoSS (v1.0) coupled an early PhysiCell version with MaBoSS to link signalling to population outcomes but suffered from design limitations that impeded maintenance and extensibility. Recent advances include PhysiCell’s tracking of internalised substrates and MaBoSS support for large models and SBML-qual. Prior studies used PhysiBoSS to explore TNF treatment dynamics and adaptive therapy, and to investigate cancer invasion modes. The paper situates PhysiBoSS 2.0 within this landscape, aiming to provide a sustainable, modular integration that leverages updated capabilities of both PhysiCell and MaBoSS, and to enable pharmacodynamics-informed simulations using dose–response data (e.g., GDSC) and established synergy frameworks (e.g., Bliss independence).

Software architecture and implementation: PhysiBoSS 2.0 is a C++ add-on to PhysiCell that provides a clean interface to MaBoSS, decoupling the two tools to allow independent updates and long-term maintainability. It uses recent PhysiCell features (e.g., tracking internalised substrates) and MaBoSS capabilities (e.g., large models, SBML-qual). Configuration is exposed via the PhysiCell XML, including per-cell-type MaBoSS model files, time step for intracellular updates, scaling, stochasticity, initial values, mutants, and parameter overrides. Custom modules standardise integration between intracellular variables and agent/environment variables. Sample projects and templates are provided; code is BSD 3-clause licensed. A nanoHUB GUI tool is available for an example model. Output processing: The PhysiCell Toolkit (PCTK) is a Python package (BSD 3-clause) that parses MultiCellDS outputs, converts them to CSV data frames, generates summary plots (e.g., time courses of alive/dead cells), and produces POV-Ray input files for 3D rendering. It can be used as a library or CLI tool. Integration patterns between scales: Continuous environmental or agent variables are mapped to Boolean inputs via transfer functions (e.g., threshold/step functions). Boolean outputs modulate cell behaviours by adjusting rates in cell cycle/death models or triggering custom rules. PhysiBoSS 2.0 reproduces prior TNF-driven spheroid results and enables alternative transport/uptake and receptor-binding mechanisms via modular submodels. TNF spheroid model replication: The authors re-implemented the 3T3 fibroblast spheroid TNF treatment experiments (from PhysiBoSS 1.0) with an extended TNF receptor binding mechanism and updated cell cycle and transport models, validating qualitative behaviours under various TNF pulse regimes (details in Supplementary Methods). Simulations used identical setups to the original work to compare outcomes. Drug-response integration and personalised models: The LNCaP prostate cancer cell line model is derived from a general prostate Boolean network, personalised via PROFILE_v2 using mutations, CNA, and RNA expression data. Drugs are mapped to target nodes (Table 2), identified via DrugBank and cross-referenced with GDSC availability. Dose–response curves per drug–cell-line pair are fitted using gdscIC50 (multilevel fixed effect model), producing normalised sigmoidal curves. Pharmacodynamics coupling: For each cell and drug, the local concentration (from nearest voxel) is used to compute a probability of target-node inhibition equal to 1 − f([X]i), where f is the normalised dose–response function and [X]i the local drug concentration. MaBoSS simulates partial inhibition by running many trajectories (e.g., 5000 Gillespie trajectories) with a specified fraction of runs having the target node inhibited vs activated, and phenotype probabilities are averaged over trajectories. In silico experiments: Spheroid simulations (~1000–1138 initial LNCaP cells) ran for 7 days in 3D domains with drugs diffusing from boundaries. Single-drug screens tested five concentrations (IC10, IC30, IC50, IC70, IC90) for six drugs targeting nodes AKT, EGFR, ERK, HSPs, MEK1/2, PI3K (see Table 2). Double-drug screens tested all pairwise combinations across concentrations. Each condition had 10 replicates; snapshots every 240 min yielded growth curves. The Growth Index (GI) was defined as GI = log2(AUC(drug)/AUC(no-drug)). Synergy was assessed using Bliss independence: Êxy = Ex + Ey − Ex·Ey and Combination Index CI = Êxy/Exy (CI<1 synergy, CI>1 antagonism). Heterogeneity analysis: To study spatial heterogeneity due to diffusion and uptake, a spheroid of radius 100 µm with boundary drug administration was simulated (10 replicates). The tumour was partitioned into concentric 50 µm-thick layers (layers 1–5), and per-layer growth curves and GIs were computed to assess positional effects on drug efficacy. Computing: Simulations ran on MareNostrum 4 (2× Intel Xeon Platinum 8160 per node, 48 CPUs total, 96 GB RAM). Individual simulations used all 48 CPUs on one node; EMEWS-based exploration used 10 nodes (three instances per node, 16 CPUs per instance). Analysis used Python 3.9 and PCTK.

- PhysiBoSS 2.0 architecture: A decoupled add-on design that cleanly interfaces PhysiCell with MaBoSS, improving maintainability, extensibility, and enabling independent upgrades. Supports XML-configured cell-type-specific Boolean models, custom modules for intracellular–agent mapping, and leverages latest PhysiCell/MaBoSS features. - Reproduction of TNF spheroid dynamics: The new implementation qualitatively reproduces PhysiBoSS 1.0 results, including resistance under continuous TNF exposure and differential effects of pulsed TNF (e.g., short pulses at 150 min frequency reduce tumour size; pulses at 600 min are ineffective). PhysiBoSS 2.0 captured more realistic exponential growth compared to a biphasic pattern previously observed, with differences attributable to updated receptor binding and cell cycle/transport models. - Single-drug screen (LNCaP): Among six tested drugs across IC10–IC90, two showed significant changes in growth behaviour: Ipatasertib (AKT inhibitor) and Pictilisib (PI3K inhibitor), with P ≤ 0.0001 (Kruskal–Wallis). For Pictilisib, Growth Index decreased monotonically from IC10 to IC90, indicating dose-dependent growth inhibition. - Double-drug synergies: Two combinations were notable—Pictilisib+Ipatasertib and Pictilisib+Luminespib—each with P ≤ 0.0001 (Kruskal–Wallis). Both exhibited stronger growth inhibition than single agents; Pictilisib+Ipatasertib achieved substantial inhibition at lower Pictilisib concentrations (IC30/IC50) compared to Pictilisib+Luminespib (which required IC70/IC90 for both drugs to reach GI ≈ −0.3). Bliss analysis revealed concentration-dependent synergy: for Pictilisib+Luminespib, synergy peaked at high concentrations with Luminespib driving the effect; for Pictilisib+Ipatasertib, synergy peaked around IC50–IC70, with higher Ipatasertib concentrations especially driving synergy. Low concentrations often yielded weak synergy or slight antagonism. - Spatial heterogeneity and shielding: Boundary drug delivery created gradients; outer layers acted as sinks, reducing drug availability centrally. In a representative case (Ipatasertib IC50 + Pictilisib IC90), inner layers (1–3) showed minimal inhibition or slight growth increases, while outer layers (4–5) showed strong inhibition and apoptosis. Layer-specific Growth Indices (median of 10 replicates) exemplified this: Layer 1: 0.05; Layer 2: 0.08; Layer 3: 0.10; Layer 4: −0.24; Layer 5: −2.57. Reduced outer-layer survival increased space for inner-layer proliferation, illustrating a competition between drug penetration and mechanical constraints. - Tooling: PCTK streamlined processing of MultiCellDS outputs, generation of summary plots, and 3D renders, facilitating reproducible analysis.

PhysiBoSS 2.0 addresses the central challenge of linking intracellular signalling decisions with multicellular dynamics in heterogeneous microenvironments through a sustainable, modular framework. By decoupling MaBoSS and PhysiCell and standardising interfaces, the tool enables reuse of arbitrary Boolean models and facilitates integration of environmental cues (e.g., diffusing drugs, receptor binding) into cell decision-making. The replication of TNF treatment dynamics validates the approach and demonstrates that updated transport/receptor-binding implementations can improve realism without compromising qualitative behaviours. The pharmacodynamics coupling via dose–response curves allows the study of drug effects and synergies in a spatial context, revealing how concentration-dependent interactions (e.g., Bliss synergy peaks at specific IC ranges) and spatial heterogeneity (outer-layer shielding and inner-layer proliferation) shape population outcomes—phenomena consistent with experimental observations of context-dependent synergy/antagonism and microenvironment-driven therapy resistance. These findings underscore the importance of multiscale, spatially explicit models when evaluating treatments and combinations and illustrate PhysiBoSS 2.0’s utility for hypothesis generation and in silico screening aligned with personalised logical models. Overall, the results demonstrate that PhysiBoSS 2.0 can reproduce prior benchmarks, extend to mechanistic or empirically driven pharmacodynamics, and provide insights into heterogeneous responses, positioning it as a foundation for scalable, personalised digital twins and virtual clinical trial simulations.

PhysiBoSS 2.0 delivers a sustainable, model-agnostic integration of stochastic Boolean and agent-based frameworks by decoupling MaBoSS from PhysiCell and exposing flexible, modular interfaces configurable via XML. It reproduces established TNF-driven spheroid behaviours, supports mechanistic or empirical drug–target coupling, and enables in silico drug screens revealing significant single-agent effects and concentration-dependent synergies, as well as spatial heterogeneity driven by diffusion and uptake. The work contributes an extensible software ecosystem, including sample projects and the PCTK analysis toolkit, to facilitate model development and reproducible analysis. Future directions include incorporating extracellular matrix interactions, vasculature and vascularisation, complex 3D patient-derived architectures (spatial omics), and explicit PK/PD modules now available in PhysiCell, to advance towards realistic, personalised digital twins for therapy design and evaluation.

- Drug–target interaction modelling often relies on fitted dose–response curves rather than detailed binding kinetics, reflecting limited mechanistic knowledge for many compounds. - Drug administration in simulations is from domain boundaries, not from explicit vasculature, leading to simplified spatial distributions; extracellular matrix and vessel structures are not yet included. - Synergy assessment uses the Bliss independence framework; other interaction models (e.g., Loewe, HSA) are not explored here. - Validation of drug effects is limited; while some single-drug treatments were compared with experimental cell survival data, broader experimental validation and multiple cell lines remain for future work. - Results are shown primarily for the LNCaP cell line and specific drugs; generalisability to other contexts requires additional personalisation and testing.