# Precision Radiotherapy via Information Integration of Expert Human Knowledge
and AI Recommendation to Optimize Clinical Decision Making

Wenbo Sun<EMAIL_ADDRESS>Dipesh Niraula<EMAIL_ADDRESS>Issam
El Naqa<EMAIL_ADDRESS>Randall K Ten Haken<EMAIL_ADDRESS>Ivo D
Dinov<EMAIL_ADDRESS>Kyle Cuneo<EMAIL_ADDRESS>Judy (Jionghua) Jin
<EMAIL_ADDRESS>Department of Industrial and Operations Engineering,
University of Michigan, Ann Arbor  Department of Machine Learning, H. Lee
Moffitt Cancer Center & Research Institute, Tampa, FL  Department of Radiation
Oncology, University of Michigan, Ann Arbor  Department of Computational
Medicine and Bioinformatics, University of Michigan, Ann Arbor

###### Abstract

In the precision medicine era, there is a growing need for precision
radiotherapy where the planned radiation dose needs to be optimally determined
by considering a myriad of patient-specific information in order to ensure
treatment efficacy. Existing artificial-intelligence (AI) methods can
recommend radiation dose prescriptions within the scope of this available
information. However, treating physicians may not fully entrust the AI’s
recommended prescriptions due to known limitations or when the AI
recommendation may go beyond physicians’ current knowledge. This paper lays
out a systematic method to integrate expert human knowledge with AI
recommendations for optimizing clinical decision making. Towards this goal,
Gaussian process (GP) models are integrated with deep neural networks (DNNs)
to quantify the uncertainty of the treatment outcomes given by physicians and
AI recommendations, respectively, which are further used as a guideline to
educate clinical physicians and improve AI models performance. The proposed
method is demonstrated in a comprehensive dataset where patient-specific
information and treatment outcomes are prospectively collected during
radiotherapy of $67$ non-small cell lung cancer patients and retrospectively
analyzed.

###### keywords:

Precision medicine , Decision making , Artificial intelligence , Computer
model calibration , Gaussian process modeling

††journal: International Journal of Medical Informatics

Develop an integrative system to help physicians design radiotherapy by
combining human knowledge and AI recommendation.

Quantify the uncertainty of the treatment outcome based on black-box AI
algorithms and physicians’ prescriptions.

Use the data analytic results to educate physicians and improve the AI
recommendations.

Demonstrate the proposed method using real patient data from radiotherapy.

##  1 Introduction

There is a growing need for precision radiotherapy since the precision
medicine era was highlighted by the Joint American Society for Radiation
Oncology (ASTRO) / National Cancer Institute (NCI) workshops in Bethesda, MD
(Benedict et al., 2016; Ashton et al., 2018). Since the treatment outcome of a
radiotherapy depends on the patient information and dose prescription, the
planned radiation dose prescription should be optimized based on the patient-
specific information to achieve the optimal treatment outcomes (El Naqa et
al., 2018) for precision radiotherapy. The radiotherapy effectiveness is often
evaluated by checking whether the local control (LC) of the tumor is assured
and the main side effects such as radiation-induced pneumonitis (RP2) is
mitigated when the patients complete all planned radiotherapy treatments. It
is desired that the tumor tissue is LC’ed while the normal tissue gets least
disturbed due to RP2. As the treatment outcome has its inherent randomness, it
is expected to maximize the probability of LC and minimize the probability of
RP2.

The radiotherapy plan is designed following 3 successful dose escalation
protocols (Luo et al., 2018) as shown in Figure 1. Each radiotherapy plan is
divided into $3$ stages, in which the patient’s state variables (simplified as
“patient variables” in the later context) $s_{1}$, $s_{2}$ and $s_{3}$ are
consistently monitored. In each stage, fractions of radiation doses with dose
per fraction $a_{1}$, $a_{2}$, and $a_{3}$ were delivered at the three stages,
respectively. In the protocols, physicians fix $a_{1}$ and $a_{2}$ as $2$
Gy/frac in the first two stages, and adjust $a_{3}$ based on the patient
variables $s_{3}$. After the three-stage radiotherapy plan, the treatment
outcomes are evaluated by two binary variables $y_{1}\in\\{0,1\\}$ and
$y_{2}\in\\{0,1\\}$, where $y_{1}=1$ indicates LC while $y_{2}=1$ indicates
RP2. In practice, $a_{3}$ is often determined based on physicians’ expert
human knowledge and may not always achieve the optimal treatment outcome. Such
a limitation motivates the development of the precision radiotherapy
techniques.

Figure 1: Illustration of the dose escalation protocol utilized in the study

The objective of the precision radiotherapy process is to adjust $a_{3}$ for
optimal treatment outcome, that is to optimize the trade-off between the high
probability of LC (denoted by $Prob[y_{1}=1]$) and low probability of RP2
(denoted by $Prob[y_{2}=1]$). The optimally planned radiation dose
prescription $a_{3}$ (simplified called prescription in latter context) is
usually searched via artificial intelligence (AI) approaches. There are two
main types of approaches: direct methods and indirect methods. The direct
methods have been proposed by Zhao et al. (2012), which search for the optimal
prescription directly based on the weighted outcomes in learning (called
O-learning) without requiring a prediction model for the potential outcomes of
future prescriptions. Thus, there is no need to fit a prediction model, and
the process is model free. In contrast, indirect methods take a two-step
analysis, i.e., first to fit a model that can accurately predict the outcomes
in the prescription space, and then to search for an optimal prescription
based on the predicted outcomes, that is a model based approach to improve
performance. Q-learning is one of the most popularly used also in these
indirect, or model-based methods (Qian and Murphy, 2011; Chakraborty, 2013;
Moodie et al., 2014). Considering the patient variables are dynamically
changing over different stages of radiotherapy, it is more appropriate to use
the indirect method for precision radiotherapy, that is, to predict the
probability of the treatment outcome
$\left(Prob[y_{1}=1],Prob[y_{2}=1]\right)$, based on which optimal
prescription $a_{3}$ is obtained by optimizing the predicted probabilities
$\left(Prob[y_{1}=1],Prob[y2=1]\right)$.

To develop the indirect methods, the treatment outcome will be predicted via
existing machine learning approaches such as penalized least square (Qian and
Murphy, 2011), Bayesian networks (Pearl et al., 2009), deep neural networks
(Goodfellow et al., 2016), etc. Since the effect of radiotherapy among
multiple stages is a complex nonlinear function on the patient variables, deep
neural networks are usually selected (Tseng et al., 2017). When the Q-learning
is integrated with the deep neural network, it is also known as deep Q
network. A major limitation of deep Q-learning is its instability due to
limited training data – a small error in prediction of the treatment outcome
may significantly affect the final AI recommendation, especially when the
sample size of patients is very limited as typically the case. Another concern
of the deep Q-learning is the interpretability – the deep neural network is
usually constructed as a black box, making it difficult to interpret the
relationship between the predictors and outcomes in a straightforward manner.

Ideally, when the Q-learning approach always provides a prescription that is
close to the oracle’s benchmark, clinical physicians feel comfortable to
follow the AI recommendation. However, when the treatment outcome is not
accurately predicted due to small and incomplete datasets, learning from the
results of practitioners’ decisions or biased inference and evaluation
processes (Rich and Gureckis, 2019), physician’s prescription may result in
preferable treatment outcomes than those are aroused by following the AI
recommendation. In this paper, instead of seeking the origin of the prediction
bias, we aim to develop a systematic approach to guide the clinical physicians
when their prescriptions are different from the AI recommendations. We propose
to use the commonly used tool for uncertainty quantification, Gaussian process
(GP) model (Rasmussen, 2003), to compensate the gap between the AI predictions
of the treatment outcomes and their true values under the computer model
calibration framework (Kennedy and O’Hagan, 2001). An integrative decision
system will be developed: when the physicians’ prescriptions are desired, the
system suggests how the AI algorithm can be improved. When the AI
recommendations provide better treatment outcomes, the system helps physicians
make better decision for future patients. The flowchart of the proposed method
is provided in Figure 2.

Figure 2: Information flow and decision support system framework

##  2 Notations and problem formulation

To start with, we would like to introduce some mathematical notations that
will be used to elaborate the proposed method. Let $T$ denote the total number
of stages of radiotherapy ($T=3$ in the dose escalation protocols (Luo et al.,
2018)). Based on Figure 1, for a specific patient, the patient variables
$s_{t}$ are observed at stage $t=1,...,T$ when the patient is in the hospital.
When the patient exits from the hospital, the treatment outcome
$(y_{1},y_{2})$, instead of the patient variables $s_{T+1}$, is recorded. At
stages $t=1,...,T-1$, radiation treatments with a known dose per fraction are
delivered to the patient. At stage $T$, the dose per fraction $a_{T}$ needs to
be adjusted based on $s_{T}$ to achieve the optimal treatment outcome.

As was mentioned in Section 1, the optimal treatment outcome refers to the
trade-off between the high LC probability and low RP2 probability. Here we
modify the “P+” reward function in Luo et al. (2018) to the smoothed version
$R$, which is defined as:

$R=-10\times\left(\left(1-Prob\left[y_{1}=1\right]\right)^{8}+\left(\frac{Prob\left[y_{2}=1\right]}{0.57}\right)^{8}\right)^{1/8}+3.281.$
(1)

Given the reward function, we define the action value function $Q$ as a Markov
Decision Process (MDP) as:

$Q_{\pi}\left(s_{T},a_{T}\right)=\mathbb{E}_{\pi}\left[R_{T+1}+\gamma
R_{T+2}+\gamma^{2}R_{T+3}+...\middle|s_{T},a_{T}\right],$ (2)

where $R_{t}$ represents the reward at time $t$ and $\gamma$ denotes the
discount factor. By setting $R_{T+1}=R$, we can further derive the action
value function as:

$Q_{\pi}\left(s_{T},a_{T}\right)=\mathbb{E}\left[R\middle|s_{T},a_{T}\right]\\\
+\gamma\int_{s_{T+1}}ds_{T+1}p\left[s_{T+1}\middle|s_{T},a_{T}\right]\int_{a_{T+1}}da_{T+1}\pi\left[a_{T+1}\middle|s_{T+1}\right]Q_{\pi}\left(s_{T+1},a_{T+1}\right),$
(3)

