Biological Factor Regulatory Neural Network
Genes are fundamental for analyzing biological systems and many recent works
proposed to utilize gene expression for various biological tasks by deep
learning models. Despite their promising performance, it is hard for deep
neural networks to provide biological insights for humans due to their
black-box nature. Recently, some works integrated biological knowledge with
neural networks to improve the transparency and performance of their models.
However, these methods can only incorporate partial biological knowledge,
leading to suboptimal performance. In this paper, we propose the Biological
Factor Regulatory Neural Network (BFReg-NN), a generic framework to model
relations among biological factors in cell systems. BFReg-NN starts from gene
expression data and is capable of merging most existing biological knowledge
into the model, including the regulatory relations among genes or proteins
(e.g., gene regulatory networks (GRN), protein-protein interaction networks
(PPI)) and the hierarchical relations among genes, proteins and pathways (e.g.,
several genes/proteins are contained in a pathway). Moreover, BFReg-NN also has
the ability to provide new biologically meaningful insights because of its
white-box characteristics. Experimental results on different gene
expression-based tasks verify the superiority of BFReg-NN compared with
baselines. Our case studies also show that the key insights found by BFReg-NN
are consistent with the biological literature.
Xinnan Dai, Caihua Shan, Jie Zheng, Xiaoxiao Li, Dongsheng Li
# Biological Factor Regulatory Neural Network

###### Abstract.

Genes are fundamental for analyzing biological systems and many recent works proposed to utilize gene expression for various biological tasks by deep learning models. Despite their promising performance, it is hard for deep neural networks to provide biological insights for humans due to their black-box nature. Recently, some works integrated biological knowledge with neural networks to improve the transparency and performance of their models. However, these methods can only incorporate partial biological knowledge, leading to suboptimal performance. In this paper, we propose the Biological Factor Regulatory Neural Network (BFReg-NN), a generic framework to model relations among biological factors in cell systems. BFReg-NN starts from gene expression data and is capable of merging most existing biological knowledge into the model, including the regulatory relations among genes or proteins (e.g., gene regulatory networks (GRN), protein-protein interaction networks (PPI)) and the hierarchical relations among genes, proteins and pathways (e.g., several genes/proteins are contained in a pathway). Moreover, BFReg-NN also has the ability to provide new biologically meaningful insights because of its white-box characteristics. Experimental results on different gene expression-based tasks verify the superiority of BFReg-NN compared with baselines. Our case studies also show that the key insights found by BFReg-NN are consistent with the biological literature.

biological factor, gene regulatory network +
Footnote †: [leftmargin=*] *The work was done during the author’s internship with Microsoft.

+
Footnote †: [leftmargin=*] *The work was done during the author’s internship with Microsoft.
the functions of genes or proteins, and ultimately give insights into the mechanism of larger living systems. Particularly, BFReg-NN is a neural network with the following architecture. For one thing, each neuron is mapped into a biological factor (e.g., a specific gene or protein), and arranged level by level based on the hierarchy of biological concepts, such as genes, proteins, pathways, biological processes, and so on. For another, since biological factors regulate each other, edges between neurons (and hyperedges among neurons) are set to reflect the existence of these regulations. In such a manner, edges also model biological meanings.

Further, we apply two different operations to simulate all the types of interactions in BFReg-NN. Inside one layer, genes regulate each other and create feedback loops to form cyclic chains of dependencies in the regulatory network. Therefore, graph neural network styled operations are suitable to build the gene reactions in the k-hop neighborhood, and obtain the "steady state" of genes. It is the same for proteins in PPI. In the layer of pathways, it can be regarded as a hypergraph where each hyperedge is a pathway containing multiple proteins. Accordingly, the hypergraph neural network is used to aggregate and balance the information of each pathway. Across the layers, we need to imitate the material transformation (e.g., genes translate protein). Thus, we adopt deep neural network styled operations to map the relations.

Moreover, BFReg-NN is flexible to explore new biological knowledge by trying to add non-existing edges in the network. The reason is that some edges are still hard to be discovered due to technological limitations and rewiring phenomenons in individual cells. If the prediction performance is improved by added edges, it is more likely that these edges exist. We illustrate BFReg-NN simulates on genome-scale cell system as an example in Figure 1.

The advantages of our proposed model include: (1) Compared with previous works, BFReg-NN merges with the structural biological knowledge in cell systems, including hierarchical relations (e.g., genes-proteins, proteins-pathways mappings), and regulatory relations among certain factors, such as GRN and PPI. Therefore, it could imitate how different biological factors work inside a cell. (2) The model of BFReg-NN is transparent and interpretable, as each neuron and edge has its corresponding biological meaning. Thus, the learned model weights give evidence of which biological parts are activated and which biological products are generated, leading to the final prediction. (3) By adding new edges between neurons, BFReg-NN not only achieves better performance in downstream tasks, but also has the potential to complete undiscovered biological knowledge. Traditional knowledge completion methods for biological domains (e.g., link prediction by knowledge bases/graphs) suffer from imbalanced data problems (Brock et al., 2017). BFReg-NN utilizes the gene expression data, which reflects the real cell states, and thus obtains more reliable results.

In the experiment, we show the effectiveness of BFReg-NN on several kinds of biological tasks, including missing gene expression value prediction, cell classification, future gene expression value forecasting and cell trajectory simulation. We also test the knowledge completing ability of BFReg-NN by the recall of the existing biological knowledge. Finally, we do case studies for newly discovered knowledge. The results demonstrate that BFReg-NN provides biologically meaningful insights.

## 2. Related Work

**Gene expression and its applications:** By RNA sequencing, it is easy to obtain gene expression which is a value to represent the amount of gene transcripts from a DNA fragment (Kraus et al., 2017). It has been used in a variety of biological applications, including single-cell analysis (Shen et al., 2017; Wang et al., 2018), disease diagnosis (Shen et al., 2018) and drug discovery (Shen et al., 2018). But most of these models lack transparency and ignore the existing biological knowledge.

