# Resonant Dynamical Casimir Effect in Waveguide-coupled Qubit Arrays

Egor S. Vyatkin  Ioffe Institute, St. Petersburg 194021, Russia  Alexander V.
Poshakinskiy  Ioffe Institute, St. Petersburg 194021, Russia  Alexander N.
Poddubny  Department of Physics of Complex Systems, Weizmann Institute of
Science, Rehovot 7610001, Israel<EMAIL_ADDRESS>
###### Abstract

We have developed a theory of parametric photon generation in the waveguides
coupled to arrays of qubits with temporally modulated resonance frequencies.
Such generation can be interpreted as a dynamical Casimir effect. We
demonstrate numerically and analytically how the emission directionality and
photon-photon correlations can be controlled by the phases of the modulation.
The emission spectrum is shown to be strongly dependent on the qubit
anharmonicity parameter. Single- and double-excited state resonances have been
identified in the emission spectrum.

##  I Introduction

Waveguide quantum electrodynamics, describing photon interaction with natural
or artificial atoms coupled to the waveguide, is now rapidly developing [1, 2,
3]. This platform allows controllable generation of quantum light [4], and
control over lifetimes [5] and entanglement [6, 7] of coupled atom-photon
excitations. Even more possibilities are opened in the structures with the
parameters dynamically modulated in time [8]. This enables Floquet engineering
and realization of synthetic dimensions [9, 10] as well as control of quantum
photon-photon correlations [11, 12]. Such time modulation has been recently
demonstrated, for example, for the superconducting transmon qubit platform
[13].

One more fundamental physical effect, that becomes possible in the dynamically
modulated structures, is the parametric generation of photon pairs. Such
generation can be also interpreted as a dynamical Casimir effect. This effect
has been first proposed for the cavity with a moving wall [14]. To the best of
our knowledge, it has been never directly observed in this setup so far,
because of the extremely low photon generation rate at realistic parameters,
see the review [15]. However, there also exist generalized dynamical Casimir
effects, where other electromagnetic properties of the medium are changing
instead of physical movement of the mirror in space. For example, parametric
photon generation due to electron-hole plasma generated and moving in a
semiconductor under the laser pulse excitation has been theoretically
considered in Ref. [16]. A seminal observation of an analogue of the dynamical
Casimir effect in a superconducting circuit was made in Ref. [17]. The
effective length of the transmission line has been modulated by changing the
inductance of a superconducting quantum interference device. A detailed
theoretical analysis of such systems has been performed in Refs. [18, 19].
Two-photon entanglement in this setup has been experimentally studied in Ref.
[20] However, the consideration has been limited to the case when the
modulated elements were not resonant for the generated photons. Dynamical
Casimir effect in the arrays of resonant modulated emitters is not yet
explored. The resonant structures could potentially allow one to selectively
enhance the photon generation and control the correlations between them.
Recently, a theory of dynamical Casimir effect in the general dynamically
modulated photonic structures has been put forward in Ref. [21]. However,
considered setup did not include two-photon interactions, essential for the
qubit platform.

Here, we consider a generalized dynamical Casimir effect for an array of
qubits with strongly dynamically modulated resonant frequencies, coupled to
the waveguide. In such a system there exist resonances for generated photons
and photon pairs near the single- and double-excited levels of the qubits. Our
goal is to explore the role of the various collective qubit resonances for the
intensity, directionality and the quantum correlations between the emitted
photons.

Figure 1: Scheme of the structure under consideration. The array of qubits
with modulated in time resonant frequencies is coupled to a waveguide.

The rest of the paper is organized as follows. Our theoretical model and
calculation approach, based on the master equation, are presented in Sec. II.
We discuss the numerical and analytical results for a single qubit in Sec.
III. The arrays with $N>1$ qubits are considered in Sec. IV and the main
results are summarized in Sec. V. Appendix A presents an equivalent
alternative Green-function-based diagrammatic approach to calculate the photon
emission intensity and the photon-photon correlation functions. Analytic
results obtained by this approach for a particular case of $N=2$ qubits are
given in Appendix B.

##  II Model

The structure under consideration is schematically illustrated in Fig. 1. It
consists of $N$ periodically spaced qubits in the one-dimensional waveguide.
The resonant frequencies of the qubits are modulated in time with the
modulation frequency $\Omega$. The system Hamiltonian can be written as
$H=H_{0}+\sum_{j}^{N}V_{j}$, where

$H_{0}=\sum_{j=1}^{N}\left(\omega_{0}a_{j}^{\dagger}a_{j}+\frac{U}{2}a_{j}^{\dagger}a_{j}^{\dagger}a_{j}a_{j}\right)+\sum_{j,k=1}^{N}\mathop{\mathrm{missing}}{Re}\nolimits(D_{jk})a_{j}^{\dagger}a_{k}$
(1)

is the unperturbed Hamiltonian and $a_{j}$ are the bosonic annihilation
operators for the qubit excitations with the frequency $\omega_{0}$, located
at the point $x_{j}=d(j-1)$. To observe the generation of photon pairs, we
consider the general case of anharmonic multilevel qubits with the finite
anharmonicity $U$. The radiative coupling and collective decay of the qubits
are described by the photon Green function

$D_{jk}=-{\rm i}\gamma_{1\rm{D}}e^{{\rm i}\omega_{0}|x_{j}-x_{k}|/c}\;,$ (2)

where $\gamma_{1\rm{D}}$ is the single qubit radiative decay rate. The
perturbation responsible for modulating the frequency of each qubit reads

$V_{j}=g_{j}(a_{j}^{\dagger}+a_{j})^{2}\cos(\Omega t+\varphi_{j})\;.$ (3)

Here, $g_{j}$ and ${\varphi}_{j}$ are the amplitude and the phase of the
modulation, respectively. The time evolution of the system is described by the
master equation [1, 11]

$\displaystyle\dot{\rho}=-{\rm
i}[H,\rho]+\sum_{j,k=1}^{N}\gamma_{jk}(2a_{j}\rho
a_{k}^{\dagger}-a_{j}^{\dagger}a_{k}\rho-\rho a_{j}^{\dagger}a_{k})\;,$ (4)

where
$\gamma_{jk}=-\mathop{\mathrm{missing}}{Im}\nolimits(D_{jk})+\delta_{jk}\gamma$
and $\gamma$ is the nonradiative decay rate.

In order to find a stationary solution for the density matrix, we represent
$\rho$ as a vector $\bm{\rho}$. Then the Eq. (4) can be written as

$\dot{\bm{\rho}}=\mathcal{L}\bm{\rho}+\frac{1}{2}\left(\mathcal{V}e^{{\rm
i}\Omega t}+\mathcal{V}^{*}e^{-i\Omega t}\right)\bm{\rho}\;,$ (5)

where the operators $\mathcal{L}$ and $\mathcal{V}$ read

$\displaystyle\mathcal{L}$ $\displaystyle=-{\rm i}(H_{0}\otimes 1-1\otimes
H_{0})$ (6) $\displaystyle+\sum_{j,k=1}^{N}\gamma_{jk}(2a_{j}\otimes
a_{k}-a_{j}^{\dagger}a_{k}\otimes 1-1\otimes a_{k}^{\dagger}a_{j})\;,$
$\displaystyle\mathcal{V}$ $\displaystyle=-{\rm i}\sum_{j=1}^{N}g_{j}e^{{\rm
i}\varphi_{j}}\left((a_{j}^{\dagger}+a_{j})^{2}\otimes
1-1\otimes(a_{j}^{\dagger}+a_{j})^{2}\right)\;,$

and the symbol $\otimes$ denotes the Kronecker product. Representing the
density matrix in the form of a Fourier series $\bm{\rho}(t)=\sum_{n}e^{-{\rm
i}n\Omega t}\bm{\rho}_{n}$, leaving only the harmonics with $n=-1,0,1$ and
substituting it into Eq. (6) we obtain a system of equations for the
stationary density matrix $\bm{\rho}_{0}$