with $p$ representing the probability density function, policy
$\pi\left[a\middle|s\right]$ representing probability density of choosing $a$,
and $s_{T+1}$ denoting the true, yet unknown patient variables after
treatment. With $Q$ defined, the optimization problem is formulated as:

$a_{T}^{*}=\pi^{*}\left(s_{T}\right)=arg\,max_{a_{T}\in\mathcal{A}}Q\left(s_{T},a_{T}\right).$
(4)

Solving the above optimization problem yields the optimal dose decision for AI
recommendation.

Equation (3) involves two unknown terms that should be modeled and estimated,
$p\left[s_{T+1}\middle|s_{T},a_{T}\right]$ and
$E\left[R\middle|s_{T+1}\right]$. The former term is named by “transition
function” $f$, which updates the transition of the patient variables from the
previous stage to the current stage under given dose per fraction $a_{T}$. The
latter term evaluates the reward function based on $Prob\left[y_{1}=1\right]$
and $Prob\left[y_{2}=1\right]$, which further requires an “evaluation
function” $g$ that predicts the probability of the treatment outcome based on
$s_{T}$. The relationship between the two functions $f$ and $g$ and the other
variables is illustrated in Figure 1. Quantifying $f$ and $g$ requires a
clinical trial dataset that collects radiotherapy records. In the dataset, the
patient variables are stored in
$S=\left\\{s_{i,t};i=1,...,n;t=1,...,T\right\\}$ where
$s_{i,t}\in\mathbb{R}^{q}$ represent the patient $i$’s patient variables at
stage $t$ with dimension of $q$. At each stage $t\in\\{1,...,T\\}$, the
physicians deliver radiation treatments with dose per fraction $a_{i,t}$ to
the patient $i$. When the radiotherapy reaches to the end, the treatment
outcome are recorded as $y_{i,1}$ and $y_{i,2}$. Suppose the estimates are
denoted by $\hat{f}$ and $\hat{g}$. The reward function $R$ in Equation 3 can
be estimated as $\hat{g}\left(\ \hat{f}\left(s_{T},a_{T}\right)\right)$. The
AI recommendation is hereby computed as

${\hat{a}}_{T}=arg\,max_{a_{T}\in\mathcal{A}}\ \hat{g}\left(\
\hat{f}\left(s_{T},a_{T}\right)\right).$ (5)

The risks of following the AI recommendations in the real practice are caused
by inaccurate $\hat{f}$ and $\hat{g}$ – when these two estimators do not well
approximate the true $f$ and $g$, the resultant ${\hat{a}}_{T}$ will deviate
from the true optimal prescription $a_{T}^{*}$, and may lead to poor treatment
outcome. To help physicians determine the prescription with a given AI
recommendation, a methodology will be developed in the following framework:
firstly, the uncertainty of the treatment outcome is quantified through the GP
modelling of $f$ and $g$. Secondly, a comparison between the treatment
outcomes is conducted via a hypothesis testing, assisting the physician to
determine whether to use their own prescription or trust in the AI
recommendation. Lastly, the comparison results will be fed back to improve
both the human expert knowledge and AI recommendation.

##  3 Methodology

###  3.1 Estimation of the transition function

The transition function $f$ provides the one-step-ahead prediction of the
patient variables given the current patient variables and the corresponding
prescription. $\left\\{\left(s_{t},a_{t}\right),t=1,...,T-1\right\\}$ are
treated as predictors and $\left\\{s_{t+1},t=1,...T-1\right\\}$ are considered
as the corresponding responses. Once $f$ is trained, the state variables at
the final stage, $s_{T+1}$, can be predicted based on $s_{T}$ and $a_{T}$. A
point estimator of $f$ can be trained as a DNN by utilizing the data generated
from a generative adversarial network (GAN) as in Goodfellow et al. (2014).
However, since DNN is a black box, there is no straightforward way to quantify
the uncertainty of $f$, which is crucial in evaluating the accuracy of the
resultant AI recommendation. To do this, we attach a GP-based bias term to the
DNN following the computer model calibration framework Kennedy and O’Hagan
(2001) as:

$s_{t+1}=f\left(s_{t}\right)=\eta_{s}\left(s_{t}\right)+\delta_{s}\left(s_{t}\right),t=1,...,T,$
(6)

where $\eta_{s}\left(s_{t}\right)$ represents the point estimator obtained
from DNN, and $\delta_{s}\left(s_{t}\right)$ represents the model bias. The
bias
$\delta_{s}\left(s_{t}\right)=\\{\delta_{s}^{\left(1\right)}\left(s_{t}\right),...,\delta_{s}^{\left(q\right)}\left(s_{t}\right)\\}$
is a $q$-dimensional random function, each dimension of which is assumed to
follow a GP as:

$\delta_{s}^{\left(j\right)}\left(\cdot\right)\sim\mathcal{GP}\left(0,\frac{1}{\lambda_{\delta}^{\left(j\right)}}c_{\delta}^{\left(j\right)}\left(\cdot,\cdot\right)\right),j=1,...,q,$
(7)

where $c_{\delta}^{\left(j\right)}\left(\cdot,\cdot\right)$ is the covariance
function and $\lambda_{\delta}^{\left(j\right)}$ is the precision parameter of
GP. The covariance function is assumed to take the squared exponential (SE)
form as:

$\displaystyle
c_{\delta}^{\left(j\right)}\left(\left(s,a\right),\left(s^{\prime},a^{\prime}\right)\right)$
$\displaystyle=$
$\displaystyle\prod_{k=1}^{q}\exp{\left(-\beta_{\delta,k,j}\left(s^{\left(k\right)}-s^{\left(k\right)\prime}\right)^{2}\right)}$
(8) $\displaystyle\times$
$\displaystyle\exp{\left(-\beta_{\delta,q+1,j}\left(a-a^{\prime}\right)^{2}\right)}.$

Here $s^{(k)}$ denotes the $k-$th dimension of the vector $s$. Following the
parameter estimation and Gaussian process prediction techniques shown in A, we
are able to generate prediction of $s_{T+1}$, denoted by $\widehat{S}_{T+1}$.
The prediction $\widehat{S}_{T+1}$ can be further used to predict the reward
function.

###  3.2 Estimation of the evaluation function

The evaluation function $g$ predicts the reward function $R$ based on the
predictions of $s_{T+1}$. Since $R$ is determined by
$Prob\left[y_{1}=1\right]$ and $Prob\left[y_{2}=1\right]$, it is sufficient to
predict the two probabilities $Prob\left[y_{1}=1\right]$ and
$Prob\left[y_{2}=1\right]$, or equivalently, the logit function of the two
probabilities, denoted by $h_{1}$ and $h_{2}$, respectively. Here we assume
the two random functions $h_{1}$ and $h_{2}$ to follow GPs as

$\displaystyle h_{1}$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}\left(0,\frac{1}{\lambda_{h_{1}}}c_{h_{1}}\left(\cdot,\cdot\right)\right),$
$\displaystyle h_{2}$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}\left(0,\frac{1}{\lambda_{h_{2}}}c_{h_{1}}\left(\cdot,\cdot\right)\right),$
(9)

where $\lambda_{h_{1}}$ and $\lambda_{h_{2}}$ are the precision parameters,
$c_{h_{1}}\left(\cdot,\cdot\right)$ and $c_{h_{1}}\left(\cdot,\cdot\right)$
are the covariance functions constructed under the squared exponential (SE)
form:

$c_{h_{j}}\left(s,s\prime\right)=\prod_{k=1}^{q}exp{\left(-\beta_{h,k,j}\left(s^{\left(k\right)}-s^{\left(k\right)^{\prime}}\right)^{2}\right)}.$
(10)

However, the two GPs cannot be directly trained because the logit of $LC$ and
$RP2$ in the dataset, denoted by $H_{1}$ and $H_{2}$, are not observed.
Instead, we only observe the binary labels $Y_{1}=\\{y_{i,1};i=1,...,n\\}$ and
$Y_{2}=\\{y_{i,2};i=1,...,n\\}$. As a solution, we utilize the Laplace
Approximation technique (Tierney and Kadane, 1986) to provide the estimation
of GP parameters and the prediction of the $LC$ and $RP2$ probabilities.
Detailed derivations can be found in B. Combination of the evaluation function
and transition function generate the prediction of $R$ for any given patient
variable $s_{T}$ and dose prescription $a_{T}$, which can be used to search
for the optimal dose prescription for AI recommendation via the deep
Q-learning algorithm.

###  3.3 Integration of physician’s prescription and AI recommendation

Now we aim to develop an integrative system that can help make prescription
decisions and improve both physician’s prescription and AI recommendation in
the future. The prerequisite is to quantify the uncertainty of $R$ for given
$s_{T}$ and $a_{T}$, that is to develop a confidence interval estimator of the
composite function $g\left(f\left(\cdot\right)\right)$. Detailed derivations
can be found in C.

Next, we generate the predictive distribution of the reward function $R$ via
Monte Carlo simulations since the reward function is not linear with respect
to the LC and RP2 probabilities. We firstly generate simulated $Prob[LC]$ and
$Prob[RP2]$ from Equation (26). Then compute the simulated reward function
based on its definition in Equation (1). Suppose the AI recommendation is
$a^{A}(s_{T-1})$ for patient variable $s_{T}$, the physician’s prescription is
$a^{P}(s_{T})$. We denote the simulated reward functions based on the AI
recommendation by

$\tilde{r}\left(s_{T},a^{A}\left(s_{T}\right)\right).$ (11)

The simulated reward functions based on the physician’s prescription is
denoted by

$\tilde{r}\left(s_{T},a^{P}\left(s_{T}\right)\right).$ (12)

Given the Monte Carlo simulated reward functions, a student-t test is
conducted to test whether

$\tilde{r}\left(s_{T},a^{A}\left(s_{T}\right)\right)$ (13)

are significantly greater than

$\tilde{r}\left(s_{T},a^{P}\left(s_{T}\right)\right).$ (14)

The system selects the AI recommendation when the null hypothesis is rejected,
or equivalently, the p-value is smaller than $0.05$. Otherwise, the system
will recommend physicians to follow their original prescription.

