# From 2D to 3D: Re-thinking Benchmarking of
Monocular Depth Prediction

Evin Pınar Örnek1,∗, Shristi Mudgal1,∗, Johanna Wald1,2, Yida Wang1, Nassir
Navab1 and Federico Tombari1,2 1:Technical University of Munich, Germany;
{evin.oernek, shristi.mudgal, yida.wang, nassir.navab,
federico.tombari}@tum.de2:Google Inc. johannawald@google.com$*$:the first two
authors contributed equally.

###### Abstract

There have been numerous recently proposed methods for monocular depth
prediction (MDP) coupled with the equally rapid evolution of benchmarking
tools. However, we argue that MDP is currently witnessing benchmark over-
fitting and relying on metrics that are only partially helpful to gauge the
usefulness of the predictions for 3D applications. This limits the design and
development of novel methods that are truly aware of - and improving towards
estimating - the 3D structure of the scene rather than optimizing 2D-based
distances. In this work, we aim to bring structural awareness to MDP, an
inherently 3D task, by exhibiting the limits of evaluation metrics towards
assessing the quality of the 3D geometry. We propose a set of metrics well
suited to evaluate the 3D geometry of MDP approaches and a novel indoor
benchmark, RIO-D3D, crucial for the proposed evaluation methodology. Our
benchmark is based on a real-world dataset featuring high-quality rendered
depth maps obtained from RGB-D reconstructions. We further demonstrate this to
help benchmark the closely-tied task of 3D scene completion.

Figure 1: State of the Art depth prediction with low absolute relative error
(absrel) often does not exhibit correct 3D structures when back-projected in
3D space (b). We suggest using 3D metrics, e.g., chamfer distance (CD), to
measure the structural quality of the prediction and propose a new depth
prediction benchmark (d) with entire 3D geometry and high-quality depth
compared to raw sensor measurements (a, c). Notice invalid depth regions of
the sensor depth maps shown in black.

##  I INTRODUCTION

Recent computer vision literature has witnessed an abundance of research on
monocular depth prediction (MDP), i.e., the task of predicting a dense depth
map from a single RGB image. Compared to traditional approaches [1, 2] which
were proposed more than a decade ago, the advent of deep learning in the last
few years offers more effective ways [3, 4, 5, 6] to regress depth maps if a
vast set of data is available for training. However, the literature is
trending towards over-fitting the development of novel ideas to the
maximization of a small subset of specific metrics on these datasets. This
poses a limit to the development of more precise and more useful MDP
approaches for two reasons: i) the main metrics tend to disregard the accuracy
for the 3D structure of the predicted geometry, so the depth maps scoring high
values for these metrics do not reflect correct 3D structures when they are
back-projected in 3D space (see Fig. 1(a, b)); current metrics compute the
accuracy globally and up to the 2D depth plane rather than in the 3D domain
ii) some datasets are supplemented with noisy 3D supervisions which are
represented by the raw depth data obtained from RGB-D sensors, such that
training and evaluating will be biased by noise in the ground truth. This is
portrayed by the examples shown in Fig. 1(a, c). We argue that improving 3D
accuracy of predictions is of utmost importance for MDP – since the vast
majority of applications using depth prediction relies on 3D geometry of the
scene and the objects therein – in order to, e.g., navigate unseen
environments through obstacle avoidance, grasp objects within clutter or
augment shapes with 3D information. However, the accuracy of predicted 3D
geometry, is not sufficiently represented in the commonly used 2D-based
metrics nor in the network designs and loss functions being proposed. This
results in continual research in a direction that does not allow the true
growth of MDP in terms of the quality of predictions.

Based on this preamble, the goal of our work is to move beyond this local
minimum by exploring novel avenues to develop, evaluate and benchmark MDP in a
more 3D-driven way, such that its applications could be well served.
Furthermore, we noticed that the current MDP datasets are not supplemented
with the necessary size, quality or semantics, and a 3D-based evaluation
requires depth maps of high quality on 3D space. Hence, we propose a large-
scale indoor dataset RIO-D3D (see Fig. 1(d))) where aside from semantic
information, a high quality, complete and denoised 3D ground truth is
available, resulting in a cleaner and finer structured, dense point cloud. Our
dataset also enables evaluating the task of semantic scene completion, which
is directly related to MDP, as both aim towards a full 3D reconstruction of a
scene. The proposed benchmark dataset and evaluation codes will be publicly
available.

In summary, our contributions are the following:

  1. 1.

We propose a comprehensive evaluation protocol based on 3D metrics, aiming at
assessing the prediction quality in terms of 3D geometry.

  2. 2.

A novel large-scale, real, indoor benchmark (see Fig. 1(d)) dataset RIO-D3D is
provided for MDP with highly accurate, complete 3D ground truth, together with
normal, occlusion and semantic annotations; allowing training and evaluation
of more suited, 3D aware approaches.

  3. 3.

We provide a comparison among state of the art MDP and 3D scene completion
methods along with the proposed evaluation criteria on RIO-D3D dataset.

  4. 4.

We bring to light the need to shift the focus of MDP research towards the
structural information in the scene and better quantify the prediction
quality. This is affirmed by experimental results, highlighting the flaws of
2D metrics and 2D-based approaches.

##  II RELATED WORK

Monocular depth prediction Initially, traditional works on MDP [1, 7, 8, 9,
10] build the inference models based on hand-crafted features tied to the
depth priors in the scenes. Deep learning provides more potential when large
scale training data is available, the very first approach [3] proposed is
based on AlexNet [11]/VGG [12]. [13] improved the performance by stacking more
layers and predicting the normals and semantic labels. Laina et al. [4]
introduced a ResNet-based fully convolutional network with novel up-sampling
blocks while other works [14, 15, 16] combined multi-scale convolutional
neural network features with a Conditional Random Field refinement. Several
methods attempted to model MDP as a relative ordering task. Fu et al. [17]
applied ordinal regression, Dhamo et al. [18] used a layer-wise scene
modeling, Lee et al. [19] proposed a model for ordinal consistency. The heavy-
tailed distribution in depth maps, where large amounts of pixels are
predominantly small depth values, leads to data imbalance problem. Reversed
Huber loss [4], binary decision trees [20], attention weighted loss [21] were
proposed for this. Recently, Bhat et al. [22] proposed a novel formulation for
predicting distance values by means of classification. Tian et al. [23]
proposed attention blocks within the decoder, while other methods adopt a
Transformer based architecture [24, 25].

Further, multi-task learning is being actively explored, focusing on jointly
learning depth prediction with tasks such as surface normals or semantic
segmentation. Some of the works [26, 27, 28] aim at learning a globally
consistent depth map but do not focus on fine local details, often resulting
in blurry boundaries and deformed planar surfaces. Several works [29, 30, 31,
32, 33] focus on accurate object boundaries. Ramamonjisoa et al. [5] aim to
improve the edges by predicting normals and edges along with depth. In
general, multi-task methods need complementary data in addition to the depth
ground truth which is often missing. Our proposed benchmark provides
semantics, normal maps, as well as occlusion boundaries for the entire
dataset.

3D-aware depth prediction Most of the aforementioned methods use 2D networks
and predict depth using a structured 2D pixel grid. These works are trending
towards over-fitting to a subset of commonly used 2D metrics. Only recently
have approaches been proposed that bring more 3D awareness to MDP. Jointly
estimating 3D surface normals [34, 6] or enforcing planar constraints [35]
improves depth predictions in indoor setups. Liu et al.reconstruct a scene
from an arbitrary number of piece-wise planar surfaces [35]. The first work to
highlight the role of depth maps in understanding 3D scene geometry is Yin et
al.[36]. They introduce a 3D training loss which backprojects the depth to 3D
and enforces consistency on 3D surface patches. While these methods present
convincing performance, they are not only limited by the quality of the
training data but also by the metrics used when evaluating.

TABLE I: Comparison of indoor depth prediction datasets.

Dataset | Sem. | Full 3D | Setting | Acquisition | Size
---|---|---|---|---|---
NYU v2 [37] | ✓ | ✗ | indoor | Kinect | 1449
SUN RGB-D [38] | ✓ | ✓ | indoor | various | 10K
DIW [39] | ✗ | ✗ | both | web (relative) | 500K
Matterport [40] | ✓ | ✓ | indoor | Matterport | 200K
2D-3D-S [41] | ✓ | ✓ | indoor | Matterport | 70K
ScanNet [42] | ✓ | ✓ | indoor | structured | 24K
ETH 3D [43] | ✗ | ✗ | both | laser | 545
DIODE [44] | ✗ | ✗ | both | Faro | 27K
RIO-D3D | ✓ | ✓ | indoor | 3D scan | 34K

Datasets and metrics Numerous datasets have been introduced for data-driven
deep learning models, see Tbl. I. The earlier benchmarks are afflicted by
either low resolution, sparse depth, small size and/or sensor limitations of
the time. NYUv2 [37] is by far the most popular indoor benchmark featuring
real-world environments but includes only a small number of semantic images,
inline with other datasets that were later introduced with semantic
annotations [38, 41, 40]. ScanNet [42] while providing 3D reconstructions,
does not provide 3D reconstructed depth maps. It is a sensor based depth
benchmark with 24K images which, due to the sensor, similarly fails on low
albedo and reflective surfaces. As opposed to these works, ETH [43] or DIODE
[44] provide finer ground truth but lack semantics and availability of the
full 3D scene including occluded regions (see Fig. 1(d)) in addition to DIODE
being relatively small.

The evaluation on these datasets is commonly carried out by using a set of
metrics: absolute relative distance (Eq. 1), squared relative distance (Eq.
2), root mean square error (Eq. 3), logged root mean square error (Eq. 4) and
accuracy under threshold (Eq. 5).

