# Modelling the effects of biological intervention in a dynamical gene network

Nicolas Champagnat, Rodolphe Loubaton, Laurent Vallat, Pierre Vallois

► To cite this version:

Nicolas Champagnat, Rodolphe Loubaton, Laurent Vallat, Pierre Vallois. Modelling the effects of biological intervention in a dynamical gene network. 2025. hal-05058235

HAL Id: hal-05058235

https://inria.hal.science/hal-05058235v1

Preprint submitted on 6 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling the effects of biological intervention in a dynamical gene network

Nicolas Champagnat<fnref n="1" />,*, Rodolphe Loubaton<fnref n="2" />,<fnref n="3" />,
Laurent Vallat<fnref n="4" />,<fnref n="5" />, Pierre Vallois<fnref n="1" />

May 6, 2025

## Abstract

Cellular response to environmental and internal signals can be modeled by dynamical gene regulatory networks (GRN). In the literature, three main classes of gene network models can be distinguished: (i) non-quantitative (or data-based) models which do not describe the probability distribution of gene expressions; (ii) quantitative models which fully describe the probability distribution of all genes co-expression; and (iii) mechanistic models which allow for a causal interpretation of gene interactions. We propose two rigorous frameworks to model gene alteration in a dynamical GRN, depending on whether the network model is quantitative or mechanistic. We explain how these models can be used for design of experiment, or, if additional alteration data are available, for validation purposes or to improve the parameter estimation of the original model. We apply these methods to the Gaussian graphical model, which is quantitative but non-mechanistic, and to mechanistic models of Bayesian networks and penalized linear regression.

*Keywords and phrases:* Dynamical gene regulatory networks, Bayesian network, Gaussian graphical model, penalized linear regression, inference, gene knockdown experiments.

*Mathematics Subject Classification:* primary: 62F10, 92C42; secondary: 62F03, 62J07, 62P10.

<footnote n="1">¹Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France
* Corresponding author: Nicolas.Champagnat@inria.fr</footnote>

<footnote n="2">²VetAgro Sup, département Territoires et Société, CS 82212, F-63370 Lempdes, France</footnote>

<footnote n="3">³INRAe, équipe CARAIBE, UMR 1213, F-63122 St Genes Champanelle, France</footnote>

<footnote n="4">⁴University of Strasbourg, CNRS, UMR-7242 Biotechnology and Cell Signaling, F-67400 Illkirch, France</footnote>

<footnote n="5">⁵Department of Molecular Genetic of Cancers, Strasbourg University Hospital, F-67200 Strasbourg, France</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of a Gaussian graphical model and prediction of alteration (knock-out) of gene 1. Each node corresponds to a gene labeled from 1 to 4. (a) Without alteration: the mean expression of gene *i* is $\mu_i$ and the covariance matrix of genes expression is $(\sigma_{ij})_{1\le i,j\le 4}$. (b) After knock-out of gene 1: $\mu_1$ and $\sigma_{11}$ were set to 0; the question of prediction amounts to determine the new values of the means and covariances $(\mu'_i)_{2\le i\le 4}$ and $(\sigma'_{ij})_{2\le i,j\le 4}$.</figcaption></figure>

# 1 Introduction

The functioning of eukaryotic cells in their basal state, and their response to changes in their environment, are dynamic phenomena involving multiple players activated in a cascade with numerous levels of regulation. The sequential activation of multiple genes (transcriptional response) following an extra-cellular signal is a crucial step in this response. Some of these genes are transcriptional factors which activate other genes in cascade, notably genes encoding proteins which are the effectors of the cellular response adapted to the initial signal.

These various interactions between genes can be represented by gene interaction networks (GRNs) [Barabasi and Oltvai, 2004, Filkov, 2005, Hecker et al., 2009, Natale et al., 2018, Huynh-Thu and Sanguinetti, 2019].

Recent developments in high-throughput sequencing techniques (RNA-seq) have made it possible to identify and quantify the various players (messenger RNAs) in these interaction networks, at different times after cell activation [Schleiss et al., 2021]. These sequencing data (dataset) can then be used to infer a mathematical model of the underlying dynamic gene interaction network.

Gene regulatory network reconstruction from gene expression data is a difficult and ill-posed problem [Ud-Dean and Gunawan, 2014]. A large variety of methods were proposed, which may give different interpretations to the edges of the inferred network [Hecker et al., 2009, Wang and Huang, 2014, Natale et al., 2018, Huynh-Thu and Sanguinetti, 2019]. These methods can be classified into three main categories: (i) *data-driven methods* which do not model gene expression and only represent interactions be-
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation of a mechanistic model and prediction of alteration (knock-out) of gene 1. Each node corresponds to a gene labeled from 1 to 4. (a) Without alteration: arrows represent direct influence of a gene on another; they encode a quantitative relation (not shown here). (b) After knock-out of gene 1: the problem of prediction consists in determining how the arrows and the quantitative relations they encode are modified.</figcaption></figure>

tween genes from statistical quantities such as correlation or mutual in-
formation; (ii) *quantitative methods* which specify a parametric family of
distributions for the genes expressions; (iii) *mechanistic methods* which are
particular cases of quantitative methods where the model quantitatively rep-
resents the effects of interactions between small subsets of genes (typically
pairs of genes). Mechanistic models may have different levels of descrip-
tion of biological interactions ranging from phenomenological, often indirect
interactions to chemical direct interactions based on the biological under-
standing of gene regulation (causal networks [Pearl, 2016]). In this work,
we focus on quantitative (mechanistic or non-mechanistic) methods. The
Gaussian graphical model and its variants are typical examples of quanti-
tative but non-mechanistic models [Chiquet et al., 2019]. In the Gaussian
graphical model, the edges in the gene network are non-directed and the
model is characterized by the vector of mean expression ($\mu_i$)$_{1\le i\le N}$ of N
genes and their covariance matrix ($\sigma_{ij}$)$_{1\le i,j\le N}$, as represented in Fig. 1(a).
Mechanistic models can range from simple models such as models based
on penalized linear regressions [Vallat et al., 2013, Omranian et al., 2016,
Ishikawa et al., 2023] to complex models, e.g. based on ordinary differential
equations [Goutsias and Lee, 2007] or piecewise-deterministic models [Bonnaffoux et al.,
2019, Ventre et al., 2023]. Other standard classes of mechanistic models
used for gene network reconstruction are Boolean networks [Ideker et al.,
1999], and Bayesian networks [Chang et al., 2011]. In mechanistic mod-
els, the gene network is usually represented using an oriented graph as in
Fig. 2(a). Another popular class of gene network models is based on ran-
dom forests [Huynh-Thu and Geurts, 2018], but they are usually difficult to
interpret as quantitative model because they only give a phenomenological
representation of interactions between genes, so we will not consider them
in this work.
Certain situations can alter the structure or expression of a gene. For example, in pathological situation, an alteration (mutation, deletion...) can occur in a gene and induce its loss of expression, leading to cell dysfunction and the development of an inflammatory disease or cancer. In experimental situations, it is possible to alter (suppress, decrease or increase) the expression of a specific gene involved in an interaction network, in order to better characterize its role in cell function. This experimental alteration can target the structure of a particular gene to modulate (decrease or increase) its expression, or directly target the expression of this gene inhibiting the messenger RNA (mRNA). These approaches are leading to the prospect of therapeutic applications in humans [de Mello et al., 2024].

The CRISPR/CAS9 approach is now one of the most widely used methods for altering gene structure. This method is derived from a natural antiviral defense mechanism found in certain bacteria. After a viral (phage) infection, these bacteria retain fragments of this virus's sequences (CRISPR sequences) in their genome. Upon re-infection with this virus, the viral DNA is recognized by a bacterial CRISPR RNA (crRNA). This crRNA then associates with a transactivating CRISPR RNA (tacrRNA) and recruits the CAS9 endonuclease, an enzyme that cleaves the viral sequence, inducing its degradation and the elimination of the virus [Boettcher and McManus, 2015]. Based on this mechanism, it is possible to synthesize a target gene-specific gRNA and associate it with the CAS9 endonuclease. After introduction of this molecule into a eukaryotic cell, the gRNA directs the endonuclease towards the target gene to induce its cleavage. During DNA repair mechanisms, a mutation altering the sequence of this gene is introduced, making it possible to specifically inhibit expression of this gene [Doudna and Charpentier, 2014, Vercauteren et al., 2024].

Some approaches directly target the expression (mRNA) of a particular gene in the GRN. These methods are also derived from physiological mechanisms of gene expression regulation and natural cell defense mechanisms based on RNA interference by microRNA (miRNA). When viral RNA is present in a cell, cellular miRNAs specifically recognize the sequence of these mRNAs. The miRNAs then bind to these viral RNAs, recruit the RISC/AGO2 molecular complex and induce their degradation [Lee et al., 1993]. The siRNA gene silencing method uses the same principle. After design of a 20-24 nucleotide synthetic siRNA specific to the gene sequence to be silenced, this siRNA is experimentally introduced into the cell. This siRNA then binds to the mRNA transcript of the target gene, recruits the RISC complex and uses the cellular machinery to induce degradation of the target RNA molecules [Fire et al., 1998, Davidson and McCray Jr, 2011]. Based on the same principle, the shRNA silencing method can also be used. In this case, a plasmid encoding a pre-shRNA is introduced into the cell and the latter matures into siRNA in the cell cytoplasm. This siRNA then recruits the RISC complex to degrade its specific target RNA [Mohr et al., 2014].
Other experimental methods can be used to increase the expression of a particular gene in a cell. For example, techniques for introducing recombinant DNA (the sequence of a gene) into a cell using a vector in order to induce expression of that gene and use the cellular machinery to produce the proteins of interest in that cell have been used for many years in the laboratory [Wolff et al., 1990] and are applicable for eukaryotic cells [Martínez-Puente et al., 2022].

It is therefore essential to predict the cellular response to a gene alteration in a network, in order to gain a better understanding of the development of pathologies caused by these alterations and acquired during life (cancer) or inherited (congenital diseases) [Chang et al., 2011, Vallat et al., 2013, Qin et al., 2014]. A model of gene alteration can also be used for design of experiment, for example to determine target genes to act on to obtain a certain biological effect or to discriminate between several possible gene networks [Ideker et al., 1999, Tegner et al., 2003, Dehghannasiri et al., 2015, Vundavilli et al., 2019, Bouvier et al., 2024]. When data of a genetic alteration experiment are available, a model of this alteration can also be used to validate or improve the statistical estimation of the initial model and/or the alteration model [Villa-Vialaneix et al., 2014, Omranian et al., 2016, Seçilmiş et al., 2022, Sarmah et al., 2022, Ishikawa et al., 2023].

Modeling studies of gene alteration have used a wide variety of methods and models that all belong to the framework of mechanistic models (Boolean networks [Azuaje et al., 2010, Ideker et al., 1999], Bayesian networks [Chang et al., 2011], linear regression models [Vallat et al., 2013, Qin et al., 2014, Villa-Vialaneix et al., 2014, Omranian et al., 2016], ODE-based models [Tegner et al., 2003, Goutsias and Lee, 2007]). Some models are static [Tegner et al., 2003, Qin et al., 2014], or focus on the stationary state of the network [Goutsias and Lee, 2007], or assume that data at different times are independent [Omranian et al., 2016]. Other studies only take into account alteration of the initial state in a dynamical model [Chang et al., 2011, Qin et al., 2014, Omranian et al., 2016]. The most advanced models are those that consider dynamic models and take into account changes in model structure over time [Vallat et al., 2013, Villa-Vialaneix et al., 2014, Dehghannasiri et al., 2015, Vundavilli et al., 2019, Ishikawa et al., 2023, Bouvier et al., 2024].

Our aim is to define general methods to construct dynamical models of gene alteration that are applicable to both quantitative and mechanistic settings and that are not limited to the simple change of initial condition. We set ourselves three goals: i) understanding how gene alteration modifies cellular behavior, ii) experiment design and iii) answering statistical issues such as estimation or validation.

We propose two different approaches, depending on whether the model is quantitative or mechanistic. In a so-called *conditional approach*, we define an alteration as a probabilistic continioning of the stochastic initial model by a certain event, so that for example the total knock-out of a gene is
taken into account by zeroing its expression. In particular, we show that this genetic alteration transforms a Gaussian graphical model into another Gaussian model whose mean values and correlations between gene expressions can be explicitly calculated, see Fig. 1(b). In a *mechanistic approach*, alteration takes the form of a modification of the model’s parameters, so that for example the total knock-out of a gene modifies the gene network as illustrated in Fig. 2(b).

