(060.3735) Fiber Bragg gratings; (050.5080) Phase shift; (060.2340) Fiber
optics components; (120.2440) Filters; (230.1480) Bragg reflectors

# Aperiodic phase masks for inscribing complex multi-notch OH-emission filters
for astronomy

Kalaga Madhav  innoFSPEC, Leibniz-Institut für Astrophysik Potsdam, An der
Sternwarte 16, Potsdam, Germany, 14482  Corresponding author<EMAIL_ADDRESS>Ziyang Zhang  innoFSPEC, Leibniz-Institut für Astrophysik Potsdam, An der
Sternwarte 16, Potsdam, Germany, 14482  Martin M Roth  innoFSPEC, Leibniz-
Institut für Astrophysik Potsdam, An der Sternwarte 16, Potsdam, Germany,
14482

###### Abstract

We demonstrate for the first time, a new type of aperiodic phase mask (APM)
for fabricating multi-channel aperiodic fiber Bragg gratings. The mask is made
of individual diffraction phase gratings with discrete unequal phase-steps
incorporated at periodic locations. The diffraction at the discrete phase-
steps in the phase mask produces corresponding half phase-steps at periodic
locations along the fiber. The accumulated phase, along with index modulation,
generates the desired multinotch reflection spectrum. Complex fiber Bragg
grating filters fabricated using APM, in a standard phase mask based
fabrication setup, can be used to simultaneously suppress multiple aperiodic
OH emission wavelengths in near infrared (NIR) existing in upper atmosphere,
and increase the sensitivity of ground based telescopes.

††journal: ol

##  1 Introduction

Observation at near infrared wavelengths (NIR) between 0.9 to 2.5$\mu$m are
critical for modern astrophysics, as they provide access to objects heavily
obscured by dust extinction, e.g. the supermassive black hole at the galactic
center, and star forming regions, to cool objects like AGB stars, to the high
redshift universe, etc. – to name but a few. Furthermore, the availability of
high sensitivity large format image sensors and the advent of adaptive optics
have made the NIR an extremely attractive wavelength range such that the new
generation of large ground-based telescopes like the ELT, TMT, or GMT must be
considered mainly NIR facilities. However, observations of faint objects in
the NIR from the ground are overwhelmed by a sky background emission line
spectrum that is typically 1000 times brighter than the NIR light from the
objects of interest. The emission occurs due to de-excitation of atmospheric
hydroxyl (OH) molecules in a cold layer of 6-10 km thickness at altitudes of
90km . At the central wavelengths of these emission lines, the signal of faint
objects is heavily affected by the photon shot noise and strong residuals
associated with the OH lines, so no reliable data can be recorded. Instead,
one has to resort to the interline continuum, a technique also known as OH
avoidance [1]. However, even when resorting to observations between the bright
OH lines, it has been discovered that the faint extended wings, that are due
to scattered light occurring inevitably within the spectrograph, are still
bright enough to affect the detection limit of faint objects in the continuum.
Therefore, it has been considered to filter out the OH lines at high
dispersion, however, first concepts have in reality not provided convincing
results (e.g. [2, 3, 4]). A radically new idea was proposed by [5] which
consists of a filter placed in front of the optical system before the light
enters the spectrograph, thus giving nowhere an opportunity to create
scattered light [6, 8]. In order to suppress or filter out the OH emission
lines, an optical filter capable of $>$30dB suppression at the emission lines
with bandwidths as small as 150nm and, with $>$0.5dB throughput outside and
between the lines will be required. Fiber Bragg gratings (FBGs) are ideal
candidates for filtering with the tight constraints and an aperiodic FBG
(AFBG) filter capable of suppressing $\sim$100 lines has been previously
demonstrated in the GNOSIS experiment [7, 9].

However, fabricating such filters, with good reproducibility is not a trivial
task and requires accurate control of the intensity, phase and exposure length
of a complex interference pattern over a moving photosensitive optical fiber.
Simple or complex gratings can be fabricated through point-by-point (PbP) [10]
or line-by-line (LbL) [11] inscription process using femtosecond lasers and
de-phasing methods [12]. Ultra-long gratings have been fabricated using
electro-optic modulators (EOMs) [13] in push-pull configuration, and complex
OH filters have been fabricated using acousto-optic modulators (AOMs) [14]. E-
or A-OMS fabrication techniques generate a “running-interference”pattern,
similar to a rack-and-pinion, that is synchronised with the velocity of
optical fiber, requiring precise control on the intensity, focus spot size,
and velocity over a long length in real time. Phase mask offers a convenience
that previously mentioned methods do not offer. By transferring the complexity
of fabrication, such as in femtosecond, EOM or AOM techniques, to the one-time
manufactured complexity in the phase mask, the convenience is preserved.
Requiring no moving parts, or stringent alignment, the complex phase mask can
be used off-the-shelf in a standard UV based FBG inscription setup. In this
paper, we introduce for the first time the design of an aperiodic phase mask
to inscribe multi-channel aperiodic filters in hydrogenated or doped
photosensitive fibers, in order to suppress the night-sky OH emission.

##  2 Aperiodic Bragg grating

The index modulation $\Delta n_{g}$ and phase $\phi_{g}$, of the complex
grating can be reconstructed from the desired reflection spectrum $|r|$, by
using layer peeling method described in [21, 22]. For the design of APM, we
selected OH sky emission lines in H-band ranging from 1400nm to 1700nm [23].
For compatibility with future tests on PRAXIS [24] system, that uses the
existing GNOSIS filters, the full-width half maximum (FWHM), transmission, and
wavelength of OH lines are selected as given in Tables.(1,2) in [8]. The
aperiodic filter is defined by [15],

$\begin{multlined}|r\big{(}\lambda\big{)}|=\sqrt{R_{i}}\sum_{i=1}^{N}{\exp\Bigg{[}-\Bigg{\\{}\frac{2\pi
n_{eff}}{p_{i}}\Big{(}\frac{1}{\lambda}-\frac{1}{\lambda_{i}}\Big{)}\Bigg{\\}}^{q_{i}}\Bigg{]}}\times\\\
\exp\Bigg{[}i2\pi
n_{eff}\Bigg{(}\frac{1}{\lambda}-\frac{1}{\lambda_{0}}\Bigg{)}g_{i}\Bigg{]}\end{multlined}|r\big{(}\lambda\big{)}|=\sqrt{R_{i}}\sum_{i=1}^{N}{\exp\Bigg{[}-\Bigg{\\{}\frac{2\pi
n_{eff}}{p_{i}}\Big{(}\frac{1}{\lambda}-\frac{1}{\lambda_{i}}\Big{)}\Bigg{\\}}^{q_{i}}\Bigg{]}}\times\\\
\exp\Bigg{[}i2\pi
n_{eff}\Bigg{(}\frac{1}{\lambda}-\frac{1}{\lambda_{0}}\Bigg{)}g_{i}\Bigg{]}$
(1)

where, $R_{i}$ is the desired reflectivity of individual OH-emission line
filter, $N$ is the number of filters, $(p_{i},q_{i})$ define the shape of the
individual filters, $n_{eff}$ is the effective index, $\lambda_{0}$ is the
seed grating, which also defines the APM’s pitch $\Lambda_{m}$, and $g_{i}$ is
the individual channel’s group delay. Since there exists an upper limit to the
index change achievable in a fiber, $g_{i}$ can be optimised to reduce the
maximum index modulation required. The right choice of $g_{i}$, or “de-
phasing”, [12, 15] effectively spreads the individual gratings over the length
of the grating, instead of crowding them in the same location spatially. We
now describe the steps required to design a mask that can be used in a
standard UV inscription setup to fabricate a grating with a reflection
spectrum defined by eq.(1).

##  3 Design of aperiodic phase mask

It is well known [16, 17, 18] that when a $\phi_{m}$-shifted phase mask is
used in side-writing technique for inscribing FBGs, the phase shift in the
mask is split into two half-phase shifts ($\frac{\phi_{m}}{2}$), separated by
$\Delta z=2y\tan\theta$ in the fiber, as shown in Fig.1, where, $y$ is the
distance between the mask and the fiber core, and $\theta$ is the angle of
diffraction of $\pm 1$ order.

Figure 1: Schematic showing the propagation of phase mask phase to fiber
grating phase.

The cumulative phase of the light propagating through the fiber gives the
desired phase $\phi_{g}$ in the grating. For example, if a $\pi$-shifted phase
mask is used to fabricate the grating, along the grating length there will be
two locations with $\phi_{g}=\pi/2$ phase. In the transmission spectrum, we
would see the characteristic single narrow transmission window at the filter
center. With increasing $\phi_{m}$, the narrow transmission window within the
Bragg stopband shifts to longer wavelength [19]. To achieve the desired phase
$\phi_{g}$ in the grating, the width of the groove, or phase-step,
$\delta_{m}$, in the phase mask is given by,

$\delta_{m}=\frac{\Lambda_{m}}{4\pi}\Big{(}2\pi+\phi_{g}\Big{)}$ (2)

For example, if $\phi_{g}=-\pi$ for the standard $\pi-$shifted FBG, we will
require a phase mask with $\delta_{m}=\frac{\Lambda_{m}}{4}$ at the center of
the mask of length $L$. If we use this mask for fabrication, in the fiber the
phase will be split, $\phi_{g}=\big{(}-\pi/2,-\pi/2\big{)}$. By tuning
$\delta_{m}$, or equivalently $\phi_{m}$, we can inscribe a desired $\phi_{g}$
in the grating at selected locations. Also, when $\delta_{m}=\Lambda_{m}/2$,
then we get $\phi_{g}=0$, and a uniform phase mask. Introducing multiple phase
shifts using high precision PZT translation was proposed in [20] for
fabricating periodic or non-periodic high-channel count FBG.

In order to design a mask that can generate a grating with the desired
aperiodic reflection spectrum defined by eq.(1), we will require the grating’s
phase ($\phi_{g}$). We use layer peeling (LP) technique to first derive the
grating’s complex coupling coefficient, $\kappa$, from which we can extract
$\phi_{g}$. By knowing $\phi_{g}$, we can then calculate the phase-steps
$\delta_{m}$ of the APM using eq.(2).