$\mathcal{L}\bm{\rho}_{0}=\frac{1}{4}\left(\mathcal{V}\left(\mathcal{L}+{\rm
i}\Omega\right)^{-1}\mathcal{V}^{*}+\mathcal{V}^{*}\left(\mathcal{L}-{\rm
i}\Omega\right)^{-1}\mathcal{V}\right)\bm{\rho}(0)\;.$ (7)

Here, a vector $\bm{\rho}(0)$ corresponds to the initial vacuum state.

Figure 2: (a) Emission spectrum $S$ of one qubit as a function of the
modulation frequency $\Omega$ and the emission frequency $\omega$. The
calculation parameters are
$g/\gamma_{1\rm{D}}=0.1,\;U/\gamma_{1\rm{D}}=10,\;\gamma=0$. Panel (b) shows
the emission spectrum obtained as the cross section of the density plot (a)
for $\Omega=2\omega_{0}+U$.

##  III Single qubit

We start by considering a single qubit coupled to a waveguide. We are
interested in the weak driving regime when $g\ll\gamma_{1\rm{D}}$. In this
case, to find the observable quantities, it is sufficient to restrict the
consideration to a three-level nonlinear qubit. The stationary density matrix
$\rho_{0}$ can be easily found from Eq. (7). It yields the photon emission
intensity

$\displaystyle I_{1}$ $\displaystyle=\langle
a^{\dagger}a\rangle=\mathop{\mathrm{missing}}{Tr}\nolimits(a^{\dagger}a\rho_{0})$
(8)
$\displaystyle=\frac{4g^{2}(\Omega^{2}+(U+2\omega_{0})^{2}+4\gamma_{\Sigma}^{2})}{[4\gamma_{\Sigma}^{2}+(\Omega-U-2\omega_{0})^{2}][4\gamma_{\Sigma}^{2}+(\Omega+U+2\omega_{0})^{2}]}\;,$

where $\gamma_{\Sigma}=\gamma_{1\rm{D}}+\gamma$. The highest intensity is
achieved when the modulation frequency is in resonance with the transition
between the ground level and the second excited level of the qubit,
$\Omega=2\omega_{0}+U$.

Another important characteristic is the spectrum of the photon emission. It
can be found using the quantum regression theorem as
$S=\frac{1}{\pi}\mathop{\mathrm{missing}}{Re}\nolimits\int_{0}^{\infty}{\rm
e}^{-{\rm
i}\omega\tau}\mathop{\mathrm{missing}}{Tr}\nolimits(a^{\dagger}e^{\mathcal{L}\tau}a\rho_{0})$.
In this way we get the spectrum

$\displaystyle
S(\omega)=\frac{I_{1}}{2\pi}\frac{\gamma_{\Sigma}}{U^{2}+4\gamma_{\Sigma}^{2}}\left(\frac{8\gamma^{2}_{\Sigma}+U[U+2(\omega-\omega_{0})]}{(\omega-\omega_{0})^{2}+\gamma_{\Sigma}^{2}}\right.$
(9)
$\displaystyle\left.+\frac{U[3U-2(\omega-\omega_{0}-U)]}{(\omega-\omega_{0}-U)^{2}+(3\gamma_{\Sigma})^{2}}\right)\>.$

The result of the calculation according to Eq. (9) is shown in Fig. 2. As
mentioned above, the maximum integral intensity is reached at the modulation
frequency $\Omega=2\omega_{0}+U$. The emission spectrum consists of two peaks,
one at the qubit frequency $\omega_{0}$, which corresponds to the transition
between the ground and first excited levels of the qubit, and the other peak
at the frequency $\omega_{0}+U$, which corresponds to the transition between
the first and second excited levels. It should be noted that in the case of
$U\gg\gamma_{1\rm D}$ the latter peak has three times smaller amplitude but is
three times wider that the former one, so the areas below the two peaks are
equal. This reflects the fact that the photons are born and emitted in pairs
with the average energies $\omega_{0}$ and $\omega_{0}+U$.

Figure 3: The photon emission intensity to the left $I_{-}=\langle
p_{-}^{\dagger}p_{-}\rangle$ for two qubits as a function of a distance
between the qubits $d$ and the modulation frequency $\Omega$. Calculations
have been performed in the case of symmetric modulation ($\varphi=0$) for (a)
$U/\gamma_{1\rm{D}}=1$, (b) $U/\gamma_{1\rm{D}}=4$, (c)
$U/\gamma_{1\rm{D}}=10$ and (d) antisymmetric modulation ($\varphi=\pi$) for
$U/\gamma_{1\rm{D}}=10$. Other calculation parameters are
$g/\gamma_{1\rm{D}}=0.1,\gamma/\gamma_{1\rm{D}}=0.1$. Gray lines represent the
double-excited eigenstates. Red dotted lines show the intensity minima. Figure
4: The photon emission intensity to the left $I_{-}$ for a pair of qubits as a
function of the distance $d$ and the modulation phase $\varphi$ for the
modulation frequency $\Omega=2\omega_{0}+U$ and different $U$ indicated in the
figure. Calculation has been performed for $\gamma/\gamma_{1\rm D}=0.05$ and
$g/\gamma_{1\rm D}=0.1$.

##  IV Qubit arrays

We now proceed to the discussion of the pair of qubits. We assume that the
coupling constants are equal, $g_{1}=g_{2}=g$, and that the phases are
$\varphi_{1}=0$, $\varphi_{2}=\varphi$. Unlike the case of a single qubit,
this system can emit directionally because of the interference between photons
from the first and second qubit. One can define the qubit lowering operators
corresponding to emission of photons in the left and right direction as
$p_{\mp}=a_{1}+a_{2}{\rm e}^{\pm{\rm i}\omega_{0}d/c}$. Then, the photon
emission intensity to the left $I_{-}=\langle p_{-}^{\dagger}p_{-}\rangle$ can
be obtained by using the stationary density matrix from Eq. (7).

Figure 5: The second order correlation function for a pair of qubits in the
case of (a) symmetric $(\varphi=0)$ and (b) antisymmetric $(\varphi=\pi)$
modulation. Calculation has been performed for the following set of
parameters: $U/\gamma_{1\rm{D}}=10$, $g/\gamma_{1\rm{D}}=0.1$,
$\gamma/\gamma_{1\rm{D}}=0.1$.

Figure 3 presents the calculation of $I_{-}$ for the three values of the
parameter $U/\gamma_{1\rm D}$. Panels (a), (b) and (c) correspond to the case
of symmetric modulation, when both qubits are excited in phase ($\varphi=0$).
In this case, we can see that the sharp intensity minima appear in the
vicinity of red dotted lines which correspond to dark single-excited states.
The highest intensity is achieved when the modulation frequency is close to
the real parts of the energies of the certain double-excited eigenstates (gray
lines). The latter are found by diagonalizing the non-Hermitian two-photon
Hamiltonian $H\otimes 1+1\otimes H+V$, where $H$ and $V$ are defined in
Appendix A. If $U/\gamma_{1\rm D}<4$ [panel (a)], the energies of the three
eigenstates are close and intersect when $\omega_{0}d/c$ is varied. As
$U/\gamma_{1\rm D}$ is increased, at a certain threshold [panel (b),
$U/\gamma_{1\rm D}=4$] the states are rearranged into two isolated bands, and
that get separated by the gap of $U$ for large $U/\gamma_{1\rm D}$ [panel
(c)]. The lower band ($\Omega\approx 2\omega_{0}$) is formed by the states
that comprise pairs of excitations in different qubits, e.g.,
$a_{1}^{\dagger}a_{2}^{\dagger}|0\rangle$. Such states cannot be excited by
the perturbation operator Eq. (3), that includes the products of operators two
$a_{j}^{\dagger}$ with the same $j$ only. Therefore, the emission intensity
near the lower band is quenched, see also Appendix B for the complete
analytical expression. The case of antisymmetric modulation ($\varphi=\pi$) is
shown in the panel (d). Here, we see no minima at the frequency of dark
single-excited states and the overall dependence on distance $d$ is rather
weak.