**Knowledge graph enhanced downstream tasks:** The emergence of knowledge bases/graphs has led to enhancing the performance in many fields of computer science, such as computer vision and natural language processing (Shen et al., 2018; Wang et al., 2018; Wang et al., 2018). Similarly, knowledge graphs also have been widely used for specific biological tasks such as cancer diagnosis in recent years (Kraus et al., 2017; Wang et al., 2018). Some methods (Shen et al., 2018; Wang et al., 2018; Wang et al., 2018) used the Gene Ontology (GO) knowledge, which defines GO terms (e.g., molecular function, cellular component, biological process) and builds the neural network architecture based on term relations. However, the network is too sketchy to simulate the complex gene/protein reactions in the cell, leading to suboptimal performance. OntoProtein (Wang et al., 2018) embedded the gene ontology knowledge in pre-training to improve the performance of several protein downstream tasks. But it limited the interpretation ability because of transforming knowledge into embeddings. P-Net (Kraus et al., 2017) modeled the relations of gene-pathway and pathway-biological process as a neural network to diagnose prostate cancer. GLRP (Gulman et al., 2017) assigned genes to the PPI network to predict breast cancer, and then explained the important genes by Layer-wise Relevance Propagation. Although they mitigate the black-box issues, they all use a small part of biological knowledge, and they cannot explore new knowledge from gene expression data.

**Knowledge complement:** Our work is related to knowledge complement (Kraus et al., 2017; Wang et al., 2018). In particular, lots of work on biological knowledge complement has been done recently (Wang et al., 2018; Wang et al., 2018). Much attention is paid to predicting the relation between specific biological factors, which can be regarded as a link prediction problem. For example, (Kraus et al., 2017) classified whether ligand and receptor proteins interact based on their molecular structures. (Wang et al., 2018) predicted drug-target interaction (DTI) by learning representations of drugs and targets from the KEGG database. (Kraus et al., 2017) discovered new drugs for diseases by embedding the drug-gene-disease database in multiple relations. However, existing biological knowledge is collected in imbalance, which causes little biological meaningful predictions (Brock et al., 2017). In our work, we expect discovered knowledge from gene expression data, which can interpret the cell state by biological factor interactions and avoid imbalanced situations.

## 3. Biological Background

The biological system is modeled by a complex network that connects many biologically relevant entities to work together to perform one or more particular functions. It could be at the organ/tissue scale, such as the nervous system, or the integumentary system. On the micro/nanoscopic scale, examples include cells, organelles, and so on. In this work, we focus on the simulation of the biological system in a cell at a genome-scale.
Thanks to the development of sequencing technologies, it is easy and cheap to obtain an amount of gene expression data to build genome-scale analysis. Starting from genes, we model the cell genome-scale system by single-omic level (intra-level) regulations and different onice level (inter-level) mappings (Becker et al., 2015). Among different ionic levels, defined as inter-level in BFReg-NN, genes are first transcribed from DNA and then translate to proteins. Then, pathways integrate individual genes or protein products to perform certain cell functions. To model single-omic level regulations, we define the intra-level interactions. Genes are regulated by each other, known as GRN, which means the gene expression values are governed to be inhabited or activated by association molecule products, such as RNA. The interactions among proteins are similar, where they inhabit, activate, or combine with others to influence the expression values (or protein abundances) in cells, called PPI. Notice that in general, the relations among factors are universal for all the cells, thus the biological knowledge could adapt to different types of cells. But sometimes the biological factor interaction would rewire in specific cells, and thus the extra knowledge discovered and verified is also important (Kang et al., 2017).

## 4. Methodology

Assume we have gene expression data \(\mathbf{x}\in\mathbb{R}^{n}\) show the state of a cell, where \(n\) is the number of genes. Besides, we also define the neural network architecture of BFReg-NN according to structural biological knowledge with the intra-level regulations \(\mathcal{A}\) and the inter-level mappings \(\mathcal{M}\). After training BFReg-NN, we obtain the hidden embedding \(\mathbf{H}_{i}=\text{BFReg-NN}(\mathbf{x},\mathcal{A},\mathcal{M})\) for each gene \(i\), and do predictions for downstream tasks.

The output format of predictions could be a single value \(y=f(\mathbf{H})\), or multiple values \(\mathbf{x}=f(\mathbf{H})\). Further, we also handle the time-series output. Given the gene expressions of cells at the \(t_{0}\) time, we aim to predict the gene expressions \(\hat{\mathbf{X}}\in\mathbb{R}^{n\times T}=f(\mathbf{H})\) in the following \(t_{1},\ldots,T\) time steps. Finally, BFReg-NN could explore new insights by adding more edges \(\mathcal{A}^{\prime}\).

In the following sections, we first extract \(\mathcal{A}\) and \(\mathcal{M}\) from the existing biological databases/graphs. Then, we propose two versions of BFReg-NN. Finally, we introduce how to apply BFReg-NN to different downstream tasks.

### Biological knowledge databases/graphs

Based on the existing biological knowledge, we first divide biological factors into different levels, \(\mathcal{L}=\{\text{Gene},\text{Protein},\text{Pathway},\ldots\}\). At each level, regulatory relations between factors in the intra-level are formulated as a matrix set \(\mathcal{M}=\{\mathbf{A}_{\text{Gene}},\mathbf{A}_{\text{Protein}},\mathbf{A}_{ \text{Pathway}},\ldots,\mathbf{A}_{L}\}\), where \(\mathbf{A}_{\text{Gene}}\) could be gene relations and defined by GRN, and \(\mathbf{A}_{\text{Protein}}\) is decided by PPI. \(\mathbf{A}_{\text{Pathway}}\) is a little special because it is a hypergraph, where each edge could connect more than two nodes, called the hyperedge. The proteins in a hyperedge propagate the information to influence each other directly.

The values in \(\mathbf{A}_{I}\) are binary to represent the existence of relations. We also identify the binary mapping matrixes \(\mathcal{M}=\{\mathbf{M}_{1},\mathbf{M}_{2},\ldots,\mathbf{M}_{L-1}\}\) from level \(l\) to its upper level \(l+1\) as the inter-level interaction. The dimension of \(\mathbf{M}_{I}\) is dependent on the factor numbers of the two levels. We set \(\mathbf{M}_{1}\) to map between genes and proteins, and \(\mathbf{M}_{2}\) as a

Figure 1. The pipeline of BFReg-NN. We build the hierarchical biological network inspired by the cell system, modeling several levels to separate mRNA, protein, pathway and phenotype, and simulating biological factor interactions both at the intra-level and inter-level. We obtain the gene expression value from RNA sequencing as input, to predict the property of the phenotype as output. Taking the cancer dataset NCI-60 as an example. We know that the regulatory activations start from gene HIF1A and TP53 to gene AURKA and CDKN1A at the gene level (Kang et al., 2017). And at the protein level, the translated proteins Q16665 and P04637 stimulate P17858 and O14965 respectively (Kang et al., 2017). We mark those relations as existing biological knowledge by the black solid lines. Multiple genes with their products consist of pathways to drive cells to different types. For example, PFKL and HIF1A are activated in the HEPG2 cell line pathway, which leads to a type of liver cancer (Kang et al., 2017). Besides, new regulatory relations would be identified at each level shown as dotted lines. After training, we learn the good representations of biological factors, and employ them for different downstream tasks with various output formats.
direct mapping between protein-level and pathway-level. \(\mathcal{A}\) and \(\mathcal{M}\) both decide the architecture of the neural network, including neurons and links between neurons.

