# Generating Multi-Modal and Multi-Attribute Single-Cell Counts with CFGen

Alessandro Palma1,2, Till Richter1,2, Hanyi Zhang1,2, Manuel Lubetzki1,
Alexander Tong3,4, Andrea Dittadi1,2,5, Fabian J. Theis1,2
1Helmholtz Munich, 2Technical University of Munich,
3Université de Montréal, 4Mila, 5MPI for Intelligent Systems, Tübingen
Correspondence to<EMAIL_ADDRESS>
###### Abstract

Generative modeling of single-cell RNA-seq data has shown invaluable potential
in community-driven tasks such as trajectory inference, batch effect removal
and gene expression generation. However, most recent deep models generating
synthetic single cells from noise operate on pre-processed continuous gene
expression approximations, ignoring the inherently discrete and over-dispersed
nature of single-cell data, which limits downstream applications and hinders
the incorporation of robust noise models. Moreover, crucial aspects of deep-
learning-based synthetic single-cell generation remain underexplored, such as
controllable multi-modal and multi-label generation and its role in the
performance enhancement of downstream tasks. This work presents Cell Flow for
Generation (CFGen), a flow-based conditional generative model for multi-modal
single-cell counts, which explicitly accounts for the discrete nature of the
data. Our results suggest improved recovery of crucial biological data
characteristics while accounting for novel generative tasks such as
conditioning on multiple attributes and boosting rare cell type classification
via data augmentation. By showcasing CFGen on a diverse set of biological
datasets and settings, we provide evidence of its value to the fields of
computational biology and deep generative models.

##  1 Introduction

Single-cell transcriptomics has enabled unprecedented exploration of cellular
heterogeneity, unveiling essential biological processes and fine-grained
cellular states [1, 2]. Advances in the single-cell RNA-seq (scRNA-seq)
technology allow parallel gene expression profiling in millions of cells,
leading to large, complex datasets that investigate cellular differentiation
[3, 4], disease progression [5, 6], and response to drug perturbation [7].
Aware of the multi-faceted complexity of the molecular state of a cell, modern
studies seek to complement the quantification of single-cell gene expression
with additional measurements, such as information on the accessibility state
of DNA regions encoding for molecular functions like gene expression and
regulation [8, 9]. This accessibility is measured with the Transposase-
Accessible Chromatin assay (ATAC-seq) [10], where open regions are referred to
as peaks. Multi-modal single-cell datasets, which offer a broader snapshot of
the biological state of cellular populations, have grown recently. Yet,
technical bias and high experimental costs still hinder the homogeneous
profiling of all possible cell states within the inspected biological process.
In light of this, generative modeling offers a promising avenue to highlight
underexplored biological states, drive scientific hypotheses, and enhance
downstream applications [11].

Generative models for single-cell data, in particular Variational Autoencoders
(VAEs), have been extensively employed in representation learning [11, 12],
perturbation prediction [13, 14, 15] and trajectory inference [16, 17, 18].
Recently, more complex approaches leveraging diffusion-based models [19] or
Generative Adversarial Networks (GAN) [20] have paved the way for the task of
synthetic data generation, demonstrating promising performance on rare cell
type augmentation and realistic single-cell data modeling. Single-cell data is
inherently discrete, as gene expression is collected as the number of
transcribed gene copies found experimentally. Due to the incompatibility of
discrete data with continuous models such as Gaussian diffusion [21], most
approaches generate data that has been pre-processed through normalization and
scaling. This limits their flexibility to support downstream tasks centered
around raw counts, such as batch correction [11], the use of differential
expression tools [22], and analyses where the total number of transcripts per
cell is important [3]. Additionally, technical and biological effects in
single-cell counts have been formalized under effective discrete noise models
[23], which should be incorporated into generative models for single-cell data
to approximate the underlying data generation process better.

Figure 1: The CFGen generative model. Gaussian noise from a prior $p_{0}$ is
converted into a latent variable $\mathbf{z}_{1}$ sampled from a target
density $p_{1}$ by a compositional flow conditioned on multiple biological and
technical attributes. Decoders for gene expression and DNA accessibility map
$\mathbf{z}_{1}$ to the parameters of negative binomial and Bernoulli
likelihood models from which single-cell gene expression and DNA accessibility
peaks are sampled.

In this work, we present Cell Flow for Generation (CFGen) (Fig. 1), a
conditional flow-based generative model designed to reproduce multi-modal
single-cell discrete counts realistically. Our approach combines the
expressiveness of recent Flow Matching techniques [24, 25, 26, 27, 28] with
modeling the statistical properties of single-cell data across multiple
modalities, each following a distinct discrete likelihood model. Moreover, we
extend the current literature on Flow Matching by introducing the concept of
compositionality, enabling the generation of cells conditioned on single
attributes or combinations thereof in a controlled setting.

We demonstrate the value of our approach in terms of generative performance
improvement and the enhancement of downstream tasks, such as data augmentation
for rare cell type classification. In summary, we introduce the following
contributions:

  * •

We propose CFGen, a generative model for discrete multi-modal single-cell
counts, which explicitly accounts for key statistical properties of single-
cell data under a specified noise model.

  * •

We extend the Flow Matching framework to incorporate guidance for
compositional generation under multiple attributes.

  * •

We demonstrate our model’s full-genome generative performance on a diverse set
of biological datasets compared to existing single-cell generative models.

  * •

We showcase potential applications of CFGen to data augmentation to improve
the generalization of cell type classification models.

##  2 Related work

The synthetic generation of single-cell datasets is a well-established
research direction pioneered by models using standard probabilistic methods to
estimate gene-wise parameters in a single modality [29, 30] or multiple
modality setting [31]. With the advent of deep generative models, VAE-based
approaches have demonstrated an extremely flexible performance, offering
popular tools for batch correction [11], modality integration [32], trajectory
inference [16] and perturbation prediction [13, 14, 15]. Despite their
relevance, most of the mentioned approaches orient towards learning meaningful
cellular representations or counterfactual predictions rather than generating
synthetic datasets from noise. Such a task has instead been extensively
explored by other works leveraging the expressive potential of diffusion
models [19], Generative Adversarial Networks (GANs) [20] and Large Language
Models (LLMs) [33] to produce realistic cells approximating the observed data
distribution. Our technical contribution builds upon Flow Matching [24, 25,
26], an efficient formulation of continuous normalizing flows for generative
modeling. Since its release, Flow Matching has been successfully applied to
optimal transport [28, 34, 35], protein generation [36, 37], interpolation on
general geometries [38, 39] and guided conditional generation [40]. Finally,
Flow Matching showed promising performance in tasks involving scRNA-seq, such
as learning cellular evolution across time [28, 39] and responses to drugs
[41].

##  3 Background

###  3.1 Deep generative modeling of single-cell data

Single-cell data is inherently discrete. At the gene expression level, cells
are high-dimensional vectors of counts, where each feature represents a gene
with, as a measurement, the number of transcripts detected in a cell.
Additionally, technical bias and biological variation result in cells
displaying unique statistical characteristics, such as sparsity and over-
dispersion. Sparsity derives from the absence of expression of genes in a
given cellular state (biological cause) or measurement dropouts in scRNA-seq
(technical cause). Over-dispersion entails that the gene-wise variance is
higher than the mean. From a probabilistic perspective, it is common to model
over-dispersed counts with a Negative Binomial (NB) distribution parameterized
by a mean $\mu$ and an inverse dispersion parameter $\theta$.

Formally, given a nonnegative count expression matrix
$\mathbf{X}\in\mathbb{N}_{0}^{N\times G}$ with $N$ cells and $G$ genes,
entries $x_{ng}$ of the expression matrix are assumed to follow the NB model:

$x_{ng}\sim\mathrm{NB}(\mu_{ng},\theta_{g})\>,$ (1)

where $\mu_{ng}\in\mathbb{R}_{\geq 0}$ is a cell-gene-specific mean and
$\theta_{g}\in\mathbb{R}_{\geq 0}$ the gene-specific inverse dispersion. In
other words, we assume each cell has an individual mean expression, while the
over-dispersion property is modeled for each gene.

When scRNA-seq is coupled with information on DNA accessibility, transcription
measurements are complemented by a binary matrix
$\mathbf{B}\in\\{0,1\\}^{N\times P}$, where $P$ is the number of DNA regions
profiled for accessibility. Here, each dimension independently follows the
Bernoulli model $b_{np}\sim\mathrm{Bern}(\pi_{np})$, with $\pi_{np}$
indicating a cell-gene-specific success probability.

In most single-cell representation learning settings, a deep latent variable
model is trained to map a latent space $\mathcal{Z}$ to the parameter space
$\mathcal{H}$ of the noise model via a decoder maximizing the log-likelihood
of the data. In practice, given a latent cell state $\mathbf{z}$, the
likelihood parameters of each modality are inferred as

$\displaystyle\bm{\mu}=l\bm{\rho}\>,\quad\bm{\rho}=\mathrm{softmax}(h_{\psi}(\mathbf{z}))\>,\quad\bm{\pi}=\mathrm{sigmoid}(h_{\phi}(\mathbf{z}))\>,$
(2)

where $h_{\psi}$ and $h_{\phi}$ are modality-specific decoders and $l$ is the
size factor, defined as the total number of counts of the generated cell.
Notably, the vector $\bm{\rho}$ represents gene expression proportions. One
can interpret the parameterization of the scRNA-seq likelihood based on the
connection between the negative binomial and Poisson-gamma distributions (see
Section B.1).

###  3.2 Continuous normalizing flows and Flow Matching

Continuous Normalizing Flows (CNF). Chen et al. [42] introduced CNFs as a
generative model to approximate complex data distributions. Given data in a
continuous domain $\mathcal{Z}\subset\mathbb{R}^{d}$, one can define a time-
dependent probability path $p:[0,1]\times\mathbb{R}^{d}\to\mathbb{R}_{\geq
0}$, transforming a tractable prior density $p_{0}$ into a more complex one
$p_{1}$, where we indicate the probability path at time $t$ as
$p_{t}:\mathbb{R}^{d}\to\mathbb{R}_{\geq 0}$ such that $\int
p_{t}(\mathbf{z})\,\textrm{d}\mathbf{z}=1$. The probability path is formally
generated by a time-dependent vector field,
$u_{t}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ with $t\in[0,1]$, if the continuity
equation $\frac{\mathrm{d}p}{\mathrm{d}t}=-\nabla\cdot(p_{t}u_{t})$ holds. The
field $u_{t}$ is the time-derivative of an invertible diffeomorphism
$\phi_{t}:[0,1]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, also called flow,
following the Ordinary Differential Equation (ODE):

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\phi_{t}(\mathbf{z})=u_{t}(\phi_{t}(\mathbf{z})),\quad\phi_{0}(\mathbf{z})=\mathbf{z}\,.$
(3)

Learning a CNF means defining a push-forward operator
$p_{t}=[\phi_{t}]_{*}p_{0}$, transforming the prior $p_{0}$ into the data
density $p_{1}$. In other words, learning the target vector field $u_{t}$
generating $p_{t}$ allows transporting samples from the prior distribution to
realistic data samples by simulating Eq. 3 in time.

Flow Matching [26]. Assume the goal is to model a complex data distribution
$q$ by learning the target vector field described in Section 3.2. One can
marginalize the probability path $p_{t}$ as $p_{t}(\mathbf{z})=\int
p_{t}(\mathbf{z}|\mathbf{z}_{1})q(\mathbf{z}_{1})\textrm{d}\mathbf{z}_{1}$,
where $\mathbf{z}_{1}$ indicates a sample from the data distribution $q$ and
$p_{t}(\cdot|\mathbf{z}_{1})$ is the conditional probability path transporting
noise to $\mathbf{z}_{1}$ under the boundary conditions
$p_{0}(\mathbf{z}|\mathbf{z}_{1})=p_{0}(\mathbf{z})$ and
$p_{1}(\mathbf{z}|\mathbf{z}_{1})\approx\delta(\mathbf{z}-\mathbf{z}_{1})$.
Here, $\delta$ represents a delta probability with mass concentrated at the
point $\mathbf{z}_{1}$. Note that, at $t=1$, the marginal distribution $p_{1}$
approximates the data distribution $q$. Moreover, $p_{t}(\mathbf{z})$ is
generated by a marginal velocity field $u_{t}(\mathbf{z})$ satisfying