Thirdly, we would like to utilize the hypothesis testing results to improve
the physician’s prescription and AI recommendation for future patients. For
physicians, they need to know the amount of dose that need to be compensated
comparing to that from their human expert knowledge. To do this, we assume the
physician’s dose bias
$\delta_{a}\left(s_{T}\right)=a^{A}\left(s_{T}\right)-a^{P}\left(s_{T}\right)$
to follow a GP with considering $s_{T}$ as inputs. The GP can be trained
following the same procedure as in subsection 3.1. Here we only utilize the
samples whose corresponding p-values are smaller than $0.05$, implying that
the AI recommendation provides better treatment outcome than the original
physician’s prescription at those samples.

When AI recommendations are not reliable, physicians should be warned. To do
this, the confidence intervals of $Prob[y_{j}=1|s_{T},a^{A}(s_{T})]$ and
$Prob[y_{j}=2|s_{T},a^{A}(s_{T})]$ are visualized for physicians. The large
confidence intervals imply that the AI recommendations are not reliable. In
these cases, it is recommended that physicians should insist on their own
prescription and be careful when delivering the dose prescription. The
detailed guidelines will be provided in Section 4.

##  4 Case study in radiotherapy of lung cancer

In this section, the proposed method is demonstrated in a dataset collecting
the radiotherapy records of $67$ patients in $4$ weeks. The dataset is part of
the clinical trial data with real patients undergoing treatment. Part of the
original dataset was published here in Kong et al. (2017). In the dataset, a
total of $250$ potential patient variables are collected as the patient
information before and during radiotherapy. They include dosimetric variables,
clinical factors, circulating microRNAs, single-nucleotide polymorphisms
(SNPs), circulating cytokines, positron emission tomography (PET) imaging
radiomics features, etc. To reduce the dimension of patient variables to a
feasible number, Bayesian Network technique (Luo et al., 2018, 2017) is
applied to search for the most related patient variables to LC and RP2 where
the selected variables are identified as the Markov blankets (MBs) of LC and
RP2. As a result, 12 out of 250 patient variables are selected as the
predictors for LC and RP2. They are il4, il10, il5, ip10, MTV, GLSZM_LZLGE,
GLSZM_ZSV, Tumor_gEUD, Lung_gEUD, Rs2234671, Rs238406, Rs1047768. To avoid
sparsity in the space of training dataset, extreme values are removed by
truncating all the patient variables at their $70\%$ quantiles. Then the
patient variables are scaled into $[0,1]$ to train the transition and
evaluation functions. Detailed descriptions of the $12$ selected patient
variables can be found in D.

After pre-processing the data, we start to train the transition function $f$
as a combination of DNN and GP. When training the DNN, GAN is implemented to
augment the sample size to ensure the fitting accuracy. A comparison between
the predictions from the original deep neural network and the Gaussian process
model is shown in Table 1. The last three patient variables are not predicted
because they are constant. The cross-validation mean square errors (MSE) of
the week 2 and week 4 data are calculated as the evaluation criterion. The
relative improvement is computed as

$RI_{j}=\frac{\sum_{t=1}^{2}\|\hat{\mu}_{t+1}^{(j)}(S_{t},A_{t})-\eta_{s}^{(j)}(S_{t},A_{t})\|^{2}}{\sum_{t=1}^{2}\|\eta_{s}^{(j)}(S_{t},A_{t})-S_{t+1}^{(j)}\|^{2}}.$
(15)

Since DNN already achieves a relatively high accuracy, the GP model does not
significantly improve the accuracy. However, the Gaussian process model
provides the predicted confidence interval of the patient variables, which can
be further used to quantify the uncertainty of the treatment outcome in the
next steps.

Table 1: Cross-validation MSE of the transition function Response | DNN | GP | Relative Improvement
---|---|---|---
il4 | 0.000287 | 0.000110 | $72.29\%$
il10 | 0.000750 | 0.000691 | $64.16\%$
il5 | 0.000394 | 0.000367 | $28.02\%$
ip10 | 0.000038 | 0.000002 | $89.79\%$
mtv | 0.010464 | 0.009343 | $44.60\%$
GLSZM_LZLGE | 0.014043 | 0.016114 | $31.86\%$
GLSZM_ZSV | 0.001932 | 0.002155 | $21.31\%$
Lung gEUD | 0.000182 | 0.000100 | $58.10\%$
Tumor gEUD | 0.000480 | 0.000442 | $26.39\%$

Next, the predicted week 6 patient variables are used to train the evaluation
function. We assign a non-informative prior to the precision parameters
$\lambda_{h_{1}}$ and $\lambda_{h_{2}}$ as

$p\left(\lambda_{h_{j}}\right)\propto\lambda_{h_{j}}^{-1},j=1,2,$ (16)

maximize the likelihood in Equation (22). through a grid search, and obtain
the two precision parameters’ values as $\lambda_{h_{1}}=\lambda_{h_{2}}=4$.
Then the other GP parameters are estimated with the optimal precision
parameters. With the parameters estimated, the LC and RP2 probabilities can be
predicted. Figure 3 compares the detailed predicted probabilities on sorted
samples. The GP model provides probabilities that are closer to the true
binary labels. The cross-entropy (Goodfellow et al., 2016) is calculated for
quantitative comparison. For LC, the GP model reduces the cross-entropy from
$35.32$ to $33.58$. For RP2, the GP model reduces the cross-entropy from
$52.19$ to $44.76$. More importantly, the GP model provides the confidence
interval estimator of the LC and RP2 probabilities, which can be used to
quantify the uncertainty of the treatment outcome based on physicians’
prescriptions and AI recommendations.

Figure 3: Comparisons between the predicted probabilities and the true binary
LC/RP2 labels. The triangular represents the DNN predictions. The cross
represents the prediction after GP calibration. The plus shows the true LC/RP2
labels. The dashed line depicts the $0.5$ probability threshold. (a)
comparisons of the predicted LC probabilities. $20$ samples have LC=0 while
$47$ samples have LC=1. (b) comparisons of the predicted RP2 probabilities.
$51$ samples have RP2=0 while $16$ samples have RP2=1.

As an illustration, the uncertainty quantification results of the treatment
outcome for two patients are visualized in Figure 4. GP enables our model to
evaluate the uncertainty of the treatment outcome at different doses. At the
observed data point, the uncertainty of the treatment outcome based on the AI
recommendation is often greater than that from the physician’s prescription.
This is because GP has higher prediction accuracy around the point where the
dose prescription and patient variables are known. For the first patient shown
in Figure 4 (a) and (b), the AI algorithm suggests a higher dose prescription,
which may result in better LC and similar RP2. To determine whether the AI
recommendation or physician’s prescription is better, we compute the p-value
following the procedure in subsection 3.3. The resultant p-value of $0$
suggests to follow the AI recommendation for the future patients with $s_{T}$.
For the second patient shown in Figure 4 (c) and (d), since the treatment
outcome from AI recommendation has a large uncertainty, we would recommend to
use physician’s prescription.

Figure 4: Uncertainty quantification of the treatment outcome for two
patients. (a) Uncertainty quantification of the treatment outcome of different
doses for the first patient. The solid, dashed, and dotted-dashed lines
represent the GP prediction, $\pm 2$ standard deviation confidence interval of
GP prediction, and DNN prediction, respectively. The triangular depicts the
physician’s prescription. The square depicts the AI recommendation. (b)
Predicted LC and RP2 probabilities for the first patient based on physician’s
prescription and AI recommendation. The plus represents the mean predictions
of LC and RP2 based on the AI recommendation. The cross shows the mean
predictions of LC and RP2 based on Physician’s prescription. The two ellipses
are the $\pm 2$ confidence region of the LC/RP2 predictions. (c) Uncertainty
quantification of the treatment outcome of different doses for the second
patient. (d) Predicted LC and RP2 probabilities for the second patient based
on physician’s prescription and AI recommendation.

We conduct the analysis on all the $67$ patients. The integrative system
suggests that the AI recommendation provides better treatment outcomes for
$35$ patients based on the p-values. We further use these samples to predict
the dose compensation, which may help physicians make better dose
prescriptions. Since the physicians mainly determine the dose prescription
based on the two patient variables “Tumor gEUD” and “Lung gEUD”, we construct
a Gaussian process to predict the dose compensation $a^{A}-a^{P}$ based on
these two patient variables, and visualize the predictions in a 2-dimensional
plot with respect to the two patient variables in Figure 5. The warm
background color implies that the physicians’ prescriptions are overall
conservative comparing to the recommended doses. It is also worth noting that
the AI algorithm suggests to deliver higher doses for small Tumor_gEUD, and
lower doses for higher Tumor_gEUD.

Figure 5: Visualization of the dose compensation for physicians’
prescriptions. The regions with warm colors depict the combination of the
patient variables that require higher doses comparing to the current
physicians’ prescriptions. The two patient variables have been transformed to
their original scales for easier interpretation. The regions with cold colors
implies that the physicians’ should reduce their dose prescription for the
patients with the corresponding patient variables. The black crosses
illustrate the training samples.

##  5 Acknowledgement:

The research is partly supported by NIH grant R01-CA233487.

##  Appendix A Parameter estimation and Gaussian process prediction for the
transition function

Denote the observed patient variables at stage $t$ by $S_{t}$, and its $j$-th
dimension by ${S_{t}}^{\left(j\right)}$. Denote the prescriptions at stage $t$
by $A_{t}$. The model parameters,
$\beta_{\delta}=\\{\beta_{\delta,k,j},k=1,...,q+1,j=1,...,q\\}$ and
$\lambda_{\delta}=\\{\lambda_{\delta}^{\left(j\right)},j=1,...,q\\}$, can be
trained by maximizing the log-likelihood as:

$\left({\hat{\beta}}_{\delta}^{\left(j\right)},{\hat{\lambda}}_{\delta}^{\left(j\right)}\right)=arg\,max\log\,p\left(\beta,\lambda\middle|S_{2}^{\left(j\right)},\ldots,S_{T}^{\left(j\right)};S_{1},\ldots,S_{T-1},A_{1},\ldots,A_{T-1}\right),$
(17)

with

