BSim

: Agent-based models (ABMs) provide a number of advantages relative to traditional continuum modeling approaches, permitting incorporation of great detail and realism into simulations, allowing in silico tracking of single-cell behaviors and correlation of these with emergent e ﬀ ects at the macroscopic level. In this study we present BSim 2.0, a radically new version of BSim, a computational ABM framework for modeling dynamics of bacteria in typical experimental environments including micro ﬂ uidic chemostats. This is facilitated through the implementation of new methods including cells with capsular geometry that are able to physically and chemically interact with one another, a realistic model of cellular growth, a delay di ﬀ erential equation solver, and realistic environmental geometries.

T he detailed simulation of cell populations is increasingly sought after in Synthetic Biology for the in silico design, testing and validation of new biological circuits and devices. 1,2gent-based models (ABMs) offer the unrivalled ability to capture complex dynamics of cell populations at different scales, by defining and taking into account specific rules for the agents and environmental characteristics.Many software packages have been developed to implement ABMs for Systems and Synthetic Biology 3 (Figure 1a).Specific tools (BSim, 4 BNSim, 5 gro 6,7 ) have been implemented to deal with the simulation of generic intracellular processes in bacteria while taking into account cells' motility as well as spatial distribution, and the diffusion of signaling molecules and chemicals in the growth environment (Figure 1a).
However, existing software packages lack some fundamental features that are crucial to realistically take into account physical interactions within the agents' population and between the cells and the simulation environment.Most notably, no existing agent-based modeling tool for bacterial populations offers the possibility of defining a realistic 3D morphology for the cells, and also of setting the physical parameters (e.g., dimensions, boundary conditions) of realistic 3D simulation domains similar to those used in experimental set-ups based on microfluidics. 8n order to overcome these drawbacks (Figure 1a), we present here BSim 2.0, a radically new ABM tool for bacterial populations which embeds all the features of BSim 4 but significantly extends its capabilities by adding the following features: (a) intracellular dynamics described in terms of delayed differential equations (DDEs) as well as ordinary differential equations (ODEs), (b) realistic 3D cell morphology implemented in terms of growth, division and mutual cell interaction, and (c) realistic physical parameters of the 3D simulation environment including microfluidic chemostat characteristics if needed (Figure 1b).

■ METHODS
Implementation.The architecture of the open source BSim 2.0 framework is written in Java, and is designed to be modular providing an extensible base of common features for modeling typical bacterial behaviors (Figure 1b; Figure S1).
Cellular Growth.Cells are implemented as individual Java objects, with per-cell methods that allow for manipulation of properties according to physical interactions and integration of internal gene regulatory network dynamics.Capsular cell growth is modeled using a per-cell ODE model of rod elongation over time as in the literature, 9 with division occurring once the mother cell has passed a set constant threshold length.Cell length is integrated according to an Euler scheme in synchrony with the global simulation time step.Initial cell length and maximal division threshold can be parametrized appropriately to correspond to the biological system of interest, furthermore cell radius, minimal and maximal length and elongation rate can be similarly chosen depending on the features of the species being simulated.Upon division, daughter cells' positions and lengths are perturbed by a small amount (randomly chosen, with an amplitude between ±5% of their length) at the location of division in order to break axial symmetry. 9ll−Cell Mechanics.Cells are represented by capsular volumes, typical of E. coli cells that are commonly used in microfluidic experiments.Cell−cell interaction forces are computed using a semirigid-body approach as employed in the literature 9−12 among others.Possible nonphysical intersections between the cells are resolved by computing, for each pair of intersecting cells, an overlap-dependent volume exclusion force.13 Positions of all cells are then simultaneously relaxed according to the computed forces in order to minimize all overlaps of all intersections.10 Details of CPU usage under a range of conditions are provided in Supporting Information and Figure S2.
Environmental Geometry.Geometry is represented internally using boxes or triangular-element meshes, e.g., a microfluidic chemostat as illustrated in Figure 1b.Meshes can be imported into the simulation domain via the Wavefront object format (this format is available as an export option in all common 3D modeling and Computer-Aided Design packages).Cell-geometry interactions are implemented using a similar rigid-body approach to the cell−cell mechanics, using methods defined in the literature. 14,15At each time step, potential intersections between cells and geometry are computed.In the case of overlap, an exclusion force proportional to the interpenetration distance is computed and applied to the cells at the position constraint relaxation stage detailed in Cell−Cell Mechanics.
GRN Implementation and Simulation.BSim 2.0 allows one to specify realistic population-level heterogeneity in system parameters, as well as in interactions between intracellular and extracellular mechanics.Differently from any other agent-based model for cell populations, gene regulatory networks can be modeled in BSim 2.0 by means of delayed differential equations.