Conditional and mechanistic methods are described in Section 2. They predict the effect of a gene alteration from an initial model inferred from a dataset obtained without gene alteration. In Section 3.1, we develop a first application of our modeling to the design of the most relevant alteration experiments (according to a certain criterion). When a gene alteration dataset is available, it is used to estimate the parameters of the model describing the alteration and/or to improve the estimation of the parameters of the initial model (Section 3.2). We then describe these methods in detail within the framework of three standard classes of gene network models: an example of quantitative but non-mechanistic model, the Gaussian graphical models, in Section 4.1, and two examples of mechanistic models in Section 4.2: Bayesian networks in Section 4.2.1, and a penalized linear regression model from [Champagnat et al., 2025] in Section 4.2.2.

# 2 Modeling gene network alteration

## 2.1 The mathematical framework

We assume that we have a gene expression dataset (e.g. RNA-seq data) that we will call in the following *initial dataset*.

$$ \mathbf{x} = (x_{i,k,p}, 1 \le i \le N, 1 \le k \le K, 1 \le p \le P), \quad (2.1) $$

where

* $P$ is the number of experiments carried out.

* $N$ is the number of genes whose expression is measured in each experiment.

* For each experiment, $K$ measurements are realized at times $t_1 = 0 < t_2 < \dots < t_K$.

* $x_{i,k,p}$ is a real number which measures the expression (after normalization) of the gene $i$ at the $k$-th time of measurement during the $p$-th experiment. In the following, we will restrict ourselves only to absolute expression data $x_{i,k,p} \in \mathbb{R}_+$, to make it simpler to model the effect of gene alteration (e.g., gene knock-out reduces the absolute expression of the gene to 0).
Our statistical model is based on the assumption that the matrices

$$ (x_{i,k,1})_{1 \le i \le N, 1 \le k \le K}, (x_{i,k,2})_{1 \le i \le N, 1 \le k \le K}, \dots, (x_{i,k,P})_{1 \le i \le N, 1 \le k \le K} $$

are realizations of independent and identically distributed (i.i.d.) experiments of the random matrix

$$ \mathbf{X} = (X_{i,k})_{1 \le i \le N, 1 \le k \le K}. \qquad (2.2) $$

We suppose that the family of laws of $\mathbf{X}$ is parametrized by $\theta$ belonging to a set $\Theta$. The statistical model is then written as :

$$ (\mathbb{P}_{\theta})_{\theta \in \Theta}, \qquad (2.3) $$

where $\mathbb{P}_{\theta}$ is the distribution of $\mathbf{X}$ when the parameter equals $\theta$. We will give in Section 4 several classes of models that have been considered in the literature. For example, in the Gaussian graphical model, the family $(\mathbb{P}_{\theta})_{\theta \in \Theta}$ is a Gaussian family (cf. Section 4.1). A model of the type (2.3) will be called *quantitative*.

Most models are associated with the notion of a gene network, which can be seen as a component of the model parameter. The interpretation of the gene network may vary from one model to another. We will come back to the notion of gene network in Section 4.

We define *mechanistic* models as those giving the expression of each gene at time $t_k$ as an explicit function of the expressions of all genes at the preceding times $t_1, \dots, t_{k-1}$ and an independent random noise. In other words,

$$ X_{ik} = F_{ik}(X_{1,[1,k-1]}, X_{2,[1,k-1]}, \dots, X_{N,[1,k-1]}; U_{ik}), \qquad (2.4) $$

for all $i \in \{1, \dots, N\}$ and $k \in \{2, \dots, K\}$, where, for all $1 \le i \le N$ and $\ell \le K$,

$$ X_{i,[1,\ell]} = (X_{ik})_{1 \le k \le \ell}, \qquad (2.5) $$

$F_{ik}$ is a function and the random variables $U_{ik}$ are independent and are also independent from the initial condition $(X_{i1})_{1 \le i \le N}$. The function $F_{ik}$ describes how the amounts of RNA expressed at times $t_1, \dots, t_{k-1}$ influence the expression of the gene $i$ at time $t_k$. Models that cannot be written in this form will be referred to as *non-mechanistic* models.

In a mechanistic model, the parameter $\theta$ can be defined as the vector formed by the family of functions $(F_{ik})_{1 \le i \le N, 1 \le k \le K}$, the law of the vector of initial expressions $(X_{i1})_{1 \le i \le N}$ and the law of the random noise $(U_{ik})_{1 \le i \le N, 2 \le k \le K}$. Such models are typical in causal inference [Pearl, 2016]. Some examples are given in section 4.
<figure><img src="image_4.png" /><figcaption>Figure 3: Prediction of gene alteration</figcaption></figure>

## 2.2 Probabilistic approach for the prediction of gene alteration

We assume that the inference of the statistical model (2.3) from the initial data $x$ has been carried out. We denote by $\hat{\theta}_P$ the estimator of $\theta$ obtained and we call $\mathbb{P}_{\hat{\theta}_P}$ the initial prediction. Our goal is to obtain a prediction of a gene alteration from the initial prediction.

Our approach consists in associating with any $\theta \in \Theta$, a model $\mathbb{P}_{\theta}^{\text{alt}}$ describing the law of expressions of genes after alteration. The prediction of gene alteration is then given by $\mathbb{P}_{\hat{\theta}_P}^{\text{alt}}$ (cf. Figure 3). Let us note that our approach has retained the same set of parameters for the prediction than the initial model, but generally the model $\mathbb{P}_{\theta}^{\text{alt}}$ model does not belong to the family $(\mathbb{P}_{\theta'})_{\theta' \in \Theta}$.

We consider two approaches to construct the laws $\mathbb{P}_{\theta}^{\text{alt}}$, one based only on the laws $(\mathbb{P}_{\theta})_{\theta \in \Theta}$, the other on a mechanistic interpretation of this model.

For reasons of simplicity, we present the principle of these approaches in the case where the alteration acts only on gene 1 and takes the form either of a gene knock-down or overexpression. We introduce a parameter $\beta \ge 0$ representing the alteration of expression of gene 1: if $\beta \in [0, 1]$, it represents the fraction of expression of gene 1 that remains after reduction of expression; if $\beta > 1$, then $\beta - 1$ represents the relative overexpression of gene 1. If $\beta = 0$, gene 1 is completely knocked-out after alteration; $\beta \in (0, 1)$ corresponds to a knock-down of gene 1. The case of alteration of a gene $j \neq 1$ or of several genes simultaneously, possibly with different values of $\beta$, is easily deduced.

### 2.2.1 Conditioning approach

We denote by $\mathbb{P}_{\theta}^{\text{alt}, \beta}$ the prediction of silencing by conditional law. This law is constructed by assuming that

* the law of $(X_{1k})_{1 \le k \le K}$ under $\mathbb{P}_{\theta}^{\text{alt}, \beta}$ is the law of $(\beta X_{1k})_{1 \le k \le K}$ under $\mathbb{P}_{\theta}$

* the conditional distribution of $(X_{ik})_{2 \le i \le N, 1 \le k \le K}$ given $(X_{1k})_{1 \le k \le K}$ is the same under $\mathbb{P}_{\theta}^{\text{alt}, \beta}$ or $\mathbb{P}_{\theta}$.

If under the probability $\mathbb{P}_{\theta}$, the vector $(X_{1k})_{1 \le k \le K}$ has a density function $p_{\theta}$ on $\mathbb{R}_{+}^K$, then for all $\beta > 0$, under $\mathbb{P}_{\theta}^{\text{alt}, \beta}$, $(X_{1k})_{1 \le k \le K}$ has for density

$$p_{\theta} \left( \frac{x_{11}}{\beta}, \dots, \frac{x_{1K}}{\beta} \right) \frac{1}{\beta^K}$$
and

$$ \begin{aligned} & \mathbb{P}_{\theta}^{\text{alt},\beta}(\mathbf{X} \in \mathrm{d}\mathbf{x}) \\ & = \mathbb{P}_{\theta} \left( X_{2,[1,K]} \in \mathrm{d}x_{2,[1,K]}, \dots, X_{N,[1,K]} \in \mathrm{d}x_{N,[1,K]} \mid X_{1,[1,K]} = x_{1,[1,K]} \right) \\ & \qquad \times p_{\theta} \left( \frac{x_{1,[1,K]}}{\beta} \right) \frac{1}{\beta^K} \mathrm{d}x_{1,[1,K]}, \end{aligned} $$

where $x_{i,[1,K]}$ stands for $(x_{i,1}, \dots, x_{i,K})$, $2 \le i \le N$.

In the case of a knock-out of gene 1 ($\beta = 0$), we obtain

$$ \begin{aligned} & \mathbb{P}_{\theta}^{\text{alt},\beta}(\mathbf{X} \in \mathrm{d}\mathbf{x}) \\ & = \mathbb{P}_{\theta} \left( X_{2,[1,K]} \in \mathrm{d}x_{2,[1,K]}, \dots, X_{N,[1,K]} \in \mathrm{d}x_{N,[1,K]} \mid X_{1,[1,K]} = 0 \right) \\ & \qquad \times \delta_0(\mathrm{d}x_{1,[1,K]}), \end{aligned} $$

where $\delta_0$ is the Dirac measure at 0 in $\mathbb{R}^K$.

### 2.2.2 Mechanistic approach

We start from the mechanistic model (2.4). The mechanistic model of gene 1 alteration is based on the assumption that the functions $F_{ik}$ and the law of the random vector $(U_{ik})_{1 \le i \le N, 1 \le k \le K}$ are not modified by the alteration. Moreover the vector $(U_{ik})_{1 \le i \le N, 1 \le k \le K}$ remains independent of the initial expressions $(X_{i1})_{1 \le i \le N}$. This leads to the model

$$ \begin{cases} X_{11}^{\text{alt},\beta} = \beta X_{11}, & X_{i1}^{\text{alt},\beta} = X_{i1}, \forall i \in \{2, \dots, N\}, \\ X_{ik}^{\text{alt},\beta} = F_{ik} \left( \beta X_{1,[1,k-1]}^{\text{alt},\beta}, X_{2,[1,k-1]}^{\text{alt},\beta}, \dots, X_{N,[1,k-1]}^{\text{alt},\beta} ; U_{ik} \right), & \end{cases} $$

for all $1 \le i \le N$ and $2 \le k \le K$. The probability $\mathbb{P}_{\theta}^{\text{alt},\beta}$ is the law of the random vector $\mathbf{X}^{\text{alt},\beta}$.

Note that the previous gene alteration model is also mechanistic. In the case of a knock-out of gene 1 ($\beta = 0$), we obtain

$$ X_{ik}^{\text{alt},0} = F_{ik} \left( 0, X_{2,[1,k-1]}^{\text{alt},0}, \dots, X_{N,[1,k-1]}^{\text{alt},0} ; U_{ik} \right), \quad (2.6) $$

for all $i \in \{2, \dots, N\}$ and $k \in \{2, \dots, K\}$.

### 2.2.3 A simple example

To illustrate the difference between the mechanistic and conditioning approaches, let us consider an example involving 3 genes and two dates ($N = 3, K = 2$), with dependency relationships between genes shown in Figure 4 (compact representation without time-dependence) and Figure 5 (full Bayesian dependency network).
<figure><img src="image_5.png" /><figcaption>Figure 4: Dependency network between the expressions of genes 1, 2 and 3</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Bayesian dependency network of $(X_{ik})_{1 \le i \le 3, 1 \le k \le 2}$ under $\mathbb{P}_{\theta}$.</figcaption></figure>

Let us consider the knock-out of gene 1 ($\beta = 0$). For simplicity, we assume that the random variables $X_{ik}$ take values in a discrete set $\mathcal{X}$ containing 0. In the formalism of Bayesian networks (cf. section 4.2.1), Figure 5 means that for all $\theta \in \Theta$, the law of the random vector $(X_{ik})_{1 \le i \le 3, 1 \le k \le 2}$ is

$$ \begin{aligned} & \mathbb{P}_{\theta}(X_{11} = x_{11}, X_{21} = x_{21}, X_{31} = x_{31}, X_{12} = x_{12}, X_{22} = x_{22}, X_{32} = x_{32}) \\ & = \mathbb{P}_{\theta}(X_{11} = x_{11})\mathbb{P}_{\theta}(X_{21} = x_{21})\mathbb{P}_{\theta}(X_{31} = x_{31})\mathbb{P}_{\theta}(X_{22} = x_{22} | X_{21} = x_{21}) \\ & \quad \times \mathbb{P}_{\theta}(X_{12} = x_{12} | X_{11} = x_{11}, X_{21} = x_{21}) \\ & \quad \times \mathbb{P}_{\theta}(X_{32} = x_{32} | X_{11} = x_{11}, X_{31} = x_{31}), \end{aligned} $$

for all $(x_{ik})_{1 \le i \le 3, 1 \le k \le 2}$ in $\mathcal{X}^6$. We suppose that $\mathbb{P}_{\theta}(X_{11} = 0, X_{12} = 0) > 0$.

**Conditioning approach** Let us denote by $\mathbb{P}_{\theta}^{\text{alt,cond}}$ the law $\mathbb{P}_{\theta}^{\text{alt},0}$ obtained using the conditioning approach. By definition:

$$ \begin{aligned} & \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{21} = x_{21}, X_{31} = x_{31}, X_{22} = x_{22}, X_{32} = x_{32}) \\ & = \mathbb{P}_{\theta}(X_{21} = x_{21}, X_{31} = x_{31}, X_{22} = x_{22}, X_{32} = x_{32} | X_{11} = X_{12} = 0) \end{aligned} $$