Figure 2: Representative 2D model of a section of the APM showing two
$\Lambda_{m}/4$ shifted grooves separated by $\Delta z$ corresponding to
($-\pi$, $+\pi$) phase. $\rho=\lambda_{uv}/2(n_{uv}-1)$ is the groove depth of
the phase mask, defined by the wavelength $\lambda_{uv}$ of the laser used for
fabrication, and the refractive index $n_{uv}$ of the mask material.

##  4 Simulation and discussion

Fig.3 shows the FDTD simulation for phase steps at two locations on a phase
mask of $60\mu m$ length, separated by 20$\mu m$, where mask period
$\Lambda_{m}=1.064\mu m$. The two phase steps
$\delta_{m}=\Big{\\{}\frac{\Lambda_{m}}{4},\frac{3\Lambda_{m}}{4}\Big{\\}}$
corresponding to mask phases $\phi_{m}=\big{\\{}-\pi,+\pi\big{\\}}$,
respectively, result in four half-phase regions in the fiber.

Figure 3: Near field of phase mask with two phase steps of $-\pi$ and $+\pi$.
Figure 4: Transmission spectrum, index modulation, phase and APM groove width
for fabricating an aperiodic filter.

For designing the APM specifically made to fabricate an aperiodic grating, we
first define the filter characteristics, such as, transmission, FWHM and
channel center. We chose $N=37$ OH-emission lines between 1500nm to 1600nm,
and $L=47.9$mm. Using LP, and the desired filter spectrum constructed using
eq.(1), we find the index modulation $\Delta n$ and $\phi_{g}$, as shown in
Fig.4. For covering the filter’s bandwidth ($\beta$), we will require a
$\phi_{g}$, or $\delta_{m}$ discretization interval or layer thickness in
layer peeling [22], $\Delta z=\pi/\beta=9.58\mu m$. We use eq.(2) to calculate
the APM’s groove width $\delta_{m}$ from $\phi_{g}$. Since
$\phi_{g}=[-\pi,+\pi]$ rad, we have $\delta_{m}=[266.32,798.96]$ nm.

$\delta_{m}$ can be incorporated in the mask as a nonlinear chirp, where the
groove width $\delta_{m}$, varies continuously over the length of the mask. To
achieve accurate continuously varying $\delta_{m}$ using e-beam process would
be challenging. Alternatively, we design the mask with a global mask pitch of
$\Lambda_{m}=1065.28$nm, corresponding to the seed grating
$\lambda_{0}=1550$nm, and at discrete locations at intervals of $\Delta z$
along the seed phase mask, we incorporate grooves of width $\delta_{m}$
defined by eq.(2). An example of a mask segment with two $\delta_{m}$ is shown
in Fig.2.

##  5 Conclusion

We have shown the steps involved in transferring the spatial structure of an
aperiodic fiber Bragg grating to the corresponding structure in an aperiodic
phase mask. Fabrication of the mask with the desired groove accuracy at
periodic intervals, or as a continuous chirp is not a trivial task. However
with recent advances in e-beam processes, the accuracy required to reproduce
$\delta_{m}$ is a reality. Fabrication of APMs based on the method described
in this paper is currently ongoing. With APMs, the complexity of alignment in
fabrication setups such as EOM, AOM or femtosecond laser, is now transferred
to the fabricated complexity in the mask, facilitating the use of standard
phase mask fabrication to inscribe complex gratings.

## Acknowledgements

This work is supported by the BMBF project “Meta-ZiK Astrooptics” (grant no.
03Z22A511).

## References

  * [1] Martini, P., & DePoy, D. L. 2000, SPIE, 4008, 695
  * [2] Piche, F., Parry, I. R., Ennico, K., et al. 1997, SPIE, 2871, 1332
  * [3] Ennico, K. A., Parry, I. R., Kenworthy, M. A., et al. 1998, SPIE, 3354, 668
  * [4] Maihara, T., Ohta, K., Tamura, N., et al. 2000, SPIE, 4008, 1111
  * [5] Bland-Hawthorn, J., Ellis, S. C., Leon-Saval, S. G., et al. 2011, Nature Communications, 2, 581
  * [6] Ellis, S. C., & Bland-Hawthorn, J. 2008, MNRAS, 386, 47
  * [7] Ellis, S. C., Bland-Hawthorn, J., Lawrence, J., et al. 2012, MNRAS, 425, 1682
  * [8] Trinh, C. Q., Ellis, S. C., Bland-Hawthorn, J., et al. 2013-1, AJ, 145, 51
  * [9] Trinh, C. Q., Ellis, S. C., Bland-Hawthorn, J., et al. 2013-2, MNRAS, 432, 3262
  * [10] Lai, Y., Zhou, K., Sugden, K., Bennion, I. 2007, Opt.Express, 15, 26, 18318
  * [11] Zhou, K., Dubov, M., Mou, C., Zhang, L., Mezentsev, V. K., Bennion, I., 2010 IEEE Photo. Tech. Lett., 22, 16, 1190
  * [12] Buryak, A. V., Kolossovski, K. Y., Stepanov, D. Y. 2003, IEEE J. of Quant. Elect., 39, 1, 91
  * [13] Loranger, S., Lambin-Iezzi, V., Kashyap, R. 2017, Optica 4, 1143
  * [14] Gbadebo, A. A., Turitsyna, E. G., Williams, J. A. R., 2018, Optics Exp., 26, 2, 1315
  * [15] Cao, H., Atai, J., Shu, X., Chen, G., 2012, Opt. Express, 20, 11, 12095
  * [16] Sheng, Y., Rothenberg, J. E., Li, H., Wang, Y., Zweiback, J. 2004, IEEE Photo. Tech. Lett., 16, 5, 1316
  * [17] Tremblay, G., Sheng, Y. 2006, J. Opt. Soc. Am. B, 23, 8, 1511
  * [18] Sheng, Y., and Sun, L. 2005 Opt. Ex., 13, 16, 6111
  * [19] Janos, M., Canning, J., Sceats, M. G. 1996, IEEE Electron. Lett., 32, 3, 245
  * [20] Dai, Y., Yao, J., 2009, IEEE J. of Quant. Elect., 45, 8, 964
  * [21] Skaar, J., Wang, L., Erdogan, T. 2001, IEEE J. of Quant. Elect., 37, 2, 165
  * [22] Skaar, J., Feced, R. 2002 J. Opt. Soc. Am. A , 19, 11, 2229
  * [23] Rousselot, P., Lidman, C., Cuby, J.-G., Moreels, G., Monnet, G. 2000, Astron. Astrophys. 354, 1134
  * [24] Horton, A., et. al., 2012, SPIE 8450, 84501V


# Identification of changes in gene expression

L. Ameis  Institute of Medical Statistics and Computational Biology (IMSB),
Faculty of Medicine, University of Cologne  K. Möllenhoff

###### Abstract

Evaluating the change in gene expression is a common goal in many research
areas, such as in toxicological studies as well as in clinical trials. In
practice, the analysis is often based on multiple t-tests evaluated at the
observed time points. This severely limits the accuracy of determining the
time points at which the gene changes in expression. Even if a parametric
approach is chosen, the analysis is often restricted to identifying the onset
of an effect. In this paper, we propose a parametric method to identify the
time frame where the gene expression significantly changes. This is achieved
by fitting a parametric model to the time-response data and constructing a
confidence band for its first derivative. The confidence band is derived by a
flexible two step bootstrap approach, which can be applied to a wide variety
of possible curves. Our method focuses on the first derivative, since it
provides an easy to compute and reliable measure for the change in response.
It is summarised in terms of a hypothesis test, such that rejecting the null
hypothesis means detecting a significant change in gene expression.
Furthermore, a method for calculating confidence intervals for time points of
interest (e.g. the beginning and end of significant change) is developed. We
demonstrate the validity of our approach through a simulation study and
present a variety of different applications to mouse gene expression data from
a study investigating the effect of a Western diet on the progression of non-
alcoholic fatty liver disease.

Keywords and Phrases: bootstrap, model-based hypothesis testing, time-response
models, gene expression data

* Corresponding author: Kathrin Möllenhoff, eMail<EMAIL_ADDRESS>
##  1 Introduction

One common objective in the field of RNA-sequencing, based on the measurement
of transcriptomes (see Wang et al. [2009]), is to identify differentially
expressed genes. Possible examples include dose-gene expression data (see,
e.g., Krug et al. [2013]) or time-gene expression data (see, e.g., Ghallab et
al. [2021]). In recent times, there have been calls among researchers that a
similar large-scale project on a par with the Human Genome Project (see
Consortium [2001b] and National Human Genome Research Institute ) for RNA is
required (see Pennisi [2024] and National Academies of Sciences, Engineering,
and Medicine and others [2024]). It is therefore of great importance to
develop new methods for the analysis of RNA sequencing data.

Anders and Huber [2010] introduced an approach to analysing gene count data
using the negative binomial distribution, that is employed in the well-known
DESeq R package. Its successor – the DESeq2 packages (see Love et al. [2014])
– has been cited over $35,000$ times as of July 2024 (see NCBI [1996]).
However, as described in the vignette of DESeq2 (see Love et al. [2024]), the
analysis of the data includes the dose/time as a discrete covariate, which
limits the precision of the analysis to the measured dose levels/time points.

One question, where this problem is particularly evident, is the
identification of an alert concentration, defined as the lowest measured
concentration where the difference between the concentration and the control
significantly exceeds a critical level of the relevant response (see
Delignette-Muller et al. [2011]) or the lowest concentration with a noticeable
effect (see Jensen et al. [2019]). Traditionally, the analysis is based on
applying multiple $t$-tests or – as recommended – Dunnett tests (see Hothorn
[2014]), focussing on the measured dose levels. Consequently, with these tests
it is not possible to identify a dosage between those levels. In order to
address this issue, Kappenberg et al. [2021] proposed a model-based approach
assuming a monotonic relationship between dose and response. This method
allows for a continuous identification of the lowest concentration at which
the response significantly exceeds the response for the control by a given
threshold. This method was further developed by Möllenhoff et al. [2022]
through the introduction of a procedure that is independent of the assumption
of monotonicity and replaces an asymptotic calculation of the variance by a
parametric bootstrap.

