Guangzhi Sun, Xianrui Zheng, Chao Zhang, Philip C. Woodland††thanks: Guangzhi
Sun and Xianrui Zheng are supported by Cambridge Trust and an Amazon
Studentship respectively. This work has been performed using resources
provided by the Cambridge Tier-2 system operated by the University of
Cambridge Research Computing Service (www.hpc.cam.ac.uk) funded by EPSRC
Tier-2 capital grant EP/T022159/1.

# Can Contextual Biasing Remain Effective with Whisper and GPT-2?

###### Abstract

End-to-end automatic speech recognition (ASR) and large language models, such
as Whisper and GPT-2, have recently been scaled to use vast amounts of
training data. Despite the large amount of training data, infrequent content
words that occur in a particular task may still exhibit poor ASR performance,
with contextual biasing a possible remedy. This paper investigates the
effectiveness of neural contextual biasing for Whisper combined with GPT-2.
Specifically, this paper proposes integrating an adapted tree-constrained
pointer generator (TCPGen) component for Whisper and a dedicated training
scheme to dynamically adjust the final output without modifying any Whisper
model parameters. Experiments across three datasets show a considerable
reduction in errors on biasing words with a biasing list of 1000 words.
Contextual biasing was more effective when applied to domain-specific data and
can boost the performance of Whisper and GPT-2 without losing their
generality.

Index Terms: contextual biasing, TCPGen, Whisper, GPT-2, large language model,
end-to-end ASR

##  1 Introduction

Recent advancements in software and hardware have facilitated the significant
expansion of end-to-end automatic speech recognition (ASR) systems through the
use of copious amounts of training data. In this regard, large-scale
supervised models with extensive training, such as Whisper [1], have been
developed in addition to numerous foundation models featuring self-supervised
training, e.g. wav2vec2.0 and GPT-2 [1, 4]. Whisper contains a series of ASR
models that have been trained in a supervised manner on hundreds of thousands
of hours of labelled speech, which can be broadly applied across a diverse
range of ASR tasks. Unlike self-supervised foundation models, Whisper aimed to
provide a readily deployable ``out-of-the-box" solution for ASR tasks in
various environments, circumventing the need for supervised fine-tuning of a
decoder for every deployment. This approach aims to mitigate the issue of
decoder brittleness when fine-tuned for specific tasks.

Notwithstanding the comprehensive training data used in Whisper and GPT-2,
domain-specific words, such as proper nouns, are rare or absent in the general
training data and may be susceptible to errors when used for a specific
domain. On the other hand, fine-tuning these models on a limited amount of
domain-specific data may potentially compromise their ability to generalise to
other tasks. Therefore, it is essential to enhance the performance of these
models on such domain-specific words without sacrificing their capacity for
generalisation, and contextual biasing is one possible solution. Contextual
biasing aims to incorporate contextual knowledge into end-to-end ASR systems.
Contextual knowledge is often represented as a biasing list, comprising words
that carry important information and are essential to downstream tasks, such
as restaurant names in an ontology for a task-oriented dialogue system. The
inclusion of these words in the biasing list has been shown to significantly
improve recognition performance.

Based on the aforementioned concerns and the solution, this paper endeavours
to answer the question: ``Can contextual biasing remain effective in Whisper
and GPT-2?" with the constraint that the parameters of both models are
unchanged. Specifically, TCPGen is chosen as the neural biasing component as
it achieves biasing via distribution-level adaptation which does not require
any modifications to the Whisper model. This paper proposes modifications
TCPGen to be compatible with Whisper, together with a dedicated training
scheme on small-scale task-specific data. Notably, the proposed method enabled
the portability of TCPGen as a neural biasing component, allowing it to
maintain the "out-of-the-box" nature of Whisper without the necessity of fine-
tuning any of its components111Code is available at
https://github.com/BriansIDP/WhisperBiasing.git.

Experiments conducted across LibriSpeech, SLURP and Dialogue State Tracking
Challenge (DSTC2) datasets demonstrated that TCPGen consistently improved the
performance of biasing words with a biasing list of 1000 biasing words.
Contextual biasing was found to be more effective when applied to domain-
specific data such as SLURP and DSTC2, especially when using a more focused
biasing list.

The rest of this paper is organised as follows: Section 2 summarises related
work. Section 3 explains the use of TCPGen in Whisper and its training scheme.
Section 4 describes the experimental setup, followed by Section 5 analysing
results. The paper concludes in Section 6.

##  2 Related Work

###  2.1 Whisper and GPT-2 for ASR

Whisper [22] is a Transformer-based [19] encoder-decoder model trained purely
on supervised data for multiple languages and multiple speech-related tasks.
The training data contained paired multilingual speech and text obtained from
the Internet. The amount of training data consisted of 680k hours of audio, of
which 117k hours were non-English data. The released Whisper models have
various sizes, ranging from 39M parameters (4 encoder blocks/4 decoder blocks)
to 1550M parameters (32 encoder blocks/32 decoder blocks). The tokeniser used
the same set of Byte-Pair Encoding (BPE) tokens as the GPT-2 [4] model for
English-only models and with extra tokens for multilingual models. Whisper
provides competitive ASR results across multiple languages and domains,
without further fine-tuning.

GPT-2 [4] is a Transformer based language model trained to predict the next
wordpiece, without explicit supervision. The training data were gathered from
a social media platform, where the resulting text was more than ten billion
words. GPT-2 and its predecessor GPT [3] have been used to improve the
performance of an existing ASR model [23, 24].

###  2.2 Contextual biasing

Contextual biasing has typically used shallow fusion with a weighted finite-
state transducer [5, 6], via deep biasing [7, 8, 9, 10, 11], or via
deliberation models [17]. However, recently, combinations of shallow fusion
and deep biasing [12, 13] or neural shortcut-based methods [16, 14] have been
proposed. The work in [13] proposed and validated a simulation of contextual
biasing on open-source data, which was also adopted in subsequent work such as
[14]. Unlike deep biasing methods which required modifications to the ASR
model structure to take additional inputs, neural shortcut-based methods
directly modified the network output distribution without changing the
structure of the ASR model.

In particular, [16] and [14] both used the neural shortcut between the biasing
list and the final output. In contrast to [16] which used a Transformer
context encoder, TCPGen [14, 21] is a light-weight component which structured
biasing lists as wordpiece prefix trees to achieve efficient processing of
biasing lists containing thousands of words. Thus, TCPGen was selected as the
biasing component for Whisper to achieve efficient biasing without fine-tuning
any Whisper parameters.

##  3 Methodology

###  3.1 Tree-constrained pointer generator

TCPGen is a neural network-based component combining the symbolic prefix-tree
search with a neural pointer generator [18] for contextual biasing, which
enables end-to-end optimisation with ASR systems. At each output step, TCPGen
calculates a distribution over all valid wordpieces constrained by a word-
piece-level prefix tree built from the biasing list (referred to as the TCPGen
distribution). TCPGen also predicts a generation probability indicating how
much contextual biasing is needed at a specific step. The final output is the
interpolation between the TCPGen distribution and the original ASR model
output distribution, weighted by the generation probability.

Specifically, a set of valid wordpieces, $\mathcal{Y}^{\text{tree}}_{i}$, is
obtained by searching the prefix tree with a given history. Then, denoting
$\mathbf{x}_{1:T}$ and $y_{i}$ as input acoustic features and output
wordpieces respectively, $\mathbf{q}_{i}$ as the query vector carrying history
and acoustic information, and $\mathbf{K}=[...,\mathbf{k}_{j},...]$ as the key
vectors, scaled dot-product attention is performed to compute the TCPGen
distribution $P^{\text{ptr}}$ and an output vector
$\mathbf{h}^{\text{ptr}}_{i}$ as shown in Eqns. (1) and (2).

$P^{\text{ptr}}(y_{i}|y_{1:i-1},\mathbf{x}_{1:T})=\text{Softmax}(\text{Mask}(\mathbf{q}_{i}\mathbf{K}^{\text{T}}/\sqrt{d}))\vspace{-0.1cm}$
(1)
$\mathbf{h}^{\text{ptr}}_{i}=\sum\nolimits_{j}P^{\text{ptr}}(y_{i}=j|y_{1:i-1},\mathbf{x}_{1:T})\,\mathbf{v}^{\text{T}}_{j}\vspace{-0.1cm}$
(2)

where $d$ is the size of $\mathbf{q}_{i}$ (see [19]), Mask$(\cdot)$ sets the
probabilities of wordpieces that are not in $\mathcal{Y}^{\text{tree}}_{i}$ to
zero, and $\mathbf{v}_{j}$ is the value vector relevant to $j$. This paper
specifically focuses on the attention-based encoder-decoder (AED) ASR model.
In AED, the query combines the context vector and the previously decoded token
embedding, while the keys and values are computed from the decoder wordpiece
embedding, with a shared projection matrix. The generation probability which
takes a value between 0 and 1, is calculated using the decoder hidden state
and the TCPGen output vector $\mathbf{h}^{\text{ptr}}_{i}$. Then, the final
output can be calculated as shown in Eqn. (3).