$\displaystyle\log\,p\left(\beta,\lambda\middle|S_{2}^{\left(j\right)},\ldots,S_{T}^{\left(j\right)};S_{1},\ldots,S_{T-1},A_{1},\ldots,A_{T-1}\right)$
(18) $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{t=1}^{T-1}\left(S_{t+1}^{\left(j\right)}-\eta_{s}^{\left(j\right)}\left(S_{t}^{\left(j\right)}\right)\right)^{T}K_{\eta}^{\left(j\right)^{-1}}((S_{t},A_{t}),(S_{t},A_{t}))$
$\displaystyle\times$
$\displaystyle\left(S_{t+1}^{\left(j\right)}-\eta_{s}^{\left(j\right)}\left(S_{t}^{\left(j\right)}\right)\right)-\frac{1}{2}\left|K_{\delta}^{\left(j\right)}((S_{t},A_{t}),(S_{t},A_{t}))\right|-\frac{n}{2}log2\pi,$

where $K^{\left(j\right)}\left(\cdot,\cdot\right)$ is the Gram matrix of the
covariance function $c_{\delta}^{\left(j\right)}\left(\cdot,\cdot\right)$.

With the parameters trained, we can obtain the prediction of the $j-$th
dimension of $s_{T+1}$, $s_{T+1}^{\left(j\right)}$. The prediction is
expressed as a normally distributed random variable, whose distribution
parameters are described as follows:

$\displaystyle s_{T+1}^{\left(j\right)}|s_{T},a_{T}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}\left({\hat{\mu}}_{T}^{\left(j\right)}(s_{T},a_{T}),{\widehat{\Sigma}}_{T}^{\left(j\right)}(s_{T},a_{T})\right),$
$\displaystyle{\hat{\mu}}_{T+1}^{\left(j\right)}(s_{T},a_{T})$
$\displaystyle=$
$\displaystyle\eta_{s}^{\left(j\right)}\left(s_{T}\right)+{K^{\left(j\right)}}^{T}\left(s_{T},\cup_{t=1}^{T-1}S_{t}\right)K_{\delta}^{\left(j\right)^{-1}}\left(\cup_{t=1}^{T-1}S_{t},\cup_{t=1}^{T-1}S_{t}\right)$
$\displaystyle\times$
$\displaystyle\left(\cup_{t=1}^{T-1}S_{t+1}^{\left(j\right)}-\eta_{s}^{\left(j\right)}\left(\cup_{t=1}^{T-1}S_{t}\right)\right),$
$\displaystyle{\widehat{\Sigma}}_{T+1}^{\left(j\right)}(s_{T},a_{T})$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda_{s}^{\left(j\right)}}\\{K_{\delta}^{\left(j\right)}\left(s_{T},s_{T}\right)-{K_{\delta}^{\left(j\right)}}^{T}\left(s_{T},\cup_{t=1}^{T-1}S_{t}\right)$
(19) $\displaystyle\times$ $\displaystyle
K_{\delta}^{\left(j\right)^{-1}}\left(\cup_{t=1}^{T-1}S_{t},\cup_{t=1}^{T-1}S_{t}\right)K_{\delta}^{\left(j\right)}\left(\cup_{t=1}^{T-1}S_{t},s_{T}\right)\\}.$

##  Appendix B Parameter estimation and Gaussian process prediction for the
transition function

To predict the binary labels $Y_{1}$ and $Y_{2}$ without the observations of
$h_{1}$ and $h_{2}$, we would like to marginalize the underlying logit of $LC$
and $RP2$ out, and hence calculate the densities of $h_{1}$ and $h_{2}$ via

$p\left(h_{j}\middle|\widehat{S}_{T+1},Y_{j},s_{T+1}\right)=\int{p\left(h_{j}\middle|H_{j},\widehat{S}_{T+1},Y_{j},s_{T+1}\right)p\left(H_{j}\middle|\widehat{S}_{T+1},Y_{j}\right)dh_{j}},j=1,2,$
(20)

where
$\widehat{S}_{T+1}=\left\\{\hat{\mu}_{T+1}^{(1)}(S_{T}),...,\hat{\mu}_{T+1}^{(q)}(S_{T})\right\\}$.
The integration cannot be computed because the term
$p\left(H_{j}\middle|\widehat{S}_{T+1},Y_{j}\right)$ is intractable.
Therefore, Laplace approximation (Tierney and Kadane, 1986) is implemented to
approximate $p\left(H_{j}\middle|\widehat{S}_{T+1},Y_{j}\right)$ with the
normal density
$\mathcal{N}\left(H_{j}\middle|{\widehat{H}}_{j},\Lambda_{j}\right)$. The
distribution parameters are calculated via maximizing the likelihood function
$\Psi_{j}$ as

$\displaystyle{\widehat{H}}_{j}$ $\displaystyle=$ $\displaystyle
arg\,max_{h_{j}}\Psi_{j}\left(h_{j}\right),$ $\displaystyle\Lambda_{j}$
$\displaystyle=$
$\displaystyle\lambda_{h_{j}}{K_{h_{j}}}^{-1}\left(\widehat{S}_{T+1},\widehat{S}_{T+1}\right)+W_{j},$
$\displaystyle W$ $\displaystyle=$
$\displaystyle\nabla\nabla\Psi_{j}\left(h_{j}\right),j=1,2.$ (21)

The likelihood function $\Psi_{j}$ is defined as

$\displaystyle\Psi_{j}\left(h_{j}\right)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\left(y_{i,j}log\frac{\exp{h_{i,j}}}{1+\exp{h_{i,j}}}+\left(1-y_{i,j}\right)log\frac{1}{1+\exp{h_{i,j}}}\right)$
(22) $\displaystyle-$
$\displaystyle\frac{\lambda_{h_{j}}}{2}h_{j}^{T}K_{h_{j}}^{-1}h_{j}+C.$

Here
$\sigma\left(h\right)=exp\left(h\right)/\left(1+exp\left(h\right)\right)$,
$K_{h_{j}}$ is the Gram matrices of the covariance function
$c_{h_{j}}\left(\cdot,\cdot\right)$, $h_{i,j}$ is the $i$-th sample of
$h_{j}$, and $C$ is some constant. Similar to the training process in
Subsection 3.1, the model parameters,
$\beta_{h}=\\{\beta_{h,k,j},k=1,...,q+1,j=1,2\\}$ and
$\lambda_{h}=\\{\lambda_{h_{1}},\lambda_{h_{2}}\\}$ can be estimated by
maximizing the approximated probabilities
$p\left(h_{1}\middle|\widehat{S}_{T+1},Y_{1},s_{T+1}\right)$ and
$p\left(h_{2}\middle|\widehat{S}_{T+1},Y_{2},s_{T+1}\right)$ by substituting
the $H_{j}$ into Equation (20) with ${\widehat{H}}_{j}$.

##  Appendix C Combination of the predictions of the transition and evaluation
functions

The confidence interval can be generated in two steps. First, we propagate the
uncertainty of the estimated $s_{T}$ to that of the estimated $Prob[y_{1}]$
and $Prob[y_{2}]$, resulting in an error-in-variable GP model (Cressie and
Kornak, 2003) for the logit functions of $h_{1}$ and $h_{2}$:

$h_{j}|s_{T},a\sim\mathcal{GP}\left({\hat{\mu}}_{j},{\widehat{\Sigma}}_{j}\right),j=1,2,$
(23)

where

$\displaystyle{\hat{\mu}}_{j}$ $\displaystyle=$
$\displaystyle{\widetilde{K}}_{\
h_{j}}^{T}\left({\hat{\mu}}_{T+1},{\widehat{S}}_{T+1}\right){\widetilde{K}}_{h_{j}}^{-1}\left({\widehat{S}}_{T+1},{\widehat{S}}_{T+1}\right){\widehat{H}}_{j},$
$\displaystyle{\widehat{\Sigma}}_{j}$ $\displaystyle=$
$\displaystyle\frac{1}{\lambda_{h_{j}}}{\widetilde{K}}_{h_{j}}^{T}\left({\hat{\mu}}_{T+1},{\hat{\mu}}_{T+1}\right)$
(24) $\displaystyle-$
$\displaystyle\frac{1}{\lambda_{h_{j}}}{\widetilde{K}}_{h_{j}}^{T}\left({\hat{\mu}}_{T+1},{\widehat{S}}_{T+1}\right){\widetilde{K}}_{h_{j}}^{-1}\left({\widehat{S}}_{T+1},{\widehat{S}}_{T+1}\right){\widetilde{K}}_{h_{j}}\left({\hat{\mu}}_{T+1},{\widehat{S}}_{T+1}\right),$

with ${\hat{\mu}}_{T+1}$ denoting the predicted mean of $s_{T+1}$ in Equation
(19). ${\widetilde{K}}_{h_{j}}$ is the Gram matrix of the covariance function:

${\tilde{c}}_{h_{j}}\left(s,s\prime\right)=\prod_{k=1}^{q}\frac{1}{1+4\beta_{k,j}\sigma_{T+1,k}^{2}}exp{\left(-\frac{1}{1/\beta_{k,j}+4\sigma_{T+1,k}^{2}}\left(s^{\left(k\right)}-s^{\left(k\right)^{\prime}}\right)^{2}\right)},j=1,2,$
(25)

where $\\{\sigma_{T+1,k}^{2},k=1,...,q\\}$ are the realization of
${\widehat{\Sigma}}_{T+1}^{\left(j\right)}$ in Equation (8) with input
$s_{T}$. Applying the delta-method to Equation (24) yields the predictive
distribution of the probabilities of LC and RP2, which are:

$\displaystyle Prob\left[y_{j}=1|s_{T},a\right]$ $\displaystyle\sim$
$\displaystyle\mathcal{N}\left(\frac{exp{\left(\hat{\mu}_{j}\right)}}{1+exp{\left(\hat{\mu}_{j}\right)}},D_{j}{\widehat{\Sigma}}_{j}D_{j}\right),$
$\displaystyle D_{j}$ $\displaystyle=$ $\displaystyle
diag\left\\{\frac{exp{\left(\hat{\mu}_{j}\right)}}{\left(1+{exp{\left(\hat{\mu}_{j}\right)}}\right)^{2}}\right\\},j=1,2.$
(26)

##  Appendix D Detailed descriptions of the 12 selected patient variables

  * 1.

il4 (interleukin 4): Th2 cytokines: (i) Regulates antibody production,
hematopoiesis, and inflammation. (ii) Promotes the differentiation of naïve
helper T cells into Th2 cells. (iii) Decreases the production of Th1 cells
(Ramírez et al., 2013; Schaue et al., 2012; Warltier et al., 2002). Here CD4 T
helper cells are lymphocytes that strongly modulate the response of the immune
system against cancer cells proliferation and tumor growth. They are
classified into Th1 (antitumor) and Th2 (protumor) cells. Th1 and Th2 cells
releases cytokines. Imbalance between Th1 and Th2 helps cancer cell.

  * 2.