When analysing time-response data, a similar problem is given by the
identification of the time frame of significant changes in gene expression.
For example, a change in gene expression is indicative of deregulation and
thus the presence of a biological process. In a two-group situation, a
difference in the time periods can indicate the deceleration or acceleration
of an underlying process. Furthermore, a gene ontology (GO) analysis of genes
that undergo a change in expression over a specified time frame allows
inference regarding the processes that occur at the organism level. In Ghallab
et al. [2021] so called rest-and-jump genes (RJG) were observed. These are
genes that showed delayed deregulation, either reaching a plateau afterwards
or continuing to deregulate. In this context, it is of interest to consider
the entire time frame during which the gene expression changes, rather than
focusing only on the beginning of an effect.

Based on this idea, in this paper we propose a model-based approach that
allows the identification of the entire time frame within which the change in
gene expression is considered significant. This is accomplished by formulating
a hypothesis test considering a test statistic based on the absolute value of
the first derivative of a selected continuous parametric model with regard to
a threshold $\lambda\geq 0$. Rejecting the null hypothesis indicates the
presence of a slope, that is significantly larger than $\lambda$ and thereby a
significant change in gene expression. The test decision is based on the
construction of a lower simultaneous confidence band for the first derivative,
inspired by the method developed in Möllenhoff et al. [2022]. Within the
period of significant change in gene expression, it may be of interest to
focus on specific time points, such as the beginning and end of the
significant change. The second goal of this paper is to construct confidence
intervals for these time points.

This paper is structured as follows: First, we introduce the method to
estimate the time frame of significant change in gene expression, based on
calculating the first derivative of the estimated time-response model and the
estimation of a corresponding lower simultaneous confidence band. Based on
this, the method is extended, allowing for the estimation of confidence
intervals for time points of interest. Second, the methodology is subjected to
a simulation study. Finally, a number of potential applications are
illustrated using a data set from the western diet mice study (see Ghallab et
al. [2021]), a toxicological study that examines the impact of a western diet
(WD) on the progression of non-alcoholic fatty liver disease (NAFLD).

##  2 Methodology

The fundamental idea of this method is to initially fit a model to the data
and subsequently assess the change in gene expression by considering its first
derivative. This is possible since the first derivative describes the slope of
the gene expression curve at any time point $t$, and thus can be regarded as a
measure for the change. Once a suitable model has been fitted, the next step
is to develop a hypothesis test that tests if and when the first derivative
significantly exceeds a chosen threshold, denoted $\lambda$. The corresponding
time point marks a change in gene expression. The first three steps of the
procedure are inspired by the method developed in Möllenhoff et al. [2022]
adjusted for the first derivative, while the fourth and last step provides the
possibility to obtain a confidence interval for the identified time points.

###  2.1 Model fit and derivation

We regard a data set with $n$ observations in total, that are conducted at $m$
different time points $t_{p}$, $p=1,\dots,m$. Let $n_{p}$ denote the number of
observations at time point $t_{p}$, i.e. $\sum_{p=1}^{m}n_{p}=n$. Furthermore,
the study period is defined as $T=\left[t_{1},t_{m}\right]$, while the study
duration is expressed as $\overline{t}=t_{m}-t_{1}$. The initial step is to
fit a parametric model to the data. Hereby, a suitable model can be selected
from a pool of possible models using a selection criteria, as, for example,
the AIC (see Fahrmeir et al. [2022]). We define

$y_{p,q}=f(t_{p},\theta)+\varepsilon_{p,q}\ \ \ p=1,\dots,m,\ \ \
q=1,\dots,n_{p},$

where $f$ denotes the selected model. We assume that the errors are
independently and identically normally distributed $\varepsilon_{p,q}\sim
N(0,\sigma^{2})$. The parameters $\theta\in B\subset\mathbb{R}^{r}$ are
estimated by applying the OLS method, i.e. $\hat{\theta}=\arg\min_{b\in
B}\sum_{p=1}^{m}\sum_{q=1}^{n_{p}}\left(y_{p,q}-f(t_{p},b)\right)^{2}$, and,
due to the assumption of normally distributed errors, coincide with the
corresponding ML estimates, so that the AIC can easily be calculated. The
corresponding variance estimate is given by
$\hat{\sigma}^{2}=\frac{1}{n-r}\sum_{p=1}^{m}\sum_{q=1}^{n_{p}}\\{y_{p,q}-f(t_{p},\hat{\theta})\\}^{2}.$
However, as we are interested in detecting significant changes, the method
developed here does not focus on the model $f$ itself, but rather on its first
derivative $f^{\prime}$. Throughout this paper, we will focus on applying the
method using the 4pLL (also called sigmoid Emax) model and the beta model,
respectively, as parameterized in Bornkamp et al. [2009]. Nevertheless, the
method is applicable to all parametric models with a first derivative that can
be calculated analytically. Exemplary, the first derivative of the 4pLL model
is given by

$f^{\prime}_{4pLL}(t,(a,b,c,h))=\frac{b\cdot c^{h}\cdot t^{h-1}\cdot
h}{(c^{h}+t^{h})^{2}}$

for $\theta=(a,b,c,h)\in\mathbb{R}^{4}$. The first derivative of the beta
model is

$\displaystyle f^{\prime}_{beta}(t,(a,b,\delta_{1},\delta_{2}))=$
$\displaystyle b\cdot
B_{\delta}\left(\frac{t^{\delta_{1}-1}\left(1-\frac{t}{scal}\right)^{\delta_{2}-1}\left(\delta_{1}scal-(\delta_{1}+\delta_{2})t\right)}{scal^{\delta_{1}+1}}\right)$

for $\theta=(a,b,\delta_{1},\delta_{2})\in\mathbb{R}^{4}$, where $scal$
denotes a fixed scaling parameter and
$B_{\delta}=\frac{(\delta_{1}+\delta_{2})^{\delta_{1}+\delta_{2}}}{\delta_{1}^{\delta_{1}}\delta_{2}^{\delta_{2}}}$.

###  2.2 Hypothesis test for the identification of significant changes

Here, we introduce a hypothesis test in which the rejection of $H_{0}$ leads
to the assertion of a significant change in gene expression. In the previously
introduced parametric framework this means that $f^{\prime}$ exceeds a
threshold $\lambda\geq 0$ at some time point t, which leads to the hypotheses

$\displaystyle H_{0}:\ \forall t\in T\ |f^{\prime}(t,\theta)|\leq\lambda\text{
vs. }H_{1}:\ \exists t\in T\ |f^{\prime}(t,\theta)|>\lambda.$ (1)

Rejecting $H_{0}$ implies the existence of at least one time point at which a
change in gene expression is detected that exceeds the critical value
$\lambda$. It is important to note that $\lambda$ should be selected
beforehand. Should any change in gene expression be of interest, $\lambda$ can
be set to $0$. However, if only a larger change may be biologically relevant,
$\lambda$ should be adjusted accordingly in consultation with experts. For
instance, regarding the absolute difference of the model fit to
$\log_{2}$-transformed counts at two time points $t_{a}$ and $t_{b}$, i.e.
$|f(t_{a},\theta)-f(t_{b},\theta)|$, the corresponding threshold is often set
to $\log_{2}(1.5)$ (see, e.g., Ghallab et al. [2021], Kappenberg et al. [2021]
and Möllenhoff et al. [2022]), investigating whether the foldchange between
$t_{a}$ and $t_{b}$ exceeds $1.5$. Based on this, a possible choice would be
given by $\lambda=\frac{\log_{2}(1.5)}{\overline{t}}$ for the test regarding
$f^{\prime}$, representing the slope of a straight line that increases by
$\log_{1}(1.5)$ between $t_{1}$ to $t_{m}$. Of note, $\overline{t}$ is
interchangeable with any other time frame, e.g. $\frac{1}{2}\overline{t}$.
Then the question is whether there is an absolute change in $\log_{2}(1.5)$
over half of the study period.

Next, let $L^{\alpha}(t,\hat{\theta})$ denote a $(1-\alpha)$ lower
simultaneous confidence band of $|f^{\prime}(t,\theta)|$, that is

$\mathbb{P}\\{\forall t\in
T:L^{\alpha}(t,\hat{\theta})\leq|f^{\prime}(t,\theta)|\\}\geq 1-\alpha.$

According to Möllenhoff et al. [2022], rejecting $H_{0}$ in (1) if
$L^{\alpha}(t,\hat{\theta})>\lambda$ for any $t\in T$, leads to an
$\alpha$-level test, as

$\displaystyle\mathbb{P}_{H_{0}}\\{\exists t\in
T:L^{\alpha}(t,\hat{\theta})>\lambda\\}$
$\displaystyle\leq\mathbb{P}_{H_{0}}\\{\exists t\in
T:L^{\alpha}(t,\hat{\theta})>|f^{\prime}(t,\theta)|\\}$
$\displaystyle=1-\mathbb{P}_{H_{0}}\\{\forall t\in
T:L^{\alpha}(t,\hat{\theta})\leq|f^{\prime}(t,\theta)|\\}$ $\displaystyle\leq
1-(1-\alpha)=\alpha.$

The time period of significant change in gene expression corresponds to the
set of all $t\in T$ for which $H_{0}$ in (1) can be rejected, i.e.
$L^{\alpha}(t,\hat{\theta})>\lambda$, which we denote by
$\tilde{S}\coloneqq\\{t|L^{\alpha}(t,\hat{\theta})>\lambda\\}.$

Figure 1: (A) Visualisation of the time period of significant change in gene
expression $\tilde{S}$ (orange) with two time points of interest. The red
dotted lines indicate the beginning and end of $\tilde{S}$. The yellow line
indicates the time point of maximum change in expression. $\tilde{S}$ is
partitioned into coherent subsets if the gene expression curve follows a 4pLL
model (B) or a beta model (C). B1 and C1 depict the fitted model, B2 and C2
the corresponding first derivative and confidence bands such that $H_{0}$ in
(1) is rejected.

Figure 1 A depicts an illustrative example of a first derivative of a gene
expression curve $f^{\prime}$ (black) and its lower confidence band (blue).
The orange line represents the time period $\tilde{S}$, during which the
confidence band exceeds the threshold $\lambda$, indicated by the red line.

