Inferring Gene Regulatory Neural Networks for Bacterial Decision Making
  in Biofilms
Bacterial cells are sensitive to a range of external signals used to learn
the environment. These incoming external signals are then processed using a
Gene Regulatory Network (GRN), exhibiting similarities to modern computing
algorithms. An in-depth analysis of gene expression dynamics suggests an
inherited Gene Regulatory Neural Network (GRNN) behavior within the GRN that
enables the cellular decision-making based on received signals from the
environment and neighbor cells. In this study, we extract a sub-network of
\textit{Pseudomonas aeruginosa} GRN that is associated with one virulence
factor: pyocyanin production as a use case to investigate the GRNN behaviors.
Further, using Graph Neural Network (GNN) architecture, we model a single
species biofilm to reveal the role of GRNN dynamics on ecosystem-wide
decision-making. Varying environmental conditions, we prove that the extracted
GRNN computes input signals similar to natural decision-making process of the
cell. Identifying of neural network behaviors in GRNs may lead to more accurate
bacterial cell activity predictive models for many applications, including
human health-related problems and agricultural applications. Further, this
model can produce data on causal relationships throughout the network, enabling
the possibility of designing tailor-made infection-controlling mechanisms. More
interestingly, these GRNNs can perform computational tasks for bio-hybrid
computing systems.
Samitha Somathilaka, Daniel P. Martins, Xu Li, Yusong Li, Sasitharan Balasubramaniam
# Inferring Gene Regulatory Neural Networks for Bacterial Decision Making in Biofilms

###### Abstract

Bacterial cells are sensitive to a range of external signals used to learn the environment. These incoming external signals are then processed using a Gene Regulatory Network (GRN), exhibiting similarities to modern computing algorithms. An in-depth analysis of gene expression dynamics suggests an inherited Gene Regulatory Neural Network (GRNN) behavior within the GRN that enables the cellular decision-making based on received signals from the environment and neighbor cells. In this study, we extract a sub-network of _Pseudomonas aeruginosa_ GRN that is associated with one virulence factor: pyocyanin production as a use case to investigate the GRNN behaviors. Further, using Graph Neural Network (GNN) architecture, we model a single species biofilm to reveal the role of GRNN dynamics on ecosystem-wide decision-making. Varying environmental conditions, we prove that the extracted GRNN computes input signals similar to natural decision-making process of the cell. Identifying of neural network behaviors in GRNs may lead to more accurate bacterial cell activity predictive models for many applications, including human health-related problems and agricultural applications. Further, this model can produce data on causal relationships throughout the network, enabling the possibility of designing tailor-made infection-controlling mechanisms. More interestingly, these GRNNs can perform computational tasks for bio-hybrid computing systems.

 Gene Regulatory Networks, Graph Neural Network, Biofilm, Neural Network.

## I Introduction

Bacteria are well-known for their capability to sense external stimuli, for complex information computations and for a wide range of responses [1]. The microbes can sense numerous external signals, including a plethora of molecules, temperatures, pH levels, and the presence of other microorganisms [2]. The sensed signals then go through the Gene Regulatory Network (GRN), where a large number of parallel and sequential molecular signals are collectively processed. The GRN is identified as the main computational component of the cell [3], which contains about 100 to more than 11000 genes (the largest genome identified so far belongs to _Sorangium cellulosum_ strain So0157-2) [4]. Despite the absence of neural components, the computational process through GRN allows the bacteria to actuate through various mechanisms, such as molecular production, motility, physiological state changes and even sophisticated social behaviors. Understanding the natural computing mechanism of cells can lead to progression of key areas of machine learning in bioinformatics, including prediction of biological processes, prevention of diseases and personalized treatment [5].

Bacterial cells are equipped with various regulatory systems, such as single/two/multi-component systems including _Quorum sensing_ (QS), to respond to environmental stimuli. The receptors and transporters on cell membranes can react and transport extracellular molecules, which subsequently interact with respective genes. In turn, the GRN is triggered to go through a complex non-linear computational process in response to the input signals. In the literature, it has been suggested that the computational process through the GRN of a bacterial cell comprises a hidden neural network (NN)-like architecture [6, 7]. This indicates that, even though bacterial cells can be categorized as non-neural organisms, they perform neural decision-making processes through the GRN. This results in recent attention towards Molecular Machine Learning systems, where AI and ML are developed using molecular systems [8]. In these systems, several neural components can be identified in GRNs, in which genes may be regarded as

Fig. 1: Illustration of the Gene Regulatory Neural Networks (GRNN) extraction and the implementation of the GNN to model the biofilm. The diffusion of molecules from one cell to another is modeled as a vector, where \(m_{q}\) represents the concentration of the \(q^{th}\) molecular signal.
computational units or neurons, transcription regulatory factors as weights/biases and proteins/second messenger Molecular Communications (MC) as neuron-to-neuron interactions. Owing to a large number of genes and the interactions in a GRN, it is possible to infer sub-networks with NN behaviors that we term Gene Regulatory Neural Networks (GRNN). The non-linear computing of genes results from various factors that expand through multi-omics layers, including proteomic, transcriptomic and metabolomic data (further explained in Section II-A). In contrast, the GRNN is a pure NN of genes with summarized non-linearity stemmed from multi-omics layers with weights/biases.