Using the independence of the random variables $X_{11}$, $X_{21}$ and $X_{31}$ we obtain after some computation

$$ \begin{aligned} & \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{21} = x_{21}, X_{31} = x_{31}, X_{22} = x_{22}, X_{32} = x_{32}) \\ & = \mathbb{P}_{\theta}(X_{31} = x_{31})\mathbb{P}_{\theta}(X_{32} = x_{32} | X_{11} = 0, X_{31} = x_{31}) \\ & \quad \times \mathbb{P}_{\theta}(X_{21} = x_{21} | X_{11} = X_{12} = 0)\mathbb{P}_{\theta}(X_{22} = x_{22} | X_{21} = x_{21}). \end{aligned} \quad (2.7) $$

We derive the following results

$$ \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{31} = x_{31}) = \mathbb{P}_{\theta}(X_{31} = x_{31}), $$

$$ \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{32} = x_{32} | X_{31} = x_{31}) = \mathbb{P}_{\theta}(X_{32} = x_{32} | X_{11} = 0, X_{31} = x_{31}), $$

$$ \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{21} = x_{21}) = \mathbb{P}_{\theta}(X_{21} = x_{21} | X_{12} = X_{11} = 0), \quad (2.8) $$

$$ \mathbb{P}_{\theta}^{\text{alt,cond}}(X_{22} = x_{22} | X_{21} = x_{21}) = \mathbb{P}_{\theta}(X_{22} = x_{22} | X_{21} = x_{21}). $$
In other words, the Bayesian dependency network associated with the model $\mathbb{P}_{\theta}^{\text{alt,cond}}$ is given by Figure 6. This network is obtained by deleting from the network in Figure 5, all edges originating from or directed towards one of the two vertices $X_{11}$ and $X_{12}$. In the sequel, we call the edges incident edges to $X_{11}$ and to $X_{12}$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Bayesian dependency network for $(X_{ik})_{1\le i\le 3, 1\le k\le 2}$ under $\mathbb{P}_{\theta}^{\text{alt,cond}}$.</figcaption></figure>

**Mechanistic approach** Let us denote by $\mathbb{P}_{\theta}^{\text{alt,meca}}$ the law $\mathbb{P}_{\theta}^{\text{alt},0}$ obtained using the mechanistic approach.

The mechanistic approach requires the model to be written as (2.4). Consider the random variables $(X_{11}^{\text{meca}}, X_{21}^{\text{meca}}, X_{31}^{\text{meca}}) = (X_{11}, X_{21}, X_{31})$ under $\mathbb{P}_{\theta}$. Let $U_{12}, U_{22}, U_{32}$ i.i.d. random variables uniformly distributed on $[0, 1]$ and independent of $(X_{11}^{\text{meca}}, X_{21}^{\text{meca}}, X_{31}^{\text{meca}})$. For any $x_{11}, x_{21}, x_{31}$, we introduce the functions:

* $F_{22}(x_{21}; \cdot)$ is the generalized inverse of the distribution function of $X_{22}$ under $\mathbb{P}_{\theta}$ conditional on $X_{21} = x_{21}$;

* $F_{12}(x_{11}, x_{21}; \cdot)$ is the generalized inverse of the distribution function of $X_{12}$ conditional on $X_{11} = x_{11}$ and $X_{21} = x_{21}$;

* $F_{32}(x_{11}, x_{31}; \cdot)$ is the generalized inverse of the distribution function of $X_{32}$ conditional on $X_{11} = x_{11}$ and $X_{31} = x_{31}$.

We now define the vector $(X_{ik}^{\text{meca}}, 1 \le i \le 3, 1 \le k \le 2)$ as follows

$$X_{22}^{\text{meca}} = F_{22}(X_{21}^{\text{meca}}; U_{22}), \quad X_{12}^{\text{meca}} = F_{12}(X_{11}^{\text{meca}}, X_{21}^{\text{meca}}; U_{12})$$

and

$$X_{32}^{\text{meca}} = F_{32}(X_{11}^{\text{meca}}, X_{31}^{\text{meca}}; U_{32}).$$

Then $(X_{ik}^{\text{meca}}, 1 \le i \le 3, 1 \le k \le 2)$ has the same distribution as $(X_{ik}, 1 \le i \le 3, 1 \le k \le 2)$ under $\mathbb{P}_{\theta}$.

According to (2.6), the mechanistic approach to model the knock-out of gene 1 leads to the model

$$(X_{21}^{\text{alt,meca}}, X_{31}^{\text{alt,meca}}) = (X_{21}^{\text{meca}}, X_{31}^{\text{meca}}), \\ X_{22}^{\text{alt,meca}} = F_{22}(X_{21}^{\text{alt,meca}}; U_{22}) \text{ and } X_{32}^{\text{alt,meca}} = F_{32}(0, X_{31}^{\text{alt,meca}}; U_{32}).$$
From the definition of functions $F_{22}$ and $F_{32}$, we deduce that

$$ \begin{aligned} & \mathbb{P}_{\theta}^{\text{alt,meca}} (X_{21}^{\text{alt,meca}} = x_{21}, X_{31}^{\text{alt,meca}} = x_{31}, X_{22}^{\text{alt,meca}} = x_{22}, X_{32}^{\text{alt,meca}} = x_{32}) \\ & = \mathbb{P}_{\theta}(X_{21} = x_{21})\mathbb{P}_{\theta}(X_{31} = x_{31})\mathbb{P}_{\theta}(X_{22} = x_{22} | X_{21} = x_{21}) \\ & \qquad \mathbb{P}_{\theta}(X_{32} = x_{32} | X_{11} = 0, X_{31} = x_{31}). \end{aligned} $$

From this we deduce

$$ \mathbb{P}_{\theta}^{\text{alt,meca}} (X_{21}^{\text{alt,meca}} = x_{21}) = \mathbb{P}_{\theta}(X_{21} = x_{21}). $$

A comparison with (2.8) shows that the two models are not identical. This difference can be explained by the fact that conditioning the expression of gene 2 at time 2 on the fact that the expressions of gene 1 are zero is a different constraint from that of the mechanistic model: conditioning has the effect of biasing the expression of genes located upstream of gene 1. More will be said on the differences between the two models in Section 4.2.1.

# 3 Applications of the gene alteration modeling

We consider two distinct settings for the application of our method. First, on the basis of the initial dataset only, we consider the problem of design of experiments, i.e. the choice of gene alteration experiments among several possibilities, in order to best achieve a given objective (Section 3.1). Second, assuming that gene alteration data were obtained in addition to the initial dataset, we can define a statistical model of these two datasets, allowing to address statistical issues such as validation or estimation (Section 3.2). We will apply these methods to typical models of gene regulatory networks in Section 4.

## 3.1 Application for the design of experiments

Let us consider the case where several possible gene alteration experiments $alt_1, \dots, alt_n$ are possible and the objective is to favor a certain cell behavior such as proliferation [Schleiss et al., 2021], or more generally maximize some utility [Venkat et al., 2017, Bouvier et al., 2024]. In a biological context different from gene regulatory networks, maximizations of measures of vulnerability [Zhu et al., 2016] or robustness [Bianconi et al., 2015] were also considered using methodological approaches. In the sequel, we consider the objective of reducing proliferation.

The search for a predictor of cellular behavior such as proliferation is a difficult task that can be addressed using various methods [Azuaje et al., 2010, Bianconi et al., 2015, Zhu et al., 2016, Venkat et al., 2017], e.g. using a learning method or known biological information on the genes involved in the
cellular proliferation program. We assume here that we have a model of (possibly nonlinear) regression of the response variable of proliferation Y in terms of the explanatory variables given by the gene expression matrix X as follows

$$Y = H(\mathbf{X}),$$

where H can be random.

We then choose the experiment i that minimizes the mean proliferation scores

$$\mathbb{E}_{\hat{\theta}_P}^{\text{alt}_1} H(\mathbf{X}), \dots, \mathbb{E}_{\hat{\theta}_P}^{\text{alt}_n} H(\mathbf{X}).$$

A simple example of predictor H is given by

$$H(\mathbf{X}) = \sum_{i \in G} \sum_{k \in T} X_{i,k},$$

where $G \subset \{1, \dots, N\}$ is a set of genes involved in the proliferation program and $T \subset \{1, \dots, K\}$ is the set of relevant dates.

Another method of selection of gene alteration experiments would be to choose the experiment that maximizes the number of down-regulated genes in G at the dates in T. The gene i is said to be up-regulated (resp. down-regulated) at time k for the experiment alt$_{\ell}$ if

$$\mathbb{E}_{\hat{\theta}_P} X_{i,k} < \mathbb{E}_{\hat{\theta}_P}^{\text{alt}_\ell} X_{i,k} \quad (\text{resp. } \mathbb{E}_{\hat{\theta}_P} X_{i,k} > \mathbb{E}_{\hat{\theta}_P}^{\text{alt}_\ell} X_{i,k}).$$

We then define the predictor of proliferation associated with experiment alt$_{\ell}$ as

$$S_{\ell} = \text{Card} \left\{ (i, k) \in G \times T : \mathbb{E}_{\hat{\theta}_P} X_{i,k} > \mathbb{E}_{\hat{\theta}_P}^{\text{alt}_{\ell}} X_{i,k} \right\},$$

where Card(A) is the number of elements in the set A, and we select the gene alteration experiment that maximizes $S_1, \dots, S_n$.

Another possibility is to choose the alteration experiment that most closely approximates gene expressions to given mean target values $m_{ik}^{\text{target}}$ (e.g. measured in non proliferative cells). In this case

$$S_{\ell} = \sum_{i \in G} \sum_{k \in T} \left| \mathbb{E}_{\hat{\theta}_P}^{\text{alt}_{\ell}} X_{i,k} - m_{ik}^{\text{target}} \right|$$

and we select the alteration experiment that minimizes $S_1, \dots, S_n$.

## 3.2 Application to the statistical modeling of initial and gene alteration datasets

We now assume that we have a second dataset from a gene alteration experiment with parameter $\beta$

$$\mathbf{x}^{\text{alt}} = (x_{i,k,p}^{\text{alt}}, 1 \le i \le N, 1 \le k \le K, 1 \le p \le P^{\text{alt}}),$$
<figure><img src="image_8.png" /><figcaption>Figure 7: Gene alteration model validation</figcaption></figure>

called *alt dataset*. Using one of the approaches described in Section 2.2, we build a statistical model of the initial and alt datasets under which the likelihood of the data writes in the discrete case as

$$ \prod_{p=1}^{P} \mathbb{P}_{\theta}(\mathbf{X} = \mathbf{x}_p) \prod_{q=1}^{P^{\text{alt}}} \mathbb{P}_{\theta}^{\text{alt}, \beta}(\mathbf{X}^{\text{alt}} = \mathbf{x}_q^{\text{alt}}), $$

where $\mathbf{x}_p = (x_{i,k,p})_{1 \le i \le N, 1 \le k \le K}$ and $\mathbf{x}_q^{\text{alt}} = (x_{i,k,q}^{\text{alt}})_{1 \le i \le N, 1 \le k \le K}$.

We can then use this statistical model to answer statistical questions. In the sequel, we focus on the validation of the *alt* model (Section 3.2.1) and on the improvement of estimation of the initial model using *alt* data (Section 3.2.2). Recall that $\hat{\theta}_P$ is an estimator of $\theta$ obtained from the *initial dataset*.

### 3.2.1 Validation of a model of gene alteration

Our aim is to check if the gene alteration model $\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}$ correctly accounts for the observations $\mathbf{x}^{\text{alt}}$. From a statistical point of view, this corresponds to a goodness-of-fit test (Figure 7). Due to the small amount of observation, it is usually impossible in practice to test if $\mathbf{x}^{\text{alt}}$ follows the distribution $\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}$. This explains why authors usually restrict to low-dimensional test statistics relevant for the biological application, typically either relating to gene expressions [Vallat et al., 2013], or to the presence or absence of edges in the network [Wagner, 2001, Markowetz et al., 2005].

Below we describe the goodness-of-fit test used in [Vallat et al., 2013], which can be adapted to any quantitative model. The idea is to test whether the number of up- and down-regulated genes correctly predicted by the model $\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}$ is significantly greater than what would be obtained if up- and down-regulated genes were chosen at random. We consider an alteration of gene 1 only, with parameter $\beta$.

More precisely, for any positive $\epsilon$, we define a partition of gene-date pairs
based on the two probabilities $\mathbb{P}_{\hat{\theta}_P}$ and $\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}$ as follows:

$$ \mathcal{U} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}}(X_{i,k}) \ge (1+\varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}, $$