The dependence of the emission on the modulation phase is shown in more detail
in Fig. 4, which presents the calculation of $I_{-}$ for the modulation
frequency $\Omega=2\omega_{0}+U$ that corresponds to the resonance for a
single qubit. In the case of $U/\gamma_{1\rm D}=1$ [panel (a)] the maximum is
observed at $\omega_{0}d/c\approx\pi$ which corresponds to the resonance with
a double-exited eigenstate, cf. Fig. 3(a). As the ratio $U/\gamma_{1\rm D}$
increases, the structure of $I_{-}$ becomes more smooth along the $d$ axis,
due to the decreasing dependence of the eigenstates on the distance $d$. In
contrast, the dependence of $I_{-}$ on $\varphi$ becomes more pronounced, see
Fig. 3(b)-(d). Note that the photon emission intensity to the right $I_{+}$
can be obtained by flipping the sign of $\varphi$. The obtained maps,
especially Fig. 3(b), are highly asymmetric with respect to this operation,
which indicates the high directivity of the emission. In the case
$U/\gamma_{1D}\gg 1$, which is shown in panel (d), the answer can be obtained
analytically. At $\Omega=2\omega_{0}+U$, we can neglect the state when
different qubits are excited simultaneously. Then, we find

$I_{-}=I_{1}\left[2+\frac{\sin\varphi\sin(2\omega_{0}d/c)}{3-\cos(2\omega_{0}d/c)}\right]$
(10)

in the absence of the nonradiative decay rate $\gamma=0$, see also Appendix B.
The second term in the right-hand side of Eq. (10) describes the interference
between the photon pairs generated by the two qubits. Importantly, the
interference term is odd in $\varphi$, thus its contribution is opposite for
$I_{-}$ and $I_{+}$ which results in a directional emission. The maximal
degree of directivity $(I_{-}-I_{+})/(I_{-}+I_{+})=\sqrt{2}/8\approx 0.18$ is
achieved for $\varphi=\pi/2$ and
$\omega_{0}d/c=\arctan(2\sqrt{2})/2\approx\pi/5$. The intensity reduction in
the corners of the panels in Fig. 4 is out of the scope of Eq. (10) and
explained by the finite nonradiative decay that was taken into account in the
calculation.

We also present in Fig. 5 the (unnormalized) photon-photon correlation
function $G^{(2)}_{--}=\langle
p_{-}^{\dagger}p_{-}^{\dagger}p_{-}p_{-}\rangle$ for symmetric [panel (a)] and
antisymmetric [panel (b)] modulation of the qubits. The calculation
demonstrates that the largest values of $G^{(2)}_{--}$ are achieved at the
two-photon resonance $\Omega=2\omega_{0}+U$. In order to explain the
calculated dependence of the correlation function on the distance between the
qubits $d$ we have obtained an analytical expression similar to Eq. (10) and
valid for $U/\gamma_{1\rm D}\gg 1$ ,

$G^{(2)}_{--}=I_{1}\left[1+\cos(2\omega_{0}d/c-\varphi)\right]\>,$ (11)

see also Appendix B. The simple form of Eq. (11) describes the interference of
the two independent coherent sources separated by distance $d$ and emitting
with the phase difference $\varphi$. The fact that the two-photon emission and
detection is considered is accounted by the factor of 2 in the phase
$2\omega_{0}d/c$ that the photon pair gains when travelling the distance $d$.
If the qubits are modulated in phase ($\varphi=0$), the maxima of Eq. (11) are
realized for the periods when $2\omega_{0}d/c=0,\pi,2\pi$ etc. The out-of-
phase modulation, when $\varphi=\pi/2$ corresponds to the maxima at
$2\omega_{0}d/c=\pi/2,3\pi/2\ldots$. This agrees with the numerical
calculations in Fig. 5. For symmetric modulation there exists also an
additional minimum in Fig. 5(a), corresponding to a strong antibunching.
According to our analytical expression Eq. (36), this minimum corresponds to
the frequency $\Omega=2\omega_{0}-2\gamma_{\rm 1D}\tan qd$. The corresponding
expression is shown by a red dotted line in Fig. 5(a) and well describes
numerical results.

Figure 6: (a) The second order correlation function for an array of four
qubits in the case of frequency modulation of the first qubit only. Panel (b)
displays the area highlighted with a green square on the panel (a). Red dotted
curves present the real parts of the energies of double-excited states.
Calculation has been performed for the following set of parameters:
$U/\gamma_{1\rm{D}}=10$, $g/\gamma_{1\rm{D}}=0.1$,
$\gamma/\gamma_{1\rm{D}}=0$.

We have also calculated the two-photon correlation function for a larger array
with $N=4$ qubits. This is the minimum qubit number required to have double-
excited subradiant states [22, 23, 5]. Such states have long radiative
lifetime because of the destructive interference in the spontaneous photon
emission processes. They exist for $\omega_{0}d/c\ll 1$ (or
$|\omega_{0}d/c-\pi|\ll 1$), and their lifetime is enhanced by the factors on
the order of $1/(\omega_{0}d/c)^{2}$ ($1/|\omega_{0}d/c-\pi|^{2}$). Hence, we
can expect appearance of additional sharp spectral features in $G^{(2)}_{--}$
for $N=4$ due to the double-excited subradiant states. The corresponding color
plots of calculated $G^{(2)}_{--}$ are shown in Fig. 6. Correlation function
$G^{(2)}_{--}$ has maxima around $\Omega=2\omega_{0}+U$, similarly to the case
of $N=2$ qubits. As expected, there also appear minima at $\Omega\approx
2\omega_{0}$ when the interqubit distance is either small or close to
$c\pi/\omega_{0}$. The map of $G^{(2)}_{--}$ in the region of small
$\omega_{0}d/c$ is shown in pannel (b) in detail. One can see that on top of
the smooth minimum two sharp maxima appear [see two lower dotted lines in Fig.
6(b)]. Their position matches the energies of the two subradiant states
$2\omega_{0}-2\gamma_{1\rm D}\omega_{0}d/c$ and
$2\omega_{0}-(14/3)\gamma_{1\rm D}\omega_{0}d/c$ [22].

##  V Summary

To summarize, we have developed a general theory of parametric photon
generation from arrays of modulated in time qubits coupled to the waveguide.
Using the two independent approaches, the master equation for the density
matrix and the diagrammatic Green-function technique, we have studied the
dependence of the photon emission spectrum and photon-photon correlation
functions on the anharmonicity of the qubit potential $U$, the distance
between the neigboring qubits and the relative modulation phase $\varphi$.

The calculated emission spectrum is very sensitive to the qubit anharmonicity
parameter. The anharmonicity controls the relative weight of the spectral
features around the single- and double-excited qubit resonances. The latter
become more prominent with the increase of the anharmonicity. When the number
of qubits is $N=4$ or larger, additional sharp spectral features,
corresponding to the double-excited subradiant states, appear in the spectrum.
We have also shown that the interference between photons emitted from
different qubits can be controlled by relative phases of their frequency
modulation. Our calculation demonstrates how this can be used to obtain
directional photon pair emission, similarly as it happens for non-parametric
quantum photon sources [24, 13, 25].

We hope that our results will be useful for engineering the parametric quantum
emission from the waveguide-coupled qubit arrays. A potentially interesting
future research direction could be the system, where the spatial position of
the light emitters, rather than their resonance frequency, oscillates in time.
This would mean generalization of our concept of optomechanical Kerker effect
[26], that is motion-induced direcitonal emission, to the quantum optics
regime.

###### Acknowledgements.

We thank E.S. Redchenko for useful discussions. This work has been funded by
the Russian Science Foundation Grant No. 20-12-00194.

##  Appendix A Diagrammatic approach

In case of weak modulation, only the states with small number of excitations
are populated. They can be described in the framework of perturbative
diagrammatic approach. The consideration generalizes the results of Refs. [27,
28, 22] and also Ref. [29] for the structures, modulated in time. Somewhat
similar consideration for time-modulated structures is also available in Ref.
[11], but that work dows not consider parametric photon generation.

The diagram representing the (not normalized) wavefunction of the generated
photon pair is shown in Fig. 7(a). The dashed line represents the modulation
that creates a pair of qubit excitations (solid lines). Then, the excitations
propagate in the structure, which is described by the Green functions (solid
lines). The excitations can either get converted to the photons directly
(first diagram) or interact and get converted to photons after that (second
diagram). For the case when both photons are emitted in the left direction
this yields