Identification of GRNNs can be used to model the decision-making process of the cell precisely, especially considering simultaneous multiple MC inputs or outputs. However, due to the limited understanding and data availability, it is still impossible to model the complete GRN with its NN-like behaviors. Therefore, this study uses a GRNN of _Pseudomonas aeruginosa_ that is associated with _PhoR-PhoB_ and _BqsS-BqsR_ two-component systems (TCSs) and three QS systems related to pyocyanin production as a use case to explore the NN-like behaviors. Although a single bacterium can do a massive amount of computing, they prefer living in biofilms. Hence, in order to understand the biofilm decision-making mechanism, we extend this single-cell computational model to an ecosystem level by designing an _in-silico_ single species biofilm with inter-cellular MC signaling as shown in Fig. 1.

The contributions of this study are as follows:

* **Extracting a GRNN:** Due to the complexity and insufficient understanding of the gene expression dynamics of the full GRN, we only focus on a sub-network associated with pyocyanin production (shown in Fig. (a)a) to investigate the NN-like computational behavior of the GRN. Further, the genes of extracted sub-network are arranged following a NN structure that comprises input, hidden and output layers, as shown in Fig. (b)b.
* **Modeling a biofilm as a GNN:** The GRNN only represents the single-cell activities. To model the biofilm-wide decision-making process, we use a Graph Neural Network (GNN). First, we create a graph network of the bacterial cell and convert it to a GNN by embedding each node with the extracted GRNN as the update function. Second, the diffusion-based MCs between bacterial cells in the biofilm are encoded as the message-passing protocol of the GNN, as shown in Fig. (c)c.
* **Exploring the behaviors of the GRNN and intracellular MC dynamics to predict cell decisions:** The output of the GRNN is evaluated by comparing it with the transcriptomic and pyocyanin production data from the literature. Finally, an edge-level analysis of the GRNN is conducted to explore the causal relationships between gene expression and pyocyanin production.

This paper is organized as follows: Section II explains the background of bacterial decision-making in two levels: cellular-level in Section II-A and population-level in Section II-B, while the background on the _P. aeruginosa_ is introduced in Section II-C. Section III is dedicated to explaining the model design of cellular and population levels. The results related to model validation and the intergenic intra-cellular signaling pattern analysis are presented in Section IV and the study is concluded in Section V.

## II Background

As the model expands through single cellular and biofilm-wide decision-making layers, this section provides the background of how a bacterium uses the GRN to make decisions and how bacterial cells make decisions in biofilms. Moreover, we briefly discuss the cellular activities of the _Pseudomonas aeruginosa_ as it is the use case species of this study.

### _Decision-Making Process of an Individual Cell_

Prokaryotic cells are capable of sensing the environment through multiple mechanisms, including TCSs that have been widely studied and it is one of the focal points of this study. The concentrations of molecular-input signals from the extracellular environment influence the bacterial activities at the cellular and ecosystem levels [9]. Apart from the extracellular signals of nutrients, it is evident that the QS input signals have a diverse set of regulative mechanisms in biofilm-wide characteristics, including size and shape [10]. These input signals undergo a computational process through the GRN, exhibiting a complex decision-making mechanism. Past studies have explored and suggested this underpinning computational mechanism in a plethora of directions, such

Fig. 2: Extraction of a GRNN considering a specific sub-network of the GRN where a) is the two-component systems (TCSs) and QS network that is associated with the pyocyanin production, b) is the derived GRNN that is equipped with hypothetical nodes (**hns**) without affecting its computation process to form a symmetric network structure and c) is the conversion of real biofilm to the suggested _in-silico_ model.
as using differential equations [11] and probabilistic Boolean networks [12] and logic circuit [13]. All of these models mainly infer that the bacterial cells can make decisions not just based on the single input-output combinations, but they can integrate several incoming signals non-linearly to produce outputs.

The studies that focus on differences in gene expression levels suggest that a hidden weight behavior controls the impact of one gene on another [6]. This weight behavior emerges through several elements, such as the number of transcription factors that induce the expression, the affinity of the transcription factor binding site, and machinery such as thermoregulatory and enhancers/silencers [14, 15]. Fig. 3 depicts a set of factors influencing the weight between genes. The weight of an edge between two genes has a higher dynamicity as it is combinedly determined by several of these factors. Based on environmental conditions, the GRN of the bacterial cell adapts various weights to increase the survivability and repress unnecessary cellular functions to preserve energy. An example of such regulation is shown in Fig. 4 where a _P. aeruginosa_ cell uses a thermoregulatory to regulate the QS behaviors. Fig. 3(a) has a set of relative weights based on cell activities in an environment at 37\(\,{}^{\circ}\)C, while Fig. 3(b) represents weights at 30\(\,{}^{\circ}\)C. The weights between the **hn21** and _rhlR_ are different in two conditions, and these cellar activities are further explained in [14].

### _Biofilm Decision-Making_