$$ \mathcal{D} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}}(X_{i,k}) \le (1-\varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}, $$

$$ \mathcal{C} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \right. \tag{3.1} $$

$$ (1 - \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) < \mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}}(X_{i,k}) < (1 + \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}. \tag{3.2} $$

The set $\mathcal{U}$ (resp. $\mathcal{D}$, resp. $\mathcal{C}$) is the set of up-regulated (resp. down-regulated, constant) gene-date pairs, i.e. those whose average expression is higher in the alt model than in the initial model at the date considered. The sets $\mathcal{U}$, $\mathcal{D}$ and $\mathcal{C}$ contain neither gene 1 nor the first date, because gene expression is identical in both models at this date. Note that these sets are constructed solely on the basis of the estimated parameter $\hat{\theta}_P$, considered here as a constant.

Then, we consider the following test statistic constructed from a alt sample $\mathbf{X}^{\text{alt}} = (X_1^{\text{alt}}, \dots, X_{P^{\text{alt}}}^{\text{alt}})$ of size $P^{\text{alt}}$:

$$ S = S(\mathbf{X}^{\text{alt}}) := \frac{\text{Card}(\mathcal{U} \cap \hat{\mathcal{U}}^{\text{alt}}) + \text{Card}(\mathcal{D} \cap \hat{\mathcal{D}}^{\text{alt}}) + \text{Card}(\mathcal{C} \cap \hat{\mathcal{C}}^{\text{alt}})}{(N-1)(K-1)}, \tag{3.3} $$

where

$$ \hat{\mathcal{U}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \bar{X}_{i,k}^{\text{alt}} \ge (1 + \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}, $$

$$ \hat{\mathcal{D}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \bar{X}_{i,k}^{\text{alt}} \le (1 - \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}, $$

$$ \hat{\mathcal{C}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \right. \\ \qquad \left. (1 - \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) < \bar{X}_{i,k}^{\text{alt}} < (1 + \varepsilon)\mathbb{E}_{\mathbb{P}_{\hat{\theta}_P}}(X_{i,k}) \right\}, $$

with

$$ \bar{X}_{i,k}^{\text{alt}} := \frac{1}{P^{\text{alt}}} \sum_{p=1}^{P^{\text{alt}}} (X_{p}^{\text{alt}})_{i,k}. $$

The S statistic counts the proportion of mean up-regulated, down-regulated or constant gene-date pairs in the alt sample that agree with the model. We consider the null and alternative hypotheses

(H0) the distribution of $X^{\text{alt}}$ is $\mathbb{P}_{\hat{\theta}_P}^{\text{alt},\beta}$,

(H1) the law of $S(\mathbf{X}^{\text{alt}})$ is $\mu_0$,
where $\mu_0$ is the distribution of

$$S' := \frac{\mathrm{Card}(\mathcal{U} \cap \mathcal{U}') + \mathrm{Card}(\mathcal{D} \cap \mathcal{D}') + \mathrm{Card}(\mathcal{C} \cap \mathcal{C}')}{(N-1)(K-1)},$$

where

$$\mathcal{U}' := \{(i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : Z_{i,k} = u\},$$

$$\mathcal{D}' := \{(i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : Z_{i,k} = d\},$$

$$\mathcal{C}' := \{(i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : Z_{i,k} = c\}$$

and $(Z_{i,k}, (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\})$ is a sequence of i.i.d. random variable with values in $\{u, d, c\}$ and distribution $(p_u, p_d, p_c)$ where

$$p_u := \frac{\mathrm{Card}(\mathcal{U})}{(N-1)(K-1)}, \quad p_d := \frac{\mathrm{Card}(\mathcal{D})}{(N-1)(K-1)}, \quad p_c := 1 - p_u - p_d.$$

Note that $S'$ is obtained from $S$ by replacing $\hat{\mathcal{U}}^{\mathrm{alt}}$, $\hat{\mathcal{D}}^{\mathrm{alt}}$ and $\hat{\mathcal{C}}^{\mathrm{alt}}$ with “purely” random sets $\mathcal{U}'$, $\mathcal{D}'$ and $\mathcal{C}'$. In particular, the random variables $\mathrm{Card}(\mathcal{U} \cap \mathcal{U}')$, $\mathrm{Card}(\mathcal{D} \cap \mathcal{D}')$ and $\mathrm{Card}(\mathcal{C} \cap \mathcal{C}')$ are independent, with binomial laws of parameters $(\mathrm{Card}(\mathcal{U}), p_u)$, $(\mathrm{Card}(\mathcal{D}), p_d)$ and $(\mathrm{Card}(\mathcal{C}), p_c)$, respectively.

Hypothesis (H0) will be rejected for low values of the statistic $S$. Hence the critical region associated with the test of significance level $\alpha$ (typically, $\alpha = 5\%$) is an interval $R = [0, u_\alpha]$, where $u_\alpha > 0$ is determined by the condition

$$\mathbb{P}_{\hat{\theta}_p}^{\mathrm{alt}, \beta} (S(\mathbf{X}^{\mathrm{alt}}) \in R) \le \alpha, \qquad (3.4)$$

If $\hat{s} := S(\mathbf{x}^{\mathrm{alt}})$ is smaller than $u_\alpha$, (H0) is rejected. Otherwise, the $p$-value of the test is $\mathbb{P}_{\hat{\theta}_p}^{\mathrm{alt}, \beta} (S(\mathbf{X}^{\mathrm{alt}}) \le \hat{s})$ and the power of the test is given by the type II error $\mu_0(R^c)$, which will be small if the sample size is large enough.

### 3.2.2 Improved estimation of initial model parameters

We now want to understand how the $\mathbf{x}^{\mathrm{alt}}$ dataset can be used to improve the estimation of the parameter $\theta$ of the initial model $(\mathbb{P}_{\theta})_{\theta \in \Theta}$. We consider here the usual cases where parameter inference is based either on likelihood maximization or on a posteriori likelihood maximization methods.

Assume that for any $\theta \in \Theta$, the probability measures $\mathbb{P}_{\theta}$ and $\mathbb{P}_{\theta}^{\mathrm{alt}, \beta}$ on $\mathbb{R}^{NK}$ are absolutely continuous with respect to a reference measure $\rho$. Let $f_\theta$ and $g_\theta^{\mathrm{alt}, \beta}$ denote respectively the densities of these laws. In the case of a posteriori maximization, we assume that the a priori law has a density $p(\theta)$ with respect to a reference measure on $\Theta$. In the case of likelihood maximization, let $p(\theta) = 1$, for any $\theta \in \Theta$.

The likelihood of the data $(\mathbf{x}, \mathbf{x}^{\mathrm{alt}})$ is then

$$p(\theta) \prod_{p=1}^{P} f_{\theta}(\mathbf{x}_p) \prod_{q=1}^{P^{\mathrm{alt}}} g_{\theta}^{\mathrm{alt}, \beta}(\mathbf{x}_q^{\mathrm{alt}}). \qquad (3.5)$$
The new estimate of parameter $\theta$ is obtained by maximizing the previous quantity with respect to $\theta \in \Theta$ (see the examples of Section 4).

# 4 Gene alteration modeling for quantitative and mechanistic models

In the following, we will consider Gaussian graphical models (which are non-mechanistic, but quantitative models) and two mechanistic models: Bayesian networks and penalized linear regressions. Our goal is to implement for each model the methods of Sections 2 and 3.

## 4.1 Gaussian graphical models

There are many variants of these models but we will concentrate on the simplest, for which we show below that the model of gene knock-out is explicit and its impact on the gene regulatory network can be fully characterized. As an application, we will implement the model validation test of Section 3.2.1. We assume that the vector $\mathbf{X}$ of gene expressions defined by (2.2) follows a multivariate normal distribution $N(m, \Sigma)$ where $m = (m_{ik})_{1 \le i \le N, 1 \le k \le K}$ is the mean and $\Sigma = (\sigma_{ikj\ell})_{1 \le i,j \le N, 1 \le k, \ell \le K}$ the covariance matrix that is supposed to be positive definite. Under these conditions, $\mathbb{P}_{\theta} = N(m, \Sigma)$ where the model parameter is $\theta = (m, \Sigma)$. The precision matrix is $\Sigma^{-1}$ and it is well-known (cf. e.g. [Kendall and Stuart, 1973, Eq. (27.6)]) that the partial correlation matrix is given by

$$C = -D\Sigma^{-1}D, \quad \text{where} \quad D = \text{diag}\left(\frac{1}{\sqrt{\Sigma_{11}^{-1}}}, \dots, \frac{1}{\sqrt{\Sigma_{NN}^{-1}}}\right).$$

The gene network is the graph $G = (E, V)$ where the set $E$ of vertices is made up of gene-time pairs $(i, k)$ where $1 \le i \le N$ and $1 \le k \le K$, and $V$ is the set of pairs $\{(i, k), (j, \ell)\}$ such that $c_{ikj\ell} \neq 0$. In other words, gene $i$ at time $t_k$ is linked to gene $j$ at time $t_\ell$ if and only $X_{ik}$ et $X_{j\ell}$ are not conditionally independent given all other gene expressions.

For simplicity, we consider the alteration of gene 1 of parameter $\beta \ge 0$. Since the model is non-mechanistic, only the conditioning approach of Section 2.2.1 is relevant. To simplify notations, we order pairs $(i, k) \in \{1, \dots, N\} \times \{1, \dots, K\}$ according to the lexicographic order and we decompose by blocks the matrix $\Sigma$ and the vector $m$ as follows

$$\Sigma = \begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{pmatrix} \quad \text{and} \quad m = \begin{pmatrix} m_1 \\ m_2 \end{pmatrix},$$

where $\Sigma_{11}$ (resp. $\Sigma_{12}$, $\Sigma_{21}$, $\Sigma_{22}$) is a matrix $K \times K$ (resp. $K \times (N-1)K$, $(N-1)K \times K$, $(N-1)K \times (N-1)K$) and $m_1$ (resp. $m_2$) is a $K$ (resp. $(N-1)K$)-dimensional vector.
**Proposition 4.1.** (i) The law $\mathbb{P}_{(m,\Sigma)}^{alt,\beta}$ is Gaussian

$$ \mathbb{P}_{(m, \Sigma)}^{alt, \beta} = \mathcal{N} \left( \begin{pmatrix} \beta m_1 \\ \mu(\beta) \end{pmatrix}, \begin{pmatrix} \beta^2 \Sigma_{11} & \beta^2 \Sigma_{12} \\ \beta^2 \Sigma_{21} & \Gamma(\beta) \end{pmatrix} \right) $$

where

$$ \mu(\beta) = m_2 - (1 - \beta)\Sigma_{21}\Sigma_{11}^{-1}m_1 \quad (4.1) $$

and

$$ \Gamma(\beta) = \Sigma_{22} - (1 - \beta^2)\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}. \quad (4.2) $$

(ii) In addition, $\Gamma(0)^{-1} = (\Sigma^{-1})_{22}$.

Proposition 4.1 will be proved in Appendix A. Note that, if $\beta = 0$, the alt model is given by

$$ \mathbb{P}_{(m, \Sigma)}^{alt, 0} = \mathcal{N} \left( \begin{pmatrix} 0 \\ \mu(0) \end{pmatrix}, \begin{pmatrix} 0 & 0 \\ 0 & \Gamma(0) \end{pmatrix} \right) $$

In particular, Item (ii) means that the gene network of the knock-out model ($\beta = 0$) is identical to the initial network in which gene 1 and all incidental edges are deleted.

The previous result allows us to explicitly determine the sets of up-regulated, down-regulated and constant gene-date pairs defined in (3.2) which take the form, after some algebra:

$$ \mathcal{U} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \varepsilon(\hat{m}_2)_{i,k} \leq -(\hat{\Sigma}_{21}\hat{\Sigma}_{11}^{-1}\hat{m}_1)_{i,k} \right\}, $$

$$ \mathcal{D} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \varepsilon(\hat{m}_2)_{i,k} \leq (\hat{\Sigma}_{21}\hat{\Sigma}_{11}^{-1}\hat{m}_1)_{i,k} \right\}, $$

$$ \mathcal{C} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : -\varepsilon(\hat{m}_2)_{i,k} < (\hat{\Sigma}_{21}\hat{\Sigma}_{11}^{-1}\hat{m}_1)_{i,k} < \varepsilon(\hat{m}_2)_{i,k} \right\}, $$

where $\varepsilon > 0$ and $\hat{\theta}_P = (\hat{m}, \hat{\Sigma})$ and the test statistic $S$ defined by (3.3) is constructed from

$$ \hat{\mathcal{U}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \bar{X}_{i,k}^{\text{alt}} \geq (1 + \varepsilon)(\hat{m}_2)_{i,k} \right\}, $$

$$ \hat{\mathcal{D}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : \bar{X}_{i,k}^{\text{alt}} \leq (1 - \varepsilon)(\hat{m}_2)_{i,k} \right\}, $$