While we can conclude significant changes in gene expression for all $t\in T$
which fulfill $L^{\alpha}(t,\hat{\theta})>\lambda$, some specific time points
might be of higher interest. An example are the beginning and end of a
significant change. Precisely, these are the first time point $t_{start}\in T$
with $L^{\alpha}(t_{start},\hat{\theta})>\lambda$ and the last time point
$t_{end}\in T$ with $L^{\alpha}(t_{end},\hat{\theta})>\lambda$. Figure 1 A
indicates both $t_{start}$ and $t_{end}$ by the red vertical lines and thus
also the period of significant change $\tilde{S}$. Another example is the time
point of maximal change, indicated by the yellow vertical line in Figure 1 A,
representing the time point $t_{max}\in T$ at which
$L^{\alpha}(t,\hat{\theta})-\lambda$, $t\in\tilde{S}$ assumes its (locally)
maximal value.

These time points of interest are not necessarily unique. In case of a non-
monotonous function $f$ (e.g. beta model), the confidence band can dip below
the threshold $\lambda$ and than exceed it again after a time period of no
significant change. Therefore, we can partition $S$ into ascending coherent
subsets $\tilde{S}=\big{(}S_{1},\dots,S_{s}\big{)},$ meaning that for all
$j=1,\dots,s-1$, it exists a $t^{\prime}\in T$ between $S_{j}$ and $S_{j+1}$
such that $L^{\alpha}(t^{\prime},\hat{\theta})<\lambda$ and
$\max(S_{j})<\min(S_{j+1})$. We refer to Figure 1 B and C for a visual
representation. Figure 1 B depicts a gene expression curve that follows a 4pLL
model on the left and its first derivative with a possible confidence band on
the right. Here, $H_{0}$ is rejected and the time period of significant change
in gene expression is coherent. Then $\tilde{S}$ equals $S_{1}$. Figure 1 C
depicts a beta model in the left and its first derivative with two possible
confidence bands on the right. The solid light-blue line represents the
scenario in which the decline in gene count in the latter part of the study
period is not statistically significant. The dotted blue line on the other
hand represents the scenario in which the aforementioned decline is
significant. Here, $\tilde{S}$ splits into the coherent sub-periods $S_{1}$
and $S_{2}$. Each $S_{j}$, $j=1,\dots,s$, has a unique set of the
aforementioned time points of interest. In the following, we will focus on
$t_{start}$ and $t_{end}$ as an example. Therefore, the set of times to search
for is given by

$\displaystyle\\{t^{(j)}_{start}=\min(S_{j})\mid
j=1,\dots,s\\}\cup\\{t^{(j)}_{end}=\max(S_{j})\mid j=1,\dots,s\\}.$ (2)

The precise estimation of these quantities depends on the accuracy of the
estimation of $L^{\alpha}(t,\hat{\theta})$, which is described in the next
section.

###  2.3 Estimating a confidence band of $\mathbf{|f^{\prime}(t)|}$

In order to estimate a lower simultaneous confidence band
$L^{\alpha}(t,\hat{\theta})$ of $|f^{\prime}(t,\theta)|$ the method developed
in Möllenhoff et al. [2022] is adjusted to the first derivative. Therefore, we
define

$L^{\alpha}(t,\hat{\theta}):=|f^{\prime}(t,\hat{\theta})|-c\hat{\sigma}_{|f^{\prime}(t,\hat{\theta})|}$

and estimate the critical value $c$, such that

$\mathbb{P}\bigg{\\{}\max_{t\in
T}\frac{|f^{\prime}(t,\hat{\theta})|-|f^{\prime}(t,\theta)|}{\hat{\sigma}_{|f^{\prime}(t,\hat{\theta})|}}\leq
c\bigg{\\}}=1-\alpha.$

To estimate both $\hat{\sigma}_{|f^{\prime}(t,\hat{\theta})|}$ and $c$, a two-
step parametric bootstrap procedure inspired by the classic bootstrap-t method
(see Efron and Tibshirani [1994]) has been proposed in Algorithm 1 by
Möllenhoff et al. [2022] and is applied here by substituting the difference
for the first derivative. In short, $\hat{\theta}$ is estimated from the
original data set and used to generate $B_{1}$ bootstrap samples. This is done
by randomly drawing iid. errors $\varepsilon_{p,q}^{*}\sim N(0,1)$ and
applying the formula

$y_{p,q}^{*}=f(t_{p},\hat{\theta})+\hat{\sigma}\varepsilon_{p,q}^{*}\ \text{
for }\ p=1,\dots,m,\ q=1,\dots,n_{p}.$

These can be considered the first level of the bootstrap. For each of the
generated samples, $\hat{\theta}^{*}_{l}$, $l=1,\dots,B_{1}$, are estimated
and used to generate another $B_{2}$ bootstrap samples, which constitute the
second level of the bootstrap. Then,
$\hat{\sigma}_{|f^{\prime}(t,\hat{\theta})|}$ is estimated from the first
bootstrap level by calculating the empirical variance of the sample
$|f^{\prime}(t,\hat{\theta}^{*}_{1})|,\dots,|f^{\prime}(t,\hat{\theta}^{*}_{B_{1}})|$.
Similarly, $\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}^{*}_{l})|}$,
$l=1,\dots,B_{1}$, is calculated for each of the generated first-level
bootstrap samples using the second-level bootstrap samples. With these values

$D^{*,l}:=\max_{t\in
T}\frac{|f^{\prime}(t,\hat{\theta}_{l}^{*})|-|f^{\prime}(t,\hat{\theta})|}{\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}_{l}^{*})|}}$

can be calculated for each $l=1,\dots,B_{1}$. The respective empirical
$(1-\alpha)$-quantile of the distribution yields an estimate of the critical
value $c$. A visual representation of this process is provided in the blue
rectangle in Figure A (see supplementary material).

###  2.4 Confidence intervals for the time points of interest

In the previous sections we obtained time points of significant change in gene
expression. However, as estimators, they are subject to uncertainty. Thus, the
construction of confidence intervals for these time points is of huge
practical importance. The method proposed is an extension to the two-step
bootstrap described in Section 2.3. Precisely, another bootstrap level is
added in advance, and the two-step procedure described above is then applied
to each sample. This provides an estimated lower confidence band for each of
the generated data sets, allowing calculation of the time frames of
significant change in gene expression and the time points of interest. In the
following analysis, only those cases with the same number of coherent subsets
as the original data set are considered. The corresponding time points of
interest from all bootstrap samples are combined (e.g. all $t^{*}_{start}$ and
all $t^{*}_{end}$, where the ’*’ denotes the bootstrap data) and used to
estimate the empirical $\frac{\alpha}{2}$\- and
$(1-\frac{\alpha}{2})$-quantiles , say $q_{\frac{\alpha}{2}}$ and
$q_{(1-\frac{\alpha}{2})}$ providing percentile confidence intervals for the
time points of interest of the original data set that can be expressed as
$[q_{\frac{\alpha}{2}},q_{(1-\frac{\alpha}{2})}]$. Details can be found in the
Supplemental material and the Appendix: Figure A (see supplementary material)
depicts a visualization of this concept, while Algorithm 1 summarizes the
procedure.

##  3 Simulation study

###  3.1 Setting and data generation

In order to evaluate the performance of our newly proposed method, a
simulation study was conducted. Motivated by Möllenhoff et al. [2022] and
Kappenberg et al. [2021] we focused on the 4pLL and the beta model. The
simulation settings were inspired by the Western diet mice trial discussed in
Section 4. Precisely, we used the parameters obtained from the expression of
the genes ’Cd163’ (4pLL model, see Section 4.1 for details) and ’Fam83a’ (beta
model) in WD-fed mice. This resulted in the six basic scenarios, which are
summarized in Table 1 and depicted in Figure 2.

Scenarios 1-3 were generated using the 4pLL model, while the latter three were
generated with the beta model. In Scenario 1, there is no relevant change in
gene expression respective to the threshold
$\lambda=\frac{\log_{2}(1.5)}{45}\approx 0.0130$. Therefore, $H_{0}$ in (1)
cannot be rejected. In all other scenarios, there is a change in expression
larger than $\lambda$, and, therefore, these scenarios correspond to the
situation under the alternative. Scenario 4 reflects the unaltered parameters
as estimated. Scenarios 2, 3, 5 and 6 were based on the real examples, but
changed to reflect a different slope, see Table 1.

Table 1: Overview of simulation scenarios. Sce. | Short name | Model | Parameters | Description
---|---|---|---|---
1 | No relevant change | 4pLL |  $a=8.791$, $b=-0.089$ | Does not exceed $\lambda=0.0130$
|  |  |  $c=17.589$, $h=10.000$ | Reflects the margin of $H_{0}$.
2 | Small Jump | 4pLL |  $a=8.791$, $b=-0.946$ | Exceeds $\lambda=0.0130$
|  |  |  $c=17.589$, $h=10.000$ | $\tilde{S}=(11.7,24.5)$
3 | Large Jump | 4pLL |  $a=8.791$, $b=-3.783$ | Exceeds $\lambda=0.0130$
|  |  |  $c=17.589$, $h=5.000$ | $\tilde{S}=(5.9,36.3)$
4 | Dip | beta |  $a=6.997$, $b=2.952$ | Exceeds $\lambda=0.0130$ twice
|  |  |  $\delta_{1}=0.506$, $\delta_{2}=0.215$ |  $S_{1}=(0,33.6)$, $S_{2}=(41.2,45)$
5 | Dip alternative | beta |  $a=6.997$, $b=2.952$ | Exceeds $\lambda=0.0130$ twice
|  |  |  $\delta_{1}=3.286$, $\delta_{2}=1.290$ |  $S_{1}=(5.7,38.1)$, $S_{2}=(39.4,45)$
6 | No Dip | beta |  $a=6.997$, $b=2.952$ | Exceeds $\lambda=0.0130$
|  |  |  $\delta_{1}=0.228$, $\delta_{2}=0.084$ | $\tilde{S}=(0,28.9)$
Figure 2: Visualization of the different scenarios. The horizontal red line
indicates the relevance threshold $\lambda=\frac{\log_{2}(1.5)}{45}$. Each
case depicts the true model on the left and the true first derivative on the
right. The vertical red lines mark the time frame of relevant change in
regards to $\lambda$.