$\displaystyle\psi_{--}(\omega_{1},\omega_{2})=$ $\displaystyle\gamma_{\rm
1D}\sum_{ij}s_{i}^{+}(\omega_{1})s_{j}^{+}(\omega_{2})[\delta_{ij}+(M\Sigma)_{ij}]g_{j}$
$\displaystyle\times 2\pi\delta(\omega_{1}+\omega_{2}-\Omega)\,.$ (12)

Here, the Green function of single excitation reads

$\displaystyle G(\omega)=(\omega-H)^{-1},\qquad
H_{ij}=\omega_{0}\delta_{ij}-D_{ij}\,,$ (13)

the outer lines of the diagrams correspond to

$\displaystyle s_{i}^{+}(\omega)=\sum_{j}G_{ij}(\omega){\rm e}^{{\rm
i}qz_{j}}=\frac{{\rm e}^{{\rm i}qz_{1}}}{{\rm i}\gamma_{\rm
1D}}\left(\frac{\omega_{0}-H}{\omega-H}\right)_{1i}\,,$ (14)

and the propagation of a pair of excitations is described by

$\displaystyle\Sigma_{ij}(\Omega)={\rm i}\int
G_{ij}(\omega)G_{ij}(\Omega-\omega)\,\frac{d\omega}{2\pi}\,.$ (15)

The dressed interaction vertex $M$ can be determined from the Dyson-like
equation shown diagrammatically in Fig. 7(b), which yields

$\displaystyle M(\Omega)=[U^{-1}-\Sigma(\Omega)]^{-1}\,.$ (16)

Substituting this into Eq. (A) we finally obtain

$\displaystyle\psi_{--}(\omega_{1},\omega_{2})=$
$\displaystyle\frac{\gamma_{\rm
1D}}{U}\sum_{ij}s_{i}^{+}(\omega_{1})s_{j}^{+}(\omega_{2})M_{ij}g_{j}$
$\displaystyle\times 2\pi\delta(\omega_{1}+\omega_{2}-\Omega)\,.$ (17)

The integration in Eq. (15) can be performed explicitly, which yields a
compact expression for the $M$ matrix:

$\displaystyle M_{ij}(\Omega)=\left[\frac{\Omega-H\otimes 1-1\otimes
H}{\Omega-H\otimes 1-1\otimes H-V}\,V\right]_{ii,jj}$ (18)

where $(A\otimes B)_{ij,kl}\equiv A_{ik}B_{jl}$, $\frac{A}{B}\equiv AB^{-1}$
and we introduced the diagonal $N^{2}\times N^{2}$ matrix
$V_{ij,kl}=\delta_{ij}\delta_{kl}\delta_{ik}U$.

Figure 7: (a) Diagrams representing the amplitude of two-photon generation by
modulated qubits. Dashed line represents the modulation, solid lines are the
Green functions of qubit excitations, wavy lines are the outgoing photons. (b)
The diagrammatic equation for the dressed vertex (solid circle) that describes
the interaction of two qubit excitations. Open circle represents the bare
vertex corresponding to the interaction amplitude $U$.

###  A.1 Correlation function

The second-order correlation function of the emitted photon pair is calculated
as

$\displaystyle G^{(2)}_{--}(\tau)=\frac{1}{\gamma_{\rm
1D}^{2}}\left|\int\psi_{--}(\omega_{1},\omega_{2}){\rm e}^{-{\rm
i}\omega_{1}t-{\rm
i}\omega_{2}(t+\tau)}\frac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\right|^{2}\>.$
(19)

At zero delay we use the representation Eq. (18) to obtain

$\displaystyle G^{(2)}_{--}(0)=\frac{1}{\gamma_{\rm
1D}^{4}}\left|\sum_{i}\left[\frac{(\omega_{0}-H)\otimes(\omega_{0}-H)}{\Omega-H\otimes
1-1\otimes H-V}\right]_{11,ii}g_{i}\right|^{2}$ (20)

In the limit $U\gg\gamma_{\rm 1D}$ and for
$|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm 1D}$, we keep in the denominator of
Eq. (20) only the diagonal terms and find

$\displaystyle G^{(2)}_{--}(0)=\dfrac{\left|\sum\limits_{i}g_{i}\,{\rm
e}^{2{\rm i}qz_{i}}\right|^{2}}{(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm
1D}^{2}}\,,$ (21)

where we used $(\omega_{0}-H)_{1i}={\rm i}\gamma_{\rm 1D}{\rm e}^{{\rm
i}q(z_{i}-z_{1})}$.

For $|\Omega-\omega_{0}|\sim\gamma_{\rm 1D}$, the result is

$\displaystyle
G^{(2)}_{--}(0)=\frac{1}{U^{2}}\left|\sum_{ij}\Sigma^{+}_{i}(\Omega)[\Sigma^{-1}(\Omega)]_{ij}g_{j}\right|^{2},$
(22)

where

$\displaystyle\Sigma_{ij}(\Omega)=\left[\frac{1}{\Omega-H\otimes 1-1\otimes
H}\right]_{ii,jj},$ (23)
$\displaystyle\Sigma^{+}_{i}(\Omega)=\frac{1}{\gamma_{\rm
1D}^{2}}\left[\frac{(\omega_{0}-H)\otimes(\omega_{0}-H)}{\Omega-H\otimes
1-1\otimes H}\right]_{11,ii}.$ (24)

###  A.2 Emission intensity

We define the intensity of emission to the left as

$\displaystyle
I_{-}=\int\left[|\psi_{--}(\omega_{1},\omega_{2})|^{2}+|\psi_{-+}(\omega_{1},\omega_{2})|^{2}\right]\frac{d\omega_{1}d\omega_{2}}{\gamma_{\rm
1D}T(2\pi)^{2}}\>,$ (25)

where $T$ is the normalization time and $\psi_{-+}$ is the wave function of a
pair of photons emitted in the opposite directions. The latter is obtained
from Eq. (A) by replacing the factor $s_{j}^{+}(\omega_{2})$ with
$s_{j}^{-}(\omega_{2})=\sum_{i}G_{ij}(\omega){\rm e}^{-{\rm i}qz_{i}}$.

Using the representation Eq. (18), we get

$\displaystyle\int|\psi_{--}(\omega_{1},\omega_{2})|^{2}\tfrac{d\omega_{1}d\omega_{2}}{T(2\pi)^{2}}=\frac{1}{{\rm
i}\gamma_{\rm 1D}^{2}}\sum_{ij}g_{i}g_{j}^{*}$ $\displaystyle\Big{[}$
$\displaystyle\frac{(\omega_{0}-H_{1})(\omega_{0}-H_{2})(\omega_{0}-H_{3}^{*})(\omega_{0}-H_{4}^{*})}{(\Omega-
H_{1}-H_{4}^{*})}$ (26) $\displaystyle\times$
$\displaystyle\left(\frac{\Omega-
H_{3}^{*}-H_{4}^{*}}{H_{1}-H_{3}^{*}}-\frac{\Omega-
H_{1}-H_{2}}{H_{4}^{*}-H_{2}}\right)$ $\displaystyle\times\frac{1}{(\Omega-
H_{1}-H_{2}-V_{12})(\Omega-H_{3}^{*}-H_{4}^{*}-V_{34})}\Big{]}_{1111,iijj}$

and the corresponding integral of $|\psi_{-+}|^{2}$ is obtained in the same
way by taking the element $[\ldots]_{1N1N,iijj}$. Here $H_{1}=H\otimes
1\otimes 1\otimes 1$, $H_{2}=1\otimes H\otimes 1\otimes 1$ etc.

In the limit $U\gg\gamma_{\rm 1D}$ and for
$|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm 1D}$, the result can be simplified:

$\displaystyle\int|$
$\displaystyle\psi_{--}(\omega_{1},\omega_{2})|^{2}\tfrac{d\omega_{1}d\omega_{2}}{T(2\pi)^{2}}=\frac{2{\rm
i}}{\gamma_{\rm 1D}^{2}[(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm 1D}^{2}]}$
$\displaystyle\times\sum_{ij}g_{i}g_{j}^{*}(\omega_{0}-H)_{1i}(\omega_{0}-H^{*})_{1j}{\mathcal{Q}}_{11,ij}\>,$
(27) $\displaystyle\int|$
$\displaystyle\psi_{-+}(\omega_{1},\omega_{2})|^{2}\tfrac{d\omega_{1}d\omega_{2}}{T(2\pi)^{2}}=\frac{{\rm
i}}{\gamma_{\rm 1D}^{2}[(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm 1D}^{2}]}$
$\displaystyle\times$
$\displaystyle\sum_{ij}g_{i}g_{j}^{*}\Big{\\{}(\omega_{0}-H)_{1i}(\omega_{0}-H^{*})_{1j}\mathcal{Q}_{NN,ij}$
(28) $\displaystyle+$
$\displaystyle(\omega_{0}-H)_{Ni}(\omega_{0}-H^{*})_{Nj}\mathcal{Q}_{11,ij}\Big{\\}}\>,$

where

$\mathcal{Q}=\frac{(\omega_{0}-H)\otimes(\omega_{0}-H^{*})}{(\omega_{0}-H)\otimes
1-1\otimes(\omega_{0}-H^{*})}\>.$ (29)

###  A.3 Optical theorem

We assume that the modulation of the qubits is performed by an external
signal. The back-action of the qubit system on that signal can be accounted by
introducing the scattering parameter $s$. Up to the second order in $g$, it
reads

$\displaystyle s=1-{\rm i}\sum_{ij}g_{i}^{*}\tilde{\Sigma}_{ij}g_{j},$ (30)
$\displaystyle\tilde{\Sigma}_{ij}=\left[\frac{1}{\Omega-H\otimes 1-1\otimes
H-V}\right]_{ii,jj}\>.$ (31)

In the absence of non-radiative losses, the optical theorem imposes
$2|s|^{2}+I_{+}+I_{-}=2$. Therefore, we get

$\displaystyle I_{+}+I_{-}$ (32) $\displaystyle=-\frac{4}{\gamma_{1D}}{\rm
Im\,}\sum_{ij}g_{i}^{*}g_{j}\left[\frac{1}{\Omega-H\otimes 1-1\otimes
H-V}\right]_{ii,jj}\>.$

In the limit $U\gg\gamma_{\rm 1D}$ , $|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm
1D}$, the result is trivial

$\displaystyle I_{+}+I_{-}=\frac{8}{(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm
1D}^{2}}\,\sum_{i}|g_{i}|^{2}\>.$ (33)

##  Appendix B Two qubits

Here we apply the results of Appendix A for the system of two qubits. The
general explicit expressions are quite bulky so consider two special cases.

###  B.1 $U\gg\gamma_{\rm 1D}$ and $|\Omega-2\omega_{0}-U|\lesssim\gamma_{\rm
1D}$

In that limit, we get

$\displaystyle G^{(2)}_{--}(0)=\frac{\left|g_{1}+g_{2}{\rm e}^{2{\rm
i}qd}\right|^{2}}{(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm 1D}^{2}}$ (34)

and

$I_{-}=\frac{1}{(\Omega-2\omega_{0}-U)^{2}+4\gamma_{\rm 1D}^{2}}\\\
\times\frac{(7-3\cos 2qd)|g_{1}|^{2}+(5-\cos 2qd)|g_{2}|^{2}+2\sin
2qd\,\text{Im\,}g_{1}g_{2}^{*}}{3-\cos 2qd}\>.$ (35)

This yields Eqs. (10) and (11) in the main text.

###  B.2 Symmetric modulation, $g_{1}=g_{2}=g$

In that case, the two-photon emission is

$\displaystyle G^{(2)}_{--}(0)=4|g|^{2}$ (36)
$\displaystyle\times\left|\frac{(\Omega-2\omega_{0})\cos qd+2\gamma_{\rm
1D}\sin qd}{(\Omega-2\omega_{0}+2{\rm i}\gamma_{\rm
1D})(\Omega-2\omega_{0}-U+2{\rm i}\gamma_{\rm 1D})+4\gamma_{\rm 1D}^{2}{\rm
e}^{2{\rm i}qd}}\right|^{2}\>.$

Note that $G^{(2)}_{--}(0)$ vanishes at $\Omega=2\omega_{0}-2\gamma_{\rm
1D}\tan qd$.

The single-photon emission intensity is found easily from the optical theorem

$\displaystyle I_{+}=I_{-}=8|g|^{2}$ (37)
$\displaystyle\times\frac{(\Omega-2\omega_{0}+\gamma_{\rm 1D}\sin
2qd)^{2}+2\gamma_{\rm 1D}^{2}(3-\cos
2qd)\sin^{2}qd}{\left|(\Omega-2\omega_{0}+2{\rm i}\gamma_{\rm
1D})(\Omega-2\omega_{0}-U+2{\rm i}\gamma_{\rm 1D})+4\gamma_{\rm 1D}^{2}{\rm
e}^{2{\rm i}qd}\right|^{2}}\>.$

We note that $I_{\pm}$ has a minimum at $\Omega\approx 2\omega_{0}-\gamma_{\rm
1D}\sin 2qd$. In particular, for $qd=0$, we have
$I_{\pm}(\Omega=2\omega_{0})=0$.