$$ \hat{\mathcal{C}}^{\text{alt}} := \left\{ (i, k) \in \{2, \dots, N\} \times \{2, \dots, K\} : (1 - \varepsilon)(\hat{m}_2)_{i,k} < \bar{X}_{i,k}^{\text{alt}} < (1 + \varepsilon)(\hat{m}_2)_{i,k} \right\}. $$

According to Proposition 4.1, the law of the random vector $(\bar{X}_{i,k}^{\text{alt}})_{2 \leq i \leq N, 2 \leq k \leq K}$ is $\mathcal{N} \left( \hat{\mu}, \frac{1}{\mathcal{P}^{\text{alt}}} \hat{\Gamma} \right)$, where $\hat{\mu}$ and $\hat{\Gamma}$ are obtained from (4.1) and (4.2) with $(m, \Sigma)$ replaced by $(\hat{m}, \hat{\Sigma})$. In particular, the $S$ statistic can be easily simulated to determine numerically the critical region.
Concerning the problem of improved parameter estimation from the ini-
tial dataset and an *alt* dataset, as considered in section 3.2.2, it seems difficult
to construct a Gaussian graphical model able to account jointly for the two
datasets, because of the lack of mechanistic interpretation of the model.

4.2 Mechanistic models

In this section, we apply our *alt* modeling approach to Bayesian networks in
Section 4.2.1 and penalized linear regression models in Section 4.2.2.

4.2.1 Dynamic Bayesian networks

First we describe the model, then we apply our *alt* scheme.

Definition of dynamic Bayesian networks

In Bayesian networks, joint law of genes expression is expressed using
an oriented acyclic graph *G* representing the dependency structure between
genes [Hecker et al., 2009, Chang et al., 2011, Venkat et al., 2017]. In this
section, we consider a particular class of dynamic Bayesian networks, where
the dependency structure between the variables **X** := (*X*<sub>*i*,*k*</sub>)<sub>1≤*i*≤*N*, 1≤*k*<*K*</sub>
is encoded by a static oriented graph *G*<sub>stat</sub> = (*V*<sub>stat</sub>, *E*<sub>stat</sub>) with *V*<sub>stat</sub> :=
{1, ..., *N*} and *E*<sub>stat</sub> the set of oriented edges of the graph, as shown in
Figure 8(a). From this static graph, the dynamic Bayesian network *G* =
(*V*, *E*) with *V* = {1, ..., *N*}×{1, ..., *K*}, where the set of edges *E* is obtained
as shown in Figure 8(b), so that a gene can only influence another gene (and
possibly itself) at the time that immediately follows and that the parents
of a gene do not change over time. Note that the acyclicity condition is
automatically satisfied.

<figure><img src="image_9.png" /><figcaption>Figure 8: Example of a dynamic Bayesian network with 4 genes and 3 times (*N* = 4 and *K* = 3): (a) static graph *G*<sub>stat</sub>; (b) dynamic Bayesian dependency network *G*.</figcaption></figure>
Usually, the values taken by gene expressions are discretized in Bayesian networks. In general, this prevents to apply our method with $\beta > 0$. For $\beta = 0$, we could consider any discretization of gene expressions containing value 0. In the sequel, for simplicity, we restrict to the case where gene expressions can take only two values 0 (for low expression) and 1 (for high expression). In this case, the model corresponds to a probabilistic Boolean network. The probabilistic model corresponding to such a dynamic Bayesian network is given by

$$ \mathbb{P}(\mathbf{X} = \mathbf{x}) = \prod_{k=1}^{K-1} \prod_{i=1}^{N} \mathbb{P}(X_{i,k+1} = x_{i,k+1} | (X_{\ell,k})_{\ell \in A(i)} = (x_{\ell,k})_{\ell \in A(i)}) \\ \times \prod_{i=1}^{N} \mathbb{P}(X_{i,1} = x_{i,1}), \qquad (4.3) $$

where $\mathbf{x} = (x_{i,k})_{1 \le i \le N, 1 \le k \le K} \in \{0, 1\}^{NK}$ and $A(i)$ is the set of parent nodes of $i \in \{1, \dots, N\}$ in the static graph.

The conditional probabilities of the dynamic Bayesian network are assumed to be independent of time. Thus, the parameters of the model are $\theta = (G, p^G)$ where G is the static graph of the model and

$$ p^G := (p_i^G(x_i, (x_\ell)_{\ell \in A(i)}), 1 \le i \le N, x_i \in \{0, 1\}, (x_\ell)_{\ell \in A(i)} \in \{0, 1\}^{\text{Card}(A(i))}), $$
where, for all $i \in \{1, \dots, N\}$,

$$ p_i^G(x_i, (x_\ell)_{\ell \in A(i)}) := \mathbb{P}((X_{i,k+1} = x_i | (X_{\ell,k})_{\ell \in A(i)} = (x_\ell)_{\ell \in A(i)})) $$

for all $k \in \{1, \dots, K-1\}$ and

$$ p_i^G(0, (x_\ell)_{\ell \in A(i)}) + p_i^G(1, (x_\ell)_{\ell \in A(i)}) = 1 $$

for all $(x_\ell)_{\ell \in A(i)} \in \{0, 1\}^{\text{Card}(A(i))}$.

## Description of conditional and mechanistic alt models

For simplicity, we consider the alteration of gene 1 of parameter $\beta = 0$, i.e. the knock-out of gene 1. From the model's distribution (4.3), it is straightforward to construct the alteration model using the conditioning approach of Section 2.2.1, denoted by $(X_{ik}^{\text{cond}})_{1 \le i \le N, 1 \le k \le K}$. The arguments developed to study the example of Section 2.2.3 can be easily extended to show that any dynamic Bayesian network is a mechanistic model. Hence the mechanistic approach of Section 2.2.2 can also be applied to construct a gene alteration model $(X_{ik}^{\text{meca}})_{1 \le i \le N, 1 \le k \le K}$ as in Section 2.2.3.

We denote by $G'_{\text{stat}}$ the graph obtained from $G_{\text{stat}}$ by removing vertex 1 and all edges incident to 1. We call a connected component of $G_{\text{stat}}$ upstream
of gene 1 (resp. downstream of gene 1) a connected component $C$ of $G'_{stat}$ such that

$$ \forall j \in C, \quad (1, j) \notin E_{\text{stat}} \quad (\text{resp. } (j, 1) \notin E_{\text{stat}}). $$

According to the acyclicity assumption of the Bayesian network, a connected component $C$ of $G'_{stat}$ is necessarily either upstream or downstream of gene 1, or both when there is no edge connecting gene 1 to $C$.

**Proposition 4.2.** Let $C$ be a connected component downstream of gene 1. Then

$$ (X_{ik}^{\text{cond}})_{i \in C, 1 \le k \le K} \stackrel{(d)}{=} (X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K} \quad (4.4) $$

and these two vectors have the same law as the random vector $(X_{ik})_{i \in C, 1 \le k \le K}$ conditionally on $(X_{1k})_{1 \le k \le K} = 0$.

In addition, $(X_{ik}^{\text{cond}})_{i \in C, 1 \le k \le K}$ is independent of $(X_{ik}^{\text{cond}})_{i \in \{2, \dots, N\} \setminus C, 1 \le k \le K}$. An identical result holds true for $X^{\text{meca}}$.

**Proposition 4.3.** Let $C$ be a connected component upstream of gene 1. Then

$$ (X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K} \stackrel{(d)}{=} (X_{ik})_{i \in C, 1 \le k \le K} \quad (4.5) $$

and $(X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K}$ is independent of $(X_{ik}^{\text{meca}})_{i \in \{2, \dots, N\} \setminus C, 1 \le k \le K}$.

The proof of these results is given in Appendix B.

We have already given in Section 2.2.3 an example where the law of genes in a connected component upstream of gene 1 is not the same in the conditional and the mechanistic approaches. In particular, the conditional model does not satisfy the relationship (4.5) in general. Similarly, the independence property of Proposition 4.3 is not true in general for the conditional model, as the following example shows.

We consider the gene expression model with $N = 3$ and $K = 2$ where $X_{11}, X_{21}, X_{22}, X_{31}$ and $X_{32}$ are i.i.d. with Bernoulli distribution with parameter $1/2$ and

$$ X_{12} = |X_{21} - X_{31}|. $$

It is easy to verify that this model is a dynamic Bayesian network, the static network of which is shown in Figure 9.

<figure><img src="image_10.png" /><figcaption>Figure 9: Example of a static network where the conditional independence property is not satisfied by the conditional alt model.</figcaption></figure>
The sets {2} are {3} are two connected components upstream gene 1. Since

$$ \{X_{12} = 0\} = \{X_{21} = X_{31} = 1\} \cup \{X_{21} = X_{31} = 0\} $$

then conditionally on $X_{12} = 0$, the random variables $X_{21}$ are $X_{31}$ not independent (they are actually equal). This shows that the independence property in Proposition 4.3 is not satisfied.

### Improved estimation of initial model parameters using alt data

Bayesian networks are usually inferred by a *posteriori* likelihood maximization on graphs. Let us recall the method. Let $(\mathbb{Q}(G), G \in \mathcal{G})$ be an *a priori* distribution over all oriented acyclic graphs and a dataset $\mathbf{x} \in D^{NKP}$ as in (2.1). We maximise the *a posteriori* log-likelihood with respect to $(G, p^G)$:

$$ \mathcal{L}(G, p^G | \mathbf{x}) := \log \mathbb{P}_{(G, p^G)}(\mathbf{X} = \mathbf{x}) + \log \mathbb{Q}(G). $$

Classical choices for the *a priori* distribution $\mathbb{Q}$ are those leading to an AIC or BIC type penalty in order to select networks with few edges [Filkov, 2005]. If an *alt* dataset $\mathbf{x}^{\text{alt}}$ corresponding to a total silencing experiment of gene 1 is available, the method of improvement of the parameter estimation of the initial model proposed in Section 3.2.2 can be applied. Using the previous notations and (3.5), the new *a posteriori* log-likelihood to be maximized is expressed as

$$ \mathcal{L}^{\text{cond}}(G, p^G | \mathbf{x}, \mathbf{x}^{\text{alt}}) := \mathcal{L}(G, p^G | \mathbf{x}) + \log \mathbb{P}_{(G, p^G)}^{\text{alt}}(\mathbf{X}^{\text{cond}} = \mathbf{x}^{\text{alt}}) $$

for the conditional model and

$$ \mathcal{L}^{\text{meca}}(G, p^G | \mathbf{x}, \mathbf{x}^{\text{alt}}) := \mathcal{L}(G, p^G | \mathbf{x}) + \log \mathbb{P}_{(G, p^G)}^{\text{alt}}(\mathbf{X}^{\text{meca}} = \mathbf{x}^{\text{alt}}). $$

for the mechanistic model.

### 4.2.2 Penalized linear regression

We consider the linear, mechanistic, Gaussian model with delay (called LiRE model in [Champagnat et al., 2025]), which generalizes the one introduced in [Schleiss et al., 2021]. This model can be seen either as the discretization of a model based on differential equations or as a regression model.

#### Definition of the LiRE model

In this model, gene expressions satisfy the following inductive relationships:

$$ X_{ik} := \mu_{ik} + \sum_{j=1}^{N} \sum_{\ell=1}^{k-1} a_{\ell} \omega_{ij} X_{j(k-\ell)} + G_{ik} \quad (4.6) $$
for all gene $i \in \{1, \dots, N\}$ and time $k \in \{1, \dots, K\}$, where $\mu_{ik}$ is the average expression of gene $i$ at time $k$, in the absence of influence from other genes, $\omega_{ij}$ is the intrinsic (time-independent) force exerted by gene $j$ on the expression of gene $i$, $a_{\ell}$ is a (gene independent) coefficient that only depends on time differences and $G_{ik}$ stands for the noise.

The coefficients $a_{\ell}$ are assumed to be nonnegative, but $\omega_{ij}$ can be positive (resp. negative) if gene $j$ activates (resp. inhibits) gene $i$ and null if gene $j$ has no influence on gene $i$. The term $a_{\ell}\omega_{ij}$ is the contribution made by gene $j$ at time $k - \ell$ to the expression of gene $i$ at time $k$. This particular form dissociates the effects of time from inter-gene effects. It is justified by the fact that the measurement dates $t_1, \dots, t_K$ are regularly spaced. Replacing the vector $(a_{\ell})_{\ell}$ by a matrix of parameters depending both on the delay $\ell$ and the measurement time $k$ allows to consider the case where the times of measurement are not regularly distributed [Schleiss et al., 2021].

The associated gene network is the directed graph in which an edge points from $j$ to $i$ if and only $\omega_{ij} \neq 0$. Finally, the random variables $(G_{ik})_{1 \le i \le N, 1 \le k \le K}$ are assumed to be Gaussian centered independent, and the variance $\sigma_i^2$ of $G_{ik}$ does not depend on $k$. In particular, the initial expressions $(X_{i1})_{1 \le i \le N} = (\mu_{i1} + G_{i1})_{1 \le i \le N}$ of the genes are independent. We set:

$$ \Sigma := \begin{pmatrix} \sigma_1^2 & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & \sigma_N^2 \end{pmatrix} $$

and

$$ \Lambda := \operatorname{diag}(\underbrace{\Sigma, \dots, \Sigma}_{K \text{ times}}). $$

In order for the model to be identifiable, we assume that $\sum_{\ell=1}^{K-1} a_{\ell} = 1$. The parameter space of the model is then $\Theta = \mathbb{R}^{NK} \times \mathcal{S}_{K-2} \times (\mathbb{R}^{N^2} \setminus \{0\}) \times (\mathbb{R}^*)^N$, where $\mathcal{S}_d$ is the simplex of dimension $d$, and

$$ \theta = ((\mu_{ik})_{1 \le i \le N, 1 \le k \le K}, (a_{\ell})_{1 \le \ell \le K-1}, (\omega_{ij})_{1 \le i,j \le N}, (\sigma_i)_{1 \le i \le N}). $$

is an element of $\Theta$.

Recall that $\mathbb{P}_{\theta}$ represents the distribution of $\mathbf{X} = (X_{ik})_{1 \le i \le N, 1 \le k \le K}$ when the parameter is $\theta \in \Theta$. In the following, we adopt the vector representation

$$ \mathbf{X} = (X_{1,1}, \dots, X_{N,1}, \dots, X_{1,K}, \dots, X_{N,K}). $$

We also define

$$ \mu = (\mu_{1,1}, \dots, \mu_{N,1}, \dots, \mu_{1,K}, \dots, \mu_{N,K}) $$
and H the matrix defined by blocks as follows

$$H := \begin{pmatrix} I & 0 & \dots & 0 & 0 \\ -a_1 \Omega & I & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ -a_{T-2} \Omega & -a_{T-3} \Omega & \dots & I & 0 \\ -a_{T-1} \Omega & -a_{T-2} \Omega & \dots & -a_1 \Omega & I \end{pmatrix},$$

where $I$ denotes the $N \times N$ identity matrix and $\Omega = (\omega_{ij})_{1\le i,j\le N}$. The following result characterizes the probability measure $P_\theta$.

**Theorem 4.4.** For any $\theta \in \Theta$, $P_\theta$ is the multivariate Gaussian distribution

$$P_{\theta} = \mathcal{N}(\Upsilon, \Gamma), \qquad (4.7)$$

where

$$\Upsilon := H^{-1}\mu \quad \text{et} \quad \Gamma := H^{-1}\Lambda(H^{-1})^{\top}. \qquad (4.8)$$

## Mechanistic model of gene alteration

Consider the alteration of gene 1 with parameter $\beta$. The LiRE model is mechanistic, so we can apply the method defined in Section 2.2.2. Let $(X_{ik}^{\text{meca}})_{1\le i\le N, 1\le k\le K}$ be the model obtained:

$$X_{ik}^{\text{meca}} = \mu_{ik} + \beta \sum_{\ell=1}^{k-1} a_{\ell} \omega_{i1} X_{1,k-\ell}^{\text{meca}} + \sum_{j=2}^{N} \sum_{\ell=1}^{k-1} a_{\ell} \omega_{ij} X_{j,k-\ell}^{\text{meca}} + G_{ik} \quad (4.9)$$

for any $i \in \{1, \dots, N\}$ and $k \in \{1, \dots, K\}$.

We recognize a LiRE model where the matrix $\Omega$ has been replaced by $\Omega D_\beta$, where $D_\beta$ is the $NK \times NK$ matrix defined by

$$D_{\beta} = \text{diag}(v), \qquad (4.10)$$

and $v = (\beta, 1, \dots, 1) \in \mathbb{R}^N$. Other parameters $\mu$, $a$ and $\sigma$ remain unchanged. Hence, we deduce from (4.9) that the law $P_\theta^{\text{alt},\beta,\text{meca}}$ of the mechanistic model of alteration of gene 1 with parameter $\beta$ is

$$P_{\theta}^{\text{alt},\beta,\text{meca}} = P_{(a,\Omega D_{\beta},\mu,\sigma)}. \qquad (4.11)$$

Relation (4.11) implies that the gene network associated with the mechanistic alteration model is the same as in the initial model, except when $\beta = 0$ (gene knock-out). In this case, all edges starting from gene 1 are deleted.
## Conditional model of gene alteration

It is convenient to operate with a vector $W = (W^{(1)}, W^{(2)})$ whose coordinates are those of $X$ permuted:

$$W^{(1)} = (X_{11}, \dots, X_{1K}) \quad \text{et} \quad W^{(2)} = (X_{21}, \dots, X_{2K}, X_{31}, \dots, X_{NK}).$$

Alteration concerns only the $W^{(1)}$ component.

According to Theorem 4.4,

$$W \sim \mathcal{N}(\Upsilon_W, \Gamma_W),$$

where

$$\Upsilon_W = (\Upsilon_W^{(1)}, \Upsilon_W^{(2)}) \quad \text{and} \quad \Gamma_W = \begin{bmatrix} \Gamma_W^{(11)} & \Gamma_W^{(12)} \\ \Gamma_W^{(21)} & \Gamma_W^{(22)} \end{bmatrix}$$

are obtained by permuting the rows and columns of $\Upsilon$ et $\Gamma$, $\Upsilon_W^{(1)} \in \mathbb{R}^K$, $\Upsilon_W^{(2)} \in \mathbb{R}^{(N-1)K}$ and $\Gamma_W^{(11)}$ (resp. $\Gamma_W^{(12)}$, $\Gamma_W^{(21)}$ and $\Gamma_W^{(22)}$) is a $K \times K$ (resp. $K \times (N-1)K$, $(N-1)K \times K$ and $(N-1)K \times (N-1)K$) matrix. Let us denote by

$$W^{\text{cond}} := (W^{\text{cond},1}, W^{\text{cond},2})$$

The random vector obtained in the alteration of gene 1 with parameter $\beta$ using the conditional model of Section 2.2.1.

**Theorem 4.5.** Let $\mathbb{P}_{\theta}^{\text{alt},\beta,\text{cond}}$ denote the law of $W^{\text{cond}}$, then

$$\mathbb{P}_{\theta}^{\text{alt},\beta,\text{cond}} = \mathcal{N}(\Upsilon_{W}^{\text{cond}}, \Gamma_{W}^{\text{cond}})$$

where

$$\Upsilon_{W}^{\text{cond}} = \left( \beta \Upsilon_{W}^{(1)}, \Upsilon_{W}^{(2)} - (1 - \beta) \Gamma_{W}^{(21)} (\Gamma_{W}^{(11)})^{-1} \Upsilon_{W}^{(1)} \right) \quad (4.12)$$

and

$$\Gamma_{W}^{\text{cond}} = \begin{pmatrix} \beta^2 \Gamma_{W}^{(11)} & \beta^2 \Gamma_{W}^{(12)} \\ \beta^2 \Gamma_{W}^{(21)} & \Gamma_{W}^{(22)} - (1 - \beta^2) \Gamma_{W}^{(21)} (\Gamma_{W}^{(11)})^{-1} \Gamma_{W}^{(12)} \end{pmatrix}. \quad (4.13)$$

We recover a similar expression as in Proposition 4.1. This result is proven in the more general context of simultaneous alteration of several genes in [Champagnat et al., 2025]. Note that the conditional *alt* model of the last theorem does not seem to be a LiRE model. It is therefore not obvious to associate a gene network with the conditional *alt* model that can be compared with that of the initial LiRE model.
## Improved estimation of initial model parameters using *alt* data

The method proposed in [Champagnat et al., 2025] for estimating the parameters of the LiRE model is based on a posteriori likelihood maximization with the prior density:

$$ \pi(\theta) = \prod_{i=1}^{N} \left( \left( \frac{\tau_i}{2\sigma_i^2} \right)^N \exp \left( -\frac{\tau_i}{\sigma_i^2} \sum_{j=1}^{N} |\omega_{ij}| \right) \right), \quad (4.14) $$

where $\tau = (\tau_1, \dots, \tau_N) \in (\mathbb{R}_+^*)^N$ is a parameter defined by the user. In other words, conditionally on $\sigma$, the $\omega_{ij}$ follow independent Laplace laws. A flat a priori is assumed for the other parameters $\mu$, $a$ and $\sigma$.

Therefore, denoting by $f_\theta$ the density of the distribution of $\mathbf{X}$ under $\mathbb{P}_\theta$, inference of the parameter $\theta$ is performed by maximizing the a posteriori log-likelihood of observations $\mathbf{x}$.

$$ \mathcal{L}(\theta | \mathbf{x}) = \sum_{p=1}^{P} \log f_{\theta}((x_{i,k,p})_{1 \le i \le N, 1 \le k \le K}) + \log \pi(\theta) \quad (4.15) $$

$$ = -\frac{1}{2} \sum_{i=1}^{N} \frac{1}{\sigma_i^2} \left( \sum_{k=1}^{K} \sum_{p=1}^{P} \left( x_{i,k,p} - \mu_{ik} - \sum_{j=1}^{N} \sum_{\ell=1}^{k-1} a_{\ell} \omega_{ij} x_{j,k-\ell,p} \right)^2 + 2\tau_i \sum_{j=1}^{N} |\omega_{ij}| \right) \\ - \left( \frac{PK + 2N}{2} \right) \sum_{i=1}^{N} \ln \sigma_i^2 + N \sum_{i=1}^{N} \ln \frac{\tau_i}{2} $$

The maximization algorithm of $\mathcal{L}(\theta | \mathbf{x})$ used in [Champagnat et al., 2025] consists in successively optimizing on the 4 parameters $\mu$, $\sigma$, $a$ and $\Omega$ while fixing the others. Optimization in $\mu$ and $\sigma$ is explicit, that in $a$ can be expressed as a least-squares problem with a positivity constraint and that in $\Omega$ can be reduced to $N$ lasso problems [Tibshirani, 1996], each involving a vector of the form $(\omega_{i1}, \dots, \omega_{iN})$. The choice of the parameters $\tau_i$ can be made by cross-validation.

If we are given an additional dataset $\mathbf{x}^{\text{alt}}$ of size $P^{\text{alt}}$ from an alteration experiment of gene 1 with parameter $\beta$, we can apply the method in Section 3.2.2. Let $f_{\theta}^{\text{alt},\beta,\text{meca}}$ stand for the density of $\mathbf{X}$ under $\mathbb{P}_{\theta}^{\text{alt},\beta,\text{meca}}$. Using (4.15), the new a posteriori log-likelihood is expressed in the case of the mechanistic approach as

$$ \mathcal{L}^{\text{meca}}(\theta | \mathbf{x}, \mathbf{x}^{\text{alt}}) = \sum_{p=1}^{P} \log f_{\theta}((x_{i,k,p})_{1 \le i \le N, 1 \le k \le K}) \\ + \sum_{p=1}^{P^{\text{alt}}} \log f_{\theta}^{\text{alt},\beta,\text{meca}}((x_{i,k,p}^{\text{alt}})_{1 \le i \le N, 1 \le k \le K}) + \log \pi(\theta). $$
Using (4.11) and Theorem 4.4, we obtain

$$
\begin{align*}
\mathcal{L}^{\text{meca}}(\theta | \mathbf{x}, \mathbf{x}^{\text{alt}}) = & -\frac{1}{2} \sum_{i=1}^{N} \frac{1}{\sigma_i^2} \left( \sum_{k=1}^{K} \sum_{p=1}^{P} \left( x_{i,k,p} - \mu_{ik} - \sum_{j=1}^{N} \sum_{\ell=1}^{k-1} a_{\ell} \omega_{ij} x_{j,k-\ell,p} \right) \right)^2 \\
& + \sum_{k=1}^{K} \sum_{p=1}^{P^{\text{alt}}} \left( x_{i,k,p}^{\text{alt}} - \mu_{ik} - \beta \sum_{\ell=1}^{k-1} a_{\ell} \omega_{i1} x_{1,k-\ell,p}^{\text{alt}} - \sum_{j=2}^{N} \sum_{\ell=1}^{k-1} a_{\ell} \omega_{ij} x_{j,k-\ell,p}^{\text{alt}} \right)^2 \\
& + 2\tau_i \sum_{j=1}^{N} |\omega_{ij}| - \left( \frac{PK + P^{\text{alt}}K + 2N}{2} \right) \sum_{i=1}^{N} \ln \sigma_i^2 + N \sum_{i=1}^{N} \ln \frac{\tau_i}{2}.
\end{align*}
$$

Maximization of this log-likelihood can be achieved using the same itera-
tive method as for inference: explicit optimization in $\mu$ and $\sigma$, optimization
in $a$ by least squares with a positivity constraint and optimization in $\Omega$ by
lasso.