The observations were simulated at seven time points ($0$, $3$, $9$, $15$,
$21$, $27$, $33$, $39$ and $45$ weeks) in accordance with the feeding times
observed in the Western diet mice trial, following the beginning of the study.
A total of $47$ data points were simulated consisting of $5$ mice at the first
$5$ time points and $4$ and $8$ mice at $39$ and $45$ weeks, respectively.

In order to gain a comprehensive understanding of the impact of the standard
deviation, five distinct levels for each of the six basic scenarios were
considered: ’small’,’mid-small’, ’medium’, ’mid-large’ and ’large’. The
’medium’ standard deviation of the realistic Scenario 4 is as estimated from
the real world data. For all other scenarios, the ’medium’ standard deviation
was calculated via linear transformation with regard to the respective maximum
absolute differences in gene expression compared to the realistic example of
the same model type. Regarding the parameterization introduced in Section 2.1,
the maximum absolute value corresponds to the parameter $b\in\theta$.
Accordingly, we designate $b_{sce}\in\theta_{sce}$ as the maximum absolute
difference of the scenario and $b_{real}\in\theta_{real}$ as the maximum
absolute difference estimated from the example. With this, we obtained the
transformed standard deviation by
$\sigma_{sce,medium}=\frac{b_{sce}}{b_{real}}\sigma_{real}$, where
$\sigma_{real}$ the estimated value from the real world data. From these
’medium’ values of standard deviation, the ’small’, ’mid-small’, ’mid-large’
and ’large’ levels were calculated by multiplying with $0.5$, $0.75$, $1.5$
and $2$, respectively, see Table 2 for details.

Table 2: Levels of standard deviation used for the simulation study. The standard deviations marked with an asterisk (∗) are approximately half the size of the maximum change in gene expression.  Scenario |  $\sigma$ small |  $\sigma$ mid-small |  $\sigma$ medium |  $\sigma$ mid-large |  $\sigma$ large
---|---|---|---|---|---
1 | $0.014$ | $0.021$ | $0.028$ | $0.042$ | $0.056^{*}$
2 | $0.149$ | $0.223$ | $0.297$ | $0.446$ | $0.595^{*}$
3 | $0.595$ | $0.892$ | $1.190$ | $1.784^{*}$ | $2.380^{*}$
4 | $0.283$ | $0.424$ | $0.566$ | $0.849$ | $1.132$
5 | $0.283$ | $0.424$ | $0.566$ | $0.849$ | $1.132$
6 | $0.283$ | $0.424$ | $0.566$ | $0.849$ | $1.132$

Given either the model function, $f$, of the beta model or the 4pLL model and
the values of $\theta$ and $\sigma$, the gene expression data was simulated by
drawing iid. errors $\varepsilon^{sim}_{p,q}\sim N(0,1)$ and applying the
formula

$y_{p,q}^{sim}=f(t_{p},\hat{\theta})+\hat{\sigma}\varepsilon^{sim}_{p,q}\text{
for }p=1,\dots,m,\ q=1,\dots,n_{p}.$

For the two-step bootstrap procedure we chose $B_{1}=500$ and $B_{2}=25$
repetitions. Each procedure was repeated in $1,000$ simulation runs. Outcomes
of the simulation study were the number of rejections of $H_{0}$ in (1) and
the identified time frame of significant change in gene expression. The
majority of the code was run on R version 4.3.2 using a Windows 11 x64 PC. The
simulation scenarios were run on a HPC-cluster employing R version 4.3.1. The
runtime of each simulated study on a commercially available Windows PC is
approximately 5 minutes, with a variability of a few minutes depending on the
number of CPU cores.

###  3.2 Results

####  3.2.1 Rejection of $\mathbf{H_{0}}$

Table 3 shows the number of rejections of $H_{0}$ in $1,000$ runs for all
configurations. The left column in Figure 3 depicts the estimated confidence
bands for a mid-small and mid-large standard deviation. We refer to the
supplementary material Figures D-F for the other cases.

Table 3: Number of rejections of $H_{0}$ in $1,000$ runs. $\sigma$ level | Sce. 1 | Sce. 2 | Sce. 3 | Sce. 4 | Sce. 5 | Sce. 6
---|---|---|---|---|---|---
small $\sigma$ | 0 | 999 | 1000 | 1000 | 1000 | 1000
mid-small $\sigma$ | 0 | 946 | 818 | 1000 | 1000 | 1000
medium $\sigma$ | 0 | 830 | 514 | 1000 | 1000 | 1000
mid-large $\sigma$ | 0 | 627 | 401 | 989 | 999 | 978
large $\sigma$ | 0 | 372 | 303 | 814 | 957 | 780

In Scenario 1, no false rejection of $H_{0}$ was observed in any of the cases,
indicating a type I error of 0. This scenario corresponds to the margin of the
null hypothesis in (1), here given by $\lambda=0.0130$. The conservatism of
the test decision even for a large standard deviation may therefore be
partially attributed to the effect being very small and an overall tendency to
slightly underestimate the time frame of significant change in expression
discussed in the next subsection.

In general, our proposed method shows high power when the beta model is the
true model. Up to a medium standard deviation, $H_{0}$ is always rejected for
Scenarios 4, 5 and 6. Even in the case of a large standard deviation, the
power remains at about $80\%$ for Scenarios 4 and 6 and above $95\%$ for
Scenario 5. Accordingly, the estimated confidence bands show a low
variability.

For the remaining scenarios generated with the 4pLL model – Scenario 2 and 3 –
we observe a large influence of the standard deviation. In particular,
significant changes in gene expression were observed for a mid-large or large
variability below $650$ and up to $303$. This influence seems to be more
pronounced for the larger jump in gene expression simulated in Scenario 3. In
this case, only $514$ of significant changes were detected for a medium
standard deviation.

It is noteworthy that the (mid-)large standard deviations (see Table 2)
employed in the generation of Scenarios 2 and 3 is considerably large compared
to the absolute change in gene expression. Naturally, this leads to a high
variability in the generated simulation studies (see Figure B in the
supplementary material). However, we observe that in the case of the 4pLL
model, the translation of a high degree of variability on the count scale
leads to an extreme degree of variability on the scale of the change in gene
expression considered when applying the first derivative (see C in the
supplementary material). It is important to interpret the results presented in
light of these considerations. It is also important to note that the same
effect occurs when generating the first and second level bootstrap runs. This
is likely to have an impact on the estimation process, especially if the data
set has a high degree of variability.

Figure 3: Results of the six simulation scenarios for a mid-small (A) and mid-
large (B) standard deviation. The left column depicts the true first
derivative and the estimated confidence bands. The right column depicts the
estimated time frames of significant change in gene expression for each run.

For Scenarios 2 and 3 we note a visible crater in the estimated confidence
bands around the week where the true first derivative reaches its maximum,
which corresponds to the inflection point of the model function. While the
extent is likely to be at least partially attributable to the size of the
standard deviation, the trend is discernible across all cases. Given that the
effect appears to be more pronounced in Scenario 3 than in Scenario 2, it is
plausible that its magnitude may correlate with the steepness of the true
model function in the area around its inflection point. This conclusion is
supported by the absence of a comparable crater for Scenario 5, the only
scenario with an inflection point generated using the beta model, which is
markedly less steep than Scenarios 2 and 3.

####  3.2.2 Estimated time frame of significant change in gene expression

While the number of rejections of $H_{0}$ is largely satisfactory, a reliable
detection of the time frame of significant change in gene expression is of
greater interest. Therefore, we considered the bias and variance of each
start/end time point of the estimated time frame of significant change in gene
expression. Furthermore, we examined the number of runs in which the true
number of coherent subsets of $\tilde{S}$ was correctly identified. We refer
to Table A in the supplementary material for a summary of all results. As
anticipated, the true number of coherent subsets was consistently detected for
Scenarios 2 and 3 if $H_{0}$ was rejected. Given the monotonicity of the 4pLL
model, the time frame of significant change in expression is, by definition,
continuous. However, as described above, a crater was observed in some
estimated lower simultaneous confidence bands. In a small number of cases,
this resulted in the identification of two disjointed subsets. Overall, the
number of these cases is less than $12$ for Scenario 2. In case of Scenario 3,
the same is observed except for a mid-small and medium standard deviation,
which respectively fall within the ranges of $200$ and $70$ cases.

While $H_{0}$ was reliably rejected for Scenarios 4 and 5, the second time
frame of significant change in expression was not always detected. The smaller
dip in Scenario 4 was only identified approximately $120$ times when the
standard deviation was small and approximately $40$ times otherwise. The more
pronounced dip in Scenario 5 was detected $999$ times for a small standard
deviation, $671$ times for a medium standard deviation and $260$ times for a
large standard deviation. However, it should be noted that in both scenarios,
the maximum value in gene count is reached around week $39$, which is the
second to last sample time point. Therefore, the detection of the dip is
contingent only on the last measurement. In Scenario 6, a dip at the end was
incorrectly identified less than $15$ times in all cases.

The right columns in Figures 3 depict the estimated time frames of significant
change in gene expression for each of the $1,000$ simulation runs for the mid-
small and mid-large standard deviations, respectively. We refer to Figures and
D-F in the supplementary material for the other cases. The true values are
indicated by the red lines. In general, we note that the estimated start of
the time frame of significant change in expression displays a positive bias
for all scenarios, whereas the estimated end displays a negative bias.
Therefore, the time frame is overall slightly underestimated. The bias depends
largely on the standard deviation and the slope of the first derivative.
Scenario 6 exhibits the smallest bias for the onset of the significant change
at the beginning of the study, never exceeding $0.375$ weeks. Scenario 4 is
similarly flat. Its bias never exceeds $2.379$. Scenario 5 is characterised by
a steeper slope and exhibits a noticeably larger bias. The range is from
$3.956$ weeks for a small standard deviation to $8.865$ weeks for a large
standard deviation. For Scenarios 2 and 3 – disregarding the crate when
detected – we observe similar results. The influence of the slope seems to
have a reverse effect on the bias of the end of the identified time frame of
significant change. In particular, flat scenarios, such as Scenario 6,
demonstrate a noticeably larger bias for all levels of standard deviation. In
addition, the bias is overall larger than for the start of the significant
change, with a range of $-1.225$ to $-23.591$ weeks.