$\text{absrel}=\frac{1}{|T|}\sum_{y\in T}\frac{|y-y^{*}|}{y^{*}},$ (1)
$\text{sqrel}=\frac{1}{|T|}\sum_{y\in T}\frac{{\|y-y^{*}\|}^{2}}{y^{*}},$ (2)
$\text{rmse (linear)}=\sqrt{\frac{1}{|T|}\sum_{y\in T}{\|y-y^{*}\|}^{2}},$ (3)
$\text{rmse (log)}=\sqrt{\frac{1}{|T|}\sum_{y\in T}{\|\log y-\log
y^{*}\|}^{2}},$ (4) $\text{accuracy at $\delta$,}=\%\text{ of }y_{i}\text{
s.t. }max(\frac{y_{i}}{y_{i}^{*}},\frac{y_{i}^{*}}{y_{i}})=\delta<t.$ (5)

Given a ground truth depth map and the corresponding prediction, metrics are
calculated iterating over all pixels in dataset $T$, comparing the ground
truth value $y^{*}$ with its predicted counterpart $y$, and averaging over all
data samples $|T|$.

In addition to standard 2D metrics, Koch et al. [45] recently suggested an
occlusion boundary metric ${acc}_{DBE}$ to measure the occlusion boundaries on
depth maps. Furthermore, they propose plane accuracy metric to measure the
connectivity and geometry of planar surfaces. Inspired from them, Ramamonjisoa
et al.[5] also emphasized the occlusion boundary metric. However, none of
these works evaluate depth maps in terms of 3D structures. To the best of our
knowledge, our work is the first to explore the usage of different 3D metrics
when evaluating MDP.

##  III EVALUATION METHODOLOGY

Figure 1(b) shows the predicted point cloud obtained using an available pre-
trained model. We can observe consistent deformations and noise, especially
around the occlusion boundaries. This can be attributed to the 2D-oriented
network architectures, datasets, 2D depth value based loss functions, and
deficient evaluation criteria. This work addresses two of these aspects:
evaluation criteria and datasets.

To that end, we redefine MDP as a single-view 3D reconstruction task. Inspired
by the related 2D-to-3D tasks, we propose the usage of 3D metrics and losses
to complete the missing 3D aspects within MDP. While most of these metrics
have already been proposed earlier [46, 47], we are the first to apply them
when evaluating depth prediction.

#### Earth Mover’s Distance (EMD)

Consider two point clouds $G,R\subseteq\mathbb{R}^{3}$ of equal size
$|G|=|R|$, where $G$ is the ground truth, $R$ is the reconstruction. EMD
between $G$ and $R$ is defined as follows:

$d_{EMD}(G,R)=\min_{\phi:G\rightarrow R}\sum_{x\in G}\|x-\phi(x)\|_{2},$ (6)

where $\phi:G\rightarrow R$ is a bijection. Practically, finding this
bijection mapping is too expensive and we follow an approximate implementation
as in [46].

#### 3D Chamfer Distance (CD)

$d_{CD}(G,R)$ between $G,R\subseteq\mathbb{R}^{3}$ is defined as:

$\displaystyle d_{CD}(G,R)=\sum_{x\in G}\min_{y\in
R}\|x-y\|^{2}_{2}+\sum_{y\in R}\min_{x\in G}\|x-y\|^{2}_{2},$ (7)

where each point from $G$ is matched with the nearest neighboring point from
$R$ and vice versa. CD is known to be sensitive to outliers and potentially
finds wrong correspondences when points are incorrectly matched.

#### 3D Completeness

This score closes the gap between EMD and CD. EMD is typically calculated from
ground truth to prediction, whereas CD covers both directions simultaneously.
Completeness measures the distances from $G$ to $R$, showing the error on the
retrieved points. It can be calculated as:

$\displaystyle d_{EMD}(R,G)=\min_{\phi:R\rightarrow G}\sum_{x\in
R}\|x-\phi(x)\|_{2},$ (8)

Completeness is necessary to understand whether the reconstructed points lie
within an acceptable the range of ground truth points.

Figure 2: Two samples from RIO-D3D a) RGB Images, corresponding b) depth map
(red is far, blue is close), c) back-projected depth, d) instance
segmentation, e) normals, and f) instance and occlusion boundaries.

#### 3D Point Cloud F-score

F-score, also known as Dice coefficient, is the harmonic mean of precision and
recall. Precision is the number of predicted points lying within a distance to
a ground truth surface, whereas recall is the number of ground truth points
lying within the prediction. The metric was introduced by Knapitsch [48] to
evaluate 3D surface quality, and further emphasized by Tatarchenko et al.[47].

The recall with a distance threshold $t$ is defined as:

$\displaystyle R(t)=\frac{1}{|G|}\sum_{g\in G}[e_{g\rightarrow R}<t].$ (9)

where the distance $e$ from a ground truth point $g$ to a reconstructed point
$r$ is $e_{g\rightarrow R}=\min_{r\in R}\|g-r\|.$

Likewise, the precision is:

$\displaystyle P(t)=\frac{1}{|R|}\sum_{r\in R}[e_{r\rightarrow G}<t].$ (10)

where the distance $e$ from a reconstructed point $r\in R$ to ground truth $G$
is defined as $e_{r\rightarrow G}=\min_{g\in G}\|r-g\|.$

$\displaystyle F(t)=2\cdot\frac{P(d)\cdot R(d)}{P(d)+R(d)}.$ (11)

#### 3D Point Cloud Intersection Over Union

IoU is commonly used to measure the volumetric quality of a 3D reconstruction.
It is also known as Jaccard Index calculated between two input sets. 3D IoU in
literature is computed on voxelized point clouds, which requires a sampling
and discretization step, losing the ability to reflect surface quality. We
instead propose an IoU calculation on non-voxelized point clouds. Inspired by
the generalized formulation of the F-score, Jaccard and the Tversky index
[49], and their relationship formulated as
$\frac{J}{F}=\frac{1}{2}+\frac{J}{2}$, we re-express the point based IoU in
terms of precision and recall:

$IoU(d)=\frac{P(d)\cdot R(d)}{P(d)+R(d)-P(d)\cdot R(d)}.$ (12)

We compare these 3D-based metrics on the proposed dataset to evaluate MDP in
the experiments section.

##  IV EVALUATION DATASET: RIO-D3D

We propose RIO-D3D, a new large-scale, highly variable, an indoor dataset with
calibrated RGB-D sequences, high-quality 3D ground truth, instance-level
semantic annotations, associated normals, and geometric as well as occlusion
boundaries. We build upon the 3RScan dataset [50] which features 1482 RGB-D
scans of 478 unique environments, camera poses, and intrinsic parameters. They
processed each sequence offline to get bundle-adjusted camera poses with loop-
closure and texture-mapped 3D reconstructions. The original dataset has 363K
real images, yet, the provided 3D reconstruction allows the generation of
unlimited semi-synthetic data from various viewpoints to create additional
training data.

The dense 3D reconstruction reduces missing surfaces and mitigates noise
(particularly on edges and in low lighting environments) typically present in
sensor-acquired depth maps. We render depth maps associated with the RGB
frames using the camera poses and the intrinsic data allowing access to a more
complete depth map with higher coverage and rich 3D structural information.
This results in a ground truth of superior quality than common sensors, making
the benchmark suitable for training and evaluating different kinds of
approaches.

A good depth prediction dataset should provide not only RGB-depth pairs in a
large-scale setting but also semantic and additional geometric information to
help propel multi-task learning. When compared to the existing datasets,
RIO-D3D provides not just high-quality ground truth but also normal maps,
occlusion boundaries, semantic and instance labels for the entirety of our
dataset. The normal maps and occlusion boundaries are generated by following a
point cloud-based approach from Tateno et al.[51].

(a) Depth dist. (b) Coverage

Figure 3: Depth distribution (in mm over number of pixels) and dataset coverage are shown.  Method | 2D Metrics | 3D Metrics
---|---|---
absrel$\downarrow$ | rmse$\downarrow$ | $\delta_{1}\uparrow$ | CD$\downarrow$ | EMD$\downarrow$ | Comp.$\downarrow$ | IoU$\uparrow$ | F-score$\uparrow$
Oracle NN | 0.097 | 0.225 | 0.891 | 0.257 | 0.130 | 0.128 | 0.328 | 0.335
Median Plane | 0.211 | 0.577 | 0.668 | 0.677 | 0.636 | 0.042 | 0.347 | 0.369
Eigen [3] | 0.217 | 0.712 | 0.637 | 0.584 | 0.529 | 0.055 | 0.254 | 0.405
FCRN [4] | 0.217 | 0.703 | 0.647 | 0.491 | 0.430 | 0.061 | 0.273 | 0.428
BTS [6] | 0.190 | 0.657 | 0.694 | 0.454 | 0.400 | 0.055 | 0.336 | 0.500
VNL [36] | 0.258 | 0.638 | 0.534 | 0.764 | 0.686 | 0.078 | 0.219 | 0.313
TABLE II: Quantitative comparison of MDP methods on standard 2D-based metrics
and the proposed 3D-based metrics. The first two lines show the metric results
for Oracle Nearest Neighbour and Median Plane experiments.

Dataset details We focus on the real images in the experiments to provide a
fair comparison to the state of the art. After preprocessing, we selected a
benchmark of 34K images (8:1:1 train, val and test). This is to provide a
content-rich subset with minimal redundancy, considering that larger datasets
are often not used in their entirety due to the limitation of train time, for
example, the KITTI-Eigen split (24K images) from the full KITTI dataset. Depth
distribution and coverage statistics can be seen in Fig. 3.

##  V EXPERIMENTS

