1. Introduction

Organismal development and pathogenesis emerge from complex genetic and cellular interactions. Understanding the developmental trajectory of multicellular organisms and predicting pathogenesis are pivotal goals of virtual cell modeling. Over the past two decades, several approaches to virtual cells have been proposed. A conventional approach is modeling the target biological variables with ordinary differential equations^[1], which is computationally intractable when the model becomes large. Another strategy is to use interpretable rule-based models that update cell states and predict cell behaviors according to predefined rules, with internal cellular states encoded as functions of environmental signals. However, these methods generally struggle with complex phenotypes involving nonlinear interactions, undefined cell types, unknown behaviors, spatial organization, temporal dynamics, and tissue-level characteristics in clinical data^[2]. In contrast, systems biology provides a holistic framework by viewing cells as systems of interacting genes and proteins^[3]. Signaling pathways and gene regulatory networks (GRNs) are inferred from gene expression profiles of perturbed cells and used as low-rank approximations of cellular states. Such models have provided valuable insights into transcriptional responses to genetic perturbations or mutations. However, several fundamental limitations make these simulations inadequate. In complex systems, perturbations often have context-dependent, non-additive effects, while conventional models are often static, based on local linear approximations, and do not consider spatial or cell-type heterogeneity^[4]. GRNs have achieved some success in holistic modeling but are often limited to being a complex and abstract representation without clear system boundaries^[5]. Recently, researchers have started to train virtual cell models by feeding vast amounts of single-cell RNA-seq data into deep generative models or agents. Although these models have shown great value in predicting cellular responses, their black-box nature makes it difficult to explicitly explain the mapping of genes to phenotypic effects and sometimes leads to uncontrolled model behaviors^[6]. After all, the detailed principles bridging individual genes and the highly organized biological structures remain poorly understood. Although efforts have been made to organize gene sets to functional modules or “pathways”, these concepts remain nebulous owing to the lack of rigorous defining criteria and validation approaches^[7]. Therefore, it is necessary to find a new anchor point for the next-generation virtual models.

2. Gene Programs as a Robust Interpretable Layer for Virtual Modeling

The advent of single-cell genomics has enabled the extraction of “gene programs” (GPs). Originating in the 1960s, this concept shares its roots with GRNs, positing that biomolecules operate modularly^[8]. The discovery of GPs, driven primarily by component analysis or non-negative matrix factorization, allows covariant gene expression patterns to naturally emerge from large datasets^[9]. These programs are highly interpretable and functionally specific representations of cellular systems. For example, Ling et al. identified a neuron-astrocyte concerted program that is associated with reduced schizophrenia risk and physiologically downregulated in older individuals^[10]. Echoing Barabási’s view of diseases as breakdowns of biological modules, this program is enriched for known schizophrenia risk genes, which suggests that aberrant program dynamics might drive pathogenesis^[11]. Furthermore, Gavish et al. found that diverse tumor categories shared common GPs that form clusters with defined biological functions^[12]. These consistent gene sets, which they called meta-programs (MPs), are de-noised forms of GPs and may reflect fundamental states of cells. In developmental biology, Farrell et al. discovered that statistically inferred GPs form tree-like structures that mirror the cell differentiation trajectories during embryogenesis^[13]. Transcription factors (TFs) of GPs can be inferred using SCENIC^[12,14], indicating that GPs are coherent functional modules emerging from noisy biological networks and establishing GPs as robust frameworks for genetic research^[14,15].

