LLM4GRN: Discovering Causal Gene Regulatory Networks with LLMs - Evaluation through Synthetic Data Generation

## Abstract

Abstract Gene regulatory networks (GRNs) represent the causal relationships between transcription factors (TFs) and target genes in single-cell RNA sequencing (scRNA-seq) data. Understanding these networks is crucial for uncovering disease mechanisms and identifying therapeutic targets. In this work, we investigate the potential of large language models (LLMs) for GRN discovery, leveraging their learned biological knowledge alone or in combination with traditional statistical methods. We develop a task-based evaluation strategy to address the challenge of unavailable ground truth causal graphs. Specifically, we use the GRNs suggested by LLMs to guide causal synthetic data generation and compare the resulting data against the original dataset. Our statistical and biological assessments show that LLMs can support statistical modeling and data synthesis for biological research.

## 1 Introduction

Single-cell RNA sequencing (scRNA-seq) is a cutting-edge technology that enables the collection of gene expression data from individual cells. This approach opens up new avenues for a wide range of scientific and clinical applications. One crucial application of scRNA-seq data is the reconstruction and analysis of gene regulatory networks (GRNs), which represent the interactions between genes. GRN analysis can deepen our understanding of disease mechanisms, identify key regulatory pathways, and provide a foundation for the development of interventional gene therapies and targeted drug discovery.

Statistical causal discovery algorithms (Scheines et al., ; Zheng et al., ; Mercatelli et al., ; Brouillard et al., ; Lippe et al., ; Yu & Welch, ; Roohani et al., ) can reveal potential causal links between TFs and their target gene. However, they often lack robustness and are prone to detecting spurious correlations, especially in high-dimensional, noisy single-cell data. Furthermore, many of these approaches rely heavily on prior knowledge from curated databases (e.g., TRANSFAC (Wingender et al., ), RegNetwork (Liu et al., ), ENCODE (de Souza, ), BioGRID (de Souza, ), and AnimalTFDB (Hu et al., )), which frequently lack essential contextual information such as specific cell types or conditions, leading to inaccuracies in the inferred regulatory relationships (Zinati et al., ).

The recent advancements in and success of large language models (LLMs) have opened up new possibilities for their use in scientific discovery (Sheth et al., ; Lu et al., ; AI4Science & Quantum, ), including causal discovery (Kıcıman et al., ; Kasetty et al., ; Vashishtha et al., ; AI4Science & Quantum, ; Abdulaal et al., ; Khatibi et al., ).
Most of the above methods involve the refinement of the statistically inferred causal graph by LLM. However, LLMs excel at synthesizing vast amounts of heterogeneous knowledge, making them well-suited for tasks that require the integration of diverse datasets, such as constructing full causal graphs based on scientific literature (Sheth et al., ). In addition, LLMs have already been shown to perform in genomic data analysis (Märtens et al., ; Jin et al., ; Tang & Koo, ; Fang et al., ; Toufiq et al., ; Elsborg & Salvatore, ), including foundation models pre-trained for genomic tasks (Wang et al., ; Cui et al., ).

Inspired by recent advancements, we harness LLMs for inferring gene regulatory networks (GRNs) from scRNA-seq data. Specifically, we utilize LLMs either to generate complete GRNs (causal graphs) directly or to provide prior knowledge in the form of a list of potential transcription factors (TFs), which can then be integrated into traditional statistical causal discovery algorithms.

We use causal synthetic data generation as a downstream task to evaluate GRNs or priors suggested by LLMs in the absence of reliable ground-truth graphs. This method embeds gene regulatory networks into the scRNA-seq data generation process and has been shown to better preserve biological plausibility (Zinati et al., ).

We assess the practical utility of the inferred causal graph by comparing the synthetic data to an oracle dataset, evaluating both its statistical similarity and biological plausibility. This approach allows testing on real-world data instead of relying solely on causal benchmark datasets. Although causal benchmark datasets are paired with ground truth graphs, their widespread use in causal research may lead to their incorporation in the training data, potentially inflating the models’ perceived causal performance. Our results highlight the potential of general-purpose LLMs for GRN inference, as evidenced by both statistical and biological metrics. The best performance is achieved by combining LLMs with statistical GRN inference, pointing to a promising direction for scRNA-seq data analysis.

Our contributions:

- 1.
We show the effectiveness of using out of the box Large Language Models (LLMs) for inferring gene regulatory networks (GRNs), showcasing their ability to capture complex biological interactions.
- 2.
Comprehensive performance evaluation based on synthetic data generation. Our methodology effectively addresses the challenge posed by the absence of ground-truth causal graphs, facilitating a more rigorous assessment of inference methods.
- 3.
Extensive biological insight on the best-performing LLM for PBMC-All dataset.

## 2 Related Works

##### LLMs and Causality.

Gene regulatory network inference from scRNA-seq data traditionally relies on statistical causal discovery methods (Pratapa et al., ; Huynh-Thu et al., ; Moerman et al., ). However, causal discovery often requires external knowledge in the form of interventions (Brouillard et al., ), expert input (Kleinegesse et al., ), or priors from curated databases (Zinati et al., ). Recent advances in large language models (LLMs) offer a promising solution, as LLMs excel at integrating diverse knowledge and providing contextual information (Wan et al., ; Ma, ). Many of the recent works leverage LLMs for causal discovery by utilizing metadata, such as variable names, to infer causal relationships (Kıcıman et al., ; Ban et al., ; Vashishtha et al., ; Ban et al., ). Further optimizations, including more advanced prompting strategies beyond pairwise variable comparisons, have been developed to enhance causal discovery (Vashishtha et al., ; Jiralerspong et al., ). Subsequently,  (Sheth et al., ) explored the effectiveness of completing a partial causal graph across diverse domains. This work, however, considers LLM as an oracle to discover causal graphs for GRN inference.

##### LLMs and Biology.

Models based on transformer architecture and trained on DNA or RNA genomic sequences are effective at prediction and generation tasks (Zhang et al., ). However, the availability of generative Gene-LLMs is limited, they lack rich contextual information and are focused on specific tasks (Zhang et al., ). To overcome these limitations, foundation models are tailored for the genomic tasks (Wang et al., ; Cui et al., ). However, models trained on the genetic data do not achieve the wide context of LLMs like GPT4 that inject information from the scientific literature and freely available genomic databases. LLMs have been used for the tasks such as gene perturbation (Märtens et al., ), protein interaction prediction (Jin et al., ), gene selection (Toufiq et al., ), analyzing biomarkers (Elsborg & Salvatore, ) and cell annotation (Fang et al., ). To the best of our knowledge, general-purpose LLMs have not been used for gene regulatory network inference.

##### Causal Synthetic Data Generation.

Causally synthetic data generation is an approach that focuses on embedding true causal relationships within the generated data. Several techniques have been proposed in which the generator is guided by the causal acyclic graph (Cinquini et al., ; Van Breugel et al., ; Wen et al., ; Wang et al., ; Bandyopadhyay & Sarkar, ). We use a recently proposed causal GAN designed for scRNA-seq data generation, shown to produce more biologically plausible results (Zinati et al., ).

## 3 Preliminaries

### 3.1 Causal Graph

A directed acyclic graph (DAG) $\mathcal{G}=(\mathbf{V},\mathcal{E})$ is composed of a set of variables/vertices $\mathbf{V}$ and a set of (directed) edges $\mathcal{E}$ between them such that no cycle is formed.
Let $P$ be the probability distribution over the same set of variables $\mathbf{V}$.
$\mathcal{G}$ and $P$ satisfy the Markov condition if every variable is conditionally independent of its non-descendants given its parents.
Assuming the Markov condition, the joint distribution of variables $V_{1},V_{2},\ldots\in\mathbf{V}$ can be factorized as:

$$ $P[V_{1},V_{2},\ldots,V_{d}]=\prod_{i}P[V_{i}|Pa(V_{i})]\ .$ (1) $$

where $Pa(V_{i})$ denotes the set of parents of $V_{i}$.

A bipartite acyclic graph $\mathcal{G}=(\mathbf{U},\mathbf{V},\mathcal{E}$) is an instance of a DAG whose vertices are divided into two subsets $\mathbf{U}$ and $\mathbf{V}$. Each edge $\mathcal{E}_{i}$ connects only two vertices from $\mathbf{U}$ and $\mathbf{V}$.

#### 3.1.1 Gene Regulatory Network (GRN)

A gene regulatory network (GRN) is a collection of molecular regulators that interact with each other and with other substances in the cell to control the gene expression levels of mRNA and proteins. GRNs describe the relationships between genes, transcription factors (TFs), RNA molecules, and other regulatory elements within a biological system, illustrating how genes are turned on or off and how their expression is modulated over time and in different conditions. GRN can be expressed as a directed acyclic graph (DAG) 3.1. In this work we are considering a bipartite DAG consisting of TFs and target genes.

### 3.2 Generative Adversarial Network (GAN)

A Generative Adversarial Network (GAN) consists of two neural networks, the generator $G$ and the discriminator $D$, which are trained simultaneously in a zero-sum game. The generator $G$ maps a random noise vector $\mathbf{z}$ sampled from a prior distribution $p_{\mathbf{z}}(\mathbf{z})$ (typically Gaussian) to the data space, producing a synthetic sample $G(\mathbf{z})$. The discriminator $D$ maps an input $\mathbf{x}$ to a scalar value $D(\mathbf{x})\in[0,1]$, representing the probability that $\mathbf{x}$ is a real sample.