This section will first lay the groundwork to establish the need for a new
evaluation protocol through two toy experiments. Then we give the qualitative
and quantitative results for our MDP and scene completion benchmarks with the
RIO-D3D and proposed evaluation methodology.

###  V-A What do deep learning based MDP networks learn?

We carry out two experiments to examine the behavior of MDP under a
hypothetical setup and gain insight into how they compare to the performance
of existing methods:

Figure 4: Sample (a) Median Plane (b) Oracle Nearest Neighbour test image
(top-left), and corresponding retrieved closest train image (top-right) with
their depth maps (below).

Median plane We argue that the existing MDP methods learn an average planar
depth due to using 2D-based distances and a global averaging. To test this, we
handcraft a set of depth maps, where the distance values are set to the ground
truth depth map’s median distance value, visualized in Fig. 4(a). We use these
planar depth maps as predictions and evaluate them with both 2D and 3D-based
metrics. These handcrafted depth maps gave an absrel of $0.211$, which is on
par with the model predictions, as shown in Tbl. II. This confirmed that if a
network were to simply learn to associate median depth values and minimize the
loss instead of trying to learn structures, they would still get reasonable 2D
metrics results. However, they would struggle on the 3D-based metrics (CD of
$0.677$ and F-score of $0.369$).

Oracle nearest neighbour This experiment measures the performance of a
retrieval-driven MDP. How would a method perform that, instead of regressing
the distance values, retrieves the structurally nearest image in the training
set? To confirm this, for each validation image, we retrieved the nearest
depth map from the train split according to the absrel value (see Fig. 4(b)).
We then calculated the overall metrics by assuming the nearest depth maps as
predictions for the images on the validation set. Results (Tbl. II) show that
OracleNN is indeed superior to the predictions, so retrieving a depth map from
the training set gives better quantitative performance than existing
regression methods.

These analysis experiments show that MDP methods are not performing as well as
they should. We do not see a true improvement in terms of the structure and
geometry of the predictions. It raises the fundamental question of what
networks can truly learn if they are correctly evaluated, and we should re-
think how MDP can be better benchmarked.

Figure 5: Qualitative depth prediction results on 3 different samples with a)
RGB input and corresponding b) ground truth depth, next to the prediction of
c) Eigen, d) FCRN, e) BTS and f) VNL.

###  V-B Monocular Depth Prediction Benchmarking

We provide a benchmark for RIO-D3D through selected MDP methods: Eigen [3],
FCRN [4], BTS [6] and VNL [36]. When training the models, we follow the
original implementations and hyperparameters in terms of batch size, learning
rate, input image size and augmentation.

Evaluation based on 2D metrics Quantitative results are provided in Tbl. II.
It can be observed that the standard 2D metrics are saturated and vary with an
insignificant difference. As can be seen in the qualitative results (Fig. 5)
that show point clouds obtained from back-projected predictions, methods
cannot predict 3D structures well, and the standard 2D metrics do not capture
this. The metric scores do not align well with the observed prediction
quality. Despite the good metric scores, BTS prediction point clouds are full
of outliers and fail to capture the 3D object surface.

This highlights the need for methods and metrics which focus on 3D aspects to
better quantify the quality of predictions that aligns with visual
observations.

Evaluation based on 3D metrics Eigen [3] and FCRN [4] have same absrel scores
with $0.21$, however, FCRN is better in terms of 3D metrics which is also
evidenced by point cloud visualizations on Fig. 5. Likewise, a recent method
BTS lies within the similar range on 2D metrics with older methods Eigen and
FCRN, yet again, as it can be seen on the point clouds, it is structurally
superior than others, which is also indicated with the best F-score on the
Tbl. II. When it comes to VNL, it performs worse than the other methods. The
standard metrics show only a slight change of $\%16$ ($0.25$ absrel compared
to $0.21$ FCRN absrel), whereas 3D metrics reflect the quality better by
showing a higher change $\%35$ ($0.76$ CD compared to $0.49$ FCRN CD).

Figure 6: Predictions with relatively good absrel $\downarrow$ (top right) and
bad F-score $\uparrow$ (bottom right). Error map (f) colors indicate
$\blacksquare$ ground truth $\blacksquare$ prediction $\blacksquare$
intersection. Figure 7: Scene completion results with 2,048 input partial
point cloud, GRNet [52], MSN [53], ForkNet [54], complete ground truth point
cloud.

Comparing 2D metrics to 3D metrics The current implementations of 3D based
metrics follow the trend of standard metrics closely (as in Tbl II). If one
method is the best on standard metrics, it is also the comparably better on 3D
metrics. Yet, when looking at the individual samples, there are numerous
conflicting examples, i.e., good absrel value but bad F-score or vice versa,
as in Figures 1(a,b) and 6. Note that the first sample in Fig. 6 has an absrel
value comparable to the test average but a poor F-score. Similarly, the second
sample with an absrel that is on-par with test average, however, has a poor 3D
structure, which is also reflected by the F-score.

The major advantage of the 3D metrics is that they align better with observed
qualitative results and quantify the quality of the structures present in the
scene more accurately. Even when the intersection between ground truth and
prediction is small, absrel values show low variance. The F-score values
however quickly get worse when structural overlap is reduced. 2D metrics
evaluate averaged per pixel depth error which can be small even if the depth
is highly inaccurate. For poor absrel, the 3D metrics are also poor. However,
for a small range of absrel values, 3D scores have a larger variety and better
reflect the quality of 3D structure.

Comparison between 3D metrics CD and EMD are commonly used for evaluating 3D
shape reconstruction. Yet, they depend on k-nearest neighbour matching which
can lead to unreliable scores due to outliers. Comparably, point cloud based
IoU and F-score are more robust, because even when the knn matching is wrong,
IoU and F-score do not magnify the effect through falsely identified
distances. They show the quality of 3D shapes more clearly, e.g., the
interpolation effect on occlusion boundaries and poor planar surfaces. The
range of IoU/F-score values depend on the selected threshold, in our
experiments set as 0.01 meters, and smaller thresholds ($\leq 10cm$) are too
strict given potential noise levels of depth sensors and the ground truth. IoU
is more sensitive than F-score and gives lower values for poor depth maps,
though beyond a threshold, results in sudden high spikes.

Recommended metric Even though we invite the community to take a critical
standpoint, F-score shows to be more advantageous over the others. Calculating
the harmonic mean between precision and recall, it measures the percentage of
3D points reconstructed correctly. It is more reliable than CD/EMD/Comp., more
robust than IoU, and considers both the structural quality and the amount of
correctly retrieved structures. Hence, F-score is the most promising 3D
indicator amongst the proposed ones for evaluating MDP.

###  V-C Scene Completion Benchmarking

RIO-D3D comes with the availability of full 3D scenes along with the partial
depth maps. This can serve as benchmark for the 3D single-view scene
completion [55, 56, 52], an active research field which recently witnessed a
lack of benchmarking tools due to the difficulty of obtaining completion
ground truth and the disappearance of one of the main datasets (SUNCG [57]).
Towards this aim, we extend our benchmark to the completion task, where the
partial scene input is based on the depth dataset, and the completed scene
input is extracted from the full 3D mesh. We provide a comparison among state
of the art scene completion methods in Tbl. III, where a 10-frame qualitative
baseline for semantic completion are shown in Fig. 7. VRCNet [58] and
SoftPoolNet [56] are applied in a cascaded way for 3D completion and semantic
segmentation separately. It can be seen that the point-cloud based methods
(GRNet [52] and MSN [53]) produce noisy results as they specifically focus on
completion tasks where samples do not vary significantly in both geometries
and poses. However, voxel based approaches (ForkNet [54]) work better both in
terms of CD and EMD, and qualitatively. In general, our scenes have a wider
variety, more diverse geometry and depths than the existing scene completion
datasets.

TABLE III: Evaluation of the scene completion with CD ($\times 10^{3}$) and EMD ($\times 10^{4}$) on RIO-D3D. Method | CD$\downarrow$ | EMD$\downarrow$
---|---|---
FoldingNet [55] | 18.64 | 11.78
PCN [59] | 13.18 | 10.63
MSN [53] | 8.72 | 5.31
SoftPoolNet [56] | 10.07 | 5.76
GRNet [52] | 7.55 | 5.20
VRCNet [58] | 6.24 | 4.92

##  VI CONCLUSION

It is crucial to rethink the MDP task regarding structural accuracy rather
than a global pixel-wise 2D error because a depth map, though represented in
2D format, actually consists of 3D information about the environment. As
current research fails to address this issue suitably, more appropriate
methods, losses, and metrics are imperative. We aim to enable novel ideas to
capture scene geometry finer than existing methods. To this end, we point out
3D metrics that directly evaluate the 3D structures of depth maps and provide
a suited depth prediction benchmark with complete 3D structures and high-
quality depth maps. Our experiments show F-score is the most promising
indicator. RIO-D3D and our evaluation codes will be publicly available. Having
brought this forth, we hope that future research on MDP will do better than
retrieving planes and images.

