Role of intercellular interactions on single cell and population level responses: considerations for multicellular bioreporter design

Introduction: Bioreporters are genetically engineered cells that produce detectable responses in the presence of specific analytes, providing a cheap, mass-producible, and accurate method of analyte detection. Most research focuses on the single cell-level, where all engineering is concentrated on the interactions within a single cell. However, intercellular communication is a well-known natural phenomenon that has been associated with sensitive responses to certain chemical stimuli, yet incorporation of intercellular communication into bioreporter design is exceedingly rare and the effect of intercellular signaling on single-cellular and population level responses has not been explicitly characterized before.

Methods: In this work, a multicellular simulator implementing the Gillespie algorithm and compatible with Virtual Cell-designed networks is created and used to demonstrate nuances to multicellular stochastic simulations. The algorithm was used to simulate multiple cells in a reaction network in which a self-promoting and membrane permeable transcription factor also induces production of a cell-bound reporter protein. A proof-of-concept bioreporter that responded to an environmental analyte while leveraging intercellular interactions for signal production was also designed and simulated for 50 s of simulation time. Simulated systems with multiple cells were compared to single-cellular simulations.

Results: Simulations for 20 s of simulated time show that while final reporter protein count per cell decreased as cell count increased, aggregate final reporter protein number across all cells significantly increased. Interestingly, 50 s simulations show final reporter protein count per cell significantly increasing as the number of cells increases. Greater number of bioreporter cells resulted in significantly greater amounts of signal protein produced in response to the same amount of starting analyte both on the population level and on the individual cellular level.

Discussion: The results show significant differences between the results of multicellular and single cellular simulations, which demonstrates the importance of simulating multiple cells to obtain nuanced results. The amplification of signal protein produced by an increasing number of simulated bioreporter cells indicates great potential for multicellular bioreporter designs to amplify response magnitude and sensitivity on the individual cellular level. These results warrant further research into the application of different simulation algorithms and multicellular bioreporter design and modeling.

GRAPHICAL ABSTRACT

Introduction

Cellular simulation is mainly concerned with accurately modeling the rates at which intracellular processes, such as transcription and translation, occur. Accurate simulation of these reactions over time allows for accurate understanding of cell behavior without costly lab expenses, saving long term costs. These reactions that occur on subcellular scales cannot be accurately modeled by differential equations due to extremely low concentrations of reactants. Algorithms that incorporate the stochastic effects dominant at such low concentrations are required. Indeed, such stochastic effects explain behavior of transcriptional networks not captured by deterministic models (Hahl and Kremling, 2016) and help elucidate biological design features to manage stochastic noise that would other not be appreciated (Boettiger, 2013). In a comparative study by Johnson et al., in 2021, they found that in biochemical networks with feedback loops produced major differences in results when simulated with stochastic dynamics when compared to deterministic simulations (Johnson et al., 2021).

A common method of simulating these processes is using the Gillespie algorithm (Gillespie, 1977). In a system with many possible reactions, the Gillespie algorithm first computes the instantaneous reaction rates of every reaction. The algorithm then chooses a reaction to occur next based on its relative reaction rate. Finally, the time it takes for that reaction to occur is drawn from an exponential distribution with λ equal to the sum of the instantaneous reaction rates of all reactions. The Gillespie algorithm is the most fundamental and exact method to predict the future evolution of a biochemical system while considering stochastic effects (Warne et al., 2019). The Gillespie algorithm, or an algorithm mathematically equivalent to it, is very commonly used in many stochastic simulation tools, most of which are designed for single-cell simulations. Examples include Virtual Cell (Resasco et al., 2011) and COPASI (Hoops et al., 2006). While these simulators can accurately describe single cell behavior, they are not equipped to simulate multicellular systems at scale. Their inability to model multicellular systems prevent possible effects of intercellular interactions from being understood. However, the effect of intercellular interactions in stochastic simulations has not been explicitly characterized before, so the knowledge lost from the inability to model multicellular systems is not well understood either. This project aims to fill this gap in knowledge by explicitly comparing the results of a single cellular simulation with the results of a multicellular simulation, both on the single cellular and population level response. The research question for this project is: do intercellular interactions significantly affect single cellular and population level responses? The hypothesis for this project is that if intercellular interactions are accounted for, then the results produced by multicellular simulations with intercellular interactions will be significantly different from the results produced by a unicellular simulation that does not account for intercellular signaling, both on the individual cellular level and on the population level. To achieve this, the results of multicellular simulations with intercellular interactions will be compared to the results of an unicellular simulation. The exact method of quantifying and comparing results of different simulations is described in the Methodology section.