$u_{t}(\mathbf{z})=\int
u_{t}(\mathbf{z}|\mathbf{z}_{1})\,\frac{p_{t}(\mathbf{z}|\mathbf{z}_{1})q(\mathbf{z}_{1})}{p_{t}(\mathbf{z})}\,\textrm{d}\mathbf{z}_{1},$
(4)

where $u_{t}(\mathbf{z}|\mathbf{z}_{1})$ is called conditional vector field.
While directly regressing $u_{t}(\mathbf{z})$ is intractable, Lipman et al.
[26] show that minimizing the Flow Matching objective

$\mathcal{L}_{\textrm{FM}}=\mathbb{E}_{t\sim\mathcal{U}[0,1],q(\mathbf{z}_{1}),p_{t}(\mathbf{z}|\mathbf{z}_{1})}\\!\left[||v_{t,\xi}(\mathbf{z})-u_{t}(\mathbf{z}|\mathbf{z}_{1})||^{2}\right]$
(5)

corresponds to learning to approximate the marginal vector field
$u_{t}(\mathbf{z})$ with the time-conditioned neural network $v_{t,\xi}$.
Defining
$p_{t}(\mathbf{z}|\mathbf{z}_{1})=\mathcal{N}(\alpha_{t}\mathbf{z}_{1},\sigma_{t}^{2}\textbf{I})$
with the functions $\alpha_{t},\sigma_{t}$ controlling the noise schedule,
$u_{t}(\mathbf{z}|\mathbf{z}_{1})$ has a closed form, and Eq. 5 is tractable
(see Section B.2 for more details). We define such a formulation Gaussian
marginal paths.

Classifier-free guidance [40]. One can guide data generation on a condition
$y$ by learning the conditional marginal field $u_{t}(\mathbf{z}|y)$ via a
neural network $v_{t,\xi}(\mathbf{z}|y)$. Given a guidance strength
hyperparameter $\omega\in\mathbb{R}$, Zheng et al. [40] show that generating
data points following the vector field
$\tilde{u}_{t}(\cdot|y)=(1-\omega)u_{t}(\cdot)+\omega u_{t}(\cdot|y)$
approximates sampling from the distribution $\tilde{q}(\mathbf{z}|y)\propto
q(\mathbf{z})^{1-\omega}q(\mathbf{z}|y)^{\omega}$, where $q(\mathbf{z})$ and
$q(\mathbf{z}|y)$ are, respectively, the unconditional and conditional data
distributions. Guidance can thus be achieved by combining conditional and
unconditional vector fields learned simultaneously during training.

##  4 CFGen

Our objective is to define a latent Flow-Matching-based generative model for
discrete single-cell data, where each cell is measured through gene expression
and DNA accessibility. Our model, CFGen, is flexible: It can handle single and
multiple modalities. Moreover, it supports guiding generation conditioned on
single or combinations of attributes without needing to train a different
model for each. In what follows, we present the assumptions and generative
process formulation in the uni-modal and multi-modal settings. We additionally
introduce our approach to compositional guidance.

###  4.1 Uni-modal and single-attribute generation

Let $\mathbf{X}\in\mathbb{N}_{0}^{N\times G}$ be a single-cell matrix where we
indicate an observed single-cell vector as $\mathbf{x}\in\mathbb{N}_{0}^{G}$,
with $G$ being the number of genes. Additionally, let
$\mathbf{y}\in\mathbb{N}_{0}^{N}$ be a categorical attribute (for example,
cell type) associated with each observation. We also define
$l=\sum_{g=1}^{G}x_{g}$ as the size factor of an individual cell $\mathbf{x}$.

Standard formulation of the generative process. When the technical bias is
negligible, we define the standard CFGen setting as the following generative
model:

$p(\mathbf{x},\mathbf{z},l,y)=p(\mathbf{x}|\mathbf{z},l)p(\mathbf{z}|y,l)p(l)p(y)\;,$
(6)

where $\mathbf{z}$ is a continuous latent variable modeling the cell state,
and we assumed that (1) $\mathbf{x}$ is independent of $y$ conditionally on
$\mathbf{z}$ and $l$, and (2) $l$ is independent of $y$. Crucially, Eq. 6
provides a standard formulation for the generative process, but the
factorization choice is flexible and can be adjusted based on the properties
of the data (e.g. when the size factor and conditioning covariate are not
independent). We provide a flexible implementation supporting said variations.

Modeling the distributions in Eq. 6. Each factor of Eq. 6 is modeled
separately: $p(y)$ is a categorical distribution
$\mathrm{Cat}(N_{y},\bm{\pi}_{y})$ where $N_{y}$ is the number of categories
and $\bm{\pi}_{y}$ a vector of $N_{y}$ class probabilities, and
$p(l)=\mathrm{LogNormal}(\mu_{l},\sigma^{2}_{l})$. The parameters of $p(y)$
and $p(l)$ can be learned as maximum likelihood estimates over the dataset
(see Section C.4). Given an attribute class $y$ and size factor $l$ sampled
from the respective distributions, $p(\mathbf{z}|y,l)$ is a conditional
normalizing flow learned via Flow Matching with Gaussian marginal paths (see
Section 3.2, Section B.2 and Section C.2 for the architecture details)
transporting samples from a standard Gaussian distribution
$\mathcal{N}(\mathbf{0},\mathbf{I})$ to latent cell representations, where
$\mathbf{I}$ is the identity matrix. Finally, given a latent variable
$\mathbf{z}$ and a size factor $l$, $p(\mathbf{x}|\mathbf{z},l)$ is a negative
binomial distribution with mean parameterized by a decoder $h_{\psi}$ as in
Eq. 2 and inverse dispersion modeled by a global parameter $\bm{\theta}$. In
practice, $h_{\psi}$ and $\bm{\theta}$ are optimized before training the flow,
together with an encoder $f_{\eta}$ that maps the data to a latent space (more
details in Section C.1).

Sampling in practice. Following the standard formulation of CFGen in Eq. 6, to
generate a cell we first sample a size factor $l$ and a condition $y$ (the
latter can be chosen to generate from a specific class). Successively, we
derive $\mathbf{z}_{0}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ and simulate a
real latent data point via

$\mathbf{z}_{1}=\mathbf{z}_{0}+\int_{0}^{1}v_{t,\xi}(\mathbf{z}_{t}|y,l)\mathrm{d}t\>;$
(7)

where $v_{t,\xi}(\mathbf{z}_{t}|y,l)$ is a parameterized vector field learned
through Flow Matching. Finally, $\mathbf{z}_{1}$ is run through the decoder
$h_{\psi}$ and a cell vector is obtained by sampling
$\mathrm{NB}(l\,\mathrm{softmax}(h_{\psi}(\mathbf{z}_{1})),\bm{\theta})$.

Size factor as a technical effect. When $l$ is influenced by technical effect
under a categorical covariate $c\in\\{1,\ldots,C\\}$, we reformulate Eq. 6 as
$p(\mathbf{x},\mathbf{z},l,y)=\frac{1}{C}\sum_{c}p(\mathbf{x}|\mathbf{z},l)p(\mathbf{z}|y,l)p(l|c)p(y)p(c)$,
where we assume that $\mathbf{z}$ is independent of $c$ given $l$ (i.e., $l$
contains all necessary technical effect information to guide the flow), and
$y$ is independent of $c$. The last assumption derives from our choice of $y$
as an attribute encoding biological identity preserved across technical
batches.

###  4.2 Multi-modal and single-attribute generation

Let $\mathbf{X}$ and $\mathbf{y}$ be defined as in Section 4.1. In the multi-
modal setting, we have additional access to a binary matrix
$\mathbf{B}\in\\{0,1\\}^{N\times P}$ representing DNA region accessibility,
with $P$ being the number of measured peaks. Each sample is, therefore, a
tuple ($\mathbf{x}$, $\mathbf{b}$, $y$), where $\mathbf{x}$ and $\mathbf{b}$
are realizations of different discrete noise models (negative binomial and
Bernoulli). Following Eq. 2, both parameters of the negative binomial and
Bernoulli likelihoods are functions of the same latent variable $\mathbf{z}$,
encoding a continuous cell state shared across modalities. We write the
likelihood term in Eq. 6 as

$p(\mathbf{x},\mathbf{b}|\mathbf{z},l)\stackrel{{\scriptstyle\scriptscriptstyle(\mkern-1.5mu1\mkern-1.5mu)}}{{=}}p(\mathbf{x}|\mathbf{z},l)p(\mathbf{b}|\mathbf{z},l)\
\stackrel{{\scriptstyle\scriptscriptstyle(\mkern-1.5mu2\mkern-1.5mu)}}{{=}}p(\mathbf{x}|\mathbf{z},l)p(\mathbf{b}|\mathbf{z})\;$
(8)

where in (1) we use the fact that the likelihood of $\mathbf{x}$ and
$\mathbf{b}$ are optimized disjointedly given $\mathbf{z}$, and in (2) that
$\mathbf{b}$ is independent of the size factor $l$ (see Eq. 2). In simple
terms, all the modalities are encoded to the same latent space $\mathcal{Z}$
used to train the conditional flow approximating $p(\mathbf{z}|y,l)$ in Eq. 6
(more details in Section C.1). During generation, separate decoders $h_{\psi}$
and $h_{\phi}$ map a sampled latent state $\mathbf{z}$ to the parameter spaces
of the negative binomial and Bernoulli distributions, representing expression
counts and binary DNA accessibility information, respectively (Fig. 1).

###  4.3 Guided compositional generation with multiple attributes

We extend the concept of classifier-free guidance (CFG) for Flow Matching [40]
to account for multiple attributes, which is useful in applications to enhance
control of the generative process towards specific data regions. This is
relevant in scRNA-seq, where datasets are usually characterized by several
covariates encoding biological and technical information. Here,
$\mathbf{Y}\in\mathbb{N}_{0}^{N\times K}$ is an integer matrix of $K$
categorical attributes. Instead of learning a conditional flow model for every
combination of attributes as a separate conditioning class, we propose a
flexible approach to generate cells from arbitrary attribute subsets by
composing multiple flow models conditioned on single attributes.

Let $q(\mathbf{z}|\mathbf{y})$ be the conditional data distribution, with
$\mathbf{y}=(y_{1},\ldots,y_{K})$ being a collection of observed categorical
attributes. In analogy with CFG in diffusion models [43], we aim to implement
a generative model to sample from $\tilde{q}(\mathbf{z}|\mathbf{y})\propto
q(\mathbf{z})\prod_{i=1}^{K}\left[\frac{q(\mathbf{z}|y_{i})}{q(\mathbf{z})}\right]^{\omega_{i}}$,
where $\omega_{i}$ is the guidance strength for attribute $i$ (see Section
B.4). Diffusion models generate data by learning to approximate the score of
the time-dependent distribution, $\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|\mathbf{y})$, with a neural network and using it to simulate
a reverse diffusion Stochastic Differential Equation (SDE) transporting noise
samples to generated data observations [44]. Importantly, the reverse
diffusion SDE corresponds to a deterministic probability flow ODE with the
same time-marginal distributions [45].

Based on the relation between score-based diffusion and energy models, Liu et
al. [46] demonstrated that compositional classifier-free guidance is
achievable through expressing the drift of the generating reverse SDE with the
compositional score:

$\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|\mathbf{y})=\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\sum\nolimits_{i=1}^{K}\omega_{i}[\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})-\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z})]\>.$ (9)

Following the direct relationship between Flow Matching and CFG provided in
[43], we build the Flow Matching counterpart to Eq. 9.

Provided that we have access to Flow Matching models for the unconditional
marginal vector field $u_{t}(\mathbf{z})$ and the single-attribute conditional
fields $\\{u_{t}(\mathbf{z}|y_{i})\\}_{i=0}^{K}$, we introduce the following
proposition:

###### Proposition 1

If the attributes $y_{1},...,y_{K}$ are conditionally independent given
$\mathbf{z}$, the vector field