## References

  * [1] D. Hoiem, A. A. Efros, and M. Hebert, “Automatic photo pop-up,” _ACM ToG_ , pp. 577––584, 2005.
  * [2] A. Saxena, M. Sun, and A. Y. Ng, “Make3D: Learning 3D Scene Structure from a Single Still Image,” in _IEEE TPAMI_ , 2009.
  * [3] D. Eigen, C. Puhrsch, and R. Fergus, “Depth map prediction from a single image using a multi-scale deep network,” in _NeurIPS_ , 2014.
  * [4] I. Laina, C. Rupprecht, V. Belagiannis, F. Tombari, and N. Navab, “Deeper depth prediction with fully convolutional residual networks,” in _3DV_ , 2016\.
  * [5] M. Ramamonjisoa and V. Lepetit, “Sharpnet: Fast and accurate recovery of occluding contours in monocular depth estimation,” _ICCV Workshops_ , 2019\.
  * [6] J. H. Lee, M.-K. Han, D. W. Ko, and I. H. Suh, “From big to small: Multi-scale local planar guidance for monocular depth estimation,” _arXiv preprint arXiv:1907.10326_ , 2019.
  * [7] C. Liu, J. Yuen, and A. Torralba, “Sift flow: Dense correspondence across scenes and its applications,” _IEEE TPAMI_ , vol. 33, no. 5, pp. 978–994, 5 2011.
  * [8] A. Oliva and A. Torralba, “Modeling the shape of the scene: A holistic representation of the spatial envelope,” _IJCV_ , vol. 42, no. 3, pp. 145–175, 2001.
  * [9] A. Saxena, S. H. Chung, and A. Y. Ng, “Learning depth from single monocular images,” in _NeurIPS_ , 2005.
  * [10] K. Karsch, C. Liu, and S. B. Kang, “Depth extraction from video using non-parametric sampling,” in _ECCV_ , 2012.
  * [11] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _NeurIPS_ , F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, Eds., vol. 25, 2012.
  * [12] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” in _ICLR_ , 2015.
  * [13] D. Eigen and R. Fergus, “Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture,” in _ICCV_ , 2015.
  * [14] Wei Zhuo, M. Salzmann, X. He, and M. Liu, “Indoor scene structure analysis for single image depth estimation,” in _CVPR_ , 2015.
  * [15] Bo Li, C. Shen, Yuchao Dai, A. van den Hengel, and Mingyi He, “Depth and surface normal estimation from monocular images using regression on deep features and hierarchical crfs,” in _CVPR_ , 2015.
  * [16] Peng Wang, Xiaohui Shen, Zhe Lin, S. Cohen, B. Price, and A. Yuille, “Towards unified depth and semantic prediction from a single image,” in _CVPR_ , 2015.
  * [17] H. Fu, M. Gong, C. Wang, K. Batmanghelich, and D. Tao, “Deep Ordinal Regression Network for Monocular Depth Estimation,” in _CVPR_ , 2018.
  * [18] H. Dhamo, N. Navab, and F. Tombari, “Object-driven multi-layer scene decomposition from a single image,” in _ICCV_ , 2019.
  * [19] J. Lee and C. Kim, “Monocular depth estimation using relative depth maps,” in _CVPR_ , 2019.
  * [20] A. Roy and S. Todorovic, “Monocular depth estimation using neural regression forest,” _CVPR_ , pp. 5506–5514, 2016.
  * [21] Y. S. Jianbo Jiao, Ying Cao and R. Lau, “Look deeper into depth: Monocular depth estimation with semantic booster and attention driven loss,” in _ECCV_ , 2018.
  * [22] S. Bhat, I. Alhashim, and P. Wonka, “Adabins: Depth estimation using adaptive bins,” _CVPR_ , 2021.
  * [23] T. Chen, S. An, Y. Zhang, C. Ma, H. Wang, X. Guo, , and W. Zheng, “Improving monocular depth estimation by leveraging structural awareness and complementary datasets,” _ECCV_ , 2020.
  * [24] R. Ranftl, A. Bochkovskiy, and V. Koltun, “Vision transformers for dense prediction,” _ICCV_ , 2021.
  * [25] G. Yang, H. Tang, M. Ding, N. Sebe, and E. Ricci, “Transformer-based attention networks for continuous pixel-wise prediction,” in _ICCV_ , 2021.
  * [26] J. Li, R. Klein, and A. Yao, “A two-streamed network for estimating fine-scaled depth maps from single rgb images,” in _ICCV_ , 2017.
  * [27] O. Hosseini Jafari, O. Groth, A. Kirillov, M. Ying Yang, and C. Rother, “Analyzing modular cnn architectures for joint depth prediction and semantic segmentation,” in _ICRA_ , 2017.
  * [28] A. Mousavian, H. Pirsiavash, and J. Kosecka, “Joint semantic segmentation and depth estimation with deep convolutional networks,” in _3DV_ , 2016.
  * [29] J. Hu, M. Ozay, Y. Zhang, and T. Okatani, “Revisiting single image depth estimation: Toward higher resolution maps with accurate object boundaries,” in _IEEE WACV_ , 2019.
  * [30] M.-L. Shih, S.-Y. Su, J. Kopf, and J.-B. Huang, “3d photography using context-aware layered depth inpainting,” in _CVPR_ , 2020.
  * [31] S. Niklaus, L. Mai, J. Yang, and F. Liu, “3d ken burns effect from a single image,” _ACM TOG_ , vol. 38, no. 6, pp. 184:1–184:15, 2019.
  * [32] M. Ramamonjisoa, Y. Du, and V. Lepetit, “Predicting sharp and accurate occlusion boundaries in monocular depth estimation using displacement fields,” in _CVPR_ , 2020.
  * [33] L. Wang, J. Zhang, Y. Wang, H. Lu, and X. Ruan, “Cliffnet for monocular depth estimation with hierarchical embedding loss,” in _ECCV_ , 2020.
  * [34] Z. Zhang, Z. Cui, C. Xu, Y. Yan, N. Sebe, and J. Yang, “Pattern-affinitive propagation across depth, surface normal and semantic segmentation.” in _CVPR_ , 2019.
  * [35] C. Liu, K. Kim, J. Gu, Y. Furukawa, and J. Kautz, “Planercnn: 3d plane detection and reconstruction from a single image,” in _CVPR_ , 2019.
  * [36] W. Yin, Y. Liu, C. Shen, and Y. Yan, “Enforcing geometric constraints of virtual normal for depth prediction,” in _ICCV_ , 2019.
  * [37] P. K. Nathan Silberman, Derek Hoiem and R. Fergus, “Indoor segmentation and support inference from rgbd images,” in _ECCV_ , 2012.
  * [38] S. Song, S. P. Lichtenberg, and J. Xiao, “Sun rgb-d: A rgb-d scene understanding benchmark suite,” in _CVPR_ , 2015.
  * [39] W. Chen, Z. Fu, D. Yang, and J. Deng, “Single-image depth perception in the wild,” in _NeurIPS_ , 2016.
  * [40] A. Chang, A. Dai, T. Funkhouser, M. Halber, M. Niessner, M. Savva, S. Song, A. Zeng, and Y. Zhang, “Matterport3D: Learning from RGB-D data in indoor environments,” _3DV_ , 2017.
  * [41] I. Armeni, S. Sax, A. R. Zamir, and S. Savarese, “Joint 2d-3d-semantic data for indoor scene understanding,” _arXiv preprint arXiv:1702.01105_ , 2017\.
  * [42] A. Dai, A. X. Chang, M. Savva, M. Halber, T. Funkhouser, and M. Nießner, “Scannet: Richly-annotated 3d reconstructions of indoor scenes,” in _CVPR_ , 2017.
  * [43] T. Schöps, T. Sattler, and M. Pollefeys, “BAD SLAM: Bundle adjusted direct RGB-D SLAM,” in _CVPR_ , 2019.
  * [44] I. Vasiljevic, N. Kolkin, S. Zhang, R. Luo, H. Wang, F. Z. Dai, A. F. Daniele, M. Mostajabi, S. Basart, M. R. Walter, and G. Shakhnarovich, “DIODE: A Dense Indoor and Outdoor DEpth Dataset,” _CoRR_ , vol. abs/1908.00463, 2019.
  * [45] T. Koch, L. Liebel, F. Fraundorfer, and M. Körner, “Evaluation of cnn-based single-image depth estimation methods,” in _ECCV Workshop_ , 2018\.
  * [46] H. Fan, H. Su, and L. J. Guibas, “A point set generation network for 3d object reconstruction from a single image,” in _CVPR_ , 2017.
  * [47] M. Tatarchenko*, S. R. Richter*, R. Ranftl, Z. Li, V. Koltun, and T. Brox, “What do single-view 3d reconstruction networks learn?” in _CVPR_ , 2019\.
  * [48] A. Knapitsch, J. Park, Q.-Y. Zhou, and V. Koltun, “Tanks and temples: Benchmarking large-scale scene reconstruction,” _ACM TOG_ , vol. 36, no. 4, 2017.
  * [49] A. Tversky, “Features of similarity,” _Psychological Review_ , vol. 84, no. 4, pp. 327–352, 1977.
  * [50] N. N. Johanna Wald, Armen Avetisyan, F. Tombari, and M. Niessner, “Rio: 3d object instance re-localization in changing indoor environments,” in _ICCV_ , 2019.
  * [51] K. Tateno, F. Tombari, I. Laina, and N. Navab, “Cnn-slam: Real-time dense monocular slam with learned depth prediction,” in _CVPR_ , 2017.
  * [52] H. Xie, H. Yao, S. Zhou, J. Mao, S. Zhang, and W. Sun, “Grnet: Gridding residual network for dense point cloud completion,” in _ECCV_ , 2020.
  * [53] M. Liu, L. Sheng, S. Yang, J. Shao, and S.-M. Hu, “Morphing and sampling network for dense point cloud completion,” _AAAI_ , 2020.
  * [54] Y. Wang, D. J. Tan, N. Navab, and F. Tombari, “Forknet: Multi-branch volumetric semantic completion from a single depth image,” in _ICCV_ , 2019\.
  * [55] Y. Yang, C. Feng, Y. Shen, and D. Tian, “Foldingnet: Point cloud auto-encoder via deep grid deformation,” in _CVPR_ , 2018.
  * [56] Y. Wang, D. J. Tan, N. Navab, and F. Tombari, “Softpoolnet: Shape descriptor for point cloud completion and classification,” _ECCV_ , 2020.
  * [57] S. Song, F. Yu, A. Zeng, A. X. Chang, M. Savva, and T. Funkhouser, “Semantic scene completion from a single depth image,” _CVPR_ , 2017.
  * [58] L. Pan, X. Chen, Z. Cai, J. Zhang, H. Zhao, S. Yi, and Z. Liu, “Variational relational point completion network,” in _CVPR_ , 2021.
  * [59] W. Yuan, T. Khot, D. Held, C. Mertz, and M. Hebert, “Pcn: Point completion network,” in _3DV_ , 2018.