A better understanding of the role of intercellular signaling will be valuable in multiple areas within biology. For example, construction of genetically engineered whole-cell biosensors, which are meant to produce detectable signals in response to the presence of certain environmental analytes, is a current area of active research (Saltepe et al., 2017). These biosensors are constructed with consideration only within a single cell and the processes within it. These biosensors can often be modeled in single-cellular cell simulators like Virtual Cell and validated for responses to their intended stimuli before being expressed in a lab, saving cost. However, a better understanding of intercellular interactions in cellular responses could allow for the construction and computational validation of biosensors that leverage intercellular signaling, expanding the possible areas for biosensor optimization. Intercellular signaling among bacteria has long been known and studied in nature, especially in the context of propagating responses to a stimulus to an entire colony of bacteria (Antunes and Ferreira, 2009). Adapting this natural process would likely yield better biosensor designs than designs that only focus on intracellular interactions due to the proven effectiveness of intercellular interactions in nature. Coordinated luminescence is already naturally achieved through quorum sensing (Li and Nair, 2012). Biosensors are cheap, onsite detection methods for soil and water contaminants (Zeng et al., 2021), so the construction of more sensitive biosensor designs could further improve the usefulness of this detection method. Cheaper detection methods could allow for higher throughput and more frequent testing for dangerous contaminants, potentially saving lives by finding sources of contamination early before the effects of these contaminants manifest severely. By stochastic simulations to model intercellular communication, a secondary research hypothesis explored in this project is as follows: a bioreporter design incorporating intercellular signaling will produce a greater signal in response to a selected analyte.

One of the most popular cellular simulation tools is Virtual Cell, created by Resasco et al., in 2011. Virtual Cell has a graphical user interface for designing reaction networks and uses a modified but mathematically equivalent version of the Gillespie algorithm to model the evolution of all species present in the user-created reaction network according to the parameters set by the user. It also allows users to define different compartments with different volumes to separate interactions within a cell and those outside. The user-friendly interface and ease of use makes Virtual Cell very useful for simulating reactions occurring in a small number of compartments, such as the reactions in a single cell and an extracellular space. However, Virtual Cell does not have any built-in capability to extend simulations to multiple cells easily. It is possible to manually create multiple compartments with duplicated reaction rules to represent multiple cells, however the process is tedious and not feasible for simulating hundreds of cells, despite many real-life situations containing hundreds or thousands of cells. This project aims to extend simulation capabilities using the Gillespie algorithm to a large number of cells while also accounting for shared chemical species between cells to represent intercellular interactions, then comparing the changes in simulation results that result from multicellular simulations that would otherwise not occur in the single cellular simulations that Virtual Cell is most used for.

Previous research on simulating intercellular interactions have taken diverse approaches. Weber and Buceta in 2013 simulated a biological “switch” comprising 1000 cells and accounting for the production and effects of intercellular quorum sensing molecules produced by each of the cells. Each cell contains the genetic components for a positive feedback loop where extracellular quorum sensing molecules diffuse into cells and stimulate the production of a reporter protein and more quorum sensing molecules. These interactions are very similar to the ones that will be used in this project, as described in the Methodology section. Weber and Buceta individually plotted the concentrations of select chemical species in each cell simulated and observed that significant variations in reporter protein concentration between cells existed. Given enough time all cells plateaued in reporter protein production. They attributed their observations to the impact of intercellular signaling, however they did not explicitly compare the simulation results of the population of cells with a simulation of only a single cell. They also did not vary the amount of cells they simulated to observe the effect of the number of cells in a population on the behavior of the individual cells or the population of cells as a whole. Therefore, the actual impact of intercellular interactions on their simulation results is not clearly defined. This project seeks to explicitly compare the difference in results caused by accounting for intercellular signaling using a similar reaction network as the one used by Weber and Buceta to determine if the presence of multiple cells and their intercellular interactions contributed significantly to the behavior they observed. This project also creates a general framework to simulate any multicellular system based off of one designed in the Virtual Cell graphic user interface. Finally, this project extends beyond the work of Weber and Buceta by analyzing generalized constructions of multicellular signaling and bioreporters for analytes.

Duso and Zechner in 2020 studied multicellular simulations from the perspective of a multi-compartmental system. Intercellular interactions were represented as cells directly acting on each other. The authors developed differential equations to accurately capture summary statistics of trajectories generated by stochastic simulation. While their work focuses solely on the time evolution of multicellular systems, this work compares multicellular systems with single cellular ones. Additionally, this work simulates the mechanism of intercellular communication, i.e., the diffusion of a signaling protein, instead of abstracting interactions as probabilistic events.

Instead of explicit simulation of multiple cells, Zhou et al., in 2005 modeled intracellular and extracellular noise of a synthetic quorum sensing circuit from the perspective of a single cell using Langevin equations. In this work, we explicitly model multiple cells sharing the same environment to determine similarities and differences between their findings and ours. Additionally, we vary the cell population size to examine its effect on cellular responses.

Smith and Grima in 2018 use stochastic simulation to compare independent, insulated cells with tissue-bound cells able to communicate with each other. Their independent cells are analogous to a single cell simulation, and their tissue-bound cells are analogous to a multicellular system. However, they do not include a shared extracellular space, and this work includes an additional simulation of a bioreporter design.

This work focuses on the explicit comparison of a given biochemical reaction system between a single cellular system and a multicellular system to establish a significant difference. Then, the findings are extended to a proof-of-concept model of a bioreporter design that leverages multicellular interactions to increase signal strength. In the process, a general framework to simulate multicellular systems based off of Virtual Cell-designed single-cellular systems is created.

Methodology

Algorithm overview