$P(y_{i})=P^{\text{mdl}}(y_{i})(1-{P}^{\text{gen}}_{i})+P^{\text{ptr}}(y_{i})P^{\text{gen}}_{i}\vspace{-0.1cm}$
(3)

where conditions, $y_{1:i-1},\mathbf{x}_{1:T}$, are omitted for clarity.
$P^{\text{mdl}}(y_{i})$ represents the output distribution from the standard
end-to-end model, and $P^{\text{gen}}_{i}$ is the generation probability.

###  3.2 Biasing Whisper with TCPGen

Figure 1: Integration of TCPGen in Whisper with corresponding terms in Eqn.
((3)). The $\bigoplus$ symbol represents linear interpolation with weights
$P^{\text{gen}}$ predicted by TCPGen. Only TCPGen part is updated during
training. Biasing list example here contains Sophia, Franklin and Francisco,
and underscore marks the word start

The integration of TCPGen into the Whisper model is illustrated in Fig. 1.
TCPGen takes as input the Whisper final decoder block hidden state
$\mathbf{h}^{\text{dec}}$ to derive the query to calculate attention as shown
in Eqn. (4)

$\mathbf{q}_{i}=\text{ReLU}(W\mathbf{h}_{i}^{\text{dec}})$ (4)

where $W$ are the parameter matrix. Key and value vectors directly used the
decoder wordpiece embeddings. The computation of the generation probability
also used the final hidden state of the decoder, as shown in Eqn. (5).

$P^{\text{gen}}_{i}=\sigma(W_{1}\mathbf{h}_{i}^{\text{dec}}+W_{2}\mathbf{h}_{i}^{\text{ptr}})$
(5)

where $W_{1}$ and $W_{2}$ are the parameter matrices. Then, through Eqn. (3),
TCPGen leveraged this distribution-level adaptation to bias the Whisper
output, by only requiring the output from the final decoder block of Whisper.
During training on task-specific data, the entire Whisper model is frozen and
only parameters in TCPGen are updated. This way, Whisper can still generalise
to other data while TCPGen as a biasing component provides task-specific
biasing information. Moreover, as the word frequencies of the training data
for Whisper are not available, a word-error-based approach was adopted to
extract a good biasing list for training instead of a frequency-based one.
This was achieved by decoding the task-specific training set and gathering
distinct word tokens with errors higher than average.

###  3.3 Rescoring with GPT-2

GPT-2, as an off-the-shelf large pretrained LM, is used to rescore the
$N$-best hypotheses from the ASR model, as the Whisper large model has a
different set of wordpiece units to the GPT-2 model. To improve the
effectiveness of using GPT-2, an ILME method was also used following [25]. The
$N$-best hypotheses were derived by the beam search algorithm during Whisper
decoding, together with the log probability scores associated with each
hypothesis. In addition to GPT-2, the log probability of each hypothesis
subtracts the internal LM score, which is estimated by forwarding the Whisper
model with the encoder output set to all-zero vectors. The total score of each
hypothesis based on which the re-ranking is performed is computed by:

$\displaystyle\text{log}P_{\text{Whisper}}(Y|X)-\lambda_{\text{ilm}}\text{log}P_{\text{ilm}}(Y)+\lambda_{\text{GPT-2}}\text{log}P_{\text{GPT-2}}(Y)$
(6)

where $P_{\text{Whisper}}(Y|X)$ is the probability of hypothesis $Y$ given
acoustic feature $X$ predicted by Whisper, $P_{\text{ilm}}(Y)$ and
$P_{\text{GPT-2}}(Y)$ are the internal and external LM probabilities of $Y$,
and $\lambda$'s are hyper-parameters to be tuned.

Table 1: WER and R-WER on LibriSpeech test sets, SLURP test set and DSTC2 test set using Whisper and TCPGen, rescored with GPT-2. Test-time biasing list selection followed the description in Sec. 4. LM weights tuned on dev sets of each data separately. | Test-clean | Test-other | SLURP | DSTC2
---|---|---|---|---
System | WER $\downarrow$ | R-WER $\downarrow$ | WER $\downarrow$ | R-WER $\downarrow$ | WER $\downarrow$ | R-WER $\downarrow$ | WER $\downarrow$ | R-WER $\downarrow$
Whisper base.en | 5.2% | 16.3% | 10.5% | 31.4% | 27.8% | 72.5% | 21.1% | 71.3%
Whisper base.en + TCPGen | 4.8% | 12.5% | 10.1% | 26.5% | 25.8% | 52.3% | 18.2% | 41.8%
Whisper base.en + GPT-2 | 4.9% | 14.9% | 9.9% | 29.4% | 26.7% | 67.7% | 20.8% | 67.9%
Whisper base.en + TCPGen + GPT-2 | 4.5% | 11.7% | 9.7% | 25.0% | 24.8% | 47.0% | 16.9% | 38.7%
Whisper large | 4.0% | 10.4% | 6.7% | 20.0% | 16.7% | 61.6% | 17.3% | 69.1%
Whisper large + TCPGen | 3.4% | 8.3% | 6.3% | 16.3% | 15.1% | 37.1% | 14.7% | 33.5%
Whisper large + GPT-2 | 3.9% | 10.1% | 6.6% | 19.6% | 16.6% | 59.3% | 17.3% | 68.5%
Whisper large + TCPGen + GPT-2 |  3.4% |  8.2% |  6.3% |  16.3% |  15.0% |  37.1% |  13.9% |  29.5%

##  4 Experimental Setup

###  4.1 Data and test-time biasing lists

Experiments were performed on three datasets, including LirbiSpeech clean
data, SLURP data and the dialogue state tracking challenge Track 2 (DSTC2).
LirbiSpeech represents a fairly generic dataset while the other two are more
specific to certain tasks. Descriptions of each data set and the extraction of
the biasing list during inference are provided below.

LibriSpeech is an audiobook dataset. The train-clean-100 set was used for
training TCPGen and the two test sets were used for evaluation. Following the
simulation in [13, 14], a full rare word list containing 200k words was first
selected by removing the 5k most frequent words from the LibriSpeech LM
training vocabulary. During inference, one biasing list for each utterance was
obtained by collecting words that appeared in the full rare word list and
adding a certain number of distractors, which, if not specified, was 1000 for
LibriSpeech and SLURP.

SLURP [26] is a dataset containing single-turn user interactions with a home
assistant, annotated with scenarios, actions and entities, where 58 hours of
real speech data were used for training. Following [20], the full biasing list
was obtained by extracting words from slot entities that appeared in the
error-based biasing list for training. Out-of-training-set words in slot
entities were also included, making the full biasing list contain 2.8k words
in total. The biasing list for each utterance during inference was organised
in the same way as LibriSpeech.

DSTC2 is a human-machine task-oriented dialogue dataset where the user-side
speech input was used for training and evaluation, which corresponds to around
10 hours of speech for training, validation and testing respectively. The
biasing list for inference was organised by collecting ontology words that
appeared in the error-based biasing list for training, and, as SLURP, adding
out-of-training-set words. The full biasing list contains 380 words due to the
small ontology and hence was used as a whole without the need for reference
texts.

Note that as the Whisper tokeniser is case sensitive, for each biasing word, a
copy of it with the first character capitalised was also included in the
biasing list, which doubled the size of each biasing list. The feature
extraction, pre-processing and tokenisation followed the Whisper pipeline [1].

###  4.2 Model specifications and evaluation metrics

Experiments were performed with both the Whisper base.en model trained on
English data and the Whisper large trained on multilingual data222Code for
Whisper models: https://github.com/openai/whisper. The base model contained a
6-block Transformer encoder and decoder with 74M parameters, and the large
model contained a 32-block Transformer encoder and decoder with 1550M
parameters. TCPGen was trained separately on each data for 30 epochs using an
Adam optimiser with a linear tri-state learning rate scheduler. The training
was performed on a single A100 GPU. As an example, the training time of the
base.en model with TCPGen on LibriSpeech clean-100 was 5 hours. GPT-2 base
model was used to rescore 50-best lists without fine-tuning, with $\lambda$'s
searched between 0 and 1 on validations sets respectively.

In addition to WER, R-WER [13, 14, 20] was also used to evaluate system
performance. R-WER is the total number of error word tokens that belong to the
biasing list divided by the total number of word tokens in the test set that
belong to the biasing list. Insertion errors were counted in R-WER if the
inserted word belonged to the biasing list. Moreover, to keep the performance
comparable and compatible with other biasing setups [13, 14], unless
indicated, text normalisation was not used for scoring, and a discussion on
this is provided in Sec. 5.2.

##  5 Results

###  5.1 Overall Evaluation