060012 2014 G. Mindlin 060012

We study the effects of multiple binding sites in the promoter of a genetic
oscillator. We evaluate the regulatory function of a promoter with multiple
binding sites in the absence of cooperative binding, and consider different
hypotheses for how the number of bound repressors affects transcription rate.
Effective Hill exponents of the resulting regulatory functions reveal an
increase in the nonlinearity of the feedback with the number of binding sites.
We identify optimal configurations that maximize the nonlinearity of the
feedback. We use a generic model of a biochemical oscillator to show that this
increased nonlinearity is reflected in enhanced oscillations, with larger
amplitudes over wider oscillatory ranges. Although the study is motivated by
genetic oscillations in the zebrafish segmentation clock, our findings may
reveal a general principle for gene regulation.

# Nonlinearity arising from noncooperative transcription factor binding
enhances negative feedback and promotes genetic oscillations

Iván M. Lengyel  [inst1] Daniele Soroldoni  [inst2, inst3] Andrew C. Oates
[inst2, inst3] Luis G. Morelli[inst1] Email<EMAIL_ADDRESS>
(3 August 2014; 28 October 2014)

††volume: 6

99 inst1 Departamento de Física, FCEyN UBA and IFIBA, CONICET, Pabellón 1,
Ciudad Universitaria, 1428 Buenos Aires, Argentina. inst2 MRC-National
Institute for Medical Research, The Ridgeway, Mill Hill, NW7 1AA London, UK.
inst3 Department of Cell and Developmental Biology, University College London,
Gower Street, WC1E 6BT London, UK.

##  1 Introduction