### Basic BFReg-NN model

Given the gene expression of a cell \(\mathbf{x}\), we first utilize an embedding layer to encode each gene independently, where \(\mathbf{H}_{i}^{0,0}=\text{Emb}(\mathbf{x}_{i})\). The embedding layer lets genes share the same parameter MLP layer to obtain the gene expression representation. As the gene expression is in a floating form with a little measurement error, the embedding layer is utilized to reduce this error and thus enhance the performance of the model. Then we obtain \(\mathbf{H}^{0,0}\in\mathbb{R}^{n\times d}\) as the input to neurons at the first level, i.e., gene neurons.

We have \(\mathbf{A}_{l}\) as intra-level relations between neurons at level \(l\). Inspired by graph neural networks, we use the message passing mechanism to update \(\mathbf{H}^{l,k}\), which is the hidden representation of the neuron in the \(k\)-hop at level \(l\). The formulation is:

\[\mathbf{H}_{i}^{l,k+1}=\text{update}\left(\sum_{j\in\mathbf{A}_{l}(i)}\text{ message}\left(\mathbf{H}_{i}^{l,k},\mathbf{H}_{j}^{l,k}\right),\mathbf{H}_{i}^{l,k} \right). \tag{1}\]

The message function is to generate the message from neuron \(j\) to neuron \(i\), where \(\mathbf{A}_{l}(i)\) decides which neurons are neighbors for neuron \(i\). Then the update function is to update the embedding of neuron \(i\) using the obtained messages and the previously hidden embedding. \(k\in[0,K-1]\) is the number of hops to determine that neuron \(i\) is influenced by other neurons in the \(K\)-hop neighborhood. Batch normalization is conducted at the end of each level. In detail, we utilize GAT-styled (graph attention) (Sutskever et al., 2017) to compute the attention weight for each message, and weighted sum up all the messages to update. As for operations in hypergraphs, we apply HGNN-styled (hypergraph neural network) (Bahdan et al., 2017) message and update functions. Specifically, the hypergraph can be represented by an incidence matrix \(\mathbf{R}\in\mathbb{R}^{|\mathbf{V}_{\mathbf{R}}|\times|\mathbf{E}_{\mathbf{R }}|}\), where \(\mathbf{R}_{ij}=1\) means the node \(i\) belongs to the hyperedge \(j\) (i.e., the protein \(i\) is contained in the pathway \(j\)). The equation is

\[\mathbf{H}^{l,k+1}=\sigma(\mathbf{D}_{\mathbf{R}}^{-1}\mathbf{R}\mathbf{B}_{ \mathbf{R}}^{-1}\mathbf{R}^{T}\mathbf{H}^{l,k}\mathbf{W}^{l,k})\]

where \(\mathbf{D}_{\mathbf{R}}\) and \(\mathbf{B}_{\mathbf{R}}\) are two corresponding degree matrices for normalization, and \(\mathbf{W}^{l,k}\) is the learnable matrix.

Due to multiple-level interactions in biological systems, the embedding is also learned level level by. Since \(\mathbf{M}_{l}\) is inter-level relations between level \(l\) and \(l+1\), we utilize the masked deep neural network to update the initial representation \(\mathbf{H}^{l+1,0}\) at next level

\[\mathbf{H}^{l+1,0}=\text{activation}\left(\left(\mathbf{M}_{l}\odot\mathbf{W }^{l,K}\right)\mathbf{H}^{l,K}+\mathbf{b}^{l}\right).\]

\(\mathbf{H}^{l,K}\) is the output after batch normalization. The element-wise multiplication \(\mathbf{M}_{l}\odot\mathbf{W}^{l,K}\) ensures that non-existing relations are not used for updating. \(\mathbf{W}^{l,K}\) and \(\mathbf{b}^{l}\) are learnable parameters.

### Enhanced BFReg-NN model

Here the enhanced version is introduced by adding new edges in \(\mathcal{A}=\{\mathbf{A}_{\text{Gene}},\mathbf{A}_{\text{Protein}},\mathbf{A}_{ \text{Pathway}},\ldots\}\). It can explore new biological insights and improve performance simultaneously.

Existing biological knowledge is detected by biological technology to reflect the implicit relations among factors. However, some knowledge is still hard to be discovered due to technological limitations and rewiring phenomenons in individual cells. Therefore, we spilt the interaction into two types. One is the universal regulation, supported by existing knowledge \(\mathcal{A}\). The other is local interaction, inferring the new biological knowledge or rewiring in individual cells, but now hidden in the non-existent edges of \(\mathcal{A}\). Instead of a binary matrix \(\mathbf{A}_{l}\) used in the basic model, we use \(\mathbf{A}_{l}\) to constrain the learnable matrix \(\mathbf{A}^{\prime}_{l}\) to discover new knowledge. As universal regulations are verified by biological methods, we use them the same as the basic model. For non-existent edges, we reweigh it by a small value \(0\leq\alpha<1\) due to it being less convincing. Thus, the edge intensity based on two types of knowledge is modified as:

\[\mathbf{A}^{\prime}_{l}=\begin{cases}\alpha^{l}_{ij^{\prime}}&\text{universal regulation that already exists in }\mathbf{A}_{l}\\ \alpha^{0}_{ij^{\prime}}&\text{local interaction that is ignored in }\mathbf{A}_{l} \end{cases} \tag{3}\]

where \(\alpha^{l}_{ij}=\sigma(\text{MLP}(\text{concat}[\mathbf{H}_{l}^{l},\mathbf{H}_{ j}^{l}]))\). In the enhanced model, we not only learn the neuron embeddings, but also utilize these embeddings with an MLP transformation to infer the intensity of the hidden interaction between neuron \(i\) and \(j\). Then we update the representations by

\[\mathbf{H}^{l,k+1}=\sigma\left(\mathbf{A}^{\prime}_{l}\mathbf{H}^{l,k}\mathbf{ W}^{l,k}\right). \tag{4}\]

After the model is trained to converge, we obtain the learned weights for non-existence edges, which provides insights for new knowledge and the rewiring phenomenon. We sort the weights \(\mathbf{w}^{l}_{ij}\) to identify the candidates which deserve verification by biological experiments. Here we use a simple edge weight modification method instead of gated edges implemented by the gumble-softmax function or advanced graph structure learning algorithms (Sutskever et al., 2017; Wang et al., 2018). The reason is that the biological knowledge is not sparse, and even dense in some core genes. Thus, it does not satisfy the sparse and low-rank requirements. Besides, we also implement a neural network where the existence of new edges is gated by the gumble-softmax function, but it does not work in the experiments. Our simple method does not add much extra computation cost while achieving great effectiveness in biological tasks.