Even though an individual cell is capable of sensing, computing, and actuating, the majority of bacterial cells live in biofilms, where the survivability is significantly increased compared to their planktonic state. Biofilm formation can cause biofouling and corrosion in water supply and industrial systems [16]. However, biofilms formation can be desirable in many situations, for example, bioreactors in wastewater treatment [17], bioremediation of contaminated groundwater [18, 19], where biofilms serve as platforms for biogeochemical reactions. A massive number of factors can influence biofilm formation, including substrant surface geometrical characteristics, diversity of species constituting the biofilm, hydrodynamic conditions, nutrient availability, and especially communication patterns [20] where the TCS and QS play significant roles. A TCS comprises a _histidine kinase_ that is the sensor for specific stimulus and a cognate response regulator that initiates expressions of a set of genes [21]. Hence, in each stage, essential functions can be traced back to their gene expression upon a response to the input signals detected by bacterial cells. For instance, in the first stage of biofilm formation, the attachment of bacteria to a surface is associated with sensing a suitable surface and altering the activities of the flagella. In the next stage, _rhramnolipids_ production is associated with ferric iron Fe\({}^{3+}\) availability in the environment, mostly sensed through _BqsS-BqsR_ TCSs. Further, Fe\({}^{3+}\) was identified as a regulator of _pqSA_, _pqsR_, and _pqsE_ gene expressions that are associated with the production of two critical components for the formation of microcolonies: eDNA and EPS [22]. Similarly, in the final stage, the dispersion process can also be traced back to a specific set of gene regulations, including _bdlA_ an _rbdA_[23, 24]. An understanding of the underlying decision-making process of bacteria may enable us to control their cellular activities.

### _Pseudomonas Aeruginosa_

The main reason for selecting _P. aeruginosa_ in this work lies in its alarming role in human health. For example, this species is the main cause of death in cystic fibrosis patients [25]. _P. aeruginosa_ is a gram-negative opportunistic pathogen with a range of virulence factors, including pyocyanin and cytotoxin secretion [26]. These secreted molecules can lead to complications such as respiratory tract ciliary dysfunction and induce proinflammatory and oxidative effects damaging the host cells [27]. The biofilms are being formed on more than 90% endotracheal tubes implanted in patients who are getting assisted ventilation, causing upper respiratory tract infections [28]. In addition, another important reason for targeting _P. aeruginosa_ is the data availability for the GRN structure [29], pathways [30], genome [31], transcriptome [32] and data from mutagenesis studies [33, 34]. Compared to the complexity of the GRN, the amount of data and information available on the

Fig. 4: Two GRNN setups with different weights associated with two environmental conditions. a) is the relative weight setup of _P. aeruginosa_ cell in 37\(\,{}^{\circ}\)C and b) is in 30\(\,{}^{\circ}\)C.

Fig. 3: Illustration of gene expression regulators that are considered the weight influencers of the edges of GRNN. Here, the \(\alpha_{(\varphi)}\), \(\alpha_{(\sim\varphi)}\), \(\alpha_{(TF)}\), \(\alpha_{(Rep)}\), \(\alpha_{(eTF)}\) and \(\alpha_{(sTF)}\) are relative concentrations of sigma factors, anti-sigma factors, transcription factors (TFs), repressors, enhancer-binding TFs and silencer-binding TFs respectively. Moreover, \(\beta_{(Prom)}\), \(\beta_{(Eph)}\), and \(\beta_{(Sil)}\) are the binding affinities of the promoter, operator, enhancer and silencers regions respectively.
gene-to-gene interactions and expression patterns is insufficient to develop an accurate full _in-silico_ model. Therefore, we chose a set of specific genes that are associated with QS, TCS, and pyocyanin production.

## III System design

This section explains the system design in two main phases, extracting a NN-like architecture from the GRN targeting the set of genes and creating a model of the biofilm ecosystem.

### _Extracting Natural Neural Network from GRN_

We first fetch the structure of the GRN graph from Reg-ulomePA [29] database that contains only the existence of interactions and their types (positive or negative). As the next step, using information from the past studies [35, 36, 37, 38], we identified the genes involved in the **Las**, **Rhl** and **PQS** QS systems, _PhoR-PhoB_ and _BqsS-BqsR_ TCSs, and pyocyanin production to derive the sub-network of GRN as shown in Fig 2a. We further explored the expression dynamics using transcriptomic data [39, 40] where we observed the non-linearity in computations that are difficult to capture with existing approaches such as logic circuits, etc. [6], making the NN approach more suitable. However, a NN model with a black box that is trained on a large amount of transcriptomic data records to do computations similar to the GRN has a number of limitations, especially in understanding the core of the computational process [41]. Our model does not use a conventional NN model; instead, we extract a NN from the interaction patterns of the GRN, which we consider a pre-trained GRNN. In this sub-network, we observed that the lengths of expression pathways are not equal. For example, the path from _PhoR-PhoB_ to the _phz2_ gene has two hops, but the path from the _BqsS-BqsR_ system to the _rhlR_ gene only has one hop. The extracted network has the structure of a random NN. Hence, we transform this GRNN to Gene Regulatory Feedforward Neural Network by introducing hypothetical nodes (**hns**) that do not affect the behaviors of the GRNN as shown in Fig 2b. In this transformation, we decide the number of hidden layers based on the maximum number of hops in gene expression pathways. In our network, the maximum number of hops is two, which determines the number of hidden layers as one, and then the number of hops of all the pathways is leveled by introducing **hns**. If a **hn** is introduced between a source and target genes, the edge weights from the source node to the **hn** and from **hn** to the target node are made "1" so that the **hn** does not have an influence on the regulation of genes. Moreover, if a gene does not induce another in the network, the weight of the edge between that pair is made "0".