$\tilde{u}_{t}(\mathbf{z}|\mathbf{y})=u_{t}(\mathbf{z})+\sum\nolimits_{i=1}^{K}\omega_{i}[u_{t}(\mathbf{z}|y_{i})-u_{t}(\mathbf{z})]$
(10)

coincides with the velocity of the probability-flow ODE associated with the
generative SDE of a diffusion model with a compositional score as in Eq. 9.

We provide a proof for 1 in Section B.4. In other words, the reversed
diffusion SDE from compositional CFG admits a deterministic probability flow
ODE with velocity as in Eq. 10. Consequently, classifier-free sampling from
compositions of attributes is achievable by simulating the probability path
$\tilde{p}_{t}(\mathbf{z}|\mathbf{y})\propto
p_{t}(\mathbf{z})\prod_{i=1}^{K}\left[\frac{p_{t}(\mathbf{z}|y_{i})}{p_{t}(\mathbf{z})}\right]^{\omega_{i}}$
from a prior Gaussian density $p_{0}$ integrating the parameterized field
$\tilde{v}_{t,\xi}(\mathbf{z}|\mathbf{y})=v_{t,\xi}(\mathbf{z})+\sum_{i=1}^{K}\omega_{i}[v_{t,\xi}(\mathbf{z}|y_{i})-v_{t,\xi}(\mathbf{z})]$.
Note that both conditional and unconditional fields are parameterized by the
same model, which is learned by providing single-attribute conditioning with a
certain probability during training (see Algorithms 1 and 2).

##  5 Experiments

In this section, we first compare CFGen with existing models in terms of uni-
modal (Section 5.1) and multi-modal (Section 5.2) generation over five
distinct datasets. The comparison is performed quantitatively in terms of
distributional proximity between real and generated cells, and qualitatively
by examining how well the compared generative models approximate real data
properties. In Section 5.3, we illustrate the effectiveness of multi-attribute
generation to guide the synthetic samples to specific combinations of
biological labels and donors. Finally, in Section 5.4 we demonstrate that
targeted data augmentation, using our generative model to increase the
representation of rare cell types enhances the performance of LLM-based cell
type classifiers.

###  5.1 Comparison with existing methods

Table 1: Quantitative performance comparison between conditional and
unconditional single-cell generative models. Evaluation is performed based on
distribution matching metrics (MMD and 2-Wasserstein distance) and the
generated cell type classification F1 score achieved by a kNN classifier
trained on real data. Results are averaged across generations performed using
three different seeds.

| PBMC3K | Dentate gyrus | Tabula Muris | HLCA
---|---|---|---|---
| MMD $(\downarrow)$ | WD $(\downarrow)$ | KNNc $(\uparrow)$ | MMD $(\downarrow)$ | WD $(\downarrow)$ | KNNc $(\uparrow)$ | MMD $(\downarrow)$ | WD $(\downarrow)$ | KNNc $(\uparrow)$ | MMD $(\downarrow)$ | WD $(\downarrow)$ | KNNc $(\uparrow)$
| Conditional
c-CFGen | 203.58 ± 7.62 | 16.94 ± 0.40 | 0.36 ± 0.15 | 427.22 ± 5.22 | 18.74 ± 0.16 | 0.21 ± 0.03 | 72.97 ± 21.49 | 7.39 ± 0.20 | 0.15 ± 0.03 | 91.04 ± 1.64 | 10.72 ± 0.08 | 0.18 ± 0.02
scDiffusion | 403.09 ± 57.49 | 22.41 ± 1.21 | 0.22 ± 0.06 | 479.89 ± 4.16 | 19.70 ± 0.09 | 0.22 ± 0.03 | 74.10 ± 22.93 | 7.89 ± 0.45 | 0.12 ± 0.01 | 218.58 ± 14.23 | 15.82 ± 0.45 | 0.08 ± 0.01
scVI | 232.05 ± 7.53 | 17.66 ± 0.29 | 0.43 ± 0.15 | 476.91 ± 8.35 | 22.61 ± 0.23 | 0.19 ± 0.03 | 86.03 ± 23.17 | 9.76 ± 0.53 | 0.10 ± 0.02 | 111.31 ± 4.39 | 11.78 ± 0.19 | 0.09 ± 0.01
| Unconditional
u-CFGen | 204.24 ± 2.01 | 16.81 ± 0.06 | - | 419.44 ± 1.90 | 21.20 ± 0.02 | - | 38.05 ± 0.33 | 7.70 ± 0.03 | - | 84.84 ± 0.41 | 10.63 ± 0.01 | -
scGAN | 241.48 ± 3.30 | 15.54 ± 0.06 | - | 616.04 ± 0.10 | 22.52 ± 0.03 | - | 17.05 ± 0.22 | 12.70 ± 0.02 | - | 30.34 ± 0.31 | 10.80 ± 0.01 | -

In Table 1, we evaluate the generative performance of CFGen conditionally
(c-CFGen) and unconditionally (u-CFGen) against three baselines on the scRNA-
seq generation task.

Baselines. We choose scVI [32] and scDiffusion [19] as conditional models and
scGAN [20] as unconditional baseline. The scVI model is based on a VAE
architecture with a negative binomial decoder and performs generation by
decoding low-dimensional Gaussian noise into parameters of the likelihood
model. Conversely, scDiffusion and scGAN operate on a continuous-space domain,
performing generation using standard latent diffusion [47] and GAN [48]
models. Thus, we train them using normalized counts (more in Appendix E).

Datasets. We evaluate the performance of the models on four datasets with
variable size: (i)
PBMC3K111https://satijalab.org/seurat/articles/pbmc3k_tutorial.html (3k
Peripheral Blood Mononuclear Cells) measures the expression of 2,700 cells
from a healthy donor clustering into 8 cell types. (ii) The Dentate gyrus
dataset [49] profiles the expression of $18,213$ cells from a developing mouse
hippocampus. (iii) Tabula Muris [50] comprises the gene expression data of
245,389 cells from the mouse hippocampus across multiple tissues. (iv) The
Human Lung Cell Atlas (HLCA) [51]–previously unconsidered for generative
modeling–is a vast cohort of 2.4 million cells from 486 individuals harvested
from 49 datasets. Here, we use a size-controlled version of the dataset
comprising 584,944 cells. We perform conditioning on the cell type covariate
on all datasets but the Tabula Muris, where we use the tissue label. A
description of dataset characteristics and pre-processing is provided in
Appendix G and Table 3.

Evaluation. As evaluation metrics, we use distribution distances (MMD [52] and
Wasserstein-2 distance) computed between the Principal Component (PC)
projections of real and generated data in 30 dimensions. The generated data is
embedded using the PC loadings of the real data for comparability. For
conditional models, we additionally evaluate the cell-type classification F1
score of a k-Nearest-Neighbors (kNN) classifier trained on real data and
evaluated on generated cells. All evaluations are performed on a held-out set
of cells, considering the whole genome, with a filtering step for low
expression genes (see Appendix G).

Results. CFGen always reaches either the top or second-best performance on the
conditional generation task and overcomes scGAN in five out of eight metrics
measured across datasets. Qualitatively, CFGen is the only conditional model
capable of generating plausible synthetic cells in terms of overlap with the
real data distribution for large datasets such as the HLCA and Tabula Muris
(see Fig. 10 and Fig. 10).

Figure 2: (a) Comparison between the gene-wise mean-variance trend in real
data and samples from generative models. (b) Number of zeroes per cell
frequency in real data and samples from generative models.

When dealing with biology, measuring realistic data generation goes beyond
distribution matching metrics. The main methodological advantage of CFGen over
diffusion or GAN-based single-cell generators is the definition of a
generative model that explicitly considers the probabilistic properties of the
data. By learning a rigorous likelihood-based model of the discrete data (Eq.
6), we find that our generated cells are faithful to the real observations in
terms of the following properties:

(1) Sparsity. Caused by technical bias in the detection of gene transcripts or
by genes being unexpressed in specific biological settings. (2) Over-
dispersion. The gene-specific mean-variance trend of scRNA-seq is nonlinear;
genes (especially when highly expressed) have higher variance than expected
under a Poisson model (where mean and variance coincide) due to technical
effects in transcript detection [53]. The inverse dispersion parameter of a
negative binomial distribution models this behavior. (3) Discreteness and
skewness. scRNA-seq data is discrete and skewed towards zero-expression values
(see sparsity) [54].

In Fig. 2, we provide qualitative evidence that using CFGen is more
appropriate at recovering (1) and (2) than scDiffusion, which assumes a
continuous data space. (3) naturally holds when modeling discrete counts with
CFGen as opposed to a continuous model. In particular, Fig. 2a illustrates
that explicitly modeling counts with gene-specific inverse dispersion ensures
better overlap in the generated gene-wise mean-variance relationship with that
of the real data. Additionally, in Fig. 2b, we demonstrate an almost perfect
recovery of the real distribution of the number of zero counts per cell. At
the same time, scDiffusion tends to produce an unrealistic shift towards
actively expressed genes.

###  5.2 Multi-modal generation

Here, we evaluate the qualitative and quantitative performance of CFGen at
generating multi-modal data comprising gene expression and binary DNA-region
accessibility.

Baselines. We compare CFGen with a VAE-based multi-modal generative model
(MultiVI) [55]. For completeness, we include as baselines two single-modality
generative models: PeakVI [56] (DNA accessibility) and scVI [11] (gene
expression).

Datasets. We introduce the multiome PBMC10K dataset, made available by 10X
Genomics 222https://www.10xgenomics.com/support/single-cell-multiome-atac-
plus-gene-expression. Here, each cell is measured both in gene expression and
peak accessibility. The dataset consists of 10,000 cells across 25,604 genes
and 40,086 peaks and was annotated with 15 cell types, with their respective
marker peaks (enriched in accessible or inaccessible points) and genes.

Table 2: Comparison between CFGen and VAE-based models on generating multiple
single-cell modalities. PeakVI and scVI only generate accessibility and gene
expression, respectively. We report distribution distance performance (MMD and
Wasserstein-2 distance) between real and generated cells across three seeds.

| RNA | ATAC
---|---|---
| MMD ($\downarrow$) | WD ($\downarrow$) | MMD ($\downarrow$) | WD ($\downarrow$)
CFGen | 175.40 ± 1.35 | 14.54 ± 0.05 | 129.50 ± 7.71 | 18.92 ± 0.36
MultiVI | 194.30 ± 8.93 | 15.92 ± 0.25 | 204.75 ± 9.39 | 21.09 ± 0.34
PeakVI | - | - | 196.75 ± 11.48 | 20.84 ± 0.45
scVI | 177.37 ± 3.48 | 14.39 ± 0.12 | - | -

Evaluation. We use the MMD and Wasserstein-2 distances in the same setting
described in Section 5.1. Before comparison, we normalize both real and
generated binary measurements of DNA accessibility via TF-IDF [57] (in analogy
to text mining). The metrics are computed in a 30-dimensional PC projection of
the data. RNA counts are treated as in Section 5.1. To obtain a more
biologically-informed comparison between CFGen and MultiVI, we also evaluate
how well the two models approximate the per-cell-type marker peak
accessibility fraction and gene expression (more in Section H.5). In short, we
subset real and generated data to marker peaks and genes taken from the
literature for each cell type. For peaks, we compute the fraction of cells
from each cell type where each marker peak is accessible. For marker genes, we
calculate the mean gene expression of markers per cell type. Finally, we
evaluate the Pearson correlation between vectors of per-cell-type marker
accessibility fraction and average expression in real and generated datasets.

Figure 3: (a) 2D UMAP overlap between real and generated cells across
modalities on the PBMC10K dataset. (b) Pearson correlation between average
cell-type-specific marker peak accessibility and marker gene expression
between real data and samples from CFGen and MultiVI.

Results. CFGen outperforms both MultiVI and PeakVI in modeling accessibility
data based on distribution matching metrics (see Table 2). It also surpasses
MultiVI in the gene expression modality, producing results on par with scVI.
Thus, CFGen is the most accurate multi-modal generative model in this setting.
Qualitatively, Fig. 9a shows substantial overlap between real and generated
modalities. Finally, Fig. 9b demonstrates that CFGen better approximates
average marker peak accessibility and gene expression, outperforming MultiVI
across all cell type categories.