In order to simulate multiple cells, a custom simulator first had to be created because none currently exist. A diagram of how multicellular simulations are constructed can be seen in Figure 1. The simulator implements the Gillespie Stochastic Simulation Algorithm, which probabilistically simulates individual reactions based on their relative rates as well as the time between reactions (Gillespie, 1977). The rate, or activity, of a given reaction is given by the elementary rate equation

$a_j=k_j*N_{Reactant1}*N_{Reactant2}*\dots N_{Reactantn}$

where N_{Reactant n} is the number of molecules of Reactant n that are present in the system, and k_j is the reaction rate corresponding to reaction j for 1 ≤ j ≤ m, where m is the total number of reactions defined for a system.

Figure 1

Figure 1. Diagram of how unicellular simulations are converted into multicellular systems. “VCell” is the abbreviation for the Virtual Cell software. VCML is the file extension used by VCell to represent simulation environments.

Simulation time is incremented by discrete time step τ seconds, which is sampled from the exponential distribution with a probability density function of

$p\left(\lambda \right)=\lambda e^{-\lambda x}\text{\hspace{0.17em}and\hspace{0.17em}}\lambda ={\displaystyle \sum _{j=1}^m}{a}_j$

Before each time increment, a reaction is chosen to occur by the algorithm by finding a number k 1 ≤ k ≤ m that satisfies the inequality

${\displaystyle \sum _{j=1}^k}{a}_j<\sim U\left(0,{\displaystyle \sum _{j=1}^m}{a}_j\right)<{\displaystyle \sum _{j=1}^{k+1}}{a}_j$

Where ∼ U(a, b) represents a sample from the Uniform Distribution with a lower bound of a and an upper bound of b. This essentially selects a reaction at random, with each reaction weighted by its activity. The kth reaction is then simulated, and the number of molecules of reactants and products of the reaction are appropriately incremented or decremented. The process is repeated until simulation time exceeds the allotted amount of time.

The algorithm was implemented in the Python (3.11.8) programming language and uses the Numba (v0.59.0) and Numpy (v1.26.0) libraries for runtime optimization. All code was written in a Jupyter notebook and run locally on a Macbook M2. Multicellularity was achieved by duplicating the intracellular cellular environment so that they all shared the same species in the extracellular environment. Thus, we reproduce a homogenous cell population within a shared environment.

All code and supporting files used to generate the results in this work are hosted on GitHub and may be accessed in Supplementary Appendix A.

Description of simulated cells and reaction networks

In this work we simulate an idealized, theoretical bacterial cell that is defined only by the biochemical reaction pathways of interest. As bioreporter pathways should be designed to be biologically orthogonal to native pathways, this simplification should be representative of real circumstances (Costello and Badran, 2021). Similarly, the environment that the cells are placed in is also defined only with the species of interest that interact with the simulated biochemical pathway. In future works, more chemical species and interactions can easily be defined within cells and the shared environment.

Any reaction that can be characterized with a set of reactants and products can be modeled with this simulator. In this work we model protein-DNA reversible binding, translation and transcription, protein degradation, ligand-protein reactions, and diffusion through membranes. Given the generality of what can be simulated, more types of reactions can be defined and simulated, e.g., small molecule-small molecule interactions, protein-protein interactions, RNA-DNA interactions, and more.

The reaction network used in the project was then designed in Virtual Cell (v7.7.0.9) and can be viewed in Figure 2.

Figure 2

Figure 2. Virtual Cell network representation of experimental reaction network. “TF” represents the reporter protein inside the cell and “TF_out” represents the protein in the extracellular space. Arrows can represent forward and reverse reactions, as is the case with TF_out diffusing into the cytoplasm of each cell and TF diffusing back out into extracellular space. null represents the products of degradation. Image produced by researcher.

In the reaction network, the expression of reporter protein “P” is increased by the concentration of a transcription factor “TF”, which also induces expression of itself. TF binds reversibly to cellular DNA to create a complex that creates mRNA for P and TF that are translated into the proteins. TF is produced within a cell and diffuses freely through the cell membrane, acting as the main mechanism of intercellular interaction. P is created in the cell but cannot diffuse out and is used as the quantification metric to test the hypothesis. This network is very similar to the one described by Weber and Buceta in 2013 and has very clear intercellular interactions, making it a good model for understanding the impact of intercellular interactions on cell responses. It assumes a well-mixed system and no changes in cell population for the duration of the simulation, which is consistent with the real-world use-case of bulk signal production by bioreporter cell populations. It is also consistent with the assumptions in related works (Weber and Buceta, 2013; Duso and Zechner, 2020).

Reaction stoichiometries, rate laws, rate constants, and starting concentrations were then defined. Rate constants were chosen as reasonable estimates rather than modeling an existing biological system. This choice was made to maintain a simple system, and because this work focuses solely on the relative differences between multicellular and unicellular dynamics, thus making the exact parameter values less critical to answering research questions. The rate at which TF diffuses out of the cell was purposefully set to be 100 times faster than the rate at which TF diffuses into the cell. This ratio is the result of a cell that has 1/100 of the volume of its environment. The starting amount of TF_out was set to be 3000 molecules and the amount of DNA was set to be 100 molecules. The starting amount of extracellular TF was estimated to be small enough so that cellular production of TF would be relatively significant compared to the amount of already existing TF, but large enough to have an impact on cellular response. The starting amounts of all other species were set to 0 molecules.