Unlike conventional GRNs, which suffer from high noise, ambiguous boundaries, and computational difficulties^[8], GPs are robust to data heterogeneity and can be easily used, once established, for downstream tasks. For example, Zhang et al. showed that established MPs can be used as robust anchors for highly heterogeneous data integration^[16]. Puram et al. found GPs directly correlated with tumor metastasis and cancer stages^[17]. Other uses of established GPs include cell state scoring (for consistent cell type identification), and identifying new cellular activities^[18,19]. These findings establish GPs as reusable and robust representations of biological systems. Crucially, Ota et al. demonstrated that GPs sometimes act as the “intermediate layer” translating genotype to phenotype^[20]. Individual traits, such as mean corpuscular hemoglobin, are found to be bidirectionally regulated by programs related to hemoglobin synthesis, cell cycle, and autophagy across various cell types. This implies that GPs may function as the intermediate layer bridging the genome and cellular phenotypes and can potentially be used as computational nodes in virtual cell computation. In this context, we use “computational nodes” to refer to computable model units that represent GP activities and can receive perturbation inputs and update downstream cellular states or phenotypic outputs. Depending on the modeling architecture, such nodes may be implemented as interpretable latent dimensions, graph nodes, or constrained neural modules. For example, linearly decoded variational autoencoders restrict the decoder to a linear configuration, where GPs serve as interpretable variables within the latent space^[21]. Similarly, but further extended to spatial modeling, the tool INSPIRE uses graph neural networks to find consistent GPs to integrate non-adjacent spatial transcriptomics (ST) data^[22], where GPs are modeled as global basis vectors jointly optimized with the network. By mathematically grounding GPs as either latent variables or structured factors, these approaches provide a blueprint for embedding biological knowledge into computational models, offering capabilities to predict cell states in response to genetic perturbations or environmental signals. Specifically, genetic perturbations can be encoded as vectors over genes, such as one-hot encoding, and linearly projected onto the GP space using gene-to-GP membership weights. Pathway-level perturbations can be encoded similarly as vectors over the corresponding genes. In the simplest linear approximation, a loss-of-function or over-expression perturbation can be modeled as a negative or positive change in the perturbed gene value, which is then projected to GP-level changes according to gene-to-GP membership weights. The updated GP states can then be used to predict the cell responses. In a neural network implementation, perturbations may also be first represented as gene embeddings or functional gene signature representations derived from Perturb-seq or drug-response datasets^[23,24]. Subsequently, perturbation effects can be propagated through GP regulatory relationships, spatial cell-cell interaction (CCI) graphs, or through learned state-update functions in different architectures. Here, propagation refers to transmitting the perturbation effects from direct molecular targets to downstream GPs, cell states, or spatial neighborhoods. For example, graph-enhanced gene activation and repression simulator (GEARS) uses graph neural networks trained on Perturb-seq datasets to predict single- or multi-gene perturbation outcomes, where the perturbation effects are propagated through a gene co-expression graph and a pathway-based gene knowledge graph^[23]. Similar perturbation encoding and propagation strategies could be adapted to a GP-centered framework.

3. Beyond Single Cells: Multi-Modal and Spatial Dimensions

Although GPs are primarily inferred from gene expression data and are sometimes viewed as alternatives to network modules or pathways, multi-omics technologies have provided evidence that certain GPs are linked to cross-layer functional modules with identifiable biophysical underpinnings^[25]. TFs, their binding motifs, and their experimentally verified regulons have been recurrently discovered by multi-omics sequencing and factor analysis^[26,27]. This is especially significant in simple systems such as E. coli, where GPs inferred by component analysis can be directly mapped to known regulons^[28]. Furthermore, existing evidence suggests that certain GPs, such as the one discovered by Puram et al., are downstream effects of specific ligand-receptor interactions^[29]. Regarding more remote molecular layers, such as metabolites and 3D chromatin structure, while there are some studies proposing their correlations with GPs within specific biosynthetic pathways and chromatin loops^[30,31], the precise and direct relationships among these structures remain an active area of investigation. These disparate molecular layers also complement the inherent noise and the temporal instability of transcriptomic state by providing modality-specific programs and cross-modality validation of top-weighted genes of GPs, serving as the cornerstone of a multi-level virtual model^[32].