In the case of the conditional approach, the log-likelihood can be com-
puted explicitely according to Theorem 4.5, but its expression lacks the struc-
ture that allows successive optimization in each parameter as above. Hence
likelihood maximization is difficult for this model.

A Conditional laws of Gaussian vectors and proof
of Proposition 4.1

Let $m$ and $d$ denote integers such that $1 \le m < d$. We consider a Gaussian vector $Y \in \mathbb{R}^d$ with mean $\mu$ and covariance matrix $\Gamma$. Write $Y = (Y_1, Y_2)$ where $Y_1 \in \mathbb{R}^m$ is the vector of the first $m$ coordinates of $Y$ and $Y_2 \in \mathbb{R}^{d-m}$ is the vector of the last $d-m$ coordinates. In a similar way, we decompose the vector $\mu$ and the matrix $\Gamma$:

$$
\mu = (\mu_1, \mu_2) \quad \text{et} \quad \Gamma = \begin{pmatrix} \Gamma_{11} & \Gamma_{12} \\ \Gamma_{21} & \Gamma_{22} \end{pmatrix} \qquad (A.1)
$$

where $\mu_1 \in \mathbb{R}^m$, $\mu_2 \in \mathbb{R}^{d-m}$, $\Gamma_{11}$ (resp. $\Gamma_{12}$, $\Gamma_{21}$, $\Gamma_{22}$) is a $m \times m$ (resp. $m \times (d-m)$, $(d-m) \times m$, $(d-m) \times (d-m)$) matrix. The following result describes the conditional distribution of $Y_1$ given $Y_2$ [Eaton, 1983, Prop. 3.13].

**Theorem A.1.** If $\Gamma$ is invertible, then for any $y_2 \in \mathbb{R}^{d-m}$,

$$
\text{Law}(Y_1 | Y_2 = y_2) = \mathcal{N}_m(\mu_{1|2}(y_2), \Gamma_{1|2}),
$$

where

$$
\begin{align*}
\mu_{1|2}(y_2) &= \mu_1 + \Gamma_{12} \Gamma_{22}^{-1} (y_2 - \mu_2) \\
\Gamma_{1|2} &= \Gamma_{11} - \Gamma_{12} \Gamma_{22}^{-1} \Gamma_{21}.
\end{align*}
$$
The case $\beta = 0$ in Point (i) of Proposition 4.1 is an immediate consequence of the last formula (exchanging the role of $Y_1$ and $Y_2$). We leave to the reader the general case.

The proof of Theorem A.1 is based on the calculation of the precision matrix $\Gamma^{-1}$ given by the following proposition.

**Proposition A.2.** If $\Gamma$ is invertible, then

$$ \Gamma^{-1} = \begin{pmatrix} \Gamma_{1|2}^{-1} & -\Gamma_{1|2}^{-1}\Gamma_{12}\Gamma_{22}^{-1} \\ -\Gamma_{22}^{-1}\Gamma_{21}\Gamma_{1|2}^{-1} & \Gamma_{22}^{-1} + \Gamma_{22}^{-1}\Gamma_{21}\Gamma_{1|2}^{-1}\Gamma_{12}\Gamma_{22}^{-1} \end{pmatrix} $$

where $\Gamma_{2|1} = \Gamma_{22} - \Gamma_{21}\Gamma_{11}^{-1}\Gamma_{12}$.

Writing as in (A.1)

$$ \Gamma^{-1} = \begin{pmatrix} (\Gamma^{-1})_{11} & (\Gamma^{-1})_{12} \\ (\Gamma^{-1})_{21} & (\Gamma^{-1})_{22} \end{pmatrix}, $$

we deduce from Proposition A.2 that $(\Gamma^{-1})_{11} = \Gamma_{1|2}^{-1}$, hence Point (ii) of Proposition 4.1 follows.

# B Proof Propositions 4.2 and 4.3

We consider only the case of a dynamic Bayesian model with discrete values (as usual for gene networks). A similar proof can be done in the case where gene expressions are continuous variables.

## B.1 Proof of Proposition 4.2

Throughout, $\mathcal{C}$ denotes a connected component of $G_{\text{stat}}$ downstream of gene 1.

*Step 1. Study of the conditional model*

As defined by the Bayesian model

$$ \begin{align*} \mathbb{P}((X_{ik})_{i \in \mathcal{C}, 1 \le k \le K} &= (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K}, (X_{1k})_{1 \le k \le K} = (x_{1k})_{1 \le k \le K}, \\ (X_{jk})_{j \notin \mathcal{C} \cup \{1\}, 1 \le k \le K} &= (x_{jk})_{j \notin \mathcal{C} \cup \{1\}, 1 \le k \le K} \\ &= \mathbb{P}((X_{ik})_{i \in \mathcal{C}, 1 \le k \le K} = (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K} | (X_{1k})_{1 \le k \le K} = (x_{1k})_{1 \le k \le K}) \\ & \times \mathbb{P}((X_{1k})_{1 \le k \le K} = (x_{1k})_{1 \le k \le K}, (X_{jk})_{j \notin \mathcal{C} \cup \{1\}}, 1 \le k \le K) \end{align*} $$

for all $(x_{1k})_{1 \le k \le K}$, $(x_{ik})_{i \in \mathcal{C}, 1 \le k \le K}$ and $(x_{jk})_{j \notin \mathcal{C} \cup \{1\}, 1 \le k \le K}$. Recall that the law of the conditional model $(X_{ik}^{\text{cond}})_{1 \le i \le N, 1 \le k \le K}$ is the conditional law of
$(X_{ik})_{2 \le i \le N, 1 \le k \le K}$ given $(X_{1k})_{1 \le k \le K} = 0$. From the previous formula, we deduce that

$$
\begin{align*}
\mathbb{P}((X_{ik}^{\text{cond}})_{2 \le i \le N, 1 \le k \le K} = (x_{ik})_{2 \le i \le N, 1 \le k \le K}) \\
&= \mathbb{P}((X_{ik})_{i \in \mathcal{C}, 1 \le k \le K} = (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K} | (X_{1k})_{1 \le k \le K} = 0) \\
&\quad \times \mathbb{P}((X_{jk})_{j \notin \mathcal{C} \cup \{1\}}, 1 \le k \le K) = (x_{jk})_{j \notin \mathcal{C} \cup \{1\}, 1 \le k \le K} | (X_{1k})_{1 \le k \le K} = 0).
\end{align*}
$$

Hence $(X_{ik}^{\text{cond}})_{i \in \mathcal{C}, 1 \le k \le K}$ is independent of $(X_{jk}^{\text{cond}})_{j \notin \mathcal{C} \cup \{1\}, 1 \le k \le K}$ and

$$
\begin{align*}
\mathbb{P}((X_{ik}^{\text{cond}})_{i \in \mathcal{C}, 1 \le k \le K} = (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K}) \\
&= \mathbb{P}((X_{ik})_{i \in \mathcal{C}, 1 \le k \le K} = (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K} | (X_{1k})_{1 \le k \le K} = 0).
\end{align*}
$$

Step 2. Study of the mechanistic model

It is easy to show, as we did in Section 2.2.3, that the dynamic Bayesian
model described in Section 4.2.1 can be written as a mechanistic model
of the form

$$
X_{i,k+1} = F_{i,k}((X_{\ell,k})_{\ell \in A(i)}; U_{i,k}),
$$

for all $i \in \{1, \dots, N\}$ and $k \in \{1, \dots, K-1\}$, where $F_{i,k}((x_{\ell,k})_{\ell \in A(i)}; \cdot)$ is
the generalized inverse distribution function of the conditional law of $X_{i,k+1}$
given $(X_{\ell,k})_{\ell \in A(i)} = (x_{\ell,k})_{\ell \in A(i)}$ and the variables $U_{i,k}$ are i.i.d. uniform on
$[0, 1]$. We deduce that the mechanistic model for total silencing of gene 1 is
given by

$$
\begin{align*}
X_{1,k+1}^{\text{meca}} &= 0, \\
X_{i,1}^{\text{meca}} &= X_{i1}, \tag{B.1}
\end{align*}
$$

$$
X_{i,k+1}^{\text{meca}} = \hat{F}_{i,k}((X_{\ell,k}^{\text{meca}})_{\ell \in A(i) \setminus \{1\}}; U_{i,k}), \quad (\text{B.2})
$$

for all $i \in \{2, \dots, N\}$ and $k \in \{1, \dots, K-1\}$, where

$$
\hat{F}_{i,k}((x_{\ell,k})_{\ell \in A(i) \setminus \{1\}}; u_{ik}) = \begin{cases} F_{i,k}(0, (x_{\ell,k})_{\ell \in A(i) \setminus \{1\}}; u_{ik}) & \text{if } 1 \in A(i) \\ F_{i,k}((x_{\ell,k})_{\ell \in A(i)}; u_{ik}) & \text{if } 1 \notin A(i). \end{cases} \tag{B.3}
$$

Since $A(i) \cap C = \emptyset$ for all $i \notin C \cup \{1\}$, the relations

$$
X_{i,k+1}^{\text{meca}} = \hat{F}_{i,k}((X_{\ell,k}^{\text{meca}})_{\ell \in A(i) \setminus \{1\}}; U_{i,k}) \quad (\text{B.4})
$$

satisfied for any $i \notin C \cup \{1\}$ and $1 \le k \le K - 1$ imply

$$
\sigma((X_{ik}^{\text{meca}})_{i \notin \mathcal{C} \cup \{1\}, 1 \le k \le K}) \subset \sigma((X_{i1})_{i \notin \mathcal{C} \cup \{1\}}, (U_{ik})_{i \notin \mathcal{C} \cup \{1\}, 1 \le k \le K-1}) \quad (\text{B.5})
$$

where $\sigma((Z_j)_{j \in J})$ stands for the $\sigma$-algebra generated of the random variables $(Z_j)_{j \in J}$.
Similarly, since for all $i \in C$ we have $A(i) \subset C \cup \{1\}$ the relations (B.4) with $i \in C$, $1 \le k \le K - 1$ imply

$$ \sigma((X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K}) \subset \sigma((X_{i1})_{i \in C}, (U_{ik})_{i \in C, 1 \le k \le K-1}). \quad (B.6) $$

Since the families of random variables $(U_{ik})_{1 \le i \le N, 1 \le k \le K}$, $(X_{i1})_{i \in C}$ and $(X_{i1})_{i \notin C \cup \{1\}}$ are independent, we deduce that $(X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K}$ is independent of $(X_{jk}^{\text{meca}})_{j \notin C \cup \{1\}, 1 \le k \le K}$.

**Step 3. Conclusion** To complete the proof of Proposition 4.2 it remains to show that

$$ \begin{aligned} \mathbb{P}((X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le K} &= (x_{ik})_{i \in C, 1 \le k \le K}) \\ &= \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le K} = (x_{ik})_{i \in C, 1 \le k \le K} | (X_{1k})_{1 \le k \le K} = 0). \end{aligned} $$

We will establish this identity by showing by induction on $\ell$ belonging to $\{1, \dots, K\}$

$$ \begin{aligned} \mathbb{P}((X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le \ell} &= (x_{ik})_{i \in C, 1 \le k \le \ell}) \\ &= \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell} | (X_{1k})_{1 \le k \le \ell} = 0). \end{aligned} \quad (B.7) $$

The property is easily verified when $\ell = 1$, using (B.1) and the fact that the random variables $(X_{i1})_{1 \le i \le N}$ are independent. Assume that (B.7) is true for $\ell \in \{1, \dots, K-1\}$ and let us prove it for $\ell+1$. Using (B.2), (B.3) and (B.7), we have for all $(x_{ik})_{i \in C, 1 \le k \le \ell+1}$,

$$ \begin{aligned} \mathbb{P}((X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le \ell+1} &= (x_{ik})_{i \in C, 1 \le k \le \ell+1}) \\ &= \mathbb{P}((X_{i,\ell+1}^{\text{meca}})_{i \in C} = (x_{i,\ell+1})_{i \in C} | (X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell}) \\ & \quad \times \mathbb{P}((X_{ik}^{\text{meca}})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell}) \\ &= \mathbb{P}(\hat{F}_{i,\ell}((X_{j,\ell}^{\text{meca}})_{j \in A(i)\setminus\{1\}}; U_{i\ell}) = x_{i,\ell+1}, \forall i \in C | (X_{i\ell}^{\text{meca}})_{i \in C} = (x_{i\ell})_{i \in C}) \\ & \quad \times \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell} | (X_{1k})_{1 \le k \le \ell} = 0) \\ &= \mathbb{P}(\hat{F}_{i,\ell}((x_{j,\ell})_{j \in A(i)\setminus\{1\}}; U_{i\ell}) = x_{i,\ell+1}, \forall i \in C) \\ & \quad \times \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell} | (X_{1k})_{1 \le k \le \ell} = 0) \\ &= \mathbb{P}(X_{i,\ell+1} = x_{i,\ell+1}, \forall i \in C | X_{1\ell} = 0, (X_{i\ell})_{i \in C} = (x_{i\ell})_{i \in C}) \\ & \quad \times \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le \ell} = (x_{ik})_{i \in C, 1 \le k \le \ell} | (X_{1k})_{1 \le k \le \ell} = 0) \\ &= \mathbb{P}((X_{ik})_{i \in C, 1 \le k \le \ell+1} = (x_{ik})_{i \in C, 1 \le k \le \ell+1} | (X_{1k})_{1 \le k \le \ell+1} = 0), \end{aligned} $$