## References

  * Sheremet  _et al._ [2023] A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, Waveguide quantum electrodynamics: collective radiance and photon-photon correlations, Rev. Mod. Phys. , in press (2023), arXiv:2103.06824 [quant-ph] .
  * Roy  _et al._ [2017] D. Roy, C. M. Wilson, and O. Firstenberg, Colloquium: strongly interacting photons in one-dimensional continuum, Rev. Mod. Phys. 89, 021001 (2017).
  * Chang  _et al._ [2018] D. E. Chang, J. S. Douglas, A. González-Tudela, C.-L. Hung, and H. J. Kimble, Colloquium: quantum matter built from nanoscopic lattices of atoms and photons, Rev. Mod. Phys. 90, 031002 (2018).
  * Prasad  _et al._ [2020] A. S. Prasad, J. Hinney, S. Mahmoodian, K. Hammerer, S. Rind, P. Schneeweiss, A. S. Sørensen, J. Volz, and A. Rauschenbeutel, Correlating photons using the collective nonlinear response of atoms weakly coupled to an optical mode, Nature Photonics 14, 719 (2020).
  * Zanner  _et al._ [2022] M. Zanner, T. Orell, C. M. F. Schneider, R. Albert, S. Oleschko, M. L. Juan, M. Silveri, and G. Kirchmair, Coherent control of a multi-qubit dark state in waveguide quantum electrodynamics, Nature Physics 18, 538 (2022).
  * Kannan  _et al._ [2020] B. Kannan, M. J. Ruckriegel, D. L. Campbell, A. F. Kockum, J. Braumüller, D. K. Kim, M. Kjaergaard, P. Krantz, A. Melville, B. M. Niedzielski, A. Vepsäläinen, R. Winik, J. L. Yoder, F. Nori, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Waveguide quantum electrodynamics with superconducting artificial giant atoms, Nature 583, 775 (2020).
  * Corzo  _et al._ [2019] N. V. Corzo, J. Raskop, A. Chandra, A. S. Sheremet, B. Gouraud, and J. Laurat, Waveguide-coupled single collective excitation of atomic arrays, Nature 566, 359 (2019).
  * Silveri  _et al._ [2017] M. P. Silveri, J. A. Tuorila, E. V. Thuneberg, and G. S. Paraoanu, Quantum systems under frequency modulation, Reports on Progress in Physics 80, 056002 (2017).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * Yuan  _et al._ [2018] L. Yuan, Q. Lin, M. Xiao, and S. Fan, Synthetic dimension in photonics, Optica 5, 1396 (2018).
  * Ilin  _et al._ [2023] D. Ilin, A. V. Poshakinskiy, A. N. Poddubny, and I. Iorsh, Frequency combs with parity-protected cross-correlations and entanglement from dynamically modulated qubit arrays, Phys. Rev. Lett. 130, 023601 (2023).
  * Jolin  _et al._ [2021] S. W. Jolin, G. Andersson, J. C. R. Hernández, I. Strandberg, F. Quijandría, J. Aumentado, R. Borgani, M. O. Tholén, and D. B. Haviland, Multipartite entanglement in a microwave frequency comb (2021).
  * Redchenko  _et al._ [2022] E. S. Redchenko, A. V. Poshakinskiy, R. Sett, M. Zemlicka, A. N. Poddubny, and J. M. Fink, Tunable directional photon scattering from a pair of superconducting qubits (2022), arXiv:2205.03293 .
  * Moore [1970] G. T. Moore, Quantum theory of the electromagnetic field in a variable-length one-dimensional cavity, J. Math. Phys. 11, 2679 (1970).
  * Dodonov [2020] V. Dodonov, Fifty years of the dynamical casimir effect, Physics 2, 67 (2020).
  * Lozovik  _et al._ [1995] Y. E. Lozovik, V. G. Tsvetus, and E. A. Vinogradov, Parametric excitation of vacuum by use of femtosecond laser pulses, Physica Scripta 52, 184 (1995).
  * Wilson  _et al._ [2011] C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, and P. Delsing, Observation of the dynamical Casimir effect in a superconducting circuit, Nature 479, 376 (2011).
  * Johansson  _et al._ [2010] J. R. Johansson, G. Johansson, C. M. Wilson, and F. Nori, Dynamical Casimir effect in superconducting microwave circuits, Phys. Rev. A 82, 052509 (2010).
  * Macrì  _et al._ [2018] V. Macrì, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Nonperturbative dynamical Casimir effect in optomechanical systems: Vacuum Casimir-Rabi splittings, Phys. Rev. X 8, 011031 (2018).
  * Schneider  _et al._ [2020] B. H. Schneider, A. Bengtsson, I. M. Svensson, T. Aref, G. Johansson, J. Bylander, and P. Delsing, Observation of broadband entanglement in microwave radiation from a single time-varying boundary condition, Phys. Rev. Lett. 124, 140503 (2020).
  * Sloan  _et al._ [2021] J. Sloan, N. Rivera, J. D. Joannopoulos, and M. Soljačić, Casimir light in dispersive nanophotonics, Phys. Rev. Lett. 127, 053603 (2021).
  * Ke  _et al._ [2019] Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, Inelastic scattering of photon pairs in qubit arrays with subradiant states, Phys. Rev. Lett. 123, 253601 (2019).
  * Zhang and Mølmer [2019] Y.-X. Zhang and K. Mølmer, Theory of subradiant states of a one-dimensional two-level atom chain, Phys. Rev. Lett. 122, 203605 (2019).
  * Kannan  _et al._ [2023] B. Kannan, A. Almanakly, Y. Sung, A. D. Paolo, D. A. Rower, J. Braumüller, A. Melville, B. M. Niedzielski, A. Karamlou, K. Serniak, A. Vepsäläinen, M. E. Schwartz, J. L. Yoder, R. Winik, J. I.-J. Wang, T. P. Orlando, S. Gustavsson, J. A. Grover, and W. D. Oliver, On-demand directional microwave photon emission using waveguide quantum electrodynamics, Nature Physics 16, 247 (2023).
  * Poddubny and Golub [2021] A. N. Poddubny and L. E. Golub, Ratchet effect in frequency-modulated waveguide-coupled emitter arrays, Phys. Rev. B 104, 205309 (2021).
  * Poshakinskiy and Poddubny [2019] A. V. Poshakinskiy and A. N. Poddubny, Optomechanical Kerker effect, Phys. Rev. X 9, 011008 (2019).
  * Fang  _et al._ [2014] Y.-L. L. Fang, H. Zheng, and H. U. Baranger, One-dimensional waveguide coupled to multiple qubits: photon-photon correlations, EPJ Quantum Technology 1, 3 (2014).
  * Poshakinskiy and Poddubny [2016] A. V. Poshakinskiy and A. N. Poddubny, Biexciton-mediated superradiant photon blockade, Phys. Rev. A 93, 033856 (2016).
  * Poddubny  _et al._ [2016] A. N. Poddubny, I. V. Iorsh, and A. A. Sukhorukov, Generation of photon-plasmon quantum states in nonlinear hyperbolic metamaterials, Phys. Rev. Lett. 117, 123901 (2016).


# Using ontology embeddings for structural inductive bias in gene expression
data analysis

Maja Trębacz Zohreh Shams Mateja Jamnik Paul Scherer
Nikola Simidjievski Helena Andres Terre Pietro Liò
Department of Computer Science and Technology
University of Cambridge, UK

###### Abstract

Stratifying cancer patients based on their gene expression levels allows
improving diagnosis, survival analysis and treatment planning. However, such
data is extremely highly dimensional as it contains expression values for over
20000 genes per patient, and the number of samples in the datasets is low. To
deal with such settings, we propose to incorporate prior biological knowledge
about genes from ontologies into the machine learning system for the task of
patient classification given their gene expression data. We use ontology
embeddings that capture the semantic similarities between the genes to direct
a Graph Convolutional Network, and therefore sparsify the network connections.
We show this approach provides an advantage for predicting clinical targets
from high-dimensional low-sample data.

##  1 Introduction

Ontologies are structured representations of semantic knowledge commonly used
to represent biological concepts. For example, the Gene Ontology (GO) [1]
unifies the representation of genes and their functions across all species.
Ontologies prove useful for developing machine learning systems operating on
biological entities [2, 3, 4], and they have an advantage of integrating data
from all the omics (genomes, transcriptomics, proteomics, metabolomics, etc.).
Recently developed ontology embeddings [5, 6, 7, 8] are structure-preserving
maps from ontologies into vector spaces. The embeddings capture the semantic
similarity between the genes. For example, genes TP53 and FOXA1 are close in
the ontology space as they are transcription factors annotated to common
molecular functions and biological processes, including terms related to
cancer. We incorporate such prior knowledge into a machine learning system
that classifies cancer subtypes from gene expressions.

The development of novel deep learning models on graphs [9, 10, 11, 12]
combined with biological motivations behind the network propagation model [13]
has created an opportunity to leverage the prior biomedical knowledge without
overly restricting the expressive power of the neural model. Rhee et al. [14]
proposed a combination of graph neural network (GNN) that uses a topology of
protein-protein interaction (PPI) with a relation network [15] to perform
convolutions on the gene expressions and classify cancer patients into PAM50
subtype [16]. A similar work of Dutil et al. [17] explores the usage of Graph
Convolutional Network (GCN) directed by PPI networks [18] on the single gene
inference tasks. The follow-up works by Bertin et al. [19] and Hashir et al.
[20] concluded that gene expression dependencies can be predicted almost as
well by using random networks as by using biological networks. One of the
reasons for this is that the degree of the nodes varies significantly and the
predictions tend to be worse for the genes with a low number of neighbours (or
no neighbours at worst). Such nodes get less signal and thus tend to lead to
lower performance. Crawford et al. [21] explored this observation and showed
that when low-degree genes are removed, then the biological networks yield
better performance than random graphs.

This work takes a different approach and is the first to use the ontology
embeddings as a similarity measure between the genes in order to impose the
inductive bias for the patient classification tasks from gene expressions. We
extend the gene expression convolution work of Dutil et al. [17] by proposing
novel graphs automatically generated from ontology embeddings, and a node
embedding method combining knowledge-based gene embeddings with the expression
values. Automatic generation of graphs expressing semantic similarity between
genes allows overcoming the problems with nodes that have a small number of
neighbours. Moreover, the ontology-based feature selection allows selecting a
biologically relevant set of features.

##  2 Methods

We built the OntoGCN (ontology-directed Graph Convolutional Network) neural
model where known similarities/relationships between the features (genes)
direct processing in the network and help the model to avoid learning spurious
correlations [17]. OntoGCN enforces convolutions on the genes related by
similarity and thus captures localised patterns of data, similarly as
convolutional neural networks capture spatial relationships of pixels in the
images [22].

Figure 1: Overview of the OntoGCN applied to gene expression data for patient
classification.