The main results were summarised in Table 1. Compared to the Whisper base.en
model, using contextual biasing achieved an average 20% relative R-WER
reduction on test-clean and test-other sets, a 28% relative R-WER reduction on
SLURP and a 42% relative R-WER reduction on DSTC2. Reductions on domain-
specific data were much larger than that on the generic LibriSpeech data,
since the test time biasing list of LibriSpeech contained mainly generic
words, whereas that for the other two test sets contained domain-specific
words such as names of podcasts or restaurants. As biasing word tokens in
LibriSpeech and SLURP occupied around 10% of the total number of word tokens
in the test sets, the reduction in overall WER was greater than the error
reduction reflected in R-WER, indicating that contextual biasing was also
beneficial to unbiased words rather than degrading their performance. This
effect was most obvious for the DSTC2 data. Contextual biasing achieved larger
improvements with the large multilingual Whisper model, with 40% relative
R-WER reduction on SLURP and over 50% R-WER reduction on DSTC2.

When GPT-2 was applied to the Whisper base.en model, relative R-WER reductions
became slightly smaller for LibriSpeech as GPT-2 served as another generic
knowledge source, where the R-WER reduction was retained for other datasets.
However, applying GPT-2 to the large model resulted in a much smaller
improvement than to the base.en model. As before, larger relative WER and
R-WER reductions were observed with the large multilingual model when GPT-2
was used for rescoring. TCPGen achieved 18% relative R-WER reduction on
LibriSpeech test sets, 37% on SLURP and 57% on DSTC compared to the Whisper
large model with GPT-2 rescoring.

###  5.2 Discussion

First, the variation of model performance against the sizes of the test-time
biasing list was examined, where SLURP was used for this purpose as a non-
generic test set, as shown in Fig. 2. Up to the full-sized biasing list of
5.6k biasing words, which achieved 60.5% R-WER which did not require any
knowledge about the reference, TCPGen still outperformed the baseline by a
relative 18% reduction in R-WER.

Figure 2: WER and R-WER on SLURP using TCPGen and Whisper base.en model. The baseline is base.en without TCPGen. Full had 5.6k biasing words. Table 2: WER on LibriSpeech test-clean using Whisper large model and TCPGen, together with text normalisation. A beam size of 10 was used.  System | WER $\downarrow$ | R-WER $\downarrow$
---|---|---
Whisper norm. | 2.5% | 8.1%
Whisper + TCPGen norm. |  2.3% |  7.0%

Meanwhile, it is necessary to discuss the influence of text normalisation used
in the Whisper paper as a post-processing stage for both reference and
hypothesis. Results with normalisation were provided in Table 2 for
LibriSpeech test-clean set.

Text normalisation had a similar effect on test-other and SLURP, with a
relative R-WER reduction of 25% on SLURP. It was much less influential to
DSTC2 and the relative R-WER improvement remained at 50%. Nevertheless, text
normalisation may obscure the improvements in a biasing word 333For instance,
changes in possessive subtractions or British or American spellings that
correspond to words in the biasing list, it was necessary to evaluate the
performance of contextual biasing without text normalisation.

Using error-based biasing list: Since the training set for Whisper was
unknown, the discrepancy in the task-specific training set may cause degraded
performance with a frequency-based biasing list (as used in [13, 14])
extracted from it. In fact, for generic datasets such as LibriSpeech train-
clean-100, using a frequency-based biasing list by including words with a
training set frequency of less than 15, similar R-WER reductions were achieved
to the error-based one. However, the error-based biasing list for training
achieved a 12% greater reduction in R-WER compared to the frequency-based one
for domain-specific sets, SLURP and DSTC, as the distribution difference to
the Whisper training data is much larger.444Unknown training set also limited
word-frequency-based analysis

Finally, as TCPGen required training on a task-specific set, it is worthwhile
examining the performance of words that did not belong to that training set.
An OOV WER [21] was used which was measured the same way as R-WER but for
words not in the task-specific training set, with results shown in Table 3.
Whisper with TCPGen biasing achieved 30% relative OOV WER reduction compared
to Whisper alone, which was much higher than that achieved on other biasing
words. As a generic dataset, words not covered by the LibriSpeech training set
had a high chance of being rare in the Whisper training set and hence TCPGen
had a much larger improvement on them. The situation was reversed on SLURP and
DSTC as two domain-specific datasets, and words not covered by the training
set may be common in the training set of Whisper, especially for DSTC.
Therefore, the OOV WER reduction was smaller on those two datasets compared to
the R-WER reduction.

Table 3: Effect of TCPGen on words not in the task-specific training set using the Whisper large model, measured by OOV WER defined the same way as R-WER. OOV WER for test-clean and -other was calculated together. System | LibriSpeech $\downarrow$ | SLURP $\downarrow$ | DSTC2 $\downarrow$
---|---|---|---
Whisper | 35.0% | 67.1% | 15.1%
Whisper + TCPGen |  24.5% |  45.3% |  11.7%

##  6 Conclusions

This paper examined the effectiveness of neural contextual biasing for Whisper
in conjunction with GPT-2. Specifically, this paper proposes a method for
adapting the TCPGen component to Whisper, along with a dedicated training
scheme without the need to fine-tune Whisper. Experiments across three
datasets revealed that the proposed method resulted in a large reduction of
R-WER when a biasing list of 1000 words was utilized. In particular, 37% and
57% relative R-WER reductions were achieved on SLURP and DSTC2 respectively.
The findings also suggested that contextual biasing is more effective when
applied to domain-specific data and has the potential to enhance the
performance of universal models without sacrificing their capacity for
generalization. Though only Whisper and GPT-2 were investigated in this paper,
as future work, the biasing scheme can be easily integrated with other
universal ASR and LMs, such as Google's USM [27] if the parameters are open-
source.

## References

  * [1] A. Radford, J.W. Kim, T. Xu, G. Brockman, C. McLeavey & I. Sutskever,  ``Robust speech recognition via large-scale weak supervision",  Proc. NeurIPS, 2022.
  * [2] A. Baevski, Y. Zhou, A. Mohamed & M. Auli,  ``wav2vec 2.0: A framework for self-supervised learning of speech representations",  Proc. ICASSP, 2023.
  * [3] A. Radford, K. Narasimhan, T. Salimans and I.Sutskever,  Improving language understanding by generative pre-training,  OpenAI blog, 2018.
  * [4] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei & I. Sutskever,  ``Language Models are Unsupervised Multitask Learners",  OpenAI blog, 1(8), 9., 2019.
  * [5] I. Williams, A. Kannan, P. Aleksic, D. Rybach & T. Sainath,  ``Contextual speech recognition in end-to-end neural network systems using beam search",  Proc. Interspeech, Hyderabad, 2018
  * [6] D. Zhao, T. Sainath, D. Rybach, P. Rondon, D. Bhatia, B. Li & R. Pang,  ``Shallow-fusion end-to-end contextual biasing",  Proc. Interspeech, Graz, 2019.
  * [7] G. Pundak, T. Sainath, R. Prabhavalkar, A. Kannan & D. Zhao,  ``Deep context: End-to-end contextual speech recognition",  Proc. ICASSP, Calgary, 2018.
  * [8] Z. Chen, M. Jain, Y. Wang, M. L. Seltzer & C. Fuegen,  ``Joint grapheme and phoneme embeddings for contextual end-to-end ASR",  Proc. Interspeech, Graz, 2019.
  * [9] M. Jain, G. Keren, J. Mahadeokar, G. Zweig, F. Metze & Y. Saraf,  ``Contextual RNN-T for open domain ASR",  Proc. Interspeech, Shanghai, 2020.
  * [10] U. Alon, G. Pundak & T. Sainath,  ``Contextual speech recognition with difficult negative training examples",  Proc. ICASSP, Brighton, 2019.
  * [11] Z. Chen, M. Jain, Y. Wang, M. Seltzer & C. Fuegen,  ``End-to-end contextual speech recognition using class language models and a token passing decoder",  Proc. ICASSP, Brighton, 2019.
  * [12] D. Le, G. Keren, J. Chan, J. Mahadeokar, C. Fuegen & M. L. Seltzer,  ``Deep shallow fusion for RNN-T personalization",  Proc. SLT, 2021.
  * [13] D. Le, M. Jain, G. Keren, S. Kim, Y. Shi, J. Mahadeokar, J. Chan, Y. Shangguan, C. Fuegen, O. Kalinli, Y. Saraf & M. L. Seltzer,  ``Contextualized streaming end-to-end speech recognition with trie-based deep biasing and shallow fusion",  arXiv: Proc Interspeech, 2021. pp. 1772–1776.
  * [14] G. Sun, C. Zhang & P. C. Woodland,  ``Tree-constrained pointer generator for end-to-end contextual speech recognition",  Proc. ASRU, Cartagena, 2021.
  * [15] G. Sun, C. Zhang & P. C. Woodland,  ``Tree-constrained pointer generator with graph neural network encodings for contextual speech recognition",  Proc. Interspeech, Incheon, 2022.
  * [16] C. Huber, J. Hussain, S. Stüker & A. Waibel,  ``Instant one-shot word-learning for context-specific neural sequence-to-sequence speech recognition",  Proc. ASRU, Cartagena, 2021.
  * [17] X. Wang, Y. Liu, J. Li, V. Miljanic, S. Zhao & H. Khalil,  ``Towards Contextual Spelling Correction for Customization of End-to-End Speech Recognition Systems",  IEEE/ACM Transactions on Audio, Speech, and Language Processing, Vol. 30, pp. 3089–3097, 2022
  * [18] A. See, P. J. Liu & C. D. Manning,  ``Get to the point: summarization with pointer-generator networks",  Proc. ACL, Vancouver, 2017.
  * [19] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser & I. Polosukhin,  ``Attention is all you need",  Proc. NIPS, Long Beach, 2017.
  * [20] G. Sun, C. Zhang & P. C. Woodland,  ``End-to-end Spoken Language Understanding with Tree-constrained Pointer Generator",  Proc. ICASSP, 2023.
  * [21] G. Sun, C. Zhang & P. C. Woodland,  ``Tree-constrained pointer generator with graph neural network encodings for contextual speech recognition",  Proc. Interspeech, Incheon, 2022.
  * [22] A. Radford and J. Kim and T.Xu and G. Brockman C.McLeavey & I. Sutskever,  ``Robust Speech Recognition via Large-Scale Weak Supervision"  arXiv preprint arXiv:2212.04356, 2022
  * [23] K. Li, Z. Liu, T. He, H. Huang, F. Peng, D. Povey & S. Khudanpur,  ``An empirical study of transformer-based neural language model adaptation"  Proc. ICASSP, 2020.
  * [24] X. Zheng, C. Zhang & P. C. Woodland,  ``Adapting GPT, GPT-2 and BERT language models for speech recognition"  Proc. ASRU, 2021.
  * [25] Z. Meng, S. Parthasarathy, E. Sun, Y. Gaur, N. Kanda, L. Lu, X. Chen, R. Zhao, J. Li & Y. Gong,  ``Internal language model estimation for domain-adaptive end-to-end speech recognition"  Proc. SLT, 2021.
  * [26] E. Bastianelli, A. Vanzo, P. Swietojanski & V. Rieser,  ``SLURP: A spoken language understanding resource package",  Proc. EMNLP, 2020.
  * [27] Y. Zhang, W. Han, J. Qin, Y. Wang, A. Bapna, Z. Chen, N. Chen, B. Li, V. Axelrod, G. Wang, Z. Meng, K. Hu, A. Rosenberg, R. Prabhavalkar, D.S. Park, P. Haghani, J. Riesa, G. Perng, H. Soltau, T. Strohman, B. Ramabhadran, T. Sainath, P. Moreno, C.-C. Chiu, J. Schalkwyk, F. Beaufays, & Y. Wu.,  ``Google USM: Scaling Automatic Speech Recognition Beyond 100 Languages",  arXiv:2303.01037, 2023.


