Novel Ground-Up 3D Multicellular Simulators for Synthetic Biology CAD Integrating Stochastic Gillespie Simulations Benchmarked with Topologically Variable SBML Models

Abstract

The elevation of Synthetic Biology from single cells to multicellular simulations would be a significant scale-up. The spatiotemporal behavior of cellular populations has the potential to be prototyped in silico for computer assisted design through ergonomic interfaces. Such a platform would have great practical potential across medicine, industry, research, education and accessible archiving in bioinformatics. Existing Synthetic Biology CAD systems are considered limited regarding population level behavior, and this work explored the in silico challenges posed from biological and computational perspectives. Retaining the connection to Synthetic Biology CAD, an extension of the Infobiotics Workbench Suite was considered, with potential for the integration of genetic regulatory models and/or chemical reaction networks through Next Generation Stochastic Simulator (NGSS) Gillespie algorithms. These were executed using SBML models generated by in-house SBML-Constructor over numerous topologies and benchmarked in association with multicellular simulation layers. Regarding multicellularity, two ground-up multicellular solutions were developed, including the use of Unreal Engine 4 contrasted with CPU multithreading and Blender visualization, resulting in a comparison of real-time versus batch-processed simulations. In conclusion, high-performance computing and client-server architectures could be considered for future works, along with the inclusion of numerous biologically and physically informed features, whilst still pursuing ergonomic solutions.

Keywords: CAD; SBML; automation; biophysics; chemical reaction networks; multicellular simulation; stochastic Gillespie; synthetic biology; systems biology; unreal engine.

Figure 1

A novel workflow was developed involving utilities, backend components/engines, as well as two multicellular spatiotemporal methodologies featuring different technical approaches. Modularity could be considered across the workflow and/or at the algorithmic level.

Figure 2

UnrealMulticell3D (UM3D) dimensionally constrained monolayer colony of bacillus cells. See also supplementary Video S1 for more UM3D capabilities and dynamic visuals.

Figure 3

Control Flow featuring NGSS integration into UM3D using a custom C++ blueprint function library. Blueprinting was also used for other cellular computations/functions. Top Left: The construction script defaults Tick Gate to 2, allowing the biochemical model to be fired. Tick Gate is set to 0. Top Right: NGSS completion is awaited and cellular functions besides physics are arrested. Tick Gate is set to 1 when the output file is detected. Lower Left: NGSS output file is cleared, and deletion is confirmed. Tick Gate is set to 2. Lower Right: The Boolean value of “FiredNGSS” is true so the “Run NGSS” node is bypassed allowing for growth, division and other phenotypic computations on the agent cell. In this prototype the biochemical model was fired only once per cell cycle, hence daughter cells would undergo the same process.

Figure 4

SynthMeshBuilder’s diverse morphology generation of mesh-based cell networks (upper left), highly scalable one million cell colonies (upper right), parallelized “stochastic chain extension” reminiscent of staphylococcus clusters [55] used for benchmarking (lower left) and an alternative on-lattice algorithm with random update order (lower right), as visualized with Blender (see supplementary Video S4).

Figure 5

600,000 “Pseudo Off-Lattice” cells clustering around an obstacle as visualized in Blender via .OBJ files output by SynthMeshBuilder. The stochastic expansion of the colony is evident from the morphology. Collision surfaces of such a nature could also be used to trigger contact signaling calculations.

Figure 6

For the low propensity models generated by SBML-Constructor, the Tau Leaping NGSS algorithm performed best using the multi-enz models. Results were generated using NGSS-Invoker. Optimized Direct Method is not shown as it was suggested retrospectively to initial benchmarking and only recommended for pathways of 2 and 3 reactions. Likely it was immediately discarded due to inferior performance and practical time constraints. SSAPredict was shown to be wrong with its prediction of performance for the Logarithmic Direct Method, i.e., Tau Leaping was 2.29 times faster than the Logarithmic Direct Method with the 48 reactions model.

Figure 7

Towards processor saturation caused by increasing “concurrent” cell target computations of NGSS activations, time “per cell” for completion of the reaction model decreased. Beyond saturation and towards cell completions, hence more .CSV file checking, time per cell increased gradually (note the uneven unit distributions on the x-axis). The above data used the NGSS settings determined in the previous section. Minimizing IO in the future would be ideal. These results were generated using the in-house NGSS-Invoker utility coupled with NGSS executions.

Figure 8

Benchmarked time performance data subset of UnrealMulticell3D on a time per reaction step basis. Note that greater cell target populations ensured processor saturation, explaining the peaks, with performance consistency upon saturation given unchanging reactions per cell cycle (RPC). There was polynomial scaling with RPC beyond the NGSS nuance up to ~7 RPC.

Figure 9

A performance profile was run for UM3D using the Unreal Insights utility with the camera turned off. FPhysScene_ProcessPhysScene (blue), processed on the CPU, required greater processing as the cell population expanded. The plots behind the blue plot can be ignored for this discussion. The lower track illustrated limited GPU activity involving intermittent SlateUI updates. PhysX calculations appear to be restricted to the CPU by UE4. The impact of physics processing on the CPU might have been missed at the relatively low cell populations tested (up to 4096 over a span of ~20 s) because of its fairly modest initial CPU usage.

Figure 10

Unlike when the camera was off (see Figure 9), the GPU had a vast assortment of functions to fulfill when the camera was turned on (GPU track, top). However, none of the functions appeared to perform physics calculations; whereas on the CPU GameThread track, FPhysScene_ProcessPhysScene continued to be processed.

Figure 11

Performance enhancement of UnrealMulticell3D without NGSS, due primarily to a move towards batch-processing and reduced rendering costs. Mesh simplification enhanced performance to generate 4096 cells by 1.33 times. Combined mesh simplification and camera redirection enhanced generation of 4096 cells by 2.63 times.

Figure 12

Despite enhancements to the UnrealMulticell3D base engine depicted in Figure 10, performance with NGSS was identical, demonstrating that NGSS was the limiting factor with respect to performance. The above depicts an average time across all benchmarked cell target populations for each reaction per cycle pathway model, used here as an indicator of the average performance.

Figure 13

Mesh generation by SynthMeshBuilder to a 252,000 vertex population was temporally enhanced up to CPU logical core saturation (12 threads) when parallelized with various numbers of launched threads demonstrating successful multithreading on the CPU architecture.

Figure 14

Both coupled with NGSS, SynthMeshBuilder scaled faster than UnrealMulticell3D due to its streamlined algorithmics, but with a similar trend. The target populations differ above.

Figure 15

The base multicellular simulation layers could generate 1000-cell populations within seconds (UM3D) or fractions of a second (SMB). However, adding a moderately sized reaction network (32 reaction steps) with parallel NGSS processes (one for each cell cycle) resulted in a drastically more time consuming performance profile operating on the order of several minutes to complete. Because NGSS could use as much as a single core (two threads) of processing power for each activation, by the time only 6 cells that reached the processor could be saturated, giving an overall linear scaling as NGSS processes were queued for completion on the rapidly saturated single processor. The same linear scaling would be expected with HPC but with a shallower gradient, at least once the HPC cores were fully saturated. For example, with the 32 reaction network UM3D proliferated 299 times slower to 1024 cells than without it (1598.20 s vs. 5.34 s).

Figure 16

SynthMeshBuilder proved much more scalable without NGSS than UnrealMulticell3D on the personal computing system. By contrast, HPC solutions from the literature could process millions [11] or even tens of millions [10,36] of cells, with thousands [23] or hundreds of thousands [11] reported on modest hardware.