il10 (interleukin 10): Th2 cytokines: (i) Inhibits synthesis of Th1 cytokines
such as IFN-$\gamma$ (interferon gamma) and IL2. (ii) Inhibits antigen-
presenting cells (Ramírez et al., 2013; Schaue et al., 2012; Warltier et al.,
2002). IFN-$\gamma$: Th1 cytokines: (i) Enhances the microbicidal function of
macrophages. (ii) Promotes the differentiation of naïve helper T cells into
Th1 cells. (iii) Activates polymorphonuclear leukocytes, cytotoxic T cells,
and NK cells. IL2: Th2 cytokines: (i) Promotes clonal expansion and
development of T and B-lymphocytes. (ii) Induces expression of adhesion
molecules. (iii) Enhances the function of NK cells.

  * 3.

il15 (interleukin 15): Th2 cytokines: (i) Induces activation and cytotoxicity
of NK (natural killer) cells. (ii) Activates macrophages. (iii) Promotes
proliferation and survival of T and B- lymphocytes and NK cells (Ramírez et
al., 2013; Schaue et al., 2012; Warltier et al., 2002).

  * 4.

ip10 (interfeon gamma-induced protein 10): IP-10 is secreted in response to
IFN-$\gamma$ by various cells including monocytes, endothelial and
fibroblasts. (i) Acts as chemoattractions for monocytes/macrophages, T cells,
NK cells, and dendritic cells. (ii) Promotes T cell adhesion to endothelial
cells. (iii) Antitumor activity (iv) Inhibition of bone marrow colony
formation (v) Angiogenesis (Luster et al., 1985; Dufour et al., 2002;
Angiolillo et al., 1995).

  * 5.

mtv: Metabolic tumor volume from PET imaging.

  * 6.

GLSZM_LZLGE: Radiomics features: the large zone low gray-level emphasis
(LZLGE) feature of a gray-level size zone matrix defined as
$\sum_{i=1^{N_{g}}}\sum_{j=1}^{L_{z}}\frac{j^{2}p(i,j)}{i^{2}}$ (Carrier-
Vallières, 2018).

  * 7.

GLSZM_ZSV: Radiomics features: the zone-size variance (ZSV) feature of a gray-
level size zone matrix defined as
$\frac{1}{N_{g}L_{z}}\sum_{i=1^{N_{g}}}\sum_{j=1}^{L_{Z}}(jp(i,j)-\mu_{j})^{2}$
(Carrier-Vallières, 2018).

  * 8.

Lung gEUD ($\alpha/\beta=4Gy,10Gy\,resp$): Generalized equivalent uniform dose
of lung converted from EQD2 dose distributions: $EQD=N_{frac}\times
d\times\left(\frac{d+\alpha/\beta}{2+\alpha/\beta}\right)$ (Bentzen and
Dische, 2000; Søvik et al., 2007).

  * 9.

Tumor gEUD ($\alpha/\beta=4Gy,10Gy\,resp$): Generalized equivalent uniform
dose of tumor converted from EQD2 dose distributions: $EQD=N_{frac}\times
d\times\left(\frac{d+\alpha/\beta}{2+\alpha/\beta}\right)$ (Bentzen and
Dische, 2000; Søvik et al., 2007).

  * 10.

Rs2234671: A SNP in the gene cxcr1 also known as Interleukin 8 receptor, alpha
(IL8RA) relation to radiation induced toxicity in non-small cell lung cancer
(Hildebrandt et al., 2010).

  * 11.

Rs238406: A SNP in the gene ercc2 known to repair DNA excision and related to
risk of lung cancer (Chang et al., 2008). SNP: Single nucleotide polymorphism
is a substitution of single nucleotide that occurs at a specific position in
the genome.

  * 12.

Rs1047768: A SNP in the gene ercc5 also known to repair DNA excision and
related to lung cancer susceptibility (Kiyohara and Yoshimasu, 2007).