Here, we summarize multiple factors of interaction into a weight that determines the transcriptional regulation of a particular gene. This regulation process occurs when the gene products get bound to the promoter region of another, influencing the transcriptional machinery. Hence, we observe this regulation process of a target gene as a multi-layered model that relies on the products of a set of source genes, the interaction between gene products, and the diffusion dynamics within the cell. Creating a framework to infer an absolute weight value using all the above factors is a highly complex task. In order to infer weight, one method is to train a NN model with the same structure as the GRN using a series of transcriptomic data. However, this approach also has numerous challenges, such as the lack of a sufficient amount of data in similar environments.

Therefore, we estimate a set of relative weights based on genomic, transcriptomic, and proteomic explanations of each interaction from the literature. The weights were further fine-tuned using the transcriptomics data. A relative weight value of an edge can be considered a summarizing of multi-layer transcriptional-translation to represent the impact of the source gene on a target gene.

In this computational process, we identify another layer of interactions that occur within the cell. The produced molecules by the considered TCs network go through a set of metabolic interactions that are crucial for the functionality of the cell. Since our primary goal is to explore the NN behaviors of GRN, we model these inter-cellular chemical reactions as a separate process, keeping the gene-to-gene interactions and metabolic interactions in two different layers. To model the complete pyocyanin production functionality of the cell, we use the inter-cellular molecular interactions shown in Fig 5. Here, RhlR is a transcriptional regulator of _P. aeruginosa_ that forms a complex by getting attached to its cognate inducer C4-HSL and then binds to the promoter regions of relevant genes [42]. Similarly, LasR transcriptional regulator protein and 3-oxo-C12-HSL (3OC), and PqsR with PQS and HHQ form complexes and get involved in the regulation of a range of genes [43, 44]. Further, C10H10O6 in the environment are converted by the _P. aeruginosa_ cells in multiple steps using the products of the GRN we consider. First, C10H10O6 is converted into phenazine-1-carboxylic using the enzymes of _pqsABCDEFG_ genes. Later, phenazine-1-carboxylic was converted into 5-Methylphenazine-1-carboxylate, and finally, 5-Methylphenazine-1-carboxylate into Pyocyanin by PhzM and PhzS, respectively [45].

Molecular accumulation within a bacterial cell can be considered its memory module where certain intra-cellular interactions occurs. Therefore, we define an internal memory matrix \(IM\) as,

\[\begin{array}{ccccc}&im_{1}&im_{2}&...&im_{J}\\ &B_{1}&\begin{pmatrix}C_{1,im_{1}}^{(t)}&C_{1,im_{2}}^{(t)}&...&C_{1,im_{J}}^{( t)}\\ C_{(2,im_{1})}^{(t)}&C_{(2,im_{2})}^{(t)}&...&C_{2,im_{J}}^{(t)}\\ \vdots&\vdots&\ddots&\vdots\\ &B_{P}&\begin{pmatrix}C_{(P,im_{1})}^{(t)}&C_{(P,im_{2})}^{(t)}&...&C_{(P,im_{J })}^{(t)}\\ \end{pmatrix},\end{array}\end{array} \tag{1}\]

where the concentration of the internal molecule \(im_{j}\) is

Fig. 5: Illustrations of intra-cellular metabolite interaction.
\(C_{(i,im_{j})}^{(t)}\).

GRNN process molecular signals from the environment and other cells. Hence, we used the approach of GNN as a scalable mechanism to model the MCs and biofilm wide decision-making process. The extreme computational power demand of modeling the diffusion-based MCs of each cell is also avoided by using this approach.

### _Graph Neural Network Modeling of Biofilm_

First, the biofilm is created as a graph network of bacterial cells where each node is a representation of a cell, and an edge between two nodes is a MC channel. We convert the graph network into a Graph Neural Network (GNN) in three steps: 1) embedding the extracted GRNN of pyocyanin production into each node as the update function, 2) encoding the diffusion-based cell-to-cell MC channels as the message passing scheme, and 3) creating an aggregation function at the reception of molecular messages by a node as shown in Fig. 6. Next, we define feature vectors of each node of the GNN to represent the gene expression profile of the individual cell at a given time. Subsequently, considering \(L\) is the number of genes in the GRNN, \(P\) is the number of bacterial cells in the biofilm and \(b_{(i,g)}^{(t)}\) is the expression of gene \(g_{l}\) by the bacteria \(B_{i}\), we derive the following matrix \(\mathbf{FV}^{(t)}\) that represents all the feature vectors of the GNN at time \(t\).

\[\mathbf{FV}^{(t)}=\begin{array}{ccccc}g_{1}&g_{2}&...&g_{L}\\ B_{1}&\begin{pmatrix}b_{(1,g_{1})}^{(t)}&b_{(1,g_{2})}^{(t)}&...&b_{(1,g_{L}) }^{(t)}\\ b_{(2,g_{1})}^{(t)}&b_{(2,g_{2})}^{(t)}&...&b_{(2,g_{L})}^{(t)}\\ \vdots&\vdots&\ddots&\vdots\\ B_{P}&\begin{pmatrix}b_{(P,g_{1})}^{(t)}&b_{(P,g_{2})}^{(t)}&...&b_{(P,g_{L}) }^{(t)}\end{pmatrix}\end{array}\end{array} \tag{2}\]