■ CASE STUDY
In order to illustrate the features of BSim 2.0, we applied it to simulate the challenging scenario of the oscillating multicellular consortium described in the literature 2 (Figure S3).This is a particularly relevant case study in Synthetic Biology as the accurate simulation of cellular consortia is becoming of increasing interest in the community. 16pecifically, the system consists of a consortium of two populations of engineered cells, activator cells and repressor cells; these are coupled together through two spatially uniform external chemical concentrations.the original system was modeled by 16 coupled delay differential equations.the varied structure of experimentally relevant microfluidic devices naturally lends itself to implementation as a 3d geometry; the effects of experimental boundary conditions such as external flow rate as well as microchemostat topology can significantly affect the outcome of synthetic circuits whose internal dynamics are coupled to external quorum signaling. 17he system dynamics were highly sensitive to the signaling chemical parameters (Figure S4) and boundary conditions (Figure S5).Furthermore, where the two populations were increasingly spatially separated from one another (Figure S6) the system did not function as desired; indeed this remained true with growth-induced mixing included in the model (Figure S7).However, in a well-mixed scenario (Figure S8), corresponding to a hypothetical "best-case" experimental condition, the spatially extended model predicted behavior equivalent to the real biological implementation.With the realistic growth and division dynamics embedded in BSim 2.0, our validation showed realistic behavior qualitatively similar to that of the real system (Movie S1, and Figure 1c, confirming the ability of BSim 2.0 to capture and predict key features of the experiments reported in the literature. 2

■ CONCLUSIONS
We have presented BSim 2.0, a new unified platform for in silico studies of cellular populations including multicellular consortia.
The key features distinguishing BSim 2.0 from any other related software package in the literature are that it allows: 1. DDE simulation that is individual to each agent and can be coupled to an arbitrary number of external diffusion/ degradation PDEs 2. Capsular cell geometry with logistic growth model coupled to division; cell contents are distributed to daughters upon division.3. Physical cell−cell interactions between capsules; physical interactions between capsules and the simulation domain.4. A constrained 3D domain with mixed boundary conditions, that can for example represent different commonly used microfluidic chemostats.The results of our case study show that BSim 2.0 is a valuable tool for the design of synthetic regulatory networks acting across multicellular consortia.
Future directions for development that are currently being explored include, in the short-to-medium term, a parallelised cell−cell collision dynamics solver, a graphical user interface, and further refinement of the growth model.The new features, and future extensions, presented here could be utilized to predict the effect of dynamically changing nutrient conditions, potentially leading to changes in cell shape and packing in a variety of 3D environments, on desired population-level behavior.Such models that couple growth conditions to intracellular gene networks, via models of cell metabolism, could permit the simulation of multispecies consortia (directly incorporating the effects of different cell shapes and growth rates) and to investigate how collective dynamics arise from intercellular interactions.
The source code and documentation for BSim 2.0, with examples including those presented here, is available at https:// github.com/bsim-bristol/bsim.

Figure 1 .
Figure 1.Implementation of BSim 2.0 features, their context with respect to the state-of-the-art, and their application to a recently characterized multicellular system.(a) A comparison of agent-based model tool features, tabulated to emphasize the new functionality that we have implemented in BSim 2.0.(b) New features of BSim 2.0.(c) Results of the case study (system described in the literature 2 ), with parameters modified to improve intercell signaling (see Case Study).Simulation performed with realistic cell morphology, growth and division with cells initialized at randomly chosen positions, in a 100 × 85 × 1 μm 3 rectangular microfluidic chemostat.A rendered animation of the simulation and corresponding quantification can be seen in Movie S1. (i) Mean population output (solid lines) is plotted, for Activator cells (blue) and Repressor cells (green); filled regions indicate mean ± SD; the red arrows indicate time points of simulation snapshots rendered in panel (iii), below.Black dashed line indicates the proportion of activator cells present in the chemostat over time.(ii) Single-cell output is plotted from the 50 longest-lived cells in the simulation (see Supporting Information for details); the normalized amplitudes of fluorescence outputs measured from Activator cells (blue) and Repressor cells (green) are indicated; time series are plotted for individual cells (rows), with simulation time increasing left-to-right.(iii) Simulation snapshots rendered at time points indicated by red arrows in (i).Individual cells are indicated by capsules, colored blue (Activator cells) and green (Repressor cells), with the intensity of each color corresponding directly to the level of output fluorescence.Simulations were initialized with 500 cells of each type; the filled chamber (snapshots) contains on average 2150 cells.