Figure 1 illustrates the workflow of our OntoGCN neural model. Each gene
contributes a node to a knowledge graph where edges represent the similarity
between the genes. We create edges that connect each gene to its K-nearest
neighbours in the ontology embedding space according to the cosine distance.
We use DL2vec ontology embeddings that are a graph-based method learning gene
representations over three biomedical ontologies (GO [1], UBERON [23], and MP
[24]). The topology of the graph is the same for all the patient samples in
the dataset, but each patient spans a new instance of the graph with different
expression values at the nodes. At each graph convolution step, the
neighbouring nodes are aggregated together based on their connectivity (Figure
1 shows the update step for the ESR1 node, but similar message-passing happens
for every node). After convolution, the gene nodes are dropped at random and
also pooled using topology-based aggregation clustering. Finally, a prediction
is made from the remaining nodes via a fully connected layer.

A possible advantage of using a graph generated based on similarity, rather
than a curated PPI network, is that one can freely control the sparsity or the
number of neighbours for each of the nodes. This overcomes the problems
regarding the nodes with few neighbours [20, 19]. Moreover, we can generate a
connected graph for any subset of genes present in the ontology embedding
space, while taking a subset of genes for PPI network may result in a
disconnected graph.

Dutil et al. [17] learn embeddings for the genes during the training of the
graph neural network. The input for the nodes in the model is these embeddings
(initialised at random) scaled by their corresponding expression level for a
particular sample. Such a method may potentially obfuscate the values of the
input features (expressions).

The simplest way to avoid this is to use the expression values directly as the
model input without scaling them via node embeddings. We propose another
method of combining knowledge-based gene embeddings with the expression values
(varying for each sample). First, the ontology embedding is concatenated with
the expression and passed through a fully connected layer. Then we concatenate
the expression value to the produced representation. Such an approach is
potentially better than a simple scaling as it adjusts the node embeddings
depending on the expression, and it also ensures the preservation of the
expression value itself. Having knowledge-driven node embeddings helps
convolution kernels to differentiate better between the genes during the graph
convolutions, while at the same time it captures the similarity between the
genes related by the biological knowledge.

##  3 Experiments and results

The experiments consider the tasks of classifying cancer patients from their
genomic data collected by Molecular Taxonomy of Breast Cancer International
Consortium (METABRIC) [25]. The classifications aim at three different
subtypes of breast cancer: PAM50 (5-class molecular cancer subtypes) [16], ER
(binary classification into immunohistochemistry subtypes) and iC10 (11
IntegrativeCluster subtypes) [25].

For the bulk of the implementation, we reused the codebase from Dutil et al.
[17]. We extended the implementation by custom graph topology generated from
ontology embeddings and novel methods of generating knowledge-driven node
embeddings. As the baselines, we use the Multi-Layer Perceptron (MLP) with
dropout and Random Forest (RF).

Table 1: Performance comparison of the methods on the PAM50, ER and iC10
patient classification tasks (accuracy$\pm$std), using all genes as input
features. Bold font marks the cases when the model obtained statistically
significantly better results than the baselines.

|  | PAM50 | ER | IC10
---|---|---|---|---
| train size | 100 | 1500 | 100 | 1500 | 100 | 1500
OntoGCN w/o node embeddings | 72.3$\pm$2.9 | 81.2$\pm$0.8 | 91.4$\pm$1.0 | 93.7$\pm$1.0 | 48.2$\pm$3.1 | 74.3$\pm$2.3
OntoGCN w/ node embeddings | 72.7$\pm$3.7 | 81.6$\pm$2.2 | 90.8$\pm$0.5 | 93.8$\pm$1.2 | 50.4$\pm 2$.3 | 73.7$\pm$3.3
Random Forest | 69.3$\pm$1.1 | 78.4$\pm$1.0 | 88.7$\pm$1.2 | 93.1$\pm$1.2 | 66.7$\pm$0.9 | 71.3$\pm$2.1
Multi-Layer Perceptron | 60.1$\pm$5.0 | 77.9$\pm$2.5 | 88.8$\pm$1.8 | 90.9$\pm$5.0 | 40.4$\pm$4.5 | 68.9$\pm$1.6

Table 1 presents the results in which all of the 24368 available genes are
used as an input to the model with no feature selection. The evaluation
explores various scenarios of data scarcity, as each experiment is performed
using the training sample of 1500 or 100 patients.

The results for PAM50 cancer subtype classification, presented in Table 1,
demonstrate a statistically significant improvement over the baseline models
not using prior biological knowledge (with a p-value of 0.035 between the
results of GCN and MLP for the training sample of 1500). The improvement
achieved by GCN over baselines is especially visible in the case of the low
data scenario of 100 patient samples. The MLP model is struggling to achieve
high scores with the scarce data, due to the vast spurious connections.
Directing the neural connections via ontology-based graph convolutions seems
to overcome these problems and performs better. On the ER task, the GCN (with
no node embeddings or the proposed ontology-based node embeddings) outperforms
the baseline models in the low training data setting of 100 patients. On the
task of iC10 classification, the GCN outperforms baselines when there is
sufficient training data used. However, it struggles when trained on scarce
data and applied to multi-label classification tasks such as iC10 (with 11
labels).

Table 2: Performance comparison of the graphs directing GCN structure on the
PAM50 classification task (accuracy $\pm$ std). The GCN used expression values
as the nodes’ input, without the node embedding mechanism.

| PAM 50
---|---
| train size=100 | train size=1500
OntoGCN | 72.3$\pm$2.9 | 81.2$\pm$0.8
Random graph | 68.9$\pm$3.9 | 78.0$\pm$4.7
GeneMANIA graph [18] | 71.7$\pm$2.6 | 80.1$\pm$2.1
STRINGdb graph [26] | 71.2$\pm$2.5 | 80.7$\pm$1.2
Random Forest | 69.3$\pm$1.1 | 78.4$\pm$1.0
Multi-Layer Perceptron | 60.1$\pm$5.0 | 77.9$\pm$2.5

In relation to the previous work on imposing bias from biological networks via
graph convolutions, we compare in our experiments the graph generated using
ontology embeddings with two graph datasets containing a mixture of protein-
protein interaction and gene co-expression data: GeneMANIA [18] and STRINGdb
[26]. We also consider a baseline of a graph with randomly generated edges
(with matching degree). Such a baseline allows determining if performance
gains come from the model itself or the underlying prior biological knowledge.
The results in Table 2 suggest that curated biological networks are a good
source of prior knowledge. The slightly improved performance in case of the
automatically generated ontology graph is probably caused by the ability to
freely tune the sparsity of the graph, which regularises the sizes of kernels
in convolution and overcomes the problems caused by sparsely connected genes
in the network [19, 21].

We further combine the proposed GCN model with ontology-based feature
selection strategy. The method combines prior knowledge of known cancer
drivers with automatic selection based on semantic similarity of genes from
the ontology embeddings. It begins by choosing a single gene that is known for
its expression levels being highly correlated with the task. For example, for
the task of predicting PAM50 subtype, such gene is ESR1, which encodes an
estrogen receptor [27, 28, 29]. Given at least one gene expression known as an
important factor in the cancer-related classification task, one can select a
set of (e.g., 1000) genes related with similar functions to the reference gene
according to the cosine distance in the ontology embedding space.

Figure 2: Accuracy obtained on the PAM50 breast cancer subtype classification
tasks. The reported accuracy is an average over 5 trials (with different
stratified split of data) and error bars represent standard error.

The results in Figure 2 for PAM50 classification show that the feature
selection is important when using a smaller set of genes of cardinality below
4000. Given the genes chosen with the proposed feature selection strategy,
both MLP and GCN outperform the baseline of MLP operating on a random set of
genes. The differences are less significant on the larger gene-set covering
most of the genes as the models can then discover signals distributed among a
variety of genes.

The proposed GCN model with structural inductive bias based on ontology prior
knowledge performs on par with MLP on the smaller selected set of genes.
However, it can outperform the MLP in the situation when the whole high-
dimensional input is considered. In such a situation the MLP (green)
performance drops while GCN (yellow) continues to improve with a larger gene-
set when given a large enough training sample. For PAM50 on scarce training
samples, the performance of GCN and MLP reaches its peak around an input set
of 1000 genes. After that MLP performance declines while GCN is still
performing well.