The simulator was then run with these conditions, simulating 1, 10, 25, 50, 75, and 100 cells, with the number of cells acting as the independent variable and the results of the single cellular simulation being the control against which the multicellular simulations will be compared. For each number of cells, two simulations were run for 20 s and 50 s for 50 trials each. 20 s and 50 s were chosen as the simulation time limits primarily to expedite experimentation, as the values provided distinct results while allowing sufficient observation of the evolution of the simulated systems.

The amount of reporter protein P was recorded for every cell after 20 s and 50 s of simulation time. To find the average amount of P present in a cell at 20 and 50 s for a given number of cells simulated, the amount of P in each cell was averaged across all cells and all trials before being recorded. The aggregate amount of P present in all cells at 20 and 50 s for a given number of cells simulated was similarly found by summing the amount of P present in each cell across the number of simulated cells and averaged across trials before being recorded. Both the average amount of reporter protein present in each cell and the total amount of reporter protein present across all cells were the dependent variables used to test the hypothesis.

2-tailed t-tests were performed to determine if multicellular simulations resulted in levels of P protein that were significantly different from the results produced by the unicellular simulation using a significance level of 0.001. A parametric t-test was determined to be appropriate due to the high number of replicates per experimental condition. A threshold of 0.001 was chosen instead of 0.05 to be more stringent on significance given the inherently noisy data being analyzed. These data points will be analyzed in the Data and Results section as well as in the Discussion section. Simulation times have also been compiled and reported in the Data and Results section.

For a proof of concept, a more complex multicellular bioreporter network was designed and simulated with different numbers of cells for 50 s to discern if a system with multiple shared chemical species across cells would behave similarly to a system with only one shared species, as in the experimental reaction network. A shared “signal” molecule is produced by the cells which is induced by the “analyte” and in turn induces the production of a “reporter” molecule. The reporter molecule is used to determine the signal and response strength of the detector cells. The reaction network is displayed as a Virtual Cell network in Figure 2.

Data and Results

20-Second simulations

After 20 s of simulation time, the amount of reporter protein molecules present in each cell was averaged across trials and across all cells simulated. Those values are represented in this table. A 2-tailed t-test was performed on the results of simulating more than 1 cell and the results of a unicellular simulation to establish a statistically significant difference in the results of multicellular simulations compared to unicellular simulations. Only the 10 cell simulation produced results with a p value greater than 0.001, with the 25, 50, 75, and 100 cell simulations producing results with a p value less than 0.001. This means the null hypothesis that multicellular simulations produce insignificantly different results than single cellular simulations can be rejected for all the multicellular simulations that simulate more than 25 cells that were tested. This strongly indicates that multicellular simulations produce statistically significantly different results than single cellular simulations at the single cell level, at least at 20 s of simulation time and for simulations with more than 25 cells. A scatter plot of the data in Table 1 is represented in Figure 3.

Table 1

Table 1. Summary statistics of average reporter protein count per cell after 20 s of simulation time.

Figure 3

Figure 3. multicellular bioreporter design with two shared chemical species, “signal” and “analyte”. Simulations were run for 50 s of simulation time. Network and reaction rate constants can be found in the public Virtual Cell model called “Bioreporter-POC”. Image produced by researcher.

The data represented in the graph demonstrates a negative, nonlinear relationship between the number of cells simulated and the amount of reporter protein present at the end of simulation time. This trend could be rationalized by considering how additional cells will dilute the intercellular signal amongst more cells, leading to less signal present in each individual cell and therefore eliciting less of a response from each cell.

At 20 s into simulation time, the amount of reporter protein molecules present in each cell was summed across all cells simulated, then averaged across trials. Those values are represented in this Table 2 tailed t-tests were similarly performed to establish a statistically significant difference in the results of multicellular simulations compared to unicellular simulations. All t-tests resulted in p-values that were multiple orders of magnitude less than the required threshold for significance, p < 0.001, which means the null hypothesis that multicellular simulations produce insignificantly different results than single cellular simulations can be rejected for all multicellular simulations tested. This strongly indicates that multicellular simulations produce statistically significantly different results than single cellular simulations at the population level, at least at 20 s of simulation time. A scatter plot of the data in Table 2 is represented in Figure 4.

Table 2

Table 2. Summary statistics of aggregate reporter protein produced across all simulated cells after 20 s of simulation time.

Figure 4

Figure 4. Scatterplot of average reporter protein count per cell at end of simulation time as a function of number of cells simulated, with 1 SD error bars. Image produced by researcher.

To show that the pattern shown in Figure 3 was robust to starting TF concentration, we re-ran the simulation outlined in Figure 3 with different starting extracellular concentrations of TF. We found the decreasing pattern held.

The data show a positive, nonlinear relationship between the number of cells simulated and the total number of reporter protein molecules produced across all cells. The positive trend appears to be a consequence of the fact that there are more cells present that can produce the reporter protein, and the nonlinearity could be attributed to the dilution effect that decreases individual cellular response shown in Figure 3. Further analysis can be found in Supplementary Appendix E.