###  5.3 Multi-attribute generation and guidance

Here, we qualitatively evaluate our approach to the multi-attribute
compositional guidance described in Section 4.3.

Datasets. We showcase guidance on datasets with extensive batch effects, as
one could combine different levels of biological and technical annotations to
either augment rare cell type and batch combinations or control for the amount
of technical effect added in simulation settings. Specifically, we consider
(i) The NeurIPS 2021 [58] bone marrow datasets, comprising 90,261 annotated
single cells from 12 healthy human donors across four sites. We use donor as a
batch attribute and cell type as a biological covariate. We also consider (ii)
the Tabular Muris dataset described in Section 5.1, using tissue and Mouse ID
as covariates.

Evaluation. The power of our guidance model is to generate data conditionally
on an arbitrary subset of attributes—including unconditional generation—using
a single trained model. For each pair of covariates $(y_{i},y_{j})$, we
evaluate generation on 500 generated cells varying the parameter $\omega_{j}$,
keeping $\omega_{i}$ fixed. The expected result is conditional generation on
$\omega_{i}$ when $\omega_{j}=0$ and generation from the intersection between
the two attributes as $\omega_{j}$ increases. We additionally test the
unconditional model, given by $\omega_{i}=\omega_{j}=0$, expected to recover
the whole single-cell dataset. In the unconditional case, we generate as many
cells as there are in the dataset to better evaluate the coverage.

Figure 4: Qualitative evaluation of guidance performance on attribute pairs in
NeurIPS 2021 and Tabula Muris datasets. Left: unconditional performance with
guidance weights at 0. Moving right: simulate 500 cells, progressively
increasing the guidance strength of one attribute while keeping the
counterpart unchanged.

Results. Qualitative guidance results are shown in Fig. 4, with examples of
double-attribute guidance between CD14+ monocytes and donor 1 for the NeurIPS
2021 dataset and tongue and mouse 18-M-52 for Tabula Muris. Unconditional
generation faithfully recreates the original data (left-hand side) for both
datasets. Setting guidance weights to zero for mouse ID and donor attributes
leads to single-attribute conditional generation. Increasing the guidance
weight steers the generation towards the intersection of the two attributes.
Interestingly, different attribute pairs in different datasets require varying
guidance strengths for optimal intersection.

###  5.4 Application example: data augmentation to improve classification of
rare cell types

Finally, we propose a use case for CFGen: enhancing the generalization
performance of cell-type classifiers by augmenting rare cell types in the
dataset. Previous studies have demonstrated performance enhancement through
leveraging large auxiliary datasets [59]. As a cell type classifier, we select
scGPT [60], an LLM-inspired single-cell transformer pre-trained on a corpus of
33 million cells.

Datasets. We evaluate our approach on two datasets: (i) the PBMC COVID dataset
[61] profiles 422,220 blood cells from 93 patients ranging across paediatric
and adult. We also consider the (ii) HLCA dataset described in Section 5.1.
Both datasets are pre-processed by selecting 2000 highly variable features and
leaving out cells from 20% of the donors from training.

Evaluation. We train CFGen on the PBMC COVID and HLCA data training sets.
Then, we augment both datasets to 800,000 samples, upsampling rare cell types.
More in detail, we assign to each cell type a probability
$\frac{1}{N_{\mathrm{ct}}}$, where $N_{\mathrm{ct}}$ is the total number of
instances of said cell type. We then generate observations to cover the gap
between the size of the dataset and 800,000 cells, sampling cell types
proportionally to the assigned probabilities.

Figure 5: Cell-type classification difference as a function of cell type
frequency before and after augmentation on the PBMC COVID and HLCA datasets. A
10-neighbor kNN classifier is trained on the feature space of a pre-trained
scGPT model. The held-out set consists of cells from 20% of the donors.

At the end of the procedure, we obtain a considerably higher number of
observations for rare cell types. However, we still do not reach uniformity,
as class imbalance may be biologically meaningful. Following the original
publication, we then train kNN cell-type classifiers on scGPT’s embeddings of
the original and augmented training sets and evaluate the generalization
performance on the held-out donors. For each cell type, we assess if
performance increases upon augmentation as a function of its frequency in the
dataset.

Results. Results are displayed in Fig. 5 for the two datasets. Remarkably,
most cell types in the held-out dataset are better classified after
augmentation, suggesting that CFGen not only generates reliable cell samples
but can be a valuable supplement to relevant downstream tasks. Moreover, the
performance difference between before and after augmentation is inversely
proportional to the frequency of the cell type in the dataset, with a Pearson
correlation of $-0.31$ on the PBMC COVID and $-0.29$ on the HCLA datasets.
Therefore, the generalization improvement is more accentuated for rare cell
types.

##  6 Discussion

We presented CFGen, a conditional flow-based generative model for single-cell
discrete data. CFGen combines state-of-the-art generative models with rigorous
probabilistic considerations, incorporating established noise models for
sampling realistic gene expression and DNA accessibility states. We
demonstrated improved performance over existing generative models for single-
cell data, showing more faithful reproduction of data sparsity and
overdispersion than continuous counterparts. Moreover, we showed enhanced
performance in comparison with established generative frameworks for multi-
modal single-cell data considering gene expression and peak accessibility. Our
core machine learning contribution extends classifier-free guidance in Flow
Matching with compositional generation with multiple attribute instances,
inspired by the necessity to control both biological and technical variation
in the generation process. We also improved rare cell type classification in
the scGPT’s representation space through realistic augmentation of the
reference set. A limitation of our framework is its reliance on multiple
assumptions, including independence assumptions in the data, which may not
hold for all biological settings and require prior consideration of the
dataset properties. In other words, exploration of the data characteristics
should be carried out before using CFGen for generation. Overall, CFGen
represents a significant advancement in the simulation and augmentation of
single-cell data, offering the research community powerful tools for
generating realistic and varied datasets for improving performance on
biological problems.

##  7 Acknowledgments

A.P. and T.R. are supported by the Helmholtz Association under the joint
research school Munich School for Data Science (MUDS). A.P. and F.J.T.
acknowledge support from the German Federal Ministry of Education and Research
(BMBF) through the HOPARL project (grant number 031L0289A). T.R. and F.J.T.
also acknowledge support from the Helmholtz Association’s Initiative and
Networking Fund via the CausalCellDynamics project (grant number
Interlabs-0029). Additionally, F.J.T. acknowledges support from the European
Union (ERC, DeepCell - grant number 101054957).