Cells can generate temporal patterns of activity by means of genetic
oscillations [1, 2]. Genetic oscillations are biochemical oscillations in the
levels of gene products [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. They can be
produced by negative feedback regulation of gene expression, in which a gene
product inhibits its own production directly or indirectly [15]. Such
autoinhibition is often performed by transcriptional repressors, proteins or
protein complexes that bind the promoter of a gene and inhibit the
transcription of new mRNA molecules [16, 17], see Fig. 1. A theoretical
description of biochemical oscillations requires a delayed negative feedback,
together with sufficient nonlinearity and a balance of the timescales of the
different processes involved [18]. Delays may occur naturally in
transcriptional regulation, since the assembly of gene products involves
intermediate steps [16, 19, 20]. Nonlinearity refers to the presence of
nonlinear terms in the equations describing the dynamics. Such nonlinear terms
may occur in the equations due to the presence of cooperative biochemical
processes, where cooperativity is understood as a phenomenon in which several
components act together to orchestrate some collective behavior [21]. Although
some processes giving rise to nonlinear terms are known, in general it is
still an open question how nonlinearity is built into genetic oscillators.

Figure 1: Delayed autoinhibition can produce genetic oscillations. (A) The
gene (light blue box) is transcribed and translated into gene products x, with
a delay $\tau$. In this example, gene products form dimers that act as
transcriptional repressors, inhibiting transcription of the gene (blunted
arrow), and decay at a rate $c$. (B) and (C) Numerical solutions of Eq. (9)
describing the oscillator in A. (B) Phase space: monomer concentration $x(t)$
vs. the delayed concentration $x(t-\tau)$ settled to a limit cycle. (C) Dimer
concentration oscillates as a function of time. Parameters in B, C: $b\,P=2$,
$\tau=1$, $c=1$, $x_{0}=1$, $N=12$, $M=6$.

A compelling model system for genetic oscillations is the vertebrate
segmentation clock [22, 23, 24, 25]. This is a tissue-level pattern generator
that controls the formation of vertebrate segments during embryonic
development [25, 26]. The spatiotemporal patterns generated by the
segmentation clock are thought to be initiated by a genetic oscillator at the
single cell level [7, 27, 28, 29, 9, 10]. In this genetic oscillator, negative
feedback is provided by genes of the Her family, which code proteins that form
dimers and can act as transcriptional repressors [30, 31, 32, 33, 34]. The
time taken from transcription, translation and splicing may introduce the
necessary feedback delays in zebrafish and mouse [19, 35, 36, 34, 37]. One
source of nonlinearity in the segmentation clock oscillator is the
dimerization of gene products that bind the promoter of cyclic genes [32].
However, this may be insufficient to generate the observed oscillations in the
levels of gene products. One way to increase nonlinearity would be cooperative
binding of repressors to regulatory binding sites at the promoter [38, 21,
39]. Cooperative binding to multiple binding sites can make the negative
feedback steeper [40]. A similar effect occurs with ultrasensitivity in
phosphorylation cascades [41]. The presence of clusters of binding sites for
transcription factors may be a common motif in gene regulation [42], and
cooperativity in transcription factor binding has been reported for some
systems [43].

In zebrafish, multiple binding sites for Her dimers have been identified in
the promoter region of Her1, Her7, and other genes of the Notch pathway [44,
32, 45]. However, there is no evidence for cooperative binding of
transcriptional regulators in the case of the zebrafish segmentation clock.
Therefore, although cooperative binding is not ruled out, this lack of
evidence raises the general question of what contribution could be expected
from the multiple binding sites reported in the promoter of segmentation clock
cyclic genes. Here we use theory to study how multiple binding sites affect
nonlinearity and biochemical oscillations in a generic description of a
genetic oscillator.

##  2 A promoter with multiple binding sites

We first evaluate the regulatory function of a promoter that contains $N$
binding sites for a transcriptional repressor, see Fig. 2. We consider
transcriptional repressors although the results in this section are more
general and would apply to other types of transcription factors. We focus on a
single transcriptional repressor for the sake of simplicity, and assume that
all binding sites are identical. That is, binding and unbinding of repressors
to the different sites occur at the same rates for all sites. Moreover, we
assume that there is no cooperativity in repressor binding. This means that
binding and unbinding rates are not affected by the presence of bound factors
to any of the other sites.

Figure 2: Schematic representation of a promoter with multiple binding sites,
numbered from $1$ to $N$. Transcriptional repressors (orange squares) bind and
unbind from binding sites (numbered platforms) at the promoter of a gene
(light blue stretch). (A) and (B) illustrate two equivalent configurations of
the promoter with identical number of bound transcriptional repressors having
the same inhibiting strength.

The state ${\cal P}_{i}$ of the promoter at any time can be characterized by
the number $i$ of bound factors, which goes from $0$ to $N$. With a rate
$k^{+}$, a free repressor binds to an empty site at the promoter, which steps
from state ${\cal P}_{i}$ to state ${\cal P}_{i+1}$; with a rate $k^{-}$, a
bound repressor falls off the promoter, which steps from state ${\cal P}_{i}$
to state ${\cal P}_{i-1}$:

$\displaystyle{\cal P}_{0}\xrightleftharpoons[k^{-}]{k^{+}}$
$\displaystyle\cdots\xrightleftharpoons[]{}{\cal
P}_{i-1}\xrightleftharpoons[k^{-}]{k^{+}}{\cal
P}_{i}\xrightleftharpoons[k^{-}]{k^{+}}{\cal P}_{i+1}\xrightleftharpoons[]{}$
$\displaystyle\cdots\xrightleftharpoons[k^{-}]{k^{+}}{\cal P}_{N}.$ (1)

Denoting by $P_{i}$ the promoter occupation probabilities which describe the
fraction of time that the promoter spends at state ${\cal P}_{i}$ in the
thermodynamic limit [46], the kinetics of binding and unbinding of
transcriptional repressors to the binding sites is given by the set of
ordinary differential equations

$\displaystyle\dot{P}_{0}$ $\displaystyle=$
$\displaystyle-k^{+}NxP_{0}+k^{-}P_{1}$ $\displaystyle\vdots$
$\displaystyle\dot{P}_{i}$ $\displaystyle=$
$\displaystyle-k^{+}(N-i)xP_{i}-k^{-}iP_{i}$
$\displaystyle+k^{+}(N-i+1)xP_{i-1}+k^{-}(i+1)P_{i+1}$ $\displaystyle\vdots$
$\displaystyle\dot{P}_{N}$ $\displaystyle=$
$\displaystyle-k^{-}NP_{N}+k^{+}xP_{N-1}\,,$

together with the conservation law

$P=\sum_{i=0}^{N}{P_{i}}\,,$ (3)

where $x$ is the repressor concentration and $P$ is proportional to the number
of gene copies.

We assume here that binding and unbinding of repressors to the promoter occur
much faster than other processes like transcription, translation, transport
and decay of molecules. This means that the promoter occupation probabilities
quickly reach equilibrium with a given concentration of transcriptional
repressors [46, 47]. This situation is described by $\dot{P}_{i}=0$ for all
$i$ in Eq. (2), and the resulting algebraic equations can be solved by
induction to obtain

$P_{i}={N\choose i}\left(\frac{x}{x_{0}}\right)^{i}P_{0}\,,$ (4)

where $x_{0}=k^{-}/k^{+}$ is the equilibrium constant for binding of factors.
Using the constraint Eq. (3) we express $P_{0}$ in terms of $P$

$P_{0}=P\left(1+\frac{x}{x_{0}}\right)^{-N}\,.$ (5)

Equation (4) and Eq. (5) describe the equilibrium occupation of the promoter
in terms of the concentration $x$ of the transcriptional repressor.

##  3 Abrupt inhibition

In the previous section, we evaluated the kinetics of noncooperative binding
to a promoter with multiple binding sites. How does the presence of bound
transcriptional repressors affect the transcription rate of the gene
downstream of the promoter? In general, the strength of inhibition will depend
on the number of bound repressors, and the regulatory function $f(x)$ will
have the form

$f(x)=\sum_{i=0}^{N}{a_{i}P_{i}}\,,$ (6)

where $a_{0}=b$ is the basal transcription rate in the absence of bound
repressors and $a_{i}$ is the transcription rate in the presence of $i$ bound
repressors. Here we assume, for the sake of simplicity, that transcription
proceeds at its basal rate $b$ while there are $M$ or less sites occupied by
the repressors, and drops to zero when the number of occupied sites is larger
than $M$, see Fig. 3. We shall consider an alternative scenario below.

In this situation, Eq. (6) becomes

$f(x)=b\sum_{i=0}^{M}{P_{i}}\,.$ (7)

Using Eq. (4) and Eq. (5), the resulting regulatory function for a promoter
with $N$ binding sites is

$\displaystyle f_{N,M}(x)$
$\displaystyle=bP\left(1+\frac{x}{x_{0}}\right)^{-N}$
$\displaystyle\times\sum_{i=0}^{M}{{N\choose
i}\left(\frac{x}{x_{0}}\right)^{i}}\,.$ (8)

The zebrafish segmentation clock may be an interesting system to evaluate
these results. In the _her1/her7_ locus of zebrafish, an estimated number of
about $12$ binding sites has been reported [32, 45]. The regulatory functions
resulting from Eq. (3) for $N=12$ are displayed in Fig. 3B. Although it is
clear that inhibition of the promoter for a given level of repressors shifts
to the right as $M$ increases, it is less obvious how the value of $M$ affects
the steepness —that is the nonlinearity— of the negative feedback. We use Hill
functions to parametrize the regulatory functions Eq. (3) in a more
transparent form, see Appendix. Hill functions are parametrized by a Hill
coefficient $h$ characterizing the steepness of the curve, and an inhibition
threshold $K$ that describes the concentration of repressors that halves the
production rate. We fit Hill functions to $f_{N,M}$ and obtain an effective
Hill coefficient $h$ and effective inhibition threshold $K$, for each value of
$N$ and $M$, Fig. 4A,B. Increasing the number of binding sites $N$ while
keeping $M$ fixed can increase the Hill coefficient. For fixed $N$, increasing
$M$ changes the Hill coefficient in a nonmonotonic way: there is an optimal
value of $M$ that maximizes the Hill coefficient and therefore nonlinearity.
The effective inhibition threshold $K$ changes in a simpler form, increasing
both with $N$ and $M$. In conclusion, multiple binding sites can effectively
increase the nonlinearity of the feedback via the regulatory function.

Figure 3: Abrupt inhibition. Full inhibition occurs when more than $M$ sites
are bound by transcriptional repressors. (A) Transcription rate as a function
of the number of bound transcriptional repressors. (B) Normalized regulatory
function, Eq. (3), as a function of repressor concentration $x$, with $N=12$
and $M=0,\ldots,11$ (dark blue to dark red).  Figure 4: Multiple binding sites
can increase nonlinearity and enhance oscillations. (A) and (B) Effective Hill
parameters for the regulatory function $f_{N,M}$, Eq. (3), with $b\,P=1$ and
$x_{0}=1$. (A) Effective Hill coefficient $h$. (B) Effective inhibition
threshold $K$. (C) and (D) Oscillations described by Eq. (9) with $b\,P=2$,
$\tau=1$, $c=1$ and $x_{0}=1$. (C) Amplitude of oscillations. (D) Period of
oscillations. In C and D, the white region represents the nonoscillatory
regime, in which the system settles to a fixed point. Color bar labels
indicate values in each panel.

As discussed above, nonlinearity is an essential ingredient in a theory of
biochemical oscillations [18]. We therefore ask how multiple binding sites in
the promoter affect a biochemical oscillator. We use the regulatory function
Eq. (3) in a generic model for genetic oscillations. We consider a gene that
encodes a protein that forms dimers, and these dimers can bind to multiple
binding sites at the promoter to inhibit transcription, Fig. 1A. We introduce
an explicit delay $\tau$ to account for transcription, translation, splicing,
and other processes involved in the assembly of the gene product and its
dimerization. We assume that dimerization is a fast reaction, with a
separation of timescales from other processes. Therefore, at any time the
dimer concentration can be approximated by that of the monomers squared. The
dynamics of the product concentration $x(t)$ is given by

$\frac{dx}{dt}=bP\,\frac{\sum\limits^{M}_{i=0}{N\choose
i}\left(\frac{x(t-\tau)}{x_{0}}\right)^{2i}}{\left(1+\left(\frac{x(t-\tau)}{x_{0}}\right)^{2}\right)^{N}}-c\,x(t),$
(9)

where $b$ is the basal production rate, $c$ is the decay rate of products, $P$
relates to the number of gene copies, $\tau$ is the total delay for product
assembly, and $x_{0}$ is the product concentration that halves the basal
transcription rate. The regulatory function Eq. (3) is parametrized by $N$ and
$M$. This genetic oscillator is a reduction of the models proposed by Lewis
[19], and other authors [48, 49]. It describes the protein concentration, but
does not include the mRNA; the duration of transcription and translation are
both included in the delay $\tau$. Furthermore, it does not describe effects
present in theories that include more than one regulator [50, 32, 51], but
here it is enough to illustrate the effects of multiple binding sites in a
simpler context.

We integrate Eq. (9) numerically and evaluate the resulting dynamics by
calculating the amplitude and period of oscillations. In all numerical
simulations, we use the function dde23 from MATLAB [52]. Scanning the values
of $N$ and $M$, we determine whether the system oscillates in steady state:
when the difference in the maxima over the last ten cycles falls below 0.01,
the simulation is stopped and we record the output.

The amplitude of oscillations grows with the number of binding sites $N$, Fig.
4C. The range where the system oscillates grows with the number of binding
sites $N$. The amplitude is nonmonotonic in $M$: for a fixed number of binding
sites $N$, there is an optimal value for $M$ that maximizes the amplitude of
oscillations, Fig.4C. The period of oscillations grows with $N$ and decreases
with $M$. These results show how the change in nonlinearity observed in the
regulatory functions is reflected in the oscillations as the number of binding
sites change.

##  4 Gradual inhibition

There is strong evidence for the products of some cyclic genes of the
zebrafish segmentation clock binding their own promoters and acting as
transcriptional repressors [27, 35, 33, 44, 32, 23]. However, we do not have
detailed knowledge of how these transcriptional repressors affect
transcription rates when bound to the promoters of cyclic genes. There is
evidence from transcriptional analysis of the _Hes1_ gene in mouse indicating
that inhibition is gradual [30]. While the wildtype _Hes1_ promoter containing
all three N Box elements that are bound by Hes1 proteins showed a 30-fold
inhibition of transcription in the presence of Hes1, mutations in one, two,
and three of the N Box elements showed impaired inhibition with 14-, 7-, and
2-fold inhibition, respectively [30].

Motivated by these results, we consider here a scenario where additional bound
repressors gradually reduce transcription rate until it drops to zero, Fig.
5A. For the sake of simplicity, we assume that transcription rate drops
linearly as a function of bound repressors, from a basal rate $b$, to zero for
$k+1$ occupied sites

$a_{i}=\left\\{\begin{array}[]{ll}b\,\left(1-i/(k+1)\right)&\textrm{if }i\leq
k+1\\\ 0&\textrm{if }i>k+1\,,\end{array}\right.$ (10)

see Fig. 5A. Using this gradual inhibition in Eq. (6) together with Eq. (4)
and (5), we obtain regulatory functions $f_{N,k}(x)$

$\displaystyle f_{N,k}(x)$
$\displaystyle=bP\left(1+\frac{x}{x_{0}}\right)^{-N}$
$\displaystyle\times\sum_{i=0}^{k+1}{\left(1-\frac{i}{k+1}\right){N\choose
i}\left(\frac{x}{x_{0}}\right)^{i}}\,,$ (11)

see Fig. 5B. As in the previous case, it is clear that the effective
inhibition threshold shifts to the right as $k$ increases, but it is not so
clear if the steepness of the regulatory function changes and, if so, how.
Performing fits to Hill functions, we find that the nonmonotonic behavior of
the effective Hill exponent $h$ is observed again as $k$ increases, Fig. 6A,B.
Oscillations are similarly affected by noncooperative binding with gradual
inhibition, Fig. 6C,D. These results show that the prediction of an optimal
value for the number of bound repressors that fully inhibits the promoter is
robust with respect to the details of how multiple bound repressors reduce the
transcription rate.

Figure 5: Gradual inhibition. Binding of multiple transcriptional repressors
gradually inhibits transcription. (A) Transcription rate as a function of the
number of bound transcriptional repressors. Transcription occurs at the basal
rate $b$ in the absence of bound repressors, and decreases linearly to zero
for more than $k$ bound repressors. (B) Normalized regulatory function Eq. (4)
as a function of repressor concentration $x$, with $N=12$ and $k=0,\ldots,11$
(dark blue to dark red).  Figure 6: Multiple binding sites can increase
nonlinearity and enhance oscillations in the presence of gradual inhibition.
(A) and (B) Effective Hill parameters for the regulatory function $f_{N,k}$,
Eq. (4), with $b\,P=1$ and $x_{0}=1$. (A) Effective Hill coefficient $h$. (B)
Effective inhibition threshold $K$. (C) and (D) Oscillations obtained using
regulatory functions Eq. (4) with $b\,P=2$, $\tau=1$, $c=1$ and $x_{0}=1$. (C)
Amplitude of oscillations. (D) Period of oscillations. In C and D the white
region represents the nonoscillatory regime, in which the system settles to a
fixed point. Color bar labels indicate values in each panel.

##  5 Discussion

We studied the effects of multiple noncooperative binding sites in the
promoter of a genetic oscillator. We evaluated the behavior of a promoter with
multiple binding sites when binding of transcriptional repressors is
noncooperative, Fig. 2. We considered two different hypotheses for how bound
transcriptional repressors affect transcription rates, Figs. 3 and 5. In both
cases, we calculated how the number of binding sites and the number of bound
repressors required to produce full inhibition affect the nonlinearity of
regulatory functions. We showed that there is an optimal value of the number
of repressors needed to fully inhibit transcription that maximizes
nonlinearity, Figs. 4AB and 6AB. This increased nonlinearity is reflected in
the behavior of a genetic oscillator controlled by such regulatory functions,
Figs. 4CD and 6CD.

Cooperative binding is a well known means to increase the nonlinearity of a
biological dynamical system [21, 39]. Here we show that the nonlinearity of a
regulatory function can be increased by multiple binding sites, even if
binding is noncooperative. This idea may have application in other biological
control systems as well. Using the same formulation as we did here, one can
also describe nonlinearity in the transcriptional activation of gene
expression, thereby creating effective on-switches in developmental and
physiological regulatory networks. A similar effect has been reported in a
theoretical study of enzymes with multiple phosphorylation sites [53, 54, 55].
It was found that nonessential phosphorylation sites give rise to an increase
in effective Hill coefficients, enhancing ultrasensitivity in signal
transduction [55].

Previous work on the segmentation clock addressed the case of three binding
sites [40], motivated by experimental observations of the _Hes7_ mouse
promoter [31]. The authors assumed that a single bound dimer inhibits
transcription completely, corresponding to the particular case $M=1$ of the
theory developed here. They considered both noncooperative and cooperative
binding, and showed that cooperativity increases the effective Hill
coefficient as expected. In a follow-up [56], the authors addressed the case
of _Hes1_ regulation based on the report of four binding sites in the _Hes1_
mouse promoter [30]. Again assuming that a single bound dimer inhibits
transcription completely, they used the data from the transcriptional analysis
of the _Hes1_ gene to estimate an effective Hill coefficient for the Mouse
Hes1 oscillator, obtaining an upper bound of about $3$. More recently, the
effect of two binding sites as compared to a single binding site was discussed
together with differential decay of the monomers [57].

Our theory predicts how changing the number of binding sites $N$, and the
number $M$ of bound repressors that produce full inhibition affects a single
cell oscillator. Although an experiment that changes the value of $M$ may
currently be challenging in a cellular system, the number of binding sites $N$
is more amenable to experimental manipulation. For example, binding sites
could be mutated [30], or deleted from the promoter, or they could be
interfered with using genome editing strategies such as TALEN [58] or CRISPR
[59] to alter or delete specific binding sites. To assess the effects of these
perturbations in experiments may also pose some challenges. Dropping the
number of binding sites from $N=12$ to $N=6$ introduces a period change of
about 2.5%, while amplitude halves over the same range, Fig. 4C,D. Experiments
will require at least such precision to reliably detect changes.

Our results suggest a possible evolutionary mechanism to increase nonlinearity
in gene regulatory systems. In this mechanism, point mutations in the promoter
that increase the number of binding sites for transcription factors may
increase the steepness of regulatory functions. If the resulting steeper
regulation performs some function better, such mutations would have a good
chance to be conserved by natural selection. In the case of the segmentation
clock, the amplitude of oscillations could increase with the number of binding
sites, possibly reducing the signal to noise ratio. Furthermore, the range for
oscillations would be wider, making the oscillatory regime less sensitive to
slow extrinsic fluctuations of parameter values. Remarkably, it may happen
that after an increase in $N$, an increase in $M$ also raises nonlinearity,
see Fig. 4A. This means that weaker repressors would result in better
oscillations. This evolutionary mechanism would provide a simple way to
gradually increase the nonlinearity of a feedback or other regulatory
function.

The theory for the zebrafish regulatory function could be refined using the
experimentally measured relative affinities of the binding sites at the _her1_
and _her7_ promoters [32]. Apart from the effects reported here, the number of
binding sites may have additional roles. For example, it could serve as a
buffer for fluctuations in gene expression [60, 61, 62, 63, 64, 20],
augmenting the precision of genetic oscillations. This will be the topic of
future work.

###### Acknowledgements.

We thank Saúl Ares, Luciana Bruno, Ariel Chernomoretz and Hernan Grecco for
helpful comments on an early version of this work, and James Ferrell for
calling our attention to Ref. [53]. Thanks to E. Aikau for inspiration. LGM
acknowledges support from MINCyT Argentina (PICT 2012 1952). A.C.O. and D.S.
were supported by the Medical Research Council UK (MC_UP_1202/3) and the
Wellcome Trust (WT098025MA).

Figure 7: Hill functions are characterized by two parameters, the Hill
coefficient $h$ and the inhibition threshold $K$. (A) Hill functions with
$K=5$ and $h=1$ (red), $h=3$ (green) and $h=7$ (blue). (B) Hill functions with
$h=3$ and $K=3$ (red), $K=5$ (green) and $K=7$ (blue).

## Appendix: Effective Hill functions

Hill functions are often used to describe the nonlinearities present in gene
regulatory networks [17]. Hill functions are sigmoidal step functions defined
by

$f_{H}(x)=\frac{1}{1+\left({x/K}\right)^{2h}},$ (12)

where the steepnes of the step is characterized by the exponent $2h$, and the
inhibition threshold $K$ is the concentration of repressor that halves the
production rate, here scaled to unity, Fig. 7. Here we include an explicit
factor $2$ in the exponent to account for the dimerization of the
transcriptional repressors. One advantage of Hill functions is that they are
very simply parametrized. In the main text, we fit Hill functions to the more
complex regulatory functions Eqs. (3) and (4). Some fits of Eq. (3) in the
case $N=12$ are displayed in Fig. 8 for illustration.

Figure 8: Effective Hill functions Eq. (12) can fit regulatory functions of
the multiple binding site regulatory function, Eq. (3) for $N=12$, with
$b\,P=1$ and $x_{0}=1$. (A-C) Examples of fits of Hill functions (solid lines)
to regulatory functions (dashed lines) for (A) $M=1$ (green), (B) $M=6$ (light
blue) and (C) $M=11$ (purple). (D) Effective Hill coefficient $h$ as a
function of $M$. (E) Effective inhibition threshold $K$ as a function of $M$.
Dots in D and E correspond to fits from panels A, B and C.

## References

  * [1] A Goldbeter,  Biochemical Oscillations and Cellular Rhythms: The Molecular Bases of Periodic and Chaotic Behaviour,  Cambridge University Press, Cambridge (1997).
  * [2] L Glass, M C Mackey,  From Clocks to Chaos: The Rhythms of Life,  Princeton University Press, Princeton (1988).
  * [3] Y Liu, N F Tsinoremas, C H Johnson, N V Lebedeva, S S G M Ishiura, T Kondo,  Circadian orchestration of gene expression in cyanobacteria,  Gene. Dev. 9, 1469 (1995).
  * [4] E Nagoshi, C Saini, C Bauer, T Laroche, F Naef, U Schibler,  Circadian gene expression in individual fibroblasts: Cell-autonomous and self-sustained oscillators pass time to daughter cells,  Cell 119, 693 (2004).
  * [5] I Mihalcescu, W Hsing, S Leibler,  Resilient circadian oscillator revealed in individual cyanobacteria,  Nature 430, 81 (2004).
  * [6] A Goldbeter, C Gérard, D Gonze, J C Leloup, G Dupont,  Systems biology of cellular rhythms,  FEBS Lett. 586, 2955 (2012).
  * [7] I Palmeirim, D Henrique, D Ish-Horowicz, O Pourquié,  Avian hairy gene expression identifies a molecular clock linked to vertebrate segmentation and somitogenesis,  Cell 91, 639 (1997).
  * [8] A Aulehla, W Wiegraebe, V Baubet, M B Wahl, C Deng, M Taketo, M Lewandoski, O Pourquié,  A $\beta$-catenin gradient links the clock and wavefront systems in mouse embryo segmentation,  Nat. Cell Biol. 10, 186 (2008).
  * [9] Y Masamizu, T Ohtsuka, Y Takashima, H Nagahara, Y Takenaka, K Yoshikawa, H Okamura, R Kageyama,  Real-time imaging of the segmentation clock: Revelation of unstable oscillators in the individual presomitic mesoderm cells,  Proc. Natl. Acad. Sci. USA 103, 1313 (2006).
  * [10] A J Krol, D Roellig, M L Dequéante, O Tassy, E Glynn, G Hattem, A Mushegian, A C Oates, O Pourquié,  Evolutionary plasticity of segmentation clock networks,  Development 138, 2783 (2011).
  * [11] H Shimojo, T Ohtsuka, R Kageyama,  Oscillations in notch signaling regulate maintenance of neural progenitors,  Neuron 58, 52 (2008).
  * [12] N Geva-Zatorsky, N Rosenfeld, S Itzkovitz, R Milo, A Sigal, E Dekel, T Yarnitzky, Y Liron, P Polak, G Lahav, U Alon,  Oscillations and variability in the p53 system,  Mol. Syst. Biol. 2, 2006.0033 (2006).
  * [13] M B Elowitz, S Leibler,  A synthetic oscillatory network of transcriptional regulators,  Nature 403, 335 (2000).
  * [14] J Stricker, S Cookson, M R Bennett, W H Mather, L S Tsimring, J Hasty,  A fast, robust and tunable synthetic gene oscillator,  Nature 456, 516 (2008).
  * [15] J J Tyson,  Computational cell biology, Chap. 9, Pag. 230,  Springer, Berlin (2002).
  * [16] B Alberts, A Johnson, J Lewis, M Raff, K Roberts, P Walter,  Molecular biology of the cell, 4th Ed.,  Garland Science, New York (2002).
  * [17] U Alon,  An introduction to systems biology: Design principles of biological circuits,  Chapman & Hall/CRC Press, Boca Raton, Florida (2006).
  * [18] B Novák, J J Tyson,  Design principles of biochemical oscillators,  Nat. Rev. Mol. Cell Biol. 9, 981 (2008).
  * [19] J Lewis,  Autoinhibition with transcriptional delay: A simple mechanism for the zebrafish somitogenesis oscillator,  Curr. Biol. 13, 1398 (2003).
  * [20] L G Morelli, F Jülicher,  Precision of genetic oscillators and clocks,  Phys. Rev. Lett. 98, 228101 (2007).
  * [21] J E Ferrell,  Q&A: Cooperativity,  J. Biol. 8, 53 (2009).
  * [22] O Pourquié,  Vertebrate segmentation: From cyclic gene networks to scoliosis,  Cell 145, 650 (2011).
  * [23] A C Oates, L G Morelli, S Ares,  Patterning embryos with oscillations: Structure, function and dynamics of the vertebrate segmentation clock,  Development 139, 625 (2012).
  * [24] Y Saga,  The synchrony and cyclicity of developmental events,  Cold Spring Harb. Perspect. Biol. 4, a008201 (2012).
  * [25] D Roellig, L G Morelli, S Ares, F Jülicher, A C Oates,  Snapshot: The segmentation clock,  Cell 145, 800 (2011).
  * [26] Y Saga,  The mechanism of somite formation in mice,  Curr. Opin. Genet. Dev. 22, 331 (2012).
  * [27] A C Oates, R K Ho,  Hairy/E(spl)-related (Her) genes are central components of the segmentation oscillator and display redundancy with the Delta/Notch signaling pathway in the formation of anterior segmental boundaries in the zebrafish,  Development 129, 2929 (2002).
  * [28] H Hirata, S Yoshiura, T Ohtsuka, Y Bessho, T Harada, K Yoshikawa, R Kageyama,  Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop,  Science 298, 840 (2002).
  * [29] S A Holley, D Jülich, G J Rauch, R Geisler, C Nüsslein-Volhard,  her1 and the notch pathway function within the oscillator mechanism that regulates zebrafish somitogenesis,  Development 129, 1175 (2002).
  * [30] K Takebayashi, Y Sasai, Y Sakai, T Watanabe, S Nakanishi, R Kageyama,  Structure, chromosomal locus, and promoter analysis of the gene encoding the mouse helix-loop-helix factor hes-1. negative autoregulation through the multiple n box elements,  J. Biol. Chem. 269, 5150 (1994).
  * [31] Y Bessho, G Miyoshi, R Sakata, R Kageyama,  Hes7: A bhlh-type repressor gene regulated by notch and expressed in the presomitic mesoderm,  Genes Cells 6, 175 (2001).
  * [32] C Schröter, S Ares, L G Morelli, A Isakova, K Hens, D Soroldoni, M Gajewski, F Jülicher, S J Maerkl, B Deplancke, A C Oates,  Topology and dynamics of the zebrafish segmentation clock core circuit,  PLoS Biol. 10, e1001364 (2012).
  * [33] A Trofka, J Schwendinger-Schreck, T Brend, W Pontius, T Emonet, S A Holley,  The Her7 node modulates the network topology of the zebrafish segmentation clock via sequestration of the Hes6 hub,  Development 139, 940 (2012).
  * [34] A Hanisch, M V Holder, S Choorapoikayil, M Gajewski, E M Özbudak, J Lewis,  The elongation rate of RNA polymerase ii in zebrafish and its significance in the somite segmentation clock,  Development 140, 444 (2013).
  * [35] F Giudicelli, E M Özbudak, G J Wright, J Lewis,  Setting the tempo in development: An investigation of the zebrafish somite clock mechanism,  PLoS Biol. 5, 1309 (2007).
  * [36] E M Özbudak, J Lewis,  Notch signalling synchronizes the zebrafish segmentation clock but is not needed to create somite boundaries,  PLoS Genet. 4(2), e15 (2008).
  * [37] Y Harima, Y Takashima, Y Ueda, T Ohtsuka, R Kageyama,  Accelerating the tempo of the segmentation clock by reducing the number of introns in the hes7 gene,  Cell Rep. 3, 1 (2013).
  * [38] J Keener, J Sneyd,  Mathematical physiology I: Cellular physiology, 2nd Ed.  Springer, Berlin (2008).
  * [39] H Qian,  Cooperativity in cellular biochemical processes: Noise-enhanced sensitivity, fluctuating enzyme, bistability with nonlinear feedback, and other mechanisms for sigmoidal responses,  Ann. Rev. Biophys. 41, 179 (2012).
  * [40] S Zeiser, H V Liebscher, H Tiedemann, I Rubio-Aliaga, G K H Przemeck, M H de Angelis, G Winkler,  Number of active transcription factor binding sites is essential for the hes7 oscillator,  Theor. Biol. Med. Model. 3, 11 (2006).
  * [41] J Gunawardena,  Multisite protein phosphorylation makes a good threshold but can be a poor switch,  P. Natl. Acad. Sci. USA 102, 14617 (2005).
  * [42] V Gotea, A Visel, J M Westlund, M A Nobrega, L A Pennacchio, I Ovcharenko,  Homotypic clusters of transcription factor binding sites are a key component of human promoters and enhancers,  Genome Res. 20, 565 (2010).
  * [43] D S Burz, R Rivera-Pomar, H Jäckle, S D Hanes,  Cooperative dna-binding by bicoid provides a mechanism for threshold-dependent gene activation in the drosophila embryo,  EMBO J. 17, 5998 (1998).
  * [44] T Brend, S A Holley,  Expression of the oscillating gene her1 is directly regulated by hairy/enhancer of split, t-box, and suppressor of hairless proteins in the zebrafish segmentation clock,  Dev. Dynam. 238, 2745 (2009).
  * [45] D Soroldoni, personal communication (2014).
  * [46] L Bintu, N E Buchler, H G Garcia, U Gerland, T Hwa, J Kondev, R Phillips,  Transcriptional regulation by the numbers: Models,  Curr. Opin. Genet. Dev. 15, 116 (2005).
  * [47] H G Garcia, A Sanchez, T Kuhlman, J Kondev, R Phillips,  Transcription by the numbers redux: experiments and calculations that surprise,  Trends Cell Biol. 20, 723 (2010).
  * [48] N A M Monk,  Oscillatory expression of Hes1, p53, and NF-$\kappa$B driven by transcriptional time delays,  Curr. Biol. 13, 1409 (2003).
  * [49] M H Jensen, K Sneppen, G Tiana,  Sustained oscillations and time delays in gene expression of protein Hes1,  FEBS Lett. 541, 176 (2003).
  * [50] O Cinquin,  Repressor dimerization in the zebrafish somitogenesis clock,  PLoS Comp. Biol. 3, e32 (2007).
  * [51] A Ay, S Knierer, A Sperlea, J Holland, E M Özbudak,  Short-lived her proteins drive robust synchronized oscillations in the zebrafish segmentation clock,  Development 140, 3244 (2013).
  * [52] L F Shampine, S Thompson,  Solving ddes in Matlab,  Appl. Numer. Math. 37, 441 (2001).
  * [53] L Wang, Q Nie, G Enciso,  Nonessential sites improve phosphorylation switch,  Biophys. J. 99, L41 (2010).
  * [54] S Ryerson, G Enciso,  Ultrasensitivity in independent multisite systems,  J. Math. Biol. 69, 977 (2014).
  * [55] G Enciso,  Nonautonomous and random dynamical systems in life sciences Lecture Notes in Mathematics (Mathematical Biosciences Subseries) No 2102, Springer Verlag, Berlin (2013).
  * [56] S Zeiser, J Müller, V Liebscher,  Modeling the hes1 oscillator,  J. Comput. Biol. 14, 984 (2007).
  * [57] M Campanelli, T Gedeon,  Somitogenesis clock-wave initiation requires differential decay and multiple binding sites for clock protein,  PLoS Comp. Biol. 6, e1000728 (2010).
  * [58] V M Bedell, Y Wang, J M Campbell, T L Poshusta, C G Starker, R G K II, W Tan, S G Penheiter, A C Ma, A Y H Leung, S C Fahrenkrug, D F Carlson, D F Voytas, K J Clark, J J Essner, S C Ekker,  In vivo genome editing using a high-efficiency talen system,  Nature 491, 114 (2012).
  * [59] E Pennisi,  The crispr craze,  Science 341, 833 (2013).
  * [60] N Barkai, S Leibler,  Biological rhythms: Circadian clocks limited by noise,  Nature 403, 267 (2000).

  * [61] M B Elowitz, A J Levine, E D Siggia, P S Swain,  Stochastic gene expression in a single cell,  Science 297, 1183 (2002).
  * [62] J Raser, E O’Shea,  Noise in gene expression: Origins, consequences, and control,  Science 309, 2010 (2005).
  * [63] B Munsky, G Neuert, A V Oudenaarden,  Using gene expression noise to understand gene regulation,  Science 336, 183 (2012).
  * [64] L S Tsimring,  Noise in biology,  Rep. Prog. Phys. 77, 026601 (2014).