Beyond cellular programs, the manifestation of biological complexity resides in the spatial architecture of tissues and intercellular communication. The integration of ST with GP inference has unveiled how GPs coordinate across spatial niches to form “multicellular programs” or “tissue programs”. For instance, Miller et al. discovered that meaningful GPs can be directly inferred from ST data by component analysis, and some of these GPs share a similar functional interpretation with single-cell GPs^[15]. After that, by mapping single-cell GPs to spatial datasets, they found that an immunosuppressive myeloid program within the tumor microenvironment is colocalized with tumor mesenchymal programs and physically entangled with angiogenesis or hypoxia niches defined by spatial GPs, serving as a critical prognostic factor for glioblastoma patients. Shi et al. discovered multicellular programs and demonstrated their strong correlation with aging and tumor progression by ST and in vivo perturbation data^[33]. Furthermore, spatial CCI analyses on ST datasets have begun to move beyond pairwise ligand-receptor interactions by incorporating sender-receiver gene co-expression and pathway-level contexts^[34-36], which motivates a GP-centered integration in which GP activities are analyzed together with spatial CCI networks to link niche-level phenotypes with candidate communication routes. These examples hint that spatial GPs may be associated with CCI, cell migration, and disease progression in processes like tumor invasion and drug resistance. Consequently, modeling virtual cells in microenvironments similar to in vivo states requires incorporating spatial context into GP extraction (Figure 1).

Display Full Size

Download

Figure 1. Schematic overview of gene program-based virtual cell modeling. Left, noisy biological networks are modularized from single-cell multi-omics data to yield robust and interpretable GPs, thereby reducing computational complexity while preserving biologically meaningful structure. Representative methods for GP extraction include NMF, ICA, and LDVAE. Right, GPs are encoded as explicit computational nodes, such as predefined latent variables in deep learning models. This framework accepts inputs from intrinsic or extrinsic perturbations and update cell states. Within this framework, intrinsic or extrinsic perturbations can be encoded as model inputs to update GP activities and downstream cell states. This enables interpretable genotype-to-phenotype mapping or in silico prediction of perturbation responses. By extending the framework to multiple interacting virtual cells and spatial cell-cell communication, it may also support simulation of tissue microenvironments in a manner closer to in vivo contexts. GPs: gene programs; NMF: non-negative matrix factorization; ICA: independent component analysis; LDVAE: linearly decoded variational autoencoder; ST: spatial transcriptomics.

4. Next-Generation Virtual Cells and Patient Digital Twins

While pioneering studies are already leveraging GPs alongside large-scale single-cell datasets such as Tahoe-100M to construct virtual cells^[37], primary bottlenecks include a scarcity of experimentally validated programs, inherent noise in perturbation profiles, and the limitations of current model architectures in encoding GPs as explicit variables^[20]. Integrating a broader array of programs as prior knowledge, along with their long-term temporal dynamics and multi-modal features, will establish a comprehensive toolkit and enable simulations of multicellularity, multi-perturbation responses, and continuous cell dynamics^[38,39]. Furthermore, how to assess the simulation quality is another critical issue that should be addressed. Technically, all criteria should be based on high-quality ground truth data, such as spatially resolved single-cell transcriptomics data for benchmarking multicellular simulations^[40], multi-gene perturbation experiments for prediction of multi-perturbation responses^[41], and temporally resolved transcriptomics for continuous cell dynamics modeling^[42]. Using these datasets, researchers can conduct multi-level assessment against the ground truth, such as the dynamics of GPs, marker genes associated with specific clinical outcomes, and cell migration, and compare model performance across different datasets and in multi-center validation. A useful model should be robust in real-world scenarios and require minimal retuning between scenarios.

The ultimate clinical promise of virtual cell technology lies in the construction of highly personalized “Patient Digital Twins”. Clinically derived multi-modal biopsy data are inherently noisy and heavily confounded by batch effects^[43]. In this context, GPs, serving as robust biological representations, emerge as the ideal framework to align macroscopic clinical phenotypes with microscopic molecular states. By leveraging GPs to extract features from an individual’s multi-omics profile, generalized virtual organ models can be rapidly fine-tuned into digital twins that accurately recapitulate the patient-specific pathophysiological states^[44]. This empowers clinicians to execute counterfactual simulations in silico prior to physical intervention, predictively evaluating the efficacies and toxicities of therapeutics. Thus, anchoring virtual biology on the foundation of GPs represents a practical route. We anticipate that a broad range of fields such as systems medicine, genetics, and oncology will be further propelled by this framework, moving towards mechanism-driven, highly individualized precision health.

Authors contribution

Wang T: Conceptualization, formal analysis, writing-original draft, writing-review & editing.

Li H, Liao J: Formal analysis, writing-review & editing.

Conflicts of interest

The authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

Not applicable.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 82522092).

Copyright

© The Author(s) 2026.