## References

  * Hwang et al. [2018] Byungjin Hwang, Ji Hyun Lee, and Duhee Bang.  Single-cell rna sequencing technologies and bioinformatics pipelines.  _Experimental & molecular medicine_, 50(8):1–14, 2018.  ISSN 2092-6413.  doi: https://doi.org/10.1038/s12276-018-0071-8.  URL https://doi.org/10.1038/s12276-018-0071-8.
  * Rozenblatt-Rosen et al. [2017] Orit Rozenblatt-Rosen, Michael J. T. Stubbington, Aviv Regev, and Sarah A. Teichmann.  The human cell atlas: from vision to reality.  _Nature_ , 550(7677):451–453, October 2017.  ISSN 1476-4687.  doi: 10.1038/550451a.  URL http://dx.doi.org/10.1038/550451a.
  * Gulati et al. [2020] Gunsagar S. Gulati, Shaheen S. Sikandar, Daniel J. Wesche, Anoop Manjunath, Anjan Bharadwaj, Mark J. Berger, Francisco Ilagan, Angera H. Kuo, Robert W. Hsieh, Shang Cai, Maider Zabala, Ferenc A. Scheeren, Neethan A. Lobo, Dalong Qian, Feiqiao B. Yu, Frederick M. Dirbas, Michael F. Clarke, and Aaron M. Newman.  Single-cell transcriptional diversity is a hallmark of developmental potential.  _Science_ , 367(6476):405–411, January 2020.  ISSN 1095-9203.  doi: 10.1126/science.aax0249.  URL http://dx.doi.org/10.1126/science.aax0249.
  * Rizvi et al. [2017] Abbas H Rizvi, Pablo G Camara, Elena K Kandror, Thomas J Roberts, Ira Schieren, Tom Maniatis, and Raul Rabadan.  Single-cell topological rna-seq analysis reveals insights into cellular differentiation and development.  _Nature Biotechnology_ , 35(6):551–560, May 2017.  ISSN 1546-1696.  doi: 10.1038/nbt.3854.  URL http://dx.doi.org/10.1038/nbt.3854.
  * Zeng and Dai [2019] Tao Zeng and Hao Dai.  Single-cell rna sequencing-based computational analysis to describe disease heterogeneity.  _Frontiers in Genetics_ , 10, July 2019.  ISSN 1664-8021.  doi: 10.3389/fgene.2019.00629.  URL http://dx.doi.org/10.3389/fgene.2019.00629.
  * Chung et al. [2017] Woosung Chung, Hye Hyeon Eum, Hae-Ock Lee, Kyung-Min Lee, Han-Byoel Lee, Kyu-Tae Kim, Han Suk Ryu, Sangmin Kim, Jeong Eon Lee, Yeon Hee Park, Zhengyan Kan, Wonshik Han, and Woong-Yang Park.  Single-cell rna-seq enables comprehensive tumour and immune cell profiling in primary breast cancer.  _Nature Communications_ , 8(1), May 2017.  ISSN 2041-1723.  doi: 10.1038/ncomms15081.  URL http://dx.doi.org/10.1038/ncomms15081.
  * Ji et al. [2021] Yuge Ji, Mohammad Lotfollahi, F. Alexander Wolf, and Fabian J. Theis.  Machine learning for perturbational single-cell omics.  _Cell Systems_ , 12(6):522–537, 2021.  ISSN 2405-4712.  doi: https://doi.org/10.1016/j.cels.2021.05.016.  URL https://www.sciencedirect.com/science/article/pii/S2405471221002027.
  * Lowe et al. [2019] Elijah K. Lowe, Claudia Cuomo, Danila Voronov, and Maria I. Arnone.  _Using ATAC-seq and RNA-seq to increase resolution in GRN connectivity_ , page 115–126.  Elsevier, 2019.  doi: 10.1016/bs.mcb.2018.11.001.  URL http://dx.doi.org/10.1016/bs.mcb.2018.11.001.
  * Baysoy et al. [2023] Alev Baysoy, Zhiliang Bai, Rahul Satija, and Rong Fan.  The technological landscape and applications of single-cell multi-omics.  _Nature Reviews Molecular Cell Biology_ , 24(10):695–713, 2023.  ISSN 1471-0080.  doi: https://doi.org/10.1038/s41580-023-00615-w.  URL https://doi.org/10.1038/s41580-023-00615-w.
  * Grandi et al. [2022] Fiorella C Grandi, Hailey Modi, Lucas Kampman, and M Ryan Corces.  Chromatin accessibility profiling by atac-seq.  _Nature protocols_ , 17(6):1518–1552, 2022.  ISSN 1750-2799.  doi: https://doi.org/10.1038/s41596-022-00692-9.  URL https://doi.org/10.1038/s41596-022-00692-9.
  * Lopez et al. [2018] Romain Lopez, Jeffrey Regier, Michael B. Cole, Michael I. Jordan, and Nir Yosef.  Deep generative modeling for single-cell transcriptomics.  _Nature Methods_ , 15(12):1053–1058, November 2018.  ISSN 1548-7105.  doi: 10.1038/s41592-018-0229-2.  URL http://dx.doi.org/10.1038/s41592-018-0229-2.
  * Palma et al. [2023] Alessandro Palma, Sergei Rybakov, Leon Hetzel, and Fabian J Theis.  Modelling single-cell rna-seq trajectories on a flat statistical manifold.  In _NeurIPS 2023 AI for Science Workshop_ , 2023.
  * Lotfollahi et al. [2019] Mohammad Lotfollahi, F Alexander Wolf, and Fabian J Theis.  scgen predicts single-cell perturbation responses.  _Nature methods_ , 16(8):715–721, 2019.  doi: 10.1038/s41592-019-0494-8.  URL https://doi.org/10.1038/s41592-019-0494-8.
  * Lotfollahi et al. [2023] Mohammad Lotfollahi, Anna Klimovskaia Susmelj, Carlo De Donno, Leon Hetzel, Yuge Ji, Ignacio L Ibarra, Sanjay R Srivatsan, Mohsen Naghipourfar, Riza M Daza, Beth Martin, Jay Shendure, Jose L McFaline-Figueroa, Pierre Boyeau, F Alexander Wolf, Nafissa Yakubova, Stephan Günnemann, Cole Trapnell, David Lopez-Paz, and Fabian J Theis.  Predicting cellular responses to complex perturbations in high-throughput screens.  _Molecular Systems Biology_ , 19(6), May 2023.  ISSN 1744-4292.  doi: 10.15252/msb.202211517.  URL http://dx.doi.org/10.15252/msb.202211517.
  * Hetzel et al. [2022] Leon Hetzel, Simon Boehm, Niki Kilbertus, Stephan Günnemann, Mohammad Lotfollahi, and Fabian Theis.  Predicting cellular responses to novel drug perturbations at a single-cell resolution.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 26711–26722. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/aa933b5abc1be30baece1d230ec575a7-Paper-Conference.pdf.
  * Gayoso et al. [2023] Adam Gayoso, Philipp Weiler, Mohammad Lotfollahi, Dominik Klein, Justin Hong, Aaron Streets, Fabian J. Theis, and Nir Yosef.  Deep generative modeling of transcriptional dynamics for rna velocity analysis in single cells.  _Nature Methods_ , 21(1):50–59, September 2023.  ISSN 1548-7105.  doi: 10.1038/s41592-023-01994-w.  URL http://dx.doi.org/10.1038/s41592-023-01994-w.
  * Chen et al. [2022] Zhanlin Chen, William C. King, Aheyon Hwang, Mark Gerstein, and Jing Zhang.  Deepvelo: Single-cell transcriptomic deep velocity field learning with neural ordinary differential equations.  _Science Advances_ , 8(48), December 2022.  ISSN 2375-2548.  doi: 10.1126/sciadv.abq3745.  URL http://dx.doi.org/10.1126/sciadv.abq3745.
  * Farrell et al. [2023] Spencer Farrell, Madhav Mani, and Sidhartha Goyal.  Inferring single-cell transcriptomic dynamics with structured latent gene expression dynamics.  _Cell Reports Methods_ , 3(9):100581, September 2023.  ISSN 2667-2375.  doi: 10.1016/j.crmeth.2023.100581.  URL http://dx.doi.org/10.1016/j.crmeth.2023.100581.
  * Luo et al. [2024] Erpai Luo, Minsheng Hao, Lei Wei, and Xuegong Zhang.  scdiffusion: conditional generation of high-quality single-cell data using diffusion model.  _arXiv preprint arXiv:2401.03968_ , 2024.
  * Marouf et al. [2020] Mohamed Marouf, Pierre Machart, Vikas Bansal, Christoph Kilian, Daniel S Magruder, Christian F Krebs, and Stefan Bonn.  Realistic in silico generation and augmentation of single-cell rna-seq data using generative adversarial networks.  _Nature communications_ , 11(1):166, 2020.  ISSN 2041-1723.  doi: 10.1038/s41467-019-14018-z.  URL https://doi.org/10.1038/s41467-019-14018-z.
  * Yang et al. [2023a] Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang.  Diffusion models: A comprehensive survey of methods and applications.  _ACM Comput. Surv._ , 56(4), nov 2023a.  ISSN 0360-0300.  doi: 10.1145/3626235.  URL https://doi.org/10.1145/3626235.
  * Michael Love [2017] Simon Anders Michael Love.  Deseq2, 2017.  URL https://bioconductor.org/packages/DESeq2.
  * Hafemeister and Satija [2019] Christoph Hafemeister and Rahul Satija.  Normalization and variance stabilization of single-cell rna-seq data using regularized negative binomial regression.  _Genome Biology_ , 20(1), December 2019.  ISSN 1474-760X.  doi: 10.1186/s13059-019-1874-1.  URL http://dx.doi.org/10.1186/s13059-019-1874-1.
  * Albergo and Vanden-Eijnden [2022] Michael S Albergo and Eric Vanden-Eijnden.  Building normalizing flows with stochastic interpolants.  _arXiv preprint arXiv:2209.15571_ , 2022.
  * Liu et al. [2022a] Xingchao Liu, Chengyue Gong, and Qiang Liu.  Flow straight and fast: Learning to generate and transfer data with rectified flow.  _arXiv preprint arXiv:2209.03003_ , 2022a.
  * Lipman et al. [2023] Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le.  Flow matching for generative modeling.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=PqvMRDCJT9t.
  * Dao et al. [2023] Quan Dao, Hao Phung, Binh Nguyen, and Anh Tran.  Flow matching in latent space.  _arXiv preprint arXiv:2307.08698_ , 2023.
  * Tong et al. [2024] Alexander Tong, Kilian FATRAS, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector-Brooks, Guy Wolf, and Yoshua Bengio.  Improving and generalizing flow-based generative models with minibatch optimal transport.  _Transactions on Machine Learning Research_ , 2024.  ISSN 2835-8856.  URL https://openreview.net/forum?id=CD9Snc73AW.  Expert Certification.
  * Zappia et al. [2017] Luke Zappia, Belinda Phipson, and Alicia Oshlack.  Splatter: simulation of single-cell rna sequencing data.  _Genome biology_ , 18(1):174, 2017.  ISSN https://doi.org/10.1186/s13059-017-1305-0.  doi: https://doi.org/10.1186/s13059-017-1305-0.  URL https://doi.org/10.1186/s13059-017-1305-0.
  * Li and Li [2019] Wei Vivian Li and Jingyi Jessica Li.  A statistical simulator scDesign for rational scRNA-seq experimental design.  _Bioinformatics_ , 35(14):i41–i50, 07 2019.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btz321.  URL https://doi.org/10.1093/bioinformatics/btz321.
  * Song et al. [2024] Dongyuan Song, Qingyang Wang, Guanao Yan, Tianyang Liu, Tianyi Sun, and Jingyi Jessica Li.  scdesign3 generates realistic in silico data for multimodal single-cell and spatial omics.  _Nature Biotechnology_ , 42(2):247–252, 2024.  ISSN 1546-1696.  doi: https://doi.org/10.1038/s41587-023-01772-1.  URL https://doi.org/10.1038/s41587-023-01772-1.
  * Gayoso et al. [2021] Adam Gayoso, Zoë Steier, Romain Lopez, Jeffrey Regier, Kristopher L Nazor, Aaron Streets, and Nir Yosef.  Joint probabilistic modeling of single-cell multi-omic data with totalvi.  _Nature methods_ , 18(3):272–282, 2021.  ISSN 1548-7105.  doi: https://doi.org/10.1038/s41592-020-01050-x.  URL https://doi.org/10.1038/s41592-020-01050-x.
  * Levine et al. [2023] Daniel Levine, Sacha Lévy, Syed Asad Rizvi, Nazreen Pallikkavaliyaveetil, Xingyu Chen, David Zhang, Sina Ghadermarzi, Ruiming Wu, Zihe Zheng, Ivan Vrkic, Anna Zhong, Daphne Raskin, Insu Han, Antonio Henrique de Oliveira Fonseca, Josue Ortega Caro, Amin Karbasi, Rahul M. Dhodapkar, and David van Dijk.  Cell2sentence: Teaching large language models the language of biology.  September 2023.  doi: 10.1101/2023.09.11.557287.  URL http://dx.doi.org/10.1101/2023.09.11.557287.
  * Eyring et al. [2023] Luca Eyring, Dominik Klein, Théo Uscidda, Giovanni Palla, Niki Kilbertus, Zeynep Akata, and Fabian Theis.  Unbalancedness in neural monge maps improves unpaired domain translation.  _arXiv preprint arXiv:2311.15100_ , 2023.
  * Pooladian et al. [2023] Aram-Alexandre Pooladian, Heli Ben-Hamu, Carles Domingo-Enrich, Brandon Amos, Yaron Lipman, and Ricky TQ Chen.  Multisample flow matching: Straightening flows with minibatch couplings.  _arXiv preprint arXiv:2304.14772_ , 2023.
  * Jing et al. [2024] Bowen Jing, Bonnie Berger, and Tommi Jaakkola.  Alphafold meets flow matching for generating protein ensembles.  _arXiv preprint arXiv:2402.04845_ , 2024.
  * Yim et al. [2023] Jason Yim, Andrew Campbell, Andrew YK Foong, Michael Gastegger, José Jiménez-Luna, Sarah Lewis, Victor Garcia Satorras, Bastiaan S Veeling, Regina Barzilay, Tommi Jaakkola, et al.  Fast protein backbone generation with se (3) flow matching.  _arXiv preprint arXiv:2310.05297_ , 2023.
  * Chen and Lipman [2023] Ricky TQ Chen and Yaron Lipman.  Riemannian flow matching on general geometries.  _arXiv preprint arXiv:2302.03660_ , 2023.
  * Kapusniak et al. [2024] Kacper Kapusniak, Peter Potaptchik, Teodora Reu, Leo Zhang, Alexander Tong, Michael Bronstein, Avishek Joey Bose, and Francesco Di Giovanni.  Metric flow matching for smooth interpolations on the data manifold, 2024.
  * Zheng et al. [2023] Qinqing Zheng, Matt Le, Neta Shaul, Yaron Lipman, Aditya Grover, and Ricky TQ Chen.  Guided flows for generative modeling and decision making.  _arXiv preprint arXiv:2311.13443_ , 2023.
  * Klein et al. [2023] Dominik Klein, Théo Uscidda, Fabian Theis, and Marco Cuturi.  Entropic (gromov) wasserstein flow matching with genot, 2023.  URL https://arxiv.org/abs/2310.09254.
  * Chen et al. [2018] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper_files/paper/2018/file/69386f6bb1dfed68692a24c8686939b9-Paper.pdf.
  * Ho and Salimans [2022] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  _arXiv preprint arXiv:2207.12598_ , 2022.
  * Song et al. [2021] Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations, 2021.  URL https://arxiv.org/abs/2011.13456.
  * Yang et al. [2023b] Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang.  Diffusion models: A comprehensive survey of methods and applications.  _ACM Computing Surveys_ , 56(4):1–39, 2023b.
  * Liu et al. [2022b] Nan Liu, Shuang Li, Yilun Du, Antonio Torralba, and Joshua B Tenenbaum.  Compositional visual generation with composable diffusion models.  In _European Conference on Computer Vision_ , pages 423–439. Springer, 2022b.
  * Rombach et al. [2022] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 10684–10695, June 2022.
  * Goodfellow et al. [2020] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  _Commun. ACM_ , 63(11):139–144, oct 2020.  ISSN 0001-0782.  doi: 10.1145/3422622.  URL https://doi.org/10.1145/3422622.
  * La Manno et al. [2018] Gioele La Manno, Ruslan Soldatov, Amit Zeisel, Emelie Braun, Hannah Hochgerner, Viktor Petukhov, Katja Lidschreiber, Maria E Kastriti, Peter Lönnerberg, Alessandro Furlan, et al.  Rna velocity of single cells.  _Nature_ , 560(7719):494–498, 2018.  ISSN 1476-4687.  doi: https://doi.org/10.1038/s41586-018-0414-6.  URL https://doi.org/10.1038/s41586-018-0414-6.
  * tab [2018] Single-cell transcriptomics of 20 mouse organs creates a tabula muris.  _Nature_ , 562(7727):367–372, October 2018.  ISSN 1476-4687.  doi: 10.1038/s41586-018-0590-4.  URL http://dx.doi.org/10.1038/s41586-018-0590-4.
  * Sikkema et al. [2023] Lisa Sikkema, Ciro Ramírez-Suástegui, Daniel C Strobl, Tessa E Gillett, Luke Zappia, Elo Madissoon, Nikolay S Markov, Laure-Emmanuelle Zaragosi, Yuge Ji, Meshal Ansari, et al.  An integrated cell atlas of the lung in health and disease.  _Nature Medicine_ , 29(6):1563–1577, 2023.  ISSN 1546-170X.  doi: https://doi.org/10.1038/s41591-023-02327-2.  URL https://doi.org/10.1038/s41591-023-02327-2.
  * Borgwardt et al. [2006] Karsten M. Borgwardt, Arthur Gretton, Malte J. Rasch, Hans-Peter Kriegel, Bernhard Schölkopf, and Alex J. Smola.  Integrating structured biological data by Kernel Maximum Mean Discrepancy.  _Bioinformatics_ , 22(14):e49–e57, 07 2006.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btl242.  URL https://doi.org/10.1093/bioinformatics/btl242.
  * Zhou et al. [2011] Yi-Hui Zhou, Kai Xia, and Fred A. Wright.  A powerful and flexible approach to the analysis of rna sequence count data.  _Bioinformatics_ , 27(19):2672–2678, August 2011.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btr449.  URL http://dx.doi.org/10.1093/bioinformatics/btr449.
  * Andrews and Hemberg [2018] Tallulah S Andrews and Martin Hemberg.  M3drop: dropout-based feature selection for scrnaseq.  _Bioinformatics_ , 35(16):2865–2867, December 2018.  ISSN 1367-4811.  doi: 10.1093/bioinformatics/bty1044.  URL http://dx.doi.org/10.1093/bioinformatics/bty1044.
  * Ashuach et al. [2023] Tal Ashuach, Mariano I Gabitto, Rohan V Koodli, Giuseppe-Antonio Saldi, Michael I Jordan, and Nir Yosef.  Multivi: deep generative model for the integration of multimodal data.  _Nature Methods_ , 20(8):1222–1231, 2023.  ISSN 1548-7105.  doi: https://doi.org/10.1038/s41592-023-01909-9.  URL https://doi.org/10.1038/s41592-023-01909-9.
  * Ashuach et al. [2022] Tal Ashuach, Daniel A. Reidenbach, Adam Gayoso, and Nir Yosef.  Peakvi: A deep generative model for single-cell chromatin accessibility analysis.  _Cell Reports Methods_ , 2(3):100182, 2022.  ISSN 2667-2375.  doi: https://doi.org/10.1016/j.crmeth.2022.100182.  URL https://www.sciencedirect.com/science/article/pii/S2667237522000376.
  * Aizawa [2003] Akiko Aizawa.  An information-theoretic perspective of tf–idf measures.  _Information Processing and; Management_ , 39(1):45–65, January 2003.  ISSN 0306-4573.  doi: 10.1016/s0306-4573(02)00021-3.  URL http://dx.doi.org/10.1016/S0306-4573(02)00021-3.
  * Luecken et al. [2021] Malte Luecken, Daniel Burkhardt, Robrecht Cannoodt, Christopher Lance, Aditi Agrawal, Hananeh Aliee, Ann Chen, Louise Deconinck, Angela Detweiler, Alejandro Granados, Shelly Huynh, Laura Isacco, Yang Kim, Dominik Klein, BONY DE KUMAR, Sunil Kuppasani, Heiko Lickert, Aaron McGeever, Joaquin Melgarejo, Honey Mekonen, Maurizio Morri, Michaela Müller, Norma Neff, Sheryl Paul, Bastian Rieck, Kaylie Schneider, Scott Steelman, Michael Sterr, Daniel Treacy, Alexander Tong, Alexandra-Chloe Villani, Guilin Wang, Jia Yan, Ce Zhang, Angela Pisco, Smita Krishnaswamy, Fabian Theis, and Jonathan M Bloom.  A sandbox for prediction and integration of dna, rna, and proteins in single cells.  In J. Vanschoren and S. Yeung, editors, _Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks_ , volume 1, 2021.  URL https://datasets-benchmarks-proceedings.neurips.cc/paper_files/paper/2021/file/158f3069a435b314a80bdcb024f8e422-Paper-round2.pdf.
  * Richter et al. [2024] Till Richter, Mojtaba Bahrami, Yufan Xia, David S Fischer, and Fabian J Theis.  Delineating the effective use of self-supervised learning in single-cell genomics.  _bioRxiv_ , 2024.
  * Cui et al. [2024] Haotian Cui, Chloe Wang, Hassaan Maan, Kuan Pang, Fengning Luo, Nan Duan, and Bo Wang.  scgpt: toward building a foundation model for single-cell multi-omics using generative ai.  _Nature Methods_ , pages 1–11, 2024.  ISSN 1548-7105.  doi: https://doi.org/10.1038/s41592-024-02201-0.  URL https://doi.org/10.1038/s41592-024-02201-0.
  * Yoshida et al. [2022] Masahiro Yoshida, Kaylee B Worlock, Ni Huang, Rik GH Lindeboom, Colin R Butler, Natsuhiko Kumasaka, Cecilia Dominguez Conde, Lira Mamanova, Liam Bolt, Laura Richardson, et al.  Local and systemic responses to sars-cov-2 infection in children and adults.  _Nature_ , 602(7896):321–327, 2022.  ISSN 1476-4687.  doi: https://doi.org/10.1038/s41586-021-04345-x.  URL https://doi.org/10.1038/s41586-021-04345-x.
  * Kingma et al. [2021] Diederik Kingma, Tim Salimans, Ben Poole, and Jonathan Ho.  Variational diffusion models.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 21696–21707. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/b578f2a52a0229873fefc2a4b06377fa-Paper.pdf.
  * Paszke et al. [2017] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.  doi: 10.1109/CVPR.2016.90.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper_files/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * [66] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, Jinkyoo Park, and Stefano Ermon.  Torchdyn: Implicit models and neural numerical methods in pytorch.
  * Van Rossum and Drake [2009] Guido Van Rossum and Fred L. Drake.  _Python 3 Reference Manual_.  CreateSpace, Scotts Valley, CA, 2009.  ISBN 1441412697.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.  URL https://proceedings.neurips.cc/paper_files/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf.
  * Wolf et al. [2018] F Alexander Wolf, Philipp Angerer, and Fabian J Theis.  Scanpy: large-scale single-cell gene expression data analysis.  _Genome biology_ , 19:1–5, 2018.  ISSN 1474-760X.  doi: https://doi.org/10.1186/s13059-017-1382-0.  URL https://doi.org/10.1186/s13059-017-1382-0.
  * Gretton et al. [2012] Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(25):723–773, 2012.  URL http://jmlr.org/papers/v13/gretton12a.html.
  * Flamary et al. [2021] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.  URL http://jmlr.org/papers/v22/20-451.html.
  * Bredikhin et al. [2022] Danila Bredikhin, Ilia Kats, and Oliver Stegle.  Muon: multimodal omics analysis framework.  _Genome Biology_ , 23(1):42, 2022.