The training objective for GANs is formulated as a minimax game defined by the following function:

$$ $\min_{G}\max_{D}\,V(D,G)=\mathbb{E}_{\mathbf{x}\sim p_{\text{data}}(\mathbf{x} )}[\log D(\mathbf{x})]+\mathbb{E}_{\mathbf{z}\sim p_{\mathbf{z}}(\mathbf{z})}[ \log(1-D(G(\mathbf{z})))]$ $$

where $p_{\text{data}}(\mathbf{x})$ is the distribution of real data, and $p_{g}(\mathbf{x})$ is the distribution induced by $G$. The discriminator tries to maximize the probability of distinguishing real data from $G$’s samples, while $G$ aims to minimize $\log(1-D(G(\mathbf{z})))$ to generate realistic data. At equilibrium, the optimal generator replicates the true data distribution, making $p_{g}(\mathbf{x})=p_{\text{data}}(\mathbf{x})$.

#### 3.2.1 GRouNdGAN

Causal GAN (CGAN) proposed by Kocaoglu et al. () is a variant of GAN that incorporates a causal directed acyclic graph (DAG) into the data generation process. GRouNdGAN (Zinati et al., ) builds on CGAN, integrating components such as a causal controller, target generators, a critic, a labeler, and an anti-labeler. The training process for GRouNdGAN occurs in two stages. First, the causal controller, which generates transcription factor (TF) expression values, is pre-trained as the generator in a Wasserstein GAN with Gradient Penalty (WGAN-GP) framework. In the second stage, the pre-trained causal controller’s TF expression values, along with randomly generated noise, are provided as input to target generators. These target generators, each an instance of WGAN-GP, produce the expressions of target genes while respecting the TF-gene relationships dictated by the input gene regulatory network (GRN). Please refer to the original paper for a detailed explanation of GRouNdGAN’s architecture and training Zinati et al. ().

## 4 LLM4GRN

In this work, we propose
a novel approach integrating Large Language Models (LLMs) for GRN inference. The overall goal is to leverage the potential of LLMs in capturing complex biological interactions and to assess their utility. Given the lack of ground-truth data, we consider causal synthetic data generation as a downstream task to perform biological (causal) and statistical evaluations.
Our methodology involves two distinct experimental settings that leverage different knowledge about the potential transcription factors (TFs) to guide gene regulatory network (GRN)-informed data generation. In each setting, we introduce LLMs into the pipeline motivated by their ability to incorporate extensive contextual information.

Figure: (a) Setting 1
Refer to caption: x1.png

- •
In the first setting ( 1(a), Setting 1.A), we use LLM to infer the GRN graph by providing the LLM with a potential list of TF candidates sourced from a human-curated database.
- •
In the second setting ( 1(b), Setting 2.A), we utilize the LLM as the knowledge base, incorporating it earlier in the pipeline to infer potential TFs and to deduce the GRN graph.

We compare with Setting 1.B and 2.B where the Human knowledge base and LLM knowledge base is used by a statistical causal inference approach, GRNBoost2, respectively.

### 4.1 GRN Graph

To model the regulatory relationships between transcription factors and target genes. In this framework, we maintain a knowledge base (KB)
that contains comprehensive information regarding the regulatory interactions, specifically indicating which genes are target genes of specific transcription factors. KB takes as input a list of genes and outputs the corresponding target genes and their regulating transcription factors:

$$ $\displaystyle\text{KB}:\text{List of Genes}\rightarrow\{(\mathbf{T}_{i}, \mathbf{R}_{j})\mid\mathbf{T}_{i}\in\mathbf{T},\mathbf{R}_{j}\in\mathbf{R}\}$ $$

We represent the gene regulatory network as a bipartite graph $\mathcal{G}=(\mathbf{T},\mathbf{R},\mathcal{E})$, where $\mathbf{T}$ represents transcription factors (TFs), $\mathbf{R}$ contains target genes regulated by these TFs. The set $\mathcal{E}$ includes directed edges, where an edge $(\mathbf{T}_{i},\mathbf{R}_{j})\in\mathcal{E}$ indicates that transcription factor $\mathbf{T}_{i}$ regulates target gene $\mathbf{R}_{j}$.

### 4.2 Setting 1: Human Knowledge Base

Given the ability of large language models (LLMs) to generate causal graphs from metadata (such as variable names) (Abdulaal et al., ), we establish the foundational components of gene regulatory networks (GRNs) using a human knowledge base, denoted as $\text{KB}^{H}$. Let $\mathbf{T}^{H}$ represent the set of transcription factors identified through the human-curated database, and let $\mathbf{R}^{H}$ denote the set of target genes regulated by these factors. Total set of genes are defined by $G$. The directed edges in the graph are represented as $\mathcal{E}$, where an edge $(\mathbf{T}^{H}_{i},\mathbf{R}^{H}_{j})\in\mathcal{E}$ signifies that transcription factor $\mathbf{T}^{H}_{i}$ regulates target gene $\mathbf{R}^{H}_{j}$.

In Setting 1A, ( 1(a)), the LLM is employed to establish causal relationships between TFs and target genes, utilizing a list of TFs transcription factor candidates sourced from a human-curated database. We model the LLM as a function $\mathcal{F}_{\text{LLM}}$ that, given a set of metadata $\mathcal{M}$, produces a bipartite graph $\mathcal{G}=(\mathbf{T},\mathbf{R},\mathcal{E})$, expressed as:

$$ $\displaystyle\mathcal{F}_{\text{LLM}}(\mathcal{M},\mathbf{T}^{H},\mathbf{R}^{H })=\mathcal{G}$ (2) $$

The input metadata $\mathcal{M}$ encompasses gene names, transcription factors (TFs), single-cell RNA sequencing (scRNA-seq) data, and relevant biological context (such as species or experimental conditions) sourced from relevant literatures about the dataset.

It is important to highlight that, unlike traditional inference methods like GRNBoost2, the LLM-based GRN inference approach in this work does not rely on observational data, ensuring that individuals’ privacy in the dataset remains uncompromised. By integrating metadata from scRNA-seq datasets, the LLM can construct GRNs that are tailored specifically for the biological context of the data, potentially capturing nuances that statistical methods cannot. Detailed descriptions of the various prompting strategies employed can be found in the Appendix A.2.

### 4.3 Setting 2: LLM Knowledge Base

In Setting 2A, ( 1(b)), we utilize LLM knowledge base, denoted as $\text{KB}^{\text{LLM}}$, to establish partition between transcription factors and target genes. In this context, the LLM is tasked with extracting relevant information directly from its knowledge base, which includes extensive biological data and relationships derived from various sources.

$$ $\displaystyle\mathcal{F}_{\text{LLM}}(\mathcal{M})=\left(\mathbf{T}^{\text{LLM }},\mathbf{R}^{\text{LLM}}\right)$ (3) $$

Here, $\mathbf{T}^{\text{LLM}}$ denotes the set of transcription factors identified from the LLM knowledge base, and $\mathbf{R}^{\text{LLM}}$ represents the corresponding target genes. The input metadata $\mathcal{M}$ includes gene names, transcription factors (TFs), and relevant biological context, leveraging the extensive knowledge embedded in the LLM.

### 4.4 GRN-induced causal synthetic data generation

We use the GRNs obtained by either of the settings, to perform synthetic causal data generation. These GRNs are fed into a causal GAN algorithm, specifically GRouNdGAN (Zinati et al., ), which uses them to generate synthetic datasets that correspond to each GRN structure in a two-stage approach.

## 5 Results

We evaluate and compare the resulting synthetic datasets based on a range of statistical and biological metrics, allowing us to assess the quality and biological relevance of the data produced in each setting. Additionally, we analyze the TFs list generated by the LLM to gain insights, comparing them to priors from human curated database.

##### Experimental Setup.

We generated synthetic data using datasets and protocols from (Zinati et al., ). Our preprocessing focused on creating train, test, and validation sets while maintaining 1,000 genes across all datasets (see AppendixA.1 for details).
For each setting, we construct three different GRNs: one derived directly from the LLM, one generated by a causal discovery algorithm that incorporates the prior information (the dataset), and a third, randomly generated graph (based on TF list extracted from $\text{KB}^{\text{H}}$ or $\text{KB}^{\text{LLM}}$).
For the GRN inference of LLM, we prompted with contextual knowledge. We used the state-of-art pretrained GPT-4 model (Achiam et al., ). Given GPT-4 strong performance for causal discovery (Abdulaal et al., ) across different domains including genomics, we prompted (see Appendix A.2) in zero-shot fashion. We also compare open-source model, Llama-3.1-70B (Dubey et al., ).

##### Evaluation.

We employed four statistical metrics: Cosine and Euclidean distances to measure the differences between centroids, maximum mean discrepancy (MMD) to assess the proximity of high-dimensional distributions without centroid creation, and Random Forest Area Under the Receiver Operating Characteristic (RF-AUROC) to determine the distinguishability of real and synthetic cells. We evaluate the biological plausibility of the datasets by conducting single-cell RNA sequencing analysis using the Scanpy Python library. The cell annotations from the original dataset were utilized to annotate the synthetic datasets. Our analysis included log1p normalization and scaling to 10,000, focusing on 1,000 highly variable features. We performed dimensionality reduction using PCA with 50 principal components, followed by Uniform Manifold Approximation and Projection (UMAP). Cell type proportion analysis was conducted with custom code, while gene expression profiling for the top markers of each cell type was visualized using dot plots. For more details, refer to Appendices A and B.