50-Second simulations

After 50 s of simulation time, the amount of reporter protein molecules present in each cell was averaged across trials and across all cells simulated. Those values are represented in this table. The full plots of the number of reporter protein molecules as a function of time per cell can be seen in Supplementary Figures B1 to B4. A 2-tailed t-test was performed on the results of simulating more than 1 cell and the results of a unicellular simulation to establish a statistically significant difference in the results of multicellular simulations compared to unicellular simulations. All t-tests produced p-values below the required threshold for significance, p < 0.01, which means the null hypothesis that multicellular simulations produce insignificantly different results than single cellular simulations can be rejected for all multicellular simulations tested. This strongly indicates that multicellular simulations produce statistically significantly different results than single cellular simulations at the single cell level, at least at 50 s of simulation time. However, the direction of difference between multicellular and single cellular simulations appears to be reversed for 50 s simulation results compared to 20 s simulation results, with the amount of reporter protein present in each cell increasing as the number of cells simulated increases. A scatter plot of the data in Table 3 is represented in Figure 5.

Table 3

Table 3. Summary statistics of average reporter protein count per cell after 50 s of simulation time.

Figure 5

Figure 5. Reproducing Figure 3 with logarithmic scan of the number of starting extracellular TF.

The data show a positive, nonlinear relationship between the number of cells and the number of reporter proteins present in each cell. The trend is initially increasing likely due to the increased cumulative production of the intercellular signal over time by more cells, counteracting the short run effects of dilution. The plateau could be caused by reaching the maximal production rate of the reporter protein, as production rate is limited by the amount of DNA in the cell that can transcribe mRNA to express both the intercellular signal and the reporter protein. The amount of intercellular signal produced by more cells may have completely saturated the available DNA in all cells.

After 50 s of simulation time, the total amount of reporter protein molecules present across all cells was averaged across trials. Those values are represented in this table. A 2-tailed t-test was performed on the results of simulating more than 1 cell and the results of a unicellular simulation to establish a statistically significant difference in the results of multicellular simulations compared to unicellular simulations. All t-tests produced p-values below the required threshold for significance, p < 0.001, which means the null hypothesis that multicellular simulations produce insignificantly different results than single cellular simulations can be rejected for all multicellular simulations tested. This strongly indicates that multicellular simulations produce statistically significantly different results than single cellular simulations at the population level, at least at 50 s of simulation time. A scatter plot of the data in Table 4 is represented in Figure 6.

Table 4

Table 4. Aggregate number of reporter protein molecules produced across all cells after 50 s of simulation time.

Figure 6

Figure 6. Scatterplot of aggregate number of reporter protein molecules present at 20 s into simulation time as a function of number of cells simulated. Image produced by researcher.

Similarly, in the 50 s simulations the pattern of increasing reporter proteins was shown to be robust.

The data show a positive, mostly linear relationship between the number of cells and the total amount of reporter protein produced across all cells. The positive relationship is likely from the increased capacity to produce reporter protein amongst more cells, while the linearity is likely the result of individual cells reaching maximal protein production rate such that the only factor influencing aggregate protein production is the total amount of DNA in the entire population, which scales linearly with the number of cells simulated. Further analysis can be found in Supplementary Appendix E.

Proof-of-concept bioreporter simulation

The results of a proof of concept bioreporter designed to leverage intercellular communication is summarized in Table 5 and Figures 7, 9. The design of the bioreporter is detailed in Figure 3. As shown, the amount of reporter protein increases as the number of detector bioreporter cells increases despite the starting analyte count remaining constant across experiments, demonstrating how a quorum of cells increases the signal strength of each individual cell and of the population as well, even in a more complex system with two shared chemical species across cells.

Table 5

Table 5. Results of multicellular bioreporter proof of concept. Simulations were run for 50 s with a starting analyte count of 10000 analyte molecules located in the extracellular environment for three replicates each.

Figure 7

Figure 7. Scatterplot of the average amount of reporter proteins present in each cell at the end of simulation time. Image produced by researcher.

Figure 8

Figure 8. Scatterplot of the aggregate amount of reporter proteins present in each cell at the end of simulation time. Image produced by researcher.

Figure 9

Figure 9. Reproducing Figure 5 with logarithmic scan of starting extracellular TF.

Runtime analysis

The time required by each simulation has been compiled and presented in this section.

The runtimes for the simulations in the work scale approximately linearly with the number of cells simulated.

Discussion

The purpose of this project was to investigate the effect of intercellular interactions on the single-cellular and population level responses in kinetic simulations. Multicellular simulations, typically in the context of quorum sensing or morphogenesis, have been explored by various researchers. However, no one has explicitly compared the difference that considering intercellular interactions has on cell- and population-level outcomes compared to only simulating a single cell. There are numerous popular applications for stochastic simulation of unicellular chemical reactions--such as Virtual Cell (Resasco et al., 2011) and COPASI (Hoops et al., 2006)--and it is important to understand how much unicellular simulation results will differ from multicellular simulations. Importantly, cellular simulations can be used for the computational validation of bioreporters, which are genetically engineered cells that produce a detectable signal in response to an environmental contaminant (Zeng et al., 2021). These bioreporters offer a low-cost method for onsite soil and water testing, however most research has focused on the single cell construction of bioreporters. Understanding the role of intercellular interactions on modulating cellular behavior as well as being able to simulate those effects opens a new dimension for the development of multicellular bioreporters to potentially achieve more sensitive detectors, enabling earlier detection of hazardous contaminants.