##  Appendix A Datasets and code

We make the code for CFGen available at https://github.com/theislab/CFGen. All
the datasets are open source and made available by the respective
publications.

##  Appendix B Theoretical supplement

###  B.1 Poisson-gamma and negative binomial distribution

A possible parameterization of the negative binomial distribution is via a
mean $\mu$ and an inverse dispersion parameter $\theta$ following the
Probability Mass Function (PMF):

$p_{\textrm{NB}}(x\mid\mu,\theta)=\frac{\gamma(\theta+x)}{x!\gamma(\theta)}\Big{(}\frac{\theta}{\theta+\mu}\Big{)}^{\theta}\Big{(}\frac{\mu}{\theta+\mu}\Big{)}^{x}\>.$
(11)

One can show that the negative binomial distribution is obtained by a
continuous mixture of Poisson distributions with gamma-distributed rate. More
formally, define a Poisson model $x\sim\mathrm{Poisson(\lambda)}$ with
$\lambda\geq 0$. The parameter $\lambda$ represents both the mean and the
variance of the distribution. Since the mean and the variance of the
distribution are equal, a Poisson model is not suited for modeling over-
dispersed counts (i.e., where the variance exceeds the mean). A way around
such a shortcoming is to model the rate of a Poisson distribution as a random
variable following a gamma distribution.

$\displaystyle x\sim\mathrm{Poisson}(\lambda)\>,$ (12)
$\displaystyle\lambda\sim\mathrm{gamma}\left(\theta,\frac{\mu}{\theta}\right)\>,$
(13)

where $\theta\geq 0$ is the shape parameter and $\mu\geq 0$ the scale
parameter. Marginalizing out $\lambda$ in the PMF of the Poisson distribution,
one retrieves the PMF of a negative binomial with mean $\mu$ and inverse
dispersion $\theta$. Notably, the variance of such a negative binomial
parameterization is $\mu+\frac{\mu^{2}}{\theta}$. As one can see, the variance
always exceeds the mean as long as $\theta$ is finite, making the negative
binomial distribution a suitable tool to model over-dispersed counts.

###  B.2 Flow Matching with Gaussian paths

Flow Matching [26] learns a time-dependent vector field $u_{t}(\mathbf{z})$,
with $t\in[0,1]$, generating the probability path $p_{t}(\mathbf{z})$, such
that $p_{0}=\mathcal{N}(\mathbf{0},\mathbf{I})$ is a standard Gaussian prior
and $p_{1}$ is a complex distribution. A common way to formulate the marginal
$p_{t}$ is via the mixture

$p_{t}(\mathbf{z})=\int
p_{t}(\mathbf{z}|\mathbf{z}_{1})q(\mathbf{z}_{1})\mathrm{d}\mathbf{z}_{1},$

where $q$ is the target data distribution. The following marginal vector field
generates such a mixture of paths [26]:

$u_{t}(\mathbf{z})=\int
u_{t}(\mathbf{z}|\mathbf{z}_{1})\frac{p(\mathbf{z}|\mathbf{z}_{1})q(\mathbf{z}_{1})}{p_{t}(\mathbf{z})}\textrm{d}\mathbf{z}_{1}\;.$

While $u_{t}$ is intractable, the conditional field
$u_{t}(\mathbf{z}|\mathbf{z}_{1})$ has a closed-form expression given an
observed data point $\mathbf{z}_{1}$ and a pre-defined choice of the
probability path $p_{t}(\mathbf{z}|\mathbf{z}_{1})$ satisfying the boundary
conditions $p_{0}(\cdot|\mathbf{z}_{1})=p_{0}$ and
$p_{1}(\cdot|\mathbf{z}_{1})=\delta(\mathbf{z}-\mathbf{z}_{1})$. Notably,
$u_{t}(\mathbf{z}|\mathbf{z}_{1})$ admits the same minimizer as
$u_{t}(\mathbf{z})$ and can be used as a regression target during training.
Following Lipman et al. [26], one can assume Gaussian probability paths
$p_{t}(\mathbf{z}|\mathbf{z}_{1})=\mathcal{N}(\mathbf{z}|\alpha_{t}\mathbf{z}_{1},\sigma_{t}^{2}\textbf{I})$,
where the tuple $(\alpha_{t},\sigma_{t})$ is called scheduler and satisfies
$\alpha_{0}=0=\sigma_{1}$ and $\alpha_{1}=1=\sigma_{0}$. In this work, we use
standard linear scheduling, where $\alpha=t$ and $\sigma=1-t$.

###  B.3 The relationship between Flow Matching and classifier-free guidance

Zheng et al. [40] draw a relationship between classifier-free guidance in
score-based models [43] and the Flow Matching vector field $u_{t}$.
Specifically, the authors show that the following relationship between the
score $\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z}|y)$ and the marginal vector
field $u_{t}(\mathbf{z}|y)$ holds:

$u_{t}(\mathbf{z}|y)=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y)\>,$ (14)

both in the conditional case and when $y=\varnothing$, with
$a_{t}=\frac{\dot{\alpha}_{t}}{\alpha_{t}}$ and
$b_{t}=(\dot{\alpha}_{t}\sigma_{t}-\alpha_{t}\dot{\sigma}_{t})\frac{\sigma_{t}}{\alpha_{t}}$.

Let the equation

$\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|y)=(1-\omega)\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\omega\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z}|y)$

be the classifier-free guidance score as formulated by Ho and Salimans [43]
with guidance strength $\omega$. Zheng et al. [40] define the vector field of
classifier-free Flow Matching as