##### Datasets.

For the direct comparison of the different settings, we focus on the dataset and genomic database information by (Zinati et al., ). Specifically, PBMC-All, PBMC-CTL and BoneMarrow data sets (details in Appendix A.1.1). Our objective is to assess the performance of large language models (LLMs) in gene regulatory network (GRN) inference. To evaluate the utility of the inferred GRNs, we employ causal synthetic data generation as a downstream task in scenarios where a reliable ground truth graph is unavailable.

### 5.1 Comparision against baseline

Evaluation of GPT-4 graphs.
For Setting 1, in the absence of ground truth for the GRN, we compare the overlap between the graph proposed by GRNBoost2, a statistical method, and the graph hypothesized by the GPT-4 . While GRNBoost2 is not the definitive ground truth, it serves as a useful reference point. Analyzing this overlap allows us to assess how the hypothesized connections align with established statistical methods and examine how these differences impact downstream synthetic data generation metrics. As a baseline, we also compared against a randomly generated causal graph. Additionally, we test the consistency of the LLM’s performance across different random seeds by measuring the overlap between graphs produced from multiple seeds. This helps us evaluate the robustness and stability of the GPT-4 -generated hypotheses.
We plot the overlaps for all of the datasets in  Figure 2. Our two main observations are that (1) the LLM-derived GRN demonstrates greater robustness compared to the GRNBoost2 (GRNB) GRN, particularly in terms of higher certainty; and (2) the overlap between the GPT-4 -generated GRN and the random GRN is smaller than between GPT-4 GRN and GRNBoost2 GRN, suggesting, that the GPT-4 could be generating meaningful graphs.

Figure: (a) PBMC-All
Refer to caption: x3.png

Evaluation of $\text{KB}^{\text{GPT4}}$.
In Setting 2, we introduced the LLM to filter TF and target genes from a list of genes. Similar to calculating overlaps in GRN evaluation, we also compute the overlap of TF produced by both approaches.
From Table 6 (Appendix C), we observe there exists around just about half an overlap between the two knowledge bases. Interestingly we observe less than $50\%$ overlap between $\text{KB}^{H}$ and $\text{KB}^{GPT4}$ for PBMC-CTL while a much higher overlap for the PBMC-ALL dataset. We also observed that LLM proposed a higher number of TFs, over 2 times in the case of the Bone Marrow dataset. As stated in the Table 3 (Appendix C), would change the density of the graphs potentially affecting the downstream tasks.

### 5.2 Statistical Evaluation of GRN Inference Methods on Synthetic Data

In Table 1, we present the statistical metrics results for the three datasets. Metrics are computed between a synthetic dataset of 1000 cells and a held-out test set of 1000 real cells for PBMC-ALL and PBMC-CTL. For BoneMarrow, 500 samples were used for synthetic and held-out test set. In GRouNdGAN’s imposed GRN, each gene is regulated by 10 transcription factors (TFs). Lower values ($\downarrow$) indicate better performance for all metrics, with the first two metrics representing the distance between the centroids of the real and synthetic cells. The “control” metrics are based on the real training dataset. The best performance values (excluding control and Stage1 which is a non-causal baseline) are highlighted in bold. Evaluations were carried out on 4 synthetic datasets, with experiments repeated using 2 cross-validation seeds.

##### Setting 1 Comparison.

The Setting 1 (using $\text{KB}^{H}$, Table 1) result shows the GPT-4-inferred GRN achieves the best performance across all metrics for PBMC-ALL. The LLM model shows the lowest Cosine distance of 0.00024 and Euclidean distance of 89, indicating that the synthetic data generated is most similar to the real data in terms of overall structure. Its lower MMD of 0.0072 compared to the GRNBoost2 model suggests it effectively captures subtle gene expression distributions. Additionally, the LLM model performs well on the Random Forest metric of 0.63, which measures binary classification accuracy in distinguishing real from synthetic data. Unsurprisingly, the Random Graph method performs the worst across all metrics, particularly with a high MMD of 0.0166 and the largest Euclidean distance of 121, highlighting the importance of informed GRN inference methods. In contrast, for the PBMC-CTL and BoneMarrow, the GRNBoost2 graph outperform the GPT-4 in this setting.

##### Setting 2 Comparison.

In Setting 2 (Table 1), where we incorporated the GPT-4 knowledge base ($\text{KB}^{GPT4}$ ), the GRNBoost2 method—combining GPT-4-proposed transcription factor (TF) lists with GRN inference—outperforms the fully GPT-4-based GRN approach (i.e., where the LLM proposes both the TF list and the connections between TFs and genes) across all datasets. GRNBoost2 also achieves the best overall performance for all datasets. For instance, in the PBMC-ALL dataset, it achieves a Euclidean distance of 83, improving from 121 in Setting 1, and an RF AUROC of 0.59, compared to 0.73 in Setting 1.
Overall, GRNBoost2 shows improvement over its performance in Setting 1, delivering better results across all metrics for both PBMC-ALL and PBMC-CTL. Notably, it surpasses the top results from Setting 1.
These findings suggest that integrating LLM-derived priors with GRNBoost2’s statistical inference yields a more accurate and robust representation of gene regulatory networks (GRNs). In contrast, BoneMarrow performed best in Setting 1 but remains competitive in Setting 2 when the $\text{KB}^{GPT4}$ is combined with GRNBoost2.

**Table 1: Performance of Different GRN Inference Methods in Simulating Realistic scRNA-seq Data across all 3 datasets. The best value for each dataset is presented in boldface and the best value in each setting is underline. The lower ($\downarrow$) the better for all metrics.**
|  |  | Cosine distance $\downarrow$ | Euclidean distance $\downarrow$ | MMD $\downarrow$ | RF AUROC $\downarrow$ |
| --- | --- | --- | --- | --- | --- |
| PBMC-ALL | Baseline |  |  |  |  |
| Control | 0.00029$\pm$0.00008 | 100$\pm$16 | 0.0051$\pm$ 0.001 | 0.49$\pm$0.017 |  |
| Stage 1 | 0.00036$\pm$0.00009 | 107$\pm$15 | 0.0057$\pm$0.005 | 0.55$\pm$0.021 |  |
| $\text{KB}^{\text{Random}}$ | 0.00132$\pm$0.00101 | 187$\pm$79 | 0.0214$\pm$0.009 | 0.86$\pm$0.080 |  |
| Setting 1 $\text{KB}^{\text{H}}$ |  |  |  |  |  |
| GPT-4 | 0.00024$\pm$0.00004 | 89$\pm$7 | 0.0072$\pm$0.001 | 0.63$\pm$0.043 |  |
| GRNBoost2 | 0.00047$\pm$0.00021 | 121$\pm$25 | 0.0139$\pm$0.006 | 0.73$\pm$0.050 |  |
| Random | 0.00045$\pm$0.00013 | 121$\pm$22 | 0.0166$\pm$0.003 | 0.85$\pm$0.019 |  |
| Setting 2 $\text{KB}^{\text{GPT4}}$ |  |  |  |  |  |
| GPT-4 | 0.00026$\pm$ 0.00009 | 90$\pm$13 | 0.0206$\pm$0.001 | 0.86$\pm$0.018 |  |
| GRNBoost2 | 0.00023$\pm$ 0.00008 | 83$\pm$17 | 0.0069$\pm$0.001 | 0.59$\pm$0.028 |  |
|  | Random | 0.00026$\pm$0.00011 | 92$\pm$14 | 0.0226$\pm$0.002 | 0.87$\pm$0.018 |
| PBMC-CTL | Baseline |  |  |  |  |
| Control | 0.00020$\pm$0.00005 | 57$\pm$7 | 0.0045$\pm$0.000 | 0.54$\pm$0.007 |  |
| Stage 1 | 0.00023$\pm$0.00003 | 63$\pm$4 | 0.0049$\pm$0.000 | 0.57$\pm$0.030 |  |
| $\text{KB}^{\text{Random}}$ | 0.00016$\pm$0.00002 | 51$\pm$3 | 0.0071$\pm$0.000 | 0.73$\pm$0.022 |  |
| Setting 1 $\text{KB}^{\text{H}}$ |  |  |  |  |  |
| GPT-4 | 0.00265$\pm$0.00253 | 176$\pm$115 | 0.0238$\pm$0.020 | 0.79$\pm$0.206 |  |
| GRNBoost2 | 0.00025$\pm$0.00004 | 65$\pm$6 | 0.0053$\pm$0.000 | 0.59$\pm$0.028 |  |
| Random | 0.00027$\pm$0.0001 | 66$\pm$7 | 0.0085$\pm$0.001 | 0.75$\pm$0.021 |  |
| Setting 2 $\text{KB}^{\text{GPT4}}$ |  |  |  |  |  |
| GPT-4 | 0.00028$\pm$0.00009 | 67$\pm$11 | 0.0067$\pm$0.001 | 0.70$\pm$0.025 |  |
| GRNBoost2 | 0.00020$\pm$0.00004 | 57$\pm$6 | 0.0049$\pm$0.000 | 0.59$\pm$0.019 |  |
|  | Random | 0.00022$\pm$0.00005 | 59$\pm$7 | 0.0080$\pm$0.000 | 0.76$\pm$0.023 |
| Bone Marrow | Baseline |  |  |  |  |
| Control | 0.00205$\pm$0.00018 | 80$\pm$4 | 0.0109$\pm$0.001 | 0.60$\pm$0.037 |  |
| Stage 1 | 0.00320$\pm$0.00115 | 101$\pm$16 | 0.0197$\pm$0.007 | 0.66$\pm$0.037 |  |
| $\text{KB}^{\text{Random}}$ | 0.00206$\pm$0.00034 | 80$\pm$7 | 0.0156$\pm$0.002 | 0.75$\pm$0.048 |  |
| Setting 1 $\text{KB}^{\text{H}}$ |  |  |  |  |  |
| GPT-4 | 0.00238$\pm$0.00052 | 86$\pm$9 | 0.0124$\pm$0.001 | 0.70$\pm$0.080 |  |
| GRNBoost2 | 0.00190$\pm$0.00023 | 77$\pm$5 | 0.0118$\pm$0.001 | 0.64$\pm$0.023 |  |
| Random | 0.00219$\pm$0.00030 | 82$\pm$6 | 0.0150$\pm$0.003 | 0.75$\pm$0.041 |  |
| Setting 2 $\text{KB}^{\text{GPT4}}$ |  |  |  |  |  |
| GPT-4 | 0.00217$\pm$0.00050 | 82$\pm$10 | 0.0137$\pm$0.003 | 0.72$\pm$0.044 |  |
| GRNBoost2 | 0.00193$\pm$0.00023 | 78$\pm$4 | 0.0119$\pm$0.001 | 0.64$\pm$0.033 |  |
|  | Random | 0.00297$\pm$0.00080 | 96$\pm$13 | 0.0172$\pm$0.004 | 0.79$\pm$0.046 |