##  4 Conclusions

We presented OntoGCN, a graph neural network model that exploits prior
biological knowledge from ontologies to analyse gene expression data in the
cancer patient classification tasks. The model uses ontology embeddings [5] to
generate a graph directing the convolution operations in the network.

The results suggest that imposing ontology-based knowledge as structural
inductive bias in the model helps to mitigate the difficulties of high-
dimensionality, by limiting the number of connections in the neural network
architecture. This can be useful for clinical tasks when little is known about
the influence of different genes to a particular target bio-marker. In such
scenarios, where all available data has to be put in use, the structural
inductive that limits spurious connections proves beneficial. However, if
there is enough knowledge about the task to do feature selection, then
limiting the number of inputs and applying a simpler machine learning model
may be a better way forward.

While in this paper we focus only on a number of clinical tasks that are
related to known gene expressions, the proposed approach is general and
applicable to other tasks. Therefore, we intend to investigate the generality
using other data from other clinical studies. Moreover, it would be
particularly interesting to evaluate the method on clinical targets derived
from heterogeneous data sources (such as image data as well as other types of
genomic data).

## References

  * [1] Michael Ashburner, Catherine A Ball, Judith A Blake, David Botstein, Heather Butler, J Michael Cherry, Allan P Davis, Kara Dolinski, Selina S Dwight, Janan T Eppig, et al.  Gene ontology: tool for the unification of biology.  Nature genetics, 25(1):25–29, 2000.
  * [2] Maxat Kulmanov, Fatima Zohra Smaili, Xin Gao, and Robert Hoehndorf.  Machine learning with biomedical ontologies.  bioRxiv, 2020.
  * [3] Jake Crawford and Casey S Greene.  Incorporating biological structure into machine learning models in biomedicine.  Current Opinion in Biotechnology, 63:126–134, 2020.
  * [4] Marinka Zitnik, Francis Nguyen, Bo Wang, Jure Leskovec, Anna Goldenberg, and Michael M Hoffman.  Machine learning for integrating data in biology and medicine: Principles, practice, and opportunities.  Information Fusion, 50:71–91, 2019.
  * [5] Jun Chen, Azza Th Althagafi, and Robert Hoehndorf.  Predicting candidate genes from phenotypes, functions, and anatomical site of expression.  Bioinformatics, 2020.
  * [6] Fatima Zohra Smaili, Xin Gao, and Robert Hoehndorf.  Opa2vec: combining formal and informal content of biomedical ontologies to improve similarity-based prediction.  Bioinformatics, 35(12):2133–2140, 2019.
  * [7] Fatima Zohra Smaili, Xin Gao, and Robert Hoehndorf.  Onto2vec: joint vector-based representation of biological entities and their ontology-based annotations.  Bioinformatics, 34(13):52–60, 2018.
  * [8] Maxat Kulmanov, Wang Liu-Wei, Yuan Yan, and Robert Hoehndorf.  El embeddings: geometric construction of models for the description logic el++.  pages 6103–6109, 2019.
  * [9] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps and spectral techniques for embedding and clustering.  In Proceedings of the 14th International Conference on Neural Information Processing Systems: Natural and Synthetic, NeurIPS’01, pages 585–591. MIT Press, 2001.
  * [10] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, NeurIPS’16, pages 3844–3852. Curran Associates Inc., 2016.
  * [11] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In Proceedings of the 5th International Conference on Learning Representations, ICLR’17, 2017.
  * [12] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, pages 1263–1272, 2017.
  * [13] Lenore Cowen, Trey Ideker, Benjamin J. Raphael, and Roded Sharan.  Network propagation: a universal amplifier of genetic associations.  Nature Reviews Genetics, 18(9):551–562, Sep 2017.
  * [14] Sungmin Rhee, Seokjun Seo, and Sun Kim.  Hybrid approach of relation network and localized graph convolutional filtering for breast cancer subtype classification.  In Proceedings of the 27th International Joint Conference on Artificial Intelligence, IJCAI’18, page 3527–3534. AAAI Press, 2018.
  * [15] Adam Santoro, David Raposo, David G Barrett, Mateusz Malinowski, Razvan Pascanu, Peter Battaglia, and Timothy Lillicrap.  A simple neural network module for relational reasoning.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30, pages 4967–4976, 2017.
  * [16] Joel S Parker, Michael Mullins, Maggie CU Cheang, Samuel Leung, David Voduc, Tammi Vickery, Sherri Davies, Christiane Fauron, Xiaping He, Zhiyuan Hu, et al.  Supervised risk predictor of breast cancer based on intrinsic subtypes.  Journal of clinical oncology, 27(8):1160, 2009.
  * [17] Francis Dutil, Joseph Paul Cohen, Martin Weiss, Georgy Derevyanko, and Yoshua Bengio.  Towards gene expression convolutions using gene interaction graphs.  In International Conference on Machine Learning (ICML) Workshop on Computational Biology (WCB), 2018.
  * [18] David Warde-Farley, Sylva L Donaldson, Ovi Comes, Khalid Zuberi, Rashad Badrawi, Pauline Chao, Max Franz, Chris Grouios, Farzana Kazi, Christian Tannus Lopes, et al.  The genemania prediction server: biological network integration for gene prioritization and predicting gene function.  Nucleic acids research, 38(suppl_2):W214–W220, 2010.
  * [19] Paul Bertin, Mohammad Hashir, Martin Weiß, Geneviève Boucher, Vincent Frappier, and Joseph Paul Cohen.  Analysis of gene interaction graphs for biasing machine learning models.  arXiv: Genomics, abs/1905.02295, 2019.
  * [20] Mohammad Hashir, Paul Bertin, Martin Weiss, Vincent Frappier, Theodore Perkins, Geneviève Boucher, and Joseph Paul Cohen.  Is graph biased feature selection of genes better than random?  Machine Learning in Computational Biology (MLCB) meeting, 2019.
  * [21] Jake Crawford and Casey S Greene.  Graph biased feature selection of genes is better than random for many genes.  BioRxiv, 2020.
  * [22] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [23] Christopher J Mungall, Carlo Torniai, Georgios V Gkoutos, Suzanna E Lewis, and Melissa A Haendel.  Uberon, an integrative multi-species anatomy ontology.  Genome biology, 13(1):R5, 2012.
  * [24] Cynthia L Smith and Janan T Eppig.  The mammalian phenotype ontology: enabling robust annotation and comparative analysis.  Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 1(3):390–399, 2009.
  * [25] Christina Curtis, Sohrab P Shah, Suet-Feung Chin, Gulisa Turashvili, Oscar M Rueda, Mark J Dunning, Doug Speed, Andy G Lynch, Shamith Samarajiwa, Yinyin Yuan, et al.  The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups.  Nature, 486(7403):346–352, 2012.
  * [26] Damian Szklarczyk, Annika L Gable, David Lyon, Alexander Junge, Stefan Wyder, Jaime Huerta-Cepas, Milan Simonovic, Nadezhda T Doncheva, John H Morris, Peer Bork, et al.  String v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.  Nucleic acids research, 47(D1):D607–D613, 2019.
  * [27] Roy RL Bastien, Álvaro Rodríguez-Lescure, Mark TW Ebbert, Aleix Prat, Blanca Munárriz, Leslie Rowe, Patricia Miller, Manuel Ruiz-Borrego, Daniel Anderson, Bradley Lyons, et al.  Pam50 breast cancer subtyping by rt-qpcr and concordance with standard clinical molecular markers.  BMC medical genomics, 5(1):44, 2012.
  * [28] Frederik Holst, Phillip R Stahl, Christian Ruiz, Olaf Hellwinkel, Zeenath Jehan, Marc Wendland, Annette Lebeau, Luigi Terracciano, Khawla Al-Kuraya, Fritz Jänicke, et al.  Estrogen receptor alpha (esr1) gene amplification is frequent in breast cancer.  Nature genetics, 39(5):655–660, 2007.
  * [29] Steffi Oesterreich and Nancy E Davidson.  The search for esr1 mutations in breast cancer.  Nature genetics, 45(12):1415–1416, 2013.