# Reply to: Assessing the precision of morphogen gradients
in neural tube development

Roman Vetter  Dagmar Iber

###### Abstract

This document supersedes the version published on 1 February 2024 in Nat.
Commun. 15, 930. The Matters Arising piece by Zagorski et al., that we were
invited to reply to, differed substantially from the version published in Nat.
Commun. 15, 929. The editorial process did not allow us to address these
changes in our reply published in Nat. Commun.

In a recent article [1], we demonstrated that single morphogen gradients in
the developing mouse neural tube (NT) can carry sufficient positional accuracy
to explain the patterning precision of progenitor domain boundaries. Zagorski
et al. had previously concluded otherwise [2], based on methodological
inconsistencies that we have revealed. The authors now comment on our work
with a Matters Arising letter. We rebut their criticism point by point in the
Supplement, and summarize the main aspects here.

The authors criticize our comparison of the different analysis methods
(FitEPM, NumEPM, DEEM) using a set of exponential functions. It is a
fundamental part of quantitative science to validate all methods used to
process data—be it part of experimental, theoretical or numerical work—against
known results. If a method fails to provide the correct results for known
problems, it is unreasonable to apply it to similar problems, for which the
answer is not known. The authors’ indirect approximation (FitEPM) leads to a
vast overestimation of the positional error in case of exponential gradients,
and is accurate only near the morphogen source [1]. The same limitations apply
to noisy gradient data, and the challenges of background subtraction and
smoothening apply to all three methods.

The authors insist that further away from the source the observed gradients
are flat, such that our arguments would not apply. First, their approximation
(FitEPM) will yield the wrong positional error whenever the mean deviates
substantially from an exponential [1] (i.e., also in flat parts of the
gradient), as it makes use of the exponential shape explicitly. Second, the
flat part of the gradients is not reliable biological data. We speculated in
our paper that the switch to a flat gradient shape may be due to insufficient
imaging depth, but the employed imaging depth remained unknown to us [1]. We
have since received confirmation that Zagorski et al. indeed employed 8-bit
imaging (A. Kicheva, personal communication, Supplementary Information). 8-bit
imaging only permits the detection of a $2^{8}=256$-fold intensity change. As
such, it is technically impossible to detect an exponential gradient beyond
$5.5$ times its decay length ($\approx 110$ µm) from the source, which
coincides with the point where Zagorski et al. find the transition from an
exponential to a flat gradient shape. The physical limits of their imaging
setup and the mathematical limitation of their approximation of the positional
error make it impossible to evaluate the positional error of the gradients at
later time points or further away from the source.

The authors nonetheless insist on the key conclusion from their paper, that
morphogen gradients are too noisy to specify cell fate in the NT beyond the
first 30 h. They argue that their analysis must be considered correct because
it 1) provides an explanation for the sensitivity of the progenitor markers to
both SHH and BMP, and 2) is consistent with their previous postulate that cell
differentiation rather than morphogen gradients defines the progenitor domain
boundary positions at later stages [3]. However, one can imagine other roles
than precision for this parallel SHH/BMP input [4], and their own work [3] and
that by others [5] showed that the key marker of the motor neuron domain,
OLIG2, remains sensitive to SHH signaling also at later stages. The limits
posed by 8-bit imaging and the inaccuracy of the chosen error approximation
are rigorous mathematical facts. They cannot be challenged biologically.

Given that the experimental data cannot be used beyond the 8-bit limit, we
developed computational approaches to estimate the gradient variability
further away from the source based on available experimental measurements.
This necessitated assumptions, which the authors now question. For one, based
on error propagation, we determined the expected gradient variability beyond
the 8-bit detection limit for the case that gradients maintain their measured
exponential shape across the domain. Secondly, we developed a cell-based
simulation framework that allows us to estimate the positional error from
measured molecular noise levels. Both approaches showed that gradients remain
sufficiently precise such that gradients can, in principle, pattern the NT
throughout its development.

Zagorski et al. claim that it is unrealistic to assume that gradients remain
exponential. However, they use exponential gradients, rather than their
measured gradients, in their own paper [2], both in their error propagation
method (FitEPM), and also as input when evaluating the potential of opposing
gradients in NT patterning via their decoding map, as the reported gradients
are too flat and noisy outside the 8-bit limit, also when both gradients are
considered simultaneously. As we showed in follow-up work, our results also
apply to non-exponential gradients [6], and gradient precision is
substantially higher still when taking into consideration that morphogens
spread at least in 2D, rather than 1D [7].

Figure 1: Gradient amplitude and decay length are uncorrelated. A Absence of
correlation between gradient amplitude and decay length in our cell-based
model with independent cell-to-cell molecular noise. Data from $n=1000$
simulated gradients obtained with independent noise between individual cells
in all three kinetic parameters, $\mathrm{CV}_{p,d,D}=0.3$. See [1] for method
details. B Fitting an exponential (orange) to an exponential gradient from
which a background was subtracted (blue, here one 256th of the maximum signal
$C_{\mathrm{max}}$ due to the 8-bit imaging depth used by Zagorski et al.)
leads to an error in the apparent gradient decay length that depends on the
gradient amplitude. For the shown example, $C_{0}/C_{\mathrm{max}}=0.2$. C We
demonstrate this with $n=1000$ synthetic noisy exponential gradients per value
of $C_{0}$, with equal $\lambda=19.26$ µm but different $C_{0}$
($C_{0}=\exp\left[(L/\text{\textmu{}m}-400)/118\right]$ as in Fig. 7D of [1]),
with independent Gaussian noise at a resolution of 1 µm. The noise level in
the gradients was set to $\mathrm{CV}_{C}=0.11$ as inferred from the gradients
with kinetic noise at $\mathrm{CV}_{p,d,D}=0.3$ from [1]. The apparent (i.e.,
fitted) decay length tends to be underestimated (orange, mean $\pm$ SD), and
correlates with the gradient amplitude relative to the maximum signal. D This
spurious correlation, as quantified by Pearson’s coefficient $R$ from the
above dataset within NT length windows of 75 µm (orange, mean $\pm$ SD),
disappears at later stages when the gradient amplitude increases. The
correlation that Zagorski et al. observe at early stages ($R=0.15$–$0.32$,
colored bars) lies well within the expected range of this fitting artifact.
The correlation coefficient they did not show for somite stages 25–30 is also
indicated ($R=-0.05$, black).