where the fourth equality follows from the definition of $\hat{F}_{i,\ell}$. This concludes the proof of Proposition 4.2.
## B.2 Proof of Proposition 4.3

Let $\mathcal{C}$ be a connected component of $G_{\text{stat}}$, upstream of gene 1. Similarly as in Step 2 of the proof of Section B.1, the $\sigma$-fields inclusions (B.5) and (B.6) still hold true in this case. Hence $(X_{ik}^{\text{meca}})_{i \in \mathcal{C}, 1 \le k \le K}$ is independent of $(X_{ik}^{\text{meca}})_{i \in \{2, \dots, N\} \setminus \mathcal{C}, 1 \le k \le K}$. Furthermore, since, for all $i \in \mathcal{C}$ and $1 \le k \le K$,

$$ \hat{F}_{i,k}((x_{\ell,k})_{\ell \in A(i) \setminus \{1\}}; u_{ik}) = F_{i,k}((x_{\ell,k})_{\ell \in A(i)}; u_{ik}), $$

we deduce

$$ \begin{aligned} \mathbb{P}((X_{ik}^{\text{meca}})_{i \in \mathcal{C}, 1 \le k \le K} &= (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K}) \\ &= \mathbb{P}((X_{ik})_{i \in \mathcal{C}, 1 \le k \le K} = (x_{ik})_{i \in \mathcal{C}, 1 \le k \le K}). \end{aligned} $$

Hence Proposition 4.3 is proved.

**Acknowledgements:** We would like to thank Alvyn Bonnet who, as a Master student, carried out a major bibliographical compilation. RL thanks Cancéropôle Est and Région Grand Est for the funding of his PhD thesis. NC, LV and PV thank INSERM / AVIESAN ITMO Cancer for the funding of the project Predi-CLL: "Quantification et prédiction de l'évolution de l'hétérogénéité clonale pour la leucémie lymphoïde chronique".

**Data availability statement:** No data are associated with this article.

## References

F. Azuaje, Y. Devaux, and D. R. Wagner. Identification of potential targets in biological signalling systems through network perturbation analysis. *Biosystems*, 100(1):55–64, 2010.

A.-L. Barabasi and Z. N. Oltvai. Network biology: understanding the cell's functional organization. *Nature Reviews Genetics*, 5(2):101–113, 2004.

F. Bianconi, E. Baldelli, V. Luovini, E. F. Petricoin, L. Crinò, and P. Valigi. Conditional robustness analysis for fragility discovery and target identification in biochemical networks and in cancer systems biology. *BMC Systems Biology*, 9:1–18, 2015.

M. Boettcher and M. T. McManus. Choosing the right tool for the job: Rnai, talen, or crispr. *Molecular Cell*, 58(4):575–585, 2015.

A. Bonnaffoux, U. Herbach, A. Richard, A. Guillemin, S. Gonin-Giraud, P.-A. Gros, and O. Gandrillon. Wasabi: a dynamic iterative framework for gene regulatory network inference. *BMC Bioinformatics*, 20:1–19, 2019.
M. Bouvier, S. Zreika, E. Vallin, C. Fourneaux, S. Gonin-Giraud, A. Bonnaffoux, and O. Gandrillon. Topodoe: a design of experiment strategy for selection and refinement in ensembles of executable gene regulatory networks. *BMC Bioinformatics*, 25(1):245, 2024.

N. Champagnat, R. Loubaton, L. Vallat, and P. Vallois. LiRE: a mechanistic Gaussian model for gene expression which can be used for prediction of biological intervention experiments. In preparation, 2025.

R. Chang, R. Shoemaker, and W. Wang. A novel knowledge-driven systems biology approach for phenotype prediction upon genetic intervention. *IEEE/ACM Transactions on Computational Biology and Bioinformatics*, 8(5):1170–1182, 2011.

J. Chiquet, G. Rigaill, and M. Sundqvist. A multiattribute gaussian graphical model for inferring multiscale regulatory networks: an application in breast cancer. *Gene Regulatory Networks: Methods and Protocols*, pages 143–160, 2019.

B. L. Davidson and P. B. McCray Jr. Current prospects for rna interference-based therapies. *Nature Reviews Genetics*, 12(5):329–340, 2011.

D. C. de Mello, J. M. Menezes, A. T. F. de Oliveira, M. M. Cristovão, E. T. Kimura, and C. S. Fuziwara. Modulating gene expression as a strategy to investigate thyroid cancer biology. *Archives of Endocrinology and Metabolism*, 68(spe 1):e240073, 2024.

R. Dehghannasiri, B.-J. Yoon, and E. R. Dougherty. Efficient experimental design for uncertainty reduction in gene regulatory networks. In *BMC Bioinformatics*, volume 16, pages 1–18. Springer, 2015.

J. A. Doudna and E. Charpentier. The new frontier of genome engineering with crispr-cas9. *Science*, 346(6213):1258096, 2014.

M. L. Eaton. *Multivariate statistics. A vector space approach*. Wiley Ser. Probab. Math. Stat. John Wiley & Sons, Hoboken, NJ, 1983.

V. Filkov. Identifying gene regulatory networks from gene expression data. *Chapman & Hall/CRC Computer & Information Science Series*, 2005.

A. Fire, S. Xu, M. K. Montgomery, S. A. Kostas, S. E. Driver, and C. C. Mello. Potent and specific genetic interference by double-stranded rna in caenorhabditis elegans. *Nature*, 391(6669):806–811, 1998.

J. Goutsias and N. Lee. Computational and experimental approaches for modeling gene regulatory networks. *Current Pharmaceutical Design*, 13 (14):1415–1436, 2007.
M. Hecker, S. Lambeck, S. Toepfer, E. Van Someren, and R. Guthke. Gene regulatory network inference: data integration in dynamic models – a review. *Biosystems*, 96(1):86–103, 2009.

V. A. Huynh-Thu and P. Geurts. dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data. *Scientific Reports*, 8(1):3384, 2018.

V. A. Huynh-Thu and G. Sanguinetti. Gene regulatory network inference: an introductory survey. *Gene Regulatory Networks: Methods and Protocols*, pages 1–23, 2019.

T. E. Ideker, V. Thorsson, and R. M. Karp. Discovery of regulatory interactions through perturbation: inference and experimental design. In *Biocomputing 2000*, pages 305–316. World Scientific, 1999.

M. Ishikawa, S. Sugino, Y. Masuda, Y. Tarumoto, Y. Seto, N. Taniyama, F. Wagai, Y. Yamauchi, Y. Kojima, H. Kiryu, et al. Renge infers gene regulatory networks using time-series single-cell rna-seq data with crispr perturbations. *Communications Biology*, 6(1):1290, 2023.

M. G. Kendall and A. Stuart. *The advanced theory of statistics. Vol. 2: Inference and relationship*. 3rd ed. London: Charles Griffin & Co., Ltd. 723 p., 1973.

R. C. Lee, R. L. Feinbaum, and V. Ambros. The c. elegans heterochronic gene lin-4 encodes small rnas with antisense complementarity to lin-14. *Cell*, 75(5):843–854, 1993.

F. Markowetz, J. Bloch, and R. Spang. Non-transcriptional pathway features reconstructed from secondary effects of rna interference. *Bioinformatics*, 21(21):4026–4032, 2005.

D. H. Martínez-Puente, J. J. Pérez-Trujillo, L. M. Zavala-Flores, A. García-García, A. Villanueva-Olivo, H. Rodríguez-Rocha, J. Valdés, O. Saucedo-Cárdenas, R. Montes de Oca-Luna, and M. d. J. Loera-Arias. Plasmid dna for therapeutic applications in cancer. *Pharmaceutics*, 14(9):1861, 2022.

S. E. Mohr, J. A. Smith, C. E. Shamu, R. A. Neumüller, and N. Perrimon. Rnai screening comes of age: improved techniques and complementary approaches. *Nature Reviews Molecular Cell Biology*, 15(9):591–600, 2014.

J. Natale, D. Hofmann, D. Hernández, and I. Nemenman. *Quantitative Biology: Theory, Computational Methods, and Models*, chapter Reverse-engineering biological networks from large data sets. The MIT Press, 2018.

N. Omranian, J. M. Eloundou-Mbebi, B. Mueller-Roeber, and Z. Nikoloski. Gene regulatory network inference using fused lasso on multiple data sets. *Scientific Reports*, 6(1):20533, 2016.
J. Pearl. *Causal inference in statistics: a primer*. John Wiley & Sons, 2016.

J. Qin, Y. Hu, F. Xu, H. K. Yalamanchili, and J. Wang. Inferring gene regulatory networks by integrating chip-seq/chip and transcriptome data via lasso-type regularization methods. *Methods*, 67(3):294–303, 2014.

D. Sarmah, G. R. Smith, M. Bouhaddou, A. D. Stern, J. Erskine, and M. R. Birtwistle. Network inference from perturbation time course data. *npj Systems Biology and Applications*, 8(1):42, 2022.

C. Schleiss, R. Carapito, L.-M. Fornecker, L. Muller, N. Paul, O. Tahar, A. Pichot, M. Tavian, A. Nicolae, L. Miguet, et al. Temporal multiomic modeling reveals a b-cell receptor proliferative program in chronic lymphocytic leukemia. *Leukemia*, 35(5):1463–1474, 2021.

D. Seçilmiş, T. Hillerton, A. Tjärnberg, S. Nelander, T. E. Nordling, and E. L. Sonnhammer. Knowledge of the perturbation design is essential for accurate gene regulatory network inference. *Scientific Reports*, 12(1): 16531, 2022.

J. Tegner, M. S. Yeung, J. Hasty, and J. J. Collins. Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling. *Proceedings of the National Academy of Sciences*, 100(10):5944–5949, 2003.

R. Tibshirani. Regression Shrinkage and Selection via the Lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, 58(1):267–288, 1996. ISSN 00359246. URL http://www.jstor.org/stable/2346178. Publisher: [Royal Statistical Society, Wiley].

S. M. Ud-Dean and R. Gunawan. Ensemble inference and inferability of gene regulatory networks. *PloS One*, 9(8):e103812, 2014.

L. Vallat, C. A. Kemper, N. Jung, M. Maumy-Bertrand, F. Bertrand, N. Meyer, A. Pocheville, J. W. Fisher III, J. G. Gribben, and S. Bahram. Reverse-engineering the genetic circuitry of a cancer cell with predicted intervention in chronic lymphocytic leukemia. *Proceedings of the National Academy of Sciences*, 110(2):459–464, 2013.

P. S. Venkat, K. R. Narayanan, and A. Datta. A bayesian network-based approach to selection of intervention points in the mitogen-activated protein kinase plant defense response pathway. *Journal of Computational Biology*, 24(4):327–339, 2017.

E. Ventre, U. Herbach, T. Espinasse, G. Benoit, and O. Gandrillon. One model fits all: combining inference and simulation of gene regulatory networks. *PLoS Computational Biology*, 19(3):e1010962, 2023.
S. Vercauteren, S. Fiesack, L. Maroc, N. Verstraeten, L. Dewachter, J. Michiels, and S. C. Vonesch. The rise and future of crispr-based approaches for high-throughput genomics. *FEMS Microbiology Reviews*, 48 (5):fuae020, 2024.

N. Villa-Vialaneix, M. Vignes, N. Viguerie, and M. S. Cristobal. Inferring networks from multiple samples with consensus lasso. *Quality Technology & Quantitative Management*, 11(1):39–60, 2014.

H. Vundavilli, A. Datta, C. Sima, J. Hua, R. Lopes, and M. Bittner. Bayesian inference identifies combination therapeutic targets in breast cancer. *IEEE Transactions on Biomedical Engineering*, 66(9):2684–2692, 2019.

A. Wagner. How to reconstruct a large genetic network from n gene perturbations in fewer than $n^2$ easy steps. *Bioinformatics*, 17(12):1183–1197, 2001.

Y. R. Wang and H. Huang. Review on statistical methods for gene network reconstruction using expression data. *Journal of Theoretical Biology*, 362: 53–61, 2014.

J. A. Wolff, R. W. Malone, P. Williams, W. Chong, G. Acsadi, A. Jani, and P. L. Felgner. Direct gene transfer into mouse muscle in vivo. *Science*, 247 (4949):1465–1468, 1990.

P. Zhu, H. M. Aliabadi, H. Uludağ, and J. Han. Identification of potential drug targets in cancer signaling pathways using stochastic logical models. *Scientific Reports*, 6(1):23078, 2016.