The hypothesis was that intercellular interactions would produce significantly different results than if they were not accounted for, both on the single cell level and on the population level. To test this hypothesis, a stochastic simulator was implemented to use the Gillespie algorithm to simulate nanoscale chemical kinetics with quantized probabilistic reactions, which more accurately describe chemical kinetics on the scale of single cells than ordinary differential equations (Gillespie, 1977). The simulator was used to simulate a hypothetical reaction network where the expression of reporter protein “P” is increased by the concentration of a transcription factor “TF”, which is also positively regulated by itself. TF is produced within a cell and diffuses freely through the cell membrane, acting as the main mechanism of intercellular interaction. P is created in the cell but cannot diffuse out, making it a useful indicator of transcriptional modulation. Levels of the reporter protein were used as indicators of the effect of intercellular interactions. Using this network, populations of 1, 10, 25, 50, 75, and 100 cells were simulated for 20 and 50 s and the amount of reporter protein present in each cell and in the entire population was recorded at the end of simulation time. The results of those simulations are discussed below.

In Figure 10, the amount of reporter protein produced in each cell decreases as more cells are simulated, with nearly all multicellular simulations producing statistically significantly different results from the unicellular simulation. This validates the hypothesis that multicellular simulations would produce statistically significantly different results from unicellular simulations at the single cellular level. The reason why reporter protein P production per cell would decrease could be explained by the faster rate of extracellular transcription factor TF depletion. The rate at which TF diffuses into a cell and thereby becoming unavailable to other cells increases as the number of cells increases, as they all share the same extracellular space. Therefore, there would be less intracellular concentrations of TF on average and less upregulation of P expression per cell due to lower TF concentrations. However, it may be possible that in the long run, more TF would be collectively produced by the increased number of cells such that it counteracts the increased rate of extracellular TF depletion.

Figure 10

Figure 10. Reporter protein production per cell at the end of 50 s of simulation time graphed against number of detector cells. Image produced by researcher.

Figure 2 indicates that the total amount of reporter proteins present in all cells at the end of 20 s of simulation time increases as the number of cells being simulated increases, with each multicellular simulation yielding statistically significant results compared to the unicellular simulation. This validates the hypothesis that multicellular simulations would yield significantly different results than unicellular simulations at the population. This trend is a natural consequence of accounting for more reporter protein to the total amount counted by including more cells in the total amount counted, and the trend is similar to how real bacteria must coordinate with each other to produce a visible glow in certain bioluminescent organisms, as single bacteria produce much less light than a colony (Smith and Grima, 2018). This graph does demonstrate that the dilution effect amongst multiple cells does not outweigh the increased reporter protein production caused by the presence of more cells and therefore greater capacity for reporter protein production. The upward trend suggests multicellular bioreporters are a promising candidate for better bioreporter signal strength. However, the graph does not increase linearly with the number of cells added, as might be naively expected of increasing the capacity for reporter production. The nonlinearity shown in the graph demonstrates the importance of actually simulating multiple cells when considering multicellular systems instead of simulating a single cell and multiplying its results by the number of cells in the system that is supposed to be simulated. If that were done here, the graph of those results would be perfectly linear, and an overestimate of the real results obtained from realistically simulating multiple cells. It is worth noting, however, that simulation of multiple cells is much more computationally intensive than approximations like scaling up the results of a single celled simulation. Future work would include optimizations to the algorithm used in this project for less prohibitive runtimes and better scalability on parallel architectures. Future work would also include comparison of different methods of approximating multicellular simulations with less complexity. For example, an analysis of how the results and runtime of a scaled single cellular simulation differs from a true multicellular simulation under different conditions could be performed. Similar analysis could also be performed on a simulation approximating a multicellular system as one large cell with modified amounts of DNA and membrane diffusion rates. Future work could even include tuning these approximations to find the best balance between accuracy of results and brevity of runtime.

The amounts of reporter protein present in each cell and across all cells were also recorded and analyzed after 50 s of simulation time. In Figure 3, the amount of reporter protein in each cell after 50 s of simulation time is plotted against the number of cells simulated. 2-Tailed t-tests show that every multicellular simulation produced results that were statistically significantly different from the results produced by the single cellular simulation, again validating the hypothesis. Interestingly, the amount of reporter protein present in each cell actually increases as the number of cells simulated increases, a reversal of the trend after 20 s of simulation time. This could be explained by the long term manifestation of the increased capacity to produce the intercellular signal due to the presence of more cells. In fact, the maximal rate of reporter protein production and intercellular signal production appears to have been reached, as indicated by the leveling off of the amount of reporter protein present in each cell for multicellular simulations of more than 25 cells. Greater long term production of the signal TF leads to greater concentrations in the extracellular space, which in turn leads to greater amounts of TF in each cell that elicit a greater response, as is shown in the graph. The effect of multiple cells positively influencing the behavior of single cells in the same population is reminiscent of the results found by Weber and Buceta in their work in simulating a similar positive feedback loop in a colony of cells in which a membrane-crossing molecule produced by each cell induces its own production (2013). Although they did not vary the amount of cells they simulated, they observed that their population of simulated cells, under the influence of an intercellular signal that they each produced, all cells increased their production of that signal up to a maximal production rate. Like in their work, a maximal rate of production was also reached in this project.