The computational output of the GRNN of each node results in the secretion of a set of molecules that are considered messages in our GNN model as illustrated in the Fig. 7.

When the number of molecular species considered in the network is \(Q\) and output \(m_{q}\) molecular message from bacterial cell \(B_{i}\) at TS \(t\) is \(msg_{(i,m_{q})}^{(t)}\), we derive the matrix

\[\mathbf{MSG}^{(t)}=\begin{array}{ccccc}m_{1}&m_{2}&...&m_{Q}\\ B_{1}&\begin{pmatrix}msg_{(1,m_{1})}^{(t)}&msg_{(1,m_{2})}^{(t)}&...&msg_{(1, m_{Q})}^{(t)}\\ msg_{(2,m_{1})}^{(t)}&msg_{(2,m_{2})}^{(t)}&...&msg_{(1,m_{Q})}^{(t)}\\ \vdots&\vdots&\ddots&\vdots\\ B_{P}&\begin{pmatrix}msg_{(P,m_{1})}^{(t)}&msg_{(P,m_{2})}^{(t)}&...&msg_{(P, m_{Q})}^{(t)}\end{pmatrix}\end{array} \tag{3}\]

Further, we use a static diffusion coefficients vector

\[\mathbf{D}=\{D_{m_{1}},D_{m_{2}},...,D_{m_{Q}}\}, \tag{4}\]

where \(D_{m_{q}}\) is diffusion coefficient of molecular species \(m_{q}\).

We define another matrix \(\mathbf{ED}\) that contains the euclidean distances between bacterial cells in the biofilm as follows

\[\mathbf{ED}=\begin{array}{ccccc}B_{1}&B_{2}&...&B_{P}\\ B_{1}&\begin{pmatrix}d_{(1,1)}&d_{(1,2)}&...&d_{(1,P)}\\ B_{2}&d_{(2,1)}&d_{(2,2)}&...&d_{(2,P)}\\ \vdots&\vdots&\ddots&\vdots\\ B_{P}&\begin{pmatrix}d_{(P,1)}&d_{(P,2)}&...&d_{(P,P)}\end{pmatrix}\end{array} \tag{5}\]

where \(d_{i,j}\) is the euclidean distance between the \(i^{th}\) and \(j^{th}\) cells.

Fig. 6: Illustration of the GNN components where a) is a snapshot of the bacterial network that has the gene expression profile as the feature vector. Further, this gene expression pattern of a cell is encoded to a message of secreted molecules where MC plays a crucial role. Moreover, b) shows the temporal behavior of the GNN, that the output of one graph snapshot influences the next.

Fig. 7: The process of one GRNN outputs reaching another GRNN as molecular messages.
The feature vector of \(i^{th}\) bacterial cell at the TS \(t+1\) is then modeled as,

\[\mathbf{FV}_{i}^{(t+1)}=GRNN_{i}(\mathbf{MSG}_{i}^{(t)}+S_{i}^{(t)}) \tag{6}\]

where \(MSG_{i}^{(t)}\) is the message generated by the same cell in the previous TS. The \(GRNN_{i}\) is the extracted GRNN that is the update function in the GNN learning process and \(\mathbf{S}_{i}^{(t)}=\mathbf{R}_{i}^{(t)}+\mathbf{K}_{i}^{(t)}\), is the aggregate function. In the aggregation component, the \(\mathbf{R}_{i}^{(t)}\) is the incoming signals from peer bacterial cells and \(K_{(i:m_{q})}^{(t+1)}\) is the external molecule input vector at the location of \(B_{i}\) and the TS \(t\) that is expressed as

\[\mathbf{K}_{i}^{(t+1)}=\big{\{}K_{i:m_{1}}^{(t+1)},K_{:i:m_{2}}^{(t+1)},...,K_{ :m_{Q}}^{(t+1)}\big{\}}. \tag{7}\]

In order to compute \(\mathbf{R}_{i}^{(t+1)}\), we use a matrix \(\mathbf{Y}_{i}\); \(\mathbf{Y}_{i}=\overset{\leftrightarrow}{1}_{[Q\times 1]}\times\)ED\({}_{i}\), where \(\overset{\leftrightarrow}{1}_{[Q\times 1]}\) is an all-ones matrix of dimension \(Q\times 1\). The \(\hat{g}\) matrix is then defined as follows,

\[\mathbf{\hat{g}}(\mathbf{D}^{\intercal},\mathbf{Y},t)= \tag{8}\] \[\begin{bmatrix}g(D_{m_{1}},d_{(i,1)},t)&g(D_{m_{1}},d_{(i,2)},t)&...&g(D_{m_{1}},d_{(i,P)},t)\\ g(D_{m_{2}},d_{(i,1)},t)&g(D_{m_{2}},d_{(i,2)},t)&...&g(D_{m_{2}},d_{(i,P)},t) \\ \vdots&\vdots&\ddots&\vdots\\ g(D_{m_{Q}},d_{(i,1)},t)&g(D_{m_{Q}},d_{(i,2)},t)&...&g(D_{m_{Q}},d_{(i,P)},t) \end{bmatrix}\]

In the above matrix, \(g(D_{m_{l}},d_{(i,j)},t)\) is the Green's function of the diffusion equation as shown below,