The authors further claim that our statistical approach was flawed because
gradient amplitude and decay length would be correlated (Pearson’s $R=0.26$,
their Fig. 2). Based on the same dataset [8], but including all developmental
stages, we concluded otherwise (Pearson’s $R=-0.0061$, Kendall’s $\tau=0.056$
[1]), and our cell-based model shows that such a correlation, while expected
in a deterministic setting, is negligibly small with independent cell-to-cell
molecular noise (Fig. 1A). The weak correlation that the authors now report
early in NT development arises as a technical artifact from the limitations of
8-bit gradient imaging (Fig. 1B–D). This highlights the pitfalls in gradient
parameter inference from imaging data and suggests that the reported gradient
variability is likely strongly overestimated, and that the true gradient
length is underestimated (Fig. 1C).

Figure 2: Effect of domain rescaling on the positional error estimate in the
tissue. Based on 50 synthetic gradients, Zagorski et al. claim that the
scaling error would be $0.15$ cell diameters (cd) at the beginning of NT
expansion. However, 50 gradients are not sufficient to obtain precise enough
statistics. A reanalysis using $n=\numprint{10000}$ gradients shows that it is
$0.32\pm 0.01$ cd (estimate $\pm$ SE) for the example picked by the authors.
Evaluating it as described by Zagorski et al., at gradient variability levels
as inferred in [1] (Supplementary Information) along the patterning axis
(red), shows that the degree of additivity is 61% (orange). For comparison,
the estimate used in [1] is shown (black line).

As we analyzed the processed gradient data that we received from the authors,
we noticed that the gradients had been binned from five somite stages and
scaled to the same domain length before determining the positional error. This
introduces an artificial positional error. We acknowledge that the simple
approximation that we used in [1] slightly overestimates the effect as the
errors are indeed only partially additive (Fig. 2). However, the difference is
not particularly important – even uncorrected, the difference between the
correctly inferred positional errors of the gradients and those reported for
the readouts is rather small and likely reflects, at least in part, remaining
technical limitations. We note that the authors overestimate the effect of
partial additivity (Fig. 2).

Zagorski et al. criticize that we did not also correct the readout data. The
method section of their Science paper [2] mentions the scaling neither for the
gradients nor for the readout, and the information remained inaccessible to us
before publication of our paper. We therefore measured the positional error of
the dorsal NKX6.1 boundary ourselves. We could reproduce their results only if
we bin, but do not scale the domains, from which we concluded that the
authors’ readout dataset was likely not scaled. During the writing of our
response to their Matters Arising letter, the authors informed us that they
did not scale the readout data (A. Kicheva, personal communication,
Supplementary Information). The raw data and scripts that would allow us to
check this remain inaccessible to us. While we remain interested in settling
the point, we consider domain scaling a minor issue, given the considerable
challenges in detecting gradients and aligning them with their readouts, in
particular in pseudostratified epithelia, where nuclei are not perfectly
aligned with their apical surface [9].

In conclusion, none of the points raised by Zagorski et al. bear relevance to
our conclusions, and remaining uncertainties could be clarified through access
to raw data and methodologies. See Supplementary Information for the detailed
rebuttal.

## References

  * Vetter and Iber [2022] R. Vetter and D. Iber.  Precision of morphogen gradients in neural tube development.  _Nat. Commun._ , 13:1145, 2022.  10.1038/s41467-022-28834-3.
  * Zagorski et al. [2017] M. Zagorski, Y. Tabata, N. Brandenberg, M. P. Lutolf, G. Tkac̆ik, T. Bollenbach, J. Briscoe, and A. Kicheva.  Decoding of position in the developing neural tube from antiparallel morphogen gradients.  _Science_ , 356:1379–1383, 2017.  10.1126/science.aam5887.
  * Kicheva et al. [2014] A. Kicheva, T. Bollenbach, A. Ribeiro, H. Pérez Valle, R. Lovell-Badge, V. Episkopou, and J. Briscoe.  Coordination of progenitor specification and growth in mouse and chick spinal cord.  _Science_ , 345, 2014.  10.1126/science.1254927.
  * Vetter and Iber [2023] R. Vetter and D. Iber.  Morphogen gradients can convey position and time in growing tissues.  _BioRxiv_ , 2023.  10.1101/2023.08.28.555133.
  * Yu et al. [2013] K. Yu, S. McGlynn, and M. P. Matise.  Floor plate-derived sonic hedgehog regulates glial and ependymal cell fates in the developing spinal cord.  _Development_ , 140:1594–1604, 2013.  10.1242/dev.090845.
  * Adelmann et al. [2023] J. A. Adelmann, R. Vetter, and D. Iber.  Patterning precision under non-linear morphogen decay and molecular noise.  _eLife_ , 12:e84757, 2023.  10.7554/eLife.84757.
  * Long et al. [2023] Y. Long, R. Vetter, and D. Iber.  2d effects enhance precision of gradient-based tissue patterning.  _iScience_ , 26:107880, 2023.  10.1016/j.isci.2023.107880.
  * Cohen et al. [2015] M. Cohen, A. Kicheva, A. Ribeiro, R. Blassberg, K. M. Page, C. P. Barnes, and J. Briscoe.  Ptch1 and Gli regulate Shh signalling dynamics via multiple mechanisms.  _Nat. Commun._ , 6:6709, 2015.  10.1038/ncomms7709.
  * Iber and Vetter [2022] D. Iber and R. Vetter.  Relationship between epithelial organization and morphogen interpretation.  _Curr. Opin. Genet. Dev._ , 75:101916, 2022.  10.1016/j.gde.2022.101916.

### Author Contributions

R.V. and D.I. jointly reviewed the Matters Arising letter and wrote the reply.

### Competing Interests

The authors declare that they have no competing interests.

### Code Availability

The source code for Figs. 1 and 2 is publicly released under the 3-clause BSD
license. It is available as a git repository at
https://git.bsse.ethz.ch/iber/Publications/2024_vetter_gradient_variability_ma.

Supplementary Information

In this supplementary document, we respond to the criticism by Zagorski et al.
on our article [1] point by point. The Matters Arising manuscript by the
authors is reprinted in bold face, and our response follows each point in
Roman.

Zagorski et al.: In the developing neural tube, pattern forms in response to
opposing BMP and Shh signaling gradients1. In a recent publication, Vetter and
Iber present theoretical analysis based on which they conclude that a single
morphogen gradient in the neural tube is sufficient to precisely position gene
expression boundaries2. Here we discuss assumptions made by Vetter and Iber
that limit the conclusions they reach, and address inaccuracies in their
analysis. Given these limitations and existing evidence, it seems likely that
both signaling gradients contribute to the precision of pattern formation in
the neural tube.

Reply: We disagree with the authors that there are inaccuracies in our
analysis, as detailed below. On the contrary, our work highlights several
inaccuracies in the authors’ 2017 Science paper [2].

Zagorski et al.: In multiple systems, morphogen gradients have been studied by
measuring fluorescent reporters of signaling activity in fixed tissues3. A
common practice is to estimate the imprecision of a gradient by assessing the
variation in fluorescent intensity (FI) between individual embryos at every
position in the tissue4. The positional error $\sigma_{x}$ of the gradient is
approximated by multiplying the variation of morphogen levels $\sigma_{C}$ by
the local gradient steepness $\absolutevalue{\frac{\partial C}{\partial
x}}^{-1}$ at that position: $\sigma_{x}\approx\absolutevalue{\frac{\partial
C}{\partial x}}^{-1}\sigma_{C}$. Vetter and Iber point out that different
methods for estimating the local gradient steepness can produce different
results. One method, numEPM, uses the spatial derivative of mean intensity at
the position of interest. Another method, fitEPM, assumes that the mean
gradient is exponential. In this case, the local steepness of the gradient is
given by the fitted mean intensity at a position divided by fitted exponential
decay length. A third method, DEEM, estimates the positional error as the
standard deviation of positions $x_{\theta,i}$ that correspond to a defined
concentration threshold: $\sigma_{x}=\operatorname{SD}\\{x_{\theta,i}\\}$. The
DEEM method is derived from the mathematical definition of positional error
and hence considered to represent the most direct measure of positional error
from an ensemble of gradients.

For low FI values, numEPM and fitEPM methods are influenced by how background
FI is estimated and subtracted and by how data is binned and smoothed along
the axis. Thus, in the tail of a gradient, the positional error estimates
generated by the two methods are inexact and may differ.

Reply: Firstly, challenges of background subtraction and smoothening apply
also to DEEM. As the authors confirmed (A. Kicheva, personal communication,
Supplementary Information), they employed 8-bit imaging, so that the flat tail
of their gradients reflects the detection limit of their chosen imaging depth.
Under such circumstances, any computational analysis is meaningless beyond the
detection limit.