Figure 5 represents the total amount of reporter protein present across all cells after 50 s of simulation time. 2-Tailed t tests again showed every multicellular simulation produced statistically significant differences from the single cellular simulation, validating the hypothesis at the population level. The graph increases monotonically as the number of cells simulated increases. When more than 25 cells are simulated, the graph appears to increase a closely linear trend, likely a result of reaching the maximal protein production rate. As a result, the actual trend shape of how much protein is produced across all cells as a function of the number of cells in a population cannot be determined from these results. Future work that aims to further analyze the dynamics of aggregate protein production in a multicellular population need to consider methods to prevent complete saturation of protein production rate. The overall positive trend of the graph is not only similar to the previously discussed results by Weber and Buceta in 2013, but also the real life phenomenon of bacterial colonies using quorum sensing signals to coordinate luminescence (Li and Nair, 2012). Bacteria in colonies coordinate amongst the individual cells to fluoresce more brightly to take advantage of greater colony sizes, producing pulses of visible color.

Finally, the proof-of-concept multicellular bioreporter simulation did exhibit significant signal amplification on the population level and on the single cellular level with two shared chemical species: an analyte, and a cell produced signal protein that also acted as a communication medium. These results validate the secondary research hypothesis that a bioreporter using intercellular interactions would create a stronger signal in response to a given analyte. The proof of concept’s success in amplifying the response of each cell shows great promise for future research in multicellular bioreporter systems.

Figure 11 through 13 show how long each simulation took to complete. As diverse results were obtained with the current experimental parameters, no simulations with more cells or for longer periods of time were executed. However, future work could focus on expanding simulation times and/or cell population sizes. These larger systems would likely necessitate multiple days’ worth of compute-time without significant improvements to algorithmic efficiency or hardware available.

Figure 11

Figure 11. Total Reporter protein production across all bioreporter cells at the end of 50 s of simulation time graphed against number of detector cells. Image produced by researcher.

Figure 12

Figure 12. Runtimes of the 20 s simulations.

Figure 13

Figure 13. Runtimes of the 50 s simulations.

Figure 14

Figure 14. Runtimes of the 50 s bioreporter simulations.

Intriguingly, all simulations resulted in a decreasing standard deviation in average reporter protein concentration per cell as the number of cells simulated increased, even despite increasing numbers of reporter proteins in each cell. This agrees with Smith and Grima’s simulation of a two-stage gene simulation network in independent cells compared to tissue-bound cells (Smith and Grima, 2018). However, the trend appears to contradict Smith and Grima’s simulation of a system in which expression products dimerize. This indicates the networks constructed in this work are more like the former system than the latter. Using the authors’ intuitive explanation, it can be inferred that the constructed networks exhibit substantially more intracellular variability than variability induced through intercellular fluctuation. Similarly, our results also agree with Zhou et al., who found extracellular noise in a quorum sensing circuit served to synchronize cell responses (Zhou et al., 2005). This intracellular response variability in environments of varying cell populations can be qualitatively observed in Supplementary Appendix B.

Overall, simulation data established a significant difference in the results of multicellular and unicellular stochastic simulations in systems where the number of shared species--such as TF in this experiment--is affected by the cell population, thereby strongly supporting the hypothesis established at the beginning of experimentation at both the single cellular and population level. Future work would include improving multicellular simulation software to determine the long-term effects of multicellular simulation compared to unicellular simulation. Future work would mainly be focused on the simulation, design, and computational validation of multicellular bioreporter designs. The preliminary results obtained from this project suggest that leveraging intercellular interactions in multicellular designs is a very promising direction for developing biosensors to emit stronger and potentially more sensitive signals in response to low concentrations of contaminants. The ability to model these interactions allow for the computational validation of novel multicellular bioreporter designs before wet lab work, saving lab expenses and time spent on faulty designs. Thus, future work involves implementing real protein cascade networks with experimentally or computationally determined rate constants for all associated reactions before simulating those networks to design novel multicellular bioreporter designs. Furthermore, given the demonstrated nuances of multicellular simulations, spatial stochastic simulations would likely uncover even more complex patterns in multicellular simulations, where the transient locations of chemical species and cells can influence the progression of a system of chemical reactions (Marquez-Lago and Burrage, 2007). The Smoldyn simulator (Andrews, 2016) is an existing program used by Virtual Cell to implement spatial stochastic modeling, and could be adapted for multicellular simulations much like how Virtual Cell capabilities were adapted for the stochastic simulations described in this work. Another direction of future work would include methods to tackle the “inverse problem” as described by Warne et al. where reaction parameters are obtained through experimental data (Warne et al., 2019). As demonstrated, multicellular systems produce distinct results on the cellular level, so solving the inverse problem will be more complex in systems where multiple cells are interacting. A summary table of challenges of this approach and avenues of future work is listed in Table 6.