\[\hat{G}(D_{m_{l}},d_{(i,j)},t)=\frac{1}{\left(4\pi D_{m_{l}}t\right)^{\frac{ 3}{2}}}\exp\left(-\frac{d_{(i,j)}^{2}}{4D_{m_{l}}t}\right). \tag{9}\]

Further, the incoming signal vector \(\mathbf{R}_{i}^{(t+1)}\) is denoted as below,

\[\mathbf{R}_{i}^{(t+1)}=diag\big{(}\mathbf{\hat{g}}(\mathbf{D}^{\intercal}, \mathbf{Y},t)\times\mathbf{MSG}^{(t)}\big{)}. \tag{10}\]

Further, we equip our model with a 3-D environment to compensate for the noise element and external molecule inputs. Environment-layer is designed as a 3-D grid of voxels that can store precise information on external nutrients (similarly to our previous model in [46]). The diffusion of nutrient molecules through the medium is modeled as a random-walk process. This layer allows us to enrich the model with the dynamics of nutrient accessibility of bacterial cells due to diffusion variations between the medium and the Extracellular Polymeric Substance (EPS).

The bacterial cells in the ecosystem also perform their own computing tasks individually, resulting in a massively parallel processing framework. Hence, we use the python-cuda platform to make our model closer to the parallel processing architecture of the biofilm, where we dedicate a GPU block for each bacterial cell and the threads of each block for the matrix multiplication of the GRNN computation associated with the particular cell. Additionally, due to the massive number of iterative components in the model, the computational power demand faces significant challenges with serial programming making parallelization the best match for the model.

## IV Simulations

In this section, we first explain the simulation setup and then discuss the results of gene expression and molecular production dynamics to prove the accuracy of the extracted GRNN, emphasizing that it works similarly to the real GRN. Later, we use computing through the GRNN to explain certain activities of the biofilm.

### _Simulation Setup_

As our interest is to investigate the NN-like computational process, we do not model the formation process of the biofilm, but we only remodel a completely formed biofilm and disregarding the maturation and dispersion stages. In this model, we consider the biofilm as a static 3-D structure of bacterial cells. Hence, we first place bacterial cells randomly in the model in a paraboloid shape using the equation, \(z<\frac{x^{2}}{5}+\frac{y^{2}}{5}+20\) where \(x\), \(y\) and \(z\) are the components of 3-D Cartesian coordinates. This paraboloid shape is chosen to make the spacial arrangement of the cells close to real biofilm while keeping the cell placement process mathematically simple. Within this 3-D biofilm region, we model the diffusivity according to \(D_{B}/D_{aq}=0.4\), which is the mean relative diffusion [47] where \(D_{B}\) and \(D_{aq}\) are the average molecular diffusion coefficients of the biofilm and pure water, respectively. Further, to start the simulation at a stage where the biofilm is fully formed and the MC is already taking place, we filled the internal memory vector of each cell with the average molecular level at the initial TS. Each bacterial cell will use the initial signals from the internal memory and use its GRNN to process and update the feature vector for the next TS. Table I presents the parameter descriptions and values used for the simulation. As shown in Table I, the model runs for 150 TSs, generating data on a range of functions for the system. For instance, this model can produce data on feature vector of each cell, MC between cells,

\begin{table}
\begin{tabular}{l r r} \hline
**Parameter** & **Value** & **Description** \\ \hline No. of cells & 2000 & The number of cells is limited due to the memory availability of the server. \\ No. of genes & 13 & The network only consists of the gene that are directly associated with QS, _PhoR-PhoB_ and _BqsS-BqsR_ TCSs, and pyeocyanin production. \\ No. internal memory molecules & 16 & The set of molecules that involved in QS, _PhoR-PhoB_ and _BqsS-BqsR_ TCSs,and ppyocyanin production. \\ No. messenger molecules & 4 & The number of molecules that were exchanged between cells in the sub network. \\ Dimensions of the environment & 20x20x20\(\mu m\) & The dimensions were fixed considering the average sizes of _P. aeruginosa_ biofilms and computational demand of the model. \\ Duration & 150 TSs & The number of TSs can be modified to explore the cellular and ecosystem level activities. For this experiment we fixed a TS to represent 30mins. \\ No. iterations per setup & 10 & Considering the stochasticity ranging from the gene expression to ecosystem-wide communications, the experiments were iterated 10 times. \\ \hline \end{tabular}
\end{table} TABLE I: Parameters utilised in the system development
molecular consumption by cells, secretion to the environment, and nutrient accessibility of cells for each TS.

In order to prove that our GRNN computes similarly to the natural bacterial cell and collective behaviors of the cells are the same as the natural biofilm, we conduct a series of experiments. We explore the GRNN computation and biofilm activities under High Phosphate (HP) and Low Phosphate (LP) levels using eight experimental setups as follows, 1) wild-type bacteria (WD) in LP, 2) _lasR_ mutant (_lasR\(\Delta\)_) in LP, 3) phoB mutant (_phoB\(\Delta\)_) in LP, 4) _lasR_ & PhoB double mutant (\(\text{LasR}\Delta\text{PhoB}\Delta\)) in LP, 5) WD in HP, 6) _lasR\(\Delta\)_ in HP, 7) _PhoB\(\Delta\)_ in HP and _LasR\(\Delta\)PhoB\(\Delta\)_ in HP. While the WD uses the full GRNN, _lasR\(\Delta\)_ is created by making the weight of the link between _hn22_ and _lasR_ as "0". Further, the GRNN of _phoB\(\Delta\)_ is created by making the weights of links from _PhoB_ to _hn23_ and _PhoB_ to _pqsABCDE_ also "0".