### Downstream Tasks

After we obtain the final embedding \(\mathbf{H}^{l,K}\), we could conduct different types of downstream tasks, whose output format could be one-dimensional, multiple-dimensional, or time-series. Here we illustrate them with three specific tasks, missing gene expression value prediction, future gene expression value forecasting, and cell classification.

**Missing gene expression value prediction:** Suppose we have a gene expression vector \(\mathbf{x}\in\mathbb{R}^{n}\) for \(n\) genes, measured by the single-cell sequencing technology. Some parts in the vector are zeros or in a low value because of dropout events (K
**Cell classification:** This task is to classify the cell types using gene expression values \(\mathbf{x}\), which is important to determine the situation of tissues or patients. Because cell type is the observable result of a biological system, we simulate each gene to pass the multiple levels of transformation to infer the property of the cell. The prediction is computed by \(\hat{y}=\mathrm{MLP}(\mathbf{H}^{L,K})=\mathrm{MLP}(\mathrm{BFReg-NN}(\mathbf{x}))\). It is a multi-class classification task and we employ a cross-entropy loss.

**Future gene expression value forecasting:** The single-cell data is efficient to analyze the cell response under different drugs. However, it is expensive to collect the data at different time steps and draw the development trend of cells. Thus, we model the cell response by future gene expression value prediction. In other words, given the gene expression data \(\mathbf{x}\in\mathbb{R}^{n}\) at \(t_{0}\), we aim to forecast the gene expression data \(\bar{\mathbf{X}}\in\mathbb{R}^{n\times T}\) in the following \(t_{1},\dots,T\) time steps. We use two backbones, MLP and LSTM. MLP predicts the gene expression in the future time steps simultaneously, and the equation is \(\hat{\mathbf{X}}=\mathrm{MLP}(\mathrm{BFReg-NN}(\mathbf{x}))\). LSTM models the dynamic values step by step, where it takes the last time output as the next step input, and the equation is \(\hat{\mathbf{x}}^{t}=\mathrm{LSTM}(\mathrm{BFReg-NN}(\hat{\mathbf{x}}^{t-1}))\). We employ MSE as the loss function.

## 5. Further Applications

BFReg-NN has the potential to be plugged into more sophisticated models, such as generative models and pre-learning.

### Cell trajectory simulation

In this task, we regard cells as a population to understand the population dynamics. In other words, our goal is to simulate the evolution of the distribution of gene expression over time. Since the single-cell sequencing technology makes the cells die after measurement, we can only know the distribution of gene expression of a part of cells at a certain time step, i.e., \(\{x_{1,t_{i}},x_{2,t_{i}},\cdots,x_{N_{i},t_{i}}\}\) where \(i=\{1,2,\cdots,T\}\). Based on previous papers (Srivastava et al., 2017; Wang et al., 2018), we also utilize continuous normalizing flows (Beng et al., 2018) as the generative model to learn an underlying differentiation landscape from scRNA-seq data. The main difference is that our proposed BFReg-NN is utilized as \(f_{0}\) instead of a traditional DNN. Because the architecture of BFReg-NN is built on biological factors and their interactions are also regarded as differential equations, BFReg-NN has a better inductive bias to model complex relationships and further improves the performance.

Specifically, assume that \(P_{t}(\mathbf{x})\) is the distribution of gene expression of cells at time step \(t\), and \(\mathbf{x}(t)\) is the gene expression of a cell drawn from \(P_{t}(\mathbf{x})\). Let \(\frac{\mathrm{d}\mathbf{x}(t)}{\mathrm{d}t}=f_{0}(\mathbf{x}(t),t)\) be a differential equation describing a continuous-in-time gene transformation of \(\mathbf{x}\). We seek to train \(f_{0}\) satisfying

\[\mathbf{x}(t_{i}) \sim P_{t_{i}}(\mathbf{x})\quad i=0,1,2,\cdots,T\] \[\frac{\mathrm{d}\mathbf{x}(t)}{\mathrm{d}t} =f_{0}(\mathbf{x}(t),t)\text{ for }t\in[t_{0},t_{T}].\]

According to the characteristic of continuous normalizing flows, we have

\[\log P_{t_{T}}(\mathbf{x}(t_{T}))=\log P_{t_{0}}(\mathbf{x}(t_{0}))-\int_{t_{0}}^{t_{T }}\mathrm{Tr}(\frac{\partial f_{0}(\mathbf{x}(t),t)}{\partial\mathbf{x}(t)})\mathrm{d}t.\]

Therefore, given a batch of empirical samples \(\{x_{0,t_{i}},x_{1,t_{i}},\cdots,x_{N_{i},t_{i}}\}\) where \(i=\{0,1,\cdots,T\}\), we can first map the distribution at the time \(t_{i}\) into the distribution at the time \(t_{i-1}\). Then we compare the distribution difference to make it as the same as possible. We use the wasserstein distance as the loss to evaluate the difference.

Then we utilize BFReg-NN to instantiate \(f_{0}\). Assume that we only have one layer, such as gene layer. We first embed \(\mathbf{x}\) by \(\mathbf{H}^{0}=\mathrm{Emb}(x)\), and change the message passing for genes

\[\mathbf{H}^{k+1}=\mathbf{H}^{k}+\int_{t_{i}}^{t_{k}+\alpha_{T}}f_{0}(\mathbf{ H}(t),t)\mathrm{d}t \tag{5}\]

where \(f_{0}(\mathbf{H}(t),t)=\sigma(A^{\prime}\mathbf{H}(t)\mathbf{W}(t))\mathbf{U}(t)\), and \(\mathbf{W}(t)\) and \(\mathbf{U}(t)\) are produced by the hypernetwork (Kang et al., 2018) to make them change over time. When we have different \(L\) layers in BRReg-NN, which mean \(f_{0}\) is a piecewise function based on the value of \(t\). One straightforward and effective approach is to define \(f_{0}\) as:

\[\begin{cases}f_{0_{1}}\leftarrow\sigma(A^{\prime}_{1}(t)\mathbf{W}^{1}(t)), \ t\in[t_{i},t_{i}+\frac{L}{L}\Delta_{t})\\ f_{0_{2}}\leftarrow\sigma((\mathbf{M}_{1}\odot\mathbf{W}^{2}(t))\mathbf{H}(t) +\mathbf{b}^{2}(t)),\ t\in[t_{i}+\frac{L}{L}\Delta_{t},t_{i}+\frac{L}{L} \Delta_{t})\\ \cdots\\ f_{0_{L}}\leftarrow\sigma(A^{\prime}_{L}\mathbf{H}(t)\mathbf{W}^{L}(t)) \mathbf{U}(t),\ t\in[t_{i}+\frac{L-1}{L}\Delta_{t},t_{i}+\Delta_{t}].\end{cases}\]

\(f_{0_{1}}\) imitates the interactions in gene layer, \(f_{0_{2}}\) simulates the mapping from genes to proteins, and then the following \(f_{0_{i}}\) imitate the reactions in higher layers.

### Pre-training

BFReg-NN is also suitable in pre-training and fine-tuning framework. For example, we utilize the missing gene expression prediction task to pre-train a BFReg-NN. Then we freeze the parameters of BFReg-NN except for the last MLP layers, and fine-tune the last MLP layers in downstream tasks, such as forecasting future gene expression and cell classification.

## 6. Experiments

In this section, we first describe the experimental setting for each task, including datasets, evaluation metrics, network architecture, baselines and training details (Section 6.1). Then we present the main results in Section 6.2. The cell trajectory simulation and pre-training are conducted in Section 6.3 and Section 6.4. Further, Section 6.5 shows the ability of BFReg-NN to discover the new biological knowledge. Finally, we did some ablation study to evaluate the each part of BFReg-NN in Section 6.6. All the experiments are run on a single A100 GPU. The code and toy datasets are in [https://github.com/DDigimon/BFReg](https://github.com/DDigimon/BFReg).

### Setup

#### 6.1.1. Datasets

Here we describe different datasets for all the tasks mentioned in Section 4 and 5.

**Missing gene expression value prediction:** In this task, we collect 10 cell lines (BT20, HS578T, LNCaP, A549, MCF7, MCF10A, MDAMB231, PC3, SKBR3 and A375) from the L1000 dataset. Because of the limitation of sequencing technology, the L1000 dataset

\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline Dataset & \#Sample & rGene & GBN knowledge & PPI knowledge & Pathways & fLabel \\ \hline GSE & 4688 & 418 & 172 & 136 & 89 & 8 \\ lung & 1642 & 382 & 148 & 332 & 121 & 11 \\ muscle & 1051 & 401 & 317 & 401 & 119 & 6 \\ trachea & 3587 & 326 & 109 & 207 & 110 & 5 \\ diaphragm & 853 & 355 & 239 & 299 & 115 & 5 \\ \hline \hline \end{tabular}
\end{table}
Table 1. Statistics for cell classification datasets
does not have the golden standard values for non-landmark genes. Thus, we employ the pre-processed data based on (Srivastava et al., 2017). Besides, we also ignore the time point information in cell lines. In each cell line, we have 1482 cell samples and 714 genes. We also collect 541 and 2305 existing knowledge edges for GRN and PPI respectively. We randomly mask the gene expression value with a 60% probability, and predict these missing values. The cell samples are randomly split by 60/20/20% as training/validation/test data. The evaluation metric is mean square error (MSE), which means the smaller the better.

**Cell classification:** We gather 5 datasets to predict cell types, and data statistics are summarized in Table 1. The first one is GSE756888 (breast cancer) collected from (Brock et al., 2017). The following four datasets are about organs, including muscle, diaphragm, lung and trachea, obtained from (Brock et al., 2017). We pre-process the data by deleting the cells and genes which have a larger ratio of zeros. The remaining data are summarized in Table 1. We also split the data same as the missing value task, and evaluate the results by macro AUC.

**Future gene expression value forecasting:** For future value prediction, we choose the breast cancer dataset from 2019 Dream Challenge (Deng et al., 2019). There are 44 cell lines and 5 treatments. Different treatments mean using various drugs on cells to test the drug effects. So we have 220 situations. In each situation, the detailed process to obtain time-series gene expression data is: all the cells are first treated and then divided into several groups. At (i)regular time step, a group of cells is killed and their gene expression is measured. There are seven timestamps, and about hundreds of cells are measured at each timestamp. To avoid noise caused by this measurement on each individual cell, we calculate the median for each gene. Thus, the format of data is \(\mathbf{X}\in\mathbb{R}^{N_{\text{treatment}}\times N_{\text{Netall,line}} \times N_{\text{mustemap}}\times N_{\text{gene}}}\).

Since the data are quantified at the proteomic level, we have 37 biomarkers as genes and 38 related PPIs as the existing knowledge, which is a small protein-protein interactions graph where nodes are phosphorylated genes. Different from the previous two tasks, we use 5 treatments as training data, and the remaining one as test data. The 5-fold cross-validation is done. We use MSE and Pearson correlation coefficient (PCC) to evaluate average performance.

**Cell trajectory simulation:** In this task, we use the same dataset as future gene expression forecasting. Instead of calculating the median value, we consider the evolution of the distribution of gene expression. Thus, we sample \(512\) cells at each timestamp. The current format of data is \(\mathbf{X}\in\mathbb{R}^{N_{\text{treatment}}\times N_{\text{Netall,line}} \times N_{\text{mustemap}}\times N_{\text{gene}}}\). In each situation (i.e., a certain combination of cell lines and treatments), we use the first two timestamps as known data, and predict the following trend for the remaining timestamps. The weessentrian distance is utilized to evaluate the similarity between two distributions. The performance is also averaged by the number of predicted timestamps.

**Pre-training:** The dataset and evaluation metric are also the same as future value forecasting.

#### 6.1.2. Neural network architecture

The architecture is totally defined by existing knowledge. From open-source databases, we select Dorothea (Deng et al., 2019) as GRN construction, and Omnipath (Marcus et al., 2019) as PPI. As the pathway is specific in different cells, we collect pathways for each dataset by Enrich (Enrich, 2017), which provides knowledge of the hyperedge connection among genes and proteins, including LINCS_L1000_Ligand_Perturbations for missing value task and WikiPathways_2019_MOUSE(HUMAN) for cell classification task. For future value forecasting and Cell trajectory simulation, we utilize its corresponding relations (38 related PPIs) in the dataset.

#### 6.1.3. Baselines

We divide baselines into three types:

(1) There is no prior knowledge integrated into the model, and the input is only gene expression data.

* **MLP** is the simplest model for static prediction.
* **LSTM** is the typical RNN-based model for dynamic value prediction.
* **Random Forest** and **XGBoost** are two classical tree-based models.
* **Transformer**(Srivastava et al., 2017) is an advanced model which achieves great performance in many NLP tasks.
* **GatedNN** is an improved version of deep neural networks, where edges between neurons are gated by the gumbel-softmax function so the network could be sparse.

(2) The second type of baselines is to first learn a gene co-expression matrix as knowledge and then do the prediction. The co-expression matrix identifies which genes have a tendency to show a coordinated expression pattern, so it can be built by the similarity between gene expressions. We follow the process described in MLA-GNN (Marcus et al., 2019) to learn the knowledge and then apply two classical graph neural networks, **GCN**(Marcus et al., 2019) and **GAT**(Srivastava et al., 2017), to aggregate the information and do the prediction.

(3) We also compare the models using the existing biological knowledge. This knowledge could be graphs, such as GRN or PPI.

* **MLP+N2V** used node2vec (Deng et al., 2019) to learn node embeddings of the prior graph, combine them with gene expression data and employ MLP to predict.
* **GCN** and **GAT** directly utilized the prior graph. As there are different prior graphs (GRN or PPI), we report the best results for GCN and GAT with the most suitable graph.
* **DCeII**(Srivastava et al., 2017) used the Gene Ontology (GO) knowledge to build the neural network architecture.
* **P-NET**(Marcus et al., 2019) utilized protein-pathway relations as prior knowledge and apply DNNs to predict.
* **BFReg-NN** is our proposed method. There are basic and enhanced versions.

#### 6.1.4. Training Details

For optimizers, we use Adam for all the models. The hyperparameters are summarized in Table 9 in Appendix. We set a small dimension of hidden embeddings to avoid overfitting. For the prediction module, we set the number of hidden units with the best value obtained from the validation data. In addition, the update function could be implemented by the sum or concatenation operation. We set \(K=1\) for all the experiments.

The performance of BFReg-NN is influenced by the value of hyperparameter \(\alpha\). While \(\alpha\) controls the acceptance of newly discovered knowledge, it is dependent on whether the existing biological knowledge is suitable for the current dataset, as shown in Table 10 in Appendix. For example, in the future value forecasting task, \(\alpha\) is larger than in other tasks to get the best result because existing knowledge is incomplete in the dataset. Besides, in the missing value task, our model performs better when \(\alpha=0\) in the protein
level, which means BFReg-NN ignores the irrelevant level when discovering knowledge. We vary it from [0, 0.00001, 0.00005, 0.0001, 0.0005, 0.001, 0.005] and select the best value when the loss is the smallest in the validation dataset.

### Main Results

Here we describe the results from two aspects, static and dynamic. Missing gene expression prediction and cell classification are static tasks, while future gene expression prediction is dynamic tasks.

**Static task results:** We report the main results for missing value prediction and cell classification, summarized in Table 2 and Table 3 respectively. Firstly, in general, compared to no prior knowledge methods, biological knowledge improves the performance of models effectively, which means that modeling specific meanings for neural networks could improve the performance in biological tasks. Although Gated NN includes the gumble-softmax function which also is able to predict edges in a probability, it is hard to generate the existing biological knowledge. Besides, co-expression knowledge does not show outstanding performance because it is extracted from a small size of data and thus brings the noise. Some complex model also achieves promising results, such as Transformer, which indicates that building dense biological relations may be needed for the cell samples.

\begin{table}
\begin{tabular}{l l c c} \hline \hline
**Knowledge** & **Models** & **MSE** & **PCC** \\ \hline \multirow{6}{*}{no prior knowledge} & MLP & 0.0761 \(\pm\) 0.0038 & 0.9704 \(\pm\) 0.0018 \\  & LSTM & 0.0822 \(\pm\) 0.0028 & 0.8888 \(\pm\) 0.0034 \\  & Transformer & 0.0743 \(\pm\) 0.0013 & 0.9717 \(\pm\) 0.0007 \\  & Gated NN & 0.0743 \(\pm\) 0.0026 & 0.9717 \(\pm\) 0.0007 \\  & Random Forest & 0.1115 \(\pm\) 0.0021 & 0.9685 \(\pm\) 0.0006 \\  & XGBoost & 0.0837 \(\pm\) 0.0033 & **0.9753**\(\pm\) 0.0008 \\ \hline \multirow{6}{*}{co-expression} & GCN & 0.3838 \(\pm\) 0.0090 & 0.8293 \(\pm\) 0.0036 \\  & GAT & 0.3447 \(\pm\) 0.0236 & 0.8455 \(\pm\) 0.0114 \\ \cline{1-1} \cline{2-4}  & MLP+N2V & 0.0760 \(\pm\) 0.0043 & 0.9706 \(\pm\) 0.0018 \\ \cline{1-1}  & GCN & 0.0920 \(\pm\) 0.0053 & 0.9646 \(\pm\) 0.0020 \\ \cline{1-1}  & GAT & 0.0869 \(\pm\) 0.0044 & 0.9661 \(\pm\) 0.0017 \\ \cline{1-1} \cline{2-4}  & BFReg-NN (MLP+B) & 0.0732 \(\pm\) 0.0026 & 0.9719 \(\pm\) 0.0009 \\ \cline{1-1} knowledge & BFReg-NN (MLP+F) & **0.0724**\(\pm\) 0.0023 & 0.9724 \(\pm\) 0.0007 \\ \cline{1-1}  & BFReg-NN (LSTM+B) & 0.0825 \(\pm\) 0.0029 & 0.8906 \(\pm\) 0.0041 \\ \cline{1-1}  & BFReg-NN (LSTM+E) & 0.0819 \(\pm\) 0.0027 & 0.8907 \(\pm\) 0.0049 \\ \hline \hline \end{tabular}
\end{table}
Table 4. Results on future gene expression value forecasting (evaluated by MSE\(\downarrow\) and PCC’).

\begin{table}
\begin{tabular}{c l c c c c c c c c c} \hline \hline
**Knowledge** & **Models** & **BT20** & **HS578T** & **LNCAP** & **A549** & **MCF7** & **MCF10A** & **MDAMB231** & **PC3** & **SKBR3** & **A375** \\ \hline \multirow{3}{*}{no knowledge} & MLP & **0.314** & 0.344 & 0.320 & 0.326 & 0.337 & 0.325 & 0.378 & 0.317 & 0.336 \\  & Transformer & 0.396 & 0.377 & 0.406 & 0.383 & 0.396 & 0.399 & 0.418 & 0.402 & 0.387 & 0.407 \\  & Gated NN & 0.709 & 0.681 & 0.644 & 0.533 & 0.647 & 0.700 & 0.667 & 0.722 & 0.594 & 0.729 \\ \hline \multirow{3}{*}{co-expression} & GCN & 0.715 & 0.714 & 0.722 & 0.653 & 0.688 & 0.675 & 0.707 & 0.731 & 0.677 & 0.702 \\  & GAT & 0.719 & 0.715 & 0.723 & 0.655 & 0.693 & 0.669 & 0.720 & 0.721 & 0.676 & 0.716 \\ \hline \multirow{6}{*}{existing knowledge} & MLP+N2V & 0.351 & 0.378 & 0.380 & 0.349 & 0.368 & 0.350 & 0.358 & 0.392 & 0.343 & 0.354 \\  & GCN & 0.331 & 0.329 & 0.330 & 0.304 & 0.340 & 0.313 & 0.319 & 0.338 & 0.331 & 0.317 \\ \cline{1-1}  & GAT & 0.324 & 0.331 & 0.332 & 0.308 & 0.321 & 0.314 & **0.317** & 0.347 & 0.317 & 0.315 \\ \cline{1-1}  & P-NET & 0.349 & 0.335 & 0.340 & 0.340 & 0.321 & 0.352 & 0.336 & 0.354 & 0.322 & 0.331 \\ \cline{1-1}  & BFReg-NN(Basic) & 0.318 & 0.336 & 0.332 & 0.329 & **0.308** & 0.328 & 0.328 & 0.344 & 0.318 & 0.328 \\ \cline{1-1}  & BFReg-NN(Enhanced) & 0.318 & **0.317** & **0.320** & **0.297** & **0.308** & **0.306** & 0.320 & **0.335** & **0.307** & **0.308** \\ \hline \hline \end{tabular}
\end{table}
Table 2. Results on missing gene expression value prediction (evaluated by MSE\(\downarrow\)).

\begin{table}
\begin{tabular}{c l c c c c} \hline \hline
**Knowledge** & **Models** & **GSE** & **muscle** & **diaphrag** & **lung** & **trachea** \\ \hline \multirow{3}{*}{no prior knowledge} & MLP & 0.9122 & 0.8586 & 0.7881 & 0.8545 & 0.9387 \\  & Random Forest & 0.9688 & 0.8771 & 0.7161 & 0.8834 & 0.9273 \\  & XGBoost & 0.8867 & 0.8534 & 0.7549 & 0.8625 & 0.9193 \\  & transformer & 0.9476 & 0.8785 & **0.8717** & 0.8900 & 0.9242 \\  & Gated NN & 0.7910 & 0.7784 & 0.7110 & 0.8240 & 0.9120 \\ \hline \multirow{3}{*}{co-expression} & GCN & 0.6916 & 0.6424 & 0.5945 & 0.5822 & 0.6612 \\  & GAT & 0.7790 & 0.6882 & 0.6305 & 0.6312 & 0.6294 \\ \hline \multirow{3}{*}{existing knowledge} & MLP+N2V & 0.9064 & 0.8772 & 0.7501 & 0.8180 & 0.9404 \\  & GCN & 0.9285 & 0.8449 & 0.7896 & 0.8581 & 0.9421 \\ \cline{1-1}  & GAT & 0.9255 & 0.8470 & 0.8039 & 0.8246 & 0.9409 \\ \cline{1-1}  & P-NET & 0.9052 & 0.8654 & 0.7973 & 0.8425 & 0.9332 \\ \cline{1-1}  & DCell & 0.9482 & 0.7158 & 0.6731 & 0.7545 & 0.9322 \\ \cline{1-1}  & BFReg-NN(Basic) & 0.9476 & 0.8808 & 0.8420 & 0.8808 & 0.9376 \\ \cline{1-1}  & BFReg-NN(Enhanced) & **0.9693** & **0.8884** & 0.8509 & **0.8903** & **0.9446** \\ \hline \hline \end{tabular}
\end{table}
Table 3. Results on cell classification (evaluated by macro AUC\(\uparrow\)).
Secondly, BFReg-NN also outperforms the baselines using the existing knowledge. MLP+N2V fails in prediction because it uses the knowledge by embedding nodes implicitly, rather than in an explicit architecture. Compared with GCN and GAT, the basic BFReg-NN(B) could merge hierarchical knowledge from different levels of the biological system. Decll and P-NET built the sketchy networks and applied DNNs to predict, which cannot simulate the complex gene/protein reactions. Thus, their performance are also worse than BFReg-NN(B).

Thirdly, the enhanced BFReg-NN(E) learns important undiscovered knowledge from the gene expression data, and enhances the performance further. Overall, BFReg-NN(E) achieves the best results in most cases.

**Dynamic task results:** Results on future value prediction are presented in Table 4. Compared with directly employing MLP/LSTM to forecast the values, BFReg-NN followed with MLP/LSTM could improve the quality of embeddings and achieve better performance. GCN and GAT with co-expression matrix perform badly due to the small data size. When using existing knowledge, MLP+N2V, GCN and GAT still cannot obtain a good result because the knowledge is incomplete. For example, there are several isolated nodes in the dataset, and these methods fail in updating the embeddings of these nodes. Transformer and Gated NN merge the gene expression densely, so they obtain better results.

The best-performing models in this dataset are tree-based ensemble methods, so we apply Random Forest and XGBoost. Random Forest has a poor performance on both MSE and PCC. XGBoost has a little higher value on PCC but also a higher MSE, which means XGBoost could predict the trend but fail on the specific values.

Finally, the enhanced BFReg-NN(MLP+E) could add the extra discovered knowledge to mitigate the problem of incomplete prior knowledge, and learn the embeddings with a suitable architecture. Therefore, it has the lowest MSE value and a higher PCC value closing to the best one.

### Results for Cell Trajectory

To show the efficiency of the BFReg-NN model, we compare it with two baselines. MLP(S) indicates that genes do not interact with each other and the value is predicted by DNNs. MLP simply applied DNNs to aggregate all the genes and then simulate the trajectory one by one. Since the number of situations is large, we show the result for five treatments averaged by different cell lines in Table 5. In detail, the treatments include stimulating cells with EGF (epidermal growth factor), and four kinase inhibitors (iEFGR, iMEK, iFXC, iPI3K) combined with EGF.

BFReg-NN obtained the small westernitan distance compared with baselines on average. Specifically, BFReg-NN achieved significant results in EGF and iEGFR treatments. The reason is that gene value is constrained by their related genes during several time steps. BFReg-NN considered that suitable relationship among genes, which helps to predict the trajectory for each gene. MLP(S) regarded each gene as independent, so it cannot use the other gene information to make a better prediction. MLP supposed that all the genes have a strong influence on each other, which is against the truth that some genes have more significant influence than others. Thus, MLP has a suboptimal performance.

### Results for Pre-training

BFReg-NN could benefit from the pre-training and fine-tuning framework, and the results are shown in Table 6. The detailed experimental steps are: (1) Utilize the missing gene value prediction task to pre-train a BFReg-NN model on the breast cancer dataset; (2) Freeze the parameters of the BFReg-NN model except for the last MLP layers, and (3) Fine-tune the last MLP layers to do the future gene expression prediction. With the pre-training, the performance increases significantly on the MSE metric. Because biological knowledge is universal and similar in the cell across different tasks, it is reasonable that pre-trained BFReg-NN improves the other tasks' performance with fine-tuning. However, the standard deviation slightly rises, which means a more stable method to fine-tune the model should be discussed in the future.

### Knowledge completion by BFReg-NN

To demonstrate the ability of BFReg-NN to complete the knowledge, we evaluate it on the future value forecasting task. We conducted two type of experiments. One is to remove the existing edges and then see if these edges can be recalled by BFReg-NN. The other is to judge whether the newly discovered biological knowledge is also proved in the recent literature.

In the first experiment, for each gene, we remove its all the edges to make it an isolated node, and conduct the model to recover the knowledge. The discovered edges come out with weights, and we rank them to obtained the top-k list. We run BFReg-NN several times to avoid randomness. The frequency of an edge is computed by (the times that the edge is in the top-k list)/(the total runs). We select the top-20 interaction pairs with the highest frequency as discovered knowledge. The results are shown in Figure 2, where BFReg-NN achieves an average 0.4172 recall score for all the nodes.

We also take p.p53 gene as an example. The detailed recalled edges are shown in Figure 3 where the purple color marks the edges existed in the current network. All the existing edges for p.p53 are recalled in the top-20 list.

Secondly, some new knowledge is discovered by biological methods in recent breast cancer research, which is also found in our proposed BFReg-NN. We take the gene p53, a well-known tumor

\begin{table}
\begin{tabular}{l l l} \hline \hline Model & MSE & PCC \\ \hline BFReg-NN (MLP+E) & 0.0724 \(\pm\) 0.0023 & 0.9724\(\pm\) 0.0007 \\ Pre-trained BFReg-NN (MLP+E) & **0.0719\(\pm\)** 0.0042 & 0.9723\(\pm\) 0.0018 \\ \hline \hline \end{tabular}
\end{table}
Table 6. Performances of pre-trained BFReg-NN (evaluated by MSE\(\downarrow\) and PCC\(\uparrow\)).

\begin{table}
\begin{tabular}{l l l l l l l} \hline \hline Model & EGF & iEGFR & iMEK & iPKC & iPI3K & Avg. \\ \hline MLP(S) & 1.2829 & 1.2087 & 1.1471 & 1.1320 & 1.1362 & 1.2265 \\ MLP & 1.0325 & 1.0592 & 1.0893 & **1.0824** & **1.0827** & 1.0692 \\ BFReg-NN (E) & **0.8490** & **0.8753** & **0.9198** & 1.0830 & 1.1077 & **0.9670** \\ \hline \hline \end{tabular}
\end{table}
Table 5. Performances of cell trajectory simulation (evaluated by WD\(\downarrow\)).
suppressor, as an example to verify the biological meaning of discovered knowledge. We show the top-10 frequent interaction pairs including p53 in Figure 4. And then we find the cues in recent biological literature. Iododeoxyuridine (Idu) is the most frequent marker, which identifies cell phases, indicating whether the cell division is continued in breast cancer (Beng et al., 2017; Chen et al., 2018). S6 takes the second place, which is regulated in mTOR signaling by p53, and the level of S6 increases when p53 is insufficient (Krishnan et al., 2017). Besides, p53 also regulates the expression level of Her2, where the interaction frequency is about 0.15 in our prediction. As a kind of Her2-positive cancer, breast cancer could be treated by inhibiting p53 (Chen et al., 2018). After that, the regulatory relation between p53 and RB is in the existing biological knowledge dataset.

Overall, we notice that most of the existing edges can be recalled by our model, while BFReg-NN also adds knowledge that does not appear in the current database.

### Ablation Study

In this section, we provide a brief description of the effectiveness of each part in BFReg-NN. The ablation analysis for the cell classification task is shown in Table 8. The accuracy is gradually improved when adding a higher level of knowledge to the model in most cases. In trachea, GRN shows the best performance which indicates PPI may lead to an overfitting problem in this organ. In addition, we present the results of _no intra-level_ and _no inter-level_ situations. _No intra-level_ version removes PPI and GRN knowledge and maintains the hierarchical network; and _no inter-level_ merges the PPI and GRN into a large graph and deletes the hierarchical structure. The results show that hierarchical structure is more important than simply linking the factors together, which obeys the phenomenon that the biological factors in a cell are produced step by step.

In addition, we could choose different levels from the biological system as layers for BFReg-NN to adapt to different tasks. For example, the cell classification task aims to capture the cell state involving a variety of factors in a cell, and thus it needs multiple levels of biological knowledge to build the cell phenotype. In contrast, shown in Table 7, the missing value task focuses on gene expression value, which means single gene level knowledge is enough while protein or pathway level knowledge may bring the noise to predictions.

## 7. Conclusion

In this paper, we propose a generic framework BFReg-NN, which offers biological meanings to neurons and links between neurons, and imitates the whole cell system as a neural network at both intra-level and inter-level. We apply BFReg-NN to different downstream tasks and our experimental results demonstrate that BFReg-NN consistently outperforms baselines and discovers new biological meaningful insights. BFReg-NN provides a novel paradigm to merge cell sequencing data and biological knowledge, and it could be extended to more types of genomics data to simulate and understand complex biological systems in the future.

In the future, we plan to apply BFReg-NN in some real and important biological problems, such as predicting drug surgery effect for a special disease, and simulating stem cell development.

\begin{table}
\begin{tabular}{l c c c c c c c c c c} \hline \hline Knowledge & BT20 & HS578T & LNCaP & A549 & MCF7 & MCF10A & MDAMB231 & PC3 & SKBR3 & A375 \\ \hline GRN & **0.318** & **0.317** & **0.320** & **0.297** & **0.308** & **0.306** & **0.320** & 0.335 & **0.307** & **0.308** \\ GRN\&PPI & 0.343 & 0.342 & 0.338 & 0.326 & 0.339 & 0.333 & 0.343 & 0.357 & 0.338 & 0.344 \\ GRN\&Pathway & 0.335 & 0.335 & 0.326 & 0.314 & 0.324 & 0.319 & 0.331 & **0.332** & 0.327 & 0.325 \\ \hline \hline \end{tabular}
\end{table}
Table 7. Ablation study results on missing value prediction (evaluated by MSE \(\downarrow\)).