## References

  * Angiolillo et al. (1995) Angiolillo, A.L., Sgadari, C., Taub, D.D., Liao, F., Farber, J.M., Maheshwari, S., Kleinman, H.K., Reaman, G.H., Tosato, G., 1995.  Human interferon-inducible protein 10 is a potent inhibitor of angiogenesis in vivo.  The Journal of experimental medicine 182, 155–162.
  * Ashton et al. (2018) Ashton, J.R., Castle, K.D., Qi, Y., Kirsch, D.G., West, J.L., Badea, C.T., 2018\.  Dual-energy ct imaging of tumor liposome delivery after gold nanoparticle-augmented radiation therapy.  Theranostics 8, 1782\.
  * Benedict et al. (2016) Benedict, S.H., Hoffman, K., Martel, M.K., Abernethy, A.P., Asher, A.L., Capala, J., Chen, R.C., Chera, B., Couch, J., Deye, J., et al., 2016\.  Overview of the american society for radiation oncology–national institutes of health–american association of physicists in medicine workshop 2015: Exploring opportunities for radiation oncology in the era of big data.  International Journal of Radiation Oncology• Biology• Physics 95, 873–879.
  * Bentzen and Dische (2000) Bentzen, S.M., Dische, S., 2000\.  Morbidity related to axillary irradiation in the treatment of breast cancer.  Acta oncologica 39, 337–347.
  * Carrier-Vallières (2018) Carrier-Vallières, M., 2018.  Radiomics: Enabling Factors Towards Precision Medicine.  Ph.D. thesis. McGill University Libraries.
  * Chakraborty (2013) Chakraborty, B., 2013.  Statistical methods for dynamic treatment regimes.  Springer.
  * Chang et al. (2008) Chang, J.S., Wrensch, M.R., Hansen, H.M., Sison, J.D., Aldrich, M.C., Quesenberry Jr, C.P., Seldin, M.F., Kelsey, K.T., Kittles, R.A., Silva, G., et al., 2008\.  Nucleotide excision repair genes and risk of lung cancer among san francisco bay area latinos and african americans.  International journal of cancer 123, 2095–2104.
  * Cressie and Kornak (2003) Cressie, N., Kornak, J., 2003\.  Spatial statistics in the presence of location error with an application to remote sensing of the environment.  Statistical science , 436–456.
  * Dufour et al. (2002) Dufour, J.H., Dziejman, M., Liu, M.T., Leung, J.H., Lane, T.E., Luster, A.D., 2002\.  Ifn-$\gamma$-inducible protein 10 (ip-10; cxcl10)-deficient mice reveal a role for ip-10 in effector t cell generation and trafficking.  The Journal of Immunology 168, 3195–3204.
  * El Naqa et al. (2018) El Naqa, I., Kosorok, M.R., Jin, J., Mierzwa, M., Ten Haken, R.K., 2018.  Prospects and challenges for clinical decision support in the era of big data.  JCO clinical cancer informatics 2, 1–12.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., Courville, A., 2016.  Deep learning.  MIT press.
  * Goodfellow et al. (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y., 2014\.  Generative adversarial nets, in: Advances in neural information processing systems, pp. 2672–2680.
  * Hildebrandt et al. (2010) Hildebrandt, M.A., Komaki, R., Liao, Z., Gu, J., Chang, J.Y., Ye, Y., Lu, C., Stewart, D.J., Minna, J.D., Roth, J.A., et al., 2010.  Genetic variants in inflammation-related genes are associated with radiation-induced toxicity following treatment for non-small cell lung cancer.  PloS one 5.
  * Kennedy and O’Hagan (2001) Kennedy, M.C., O’Hagan, A., 2001\.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63, 425–464.
  * Kiyohara and Yoshimasu (2007) Kiyohara, C., Yoshimasu, K., 2007\.  Genetic polymorphisms in the nucleotide excision repair pathway and lung cancer risk: a meta-analysis.  International journal of medical sciences 4, 59.
  * Kong et al. (2017) Kong, F.M., Ten Haken, R.K., Schipper, M., Frey, K.A., Hayman, J., Gross, M., Ramnath, N., Hassan, K.A., Matuszak, M., Ritter, T., et al., 2017\.  Effect of midtreatment pet/ct-adapted radiation therapy with concurrent chemotherapy in patients with locally advanced non–small-cell lung cancer: a phase 2 clinical trial.  JAMA oncology 3, 1358–1365.
  * Luo et al. (2017) Luo, Y., El Naqa, I., McShan, D.L., Ray, D., Lohse, I., Matuszak, M.M., Owen, D., Jolly, S., Lawrence, T.S., Ten Haken, R.K., et al., 2017\.  Unraveling biophysical interactions of radiation pneumonitis in non-small-cell lung cancer via bayesian network analysis.  Radiotherapy and Oncology 123, 85–92.
  * Luo et al. (2018) Luo, Y., McShan, D.L., Matuszak, M.M., Ray, D., Lawrence, T.S., Jolly, S., Kong, F.M., Ten Haken, R.K., El Naqa, I., 2018.  A multiobjective bayesian networks approach for joint prediction of tumor local control and radiation pneumonitis in nonsmall-cell lung cancer (nsclc) for response-adapted radiotherapy.  Medical physics 45, 3980–3995.
  * Luster et al. (1985) Luster, A.D., Unkeless, J.C., Ravetch, J.V., 1985.  $\gamma$-interferon transcriptionally regulates an early-response gene containing homology to platelet proteins.  Nature 315, 672–676.
  * Moodie et al. (2014) Moodie, E.E., Dean, N., Sun, Y.R., 2014.  Q-learning: Flexible learning about useful utilities.  Statistics in Biosciences 6, 223–243.
  * Pearl et al. (2009) Pearl, J., et al., 2009.  Causal inference in statistics: An overview.  Statistics surveys 3, 96–146.
  * Qian and Murphy (2011) Qian, M., Murphy, S.A., 2011\.  Performance guarantees for individualized treatment rules.  Annals of statistics 39, 1180\.
  * Ramírez et al. (2013) Ramírez, M.F., Huitink, J.M., Cata, J.P., 2013.  Perioperative clinical interventions that modify the immune response in cancer patients .
  * Rasmussen (2003) Rasmussen, C.E., 2003.  Gaussian processes in machine learning, in: Summer School on Machine Learning, Springer. pp. 63–71.
  * Rich and Gureckis (2019) Rich, A.S., Gureckis, T.M., 2019\.  Lessons for artificial intelligence from the study of natural stupidity.  Nature Machine Intelligence 1, 174–180.
  * Schaue et al. (2012) Schaue, D., Kachikwu, E.L., McBride, W.H., 2012.  Cytokines in radiobiological responses: a review.  Radiation research 178, 505–523.
  * Søvik et al. (2007) Søvik, Å., Ovrum, J., Olsen, D.R., Malinen, E., 2007\.  On the parameter describing the generalised equivalent uniform dose (geud) for tumours.  Physica Medica 23, 100–106.
  * Tierney and Kadane (1986) Tierney, L., Kadane, J.B., 1986\.  Accurate approximations for posterior moments and marginal densities.  Journal of the american statistical association 81, 82–86.
  * Tseng et al. (2017) Tseng, H.H., Luo, Y., Cui, S., Chien, J.T., Ten Haken, R.K., El Naqa, I., 2017.  Deep reinforcement learning for automated radiation adaptation in lung cancer.  Medical physics 44, 6690–6705.
  * Warltier et al. (2002) Warltier, D.C., Laffey, J.G., Boylan, J.F., Cheng, D.C., 2002\.  The systemic inflammatory response to cardiac surgeryimplications for the anesthesiologist.  Anesthesiology: The Journal of the American Society of Anesthesiologists 97, 215–252.
  * Zhao et al. (2012) Zhao, Y., Zeng, D., Rush, A.J., Kosorok, M.R., 2012.  Estimating individualized treatment rules using outcome weighted learning.  Journal of the American Statistical Association 107, 1106–1118.


# Learning the helix topology of musical pitch

###### Abstract

To explain the consonance of octaves, music psychologists represent pitch as a
helix where azimuth and axial coordinate correspond to pitch class and pitch
height respectively. This article addresses the problem of discovering this
helical structure from unlabeled audio data. We measure Pearson correlations
in the constant-$Q$ transform (CQT) domain to build a $K$-nearest neighbor
graph between frequency subbands. Then, we run the Isomap manifold learning
algorithm to represent this graph in a three-dimensional space in which
straight lines approximate graph geodesics. Experiments on isolated musical
notes demonstrate that the resulting manifold resembles a helix which makes a
full turn at every octave. A circular shape is also found in English speech,
but not in urban noise. We discuss the impact of various design choices on the
visualization: instrumentarium, loudness mapping function, and number of
neighbors $K$.

Index Terms—  Continuous wavelet transforms, distance learning, music, pitch
control (audio), shortest path problem.

##  1 Introduction

Listening to a sequence of two pure tones elicits a sensation of pitch going
“up” or “down”, correlating with changes in fundamental frequency ($f_{0}$).
However, contrary to pure tones, natural pitched sounds contain a rich
spectrum of components in addition to $f_{0}$. Neglecting inharmonicity, these
components are tuned to an ideal Fourier series whose modes resonate at
integer multiples of the fundamental: $2f_{0}$, $3f_{0}$, and so forth. By the
change of variable $f_{0}^{\prime}=2f_{0}$, it appears that all even-numbered
partials $2f_{0}=f_{0}^{\prime}$, $4f_{0}=2f_{0}^{\prime}$,
$6f_{0}=3f_{0}^{\prime}$, and so forth make up the whole Fourier series of a
periodic signal whose fundamental frequency is $f_{0}^{\prime}$. Figure 1
illustrates, in the case of a synthetic signal with perfect harmonicity, that
filtering out all odd-numbered partials $(2p+1)f_{0}$ for integer $p\geq 0$
results in a perceived pitch that morphs from $f_{0}$ to $2f_{0}$, i.e., up
one octave [1].

Such ambivalence brings about well-known auditory paradoxes: tones that have a
definite pitch class but lack a pitch register [2]; glissandi that seem to
ascend or descend endlessly [3]; and tritone intervals whose pitch
directionality depends on prior context [4]. To explain them, one may roll up
the frequency axis onto a spiral or helix which makes a full turn at each
octave, thereby aligning power-of-two harmonics onto the same radii.

The consonance of octaves has implications in several disciplines. In music
theory, it allows for pitches to be grouped into pitch classes: e.g., in
European solfège, _do re mi_ etc. eventually “circle back” to _do_. In
ethnomusicology, it transcends the boundaries between cultures, to the point
of being held as a nearly universal attribute of music [5]. In
neurophysiology, it explains the functional organization of the central
auditory cortex [6]. In music information research, it motivates the design of
chroma features, i.e., a representation of harmonic content that is meant to
be equivariant to parallel chord progressions but invariant to chord
inversions [7, chapter 5].

Despite the wealth of evidence for the crucial role of octaves in music, there
is, to this day, no data-driven criterion for assessing whether a given audio
corpus exhibits a property of octave equivalence. Rather, the disentanglement
of pitch chroma and pitch height in time–frequency representations relies on
domain-specific knowledge about music [8]. More generally, although the
induction of priors on the topology of pitch is widespread in symbolic music
analysis [9], few of them directly apply to audio signal processing [10].

Yet, in recent years, the systematic use of machine learning methods has
progressively reduced the need for domain-specific knowledge in several other
aspects of auditory perception, including mel-frequency spectrum [11] and
adaptive gain control [12]. It remains to be known whether octave equivalence
can, in turn, be discovered by a machine learning algorithm, instead of being
engineered ad hoc. Furthermore, it is unclear whether octave equivalence is an
exclusive characteristic of music or whether it may extend to other kinds of
sounds, such as speech or environmental soundscapes.

Fig. 1: Two continuous trajectories in pitch space, either by octave glissando
(left) or by attenuation of odd-numbered partials (right). Darker shades
indicate larger magnitudes of the constant-$Q$ transform. Axes $t$ and
$\gamma$ respectively correspond to time and log-frequency. Vertical ticks
denote octaves. See Section 1 for details.

In this article, we conduct an unsupervised manifold learning experiment,
inspired by the protocol of “learning the 2-D topology of images” [13], to
visualize the helix topology of constant-$Q$ spectra. Starting from a dataset
of isolated notes from various musical instruments, we run the Isomap
algorithm to represent each of these frequency subbands as a dot in a 3-D
space, wherein spatial neighborhoods denote high correlation in loudness.
Contrary to natural images, we find a mismatch between the physical
dimensionality of natural acoustic spectra (i.e., 1-D) and their statistical
dimensionality (i.e., 3-D or greater).

The companion website of this paper111Companion website:
https://github.com/BirdVox/lostanlen2020icassp contains a Python package to
visualize octave equivalence in audio data.

Fig. 2: Functional diagram of the proposed method, comprising: constant-$Q$
transform $\mathbf{W}$, data matrix $\mathbf{X}$, extraction of Pearson
correlations $\boldsymbol{\rho}^{2}$, shortest path distance matrix
$\mathbf{D}_{\mathcal{G}}$, and Isomap eigenbasis $\mathbf{e}$ for the set of
vertices $\mathcal{V}$. Darker shades in $\boldsymbol{\rho}^{2}$ and
$\mathbf{D}_{\mathcal{G}}$ indicate larger absolute values of the Pearson
correlation and distance respectively. The hue of colored dots and the solid
grey line respectively denote pitch chroma and pitch height. The first three
dimensions in the Isomap embedding explain 36%, 35%, and 9% of the total
variance in $\mathbf{D}_{\mathcal{G}}$ respectively. Note that our method is
unsupervised: neither pitch chroma nor pitch height are directly supplied to
the Isomap manifold learning algorithm. See Section 2 for details.

##  2 Isomap embedding of subband correlations

###  2.1 Constant-$Q$ transform and loudness mapping

Given a corpus of audio signals
$\boldsymbol{x}_{1}\ldots\boldsymbol{x}_{N}\in\mathbf{L}^{2}(\mathbb{R})$ in
arbitrary order, we define their constant-$Q$ transforms (CQT) as the
convolutions with wavelets $t\mapsto
2^{\gamma}\boldsymbol{\psi}(2^{\gamma}t)$, where $2^{-\gamma}$ is a scale
parameter:

$\mathbf{W}\boldsymbol{x}_{n}:(t,\gamma)\mapsto
2^{\gamma}\int_{-\infty}^{+\infty}\boldsymbol{x}_{n}(t^{\prime})\boldsymbol{\psi}(2^{\gamma}(t-t^{\prime}))\,\mathrm{d}t.$
(1)

The wavelet filterbank in the operator $\mathbf{W}$ covers a range of $J$
octaves. In the case of musical notes, we define a region of interest $t_{n}$
in each $\boldsymbol{x}_{n}$ as the frame of highest short-term energy at a
scale of $93\text{\,}\mathrm{ms}$. We compute their scalogram representation
as the vector of CQT modulus responses at some $t=t_{n}$, for discretized
values $u=1+Q\gamma$:

$\mathbf{X}[n,u]=\left|\mathbf{W}\boldsymbol{x}_{n}\left(t_{n},\frac{u-1}{Q}\right)\right|,$
(2)

where the integer $u$ ranges from $1$ to $QJ$. Then, we apply a pointwise
logarithmic compression $\Lambda$ to map each magnitude coefficient
$\mathbf{X}[n,u]$ onto a decibel-like scale, and clip it to
$-100\text{\,}\mathrm{dB}$:

$\Lambda(\mathbf{X})[n,u]=\max(-100,10\log_{10}\mathbf{\mathbf{X}}[n,u]).$ (3)

Unless stated otherwise, we set $|\boldsymbol{\psi}|$ to a Hann window, the
quality factor $Q$ to $24$, and the number of octaves $J$ to 3 in the
following. We compute constant-$Q$ transforms with librosa v0.6.1 [14].
Section 3.4 will discuss the effect of alternative choices for $\Lambda$.

###  2.2 Pearson autocorrelation of log-magnitude spectra

In accordance with [13], we express the similarity between two features $u$
and $v$ in terms of their squared Pearson correlation
$\boldsymbol{\rho}^{2}[u,v]$. We begin by recentering each feature to null
mean, yielding the matrix

$\mathbf{Y}[n,u]=\Lambda(\mathbf{X})[n,u]-\dfrac{1}{N}\sum_{n=1}^{N}\Lambda(\mathbf{X})[n,u],$
(4)

and then compute squared cosine similarities on all pairs $(u,v)$:

$\boldsymbol{\rho}^{2}[u,v]=\dfrac{(\sum_{n=1}^{N}\mathbf{Y}[n,u]\mathbf{Y}[n,v])^{2}}{(\sum_{n=1}^{N}\mathbf{Y}^{2}[n,u])\times(\sum_{n=1}^{N}\mathbf{Y}^{2}[n,v])}.$
(5)

Let $\mathcal{V}$ be the set of all $QJ$ features $u$ in $\mathbf{Y}$.
Adopting a manifold learning perspective, we may regard
$\boldsymbol{\rho}^{2}$ as the values of a standard Gaussian kernel
$\kappa:\mathcal{V}\times\mathcal{V}\longrightarrow\mathbb{R}$. Inverting the
identity
$\boldsymbol{\rho}^{2}[u,v]=\kappa(u,v)=\exp(-2\mathbf{D}_{\boldsymbol{\rho^{2}}}^{2}[u,v])$
leads to a pseudo-Euclidean distance

$\mathbf{D}_{\boldsymbol{\rho}^{2}}[u,v]=\sqrt{-\frac{1}{2}\log\boldsymbol{\rho}^{2}[u,v]}.$
(6)

###  2.3 Shortest path distance on the $K$-nearest neighbor graph

Following the Isomap manifold learning algorithm [15], we compute the $K$
nearest neighbors of each feature $u\in\mathcal{V}$ as the set
$\mathcal{N}_{K}(u)$ of $K$ features $v\neq u$ that minimize the distance
$\mathbf{D}_{\boldsymbol{\rho^{2}}}[u,v]$, i.e., maximize the squared Pearson
correlation $\boldsymbol{\rho}^{2}[u,v]$. We construct a $K$-nearest graph
$\mathcal{G}$ whose vertices are $\mathcal{V}$ and whose adjacency matrix
$\mathbf{A}[u,v]$ is equal to $\mathbf{D}_{\boldsymbol{\rho}}[u,v]$ if
$v\in\mathcal{N}_{K}(u)$ or if $u\in\mathcal{N}_{K}(v)$, and infinity
otherwise. We run Dijkstra’s algorithm on $\mathcal{G}$ to measure geodesics
on the manifold induced by $\mathbf{A}$. These geodesics yield a shortest path
distance function over $\mathcal{V}^{2}$:

$\mathbf{D}_{\mathcal{G}}[u,v]=\left\\{\begin{array}[]{c}\mathbf{D}_{\boldsymbol{\rho}^{2}}[u,v]\text{
if }v\in\mathcal{N}_{K}(u)\text{ or }u\in\mathcal{N}_{K}(v)\\\
\displaystyle{\min_{z}}\;(\mathbf{D}_{\mathcal{G}}[u,z]+\mathbf{D}_{\mathcal{G}}[z,v])\text{
otherwise.}\end{array}\right.$ (7)

If $\mathcal{G}$ has more than one connected component, then
$\mathbf{D}_{\mathcal{G}}$ is infinite over pairs of mutually unreachable
vertices. On some datasets, $\boldsymbol{\rho}$ may lead to a disconnected
$K$-nearest neighbor graph $\mathcal{G}$, especially for small $K$, thereby
causing numerical aberrations in Isomap. However, we make sure that the
effective bandwidth of the wavelet $\widehat{\boldsymbol{\psi}}$ is large
enough in comparison with $Q$ to yield strong correlations
$\boldsymbol{\rho}[u,u+1]$ for all $u$, and thus $u\in\mathcal{N}_{K}(u+1)$.
With this caveat in mind, we postulate that $\mathbf{D}_{\mathcal{G}}$ is
finite in all of the following. We set $K=3$ in the following unless stated
otherwise.

###  2.4 Classical multidimensional scaling and 3-D embedding

Let $\mathbf{S}_{\mathcal{G}}[u,v]=\mathbf{D}_{\mathcal{G}}[u,v]^{2}$.
Classical multidimensional scaling (MDS) diagonalizes
$\boldsymbol{\tau}_{\mathcal{G}}=-\frac{1}{2}(\mathbf{H}\mathbf{S}_{\mathcal{G}}\mathbf{H})$
where $\mathbf{H}[u,v]=\boldsymbol{\delta}(u-v)-\frac{1}{QJ}$ [16]. We denote
by $\mathbf{e}_{m}$ and $\lambda_{m}$ the respective eigenvectors and
eigenvalues of $\boldsymbol{\tau}_{\mathcal{G}}$, satisfying
$\boldsymbol{\tau}_{\mathcal{G}}\mathbf{e}_{m}=\lambda_{m}\mathbf{e}_{m}$. We
rank eigenvalues in decreasing order, without loss of generality. Lastly, we
display the Isomap embedding as a 3-D scatter plot with Cartesian coordinates
$\mathbf{e}[u]=(\mathbf{e_{1}}[u],\mathbf{e_{2}}[u],\mathbf{e_{3}}[u])$ for
every vertex $u\in\mathcal{V}$. We compute Isomap using scikit-learn v0.21.3
[17].

(a) Trumpet in C only.

(b) Harp only.

Fig. 3:  Pearson correlation matrices $\boldsymbol{\rho}^{2}$ (left) and
Isomap embeddings (right) for two instruments in the SOL dataset: Trumpet in C
(a) and Harp (b). In both cases, the number of neighbors is set to $K=3$ and
the loudness mapping $\Lambda$ is logarithmic. See Section 3.2 for details.

##  3 Experiments with musical sounds

###  3.1 Main protocol

We extract $\mathbf{X}$ from the isolated musical notes played by eight
instruments in the SOL dataset [18]: accordion, alto saxophone, bassoon,
flute, harp, trumpet in C, and cello. For each of these instruments, we
include three levels of intensity dynamics: _pp_ , _mf_ , and _ff_. We include
all available pitches in the tessitura ($\min=\mathrm{B_{0}}$,
$\operatorname{median}=\textrm{E}_{4}$,
$\max=\textrm{C}\mbox{\musFontLarge\kern 0.90005pt\raisebox{1.93748pt}{\char
92\relax}\kern-0.1pt}_{8}$). We exclude extended playing techniques, because
some of them may lack a discernible pitch class [19]; thus resulting in a
total of $N=1212$ audio recordings for the _ordinario_ technique. We
hypothesize that cross-correlations across octaves are weaker than cross-
correlations along the log-frequency axis. Because the helix topology of
musical pitch relies on both kinds of correlation, we set the number of
neighbors to $K=3$.

Figure 2 illustrates our protocol and main finding. Isomap produces a quasi-
perfect cylindrical manifold in dimension three. We color each dot
$\mathbf{e}[u]$ according to a hue of
$\boldsymbol{\theta}(u)=\frac{2\pi}{Q}(u\mod Q)$ radians. Furthermore, we draw
a segment between each $\mathbf{e}[u]$ and its upper adjacent subband
$\mathbf{e}[u+1]$. Once these visual elements are included in the display of
the scatter plot $(\mathbf{e}[u])_{u}$, the cylindrical manifold appears to
coincide with the Drobisch-Shepard helix in music psychology [20, 21]. Indeed,
hues $\boldsymbol{\theta}[u]$ appear to align on the same radii, whereas the
grey line $\gamma$ grows monotonically with $\mathbf{e_{3}}$. This result
demonstrates that, even without prior knowledge about the perception of
musical pitch, it is possible to discover octave equivalence in a purely data-
driven fashion, by computing the graph of greatest cross-correlations between
CQT magnitudes in a corpus of isolated musical notes.

###  3.2 Varying the instrumentarium

We reproduce the main protocol on subsets of the SOL dataset, involving a
single instrument at once. The bright timbre of brass instruments correlates
with relatively loud high-order partials [22], resulting in large octave
equivalence and a helical topology (see Figure 3 (a)). In contrast, harp tones
carry little energy at twice the fundamental; yet, they induce sympathetic
resonance along the soundboard, which affects nearby strings predominantly
[23]. These two phenomena in combination favor semitone correlations over
octave correlations, resulting in a rectilinear topology (see Figure 3 (a)).

###  3.3 Varying the number of neighbors $K$

We reproduce the main protocol with varying $K$-nearest neighbor graphs.
Setting $K=2$ results in multiple lobes, each corresponding to an octave, and
connected at a single pitch class (see Figure 4). This confirms our hypothesis
that large semitone correlations outnumber large octave correlations.
Conversely, setting $K=5$ results in a topology that is more intricate than
the Drobisch-Shepard helix, involving correlations across perfect fourths and
fifths..

(a) $K=2$ neighbors. (b) $K=5$ neighbors.

Fig. 4:  Isomap embeddings for varying $K$-nearest neighbor graphs: $K=2$
(left) and $K=5$ (right). In both cases, the Pearson correlation matrix
$\boldsymbol{\rho}^{2}$ results from all instruments in the SOL dataset and
the loudness mapping $\Lambda$ is logarithmic. See Section 3.3 for details.

(a) Linear loudness mapping function.

(b) Cubic root loudness mapping function.

Fig. 5:  Pearson correlation matrices $\boldsymbol{\rho}^{2}$ (left) and
Isomap embeddings (right) for cubic root (a) and linear (b) loudness mappings.
In both cases, $\boldsymbol{\rho}^{2}$ results from all instruments in the SOL
dataset and the number of neighbors is set to $K=3$. See Section 3.4 for
details.

###  3.4 Varying the loudness mapping function $\Lambda$

We reproduce the main protocol with varying loudness mappings $\Lambda$ (see
Figure 5). On one hand, setting $\Lambda$ to the identity function yields a
$K$-nearest neighbor graph in which octave correlations are numerically
negligible, except in the lower register. This results in an Isomap embedding
which is circular in the bottommost octave and irregular in the topmost
octaves. On the other hand, setting $\Lambda$ to the cubic root function
yields a helical topology. This experiment demonstrates the need for nonlinear
loudness compression in the protocol; in contrast with [13], which relied on
raw grayscale intensities in handwritten digits.

##  4 Extension to speech and urban sounds

###  4.1 Experiment with speech data

We analyze the North Texas vowel database (NTVOW), which contains utterances
of 12 English vowels from 50 American speakers, including children aged three
to seven as well as male and female adults [24]; resulting in a total of
$N=9570$ audio recordings. As seen in Figure 6 (a), which includes the data
from all age groups, there is some notion of pitch circularity seen in the
mid-frequency range, but not so in the low and high-frequency ranges. This is
because the distribution of $f_{0}$ in human speech is polarized around
certain pitch classes, rather than uniform over the chromatic scale.

###  4.2 Experiment with environmental audio data

We analyze a portion of the SONYC Urban Sound Tagging dataset (SONYC-UST
v0.4), which contains a collection of 3068 acoustic scenes from a network of
autonomous sensors in various locations of New York City [25]. We restrict
this collection to the acoustic scenes in which the consensus of expert
annotators has confirmed the absence of both human speech and music. As a
result of this preliminary curation step, we obtain $N=233$ audio recordings
from eight different sensor locations. Each of these scenes contains one or
several sources of urban noise pollution, among which: engines, machinery and
non-machinery impacts, powered saws, alert signals, and dog barks. Figure 6
(b) shows that no discernible correlations across octaves are observed in this
dataset. This finding confirms the conclusions of a previous publication [26],
which stated that “music audio signal processing techniques must be informed
by a deep and thorough insight into the nature of music itself”.

(a) Spoken English vowels (NTVOW dataset).

(b) Urban sounds (SONYC-UST excluding speech and music).

Fig. 6: Extension of the proposed method to speech (a) and environmental
soundscapes (b). In both cases, the number of neighbors is set to $K=3$ and
the loudness mapping $\Lambda$ is logarithmic. See Section 4 for details.

##  5 Conclusion

The Isomap manifold learning algorithm offers an approximate visualization of
nearest neighbor relationships in any non-Euclidean metric space. Thus, a
previous publication [13] proposed to apply Isomap onto cross-correlations
between grayscale intensities in natural images. In this article, we have
borrowed from the protocol of this publication to apply it on music data.

Despite their methodological resemblance, the two studies lead to different
insights. While [13] recovered a quasi-uniform raster, we do not recover a
straight line from cross-correlations along the log-frequency axis. Instead,
we obtain a cylindrical lattice in dimension three. Assigning pitch classes to
subbands reveals that this lattice is akin to a Drobisch-Shepard helix [20],
which makes a full turn at each octave. Thus, whereas [13] learned a 2-D
topology from 2-D data, we learned a 3-D topology from 1-D data. Furthermore,
after benchmarking several design choices, we deduce that the most regular
helical shape results from: a diverse instrumentarium; a graph of $K=3$
nearest neighbors; and a logarithmic loudness mapping. Lastly, we have
discussed the limitations of our findings: although spoken vowels also exhibit
a quasi-helical topology in subband neighborhoods, the same cannot be said of
urban noise.

Beyond the realm of manifold learning, the present article motivates the
development of weight sharing architectures that foster octave equivalence in
deep representations of music data. Three examples of such architectures are:
spiral scattering transform [27]; spiral convolutional networks [28]; and
harmonic constant-$Q$ transform [29]. Future work will extend this protocol to
bioacoustic data, and comparing the influence of species-specific
vocalizations onto the empirical topology of the frequency domain.

##  6 Acknowledgment

The authors wish to thank Stéphane Mallat for helpful discussions.

## References

  * [1] Diana Deutsch, Kevin Dooley, and Trevor Henthorn,  “Pitch circularity from tones comprising full harmonic series,”  Journal of the Acoustical Society of America, vol. 124, no. 1, pp. 589–597, 2008.
  * [2] Diana Deutsch,  “The paradox of pitch circularity,”  Acoustics Today, vol. 7, pp. 8–15, 2010.
  * [3] Jean-Claude Risset,  “Pitch control and pitch paradoxes demonstrated with computer-synthesized sounds,”  Journal of the Acoustical Society of America, vol. 46, no. 1A, pp. 88–88, 1969.
  * [4] Claire Pelofi, Vincent de Gardelle, Paul Egré, and Daniel Pressnitzer,  “Interindividual variability in auditory scene analysis revealed by confidence judgements,”  Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 372, no. 1714, pp. 20160107, 2017.
  * [5] Edward M. Burns,  “Intervals, scales, and tuning,”  in The Psychology of Music, pp. 215–264. Elsevier, 1999.
  * [6] Jason D. Warren, Stefan Uppenkamp, Roy D. Patterson, and Timothy D. Griffiths,  “Separating pitch chroma and pitch height in the human brain,”  Proceedings of the National Academy of Sciences, vol. 100, no. 17, pp. 10038–10042, 2003.
  * [7] Meinard Müller,  Fundamentals of music processing: Audio, analysis, algorithms, applications,  Springer, 2015.
  * [8] Dmitri Tymoczko,  “The geometry of musical chords,”  Science, vol. 313, no. 5783, pp. 72–74, 2006.
  * [9] Louis Bigo and Moreno Andreatta,  “Topological structures in computer-aided music analysis,”  in Computational Music Analysis, David Meredith, Ed., pp. 57–80. Springer International Publishing, Cham, 2016.
  * [10] Ching-Hua Chuan and Elaine Chew,  “Polyphonic audio key finding using the spiral array ceg algorithm,”  in 2005 IEEE International Conference on Multimedia and Expo. IEEE, 2005, pp. 21–24.
  * [11] Neil Zeghidour, Nicolas Usunier, Iasonas Kokkinos, Thomas Schaiz, Gabriel Synnaeve, and Emmanuel Dupoux,  “Learning filterbanks from raw speech for phone recognition,”  in Procedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 5509–5513.
  * [12] Yuxuan Wang, Pascal Getreuer, Thad Hughes, Richard F. Lyon, and Rif A. Saurous,  “Trainable frontend for robust and far-field keyword spotting,”  in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2017, pp. 5670–5674.
  * [13] Nicolas Le Roux, Yoshua Bengio, Pascal Lamblin, Marc Joliveau, and Balázs Kégl,  “Learning the 2-D topology of images,”  in Advances in Neural Information Processing Systems (NeurIPS), 2008, pp. 841–848.
  * [14] Brian McFee, Matt McVicar, Stefan Balke, Carl Thomé, Colin Raffel, Dana Lee, Oriol Nieto, Eric Battenberg, Dan Ellis, Ryuichi Yamamoto, Josh Moore, Rachel Bittner, Keunwoo Choi, Pius Friesch, Fabian-Robert Stöter, Vincent Lostanlen, Siddhartha Kumar, Simon Waloschek, Seth Kranzler, Rimvydas Naktinis, Douglas Repetto, Curtis “Fjord” Hawthorne, CJ Carr, Waldir Pimenta, Petr Viktorin, Paul Brossier, Joao Felipe Santos, Jackie Wu, Erik Peterson, and Adrian Holovaty,  “librosa/librosa: 0.6.1,” May 2018.
  * [15] Joshua B. Tenenbaum, Vin De Silva, and John C. Langford,  “A global geometric framework for nonlinear dimensionality reduction,”  Science, vol. 290, no. 5500, pp. 2319–2323, 2000.
  * [16] Warren S. Torgerson,  “Multidimensional scaling: I. theory and method,”  Psychometrika, vol. 17, no. 4, pp. 401–419, 1952.
  * [17] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay,  “Scikit-learn: Machine learning in Python,”  Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011\.
  * [18] Guillaume Ballet, Riccardo Borghesi, Peter Hoffmann, and Fabien Lévy,  “Studio Online 3.0: An Internet “killer application” for remote access to Ircam sounds and processing tools,”  in Actes des Journées d’Informatiques Musicales (JIM), 1999\.
  * [19] V. Lostanlen, J. Andén, and M. Lagrange,  “Extended playing techniques$:$ the next milestone in musical instrument recognition,”  in Proceedings of the International Conference on Digital Libraries for Musicology (DLfM), 2018.
  * [20] Roger N. Shepard,  “Circularity in judgments of relative pitch,”  Journal of the Acoustical Society of America, vol. 36, no. 12, pp. 2346–2353, 1964.
  * [21] Fred Lerdahl,  Tonal pitch space,  Oxford University Press, 2004.
  * [22] Emilie Poirson, Jean-François Petiot, and Joël Gilbert,  “Study of the brightness of trumpet tones,”  Journal of the Acoustical Society of America, vol. 118, no. 4, pp. 2656–2666, 2005.
  * [23] Jean-Loïc Le Carrou, François Gautier, and Roland Badeau,  “Sympathetic string modes in the concert harp,”  Acta Acustica united with Acustica, vol. 95, no. 4, pp. 744–752, 2009.
  * [24] Peter F. Assmann and William F. Katz,  “Time-varying spectral change in the vowels of children and adults,”  Journal of the Acoustical Society of America, vol. 108, no. 4, pp. 1856–1866, 2000.
  * [25] Mark Cartwright, Ana Elisa Mendez Mendez, Graham Dove, Jason Cramer, Vincent Lostanlen, Ho-Hsiang Wu, Justin Salamon, Oded Nov, and Juan Pablo Bello,  “SONYC Urban Sound Tagging (SONYC-UST): A multilabel dataset from an urban acoustic sensor network,”  in Proceedings of the International Workshop on Detection and Classification of Acoustic Scenes and Events (DCASE), 2019.
  * [26] Meinard Müller, Daniel P. W. Ellis, Anssi Klapuri, and Gaël Richard,  “Signal processing for music analysis,”  IEEE Journal of Selected Topics in Signal Processing, vol. 5, no. 6, pp. 1088–1110, 2011.
  * [27] Vincent Lostanlen and Stéphane Mallat,  “Wavelet scattering on the pitch spiral,”  in Proceedings of the International Conference on Digital Audio Effects (DAFX), 2016, pp. 429–432.
  * [28] Vincent Lostanlen and Carmine Emanuele Cella,  “Deep convolutional networks on the pitch spiral for musical instrument recognition,”  in Proceedings of the International Society on Music Information Retrieval (ISMIR) Conference, 2016.
  * [29] Rachel Bittner, Brian McFee, Justin Salamon, Peter Li, and Juan Pablo Bello,  “Deep salience representations for $f_{0}$ estimation in polyphonic music,”  in Proceedings of the International Society on Music Information Retrieval (ISMIR) Conference, 2017.