Secondly, the differences between the computational methods exist
independently of any measurement consideration. FitEPM yields the wrong
positional error whenever the mean of the gradients deviates from an
exponential function. The gradients that have been reported by Zagorski et al.
are fitted well by exponential functions close to the source, and are roughly
constant in the center of the domain. In both cases, FitEPM yields the wrong
result further away from the source as neither the mean of exponential
functions nor the mean of constant functions is an exponential function.

Zagorski et al.: Vetter and Iber claim they can determine which of the two
methods is correct by testing which method gives the result closest to
estimating the precision of an artificial dataset consisting of an ensemble of
exponential gradients using the DEEM method. This leads them to conclude that
NumEPM is correct while FitEPM overestimates the positional error. However,
this conclusion depends on the assumption that experimental gradients are
perfectly exponential.

Reply: It is an indispensable pillar of quantitative science to test the
accuracy of methods by applying them to problems for which the result is
known. If a method fails to provide the correct results for known problems, it
is unreasonable to apply it to similar problems, for which the correct answer
is not known. In this spirit, we are comparing the three methods (FitEPM,
NumEPM, DEEM) using a set of exponential functions, as for those, the result
is known. We are not claiming that morphogen gradients are necessarily
perfectly exponential. What we are showing is that for exponential functions,
FitEPM fails to yield accurate results (except close to the source), and thus
it cannot be expected to work for similarly shaped noisy gradients.

To compare the positional errors of morphogen gradients and of their readout,
the same method must be used to calculate both of them. The authors chose DEEM
to quantify the positional error of the readout. DEEM is unrelated to the
shape of the gradients or their readouts. Unlike DEEM, the two other methods,
NumEPM and FitEPM, estimate the standard deviation of the position indirectly,
via error propagation. We can therefore evaluate whether NumEPM and FitEPM are
appropriate by comparing their results to DEEM, when applied to data. As we
show, NumEPM gives very similar results as DEEM, while FitEPM greatly
overestimates the positional error for sets of exponential gradients that have
the same variability in their two parameters (amplitude and gradient length)
as reported by Zagorski et al., based on their own fitting of their gradient
data with exponential gradients. In fact, the positional error that we obtain
with FitEPM for the perfectly exponential gradients is very similar to the
results reported by Zagorski et al. based on the measured, noisy data (Fig.
1D).

Further away from the source, the reported gradient shapes lose their
exponential shape and become flat. Here, FitEPM performs even worse, as it is
explicitly based on the assumption of a continued exponential shape.

Zagorski et al.: The cellular response to the signal and tissue
heterogeneities result in gradient shapes that deviate from a perfectly
exponential curve5,6.

Reply: We agree with this statement, and for non-exponential gradients, the
mean gradient will potentially deviate even further from an exponential
function, introducing an even greater error when using FitEPM.

Zagorski et al.: The poor signal-to-noise ratio in the gradient tail means
that the real shape of gradients in this region cannot be reliably measured.

Reply: It is not just that the gradient cannot be reliably measured, but
rather not be detected at all, except close to the source. Given the author’s
use of 8-bit imaging, which limits the detectable signal range to 256-fold,
the reported flat gradient “signal” is likely just technical noise.

The fact that the measured gradients turn from an exponential shape to a flat
noisy shape near the 8-bit limit led us to suspect that the reported poor
signal-to-noise ratio is most likely the consequence of an inappropriate
imaging regime. The imaging depth is not provided in the methods section of
[2], and was inaccessible to us prior to the publication of our article. The
authors meanwhile confirmed that they had used 8-bit imaging (A. Kicheva,
personal communication, Supplementary Information). It is unclear to us why
the authors nonetheless insist that after 30 hours, morphogen gradients are
too imprecise to inform patterning in the center of the neural tube — even
though the employed 8-bit imaging depth does not allow to image the gradient.

Zagorski et al.: Thus, judging the two methods by comparison to an artificial
dataset, which may not represent the true shape of gradients, is misleading.
In other words, the performance of a method on an idealized dataset does not
determine whether this method will work well on real data which may differ
from the idealized dataset.

Reply: As we explained above, the same mathematical analysis method needs to
be used when evaluating the positional error of gradient and readout. This was
not done by Zagorski et al., but a method was used that even with perfectly
exponential gradients wrongly yields huge positional errors unless the
gradients are identical. Our analysis shows that when the same methods are
used for both the gradient and readout, their positional errors are compatible
with each other.

Zagorski et al.: More importantly, the analysis of Vetter and Iber indicate
that there is in fact very good agreement between the precision estimated by
the different methods during the relevant stages of neural tube development
(0–15 ss, corresponding to 0–30 h). An examination of their Fig 1E shows that
the two methods produce identical precision estimates for time points 0–5 ss.
For 10–15 ss, the estimates are also very similar and diverge only in the
gradient tail: DEEM and numEPM estimate 5–6 cell diameters, fitEPM results in
6–8 cell diameters. These positional errors occur at distances $>60\%$ tissue
length from the morphogen source for GBS-GFP and $>45\%$ for pSmad.

Reply: As we pointed out in our article, it lies in the nature of the method
that the mean of exponential gradients is fitted well by an exponential
function close to the source, but not further away. Our paper does not claim
that FitEPM is inaccurate early on or very near the source, but later on and
away from the source.

Zagorski et al. were not able to image the exponential gradients beyond somite
stage 15 because they employed 8-bit imaging, which allowed them to detect
fluorescence intensity changes only over a 256-fold range. Beyond the 8-bit
limit ($\approx 110$ µm), their imaging returns technical noise, likely
explaining the sharp transition from an exponential to a flat gradient shape.
Given this technical detection limit, it is impossible to know whether the
reported flat part of the gradient reflects biological reality or the
technical limits of 8-bit imaging (although the latter seems more likely to
us). No matter which of the two scenarios applies, FitEPM will return a wrong
positional error in the flat part of the reported gradients. FitEPM relies on
an exponential fit to the mean signal, and this fit is dominated by the higher
signal close to the source. As such, FitEPM neither provides a correct
estimate of the positional error if the gradients were indeed flat, nor for
the case that the exponential gradients continued with the same gradient shape
as near the source (but cannot be detected by 8-bit imaging). We analyzed the
latter scenario via error propagation methods in our article and showed that
the expected positional error of the gradients would be largely consistent
with that of the readouts if the gradients indeed remained exponential.

In summary, by employing 8-bit imaging and using FitEPM over the entire
spatial and temporal range for the gradients, but DEEM for the readouts, the
2017 Science paper wrongly concludes that single gradients become too
imprecise to pattern the center of the neural tube.

Zagorski et al.: This similarity in the estimates at early stages is relevant,
because, as we show1, early (before 15 ss) but not late stage gradients are
used to establish pattern.

Reply: The authors provide no such evidence in their paper. In fact, in a
previous paper [3], Kicheva and Briscoe show that late interference with SHH
signaling still impacts OLIG2 levels, in agreement with previous work by
others [4].

Zagorski et al.: In 1, we derive a decoding map of Shh and BMP signaling using
the profiles measured at 5 ss (Fig. 2A and S3A therein). We validate this map
with experiments that are independent of how the morphogen signaling gradients
were imaged.

Reply: The part of our paper that addresses the authors’ 2017 Science paper is
only concerned with the correct method to compare the positional errors of
gradients and their readouts, not with the decoding map. But as the authors
bring this up, we emphasize that the decoding map would not work to explain
the position of domain boundaries were the authors to use their measured
gradients at time points other than the very first. At these later time
points, at least one of the two gradient profiles reported by Zagorski et al.
is flat in the part of the domain where the NKX6.1 and PAX3 boundaries are
observed. Consequently, their reported gradients could not define the readout
position, whether they use a decoding map based on opposing gradients or a
direct threshold-based readout of a single gradient. The presented experiments
are also consistent with alternative explanations and provide no proof for the
authors’ decoding map.

Zagorski et al.: We demonstrate that the downstream transcriptional network
requires morphogen input for $<30$ h to generate the pattern and this
mechanism is sufficient to maintain gene expression in the absence of ongoing
signaling at late stages. This reinforces previous experimental evidence,
based on growth rate measurements, lineage tracing and perturbation
experiments, that indicates that the temporal window for morphogen-dependent
cell fate specification is during the first 30 h of mouse neural tube
development7. Thus, for the time interval that is relevant for pattern
formation, the fitEPM, NumEPM and DEEM methods produce similar estimates of
positional error.

Reply: Despite all their quantitative work, the authors have never shown that
this transcriptional network yields the correct progenitor domain boundaries
when the measured gradients are used as input. Even though Zagorski et al.
claimed above that our use of exponential gradients was unrealistic, they used
an idealized exponential gradient, rather than their measured gradients, as
input, and did not quantitatively match their boundary data over time.

As to the claim that morphogen gradients specify patterns only during the
first 30 h of NT development, genetic experiments have shown that removal of
the Shh source at later stages affects SHH-dependent NT patterning [4], which
should not happen if patterning became independent of morphogen input after 30
h. The claim that the gradients were too imprecise beyond the first 30 h
hinges on the use of a mathematical flawed method, and on 8-bit imaging that
makes it impossible to detect exponential gradients beyond this point, as
detailed above.