$\tilde{u}_{t}(\mathbf{z}|y)=(1-\omega)u_{t}(\mathbf{z})+\omega\,u_{t}(\mathbf{z}|y)$
(15)

which is related to the classifier-free guidance score by

$\tilde{u}_{t}(\mathbf{z}|y)=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|y)\>,$

derived by substituting Eq. 14 into Eq. 15.

###  B.4 Proof of 1

###### Proposition 1

If the attributes $y_{1},...,y_{K}$ are conditionally independent given
$\mathbf{z}$, the vector field

$\tilde{u}_{t}(\mathbf{z}|\mathbf{y})=u_{t}(\mathbf{z})+\sum_{i=1}^{K}\omega_{i}[u_{t}(\mathbf{z}|y_{i})-u_{t}(\mathbf{z})]$

coincides with the velocity of the probability-flow ODE associated with the
generative SDE of a diffusion model with the compositional score as in Eq. 9.

Proof. (1) We first justify the conditional independence assumption and
successfully show the described equality.

#### Conditional independence assumption.

Given a variable $\mathbf{z}$ and a set of attributes
$\mathbf{y}=y_{1},...,y_{K}$, our aim is to sample conditionally from the
marginal distribution $p_{t}(\mathbf{z}|y_{1},...,y_{K})$. To obtain the
compositional score formulation in Eq. 9, one must first assume that the
attributes are conditionally independent given $\mathbf{z}$. Then:

$p_{t}(\mathbf{z}|y_{1},...,y_{K})\propto
p_{t}(\mathbf{z},y_{1},...,y_{K})=p_{t}(\mathbf{z})\prod_{i=1}^{K}p_{t}(y_{i}|\mathbf{z})\propto
p_{t}(\mathbf{z})\prod_{i=1}^{K}\frac{p_{t}(\mathbf{z}|y_{i})}{p_{t}(\mathbf{z})}\;.$
(16)

Taking the logarithm and then the gradient with respect to $\mathbf{z}$ on
both sides in Eq. 16, we obtain:

$\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{1},...,y_{K})=\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\sum_{i=1}^{K}\left[\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})-\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z})\right]\;.$
(17)

Additionally, if the goal is to sample with attribute-specific guidance
strengths $\\{\omega_{i}\\}_{i=1}^{K}$ according to a modified conditional
distribution

$\tilde{p}(\mathbf{z}|y_{1},...,y_{K})\propto
p_{t}(\mathbf{z})\prod_{i=1}^{K}\left[\frac{p_{t}(\mathbf{z}|y_{i})}{p_{t}(\mathbf{z})}\right]^{\omega_{i}}\;,$
(18)

the score in Eq. 17 becomes:

$\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|y_{1},...,y_{K})=\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\sum_{i=1}^{K}\omega_{i}\left[\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})-\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z})\right]\;.$
(19)

The formulation in Eq. 19 is used to parameterize the drift of the reverse-
time SDE that generates data points conditionally on the attributes
$y_{1},...,y_{K}$ with guidance strengths $\\{\omega_{i}\\}_{i=1}^{K}$.

#### Proof of equality.

Following the standard theory of score-based models [45] and their
compositional version [46], we first note that one can use the compositional
classifier-free guidance score

$\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|\mathbf{y})=\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\sum_{i=1}^{K}\omega_{i}\left[\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})-\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z})\right]\;$

to simulate the probability-flow ODE

$\dot{\mathbf{z}}=f_{t}\mathbf{z}-\frac{1}{2}g_{t}^{2}\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|\mathbf{y})\;.$
(20)

Given a scheduling pair $(\alpha_{t},\sigma_{t})$, Kingma et al. [62] show
that

$f_{t}=\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t},\quad
g_{t}^{2}=\frac{\mathrm{d}\sigma_{t}^{2}}{\mathrm{d}t}-2\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t}\sigma_{t}\>,$
(21)

which yield

$\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t}=\frac{\dot{\alpha}}{\alpha_{t}}=a_{t},\quad-\frac{1}{2}\frac{\mathrm{d}\sigma_{t}^{2}}{\mathrm{d}t}+\frac{\mathrm{d}\log\alpha_{t}}{\mathrm{d}t}\sigma_{t}=(\dot{\alpha}\sigma_{t}-\alpha_{t}\dot{\sigma}_{t})\frac{\sigma_{t}}{\alpha_{t}}=b_{t}\>.$
(22)

Hence, the probability-flow ODE is written as

$\dot{\mathbf{z}}=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|\mathbf{y})\;.$
(23)

Given the score in Eq. 10, we first use results from [40] explained in Section
B.3 to justify the following equalities:

$\displaystyle u_{t}(\mathbf{z})=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})\>$ (24) $\displaystyle
u_{t}(\mathbf{z}|y_{i})=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})\>$ (25)

Plugging Eq. 25 and Eq. 24 into Eq. 10 yields:

$\displaystyle\tilde{u}_{t}(\mathbf{z}|\mathbf{y})$
$\displaystyle=a_{t}\mathbf{z}+b_{t}\left[\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z})+\sum_{i=1}^{K}\omega_{i}(\nabla_{\mathbf{z}}\log
p_{t}(\mathbf{z}|y_{i})-\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z}))\right]$
(26)
$\displaystyle=a_{t}\mathbf{z}+b_{t}\nabla_{\mathbf{z}}\log\tilde{p}_{t}(\mathbf{z}|\mathbf{y})\>.$
(27)

completing the proof.

###  B.5 Guidance algorithm

Algorithm 1 Train CFGen with multiple attributes on scRNA-seq

0: Probability of unconditional generation $p_{\mathrm{uncond}}$, trained
encoder $f_{\eta}$, scheduling $(\alpha_{t},\sigma_{t})$.

1: Initialize $v_{t,\xi}$

2: while not converged do

3: Sample $(\mathbf{x}_{1},y_{1},...,y_{K})$ from the data

4: $\mathbf{z}_{1}\leftarrow f_{\eta}(\mathbf{x}_{1})$

5: Sample $t$ from $\mathcal{U}[0,1]$

6: $l\leftarrow\textrm{Sum of entries of }\mathbf{x}_{1}$

7: Sample $b$ from Bernoulli($p_{\mathrm{uncond}}$)

8: if $b$ = 1 then

9: $y\leftarrow\emptyset$

10: else

11: $y\leftarrow\textrm{random sample among}\>y_{1},...,y_{K}$

12: end if

13: $\mathbf{z}_{0}\sim p(\mathbf{z}_{0})$ {sample noise}

14:
$\mathbf{z}_{t}\leftarrow\alpha_{t}\mathbf{z}_{1}+\sigma_{t}\mathbf{z}_{0}$
{noisy data point}

15:
$\dot{\mathbf{z}}_{t}\leftarrow\dot{\alpha}_{t}\mathbf{z}_{1}+\dot{\sigma}_{t}\mathbf{z}_{0}$

16: Take gradient step on
$\nabla_{\xi}||v_{t,\xi}(\mathbf{z}_{t}|y,l)-\dot{\mathbf{z}}_{t}||^{2}$

17: end while

Output: $v_{t,\xi}$

Algorithm 2 Sampling from multi-attribute guided CFGen for scRNA-seq

0: Trained velocity field $v_{t,\xi}$, conditions $y_{1},...,y_{K}$, guidance
parameters $\omega_{1},...,\omega_{K}$, size factor distribution parameters
$(\mu_{l},\sigma_{l})$, number of ODE steps $n_{\textrm{ode}}$,trained decoder
$h_{\psi}$, trained inverse dispersion parameter $\bm{\theta}$.

1: Sample size factor $l$ from $\mathrm{LogNormal}(\mu_{l},\sigma_{l})$

2: $\mathbf{z}_{0}\sim p(\mathbf{z}_{0})$ {sample noise}

3: $n\leftarrow 1/n_{\textrm{ode}}$ {step size}

4: $\tilde{u}_{t}(\cdot)\leftarrow
v_{t,\xi}(\cdot|\emptyset,l)+\sum_{i=1}^{K}\omega_{i}\left[v_{t,\xi}(\cdot|y_{i},l)-v_{t,\xi}(\cdot|\emptyset,l)\right]$
{guided velocity}

5: for $t=0,n,\ldots,1-n$ do

6: $\mathbf{z}_{t+n}\leftarrow\textrm{ODEStep}(\tilde{u}_{t},\mathbf{z}_{t})$
{ODE solver step}

7: end for

8: $\mathbf{x}_{1}\leftarrow\textrm{Sample from
NB}(l\,\mathrm{softmax}(h_{\psi}(\mathbf{z}_{1})),\bm{\theta})$

Output: $\mathbf{x}_{1}$

##  Appendix C Model details

The CFGen model is implemented in PyTorch [63], version 2.1.2.

###  C.1 The CFGen autoencoder

Before training the flow generating noise from data, we embed the data with an
autoencoder model trained with maximum likelihood optimization.

#### Encoder.

The encoder is a multi-layer perceptron (MLP) with two hidden layers of
dimension [512, 256]. The last layer maps to a latent space of variable size.
In our experiment, we use 50 latent dimensions for all datasets excluding the
Human Lung Cell Atlas (HLCA) and the Tabula Muris datasets, where we set the
latent dimension to 100 for higher expressivity. In the multi-modal setting,
the different modalities are embedded into the same latent space. Both RNA and
ATAC inputs initially pass through a modality-specific MLP encoder. For ATAC,
we fix the dimension of the hidden layers to [1024, 512] due to its higher
dimensionality. The outputs of the two modality-specific encoders are then
concatenated and further encoded by a shared encoder layer, mapping to a
100-dimensional latent space.

#### Decoder.

The decoder maps the latent space to the parameter space of a likelihood
model. When dealing with multi-modal data, one decoder per modality is used to
sample from the noise model corresponding to the specific modality.

  * •

For scRNA-seq, we map the latent representation to the mean of the Negative
Binomial likelihood, $\bm{\mu}$, with one dimension per gene. Following [11],
we compute the softmax transformation of the output across the gene dimension.
The resulting probabilities are multiplied by the size factor, which is the
total number of transcripts per cell. The inverse dispersion is a learned
model parameter implemented via torch.nn.Parameter. Crucially, we offer the
option to model the inverse dispersion per gene and attribute, in case this
reflects the properties of the datasets.

  * •

When the DNA accessibility modality is present, the decoder is used to map the
latent space to continuous logit values. These are processed with a sigmoid
function for each dimension, with no need for a size factor.

#### Additional training details.

We train the encoder and decoder networks jointly via likelihood optimization.
Therefore, the weights of the networks are tweaked (together with the inverse
dispersion parameter) to produce the parameters that maximize the likelihood
of the data under a pre-defined noise model. For scRNA-seq, we employ the
negative binomial distribution. For ATAC, we use Bernoulli likelihood. The
losses from different modalities are summed before applying backpropagation.
Notice that scRNA-seq data are provided by the encoder in their
$\log$-transformed version for training stability. Intuitively, the loss is
evaluated on the original counts.

For all settings, we keep the learning rate to $0.001$, while all couples of
layers are interleaved with 1-dimensional batch normalization layers. We use a
standard AdamW optimizer and the ELU activation function as a non-linearity.

###  C.2 The flow model

#### The flow architecture.

The flow model inputs the latent representation computed by the encoder and
produces a vector field used to simulate paths generating data from noise. The
architecture is essentially a deep dimensionality-preserving ResNet [64]. The
flow architecture is made of the following modules:

  * •

A linear projection layer from the input dimension to the hidden dimension of
the flow model.

  * •

Three middle ResNet blocks stacked on top of each other performing
representation learning and conditioning.

  * •

An output layer with one level of non-linearity, given by a SiLU activation
function.

  * •

A time embedder, inputting sinusoidal multi-dimensional encoding [65] of the
time with 1e4 frequency value. It is represented by an MLP with two layers and
SiLU non-linearity.

  * •

A size factor embedder, inputting a sinusoidal embedding [65] of the size
factor to approximately guide generation to a pre-defined number of
transcripts. This is used only if the size factor is a conditioner for the
model, that is, we do not assume $p(\mathbf{z}|y,l)=p(\mathbf{z}|y)$. Before
being passed through sinusoidal embeddings, the $\log$ size factor is
normalized to a value lying approximately between 0 and 1 using the maximum
and minimum $\log$ size factors from the dataset.

  * •