In Scenarios 4 and 5, a second time frame of significant change in gene
expression was simulated. Since in both cases it is cut off at the end of the
study, we consider only the starts of the second coherent subsets of
$\tilde{S}$, respectively. Again, we note a positive bias, which ranges from
almost $2$ weeks to nearly $3.5$ weeks. Finally, the variance shows a similar
behavior to the bias.

##  4 Case study

In this section, three possible applications of our method are presented. For
illustration, we use data from the Western diet mice study (see Ghallab et al.
[2021]). A total of $79$ male mice were fed a Western (WD) or standard (SD)
diet for a maximum of $48$ weeks to assess the effect of a high fat diet on
the development of non-alcoholic fatty liver disease. Up to $8$ mice per diet
were sacrificed at $3,\ 6,\ 30,\ 36,\ 42$ and $48$ weeks, and WD-fed mice were
additionally sacrificed at $12,\ 18$ and $24$ weeks. RNA was extracted from
frozen liver tissue and RNA-seq analysis was performed. Here, we use the pre-
processed data set, that was generated using Salmon (see Patro et al. [2017])
and the R packages tximeta (see Love et al. [2020]) and DESeq2 (see Love et
al. [2014]). The count data was normalized and $\log_{2}$-transformed using
the vst() function (see Tibshirani [1988], Huber et al. [2003] and Anders and
Huber [2010]) from the DESeq2 package. As the initial samples were collected
at $3$ weeks, the beginning of the study at $t=0$ in the subsequent analysis
represents these initial data points.

The first example is an application to a single gene. In the second example,
our method is used to assess the effect of SD vs. WD on gene expression, using
a gene that shows differential expression in both cases as an example. In the
third example, we select a pool of nearly $10000$ genes and use our method to
identify those that change significantly between 10 and 25 weeks. A Gene
Ontology (GO) enrichment analysis is then performed on the identified genes.
In all cases, the two-step bootstrap procedure was implemented with
$B_{1}=500$ and $B_{2}=25$ bootstrap runs, while the three-step bootstrap
procedure was executed with an additional $B_{3}=500$.

###  4.1 Application to an individual gene

’Cd163’ is a gene found in both mice (see National Center for Biotechnology
Information [b]) and humans (see National Center for Biotechnology Information
[a]). The soluble form of the encoded molecule was found to be a predictor for
both fibrosis and hepatocellular carcinoma development in nonalcoholic
steatohepatitis in humans, as reported in Kawanaka et al. [2023]. In
particular, decreased levels were observed in patients with NAFLD who
exhibited alleviated fibrosis and inflammation. In Ghallab et al. [2021], the
gene has been identified as a RJG. This implies that the gene in question is
only deregulated after a specific period of time, in this case, after $24$
weeks. In a preceding model selection step, both the 4pLL and the beta model
were fitted to the gene expression data of the WD-fed mice. The AIC for the
4pLL model was $90.4$, while that for the beta model was $98.2$. Therefore,
the 4pLL model was used for further analyses.

Figure 4: On the left, the fit of the 4pLL model to the normalized gene count
of the gene ”Cd163” over time is depicted. On the right, the first derivative
and the estimated confidence band of the model fit are displayed. The
horizontal red line indicates $\lambda=\frac{\log_{2}(1.5)}{45}$. The vertical
dashed red lines indicate the corresponding time frame of significant change
in gene expression. The yellow colored lines correspond to an alternative
choice of $\lambda=\frac{\log_{2}(1.5)}{27.5}$. The black lines represent the
case $\lambda=0$.

The left panel of Figure 4 depicts the estimated gene expression curve. In
accordance with the identification as a RJG, it appears that approximately
$10-12$ weeks pass, before the gene becomes down-regulated. Following an
additional $10-12$ weeks, the deregulated expression appears to reach a
constant value. Overall the gene changes by a value of $-1.891$. In the right
panel of Figure 4, the corresponding first derivative and the estimated
confidence band are depicted. The various horizontal lines represent different
potential choices for the threshold $\lambda$, while the corresponding
vertical lines represent the respective time frames of significant change in
gene expression. The black line represents the case where $\lambda=0$, while
the red line corresponds to the aforementioned case
$\lambda=\frac{\log_{2}(1.5)}{45}=0.0130$ (see Section 2.2). Similarly, a
third threshold was calculated, indicated by the yellow line, which
corresponds to a change of $\log_{2}(1.5)$ in gene expression over the span of
half the study. This third threshold is given by
$\lambda=\frac{\log_{2}(1.5)}{27.5}=0.0213$.

We observe that $H_{0}$ is rejected for all thresholds. In the following, we
will present the analysis for $\lambda=\frac{\log_{2}(1.5)}{45}$ as an
illustrative example. The period during which the estimated significant
changes in gene expression occur begins at $15.1$ weeks ($95\%$ CI
$\left[13.2,17.8\right]$) and ends at $20.3$ weeks ($95\%$ CI
$\left[17.4,23.5\right]$). Consequently, we can conclude the classical RJG
form of a delayed deregulation. The change is observed over a period of $5.2$
weeks, with a plateau following the initial ’jump’ in activity. Over the
course of the $5.2$ week period, the expression undergoes a change of $1.189$.
Therefore, we can account for $62.883\%$ of the observed maximal change via a
significant effect.

For comparison, the other cases indicate a time period of significant change
in gene expression between weeks $14.8$ to $20.7$ ($\lambda=0$, black) and
$15.3$ to $20.1$ ($\lambda=0.0213$, yellow), respectively. Consequently, while
the ’Cd163’ gene was also identified as a RJG, the estimated time frame of
deregulation is set earlier than reported in Ghallab et al. [2021] for all
thresholds.

###  4.2 Aging

For this application we consider a gene that deregulates in both diet groups
and compare the estimated time frames of significant change in gene
expression. The biological question inspiring this example is whether the WD
accelerates or decelerates, and thereby ’ages’ biological processes, which
also occur in SD-fed mice. The chosen example, ’Dbp’, is a protein encoding
gene detected in humans National Center for Biotechnology Information [c] and
mice National Center for Biotechnology Information [d] associated with the
circadian rhythm. A preliminary model selection step indicated a superior fit
of the 4pLL model to the gene expression data for both diet groups.
Consequently, it was utilised in the subsequent analysis.

Figure 5: A comparison of the SD versus the WD, exemplary for the gene ”Dbp”.
The gene counts, as modelled by the 4pLL model, are depicted on the left. On
the right, the first derivatives and their estimated confidence bands are
shown. The red dotted line outlines the time period of significant change in
gene expression for the WD, while the dot-dashed line outlines the analogous
period for the SD.

The left panel of Figure 5 depicts the modeled gene expression of both the SD
and the WD groups. The gene count in SD-fed mice initially appears to be
smaller and both genes seem to undergo similar deregulation at approximately
the same time and converge to the same value. On the right, the corresponding
first derivatives and the estimated confidence bands are displayed. The dotted
line represents the time period of significant change in gene expression for
the WD-fed group, while the dot-dashed line outlines the time period for the
SD-fed mice. Although the first derivative of the WD-fed group appears to
increase earlier, the confidence band of the SD group exceeds the threshold
$\lambda=\frac{\log_{2}(1.5)}{45}=0.0130$ earlier. The estimated time frame of
significant change in gene expression in WD-fed mice ranges from $36.9$ weeks
($95\%$ CI $\left[33.4,38.4\right]$) to the end of study at $45$ weeks ($95\%$
CI $\left[40.5,45.0\right]$) , indicated by the red dotted lines. Over the
span of $8.1$ weeks the gene count undergoes an absolute change of $2.213$.
Similarly, the estimated time period during which significant change in gene
expression in SD-fed mice occurs begins at $35.5$ weeks ($95\%$ CI
$\left[32.2,36.3\right]$) and ends with the end of study at $45$ weeks ($95\%$
CI $\left[43.2,45.0\right]$). The gene count undergoes a change of $2.996$
over the span of $9.5$ weeks. Overall, both genes show very similar results.
The estimated confidence band for the SD-fed mice exceeds the threshold
$\lambda=0.0130$ about $-1.4$ weeks earlier and the estimated $95\%$
confidence interval for this difference is given by $\left[-4.9,2.0\right]$.

###  4.3 GO enrichment analysis

For this example we applied our proposed method to identify a subset of
differentially expressed genes that deregulate during a pre-selected time
period ($10$ to $25$ weeks). In a subsequent step, we conducted a GO
enrichment analysis on these genes. This means, we examined if a significant
number of the identified genes were assigned to pre-specified GO groups (see
Consortium [2001a]), in which the genes were grouped according to their
connection to biological functions. This approach would permit a practitioner
to draw inferences regarding the time frame within which specific biological
processes occur.

To achieve this, we used the differential gene expression analysis as
performed in Ghallab et al. [2021] using the DESeq2 package to pre-select the
subset of genes that are estimated to undergo significant changes in
expression between week $0$ (time on study) in SD-fed mice (SD0) and week $45$
in WD-fed mice (WD45). Precisely, the pre-processed data set was created by
applying a generalized linear model with one factor, a combination of the diet
and the week. The SD0 group was used as reference. To account for the
potentially high variability of the counts, we used this as basis and
performed a shrinkage of the effect size WD45 vs. SD0 as described in Zhu et
al. [2018] by applying the lfcShrink() function from the DESeq2 R package.
Afterwards we identified all genes with an adjusted p-value (see Benjamini and
Hochberg [1995]) smaller than $0.05$ when comparing the expression of SD-fed
mice at $0$ weeks to the expression of WD-fed mice at $48$ weeks.

Following the exclusion of $22$ genes due to measurement issues, the total
number of identified genes was $9881$. For all the aforementioned genes, we
assumed that the 4pLL model was the most appropriate for the analysis and
applied our method to identify the time periods of significant change in gene
expression. Figure 6 depicts the estimated confidence bands.