### _Model Validation_

First, we show the nutrient accessibility variation in the biofilm through Fig 8. The cells in the biofilm core have less accessibility while the cells closer to the periphery have more access to nutrients due to variations in diffusion between the environment and the EPS. Fig 8 shows that when a low phosphate concentration is introduced to the environment, the direct access to the nutrient by the cells is limited. After the \(TS=10\), around 60% of cells have access to 20% of the nutrient concentration. Further, Fig 8 shows that the increased nutrient introduction to the environment positively reflects on the accessibility. This accessibility plays a role mainly in the deviation of gene expression patterns resulting in phenotypic differentiations that is further analyzed in Section IV-C.

Comparing the predictions of molecular production through GRNN computing with the wet-lab experimental data from the literature, we are able to prove that the components of the GRN work similarly to a NN. Fig. 9 shows the pyocyanin accumulation variations of the environment in the eight setups mentioned earlier as results of decision-making of the GRNN. Production of pyocyanin of the WD _P. aeruginosa_ biofilms is high in LP, compared to the HP environments as shown in Fig: 9. Further, the same pattern can be observed in the _lasR\(\Delta\)_ biofilms, but with a significantly increased pyocyanin production in LP as shown in Fig. 9. The _phob\(\Delta\)_ and _LasR\(\Delta\)phob\(\Delta\)_ biofilms produce a reduced level of pyocyanin compared to WD and _LasR\(\Delta\)_ that are shown in Fig. 9 and Fig. 9 respectively. We then present a comparison between GRNN prediction and wet-lab experimental data [48] as ratios of HP to LP in Fig. 10. The differences between pyocyanin production through GRNN in HP and LP condition for all the four setups in Fig. 10 are fairly close to the wet-lab data. In the WD setup, the difference between the GRNN model and wet-lab data only has around 5% difference, while deviations around 10% can be observed in _lasR\(\Delta\)_ and _phoB\(\Delta\)_. The most significant deviation around 20% of pyocyanin production difference is visible in _lasR\(\Delta\)phoB\(\Delta\)_ that is caused by the lack of interaction from other gene expression pathways, as we only extracted a sub-network portion of the GRN. Therefore, these results prove that the extracted GRNN behaves similarly to the GRN dynamics.

We further prove that the GRNN computing process performs similarly to the GRN by comparing the gene expression behaviors of the model with the wet-lab data [48] as shown in Fig. 11. First, we show the expression dynamics of genes _lasI_, _pqsA_ and _rhR_ of WD in LP in Fig. 11, Fig. 11 and Fig. 11 respectively. All the figures depict that gene expression levels are higher in LP compared to HP until around \(TS=100\). Beyond that point, relative gene expression levels are close to zero as the the environment run out of nutrients. Moreover, the differences in gene expression levels predicted by the GRNN

Fig. 8: The nutrient accessibility variations of cells is expressed in two different environment conditions: a) low phosphate and b) high phosphate concentrations.

Fig. 10: Evaluation of the model accuracy by comparing HP to LP pyocyanin production ratio with wet-lab data from [48].

Fig. 9: Relative Pyocyanin accumulation of four different biofilms of a) WD, b) lasR\(\Delta\), c) phoB\(\Delta\) and d) lasR\(\Delta\)phob\(\Delta\) in both low and high phosphate levels.
computing for LP and HP are also compared with the wet-lab data in Fig. (d)d. In this comparison, it is evident that the predicted gene expression differences of all three genes are close to the wet-lab data with only around 10% variation. The performance similarities between the GRNN and real cell activities once again prove that the GRN has underpinning NN-like behaviors.

### _Analysis of GRNN Computing_

Fig. 12 and Fig. 13 are used to show the diverse information flow of the GRNN that cause the variations in pyocyanin

Fig. 11: Expression levels of three different genes to that were used to prove the accuracy of the GRNN: a) _lasI_, b) _pqSA_, c) _rhIR_ expression levels in LP and HP and d) comparison between GRNN computing results and wet-lab data.

Fig. 12: Gene expression and associated information flow variations in GRNNs of a) WD b) _lasR\(\Delta\)_ and c) phoB\(\Delta\) in LP.