Table 6

Table 6. Summary of limitations and avenues of future improvement.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

DL: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing. MM: Conceptualization, Project administration, Supervision, Writing – review and editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2025.1595363/full#supplementary-material

References

Andrews, S. S. (2016). Smoldyn: particle-based simulation with rule-based modeling, improved molecular interaction and a library interface. Bioinformatics 33 (5), 710–717. doi:10.1093/bioinformatics/btw700

PubMed Abstract | CrossRef Full Text | Google Scholar

Antunes, L. C., and Ferreira, R. B. (2009). Intercellular communication in bacteria. Crit. Rev. Microbiol. 35 (2), 69–80. doi:10.1080/10408410902733946

PubMed Abstract | CrossRef Full Text | Google Scholar

Boettiger, A. N. (2013). Analytic approaches to stochastic gene expression in multicellular systems. Biophysical J. 105 (12), 2629–2640. doi:10.1016/j.bpj.2013.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Costello, A., and Badran, A. H. (2021). Synthetic biological circuits within an orthogonal central dogma. Trends Biotechnol. 39 (1), 59–71. doi:10.1016/j.tibtech.2020.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowan, A. E., Moraru, I. I., Schaff, J. C., Slepchenko, B. M., and Loew, L. M. (2012). Spatial modeling of cell signaling networks. Methods Cell Biol. 110, 195–221. doi:10.1016/B978-0-12-388403-9.00008-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Duso, L., and Zechner, C. (2020). Stochastic reaction networks in dynamic compartment populations. Proc. Natl. Acad. Sci. 117 (37), 22674–22683. doi:10.1073/pnas.2003734117

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81 (25), 2340–2361. doi:10.1021/j100540a008

CrossRef Full Text | Google Scholar

Gross, M. (2017). Roots of Mediterranean civilisations. Curr. Biol. 27 (24), R979–R982. doi:10.1016/j.cub.2017.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hahl, S. K., and Kremling, A. (2016). A comparison of deterministic and stochastic modeling approaches for biochemical reaction systems: on fixed points, means, and modes. Front. Genet. 7, 157. doi:10.3389/fgene.2016.00157

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., et al. (2006). COPASI—a complex pathway simulator. Bioinformatics 22 (24), 3067–3074. doi:10.1093/bioinformatics/btl485

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, M. E., Chen, A., Faeder, J. R., Henning, P., Moraru, I. I., Meier-Schellersheim, M., et al. (2021). Quantifying the roles of space and stochasticity in computer simulations for Cell Biology and cellular biochemistry. Mol. Biol. Cell 32 (2), 186–210. doi:10.1091/mbc.e20-08-0530

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., and Nair, S. K. (2012). Quorum sensing: how bacteria can coordinate activity and synchronize their response to external signals? Protein Sci. 21 (10), 1403–1417. doi:10.1002/pro.2132

PubMed Abstract | CrossRef Full Text | Google Scholar

Marquez-Lago, T. T., and Burrage, K. (2007). Binomial Tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics. J. Chem. Phys. 127 (10), 104101. doi:10.1063/1.2771548

PubMed Abstract | CrossRef Full Text | Google Scholar

Resasco, D. C., Gao, F., Morgan, F., Novak, I. L., Schaff, J. C., and Slepchenko, B. M. (2011). Virtual cell: computational tools for modeling in cell Biology. WIREs Syst. Biol. Med. 4 (2), 129–140. doi:10.1002/wsbm.165

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltepe, B., Kehribar, E. Ş., Su Yirmibeşoğlu, S. S., and Şafak Şeker, U. Ö. (2017). Cellular biosensors with engineered genetic circuits. ACS Sensors 3 (1), 13–26. doi:10.1021/acssensors.7b00728

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaff, J., Fink, C. C., Slepchenko, B., Carson, J. H., and Loew, L. M. (1997). A general computational framework for modeling cellular structure and function. Biophysical J. 73, 1135–1146. doi:10.1016/S0006-3495(97)78146-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S., and Grima, R. (2018). Single-cell variability in multicellular organisms. Nat. Commun. 9 (1), 345. doi:10.1038/s41467-017-02710-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Warne, D. J., Baker, R. E., and Simpson, M. J. (2019). Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art. J. R. Soc. Interface 16 (151), 20180943. doi:10.1098/rsif.2018.0943

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, M., and Buceta, J. (2013). Dynamics of the quorum sensing switch: stochastic and non-stationary effects. BMC Syst. Biol. 7 (1), 6. doi:10.1186/1752-0509-7-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, N., Wu, Y., Chen, W., Huang, Q., and Cai, P. (2021). Whole-cell microbial bioreporter for soil contaminants detection. Front. Bioeng. Biotechnol. 9, 622994. doi:10.3389/fbioe.2021.622994

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, T., Chen, L., and Aihara, K. (2005). Molecular communication through stochastic synchronization induced by extracellular fluctuations. Phys. Rev. Lett. 95 (17), 178103. doi:10.1103/physrevlett.95.178103

PubMed Abstract | CrossRef Full Text | Google Scholar