Figure 6: Estimated confidence bands of the $9881$ preselected genes. The
darker blue lines correspond to the genes for which the estimated period of
significant change in expression was determined to be between $10$ and $25$
weeks.

The results were then employed to identify the genes with an estimated time
period of significant change in expression between weeks $10$ and $25$. A
total of $569$ genes were identified.

A GO enrichment analysis was conducted on the aforementioned genes using the
topGo R package (see Alexa and Rahnenführer [2024] and Alexa et al. [2006]).
In essence, Fisher’s exact test is employed to verify whether a number of
interesting genes is annotated to a specific functional group that exceeds
what would be expected by chance. The $20$ gene ontology groups with the
smallest adjusted p-value are presented in Table B in the supplementary
material. For instance, the majority of the genes ($525$) were assigned to the
GO group ’inflammatory response’. This is in accordance with the findings that
inflammation is a pivotal event in the progression of NAFLD (see Ghallab et
al. [2021] and Sutti and Albano [2020]). In Ghallab et al. [2021] few
inflammatory foci were observed in week $3$ and the number strongly increased
after week $36$. Lipogranulomas, single or multiple fat globules surrounded by
chronic inflammatory cells and Kupffer cells (see Brunt [2002]), were observed
as early as week $6$ and increased at week $18$.

##  5 Discussion and conclusion

In this paper, we proposed a parametric method to estimate the time frame
during which the expression of a gene changes significantly. Thereby, we have
provided a tool, that allows for an interpretation that is not limited to the
onset of an effect or the observed time points and produces results based on a
hypothesis test to the $\alpha$-level. As demonstrated by the 4pLL model and
the beta model, the approach is applicable to a wide range of potential
models, with the singular restriction that the first derivative can be derived
analytically from the model function. It can be employed to analyse the
expression of a single gene or, when used in conjunction with a GO enrichment
analysis, to infer biological function. Moreover, the methodology can be
applied to other forms of time-response or dose-response data. By extending
the method to a three-step bootstrap approach, we additionally provided the
means to estimate confidence intervals for specific time points of interest.
This provides further evidence for the results of the statistical analysis and
enables a comparison of multiple gene curves.

Of course, there are limits to the method. The simulation results discussed in
Section 3.2 show a clear tendency to underestimate the time frame of
significant change in expression. The effect is relatively modest for medium
levels of data variability. However, it becomes more pronounced with greater
variability. Overall, the method would benefit from a built-in stabilization
procedure against the predictable effects of a large standard deviation.
Moreover, a potential solution for a real data set with considerable
variability could be the replacement of the model selection phase with model
averaging (see Claeskens and Hjort [2008] and Fletcher [2019]). It has been
demonstrated that model averaging offers advantages over model selection, as
it circumvents the potential bias introduced by model misspecifications (see
Bornkamp [2015], Breiman [1996] and, recently, Hagemann and Möllenhoff
[2024]).

Often, a statistical conclusion on more than one gene expression curve is of
interest. In Section 4.3, we presented a possible application to a large
number of genes in combination with a GO enrichment analysis. Nevertheless, a
multitude of potential biological questions may necessitate alternative
approaches. This frequently gives rise to to a multiple testing problem.
Hence, our proposed method would highly benefit from the development of a
method, that corrects for this issue beyond a simple correction of the
p-values.

In conclusion, we posit that our proposed method is a valuable addition to the
existing tools for the analysis of gene expression data. It allows for an
comprehensive analysis over the entire study period and is applicable to a
wide range of expression curves.

## Supplementary material

Supplementary material providing additional tables and figures can be found
online. Corresponding R code, which can be used to reproduce the analysis of
the case study and the simulation results, is available at
https://github.com/LuciaAmeis/Identification-of-changes-in-gene-expression.

## Funding

This work has been supported by the Research Training Group ”Biostatistical
Methods for High-Dimensional Data in Toxicology” (RTG 2624, P7) funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation - Project
Number 427806116).

## References

  * Alexa and Rahnenführer [2024] A. Alexa and J. Rahnenführer.  topGO: Enrichment Analysis for Gene Ontology.  2024\.  doi: 10.18129/B9.bioc.topGO.  R package version 2.56.0.
  * Alexa et al. [2006] A. Alexa, J. Rahnenführer, and T. Lengauer.  Improved scoring of functional groups from gene expression data by decorrelating GO graph structure.  _Bioinformatics_ , 22(13):1600 – 1607, 2006.  doi: 10.1093/bioinformatics/btl140.
  * Anders and Huber [2010] S. Anders and W. Huber.  Differential expression analysis for sequence count data.  _Genome Biology_ , 2010.  doi: 10.1186/gb-2010-11-10-r106.
  * Benjamini and Hochberg [1995] Y. Benjamini and Y. Hochberg.  Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 57(1):289 – 300, 1995.  doi: 10.1111/j.2517-6161.1995.tb02031.x.
  * Bornkamp [2015] B. Bornkamp.  Viewpoint: model selection uncertainty, pre-specification, and model averaging.  _Pharmaceutical statistics_ , 14(2):79 – 81, 2015.  doi: 10.1002/pst.1671.
  * Bornkamp et al. [2009] B. Bornkamp, J. Pinheiro, and F. Bretz.  MCPMod: An R Package for the Design and Analysis of Dose-Finding Studies.  _Journal of Statistical Software_ , 29(7):1 – 23, 2009.  doi: 10.18637/jss.v029.i07.
  * Breiman [1996] L. Breiman.  Heuristics of instability and stabilization in model selection.  _The Annals of Statistics_ , 24(6):2350 – 2383, 1996.  doi: 10.1214/aos/1032181158.
  * Brunt [2002] E.M. Brunt.  Alcoholic and nonalcoholic steatohepatitis.  _Clinics in Liver Disease_ , 6(2):399 – 420, 2002.  doi: doi.org/10.1016/S1089-3261(02)00002-8.
  * Claeskens and Hjort [2008] G. Claeskens and N.L. Hjort.  _Model Selection and Model Averaging_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2008.  doi: 10.1017/CBO9780511790485.
  * Consortium [2001a] G. O. Consortium.  Creating the gene ontology resource: design and implementation.  _Genome Research_ , 11(8), 2001a.  doi: 10.1101/gr.180801.
  * Consortium [2001b] International Human Genome Sequencing Consortium.  Initial sequencing and analysis of the human genome.  _Nature_ , 409:860 – 921, 2001b.  doi: 10.1038/35057062.
  * Delignette-Muller et al. [2011] M.-L. Delignette-Muller, C. Forfait, E. Billoir, and S. Charles.  A new perspective on the Dunnett procedure: Filling the gap between NOEC/LOEC and ECx concepts.  _Environmental Toxicology and Chemistry_ , 30(12):2888–2891, 2011.  doi: 10.1002/etc.686.
  * Efron and Tibshirani [1994] B. Efron and R.J. Tibshirani.  _An Introduction to the Bootstrap_.  CRC Press, 1st ed. edition, 1994.  doi: 10.1201/9780429246593.
  * Fahrmeir et al. [2022] L. Fahrmeir, T. Kneib, S. Lang, and B. D. Marx.  _Regression: Models, Methods and Applications_.  Springer Berlin, Heidelberg, 2nd edition, 2022.  doi: 10.1007/978-3-662-63882-8.
  * Fletcher [2019] D. Fletcher.  _Model Averaging_.  SpringerBriefs in Statistics. Springer Berlin, 2019.  doi: 10.1007/978-3-662-58541-2.
  * Ghallab et al. [2021] A. Ghallab, M. Myllys, A. Friebel, J. Duda, K. Edlund, E. Halilbasic, M. Vucur, Z. Hobloss, L. Brackhagen, B. Begher-Tibbe, R. Hassan, M. Burke, E. Genc, L.J. Frohwein, U. Hofmann, C.H. Holland, D. Gonzáez, M. Keller, A. Seddek, T. Abbas, E.S.I. Mohammed, A. Teufel, T. Itzel, S. Metzler, R. Marchan, C. Cadenas, C. Watzl, M.A. Nitsche, F. Kappenberg, T. Luedde, T. Longerich, J. Rahnenführer, S. Hoehme, M. Trauner, and J.G. Hengstler.  Spatio-Temporal Multiscale Analysis of Western Diet-Fed Mice Reveals a Translationally Relevant Sequence of Events during NAFLD Progression.  _Cells_ , 10(10), 2021.  doi: 10.3390/cells10102516.
  * Hagemann and Möllenhoff [2024] N. Hagemann and K. Möllenhoff.  Overcoming model uncertainty–how equivalence tests can benefit from model averaging.  _arXiv preprint arXiv:2405.00827_ , 2024.  doi: 10.48550/arXiv.2405.00827.
  * Hothorn [2014] L.A. Hothorn.  Statistical evaluation of toxicological bioassays - a review.  _Toxicology Research_ , 3(6):418–432, 08 2014.  doi: 10.1039/c4tx00047a.
  * Huber et al. [2003] W. Huber, A. von Heydebreck, H. Sueltmann, A. Poustka, and M. Vingron.  Parameter estimation for the calibration and variance stabilization of microarray data.  _Statistical applications in genetics and molecular biology_ , 2003.  doi: 10.2202/1544-6115.1008.
  * Jensen et al. [2019] S.M. Jensen, F.M. Kluxen, and C. Ritz.  A review of recent advances in benchmark dose methodology.  _Risk Analysis_ , 39(10):2295–2315, 2019.  doi: 10.1111/risa.13324.
  * Kappenberg et al. [2021] F. Kappenberg, M. Grinberg, X. Jiang, A. Kopp-Schneider, J.G. Hengstler, and J. Rahnenführer.  Comparison of observation-based and model-based identification of alert concentrations from concentration-expression data.  _Bioinformatics_ , 37(14):1990 – 1996, 01 2021.  doi: 10.1093/bioinformatics/btab043.
  * Kawanaka et al. [2023] M. Kawanaka, K. Nishino, M. Kawada, K. Ishii, T. Tanikawa, R. Katsumata, N. Urata, J. Nakamura, M. Suehiro, K. Haruma, and H Kawamoto.  Soluble CD163 is a predictor of fibrosis and hepatocellular carcinoma development in nonalcoholic steatohepatitis.  _BMC gastroenterology_ , 23(10), 2023.  doi: 10.1186/s12876-023-02786-4.
  * Krug et al. [2013] A.K. Krug, R. Kolde, J.A. Gaspar, E. Rempel, N.V. Balmer, K. Meganathan, K. Vojnits, M. Baquié, T. Waldmann, R. Ensenat-Waser, et al.  Human embryonic stem cell-derived test systems for developmental neurotoxicity: a transcriptomics approach.  _Archives of toxicology_ , 87:123–143, 2013.  doi: 10.1007/s00204-012-0967-3.
  * Love et al. [2014] M.I. Love, W. Huber, and S. Anders.  Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.  _Genome Biology_ , 15:550, 2014.  doi: 10.1186/s13059-014-0550-8.
  * Love et al. [2020] M.I. Love, C. Soneson, P.F. Hickey, L.K. Johnson, N.T. Pierce, L Shepherd, M. Morgan, and R Patro.  Tximeta: Reference sequence checksums for provenance identification in RNA-seq.  _PLOS Computational Biology_ , 16:e1007664, 2020.  doi: 10.1371/journal.pcbi.1007664.
  * Love et al. [2024] M.I. Love, W. Huber, and S. Anders.  Analyzing RNA-seq data with DESeq2, 2024.  URL https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html.  (acessed 17.07.2024).
  * Möllenhoff et al. [2022] K. Möllenhoff, K. Schorning, and F. Kappenberg.  Identifying Alert Concentrations Using a Model-Based Bootstrap Approach.  _Biometrics_ , 79(3):2076 – 2088, 11 2022.  doi: 10.1111/biom.13799.
  * National Academies of Sciences, Engineering, and Medicine and others [2024] National Academies of Sciences, Engineering, and Medicine and others.  Charting a Future for Sequencing RNA and Its Modifications: A New Era for Biology and Medicine.  2024\.  Prepublication.
  * National Center for Biotechnology Information [a] National Center for Biotechnology Information.  CD163 molecule [Homo sapiens (human)], a.  URL https://www.ncbi.nlm.nih.gov/gene/93671.  (accessed 13th June 2024).
  * National Center for Biotechnology Information [b] National Center for Biotechnology Information.  Cd163 antigen [Mus musculus (house mouse)], b.  URL https://www.ncbi.nlm.nih.gov/gene/93671.  (accessed 13th June 2024).
  * National Center for Biotechnology Information [c] National Center for Biotechnology Information.  DBP D-box binding PAR bZIP transcription factor [Homo sapiens (human)], c.  URL https://www.ncbi.nlm.nih.gov/gene/1628.  (accessed 14th June 2024).
  * National Center for Biotechnology Information [d] National Center for Biotechnology Information.  Dbp D site albumin promoter binding protein [Mus musculus (house mouse)], d.  URL https://www.ncbi.nlm.nih.gov/gene/13170.  (accessed 14th June 2024).
  * [33] National Human Genome Research Institute.  The Human Genome Project.  URL https://www.genome.gov/human-genome-project.  (acessed 16.07.2024).
  * NCBI [1996] NCBI.  Pubmed, 1996.  URL https://pubmed.ncbi.nlm.nih.gov/?linkname=pubmed_pubmed_citedin&from_uid=25516281.  (acessed 15.07.2024).
  * Patro et al. [2017] R. Patro, G. Duggal, M. I. Love, R. A. Irizarry, and C. Kingsford.  Salmon provides fast and bias-aware quantification of transcript expression.  _Nature Methods_ , 14:417–419, 2017.  doi: 10.1038/nmeth.4197.
  * Pennisi [2024] E. Pennisi.  RNA deserves its own massive counterpart to the Human Genome Project, researchers argue.  _ScienceInsider_ , 2024.  doi: 10.1126/science.zj5d4ag.
  * Sutti and Albano [2020] S. Sutti and E. Albano.  Adaptive immunity: an emerging player in the progression of NAFLD.  _Nature Reviews Gastroenterology & Hepatology_, 17:81 – 92, 2020.  doi: 10.1038/s41575-019-0210-2.
  * Tibshirani [1988] R. Tibshirani.  Estimating Transformations for Regression via Additivity and Variance Stabilization.  _Journal of the American Statistical Association_ , 83(402):394 – 405, 1988.  doi: 10.1080/01621459.1988.10478610.
  * Wang et al. [2009] Z. Wang, M. Gerstein, and M. Snyder.  RNA-Seq: a revolutionary tool for transcriptomics.  _Nature reviews genetics_ , 10(1):57–63, 2009.  doi: 10.1038/nrg2484.
  * Zhu et al. [2018] A. Zhu, J. G. Ibrahim, and M. I. Love.  Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences.  _Bioinformatics_ , 35(12):2084 – 2092, 2018.  doi: 10.1093/bioinformatics/bty895.