Fig. 13: Gene expression and associated information flow variations in GRNNs of a) WD b) _lasR\(\Delta\)_ and c) phoB\(\Delta\) in HP.
production in LP and HP conditions, respectively. Here, we use gene expression profiles extracted from one bacterial cell located at (7, 9, 2)\(\mu\)m in the Cartesian coordinates that is in the middle region of the biofilm with limited access to the nutrients. First, the gene expression variations of WD, _lasR\(\Delta\)_, and _phob\(\Delta\)_ bacterial cells in LP (Fig. 12) and HP (Fig. 13) are shown for \(TS<50\). Next, the information flow through the GRNN is illustrated above each expression profile at time \(TS=20\), where the variations will be discussed. In Fig. (a)a, impact of the inputs 3OC-LasR and phosphate cause higher expression levels of the nodes **hn12** and _phoB_ in the input layer that cascade the nodes _phZ1_, _phZ2_, _pqsR_, _lasR_, 3OC, _rhlR_ and _PqsH_ in the output layer at \(TS=20\). Fig. (b)b has significantly higher _pqsA_ operon expression levels compared to HP conditions (Fig. (b)b), reflecting higher yyocanin production that can be seen in Fig. (b)b. Nevertheless, the reduced gene expression levels, except _pqsA_ operon, of _lasR\(\Delta\)_ biofilm in both LP (Fig. (b)b) and HP (Fig. (b)b) conditions compared to the other setups emphasize that the inputs via inter-cellular MC significantly alter GRNN computing outputs. In contrast, only a smaller gene expression difference can be observed between the two setups of _phob\(\Delta\)_ in LP (Fig. (c)c and _phob\(\Delta\)_ in HP (Fig. (c)c) resulting in minimized myocardium production differences as shown earlier in Fig. (c)c.

The GRNN model supports the understanding of the gene expression variations due to factors such as nutrient accessibility, where in our case is a single species biofilm. Fig. 14 depicts the variability in the gene expression levels for four different locations of the biofilm at \(TS=3\). Fig. (a)a and Fig. (b)b are the gene expression profiles and the signal flow through GRNN pairs of two cells located close to the attached surface and the center of the biofilm. The phosphate accessibility for these two locations is limited. Hence, edges from _phob_ have a higher information flow compared to the other two cells near the periphery of the biofilm, which can be observed in Fig. (c)c and Fig. (d)d. The microbes in the center (Fig. (a)a) and the bottom (Fig. (b)b) mainly have access to the inter-cellular MCs, while the other two bacteria have direct access to the extracellular phosphate.

This GRNN produced data can further be used to understand the spatial and temporal dynamics of phenotypic clustering of gene expressions which is important in predicting and diagnosis of diseases [49]. Fig. 15 shows the phenotypic variation of WD biofilm in LP. Fig. (a)a shows the number of cluster variations over the first 30 \(TSs\) when the significant phenotypic changes of the biofilm is evident. At around \(TS=9\) and \(TS=10\), the bacterial cells have the most diverse expression patterns due to the highest extracellular nutrient penetration (can be seen in Fig. (a)a) to the biofilm and inter-cellular communications. Here we use four TSs (\(TS=5\) - Fig. (b)b, \(TS=15\) - Fig. (c)c, \(TS=23\) - Fig. (d)d and \(TS=30\) - Fig. (e)e) to analyze this phenotypic differentiation. Each pair of Uniform Manifold Approximation and Projection (UMAP) plot and diagram of cell locations of each cluster explain how nutrient accessibility contribute to the phenotypic clustering. Although at \(TS=5\) (Fig. 15) the average number of clusters is over four, there are only two major clusters that can be observed with higher proportions, as shown in the pie chart. Among the two major clusters (blue and green) of Fig. (b)b, the bacteria in the blue cluster can mostly be found in the

Fig. 14: Illustration of GRNN information flow variations concerning the particular positions of cells within the biofilm. We selected four cells at a) [10, 10, 0] – close to the attached surface, b) [10, 10, 5]- close to the periphery, c) [7, 10, 13] – at the center and d) [3, 15, 0] – close to the attached surface and the periphery of the biofilm.
center of the biofilm, while the green cluster cells are close to the periphery. Fig. 15c and Fig. 15d have more clusters as the nutrient accessibility among cells is high. In contrast, due to the lack of nutrients in the biofilm, a limited number of clusters can be seen in the biofilm after around \(TS=30\), which can be observed Fig. 15e.

## V Conclusion

The past literature has captured the non-linear signal computing mechanisms of Bacterial GRNs, suggesting underpinning NN behaviors. This study extracts a GRNN with summarized multi-omics gene expression regulation mechanisms as weights that can further analyze gene expression dynamics, design predictive models, or even conduct _in-vivo_ computational tasks. We used _P. aeruginosa_ single species biofilm as a use case and extracted relevant gene expression data from databases such as RegulomePA and transcriptomic data from databases including GEO. Due to the complexity of the GRN and expression dynamics, we only considered a smaller sub-network of the GRN as a GRNN that is associated with QS, iron and phosphate inputs, and pyocyanin production. Considering this GRNN, we modeled the computation process that drives cellular decision-making mechanism. As bacteria live in ecosystems in general where intra-cellular communication play a significant role in cellular activities, an _in-silico_ biofilm is modeled using GNN to further analyze the biofilm-wide decision-making. A comparison between the GRNN generated data and the transcriptomic data from the literature exhibits that the GRN behaves similarly to a NN. Hence, this model can explore the causal relationships between gene regulation and cellular activities, predict the future behaviors of the biofilm as well as conduct bio-hybrid computing tasks. Further, in the GRNN extraction phase, we were able to identify the possibility of modeling more network structures with various number of input nodes, hidden layers, and output nodes. In addition, GRN components including auto regulated genes and bidirectional intergenic interactions hints the possibility of extracting more sophisticated types of GRNNs such as Recurrent NN and Residual NN in the future. The idea of extracting sub-networks as NNs can lead to more intriguing intra-cellular distributed computing. Further, this model can be extended to multi-species ecosystems for more advanced predictive models as well as distributed computing architectures combining various NNs.