A covariate embeddings for all different conditioning attributes.

#### Additional technical details.

During training the covariate embeddings are elementwise summed to the time
embedding and the size factor embedding (if applicable). Hence, all embeddings
are either designed to share the same size or are transformed to a common
dimensionality. The sum of such covariate representations is passed to the
ResNet middle blocks as a single vector.

The summed conditioning embedding and the down-projected input are provided to
the ResNet blocks. The ResNet blocks consist of:

  * •

A non-linear input transformation of the state embedding.

  * •

A linear encoder for the covariate embedding.

  * •

A non-linear output transformation.

  * •

A skip connection.

The results of the non-linear input transformation and the covariate encoder
are summed and passed to the output transformation. The result is summed to
the input of the ResNet via the skip connection as in traditional residual
blocks [64]. All non-linear transformations are simple [siLU, Linear] stacks.

In the standard setting, we train the flow for $1,000$ epochs, using AdamW as
an optimizer, a learning rate of $0.001$ and batch size $256$.

###  C.3 Covariate embeddings

Covariate embeddings are trainable torch.nn.Embdding layers of pre-defined
size. In our experiments, we use size $100$ in most of the settings.

###  C.4 Sampling from noise

To generate discrete observations from noise, we first draw a covariate from
the associated categorical distribution with proportions obtained from the
observed data and a size factor from the $\mathrm{LogNormal}$ distribution
(with mean and standard deviation as the Maximum Likelihood Estimates (MLE)
from the whole dataset or conditional on a technical effect covariate). Then,
we sample Gaussian noise. We simulate a latent observation from the real
datasets conditionally by integrating the vector field computed by the neural
network in Section C.2 starting from Gaussian noise and using the dopri5
solver with adjoint sensitivity and 1e-5 tolerance from the torchdyn [66]
package in Python3 [67]. We integrate over the $(0,1)$ time interval. The
generated latent vector is decoded to the parameter space of the data
likelihood (negative binomial for scRNA-seq or Bernoulli for ATAC-seq) and
single cells are sampled from the parameterized noise model.

##  Appendix D Scheduling

We use linear scheduling, following the original formulation advanced by
Lipman et al. [26]. For more details, we refer to Section B.2.

##  Appendix E Baseline description

###  E.1 scVI, MultiVI, PeakVI

scVI [11], MultiVI [55] and PeakVI [56] are all VAE-based generative models
for single-cell discrete data. Following the standard VAE setting, such models
learn a Gaussian latent space which is decoded to the parameters of the
discrete likelihood models describing different single-cell modalities. While
PeakVI and scVI generate single modalities (respectively, ATAC and scRNA-seq
data), MultiVI learns a common latent space between modalities, while sampling
from different discrete decoders.

###  E.2 scGAN [20]

The scGAN model is a Generative Adversarial Network (GAN) [68] tailored for
realistic scRNA-seq data generation. It minimizes the Wasserstein distance
between distributions of real and generated cells, utilizing a generator
network to produce synthetic samples and a critic network for discrimination.
The model employs fully connected layers, incorporates a custom library-size
normalization (LSN) layer for stable training, and extends to conditional
scGAN (cscGAN) for type-specific cell generation. Evaluation involves metrics
like t-SNE and marker gene correlation to gauge generated cell quality.
Notably, scGAN has been explored for conditional generation as well. However,
we found the results conditioning the model on cell type to be way worse than
when not providing the label. In the latter case, the model is trained
conditioned on data-dependent Leiden cluster labels. We name such version
unconditional as it does not exploit real labels but data-driven attributes.

###  E.3 scDiffusion [19]

The scDiffusion model comprises an autoencoder, a diffusion backbone network,
and a conditional classifier. The autoencoder transforms gene expression
profiles into latent space embeddings, the diffusion backbone network learns
the reverse diffusion process, and the conditional classifier guides cell
generation under specific conditions. During training, the autoencoder creates
embeddings from real data, followed by diffusion to generate noisy embeddings
for backbone training, while the classifier predicts labels. During inference,
the diffusion backbone denoises embeddings to produce new ones for gene
expression data generation. The autoencoder addresses high-dimensional data
and a non-Gaussian distribution, while the diffusion backbone network utilizes
fully connected layers and a residual structure. In the diffusion process,
noise is iteratively added to embeddings, and during inference, noise is
iteratively removed to generate new embeddings for final data generation.

##  Appendix F Computational Resources

For the implementation of our model, we utilized Python 3.10 [67] for the deep
learning components. The experiments were executed on a variety of GPU
servers, each possessing unique specifications:

  * •

GPU: 16 Tesla V100 GPUs, each with 32GB of RAM

  * •

GPU: 2 Tesla V100 GPUs, each with 16GB of RAM

  * •

GPU: 8 A100-SXM4 GPUs, each with 40GB of RAM

##  Appendix G Data preprocessing and description

Single cells were pre-processed via the Scanpy [69] software. Count
normalization was applied only to baselines requiring real-valued data. In
such settings, cells were normalized to sum to 1e4 and $\log$-transformed.
Since CFGen, MultiVI, PeakVI and scVI work in discrete spaces, we did not
normalize the data to train them. Additionally, we filter out genes that are
expressed in less than 20 cells in all the datasets. In each dataset, we
allocate 80% of the observations to training and 20% to testing.

Table 3: List of datasets considered in this work with the associated number
of genes, cells and cell types.

Dataset name | Number of cells | Number of genes | Number of cell types
---|---|---|---
PMBC3K | 2,638 | 8,573 | 8
Dentate gyrus | 18,213 | 17,002 | 14
Tabula Muris | 245,389 | 19,734 | 123
HLCA | 584,944 | 27,997 | 50
PBMC10k | 10,025 | 25,604 | 14
NeurIPS | 90,261 | 14,087 | 45
PBMC COVID | 344,820 | 2,000 | 29

##  Appendix H Experiment description and evaluation metrics

###  H.1 Wasserstein-2 distance and MMD

We use the Wasserstein-2 distance and the linear Mean Maximum Discrepancy
(MMD) [70] to measure the overlap between real and generated data. To
implement the former we use the Python Optimal Transport (POT) [71] package.
For the linear MMD metric, we resort to the implementation proposed in
333https://github.com/atong01/conditional-flow-matching.

###  H.2 KNNc

Similarly to [33], we train a 10-NN classifier for cell type on the real test
data and evaluate its F1 performance on the generated set. This indicates
whether generated cell types are proximal to the real ones in the real
dataset. Notice that a low score is not necessarily a symbol of poor
generation performance, as it may be caused by the F1 scores being low on the
real data. Therefore, such a metric should be considered comparatively.

###  H.3 Distribution metrics comparisons

To compute the metrics in Table 1, we use all competing models to generate
three datasets of the same size as the original. In the conditional case,
distribution metrics are computed per cell type, namely, each cell type in the
real and generated datasets is a subset and used for comparison. In the
unconditional case, we compare batches of $5,000$ cells sampled from the whole
data distribution. The MMD and Wasserstein-2 distances are computed in the
Principal Component (PC) projection of the data onto 30 dimensions. This is
done to face the curse of dimensionality, which hinders the reliability of
distances in high dimensions. To make PC embeddings comparable between real
and generated cells, we project generated cells using the PC loadings of the
real cells. Crucially, scDiffusion and scGAN generate normalized data, while
CFGen and scVI produce discrete count data. To make the obtained numbers
comparable, we first normalize the output of CFGen and scVI to ensure that
gene counts sum to 1e4, then we $\log$-transform the results of all generative
models to obtain a better range for the distances. The same is performed on
the real data. After pre-processing, all generative models and the real data
represent the same quantity. All metric measures are reported on the test set.

###  H.4 Variance-mean trend plot and sparsity histograms

We select scDiffusion as a continuous baseline for demonstrative purposes.
Given generated cells across datasets (after pre-processing as explained in
Section H.3), we compute the mean and variance expression across cells per
gene and the frequency of unexpressed genes per cell. These values are
compared to each other. Note that the mean-variance relationship is expected
to be quadratic when considering raw counts (see Section B.1). However, after
the normalization and $\log$-transformation needed to compare with
scDiffusion, the trend is no longer quadratic. Nevertheless, the relationship
between mean and variance gene expression is still worth investigating, as it
should follow the empirical behavior of real data.

###  H.5 Multi-modal evaluation

We generate multi-modal data and perform an unconditional comparison with the
ground truth as described in Section H.3. When dealing with ATAC data, we
normalize both real and generated cells using the TF-IDF algorithm implemented
in the MUON package [72].

To generate Fig. 9, we perform the following steps:

  1. 1.

Aggregate average gene expression and peak accessibility (i.e, the fraction of
accessible regions) per cell type per marker gene/peak as described in
444https://muon-tutorials.readthedocs.io/en/latest/single-cell-rna-
atac/pbmc10k/3-Multimodal-Omics-Data-Integration.html. The result is, for both
real and all generated datasets, a matrix cell_type x marker with, as value,
the average expression or accessibility of such marker in such cell type.

  2. 2.

We correlate the rows of such matrix between real and generated datasets. A
high correlation signifies that the generative model correctly captures the
mean marker expression and accessibility across markers per cell type.

###  H.6 Guidance strength experiments

In Fig. 4, first train CFGen on each dataset following Algorithm 1. Upon
successful training, we show the guidance performance qualitatively by
sampling 500 cells for an array of guidance strength combinations between
attributes, keeping the guidance weight of an attribute fixed while varying
the other as shown in the figure. For the unconditional generation (hence,
with guidance strength equal to 0 for both attributes) we generated as many
cells as there are in the dataset to better show the overlap between real and
synthetic cells. When using guidance we train the guided CFGen model using a
probability $p_{\mathrm{uncond}}=0.2$ (see Algorithm 1).

###  H.7 scGPT [60] generalization performance enhancement

We split the PBMC COVID and HLCA datasets into a training set and a held-out
set. To make the generalization task more challenging, we leave out all cells
from 20% of the donors in both datasets. This makes 80 training and 27 test
donors for HLCA and 60 training and 15 test donors for PBMC COVID. After
augmenting the training set (see Section 5.4) we use a pre-trained scGPT model
to embed training and validation sets, for both original and augmented data.
Then we fit a cell-type kNN classifier on the training embeddings and evaluate
it on the held-out set. The results displayed in Fig. 5 show to what extent
the held-out classification performance on a cell type varies as a function of
its frequency in the dataset. A performance improvement indicates that a cell
type better separates from the rest in the LLM representation space after the
model sees additional synthetic examples.

##  Appendix I Additional results

###  I.1 Example of synthetic generation by CFGen

Figure 6: Uni-modal generation of scRNA-seq by CFGen on the PBMC3K dataset.
Real and generated cells are embedded together and visualized as 2D UMAP
coordinates. Figure 7: Uni-modal generation of scRNA-seq by CFGen on the
Dentate gyrus dataset. Real and generated cells are embedded together and
visualized as 2D UMAP coordinates.

Figure 8: Uni-modal generation of scRNA-seq by CFGen on the Tabula Muris
dataset. Real and generated cells are embedded together and visualized as 2D
UMAP coordinates. Figure 9: Uni-modal generation of scRNA-seq by CFGen on the
HLCA dataset. Real and generated cells are embedded together and visualized as
2D UMAP coordinates.

###  I.2 Comparison between CFGen, scVI and scDiffusion on the HLCA and Tabula
Muris datasets

Figure 10: Qualitative comparison of the generation results of CFGen, scVI and
scDiffusion on the HLCA and Tabula Muris datasets. Comparison is performed by
evaluating the similarity of the generated results to real cells. Real and
generated cells for all models are embedded together and visualized as 2D UMAP
coordinates.

###  I.3 Marker gene expression and peak accessibility fraction in real and
generated datasets

Figure 11: Average marker expression per cell type in real and generated data
on the PBMC10k dataset. x-axis \- marker genes. y-axis \- cell types. Figure
12: Average number of cells with accessible marker peaks per cell type in real
and generated data on the PBMC10k dataset. x-axis \- marker peaks. y-axis \-
cell types.