The conclusion that the rate of differentiation changes rapidly at around 30 h
in the cited publication [3] hinges on a mathematical analysis that is
inaccessible to us and that cannot be reproduced by us with their published
data.

But again, this is of no concern regarding the topic of our analysis, which
focuses on a consistent way to compare the positional errors of gradients and
their readouts.

Zagorski et al.: Vetter and Iber also argue that imprecision of the signaling
gradients is overestimated by grouping signaling profiles into temporal bins
that correspond to 10h of developmental time. For a given bin, all signaling
profiles are assumed to have the same DV length. Vetter and Iber suggest that
this introduces a “scaling error”. To define it, they assume that all profiles
in the temporal bin have equal amplitudes and decay lengths, but different
absolute lengths. They reason that any resulting positional error is therefore
the product of the differing lengths, rather than actual variability in the
amplitude and decay length.

Reply: This is a misinterpretation of our work. In the quantitative part (Fig.
2C–E,H–J), we are using the variability in the amplitude and gradient length
that the authors reported in their 2017 Science paper. Only to illustrate the
concept of the scaling error in the schematics of Fig. 2F,G, the same $C_{0}$
and $\lambda$ are used for visual clarity. This variability was determined
after the gradients had been scaled. We are correcting for this effect by
exploiting the reported uniform NT growth behavior. But we agree that this
correction can, unfortunately, only be approximate.

Zagorski et al.: This reasoning is problematic. First, if the signaling
gradient profiles are corrected in this way, so should the gene expression
boundaries of Pax3 and Nkx6.1. Vetter and Iber did not do this. Instead they
compare the corrected signaling gradients to the imprecision of Pax3 and
Nkx6.1 as reported in Zagorski et al., that is without correction.

Reply: The information whether Zagorski et al. also scaled the NT length in
each 5-somite bin before evaluating the PAX3 and NKX6.1 boundaries, and
whether they would therefore need to be corrected as well, was inaccessible to
us before publication of our paper. When we determined the positional error of
the PAX3 and NKX6.1 domain boundaries ourselves using our own imaging data, we
could reproduce the reported positional error only when we binned, but did not
scale the domains. Having said this, this could, of course, be the result of
differences between mouse lines. Based on this finding and the different
methods that Zagorski et al. used to calculate the positional error for
gradients and readouts, we concluded that it is most plausible that they
measured the PAX3 and NKX6.1 boundaries on unscaled domains.

During the writing of our response to their Matters Arising letter, the
authors informed us that they did not scale the readout data (A. Kicheva,
personal communication, Supplementary Information). The raw data and scripts
that would allow us to check this remain inaccessible to us.

The remaining difference is in any case small (at most one cell diameter), and
we expect that additional effects from the epithelial structure
(pseudostratification) amplify the measured positional error [5], and that,
when carefully remeasured, the positional errors of both gradients and
readouts will turn out to still be lower. As we show in follow-up work, the
expected positional error of the morphogen gradients is even lower when
considering that they are not spreading on a 1D domain, but at least in 2D
[6].

Zagorski et al.: Furthermore, by subtracting the scaling error, Vetter and
Iber assume that it has an additive contribution to the overall profile
variability. This excludes the possibility that variability in decay length
and amplitude could dominate any scaling variability. In such a scenario,
subtracting the scaling error would lead to unrealistic underestimation of the
actual error (Fig. 1). Taken together, the proposed “scaling error” correction
is applied inconsistently and might underestimate the actual variability.

Reply: Given our lack of access to the authors’ raw imaging data, the degree
of additivity of the scaling error in the authors’ methodology on the
biological positional error remains unknown to us. However, as we show in our
article, the positional error of the gradients and the readout are comparable
even if the errors are only partially additive.

Zagorski et al.: Vetter and Iber suggest that gene expression boundaries in
the neural tube are positioned by a single morphogen gradient, rather than the
combined interpretation of both signaling pathways.

Reply: Our work shows that single gradients can provide sufficient positional
information in the neural tube. We do not claim that other precision-enhancing
effects are excluded, nor that gene expression boundaries are positioned by
just one of the morphogen gradients in the neural tube.

Zagorski et al.: Implicit in this idea is that cells somehow distinguish which
of two independent gradients is the most precise and use that to determine
their identity.

Reply: We do not make such an assumption. It appears to us that the authors
might have misinterpreted our Fig. 6A,D,E to arrive at this impression. In
Fig. 6B,D,E we explicitly show that such an assumption is not necessary.

Zagorski et al.: This interpretation also misses a crucial point: there is
experimental evidence that neural progenitors respond to combinations of
signaling factors. Consistent with prior studies8, we1 show that neural
progenitor identities depend on the levels of both BMP and Shh signaling.

Reply: We do not challenge these experimental observations. However, there is
no evidence that they serve to increase precision. In fact, as shown by others
[7], also the ventral-most NKX2.2 domain responds to BMP, even though
according to the analysis by Zagorski et al., the BMP gradient would be way
too noisy to yield any useful positional information there.

Zagorski et al.: Vetter and Iber further suggest that gradient variability can
be accurately inferred from “summary statistics of exponential gradients”.
This necessitates several assumptions. First, gradients are assumed to be
exponential. However, diffusion and degradation often depend on feedback from
morphogen signaling, which can lead to deviations from exponential shape9.

Reply: In this second part of our paper, we developed an error propagation
formula to estimate the positional error far from the source (where
measurements have so far reached technical limits) based on the variability
that can be measured near the source. Whether or not gradients are exponential
is indeed an open question. For those cases, in which gradients do remain
exponential, our formula can be used. Recognising that they may not
necessarily need to do so in other systems of gradient-based patterning, we
also presented a simulation framework in the final part of our paper that can
be used to estimate the positional error also for more complex cases, as the
authors refer to. Using this simulation framework, we have recently shown that
non-linear decay, which results in non-exponential gradients, and which had
previously been proposed to increase the robustness to variability in the
source [8], yields very similar results as linear decay, and does not result
in a relevant increase in precision [9]. Morphogen spreading in 2D rather than
in 1D, however, substantially reduces gradient variability and therefore
permits subcellular precision for physiological levels of molecular noise [6].
Together, this reinforces our previous findings that single morphogen
gradients can be sufficiently precise to pattern the mouse neural tube, also
when considering additional physiological aspects.

Zagorski et al.: Second, ligand and signaling gradients are assumed to have
comparable variability and any discrepancy results from technical measurement
errors. This ignores the possibility that the signal transduction mechanisms
alter the noise properties of a signal10.

Reply: We do not make such an assumption. However, if the readout was more
precise than the gradients, then the missing information would have to be
introduced somehow. Zagorski et al. proposed that this missing information is
obtained by reading out both gradients simultaneously in the center of the
domain. As we show, single gradients can be sufficiently precise to define the
progenitor domain boundaries also in the center of the NT, offering a simple
explanation of how patterning is controlled in the NT.

Zagorski et al. themselves make the assumption that the GBS-GFP and pSMAD
gradients can serve as proxies for the SHH and BMP gradients, even though GBS-
GFP may respond very differently from more dorsal SHH-dependent genes, as we
explained in detail in our paper.

Zagorski et al.: Third, variables, such as $C_{0}$ and $\lambda$, are assumed
to be independent and uncorrelated. Given that both $C_{0}$ and $\lambda$
depend on the diffusion coefficient and degradation rate, this assumption can
easily be violated. Indeed assessing the correlation between $C_{0}$ and
$\lambda$ for measurements taken from 5–25 ss embryos reveals a modest but
significant correlation of $R=0.26$ (Pearson correlation coefficient;
$p=0.001$) (Fig. 2). This is inconsistent with the assumption that $C_{0}$ and
$\lambda$ vary independently.

Reply: We have checked this in our paper using the same dataset [10] but over
the entire time course rather than limited to the first 25 somite stages, and
concluded otherwise (Fig. 4G in [1]). As the data was plotted using vector
graphics, we could extract the plotted point pairs, from which we determined
Pearson’s $R=-0.0061$ ($p=0.94$) and Kendall’s $\tau=0.056$ ($p=0.26$),
suggesting that any correlation between $C_{0}$ and $\lambda$, if it exists,
is negligible.

This is confirmed by our cell-based simulations (Fig. 1A): The values of
$C_{0}$ and $\lambda$, that we determined from fitting the numerically
simulated gradients with independent noise in the three kinetic parameters,
are largely uncorrelated: $R=-0.05$ ($p=0.11$) and $\tau=-0.03$ ($p=0.12$).
This confirms the validity of our assumption.

The correlation that the authors observe at early stages ($R=0.15$–$0.32$)
lies within the expected range of a fitting artifact arising from the
technical limitations of 8-bit imaging, as detailed in our main reply (Fig.
1B–D). Thus, the author’s data provides no evidence for real correlation.

Zagorski et al.: In conclusion, the assumptions inherent to the work of Vetter
and Iber and their decision not to take into account experimental evidence
make their conclusion, that gene expression boundaries in the neural tube are
accurately positioned by a single morphogen gradient, unconvincing.