##### Setting 2 Comparison using Llama-3.1-70B as Knowledge Base.

Table 2 presents the results of Setting 2 for Llama-3.1-70B Knowledge Base ($\text{KB}^{Llama}$) . The findings indicate that Llama, when queried for knowledge, achieves the best overall performance in combination with GRNBoost2, outperforming $\text{KB}^{GPT4}$ (Table 1) across all metrics. These results are promising, suggesting that integrating LLM knowledge bases with causal inference approaches like GRNBoost2 offers a highly promising direction for scRNA causal synthetic data generation.

**Table 2: Comparison of $\text{KB}^{\text{Llama}}$ - performance GRN Inference Methods in Simulating Realistic scRNA-seq Data for PBMC-ALL dataset. The best value is presented in boldface and the best value in each setting is underline. The lower ($\downarrow$) the better for all metrics.**
|  | Cosine distance $\downarrow$ | Euclidean distance $\downarrow$ | MMD $\downarrow$ | RF AUROC $\downarrow$ |
| --- | --- | --- | --- | --- |
| Baseline |  |  |  |  |
| Control | 0.00029$\pm$0.00008 | 100$\pm$16 | 0.0051$\pm$ 0.001 | 0.49$\pm$0.017 |
| Stage 1 | 0.00036$\pm$0.00009 | 107$\pm$15 | 0.0057$\pm$0.005 | 0.55$\pm$0.021 |
| $\text{KB}^{\text{Random}}$ | 0.00132$\pm$0.00101 | 187$\pm$79 | 0.0214$\pm$0.009 | 0.86$\pm$0.080 |
| Setting 1 $\text{KB}^{\text{H}}$ |  |  |  |  |
| GPT-4 | 0.00024$\pm$0.00004 | 89$\pm$7 | 0.0072$\pm$0.0012 | 0.63$\pm$0.04 |
| Llama-3.1 | 0.00114$\pm$0.00048 | 193$\pm$52 | 0.0244$\pm$0.0075 | 0.84$\pm$0.04 |
| GRNBoost2 | 0.00047$\pm$0.00021 | 121$\pm$25 | 0.0139$\pm$0.0058 | 0.73$\pm$0.05 |
| Random | 0.00045$\pm$0.00013 | 121$\pm$22 | 0.0166$\pm$0.0031 | 0.85$\pm$0.02 |
| Setting 2 $\text{KB}^{\text{GPT4}}$ |  |  |  |  |
| GPT-4 | 0.00026$\pm$0.00009 | 90$\pm$13 | 0.0206$\pm$0.0025 | 0.86$\pm$0.02 |
| GRNBoost2 | 0.00023$\pm$0.00008 | 83$\pm$17 | 0.0069$\pm$0.0012 | 0.59$\pm$0.03 |
| Random | 0.00026$\pm$0.00011 | 92$\pm$14 | 0.0226$\pm$0.0020 | 0.87$\pm$0.02 |
| Setting 2 $\text{KB}^{\text{Llama}}$ |  |  |  |  |
| Llama-3.1 | 0.00029$\pm$0.00005 | 97$\pm$10 | 0.0100$\pm$0.0007 | 0.77$\pm$0.030 |
| GRNBoost2 | 0.00022$\pm$0.00006 | 83$\pm$11 | 0.0067$\pm$0.0011 | 0.58$\pm$0.023 |
| Random | 0.00035$\pm$0.00009 | 105$\pm$17 | 0.0152$\pm$0.0011 | 0.86$\pm$0.028 |

### 5.3 Biological Plausibility of Causal Synthetic Data

We conduct cell type annotation based on the original dataset, gene expression analysis to identify top markers for each cell type and cell-type proportion analysis on the most statistically robust datasets. Full analysis can be found in Appendix D (subsection C.4).

Figure: (a) Setting 2. $\text{KB}^{\text{Llama}}$, GRNBoost2 Top Markers
Refer to caption: extracted/5937713/Images/bio_eval/llama_dotplot.png

#### 5.3.1 Gene Expression Profiling Per Cell Type

##### General Performance of KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2 in Cell-Specific Expression Profiles.

The $\text{KB}^{\text{Llama}}$ GRNBoost2 model shows that specific cell types, such as CD8+/CD45R+ naive cytotoxic T cells and CD4+ T helper cells, exhibit significantly higher mean expression levels for the same markers across multiple cells, while maintaining similar cell fraction of expression. Additionally, the $\text{KB}^{\text{Llama}}$ GRNBoost2 model demonstrates superior performance in elucidating expression patterns in certain cases, such as with dendritic cells. Although $\text{KB}^{\text{Llama}}$ GRNBoost2 synthetic data 4(b) introduces some noise, it generally improves cell-specific expression profiles, which could be further optimized for even cleaner cell type differentiation. A key observation was that when mean expression was higher in specific cell types, it tended to be lower in others, with fewer than $30\%$ of cells displaying similar expression fractions. In contrast, when expression was noisy across multiple cell types, both mean expression and cell fraction tended to be similar across these groups.

##### KB GPT4 superscript KB GPT4 \text{KB}^{\text{GPT4}} KB start_POSTSUPERSCRIPT GPT4 end_POSTSUPERSCRIPT GRNBoost2 Dataset Performance and Noise Patterns.

The $\text{KB}^{\text{GPT4}}$ GRNBoost2 model identifies top marker genes, revealing notable differences in expression profiles. While expression levels among some markers are noisy and indistinguishable, this model effectively highlights markers that are more discriminative for certain cell types. The $\text{KB}^{\text{GPT4}}$ GRNBoost2 dataset 4(b) exhibited more noise than the $\text{KB}^{\text{Llama}}$ dataset 4(b), with multiple markers expressed across various cell types at comparable mean expression levels and with higher cell fractions. For the $\text{KB}^{\text{H}}$ GRNBoost2 dataset ( 4(c)), noise was primarily observed in a few cell types, including CD4+/CD45A+/CD25– Naïve T cells and CD8+/CD45RA+ Naïve cytotoxic T cells. Some markers from other cell types also showed similar expression across multiple cell types. These noisy patterns were not observed in the Stage 1 or random datasets and were less pronounced in the $\text{KB}^{\text{Llama}}$ and $\text{KB}^{\text{GPT4}}$ models, though $\text{KB}^{\text{Llama}}$ GRNBoost2 outperformed $\text{KB}^{\text{GPT4}}$ GRNBoost2 in reducing noise for these specific cell types.

#### 5.3.2 Differential Cell Type Proportion Analysis

Figure: (a) Setting 2, $\text{KB}^{\text{GPT4}}$, GRNBoost2 graph.
Refer to caption: extracted/5937713/Images/bio_eval/kbgpt_prop.png

##### Variations in Cell Type Proportions Across Datasets.

A notable observation is that the cell type proportions in each generated dataset differ significantly from those in the original dataset, indicating a noisy overall expression profile that can alter the distribution of cell types. In the $\text{KB}^{\text{Llama}}$ GRNBoost2 dataset ( 5(b)), CD8+/CD45RA+ Naïve Cytotoxic T cells were the most abundant at $34.1\%$, followed closely by CD56+ NK cells, CD8+ Cytotoxic T cells, and CD4+/CD25+ T regulatory cells. This trend was similarly observed $\text{KB}^{\text{GPT4}}$ GRNBoost2 ( 5(a)) and other datasets, albeit with slightly varying percentages.