## Appendix

######  Algorithm 1.

Let $\hat{\theta}$, $\hat{\sigma}^{2}$, $\tilde{S}$ and the time points of
interest be calculated from the original data set as described in Section 2.1.

  * Step 1 Generate bootstrap data at the $m$ different time points using $\hat{\theta}$ and $\hat{\sigma}^{2}$. Do so by drawing iid. errors $\varepsilon_{p,q}^{*}\sim N(0,1)$ and applying the formula

$y_{p,q}^{*}=f(t_{p},\hat{\theta})+\hat{\sigma}\varepsilon_{p,q}^{*}\text{ for
}p=1,\dots,m,\ q=1,\dots,n_{p}.$

Repeat this step $B_{3}$ times. This returns
$y_{p,q}^{*,1},\dots,y_{p,q}^{*,B_{3}}$.

Next, calculate the estimates $\hat{\theta}_{k}^{*}$ and
$\hat{\sigma}_{k}^{2,*}$, $k=1,\dots,B_{3}$, using the OLS methodology as
descried above for each of the generated data sets. This yields
$|f^{\prime}(t,\hat{\theta}_{k}^{*})|$.

  * Now for $k=1,\dots,B_{3}$ perform the two-step bootstrap procedure to estimate the confidence band:

  * Step 2 For $k=1,\dots,B_{3}$: Generate bootstrap data at the $m$ different time points using $\hat{\theta}^{*}_{k}$ and $\hat{\sigma}_{k}^{2,*}$. Do so by drawing iid. errors $\varepsilon_{p,q}^{**}\sim N(0,1)$ and applying the formula

$y_{p,q}^{**}=f(t_{p},\hat{\theta}^{*}_{k})+\hat{\sigma}_{k}^{*}\varepsilon_{p,q}^{*}\text{for
}p=1,\dots,m,\ q=1,\dots,n_{p}.$

Repeat this step $B_{1}$ times. This returns
$y_{{p,q}}^{**,1},\dots,y_{{p,q}}^{**,{B_{1}}}$.

For each combination of $k=1,\dots,B_{3}$ and $l=1,\dots,B_{1}$ estimate
$\hat{\theta}^{**}_{k,l}$ and $\hat{\sigma}^{2,**}_{k,l}$. This yields
$|f^{\prime}(t,\hat{\theta}^{**}_{k,l})|$.

The replicates
$|f^{\prime}(t,\hat{\theta}^{**}_{k,1})|,\dots,|f^{\prime}(t,\hat{\theta}^{**}_{{k,B_{1}}})|$
can be interpreted as the second level of the bootstrap. The standard error of
this sample yield the estimate
$\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}^{*}_{k})|}$.

  * Step 3 For each combination of $k=1,\dots,B_{3}$ and $l=1,\dots,B_{1}$: Generate bootstrap data at the $m$ different time points using $\hat{\theta}^{**}_{k,l}$ and $\hat{\sigma}_{k,l}^{2,**}$. Do so by drawing iid. errors $\varepsilon_{p,q}^{***}\sim N(0,1)$ and applying the formula

$y_{p,q}^{***}=f(t_{p},\hat{\theta}^{**}_{k,l})+\hat{\sigma}_{k,l}^{**}\varepsilon_{p,q}^{***}\text{for
}p=1,\dots,m,\ q=1,\dots,n_{p}.$

Repeat this step $B_{2}$ times.

Estimate $\hat{\theta}^{***}$ (indices omitted). This yields $B_{2}$
replicates of the form $|f^{\prime}(t,\hat{\theta}^{***})|$. The standard
errors of the $B_{2}$ replicates yield estimates
$\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}^{**}_{k,l})|}$.

  * Step 4 For each combination of $k=1,\dots,B_{3}$ and $l=1,\dots,B_{1}$: Calculate

$D^{*,k,l}:=\max_{t\in
T}\frac{|f^{\prime}(t,\hat{\theta}_{k,l}^{**})|-|f^{\prime}(t,\hat{\theta^{*}_{k}})|}{\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}_{k,l}^{**})|}}.$

The respective empiric $(1-\alpha)$-quantiles of the distribution yield
estimates of the critical values $c_{k}^{*}$.

  * Step 5 For $k=1,\dots,B_{3}$: Calculate the lower simultaneous confidence band

$L^{\alpha}(t,\hat{\theta}^{*}_{k})=|f^{\prime}(t,\hat{\theta}_{k}^{*})|+c_{k}^{*}\hat{\sigma}_{|f^{\prime}(t,\hat{\theta}^{*}_{k})|}.$

Calculate the time frames of significant change in gene expression:

$\tilde{S}_{k}\coloneqq\\{t|L^{\alpha}(t,\hat{\theta}^{*}_{k})>\lambda\\}$

Partition the $\tilde{S}_{k}$ into coherent subsets:

$\tilde{S}_{k}=(S_{k,1},\dots,S_{k,s})$

  * Step 6 For $k=1,\dots,B_{3}$: Dismiss the cases with more or less coherent subsets then the original data set. Identify the time points of interest for the remaining runs. Now, calculate the empiric $\frac{\alpha}{2}$\- and $(1-\frac{\alpha}{2})$-quantiles ($Q_{\frac{\alpha}{2}}$ and $Q_{(1-\frac{\alpha}{2})}$) for each time point of interest, respectively. The resulting confidence intervals are:

$\left[Q_{\frac{\alpha}{2}},Q_{(1-\frac{\alpha}{2})}\right]$