Reply: None of the points of criticism raised by Zagorski et al. bears
relevance to the conclusions of our article, and as pointed out above, we
clearly state that our work shows that single gradients can, in principle, be
precise enough to encode patterns in the neural tube. We will discuss how the
progenitor domains are defined in forthcoming work.

In summary, we demonstrated that the conclusion of Zagorski et al., that the
positional error of gradients is much higher than that of their readouts, is
the consequence of using different methods to calculate them, and of using
8-bit imaging that can detect fluorescent signals only over a 256-fold range.
When using consistent methods, the positional errors are very similar within
the distance from the source where 8-bit imaging can yield technically sound
results.

### Figure Legends:

Zagorski et al.: Fig. 1. A numerical example of the rescaling error for
exponential gradients with variable $C_{0}$ and $\lambda$.

A. Left: Set of 50 randomly generated exponential morphogen profiles
$C(x)=C_{0}e^{-x/\lambda}$. Mean $\lambda=20$ µm and $C_{0}=1$. $\lambda$ and
$C_{0}$ were varied by adding Gaussian noise with $\mathrm{CV}_{\lambda}=0.2$
and $\mathrm{CV}_{C_{0}}=0.2$. The domain length was randomly selected from
uniform distribution between min $L=100$ µm and max $L=150$ µm. 1 cell
diameter (cd) $=4.9$ µm (as in 1). Dashed horizontal line indicates a
concentration threshold $C_{\theta}=0.1$. At this threshold, the histogram of
positions is shown, the mean position $\operatorname{mean}\\{x_{i,\theta}\\}$
is 48.4 µm from the source, and the positional error is
$\operatorname{std}\\{x_{i,\theta}\\}=9.3\,\text{\textmu{}m}=1.90$ cd.

Right: The profiles are rescaled to the average length
$\operatorname{mean}\\{L_{i}\\}=125.7$ µm by rescaling each $\lambda_{i}$ by a
factor $\operatorname{mean}\\{L_{i}\\}/L_{i}$. In this set,
$\operatorname{mean}\\{x_{i,\theta}\\}=48.8$ µm, and
$\operatorname{std}\\{x_{i,\theta}\\}=10.0\,\text{\textmu{}m}=2.05$ cd. This
indicates that rescaling changed the positional error estimate
$\operatorname{std}\\{x_{i,\theta}\\}$ by 0.15 cd.

B. In Vetter and Iber, a scaling correction is estimated for exponential
profiles without variability. A set of 50 such profiles is shown without and
with rescaling (left and right, respectively). Assuming a uniform distribution
of values at any given concentration threshold, the scaling error increases
with distance to the source and reaches a maximum of 3 cd. Thus, the scaling
error corresponds to $3\xi$, where $\xi$ denotes the relative position of the
bin from the source. For the mean position at $C_{\theta}$ considered here
($\operatorname{mean}\\{x_{i,\theta}\\}=46.7$ µm), $\xi=46.7/125.7=0.37$,
hence the implied scaling correction is 1.11 cd ($=5.4$ µm). This is much
higher than the rescaling error of 0.15 cd that we obtained for the dataset in
A which incorporates realistic variability in $C_{0}$ and $\lambda$. Note that
for an opposing gradient using the same coordinate system, the scaling
correction should be $3(1-\xi)$. Yet, in Fig. 2C and E of Vetter and Iber2,
the same correction of $3(1-\xi)$ is incorrectly applied to both the GBS-GFP
and pSmad gradients. Had we used $3(1-\xi)$ as in 2, the implied correction
would be 1.8 cd and be even more overestimated compared to the actual
rescaling error.

Reply: 50 gradients are not sufficient to obtain precise enough statistics. We
repeated the authors’ analysis with $n=\numprint{10000}$ gradients and
obtained a scaling correction of $0.32\pm 0.01$ cd (estimate $\pm$ SE) instead
of 0.15 cd as purported. The scaling error is thus greater than the authors
claim. Moreover, the example mentioned by the authors is specific to
variability in the gradient parameters that is likely dominated by technical
limitations ($\mathrm{CV}_{\lambda,C_{0}}=0.2$), as we showed in our reply
here (Fig. 1C) and in [1]. Repeating the same analysis along the entire
domain, but taking as gradient variability the relationships inferred for
molecular noise at $\mathrm{CV}_{p,d,D}=0.3$ (Fig. 8F,I in [1]),

$\mathrm{CV}_{\lambda}=\sqrt{\frac{0.270\,\text{\textmu{}m}}{L}}$

and

$\mathrm{CV}_{C_{0}}=0.1773+\frac{L}{\numprint{12871}\,\text{\textmu{}m}}+\left(\frac{L}{5659\,\text{\textmu{}m}}\right)^{2},$

yields domain scaling errors whose purely additive contribution is about 61%
of the $3\xi$ approximation we had used (Fig. 2). In light of this new data,
the scaling correction to the positional error estimated in [2] is likely
somewhat smaller than shown in Fig. 2 in [1]. However, the corrected
positional errors were previously even lower than the reported ones. We had
included the analysis with a simple rough estimate for the scaling correction
to address the remaining minor gap between the positional errors inferred for
the gradients and the one that had been reported for the readouts. Domain
scaling is a subordinate issue in comparison to the other limitations in [2].

Zagorski et al.: Fig. 2. Correlation between the amplitude $C_{0}$ and decay
length $\lambda$ of measured Shh gradients. $C_{0}$ and $\lambda$ are obtained
from exponential fits to the measured Shh ligand gradients from Cohen et al11.
Here, Shh profiles were assigned to developmental stages (designated ss for
somite stage) based on their DV length as described in Zagorski et al. Dashed
lines are linear fits to the data. For each stage, the Pearson correlation
coefficient R is shown in the plot. For the pooled set of profiles between ss5
and ss25, the correlation coefficient is 0.26. Only stages up to 25 ss are
shown.

Reply: The apparent correlation between gradient amplitude and decay length
can be accounted for by technical limitations when the gradient amplitude is
still small. For details, see our reply in the main text and Fig. 1.

## Supplementary References

## References

  * Vetter and Iber [2022] R. Vetter and D. Iber.  Precision of morphogen gradients in neural tube development.  _Nat. Commun._ , 13:1145, 2022.  10.1038/s41467-022-28834-3.
  * Zagorski et al. [2017] M. Zagorski, Y. Tabata, N. Brandenberg, M. P. Lutolf, G. Tkac̆ik, T. Bollenbach, J. Briscoe, and A. Kicheva.  Decoding of position in the developing neural tube from antiparallel morphogen gradients.  _Science_ , 356:1379–1383, 2017.  10.1126/science.aam5887.
  * Kicheva et al. [2014] A. Kicheva, T. Bollenbach, A. Ribeiro, H. Pérez Valle, R. Lovell-Badge, V. Episkopou, and J. Briscoe.  Coordination of progenitor specification and growth in mouse and chick spinal cord.  _Science_ , 345, 2014.  10.1126/science.1254927.
  * Yu et al. [2013] K. Yu, S. McGlynn, and M. P. Matise.  Floor plate-derived sonic hedgehog regulates glial and ependymal cell fates in the developing spinal cord.  _Development_ , 140:1594–1604, 2013.  10.1242/dev.090845.
  * Iber and Vetter [2022] D. Iber and R. Vetter.  Relationship between epithelial organization and morphogen interpretation.  _Curr. Opin. Genet. Dev._ , 75:101916, 2022.  10.1016/j.gde.2022.101916.
  * Long et al. [2023] Y. Long, R. Vetter, and D. Iber.  2d effects enhance precision of gradient-based tissue patterning.  _iScience_ , 26:107880, 2023.  10.1016/j.isci.2023.107880.
  * Mizutani et al. [2006] C. M. Mizutani, N. Meyer, H. Roelink, and E. Bier.  Threshold-Dependent BMP-Mediated Repression: A Model for a Conserved Mechanism That Patterns the Neuroectoderm.  _PLoS Biol._ , 4:e313, 2006.  10.1371/journal.pbio.0040313.
  * Eldar et al. [2003] A. Eldar, D. Rosin, B.-Z. Shilo, and N. Barkai.  Self-enhanced ligand degradation underlies robustness of morphogen gradients.  _Dev. Cell_ , 5:635–646, 2003.  10.1016/S1534-5807(03)00292-2.
  * Adelmann et al. [2023] J. A. Adelmann, R. Vetter, and D. Iber.  Patterning precision under non-linear morphogen decay and molecular noise.  _eLife_ , 12:e84757, 2023.  10.7554/eLife.84757.
  * Cohen et al. [2015] M. Cohen, A. Kicheva, A. Ribeiro, R. Blassberg, K. M. Page, C. P. Barnes, and J. Briscoe.  Ptch1 and Gli regulate Shh signalling dynamics via multiple mechanisms.  _Nat. Commun._ , 6:6709, 2015.  10.1038/ncomms7709.

See pages - of anc/personal_communication