##### Overall Performance of KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2.

We observe that CD4+/CD45RA+/CD25– Naïve T Cells and CD8+/CD45RA+ Naïve Cytotoxic T Cells exhibited noisy expression patterns consistently across all datasets. The original dataset performed better for CD8+/CD45RA+ Cytotoxic T Cells, showing lower noise in their expression compared to the $\text{KB}^{\text{Llama}}$ GRNBoost2 ( 5(b)) dataset. Among the various models, the $\text{KB}^{\text{Llama}}$ GRNBoost2 model emerged as the best performer, providing clearer segregation of cell types and reduced noise, particularly for CD4+/CD45RA+/CD25– Naïve T Cells and dendritic cells. In contrast, the $\text{KB}^{\text{GPT4}}$ GRNBoost2 ( 5(a)) model exhibited more noise than $\text{KB}^{\text{Llama}}$ GRNBoost2, with similar markers expressed across various cell types at comparable mean expression levels.

### 5.4 Discussion and Limitations

We observe promising results in using LLM for GRN discovery, especially on the PBMC dataset. The hybrid approach incorporating TFs suggested by LLM ($\text{KB}^{\text{LLM}}$) and GRNBoost2 causal discovery yield the best overall performance (for PBMC-ALL and PBMC-CTL data sets). Surprisingly, a smaller open-source model Llama has shown the best result in this setting. One possible explanation could be the number of TFs proposed in Setting 2. The number of TFs in $\text{KB}^{\text{Llama}}$ is higher than those in $\text{KB}^{\text{H}}$ or $\text{KB}^{\text{GPT4}}$. Our additional experiments (subsection C.2) indicate that the number of transcription factors (TFs) may be significant, as it provides a broader set of variables for the GRNBoost2 algorithm to evaluate during causal discovery.
Nevertheless, although further evaluation is needed, new Llama models might perform as good or better than state-of-the-art GPT for specific tasks Valero-Lara et al. (). However, Llama performs worse on the more challenging, GRN construction task.
The lower performance of LLMs on CTL and BoneMarrow data sets is not surprising. Recent studies suggest that genomic LLMs under-perform on cell-specific data (Tang & Koo, ). In addition, human transcription factors outnumber mice transcriptomic information in the databases (Members & Partners, ) making the information on mice genes scarcer in the LLMs training data.

Biological plausibility analysis reveals that the Setting 2, $\text{KB}^{\text{Llama}}$ dataset significantly improves cell type differentiation and reduces noise compared to other generated datasets, particularly for CD4+/CD45RA+/CD25– Naïve T cells and dendritic cells, which indicates better performance of the model due to its ability to handle large contexts (Xiong et al., ).

In contrast, while the Setting 2, $\text{KB}^{\text{GPT4}}$ model exhibited higher noise and less specificity, the Setting 1, $\text{KB}^{\text{H}}$, GRNBoost2 model also showed noisy patterns, especially in Naïve T and cytotoxic T cells, emphasizing the need for models that clearly delineate cell type-specific expression. Notably, discrepancies in cell type proportions across generated datasets, particularly the inflated presence of CD8+/CD45RA+ Naïve Cytotoxic T cells, raise concerns about the biological relevance of these models. As reported previously, large-scale single-cell studies tend to be more noisy which can lead to sub-optimal biological inferences (Kavran & Clauset, ). Therefore, further refinements are essential across all approaches to enhance specificity and reliability in representing true cellular compositions, which is crucial for accurate downstream biological analyses and interpretations. Finally, the $\text{KB}^{\text{H}}$ GPT-4 model exhibited highly noisy or non-specific expression profiles, particularly for CD4+/CD45A+/CD25– Naïve T cells, CD8+/CD45RA+ Naïve cytotoxic T cells, CD4+ T helper cells, CD4+/CD25+ T regulatory cells, and CD4+/CD45RO+ memory cells. For other cell types, the model either failed to express most top markers or showed high cell fractions with reduced mean expression across different cell types.

One of the limitations of the study is the unifying constraints that are imposed on the GRN discovery due to the parametric requirements of GRouNdGan (Zinati et al., ). Namely, all GRN graphs are set to be bipartite graphs with the same number of TFs and same number of target genes for each TF. These constraints restrict the diversity of the generated graphs, resulting in fairly similar performance metrics among different GRN discovery approaches. In addition, multi-layer graphs with transcription factors, cofactors, and target genes are more realistic and can better reflect biological complexity (Karlebach & Shamir, ).

Despite the promising results, we note, that LLMs should be used cautiously in any high-stakes decision-making. The models are prone to reporting false information with confidence (Ji et al., ; Farquhar et al., ). In addition, they are susceptible to biases in the training data.
One of the most common sources of bias in machine learning is the under-representation of minority populations in the training data. Transcription factor databases and literature often lack balanced representation across ethnicity, ancestry, gender, and age groups. The vast majority of genetic data and studies are based on individuals of European ancestry (Sirugo et al., ; Bentley et al., ).
The underrepresentation of non-European ancestry populations in genomic databases can obscure gene-disease associations that are rare in European groups, leading to treatments that may not be effective for these communities (Landry et al., ; Tawfik et al., ; Barral-Arca et al., ). Age groups, specifically children and elderly, females, low-income populations and individuals from remote rural areas are also underrepresented in clinical trials and data, making them less likely to benefit from the achievements of precision medicine (Mosenifar, ; Davis et al., ; Steinberg et al., ). LLMs are shown to exacerbate some existing health disparities (Pfohl et al., ), and are likely to be influenced by the under-representation bias in the genomic data.

## 6 Conclusion

In conclusion, our analysis highlights the potential of using Large Language Models (LLMs) for gene regulatory network (GRN) inference, as supported by both statistical and biological performance metrics. Notably, the best performing approach combines LLM-derived TF priors with GRNBoost2 for statistical inference. In addition, we see good potential for using open-source models such as Llama for GRNs, and we aim to explore their utility further, either directly or with fine-tuning. Our biological results suggest that additional refinement is needed. For this we foresee developing GRNs tailored to individual cell types, as evidence indicates that GRNs are often cell-type-specific.

## Acknowledgements

This work is funded in part by Bundesministeriums fur Bildung
und Forschung (PriSyn), grant No. 16KISAO29K. The work is also supported by Medizininformatik-Plattform "Privatsphären-schutzende
Analytik in der Medizin" (PrivateAIM), grant No. 01ZZ2316G, and ELSA – European
Lighthouse on Secure and Safe AI funded by the European Union under grant agreement No.
101070617. Moreover, the computation resources used in this work are supported by the Helmholtz
Association’s Initiative and Networking Fund on the HAICORE@FZJ partition. Views and opinions
expressed are, however, those of the authors only and do not necessarily reflect those of the European
Union or European Commission. Neither the European Union nor the European Commission can be
held responsible for them.

## Appendix A Experimental setup

### A.1 Data Preparation and Modifications

In synthetic data generation, we use the same data sets and follow the protocol described by Zinati et al. Zinati et al. ().

We downloaded the three datasets using the link provided at [https://emad-combine-lab.github.io/GRouNdGAN/tutorial.html#demo-datasets](https://emad-combine-lab.github.io/GRouNdGAN/tutorial.html#demo-datasets).

To ensure consistent results across different runs, we modified the code in two key ways: (1) we seeded the randomness in the ‘main.py’ file. While the GRNBoost2 algorithm is already seeded by default, ensuring consistent output, the rest of the code was not. Given that we are testing different GRN graphs, it was critical to control for randomness to prevent it from influencing our results. (2) We addressed additional bugs to ensure smooth execution, such as handling sparse data loading and ensuring that the code was properly assigned to the correct device for training. All modifications have been incorporated into our forked version of the repository, available upon publication. Some of the bugs we identified were also reported and fixed by the original authors.

After downloading the raw data, we preprocessed it to create separate train, test, and validation files. Following the paper’s recommendations, we used a test set of 1,000 samples for the PBMC and CTL datasets, and 500 samples for the BoneMarrow dataset. For validation, we set aside 1,000 samples each for PBMC and CTL, and 500 for BoneMarrow. The subsets were disjoint. We refer to these preprocessed datasets as “Real Data” throughout the paper. The preprocessed data will be made available upon publication.

During preprocessing, we observed that the BoneMarrow dataset had two genes absent in the training data, as they were not expressed in any cells. This issue did not occur in the full dataset, where genes expressed in fewer than three cells and cells expressing fewer than 10 genes had been filtered out. However, after splitting the data into train, test, and validation sets, some genes in specific splits were not expressed at all.

To address this, we modified the code to filter out genes expressed in fewer than one cell for each data split. This reduced the number of genes in the BoneMarrow dataset from 1,000 to 910. In the PBMC dataset, we found that 493 genes were not expressed in the test set, requiring extensive filtering.

To resolve this, we revised our approach. Instead of applying the original threshold, which discarded genes expressed in fewer than three cells, we increased this threshold to 680 cells for the PBMC dataset, equivalent to 0.01% of the full dataset. This change allowed us to retain the full 1,000 highly variable genes. However, the resulting set of 1,000 genes may differ slightly from those identified with the initial criteria, given that the full PBMC dataset contains over 32,738 genes, while BoneMarrow contains 3,451 genes.

This revised strategy allowed us to retain the full 1,000 genes for the PBMC and CTL datasets but not for BoneMarrow. To address this, we further adjusted the strategy by filtering out genes expressed in fewer than the number of cells in the test set for each dataset (1,000 for PBMC and CTL, and 500 for BoneMarrow). This approach successfully retained 1,000 genes across all datasets.

After preprocessing, we followed the rest of the process as detailed in the tutorial provided by the authors.

#### A.1.1 Data Sets

- •
Human peripheral blood mononuclear cell (PBMC-Al). 68579 samples corresponding to 11 cell types
- •
Human peripheral blood mononuclear CD8+ Cytotoxic T-cells (PBMC-CTL). 20773 samples from the most common cell type in PBMC-All
- •
Mouse bone marrow Hematopoietic stem cells lineage differentiation (BoneMarrow). 2730 cells.

**Table 3: The density of the causal graphs.**
| Dataset | # TFs | # Targets | # Genes | # Possible Edges | # Imposed Edges | GRN density Edges |
| --- | --- | --- | --- | --- | --- | --- |
| PBMC | 75 | 925 | 1000 | 69375 | 9250 | 0.1333 |
| CTL | 65 | 935 | 1000 | 60775 | 9350 | 0.1538 |
| BoneMarrow | 68 | 932 | 1000 | 63376 | 9320 | 0.1471 |

### A.2 LLM Prompting Strategies

In this study, we utilize prompting techniques to guide the behavior of a Large Language Model (LLM) for gene regulatory network (GRN) inference and other downstream analyses. The model employed is based on a pretrained large language model, specifically GPT-4, which has not been fine-tuned for this task. Instead, we rely on advanced prompting strategies, including Chain of Thought (CoT) reasoning and context provision, to enhance the performance of the LLM in generating biologically plausible results.

#### A.2.1 Chain of Thought (CoT) Prompting

Chain of Thought (CoT) prompting is a method that encourages the LLM to reason through intermediate steps, producing a more transparent and logical progression towards its final answer. By guiding the model to provide step-by-step explanations before arriving at a conclusion, we aim to improve the interpretability and accuracy of its predictions. CoT prompting is particularly valuable for tasks that require complex reasoning, such as GRN inference, where multiple factors, such as transcription factor interactions and gene expression patterns, must be considered.

#### A.2.2 Context Provision

To further enhance the performance of the LLM, we utilize context provision by supplying the model with relevant background information prior to prompting. This technique is especially important when the task requires domain-specific knowledge, such as GRN inference from single-cell RNA sequencing (scRNA-seq) data. By embedding relevant context in the prompt, we can better align the model’s responses with the biological characteristics of the data being analyzed. Specifically, we provide excerpts from the original articles where the analyzed data sets were introduced Paul et al. (); Zheng et al. ().

#### A.2.3 Extracting Potential TFs From The Gene List

we provide the LLM with a curated list of genes, accompanied by relevant contextual information such as biological function, tissue type, and experimental conditions derived from scRNA-seq datasets. The prompt explicitly requests the LLM to identify and propose TFs that are known to regulate the expression of each gene within the list.
In our approach, we start with a total of 1000 genes and sequentially query the Large Language Model (LLM) about subsets of 20 genes at a time to extract potential transcription factors (TFs). For each initial query, the LLM identifies and proposes TFs that may regulate the selected genes, leveraging its extensive contextual knowledge of gene regulatory mechanisms. In subsequent prompts, we maintain a 50$\%$ overlap with the previous set of genes, ensuring that 10 of the genes are revisited while introducing 10 new genes. This iterative process allows us to refine the extraction of TFs, incorporating insights from the previous queries while continuously expanding our understanding of the regulatory landscape. By doing so, we aim to capture a comprehensive set of potential TFs that interact with the entire gene pool, facilitating a more robust inference of the Gene Regulatory Network (GRN).

## Appendix B Metrics

### B.1 Statistical Evaluation Metrics

##### Euclidean Distance:

To compute the Euclidean distance between the centroids of the real and simulated cells, we first calculate the centroid by finding the mean along the gene axis across all simulated and real cells. The Euclidean distance
$d(r,s)$ is then given by:

$$ $d(r,s)=\|\mu(R)-\mu(S)\|_{2}$ $$

Where $R$ amd $S$ are matrix of real and simulated cells, with elements $R_{i,j}$ and $S_{i,j}$ where $i$ indexes the cells and $j$ indexes the genes. $\mu(R)$ and $\mu{S}$ is the mean vector of the real and synthetic cells along the gene axis (columns) respectively.

$$ $\mu(R)=\left(\frac{1}{n}\sum_{i=1}^{n}R_{i,j}\right)_{j=1,\ldots,m}~{}\text{ and}~{}\mu(S)=\left(\frac{1}{n}\sum_{i=1}^{n}S_{i,j}\right)_{j=1,\ldots,m}$ $$

##### Cosine Distance:

This computes the cosine distance between the centroids of the real and simulated cells. The centroid is obtained by calculating the mean along the gene axis across all simulated and real cells.

$$ $d_{\cos}(r,s)=1-\frac{\mu(R)\cdot\mu(S)}{\|\mu(R)\|_{2}\|\mu(S)\|_{2}}$ $$

Cosine distance measures the difference in orientation of the two centroids while the Euclidean distance measures the absolute difference between the two centroids. Since cosine focus on angle between the vectors, it is useful in comparing the shape of data distributions rather than their scale.

##### Maximum Mean Discrepancy (MMD):

Maximum Mean Discrepancy (MMD) serves as a non-parametric two-sample test to determine if samples are drawn from the same distribution. MMD metric is identified as a particularly convenient method for assessing the similarity of real data. We followed the description in Zinat et al. Zinati et al. ().

##### Discriminative Metric (RF AUROC):

This metric evaluates a model’s ability to distinguish between real and synthetic datasets. A random forest (RF) classifier is employed, and the area under the receiver operating characteristic (AUROC) curve is used to assess whether real and simulated cells can be effectively differentiated.

## Appendix C Additional Results

### C.1 Statistical Metrics

##### Comparison of synthetically generated dataset with known causal graph approach.

We include results for a synthetically generated dataset, LinearUniform (Table 4), to evaluate the impact of causality. Since the LinearUniform dataset was generated with a known causal graph, we applied GRNBoost2 as the causal inference method. For comparison, the non-causal model is represented in Stage 1. The results indicate that the causal model (GRNBoost2) outperforms across all three metrics.

**Table 4: Evaluation of synthetically generated dataset.**
|  | Cosine distance | Euclidean distance | MMD | RF AUROC |
| --- | --- | --- | --- | --- |
| LinearUniform |  |  |  |  |
| Control | 0.00082$\pm$0.00012 | 108$\pm$8 | 0.0090$\pm$0.000 | 0.57$\pm$0.016 |
| Stage 1 | 0.00097$\pm$0.00010 | 119$\pm$6 | 0.0139$\pm$0.001 | 0.64$\pm$0.020 |
| GRNBoost2 | 0.00061$\pm$0.00011 | 93$\pm$9 | 0.0152$\pm$0.001 | 0.63$\pm$0.028 |

### C.2 Llama-3.1 Ablation

We carried out ablation of using Llama-3.1 as knowledge base to study the impact of size of TFs used for GRN creation.The result shown in Table 2 indicate that the number of transcription factors (TFs) may be significant, as it provides a broader set of variables for the GRNBoost2 algorithm to evaluate during causal discovery.

**Table 5: Comparison of $\text{KB}^{\text{Llama}}$ with varying number of TFs - performance of GRNBoost2 Inference Method in Simulating Realistic scRNA-seq Data for PBMC-ALL dataset. The best value is presented in boldface. The lower ($\downarrow$) the better for all metrics.**
|  | Cosine distance $\downarrow$ | Euclidean distance $\downarrow$ | MMD $\downarrow$ | RF AUROC $\downarrow$ |
| --- | --- | --- | --- | --- |
| Setting 2 $\text{KB}^{\text{Llama}}$ |  |  |  |  |
| GRNBoost2 (Tf=266) | 0.00022$\pm$0.00006 | 83$\pm$11 | 0.0067$\pm$0.0011 | 0.58$\pm$0.023 |
| GRNBoost2 (Tf=95) | 0.00032$\pm$0.00008 | 103$\pm$17 | 0.0083$\pm$0.0006 | 0.65$\pm$0.024 |
| GRNBoost2 (Tf=75) | 0.00034$\pm$0.00023 | 102$\pm$38 | 0.0097$\pm$0.0029 | 0.67$\pm$0.043 |

#### C.2.1 Overlaps between Human KB and GPT-4 KB

**Table 6: The overlaps of TFs between of $\text{KB}^{\text{H}}$ and $\text{KB}^{\text{GPT4}}$.**
| Dataset | # TFs | % Overlap |
| --- | --- | --- |
| PBMC | 95 | 64.00 |
| CTL | 135 | 49.23 |
| BM | 157 | 55.82 |

#### C.2.2 Density of graphs

**Table 7: The density of different causal graphs wrt different KB.**
| Dataset | # TFs | # Targets | # Genes | # Possible Edges | # Imposed Edges | GRN density Edges |
| --- | --- | --- | --- | --- | --- | --- |
| KB^H |  |  |  |  |  |  |
| PBMC | 75 | 925 | 1000 | 69375 | 9250 | 0.1333 |
| CTL | 65 | 935 | 1000 | 60775 | 9350 | 0.1538 |
| BoneMarrow | 68 | 932 | 1000 | 63376 | 9320 | 0.1471 |
| KB^GPT |  |  |  |  |  |  |
| PBMC | 95 | 905 | 1000 | 85975 | 9050 | 0.1052 |
| CTL | 135 | 865 | 1000 | 116775 | 8650 | 0.0740 |
| BoneMarrow | 157 | 843 | 1000 | 132351 | 8430 | 0.0639 |
| KB^Llama |  |  |  |  |  |  |
| PBMC | 266 | 843 | 1000 | 224238 | 8430 | 0.0375 |

#### C.2.3 Qualitative Evaluation of LLM Priors

The TFs that were misssing from LLM list as compared to the human curated TFs list for the PBMC data set:"PLEK", "SNAI3", "ZNF264", "ZNF708", "ZNF277","SMARCC1", "ZBTB38","BAZ2A","ARID5A","TGIF1","LYAR","DNAJC1","HOPX".
These transcription factors (TFs) play diverse roles in regulating immune function, cell differentiation, and responses to stress or inflammation. Their expression in PBMCs highlights their importance in maintaining immune homeostasis, responding to infections, and their potential involvement in disease states such as cancer, autoimmune diseases, and inflammatory disorders.

Extra genes proposed as TFs by LLM: "TNFRSF18","IL2RB","NCOA4","CHD8","CDK9","LMO4", "HEXIM1","TAF10", "CCNK","POLR2A","BRD4","ANKRD11".
Notably the extra transcription factors (TFs) proposed by LLM play key roles in regulating immune cell function, influencing various aspects of immune responses, disease states (like cancer and autoimmunity), and cell differentiation pathways.

The list with the missing TFs has more genes that are not as well-defined in the literature. Specifically, genes like SNAI3, ZNF264, ZNF708, and ZNF277 have limited characterization and fewer studies addressing their specific roles in cellular processes, especially in the context of immune cells and PBMCs.
In contrast, many genes in the extra proposed TFs list (e.g., TNFRSF18, CDK9, BRD4) are more extensively studied and have clearer functions, particularly in transcriptional regulation and immune processes.

### C.3 Qualitative Evaluation of Synthetic Data Sets

In addition we perform a qualitative evaluation of the synthetic data sets by visualizing them using TSNE projections (Figure 6). We observe that using random GRN induces "hallucinated" extra clusters in the data, whereas both GRNB2 and LLM proposed graphs stay faithful to the original distribution.

Figure: Figure 6: TSNE projections of synthetic vs. real data for different GRN graphs (Setting 1). “Control” corresponds to the projection of training and testing data, “GRNB2” -the synthetic data based on GRNBoost2 graph, “LLM” - synthetic data based on LLM graph, and “RAND” - the data based on the random graph. Red dots correspond to the real data points and blue ones - to the synthetic data points. “Hallucinated” extra blue clusters in the RAND graph are marked with red circles.
Refer to caption: extracted/5937713/Images/TSNESynthetix.png

### C.4 Biological evaluation

##### KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2 Dataset and Improved Cell Type Segregation

Compared to the original dataset, it was observed that the $\text{KB}^{\text{Llama}}$ GRNBoost2 dataset segregates cell types more effectively. In the original dataset, four of the top five marker genes for CD4+/CD45RA+/CD25– Naïve T cells exhibited similar expression across multiple cell types, introducing noise that could complicate biological analysis. In contrast, the $\text{KB}^{\text{Llama}}$ dataset demonstrated more distinct expression for CD4+/CD45RA+/CD25– Naïve T cells, with the exception of a single marker, LTB, which was expressed in other cell types as well.

##### Reduction of Noise in Other Cell Types in KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2Dataset

Additionally, dendritic cells in the $\text{KB}^{\text{Llama}}$ GRNBoost2 model exhibited a less noisy expression pattern compared to the original dataset, a trend also observed in CD8+ cytotoxic T cells and CD4+/CD25+ T regulatory cells. A notable difference was found in CD8+/CD45RA+ cytotoxic cells, where the original dataset showed little to no expression of top markers, while the $\text{KB}^{\text{Llama}}$ dataset had a more widespread, noisy expression.

##### General Performance of KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2 in Cell-Specific Expression Profiles

Although $\text{KB}^{\text{Llama}}$ GRNBoost2 synthetic data 4(b) introduces some noise, it generally improves cell-specific expression profiles, which could be further optimized for even cleaner cell type differentiation. A key observation was that when mean expression was higher in specific cell types, it tended to be lower in others, with fewer than $30\%$ of cells displaying similar expression fractions. In contrast, when expression was noisy across multiple cell types, both mean expression and cell fraction tended to be similar across these groups.

##### Comparison with KB GPT4 superscript KB GPT4 \text{KB}^{\text{GPT4}} KB start_POSTSUPERSCRIPT GPT4 end_POSTSUPERSCRIPT B GRNBoost2 and Random Datasets

In comparison, the $\text{KB}^{\text{GPT4}}$ GRNBoost2 dataset 4(b) exhibited more noise than the $\text{KB}^{\text{Llama}}$ dataset 4(b), with multiple markers expressed across various cell types at comparable mean expression levels and with higher cell fractions. The random GRN dataset resembled the Stage 1 non-causal dataset, with mean expression profiles being more clearly defined for each specific cell type. However, cell fractions for the same markers in other cell types were still significant, often exceeding $80\%$, indicating suboptimal gene expression specificity.

##### Stage 1 Dataset and Noise Reduction in Markers

Stage 1 data demonstrated a less noisy expression pattern, where markers were primarily confined to their specific cell types, showing lower mean expression in others. However, even though mean expression in non-specific cell types were low, cell fractions remained similar to those of the main cell type, which is not ideal for clear differentiation.

##### KB H superscript KB H \text{KB}^{\text{H}} KB start_POSTSUPERSCRIPT H end_POSTSUPERSCRIPT GRNBoost2 Dataset Performance and Noise Patterns

For the $\text{KB}^{\text{H}}$ GRNBoost2d dataset, noise was primarily observed in a few cell types, including CD4+/CD45A+/CD25– Naïve T cells and CD8+/CD45RA+ Naïve cytotoxic T cells. Some markers from other cell types also showed similar expression across multiple cell types. These noisy patterns were not observed in the Stage 1 or random datasets and were less pronounced in the $\text{KB}^{\text{Llama}}$ and $\text{KB}^{\text{GPT4}}$ models, though $\text{KB}^{\text{Llama}}$ GRNBoost2 outperformed $\text{KB}^{\text{GPT}}$ in reducing noise for these specific cell types. Beyond these examples, most cell types in the $\text{KB}^{\text{H}}$ GRNBoost2 dataset did not exhibit significant noise.

##### KB H superscript KB H \text{KB}^{\text{H}} KB start_POSTSUPERSCRIPT H end_POSTSUPERSCRIPT LLM Model’s Noisy Expression Profiles

Finally, the $\text{KB}^{\text{H}}$ LLM model exhibited highly noisy or non-specific expression profiles, particularly for CD4+/CD45A+/CD25– Naïve T cells, CD8+/CD45RA+ Naïve cytotoxic T cells, CD4+ T helper cells, CD4+/CD25+ T regulatory cells, and CD4+/CD45RO+ memory cells. For other cell types, the model either failed to express most top markers or showed high cell fractions with reduced mean expression across different cell types.

##### Variations in Cell Type Proportions Across Datasets

A notable observation is that the cell type proportions in each generated dataset differ significantly from those in the original dataset, indicating a noisy overall expression profile that can alter the distribution of cell types. In the $\text{KB}^{\text{Llama}}$ GRNBoost2 dataset, CD8+/CD45RA+ Naïve Cytotoxic T cells were the most abundant at $34.1\%$, followed closely by CD56+ NK cells, CD8+ Cytotoxic T cells, and CD4+/CD25+ T regulatory cells. This trend was similarly observed in the Stage 1, $\text{KB}^{\text{H}}$GRNBoost2, $\text{KB}^{\text{H}}$ LLM, and $\text{KB}^{\text{GPT4}}$ GRNBoost2 datasets, albeit with slightly varying percentages.

##### Discrepancies in Cell Type Proportions

The most considerable variation in proportions was noted for CD8+/CD45RA+ Naïve Cytotoxic T cells, with $\text{KB}^{\text{Llama}}$ GRNBoost2 reporting the highest percentage at $34.1\%$ and GPT-4 the lowest at $28.2\%$. Importantly, all generated datasets displayed significant differences when compared to the original dataset, where CD8+/CD45RA+ Naïve Cytotoxic T cells constituted only $0.1\%$ of the population, and CD4+/CD45RA+ Naïve T cells accounted for $65.1\%$. These findings suggest that the models manipulate the expression profiles in such a way that the proportions of CD8+/CD45RA+ Naïve Cytotoxic T cells, CD4+/CD45RA+ Naïve T cells, and other cell types are elevated in the generated datasets.

##### Overall Performance and Future Improvements for KB Llama superscript KB Llama \text{KB}^{\text{Llama}} KB start_POSTSUPERSCRIPT Llama end_POSTSUPERSCRIPT GRNBoost2

Therefore, it was observed that CD4+/CD45RA+/CD25– Naïve T Cells and CD8+/CD45RA+ Naïve Cytotoxic T Cells exhibited noisy expression patterns consistently across all datasets. The original dataset performed better for CD8+/CD45RA+ Cytotoxic T Cells, showing lower noise in their expression compared to the $\text{KB}^{\text{Llama}}$ GRNBoost2 generated dataset. Among the various models, the Llama model emerged as the best performer, providing clearer segregation of cell types and reduced noise, particularly for CD4+/CD45RA+/CD25– Naïve T Cells and dendritic cells.

##### Conclusion: Challenges and Future Refinements

In contrast, the $\text{KB}^{\text{GPT4}}$ GRNBoost2 model exhibited more noise than $\text{KB}^{\text{Llama}}$ GRNBoost2, with similar markers expressed across various cell types at comparable mean expression levels. The Stage 1 model displayed less noisy expression patterns, but some markers were still expressed significantly across non-specific cell types, while the random-generated dataset had poorer cell-specific profiling. In conclusion, the discrepancies in cell-type proportions between generated datasets and the original dataset highlight the challenges in achieving accurate cell-type specificity. These variations can significantly impact downstream analyses and interpretations in biological research. Consequently, refining these models to enhance the specificity of marker expression is essential for producing more reliable datasets that reflect the true cellular composition observed in biological samples.

Figure: (a) Setting 1, $\text{KB}^{\text{H}}$, LLM graph.
Refer to caption: extracted/5937713/Images/bio_eval/gpt_dotplot.png

Figure: (a) Original data.
Refer to caption: extracted/5937713/Images/bio_eval/originaldata_dotplot.png

Figure: (a) Setting 1, $\text{KB}^{\text{H}}$, LLM graph
Refer to caption: extracted/5937713/Images/bio_eval/gpt_prop.png

Figure: (a) Original data.
Refer to caption: extracted/5937713/Images/bio_eval/originalProp.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 1
questions: - Please consider adding other orthogonal direct ways of validating GRNs, such as the one suggested above by testing the exact same probing strategy on whether it can identify context specific variations of a well known and validated TF GRN
- Please include an orthogonal method to validate the context-specific TF list, such as distinguishing TFs based on uncorrelated expression patterns across distinct cell types, as mentioned in the Weaknesses section
soundness: 2
strengths: - The study addresses a critical need for context-specific GRNs, as many practitioners rely on generic knowledge-base GRNs, which lack context specificity and often lead to less accurate inferences.
- A thorough analysis, evaluating both the TF lists generated by foundation models and the resulting GRNs with controls at each step
summary: This paper explores the application of general-purpose foundation models (GPT-4, LLaMA) for identifying context-specific transcription factors (TFs) and constructing corresponding gene regulatory networks (GRNs). The authors evaluate the GRNs by contrasting them with knowledge base-derived GRNs and those inferred using GRNBoost2, then apply causal synthetic data generation based on the method of Zinati et al. (2024) to assess the biological plausibility of the generated networks.
weaknesses: The primary limitation is the lack of objective metrics for validating the inferred GRNs. While the authors use synthetic data generation as an indirect evaluation method, results vary significantly across datasets, with the GPT-based GRNs performing well on PBMC_ALL but showing weaker, near-random performance on other datasets (as measured by MMD and RF AUROC). These results indicate low systematic performance on the author's suggested standards.

The above may also be a result of knock-effects caused by the Causal GAN synthetic data generation. To address these issues, it would be helpful to incorporate additional validation methods that test the sensitivity of the probing strategy to known, context-dependent variations within well-studied GRNs. For instance, starting with a validated GRN and introducing experimentally confirmed context-specific variations could help determine if the models correctly adapt the network structure. Although this approach still has limitations (e.g., it could have been in the training set of GPT4), it could provide orthogonal evidence of the probing strategy’s effectiveness in identifying/retrieving context-sensitive elements.

Similarly, for the initial TF lists generated by foundation models, I recommend testing their context specificity by comparing TF lists for pairs of contexts with low TF overlap based on expression in resources like the Human Cell Atlas. This could reveal whether the models distinguish between TFs active in distinct biological contexts.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. How might the LLM-based approach handle incomplete pathway information or limited TF-gene knowledge in certain datasets?
2. How does pathway information specifically improve GRN inference compared to simpler approaches that omit it?
3. What criteria were used to select the three datasets, and how might the results generalize to datasets with different biological characteristics?
4. How does the complexity of LLM4GRN impact reproducibility, and are there ways to simplify the model without losing performance?
5. Can you discuss how robust the model is to potential biases in the LLM-based knowledge base, especially given the dependency on prior information?
6. Why do you believe Llama performs better than GPT-4 in certain settings, and what does this imply for LLM-based GRN inference?
soundness: 2
strengths: 1. Innovative Use of LLMs: Leveraging LLMs in GRN discovery demonstrates novel application potential in biological data analysis.
2. Comparison of Multiple LLMs: Provides comparative insights into GPT-4 and Llama-3.1 performances, showing valuable differences that inform model selection for scRNA-seq analysis.
3. Promising Hybrid Approach: Results show that combining LLM-derived TF lists with statistical models like GRNBoost2 can improve causal inference quality.
summary: This paper presents LLM4GRN, a method that integrates Large Language Models (LLMs) to discover causal Gene Regulatory Networks (GRNs) from single-cell RNA sequencing (scRNA-seq) data. The study employs a causal synthetic data generation approach, comparing two experimental setups where LLMs either generate GRNs directly or assist statistical methods (e.g., GRNBoost2) by providing transcription factor (TF) candidates. The authors test the approach on three datasets, highlighting LLMs' effectiveness in biological plausibility and similarity to real data. Notably, Llama-3.1 outperformed GPT-4 in some aspects, suggesting potential advantages for using certain open-source models over more advanced commercial ones.
weaknesses: 1. Limited Dataset Variety: Only three datasets (PBMC-All, PBMC-CTL, and Bone Marrow) are used, which limits generalizability to broader biological contexts.
2. Dependency on Pathway Completeness: The method’s performance may suffer in cases with incomplete pathway information, which isn’t fully addressed.
3. Biological Relevance of Generated Data: Some results indicate noise and inconsistency in cell-type differentiation, which could impact downstream analyses and applications.
4. Complexity of Methodology: The use of multiple knowledge bases and LLM models creates a complex pipeline that may be difficult for other researchers to reproduce.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Could the authors elaborate on how the performance of the LLM-based GRN inference might be affected by the quality and representativeness of the data used to train the LLMs?
2. It would be helpful if the authors could discuss potential future work on refining the GRN inference process, especially in addressing the noise and specificity issues mentioned in the biological plausibility analysis.
3. The paper mentions the use of GPT-4 and Llama-3.1-70B models. Are there any plans to compare the approach with other LLMs or to investigate the impact of model size on the performance of GRN inference?
soundness: 2
strengths: 1. The methodology involving the use of causal synthetic data generation to evaluate GRNs is well-designed and addresses a significant challenge in the field, namely the lack of ground truth data for validation.
2. The paper is well-organized, with clear explanations of the methodology, experimental setup, and results. The use of visual aids such as figures and tables enhances the understanding of the presented material.
3. The potential of LLMs to capture complex biological interactions and support statistical modeling in genomics is a significant contribution, with broad implications for disease mechanism understanding and therapeutic target identification.
summary: The paper titled "LLM4GRN: Discovering Causal Gene Regulatory Networks with LLMs – Evaluation through Synthetic Data Generation" presents a novel approach leveraging large language models (LLMs) for inferring gene regulatory networks (GRNs) from single-cell RNA sequencing (scRNA-seq) data. The authors propose using LLMs to generate complete GRNs or to provide prior knowledge to traditional statistical methods. A key contribution is the use of causal synthetic data generation as an evaluation strategy, which allows for the assessment of the inferred GRNs in the absence of ground truth causal graphs. The study demonstrates the potential of LLMs in supporting statistical modeling and data synthesis for biological research, with a particular emphasis on the combination of LLMs and statistical inference methods showing promise for scRNA-seq data analysis.
weaknesses: 1. The paper could benefit from a more detailed discussion on the limitations of the current approach, particularly regarding the assumptions made in the GRN structure due to the constraints of the GRouNdGAN framework.
2. While the study highlights the potential of LLMs, it does not explore the reasons behind the performance differences observed between the LLM-based and statistical methods, which could provide further insights into the underlying mechanisms.
3. The paper might consider discussing the scalability of the proposed method, especially how it might perform with larger and more complex datasets.
4.The paper does not seem to have much theoretical innovation. It simply uses LLM as a plug-in to replace previous human prior knowledge or statistical analysis methods.
5. The results are not very convincing: Figure 2 aims to illustrate the robustness and superiority of the method by overleaping different methods. However, the small gap cannot reflect the advantages of the proposed method.

Meta-review: This paper investigates leveraging LLMs to discover Gene Regulatory Network (GRN) from scRNA-seq data, evaluated via synthetic data generation. While reviewers appreciated the novelty and potential impact, concerns were raised about limited datasets, lack of robustness to noise and biases, minimal benchmarking against validated GRNs, and marginal performance improvements. The authors addressed some concerns by clarifying the theoretical basis, adding ablation studies, and outlining future directions, but many reviewers found these responses insufficient. Overall, the paper presents an interesting idea but requires further evaluation and refinement.

Decision: Reject