#  CH3-Terminated Carbon Chains in the GOTHAM Survey of TMC-1: Evidence of
Interstellar CH3C7N

Mark A. Siebert  Department of Astronomy, University of Virginia,
Charlottesville, VA 22904, USA  Kin Long Kelvin Lee  Department of Chemistry,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Anthony J.
Remijan  National Radio Astronomy Observatory, Charlottesville, VA 22903, USA
Andrew M. Burkhardt  Center for Astrophysics $\mid$ Harvard & Smithsonian,
Cambridge, MA 02138, USA  Department of Physics, Wellesley College, 106
Central Street, Wellesley, MA 02481, U.S.A.  The GOTHAM Collaboration  Ryan A.
Loomis  National Radio Astronomy Observatory, Charlottesville, VA 22903, USA
Michael C. McCarthy  Center for Astrophysics $\mid$ Harvard & Smithsonian,
Cambridge, MA 02138, USA  Brett A. McGuire  Department of Chemistry,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  National
Radio Astronomy Observatory, Charlottesville, VA 22903, USA  Center for
Astrophysics $\mid$ Harvard & Smithsonian, Cambridge, MA 02138, USA Mark A.
Siebert, Brett A. McGuire<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted October 20, 2021)

###### Abstract

We report a systematic study of all known methyl carbon chains toward TMC-1
using the second data release of the GOTHAM survey, as well as a search for
larger species. Using Markov-Chain Monte Carlo simulations and spectral line
stacking of over 30 rotational transitions, we report statistically
significant emission from methylcyanotriacetylene (CH3C7N) at a confidence
level of $4.6\sigma$, and use it to derive a column density of ${\sim}10^{11}$
cm-2 . We also searched for the related species, methyltetraacetylene
(CH3C8H), and place upper limits on the column density of this molecule. By
carrying out the above statistical analyses for all other previously detected
methyl-terminated carbon chains that have emission lines in our survey, we
assess the abundances, excitation conditions, and formation chemistry of
methylpolyynes (CH3C2_nH) and methylcyanopolyynes (CH3C2_n-1N) in TMC-1, and
compare those with predictions from a chemical model. Based on our observed
trends in column density and relative populations of the $A$ and $E$ nuclear
spin isomers, we find that the methylpolyynes and methylcyanopolyynes families
exhibit stark differences from one another, pointing to separate interstellar
formation pathways, which is confirmed through gas-grain chemical modeling
with nautilus.

Astrochemistry, ISM: molecules

††journal: ApJ

##  1 Introduction

One of the most well-studied sites of cold, interstellar carbon chemistry is a
prestellar core within the Taurus Molecular Cloud complex commonly referred to
as Taurus Molecular Cloud 1 (TMC-1) (Churchwell et al., 1978). The chemistry
occurring in TMC-1 produces a suite of exotic molecules (by terrestrial
standards) including unsaturated cyanopolyyne (HC2_n+1N; $n=1,2,3,...$) and
acetylenic (C2_nH2) linear carbon chains. Among the variations on these
abundant molecules are the symmetric top methylpolyynes (MPs) and
methylcyanopolyynes (MCPs), which have the form CH3C2_nH and CH3C2_n-1N,
respectively. These species have been known in the interstellar medium (ISM)
since the discovery of methyl cyanide (CH3CN) (Solomon et al., 1971) and
methyl acetylene (CH3CCH) (Irvine et al., 1981). Toward TMC-1, methylpolyynes
and methylcyanopolyynes as large as CH3C6H and CH3C5N have been detected and
characterized (Matthews & Sears, 1983; Irvine et al., 1981; Broten et al.,
1984; Walmsley et al., 1984; MacLeod et al., 1984; Loren et al., 1984; Snyder
et al., 2006; Remijan et al., 2006).

Currently, the chemical formation of MPs and MCPs in TMC-1 is unconstrained.
While it is possible that simple carbon addition reactions occurring on the
surfaces of dust grains may account for the growth of these species toward
cold, dark clouds (Turner et al., 2000), it is also plausible that ion-
molecule reaction pathways in the gas phase can account for their formation
(Quan & Herbst, 2007). Because molecules with C3v symmetry have unique
excitation properties that do not allow for easy interconversion between $A$
and $E$ symmetry states, the relative populations of their $K$-ladders can
offer additional clues to their formation conditions (see e.g. Minh et al.
(1993) Askne et al. (1984), Willacy et al. (1993), (Mendoza et al., 2018)),
but to date, analysis of these populations has only been done for the shortest
MPs and MCPs. Furthermore, comparing the relative column densities of the
similarly structured cyanopolyynes has proven to be a useful tool in testing
chemical models and assessing the mechanisms through which carbon chain
species grow in the ISM (Loomis et al., 2021; Agúndez et al., 2017; Burkhardt
et al., 2018).

Obtaining a more complete understanding of their chemistry and setting
constraints on the abundances and excitation conditions of the various methyl-
terminated carbon chains requires exploring the extent to which these types of
species can grow (i.e. how long of a C-chain can form). Such observations
would enable the direct comparison of measured abundances with both grain-
surface and gas-phase reaction pathways, especially when combined with a fully
self-consistent excitation analysis of MPs and MCPs in TMC-1. In addition,
these studies will set the overall limit to the detectability of longer chain
MPs and MCPs based on both an improved understanding of the chemical formation
pathway and on a new computational methodology to detect the weak signals
coming from the larger C-chain species.

The ongoing GOTHAM survey (GBT Observations of TMC-1: Hunting for Aromatic
Molecules) along with a deep Q-band survey being conducted on the Yebes 40 m
radio telescope (Cernicharo et al., 2020) have thus far proven to be very
useful for rigorous studies of exotic molecules in TMC-1. With many new
detections of various carbon-bearing molecules, their isomers, and
cyclic/polycyclic species, these projects have effectively expanded our
knowledge of the interstellar molecular carbon reservoir before star formation
occurs (Xue et al., 2020; Loomis et al., 2021; McGuire et al., 2018; McCarthy
et al., 2020; McGuire et al., 2020, 2021; Cabezas et al., 2021; Cernicharo et
al., 2021). In this work, we present a detailed analysis of CH3-terminated
carbon chains toward TMC-1 as well as the first detection of
methylcyanotriacetylene (CH3C7N, Chen et al. (1998)) using the second data
release of GOTHAM. We also search for methyltetraacetylene (CH3C8H, Travers et
al. (1998a)) in our data, but do not find significant emission and thus report
upper limits on its abundance.

In Section 2, we outline the specifics of our observations and reduction
methods. In Section 3 we describe our spectroscopic calculations for the
targeted molecules and statistical analysis procedures. In Section 4 we
present our results for all targeted molecules and discuss their relative
column densities and excitation physics. Finally, in Section 5 we discuss how
our results fit in to the current understanding of chemistry in TMC-1, and
compare them with predictions of our chemical model.

##  2 Observations

Observations obtained in this study were collected as part of the GOTHAM
Survey. GOTHAM is a large project on the 100m Green Bank Telescope (GBT)
currently carrying out dedicated spectral line observations of TMC-1 covering
almost 30 GHz of radio bandwidth at high sensitivity and ultra-fine spectral
resolution. This work uses the second data release (DR2) of GOTHAM, which
includes over 600 hours of observations targeting the cyanopolyyne peak (CP)
of TMC-1, centered at $\alpha_{\mathrm{J2000}}$ = 04$\fh$41$\fm$42.5$\fs$,
$\delta_{\mathrm{J2000}}$ = +25°41′26.8″. A full description of DR2
specifications and our reduction pipeline is described (McGuire et al., 2020),
but put briefly, the spectra in this data set cover the entirety of the X-,
K-, and Ka-receiver bands with nearly continuous coverage from 8.0 to 11.6 GHz
and 18.0 to 33.5 GHz (25.6 GHz of total bandwidth). All spectra have a uniform
frequency resolution of 1.4 kHz (0.05–0.01 km/s in velocity) and an RMS noise
of ${\sim}2$–20 mK, with the RMS gradually increasing toward higher frequency
because of the lower integration times. Data reduction involved removal of RFI
and artifacts, baseline continuum fitting, and flux calibration using
complementary VLA observations of the source J0530+1331. Uncertainty from to
this flux calibration is estimated at $\sim$20%, and is factored into our
statistical analysis described below (McGuire et al., 2020).

##  3 Methods & Analysis

###  3.1 Molecular spectroscopy

Table 1: Spectroscopic properties of stacked CH3C7N lines covered by the GBT
X-band receiver

Transitions | Symm. | Frequency | $E_{up}$ | $S_{ij}\mu^{2}$
---|---|---|---|---
$J^{\prime}\ \rightarrow\ J^{\prime\prime}$ | $K$ | $F^{\prime}\ \rightarrow\ F^{\prime\prime}$ |  | (MHz) | (K) | (D2)
$11\ \rightarrow\ 10$ | $1$ | $10\ \rightarrow\ 9$ | $E$ | 8243.824 | 9.92 | 355.325
|  | $11\ \rightarrow\ 10$ | $E$ | 8243.824 | 9.92 | 389.482
|  | $12\ \rightarrow\ 11$ | $E$ | 8243.836 | 9.92 | 426.877
| $0$ | $10\ \rightarrow\ 9$ | $A$ | 8243.851 | 2.37 | 358.286
|  | $11\ \rightarrow\ 10$ | $A$ | 8243.859 | 2.37 | 392.727
|  | $12\ \rightarrow\ 11$ | $A$ | 8243.865 | 2.37 | 430.435
$12\ \rightarrow\ 11$ | $1$ | $11\ \rightarrow\ 10$ | $E$ | 8993.262 | 10.35 | 391.696
|  | $12\ \rightarrow\ 11$ | $E$ | 8993.263 | 10.35 | 426.021
|  | $13\ \rightarrow\ 12$ | $E$ | 8993.272 | 10.35 | 463.320
| $0$ | $11\ \rightarrow\ 10$ | $A$ | 8993.293 | 2.81 | 394.435
|  | $12\ \rightarrow\ 11$ | $A$ | 8993.299 | 2.81 | 429.000
|  | $13\ \rightarrow\ 12$ | $A$ | 8993.304 | 2.81 | 466.560
$13\ \rightarrow\ 12$ | $1$ | $12\ \rightarrow\ 11$ | $E$ | 9742.699 | 10.82 | 428.012
|  | $13\ \rightarrow\ 12$ | $E$ | 9742.700 | 10.82 | 462.478
|  | $14\ \rightarrow\ 13$ | $E$ | 9742.708 | 10.82 | 499.692
| $0$ | $12\ \rightarrow\ 11$ | $A$ | 9742.733 | 3.27 | 430.560
|  | $13\ \rightarrow\ 12$ | $A$ | 9742.739 | 3.27 | 465.231
|  | $14\ \rightarrow\ 13$ | $A$ | 9742.743 | 3.27 | 502.666
$14\ \rightarrow\ 13$ | $1$ | $13\ \rightarrow\ 12$ | $E$ | 10492.136 | 11.33 | 464.286
|  | $14\ \rightarrow\ 13$ | $E$ | 10492.137 | 11.33 | 498.870
|  | $15\ \rightarrow\ 14$ | $E$ | 10492.144 | 11.33 | 536.010
| $0$ | $13\ \rightarrow\ 12$ | $A$ | 10492.173 | 3.78 | 466.667
|  | $14\ \rightarrow\ 13$ | $A$ | 10492.178 | 3.78 | 501.429
|  | $15\ \rightarrow\ 14$ | $A$ | 10492.181 | 3.78 | 538.758
$15\ \rightarrow\ 14$ | $1$ | $14\ \rightarrow\ 13$ | $E$ | 11241.572 | 11.86 | 500.524
|  | $15\ \rightarrow\ 14$ | $E$ | 11241.573 | 11.86 | 535.211
|  | $16\ \rightarrow\ 15$ | $E$ | 11241.579 | 11.86 | 572.284
| $0$ | $14\ \rightarrow\ 13$ | $A$ | 11241.612 | 4.32 | 502.759
|  | $15\ \rightarrow\ 14$ | $A$ | 11241.616 | 4.32 | 537.600
|  | $16\ \rightarrow\ 15$ | $A$ | 11241.619 | 4.32 | 574.839

Note. — Only $K=0$ and $K=1$ are included here since they contribute the
brightest transitions; however we still considered states as high as $K=13$ in
our MCMC model and stacking procedure. Similarly, lines ranging from
$J^{\prime}=27$ to $J^{\prime}=44$ are covered by the GBT K- and Ka-bands and
included in our analysis, but not shown here. For the full list of transitions
used for CH3C7N and all other methyl carbon chains treated in this work,
please refer to the linked Dataverse repository,10.7910/DVN/K9HRCK (catalog
10.7910/DVN/K9HRCK).

All the MPs and MCPs are prolate symmetric top molecules with C3v symmetry,
meaning their rotational energy states are defined by the total angular
momentum state $J$ in addition to its projection along the unique axis of
rotation, denoted by the quantum number $K$. Depending on the value of $J$ and
$K$, symmetric top molecules will either exhibit $A$ or $E$ symmetry due to
the nuclear spin on the methyl group; and since radiative and collisional
transitions between these symmetry states are strictly limited, their relative
populations can be far from the expected thermal distribution. Because of
this, we treat transitions from non-degenerate ($A_{1},A_{2}$ symmetry) and
doubly degenerate ($E$) levels separately by running independent MCMC fits. As
such, we can compare the relative abundances of carbon chains in these energy
states, similar to the separate treatment of the $K=0$ and $K=1$ components
employed by Snyder et al. (2006) and Remijan et al. (2006).

For molecules with C3v symmetry, the total statistical weight of rotational
energy levels is the product of the typical $J$-degeneracy $g_{J}=2J+1$, the
$K$-level degeneracy ($g_{k}=1$ for $K=0$ and $g_{k}=2$ for $K\neq 0$), and
the nuclear spin degeneracy $g_{I}$ (Gordy & Cook, 1984). The reduced
statistical weight contributed by nuclear spin statistics of three identical
hydrogen nuclei in a methyl group (each with spin $I=1/2$) can be computed
with the following:

$g_{I}=\frac{1}{3}\frac{\left(4I^{2}+4I+3\right)}{\left(2I+1\right)^{2}}\mathrm{,}\quad\mathrm{for}\quad
K=0,3,6,9,...$ (1)
$g_{I}=\frac{1}{3}\frac{\left(4I^{2}+4I\right)}{\left(2I+1\right)^{2}}\mathrm{,}\quad\mathrm{for}\quad
K\neq 0,3,6,9,...$ (2)

After combining all degeneracies, the states with $A$ symmetry are weighted
2:1 relative to those with $E$ symmetry, except in the case of $K=0$ which is
weighted equally because of the lower factor of $g_{k}$. The rotational
constants for the largest MCPs and MPs were measured by Chen et al. (1998) and
Travers et al. (1998b). To generate the spectroscopic data used for this work
from these measurements, we used PGopher (Western, 2017; Western &
Billinghurst, 2019), which is able to account for the necessary symmetry and
statistical considerations. The input files and line lists for all six species
are available on the Harvard DataVerse10.7910/DVN/K9HRCK (catalog
10.7910/DVN/K9HRCK).

As an example, Table 1 summarizes the spectroscopic properties of the
rotational transitions of CH3C7N in X-band that were a part of this
investigation. The transition quantum numbers, calculated rest frequencies
(MHz), upper state energy level ($K$), transition line strengths S($J,K$) are
presented. In the case of the MCPs, the nuclear spin on the nitrogen nucleus
contributes hyperfine splitting to each rotational energy level.

Methyl acetylene (CH3CCH) was not included in our analysis as it does not have
any transitions covered by the GOTHAM survey. Similarly, although the
fundamental transition of methyl cyanide (CH3CN) is included in our spectra, a
lack of additional lines prevents us from performing a full characterization
of this molecule. For both CH3CN and CH3CCH, we refer to previous works
studying these molecules in TMC-1 (e.g. Askne et al. (1984); Minh et al.
(1993)), and instead base our analysis on the longer species.

###  3.2 MCMC modeling

In order to derive physical characteristics for the CH3- polyyne and
cyanopolyyne species in TMC-1, we use the same Markov-Chain Monte Carlo (MCMC)
model employed in previous publications from the GOTHAM collaboration McGuire
et al. (2020); Xue et al. (2020); McCarthy et al. (2020). This procedure is
discussed at length in Loomis et al. (2021), but we will summarize it here as
well. In short, the MCMC model calculates probability distributions and
covariances for parameters describing the physical and excitation conditions
of a molecule in TMC-1. Based on recent observations performed with the 45 m
telescope at Nobeyama Radio Observatory (Dobashi et al., 2018, 2019), as well
as our data (Loomis et al., 2021), emission from molecules in TMC-1 toward the
cyanopolyyne peak display at least four individual velocity components. In our
current model, we make the assumption that the emission from each velocity
component is cospatial, similar to the approach adopted in Loomis et al.
(2021) for the cyanopolyynes. To summarize, the model describing a given
methyl chain comprises a source size, four radial velocities, eight column
densities and two excitation temperatures describing A/E components, and a
linewidth parameter, giving a total of 16 parameters.

Line profile simulations were performed using molsim (Lee & McGuire, 2020).
The MCMC simulations used wrapper functions in molsim to arviz (Kumar et al.,
2019) and emcee (Foreman-Mackey et al., 2013); the former for analyzing the
results of sampling, and the latter implements an affine-invariant MCMC
sampler. As prior parameters, we used the marginalized posteriors from the
corresponding cyanopolyyne chain (e.g. HC3N and CH3C3N). The prior
distributions are approximated as normally distributed (i.e. $p(\theta)\sim
N(\mu_{\theta},\sigma_{\theta})$ for parameter $\theta$) with modifications to
the variance $\sigma_{\theta}$ as to avoid overly constrictive/influential
priors. Convergence of the MCMC was confirmed using standard diagnostics such
as the Gelman & Rubin (1992) $\hat{R}$ statistic, and by visually inspecting
the posterior traces. With each of the six species studied in this work, we
compute separate MCMC model fits for the $A$ and $E$ symmetry states so their
relative populations could be derived directly through comparison of the
posterior distributions.

###  3.3 Spectral line stacking and matched filter

Following the formalism of Loomis et al. (2018, 2021), we perform a combined
velocity stack and matched filter analysis in order to verify that the results
of the MCMC model are consistent with the data, and determine the statistical
significance of molecular detections. The first step involves a noise-weighted
sum in velocity space of all expected transitions in our spectra, excluding
any interloping transitions from known molecules in TMC-1; however, this was
not necessary for any of the transitions in this analysis as there were no
such interlopers. Next, the resulting composite spectrum is passed through a
matched filter, using the best-fit parameters of the MCMC to create a model
stacked spectrum. The impulse response of the matched filter, measured in
signal-to-noise ratio $\sigma$, is a representation of how well the MCMC model
reproduces GOTHAM spectra. Generally, we adopt a lower limit of $5\sigma$ as a
threshold for confirmation of a molecule in TMC-1.

##  4 Results & Discussion

###  4.1 Column densities and physical characteristics

Figure 1: (First and third rows) Velocity stacks of the observational (blue) and MCMC simulated (red) spectra. The cross-correlation of the two velocity stacks for each molecule corresponds to the matched filter (MF) spectrum (second and fourth rows), which are shown in the second and fourth rows. Though the transitions belonging to $A$ and $E$ symmetry states were treated with separate MCMC fits, their results are combined here. The peak impulse response are annotated on the MF spectra. Table 2: MCMC model results for all observed CH3-terminated carbon chains. Uncertainties are listed as the 95% confidence level.  Molecule | $N_{T}$ | $N_{A}$ | $N_{E}$ |  $T_{ex}$ ($A$ sym.) |  $T_{ex}$ ($E$ sym.) | Source Size |  $A/E$ Ratio
---|---|---|---|---|---|---|---
- | ($10^{11}$ cm-2) | ($10^{11}$ cm-2) | ($10^{11}$ cm-2) | (K) | (K) | $"$ | -
CH3C3N | 8.66 (.46) | 5.00 (.40) | 3.66 (.24) | 3.36 (.40) | 5.70 (.26) | 481.0 (8.7) aaHC3N source size adopted as prior. See Loomis et al. (2021) for more details. | 1.37 (.15)
CH3C5N | 2.86 (.30) | 2.01 (.29) | 0.85 (.07) | 3.48 (.34) | 3.97 (.48) | 128.7 (4.9) bbHC5N source size adopted as prior. | 2.40 (.40)
CH3C7N | 0.86 (.19) | 0.69 (.18) | 0.17 (.06) | 3.82 (.25) | 3.90 (.44) | 53.9 (1.0) ccHC7N source size adopted as prior. | 4.20 (1.8)
CH3C4H | 100.8 (5.7) | 60.0 (5.0) | 40.8 (2.8) | 4.31 (.84) | 5.50 (.37) | 481.0 (8.6) aaHC3N source size adopted as prior. See Loomis et al. (2021) for more details. | 1.46 (.15)
CH3C6H | 10.4 (.72) | 6.10 (.60) | 4.30 (.40) | 7.19 (1.2) | 6.41 (1.1) | 128.3 (4.8) bbHC5N source size adopted as prior. | 1.41 (.20)
CH3C8H | $<0.98$ | $<0.8$ | $<0.18$ | - | - | - | -

Figure 1 shows the stacked data, the stacked MCMC model, as well as the
matched filter response of all the methyl carbon chains considered in this
work. We clearly detect CH3C3N, CH3C5N, CH3C4H, and CH3C6H at high
significance, and the model is able to fit all components (both velocity and
hyperfine) of these molecules. In the case of CH3C7N, though no individual
lines are present above the current noise level of the survey, the stacked
emission in the top right of Figure 1 exhibits noteworthy signal. Furthermore,
its matched filter has a central peak at 4.6$\sigma$, indicating evidence that
this molecule is present in our data, and that the MCMC model converged to a
set of parameters that reproduce its emission. The matched filter response is
weaker than previous molecules that were found using this same method (Loomis
et al., 2021; Lee et al., 2021; McGuire et al., 2021) and falls just below our
desired 5.0$\sigma$ threshold for a definitive discovery, but it is sufficient
to consider a tentative detection and adopt the parameters of its fit in our
analysis. For more information on this detection as it is indicated in the
Bayesian fit to the spectra (prior to any stacking of emission lines), in
Appendix A, we present corner plots illustrating the full results of the MCMC
model for both CH3C7N and CH3C8H, and discuss the posterior distributions and
covariances in the parameter spaces for each. In contrast to CH3C7N, CH3C8H
shows a much higher degree of uncertainty in its fit, in addition to no signal
in the stack nor its matched filter response (bottom-right panels of Fig. 1),
so we therefore place upper limits on its column density.

A summary of parameters derived from the posterior distributions is provided
in Table 2. Here, we note that the excitation temperatures of MPs (4.3 – 7.2
K) are systematically higher than those of the MCPs (3.4 – 3.9 K), which is in
agreement with previous studies of these species in TMC-1 (Broten et al.,
1984; Askne et al., 1984; Snyder et al., 2006). However, one exception to this
is seen for the $E$ state of CH3C3N, which has $T_{ex}{\sim}5.7$ K. At all
derived excitation temperatures, we find that less than 4% of the total
integrated line flux is contributed by transitions with $K>1$. Furthermore,
the population of the $K=2$ state at these low temperatures is less than 1%
for all detected chains, which is well within the uncertainty of our measured
column densities. In that sense, although states up to $K=13$ were considered,
the column densities of the $A$ and $E$ components reported in Table 2 can be
treated as equivalent to the column densities of the $K=0$ and $K=1$
components, respectively.

The total column densities (combining both $A$ and $E$ symmetries) decrease
with increasing carbon chain length, much like the similar linear
cyanopolyynes (Loomis et al., 2021). This is shown pictorially in Figure 2,
with a comparison to the column densities found by Remijan et al. (2006).
There is small but notable disagreement (less than a factor of four) between
the column densities we find and those derived in Remijan et al. (2006). We
attribute this discrepancy to the intrinsic line strengths ($S_{ij}$) adopted
from our numerical treatment of these molecules (see Section 3.1) as they are
factors of 2–3 larger than the analytic values employed in Remijan et al.
(2006) and Snyder et al. (2006). Another important factor is that the column
densities we derive take into account decreasing source sizes, whereas
previous studies of TMC-1 assume no effects due to beam dilution.

CH3C7N has a total best-fit column density of $8.60\pm 1.9\times 10^{10}$
cm-2, which matches well with the log-linear slope set by CH3C3N and CH3C5N
(Fig. 2). Furthermore, this slope is very similar to that observed for the
cyanopolyynes, and would imply the next chain (CH3C9N) to have a column
density of about $3\times 10^{10}$ cm-2. To detect stacked emission from this
molecule at such an abundance, the RMS noise of the survey would need to be
reduced by at least factor of five. This would require several hundred more
hours of integration time, primarily at lower frequencies (e.g. the GBT X-band
receiver) because the brightest transitions of this molecule are simulated to
be present in the 8–10 GHz range.

In contrast to the cyanopolyynes and MCPs, the MPs decrease in abundance at a
much faster rate than the other carbon chains in TMC-1 (about 1 dex for each
subsequent species). Furthermore, the upper limit we place on CH3C8H appears
to break the trend seen in the shorter methylpolyynes. A linear extrapolation
on the column densities of CH3C4H and CH3C6H (the blue points in Fig. 2) would
imply the next chain length has a column density ${\sim}2\times 10^{11}$ cm-2,
but in this case we would have detected emission from CH3C8H at the
sensitivity of our survey. Instead, it appears that the methylpolyynes
experience a sharp drop-off in column density after CH3C6H, which indicates
that the production of these carbon chains is severely hindered for these
longer species.

A comparison to analogous linear species in TMC-1 is more difficult for the
MPs, because the acetylene chains H2C2n have no dipole moment, and thus no
rotational transition lines with which we could apply a similar analysis.
However, the hydrocarbon radicals (C2nH) do exhibit transitions in the
radio/sub-mm and have been characterized in TMC-1 up to a chain length of
eight by Brünken et al. (2007) (green points in Fig. 2). Here, it is apparent
that those species also exhibit a non-linear decrease in column density after
a carbon chain length of six. In other words, a similar ratio of
$\frac{[\mathrm{\ce{CH3C8H}}]}{[\ce{CH3C6H}]}$ to what is observed for
$\frac{[\mathrm{\ce{C8H}}]}{[\ce{C6H}]}$ would be in agreement with the upper
limits we set in this work. This may suggest that similar chemical routes
govern the abundances of both the MPs and the hydrocarbon radicals. However,
since our upper limit on CH3C8H is still very close to the log-linear
extrapolated column density, it will be imperative to continue studying this
molecule with future more sensitive data releases from GOTHAM survey to
identify how steep this drop-off is, and how similar it is to the relative
column density of C8H.

Figure 2: Derived total column densities of methyl-terminated carbon chains in
TMC-1 as a function of chain length. Methylcyanopolyynes are shown in red
while methylpolyynes are shown in blue. Contributions from the A and E
symmetries of each species are summed, with their errors added in quadrature.
CH3C8H is denoted with a triangle as it is an upper limit. Results from
Remijan et al. (2006) are overplotted in their respective colors with no
connecting line. Additionally, column densities of the linear cyanopolyynes
and hydrocarbons radicals published in Loomis et al. (2021) and Brünken et al.
(2007) are shown in pink and green, respectively.

###  4.2 A/E symmetry state populations

Figure 3: Column density ratios toward TMC-1 of $A$ and $E$ symmetry states
for each of the chains studied in this work (colored points), as well as those
derived for methyl cyanide and methyl acetylene by Minh et al. (1993) and
Askne et al. (1984), respectively (black points). Also plotted in black is the
expected ratio of these two populations if the molecules are thermalized to a
Boltzmann distribution at a range of kinetic temperatures. The vertical error
on the points represents the 95% confidence level from our MCMC modeling
results, and the horizontal error bars denote the range of kinetic
temperatures observed for gas phase molecules toward the CP region TMC-1
(Fehér et al., 2016).

Askne et al. (1984) found that in TMC-1, methyl acetylene (CH3CCH) has about
equal abundances of $A$ and $E$ symmetry states, which is not expected under
LTE conditions and could indicate this molecule formed at a higher
temperature. A similar result was found for CH3C4H by Walmsley et al. (1984).
In contrast, Minh et al. (1993) found that for CH3CN in TMC-1,
$N_{E}/N_{A}=0.76\pm 0.09$, which suggests that these states are equilibrated
to the kinetic temperature. With the sensitivity and bandwidth of GOTHAM DR2,
we are able to revisit some of these population studies with a larger data set
and expand them to the longer molecules CH3C6H, CH3C3N, CH3C5N, and CH3C7N.
Because we performed separate MCMC analyses for the $A$ and $E$ states of each
symmetric top carbon chain, their ratios can be computed directly and are
shown in the last column of Table 2.

Figure 3 shows a comparison of all $A/E$ ratios derived in this work plotted
alongside the expected ratio in a thermalized distribution for a range of
excitation temperatures. The gas phase kinetic temperature at the CP of TMC-1
has recently been placed at $11\pm 1.0$ K by Fehér et al. (2016), so methyl
chains with $A/E\approx 1.75$ can be considered as having equilibrated to this
temperature.

For the MCPs (left plot in Fig. 3), we note that the $A/E$ ratio increases for
longer chain lengths, starting at a superthermal value of 1.33 for CH3CN (as
measured by Minh et al. (1993)) and reaching $4.2\pm 1.8$ for CH3C7N which
corresponds to an excitation temperature of $6.1\pm 2$ K, which is much lower
than the kinetic temperature in TMC-1. In contrast, the MPs (right plot in
Fig. 3) have nearly constant $A/E$ ratios between 1.4 and 1.5. While this is
only slightly smaller than the equilibrium value of 1.75, it is indicative
that all MPs studied here have slightly larger $E$ symmetry populations than
what is expected given the physical conditions in TMC-1.

As noted in Section 3.1, both radiative and collisional interconversion
between the $A$ and $E$ symmetry states are forbidden for both the MCPs and
MPs (unless proton exchange collisions are non-negligible). Willacy et al.
(1993) investigated the possibility that $E\rightarrow A$ conversion of CH3CN
occurs through $E$-state species adsorbing to and subsequently desorbing from
grain surfaces in TMC-1. They found that this is an efficient process so long
as the desorption rate is high enough. This may suggest that an important
distinction may be happening in TMC-1. Highly exothermic reactions in the gas
phase (e.g. dissociative recombination) may drive the formation of the MPs,
forcing them into superthermal $A/E$ distributions, where as the largest MCPs
may preferentially form on (or stick to) the surfaces of dust grains. The non-
thermal desorption of large complex molecules like the MCPs in TMC-1 has been
a topic of numerous recent computation efforts (Garrod et al., 2007; Herbst &
Cuppen, 2006; Hoang & Tung, 2019; Minissale et al., 2016; Shingledecker et
al., 2021), but it is still difficult to say exactly how efficient this
process is for the species explored here. The investigation of the potential
bifurcation of the formation pathways between the MPs and the MCPs using
astrochemical models might shed further light on this intriguing situation.

##  5 Chemical modeling of methyl chains

To further explore the chemistry of methyl-terminated carbon chains in TMC-1,
we utilized the adapted three-phase gas-grain chemical network model nautilus
v1.1 code (Ruaud et al., 2016) discussed in previous analysis of GOTHAM data,
which has previously been used to successfully study the formation of carbon-
chain molecules (Xue et al., 2020; McGuire et al., 2020; Shingledecker et al.,
2021). The physical conditions of the model are equal to what has been
previously used with this network ($T_{\mathrm{gas}}=T_{\mathrm{grain}}=10$ K,
$n_{\text{H}_{2}}$=2$\times$104 cm-3, $A_{V}$=10, and
$\zeta_{\text{CR}}$=1.3$\times$10-17 s-1; Hincelin et al. (2011)) as are the
elemental abundances (Loomis et al., 2021). Based originally off of the KIDA
network, the network already contained the formation of MCPs up to $n=7$ and
MPs up to $n=6$. Here, we expanded on this to include the formation of CH3C8H,
as well as additional reaction to better constrain the formation of the
smaller members of this family.

To simulate the formation and destruction of CH3C8H, we utilized the analogous
pathways for the existing network of CH3C6H, whose primary formation routes
are from the dissociate recombination of C7H5+. As such, we expanded out the
formation of C9H5+ from ion-neutral reactions with semi-saturated carbon
chains, whose rates were estimated using the Langevin formula (Woon & Herbst,
2009). It should be noted that the rates for C7H5+ and CH3C6H are calculated
using the Modified Arrhenius formula, which may result in some potential
discrepancies between this and CH3C8H. As a source of destruction of C9H5+, we
adapted existing dissociate recombination of similar cations. The rates and
branching ratios C9H5+ and C7H5+ were estimated and updated, respectively,
based on the dissociate recombination of C5H5+ whose rates are taken from
Herbst & Leung (1989) and branching ratios are discussed in detail on the KIDA
website111http://kida.astrophy.u-bordeaux.fr/datasheet/datasheet_2761_C5H5++e-_V1.pdf.
In addition to dissociative recombination, CnH5+ ions present in KIDA were
also destroyed by reactions with anions. As such, the rates of the anion
destruction of C9H5+ were based off analogous rates estimated for C7H5+ by
Harada & Herbst (2008). Due to a lack of robust isomerization for large
molecules within existing chemical networks, we only consider the C7H5+ and
C9H5+ isomers with KIDA entries (i.e. CH3CnCH2+)

The analogous destruction of CH3C6H was also adapted for CH3C8H. This included
1) dissociation due to photons and cosmic rays 2) ion-neutral reactions with
abundant ions (i.e. H+, H3+, C+, HCO+, and He+) which were estimated with
Langevin formula and 3) reaction with elemental carbon to form C10H2 and H2.
To properly study C9H5+, we also added in several neutral-neutral reactions to
better account for the formation of C5H3, a C9H5+-precursor, which were
obtained originally from Hébrard et al. (2009). Finally, we included an
additional formation pathway for the MPs by reactions of the C2nH family with
methane based on work by Quan & Herbst (2007) and formation of CH3C4H from
CH3CHCH2 by Berteloite et al. (2010).

For the MCPs, no additional reactions or rates were added from the KIDA
network beyond the related adaptations from other recent molecular detections
(e.g. the expanded semi-saturated carbon-chain networks for HC4NC (Xue et al.,
2020), C6H5CN (Burkhardt et al., 2021a), H2CCCHC3N (Shingledecker et al.,
2021), and indene (Burkhardt et al., 2021b)).

Figure 4: Simulated abundances and column densities of the four carbon chain
families studied here at a model time of 2.5$\times$105 years, with the same
colors and markers as Figure 2. Observed values are shown as dotted lines and
the CH3C8H upper limit is denoted as an unfilled triangle.

The results of this model can be seen in Figure 4 at a time of 2.5$\times$105
yrs, assuming a TMC-1 hydrogen column density of $N_{\text{H}_{2}}$=1022 cm-2.
At this time in the model, the simulated abundances of all studied species are
within an order of magnitude of their observed values. For the MCPs (solid red
line in Fig. 4), the abundance of CH3C3N is higher than the observed value,
while the remainder of the family are somewhat below their observed value. As
such, the log-linear trend is less well constrained here. The relative trend
line was also fairly independent of time, with a slight increase of the longer
chains at later times, as seen in Fig. A3. This is expected as longer chains
typically have their peak abundance several 105 years after similar smaller
chains, as seen by the cyanopolyynes (Loomis et al., 2021). The dominant
formation pathways for the MCPs are the dissociative recombination of the
Cn+1H4N+ ion family, which are in turn primarily produced by reactions between
the cyanopolyynes and CH3+, with minor contributions from Cn+1H3N+ and N +
Cn+1H5+ after several 105 years. In particular, CH3C3N has an additional
production from H3C4N+ + H2 that may account for it deviation from the
trendline of the rest of the MCP family. The MCPs are primarily destroyed by
ion-neutral reactions with abundant ions (i.e. H+, H3+, C+, HCO+, and He+).
The dependence on the cyanopolyynes (solid pink line in Fig. 4) is consistent
with the similar observed abundance trends seen in Figure 2.

For the MPs (solid blue line in Fig. 4), the relative abundance trends of both
the simulated and observed values are in good agreement at this time. This
slope of this trend is also very consistent with the simulated abundance trend
of the CnH family (solid green line in Fig. 4). Within the general
uncertainties of the source age of TMC-1 ($\sim$2-5$\times$105 years), the
abundance of CH3C4H was found to be consistently in worst agreement with the
observed values but improving at later times by a $\sim$60% abundance
increase, while the longer chains increased up to nearly an order of magnitude
(See Fig. A3). As a result, the trend line was in less consistently in
agreement with the observed trend at later times, while reproducing the
observed abundances better. Due to the overall depletion of CH3C4H, we have
chosen this as a representative time as the simulated abundances are still
within typical uncertainties for kinetic chemical models and we suspect the
trend at this model time would remain consistent if a more efficient
production is found either through new proposed pathways or improved reaction
rate measurement/calculations. More generally, this strong dependence on the
source age contrasts with other carbon chain families such as the
cyanopolyynes and MCPs, whose observed trends are fairly consistent across the
possible source ages of TMC-1.

The dominant formation route for the MPs are from the dissociative
recombination of the Cn+1H5+ family, with a significant contribution from CnH
reactions with methane after $\sim$105 years. The Cn+1H5+ family are produce
by various ion neutral reactions, notably CnH2+ with methane and Cn+1H2 with
C2H4+ ($<10^{5}$ years). These dominant pathways provide strength for the
observed relation between the MPs and CnH families. There is also an increased
importance after 105 years of C4H3+ + methane for CH3C4H and C3H3 + C4H2 for
CH3C6H. MPs are mostly destroyed by ion-neutral reactions and reactions with
carbon atoms.

Overall, the predictions of these chemical models are within an order of
magnitude of agreement with the observed values of TMC-1, indicating that
carbon chain chemistry can be fairly well understood. In addition, the
observed MP trend line can be reproduced well and is also analogous to what is
observed for CnH. The estimated branching ratios and reaction rates still
remain a major source of uncertainties in the models, which provide a strong
motivation for further laboratory and theoretical studies. As discussed in
Burkhardt et al. (2021b), detailed study of increasing larger molecules (e.g.
molecules with 7+ heavy atoms and/or aromatic rings) will likely require a
rigorous study of the formation of semi-saturated carbon chains that have
historically been performed piecewise. The recent expanded inventory of known
carbon chains in TMC-1 will provide key constraints in this study.

##  6 Conclusions

We performed a rigorous, self-consistent study of methylpolyynes and
methylcyanopolyynes in TMC-1 using the second data release of GOTHAM. Through
MCMC modeling and matched filter analysis of stacked of emission lines, we
derived column densities, excitation temperatures, and $A/E$ symmetry ratios
for all previously found species, and discovered evidence for a new
interstellar symmetric top, methylcyanotriacetylene (CH3C7N), at a confidence
level of 4.6$\sigma$. We also searched for methyltetraacetylene CH3C8H) in our
spectra and place upper limits on its column density in TMC-1.

In our analysis, we found two important divisions between the different
families of molecules explored in this work:

  1. 1.

The column densities of methylcyanopolyynes decrease in a log-linear manner
with increasing carbon chain length, and the slope of this trend is similar to
the cyanopolyynes. In contrast, the methylpolyynes experience a drop-off in
column density for species larger than CH3C6H that is not consistent with the
trend set by the smaller species.

  2. 2.

The $A/E$ ratios of methylcyanopolyynes increase with carbon chain length, and
are subthermal for the larger species. The detected methylpolyynes have
systematically smaller $A/E$ ratios which are not equilibrated to the kinetic
temperature in TMC-1.

Given the structural similarity of the two classes of molecules, these
dichotomies are striking and point to separate formation and carbon chain
growth mechanisms in TMC-1. Whether these differences exist between the
analogous cyanopolyynes (HC2n+1N) and the pure hydrocarbon polyynes (H2C2n) is
unclear, as the latter have no dipole moment. To understand this, future
infrared observations of vibrational bands from these molecules toward TMC-1
would be instrumental in constraining their abundances and understanding their
formation chemistry.

Finally, utilizing a gas-grain chemical network, we modeled the formation and
destruction of CH3-terminated carbon chains in TMC-1. The model is able to
reproduce our observed column densities to within an order of magnitude, and
we see very similar trends in the exponentially decreasing abundances with
carbon chain length. However, we are not able to predict the observed drop-off
in column density for CH3C8H or C8H using this model.

The results of this analysis offer important insights into the formation of
carbon chain molecules in cold pre-stellar conditions, but it is important to
note that many other families of unsaturated carbon chain species are known to
form in TMC-1 (e.g. CnO, CnS, HCnO). Similar in-depth analyses of these groups
of molecules would be instrumental in improving our understanding the extent
and efficiency of carbon chemistry in TMC-1.

We thank the reviewers for their thoughtful notes on this manuscript. The
National Radio Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc. The Green Bank Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc. M. A. S. acknowledges support from the Virginia Space Grant Consortium
for this work.

## References

  * Agúndez et al. (2017) Agúndez, M., Cernicharo, J., Quintana-Lacaci, G., et al. 2017, A&A, 601, A4, doi: 10.1051/0004-6361/201630274
  * Askne et al. (1984) Askne, J., Hoglund, B., Hjalmarson, A., & Irvine, W. M. 1984, A&A, 130, 311
  * Berteloite et al. (2010) Berteloite, C., Le Picard, S. D., Balucani, N., Canosa, A., & Sims, I. R. 2010, Physical Chemistry Chemical Physics, 12, 3666
  * Broten et al. (1984) Broten, N. W., MacLeod, J. M., Avery, L. W., et al. 1984, ApJ, 276, L25, doi: 10.1086/184181
  * Brünken et al. (2007) Brünken, S., Gupta, H., Gottlieb, C. A., McCarthy, M. C., & Thaddeus, P. 2007, ApJ, 664, L43, doi: 10.1086/520703
  * Burkhardt et al. (2018) Burkhardt, A. M., Herbst, E., Kalenskii, S. V., et al. 2018, Monthly Notices of the Royal Astronomical Society, 474, 5068, doi: 10.1093/mnras/stx2972
  * Burkhardt et al. (2021a) Burkhardt, A. M., Loomis, R. A., Shingledecker, C. N., et al. 2021a, Nature Astronomy, 5, 181, doi: 10.1038/s41550-020-01253-4
  * Burkhardt et al. (2021b) Burkhardt, A. M., Long Kelvin Lee, K., Bryan Changala, P., et al. 2021b, The Astrophysical Journal, 913, L18, doi: 10.3847/2041-8213/abfd3a
  * Cabezas et al. (2021) Cabezas, C., Tercero, B., Agúndez, M., et al. 2021, A&A, 650, L9, doi: 10.1051/0004-6361/202141274
  * Cernicharo et al. (2021) Cernicharo, J., Agúndez, M., Cabezas, C., et al. 2021, A&A, 649, L15, doi: 10.1051/0004-6361/202141156
  * Cernicharo et al. (2020) Cernicharo, J., Marcelino, N., Agúndez, M., et al. 2020, A&A, 642, L17, doi: 10.1051/0004-6361/202039351
  * Chen et al. (1998) Chen, W., Grabow, J.-U., Travers, M., et al. 1998, Journal of Molecular Spectroscopy, 192, 1, doi: https://doi.org/10.1006/jmsp.1998.7665
  * Chen et al. (1998) Chen, W., Grabow, J. U., Travers, M. J., et al. 1998, Journal of Molecular Spectroscopy, 192, 1, doi: 10.1006/jmsp.1998.7665
  * Churchwell et al. (1978) Churchwell, E., Winnewisser, G., & Walmsley, C. M. 1978, A&A, 67, 139
  * Dobashi et al. (2018) Dobashi, K., Shimoikura, T., Nakamura, F., et al. 2018, The Astrophysical Journal, 864, 82, doi: 10.3847/1538-4357/aad62f
  * Dobashi et al. (2019) Dobashi, K., Shimoikura, T., Ochiai, T., et al. 2019, The Astrophysical Journal, 879, 88, doi: 10.3847/1538-4357/ab25f0
  * Fehér et al. (2016) Fehér, O., Tóth, L. V., Ward-Thompson, D., et al. 2016, A&A, 590, A75, doi: 10.1051/0004-6361/201424385
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Garrod et al. (2007) Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103, doi: 10.1051/0004-6361:20066704
  * Gelman & Rubin (1992) Gelman, A., & Rubin, D. B. 1992, in Bayesian Statistics (Oxford University Press), 625 – 631
  * Gordy & Cook (1984) Gordy, W., & Cook, R. 1984, Microwave Molecular Spectra (John Wiley & Sons)
  * Harada & Herbst (2008) Harada, N., & Herbst, E. 2008, The Astrophysical Journal, 685, 272, doi: 10.1086/590468
  * Hébrard et al. (2009) Hébrard, E., Dobrijevic, M., Pernot, P., et al. 2009, Journal of Physical Chemistry A, 113, 11227, doi: 10.1021/jp905524e
  * Herbst & Cuppen (2006) Herbst, E., & Cuppen, H. M. 2006, Proceedings of the National Academy of Science, 103, 12257, doi: 10.1073/pnas.0601556103
  * Herbst & Leung (1989) Herbst, E., & Leung, C. M. 1989, The Astrophysical Journal Supplement Series, 69, 271, doi: 10.1086/191314
  * Hincelin et al. (2011) Hincelin, U., Wakelam, V., Hersant, F., et al. 2011, Astronomy & Astrophysics, 530, A61, doi: 10.1051/0004-6361/201016328
  * Hoang & Tung (2019) Hoang, T., & Tung, N.-D. 2019, ApJ, 885, 125, doi: 10.3847/1538-4357/ab4810
  * Irvine et al. (1981) Irvine, W. M., Hoglund, B., Friberg, P., Askne, J., & Ellder, J. 1981, ApJ, 248, L113, doi: 10.1086/183637
  * Kumar et al. (2019) Kumar, R., Carroll, C., Hartikainen, A., & Martin, O. 2019, Journal of Open Source Software, 4, 1143, doi: 10.21105/joss.01143
  * Lee & McGuire (2020) Lee, K. L. K., & McGuire, B. A. 2020, molsim, Zenodo.  https://zenodo.org/record/4122749
  * Lee et al. (2021) Lee, K. L. K., Loomis, R. A., Burkhardt, A. M., et al. 2021, ApJ, 908, L11, doi: 10.3847/2041-8213/abdbb9
  * Loomis et al. (2021) Loomis, R. A., Burkhardt, A. M., Shingledecker, C. N., et al. 2021, Nature Astronomy, 5, 188, doi: 10.1038/s41550-020-01261-4
  * Loomis et al. (2018) Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, The Astronomical Journal, 155, 182, doi: 10.3847/1538-3881/aab604
  * Loren et al. (1984) Loren, R. B., Wootten, A., & Mundy, L. G. 1984, ApJ, 286, L23, doi: 10.1086/184377
  * MacLeod et al. (1984) MacLeod, J. M., Avery, L. W., & Broten, N. W. 1984, ApJ, 282, L89, doi: 10.1086/184312
  * Matthews & Sears (1983) Matthews, H. E., & Sears, T. J. 1983, The Astrophysical Journal, 272, 149, doi: 10.1086/161271
  * McCarthy et al. (2020) McCarthy, M. C., Lee, K. L. K., Loomis, R. A., et al. 2020, Nature Astronomy, 1, doi: 10.1038/s41550-020-01213-y
  * McGuire et al. (2018) McGuire, B. A., Burkhardt, A. M., Kalenskii, S., et al. 2018, Science, 359, 202, doi: 10.1126/science.aao4890
  * McGuire et al. (2020) McGuire, B. A., Burkhardt, A. M., Loomis, R. A., et al. 2020, The Astrophysical Journal Letters, 900, L10, doi: 10.3847/2041-8213/aba632
  * McGuire et al. (2021) McGuire, B. A., Loomis, R. A., Burkhardt, A. M., et al. 2021, Science, 371, 1265, doi: 10.1126/science.abb7535
  * Mendoza et al. (2018) Mendoza, E., Bronfman, L., Duronea, N. U., et al. 2018, ApJ, 853, 152, doi: 10.3847/1538-4357/aaa1ec
  * Minh et al. (1993) Minh, Y. C., Irvine, W. M., Ohishi, M., et al. 1993, A&A, 267, 229
  * Minissale et al. (2016) Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24, doi: 10.1051/0004-6361/201525981
  * Quan & Herbst (2007) Quan, D., & Herbst, E. 2007, A&A, 474, 521, doi: 10.1051/0004-6361:20078246
  * Remijan et al. (2006) Remijan, A. J., Hollis, J. M., Snyder, L. E., Jewell, P. R., & Lovas, F. J. 2006, ApJ, 643, L37, doi: 10.1086/504918
  * Ruaud et al. (2016) Ruaud, M., Wakelam, V., & Hersant, F. 2016, Monthly Notice of the Royal Astronomical Society, 459, 3756, doi: 10.1093/mnras/stw887
  * Shingledecker et al. (2021) Shingledecker, C. N., Lee, K. L. K., Wandishin, J. T., et al. 2021, A&A, 652, L12, doi: 10.1051/0004-6361/202140698
  * Snyder et al. (2006) Snyder, L. E., Hollis, J. M., Jewell, P. R., Lovas, F. J., & Remijan, A. 2006, ApJ, 647, 412, doi: 10.1086/505323
  * Solomon et al. (1971) Solomon, P. M., Jefferts, K. B., Penzias, A. A., & Wilson, R. W. 1971, ApJ, 168, L107, doi: 10.1086/180794
  * Travers et al. (1998a) Travers, M. J., Chen, W., Grabow, J. U., McCarthy, M. C., & Thaddeus, P. 1998a, Journal of Molecular Spectroscopy, 192, 12, doi: 10.1006/jmsp.1998.7666
  * Travers et al. (1998b) —. 1998b, Journal of Molecular Spectroscopy, 192, 12, doi: 10.1006/jmsp.1998.7666
  * Turner et al. (2000) Turner, B. E., Herbst, E., & Terzieva, R. 2000, ApJS, 126, 427, doi: 10.1086/313301
  * Walmsley et al. (1984) Walmsley, C. M., Jewell, P. R., Snyder, L. E., & Winnewisser, G. 1984, A&A, 134, L11
  * Western (2017) Western, C. M. 2017, Journal of Quantitative Spectroscopy and Radiative Transfer, 186, 221, doi: 10.1016/j.jqsrt.2016.04.010
  * Western & Billinghurst (2019) Western, C. M., & Billinghurst, B. E. 2019, Physical Chemistry Chemical Physics, 21, 13986, doi: 10.1039/C8CP06493H
  * Willacy et al. (1993) Willacy, K., Williams, D. A., & Minh, Y. C. 1993, MNRAS, 263, L40, doi: 10.1093/mnras/263.1.L40
  * Woon & Herbst (2009) Woon, D. E., & Herbst, E. 2009, The Astrophysical Journal Supplement Series, 185, 273, doi: 10.1088/0067-0049/185/2/273
  * Xue et al. (2020) Xue, C., Willis, E. R., Loomis, R. A., et al. 2020, The Astrophysical Journal Letters, 900, L9, doi: 10.3847/2041-8213/aba631

##  Appendix A MCMC posterior distributions

Figure A1 and A2 shows the MCMC posterior distributions for CH3C7N and CH3C8H
respectively. Based on the off-diagonal heatmaps, we see that the majority of
modeling parameters do not demonstrate significant covariance; those that do
typically pertain to radial velocities and column densities of each velocity
component, often those adjacent to one another in velocity and within the same
symmetry group (i.e. $E$ state column densities in neighboring velocity
components). We do not expect significant covariance between the $A/E$ states
for both species, as they are by in large spectroscopically separate.

Regarding the treatment of source sizes, we observe quite different posteriors
in comparison to the earlier work done by Loomis et al. (2021), which
highlighted significant covariance between source sizes and column densities.
We do not observe this for CH3C7N most likely due to an overly constrained
prior placed on the source size, and we likely underestimate its uncertainty
and covariance. In the case of CH3C8H, we fixed the mean value of HC7N Loomis
et al. (2021), as we were interested in estimating upper limits to the column
densities, and thus not shown in Figure A2.

Comparison of the posterior distributions for the column densities between
CH3C7N and CH3C8H provides a margin for detection and non-detection cases. In
the former, where the matched filtering shown in Figure 1 indicates a
$4.6\sigma$ likelihood for CH3C7N, corroborating with smaller uncertainties in
the column densities in comparison to CH3C8H, particularly for the first
velocity component. We interpret this as a smaller range of column density
values provide evidential support for the tentative detection of CH3C7N: even
with the absence of individual lines, the likelihood-based sampling is able to
place relatively large constraints on the possible values.

Figure A1: Corner plot for CH3C7N. The diagonal traces correspond to the
marginalized, cumulative posterior distributions for each parameter, and off-
diagonal heatmaps represent parameter covariances. Vertical lines in the
cumulative distribution plots represent the first, second, and third
quantiles. Figure A2: Corner plot for CH3C8H. The diagonal traces correspond
to the marginalized, cumulative posterior distributions for each parameter,
and off-diagonal heatmaps represent parameter covariances. Vertical lines in
the cumulative distribution plots represent the first, second, and third
quantiles.

##  Appendix B Chemical Model Time Dependence

Figure A3 shows the time-dependence of the simulated abundances within the
nautilus chemical models in relation to the observed values. As it can be
seen, the relative trend lines of these species can be quite time dependent,
the time of peak abundance is strongly dependent on carbon chain length.

Figure A3: Simulated gas-phase abundance and column densities of the CH3CnN
(left) and CH3CnH (right) families from nautilus chemical models in comparison
to the observed values with uncertainties as a horizontal bars. The upper
limit of the CH3C8H is shown as green dotted line. The time used in Figure 4
is shown as a vertical dashed gray line.


11institutetext:  Biological and Environmental Sciences and Engineering
Division, Computer, Electrical and Mathematical Sciences and Engineering
Division, King Abdullah University of Science and Technology, Thuwal, Saudi
Arabia 11email: {chao-han.yang, jiahong.huang, hiromasa.morikawa,
<EMAIL_ADDRESS>
22institutetext: Georgia Institute of Technology, GA, USA
33institutetext: University of Waterloo, Canada 33email:
<EMAIL_ADDRESS>
44institutetext: Department of Ophthalmology, Bern University Hospital, Bern,
Switzerland 55institutetext: Department of Ophthalmology, Tri-Service General
Hospital, Taiwan66institutetext: National Taiwan University, Taiwan
77institutetext: Unit of Computational Medicine, Center for Molecular
Medicine, Department of Medicine, Karolinska Institutet, Sweden

# Auto-Classification of Retinal Diseases in the Limit of Sparse Data Using a
Two-Streams Machine Learning Model

C.-H. Huck Yang* 1122 Fangyu Liu 33 Jia-Hong Huang 11 Meng Tian 44 Hiromasa
Morikawa 11 I-Hung Lin M.D 55 Yi Chieh Liu 66 Hao-Hsiang Yang 66 Jesper Tegnèr
1177

###### Abstract

Automatic clinical diagnosis of retinal diseases has emerged as a promising
approach to facilitate discovery in areas with limited access to specialists.
Based on the fact that fundus structure and vascular disorders are the main
characteristics of retinal diseases, we propose a novel visual-assisted
diagnosis hybrid model mixing the support vector machine (SVM) and deep neural
networks (DNNs). Furthermore, we present a new clinical retina labels
collection sorted by the professional ophthalmologist from the educational
project Retina Image Bank, called EyeNet, for ophthalmology incorporating 52
retina diseases classes. Using EyeNet, our model achieves 90.40% diagnosis
accuracy, and the model performance is comparable to the professional
ophthalmologists.

##  1 Introduction

Computational retinal disease methods [1, 2] has been investigated extensively
through different signal processing techniques. Retinal diseases are
accessible to machine learning techniques due to their visual nature in
contrast to other common human diseases requiring invasive techniques for
diagnosis or treatments. Typically, the diagnosis accuracy of retinal diseases
based on the clinical retinal images is highly dependent on the practical
experience of a physician or ophthalmologist. However, training highly-skilled
ophthalmologists usually take years and the number of them, especially in the
less-developed area, is still far from enough. Therefore, developing an
automatic retinal diseases detection system is important, and it will broadly
facilitate diagnostic accuracy of retinal diseases. Moreover, for remote rural
areas, where there are even no ophthalmologists locally to screen retinal
disease, the automatic retinal diseases detection system also helps non-
ophthalmologists find the patient of the retinal disease, and further, refer
them to the medical center for further treatment.

Figure 1: This figure represents our proposed two-streams model. A raw retinal
image as an input of DNNs, U-Net, and as the other input to a contrast
enhancement algorithm. Then we pass the output of U-Net to two separated PCA
processing. Finally, the output of these two PCA modules is sent as inputs to
the retina disease classifier, SVM, which give the outcome of predicted retina
disease.

The development of automatic diseases detection (ADD) [3] alleviates enormous
pressure from social healthcare systems. Retinal symptom analysis [4] is one
of the important ADD applications given that it offers a unique opportunity to
improve eye care on a global stage. The World Health Organization estimates
that age-related macular degeneration (AMD) and Diabetic Retinopathy, which
are two typical retinal diseases, are expected to affect over 500 million
people worldwide by 2020 [5].

Moreover, the increasing number of cases of diabetic retinopathy globally
requires extending efforts in developing visual tools to assist in the
analytic of the series of retinal disease. These decision support systems for
retinal ADD, as [6] for non-proliferative diabetic retinopathy have been
improved from recent machine learning success on the high dimensional images
processing by featuring details on the blood vessel. [7] demonstrated an
automated technique for the segmentation of the blood vessels by tracking the
center of the vessels on Kalman Filter. However, these pattern recognition
based classification still rely on hand-crafted features and only specify for
evaluating single retinal symptom. Despite extensive efforts using wavelet
signal processing, retinal ADD remains a viable target for improved machine
learning techniques applicable for point-of-care (POC) medical diagnosis and
treatment in the aging society [8].

To the best of our knowledge, the amount of clinical retinal images are less
compared to other cell imaging data, such as blood cell and a cancer cell.
However, a vanilla deep learning based diseases diagnosis system requires
large amounts of data. Therefore, we propose a novel visual-assisted diagnosis
algorithm which is based on an integration of the support vector machine and
deep neural networks. The primary goal of this work is to automatically
classify 52 specific retinal diseases for human beings with the reliable
clinical-assisted ability on the intelligent medicine approaches. To foster
the long-term visual analytics research, we also present a visual clinical
label collection, EyeNet, including several crucial symptoms as AMD, DR,
uveitis, BRVO, BRAO.

Contributions.

  * •

We design a novel two-streams-based algorithm on the support vector machine
and deep neural networks to facilitate medical diagnosis of retinal diseases.

  * •

We present a new clinical labels collection, EyeNet, for Ophthalmology with 52
retina diseases classes as a crucial aid to the ophthalmologist and medical
informatics community.

  * •

Finally, we visualize the learned features inside the DNNs model by heat maps.
The visualization helps in understanding the medical comprehensibility inside
our DNNs model.

##  2 Related Work

In this section, we review some works related to our proposed method. We
divide the related works into three parts including medical dataset
comparison, dimension reduction by feature extraction, and image segmentation
by neural networks.

2.1 Medical Dataset Comparison

Large-scale datasets help the performance of deep learning algorithms
comparable to human-level on the tasks of speech recognition [9], image
classification and recognition [10], and question answering [11, 12, 13, 14,
15]. In the medical community, large scale medical datasets also help
algorithms achieve expert-level performance on detection of skin cancer [16],
diabetic retinopathy [17], heart arrhythmias [18], pneumonia [19], brain
hemorrhage [20], lymph node metastases [21], and hip fractures [22].

Recently, the number of openly available medical datasets is growing. In Table
1, we try to provide a summary of the publicly available medical image
datasets related to ours. According to Table 1, we notice that the recently
released ChestXray14 [23] is the largest medical dataset containing 112,120
frontal-view chest radiographs with up to 14 thoracic pathology labels.
Moreover, the smallest medical dataset is DRIVE [24] containing 40 retina
images. Regarding the openly available musculoskeletal radiograph databases,
the Stanford Program for Artificial Intelligence in Medicine and Imaging has a
medical dataset containing pediatric hand radiographs annotated with skeletal
age (AIMI). The Digital Hand Atlas [25] includes the left-hand radiographs
which are from children of various ages labeled with radiologist readings of
bone age. Then, our proposed EyeNet contains 52 classes of diseases and 1747
images.

Name of Dataset | Study Type | Label | Number of Images
---|---|---|---
EyeNet |  Retina |  Labels mining of Retinal Diseases |  1747
DRIVE [24] |  Retina |  Retinal Vessel Segmentation |  40
MURA [26] |  Musculoskeletal (Upper Extremity) |  Abnormality |  40,561
Digital Hand Atlas [25] |  Musculoskeletal (Left Hand) |  Bone Age |  1,390
ChestX-ray14 [23] |  Chest |  Multiple Pathologies |  112,120
DDSM [27] |  Mammogram |  Breast Cancer |  10,239

Table 1: Overview of available different types of medical label collection and
image datasets.

2.2 Dimension Reduction by Feature Extraction

Feature extraction is a method to make the task of pattern classification or
recognition easier. In image processing and pattern recognition, feature
extraction is one of the special forms of dimensionality reduction [28] in
some sense. The purpose of feature extraction is to exploit the most relevant
information based on the original data and describe the information in a space
with lower dimensions. For example, typically the size of original medical
image data, such as functional magnetic resonance imaging (fMRI) scans, is
very large and it causes algorithms computationally inefficient. In this case,
we will transform the original data into a reduced representation set of
features. That is, we exploit a set of feature vectors to describe the
original data and the process is called image feature extraction. In [29], the
authors mention that the representation by extracted feature vectors should
have a dimensionality that corresponds to the intrinsic dimensionality of the
original data. Then, the intrinsic dimensionality of data is the minimum
number of parameters required to account for the properties of the original
data. Moreover, the authors of [30] claim that dimensionality reduction
mitigates the curse of dimensionality and the other undesired properties of
spaces with high dimensions. The dimensionality reduction by feature
extraction method has been used in many different application fields such as
document verification [31], character recognition [32], extracting information
from sentences [14, 33, 13], machine translation [34, 35] and so on.

2.3 Image Segmentation by Neural Networks

Typically, researchers exploit the convolutional neural networks to do image
classification tasks with a single class output label. However, in the
biomedical image processing tasks, the output should contain localization.
That is to say, a class label is assigned to each pixel. Furthermore,
thousands of images in training set are typically beyond reach in the
biomedical tasks. Therefore, the authors of [36] train a neural networks
model, with sliding-window, to predict the output class label of each pixel by
providing a sub-region, small patch, around that pixel as input.

In [36], we know that the neural network model can do localization and the
number of training data, in the sense of patches, is much larger than the
training images. Apparently, [36] has two drawbacks. First, there exists some
trade-off between the use of context and localization accuracy. Then, since
the model runs separately for each small patch and there is much redundancy
due to overlapping patches, it is not efficient in the sense of computational
speed. Recently, the authors of [37, 38] have proposed an approach which can
do the good localization and use of context at the same time.

In the U-Net paper [39], the authors build upon an even more elegant neural
network architecture, the so-called fully convolutional network [40]. The
authors modify the architecture such that it works with very few training
images and produces even more accurate image segmentation. The main idea of
[40] is to supplement a usual contracting neural network by the successive
layers. Then, the authors exploit upsampling operators to replace pooling
operators, so the resolution of output is enhanced by these layers. In order
to do localization, the authors combine the upsampled output and high-
resolution features from the contracting path. Furthermore, a successive
convolutional layer learns to assemble a more accurate output based on this
information. Due to the advantages of U-Net mentioned above, we modify and
incorporate the U-Net to our proposed method.

Figure 2: The figure shows the result of U-Net effects on (a), unseen eyeball
clinical images with different morphological shapes. (b) is the ground truth
and (c) is the generated result of vessel-subtracted U-Net. Based on (b) and
(c), we discover that the results are highly similar to the ground truth.

##  3 Methodology

In this section, we present the workflow of our proposed model, referring to
Figure 1.

3.1 U-Net

DNNs has greatly boosted the performance of image classification due to its
power of image feature learning [41]. The active retinal disease is
characterized by exudates around retinal vessels resulting in cuffing of the
affected vessels [42]. However, ophthalmology images from clinical microscopy
are often overlayed with white sheathing and minor features. Segmentation of
retinal images has been investigated as a critical [43] visual-aid technique
for ophthalmologists. U-Net [39] is a functional DNNs, especially for
segmentation. Here, we proposed a modified version of U-Net by reducing the
copy and crop processes with a factor of two. The adjustment could speed up
the training process and have been verified as an adequate semantic effect on
small size images. We use cross-entropy for evaluating the training processes
as:

$\ E=\sum_{x\in\Omega}w(x)log(p_{l}(x))\hskip 35.56593pt(1)$

where $p_{l}$ is the approximated maximum function, and the weight map is then
computed as:

$\ w(x)=w_{c}(x)+w_{0}\cdot
exp(\frac{-(d_{x1}+d_{x2})^{2}}{2\sigma^{2}})\hskip 14.22636pt(2)$

$d_{x1}$ designates the distance to the border of the nearest edges and
$d_{x2}$ designates the distance to the border of the second nearest edges. LB
score is shown as [8]. We use the deep convolutional neural network (CNN) of
two $3\times 3$ convolutions. Each step followed by a rectified linear unit
(ReLU) and a $2\times 2$ max pooling operation with stride 2 for downsampling;
a layer with an even x- and y-size is selected for each operation. For the
U-Net model, we use existing DRIVE [24] dataset as the training segmentation
mask. Then, we use Our proposed model converges at the 44th epoch when the
error rate of the model is lower than $0.001$. The Jaccard similarity of our
U-Net model is 95.59% by validated on a 20% test set among EyeNet shown in
Figure 2. This model is robust and feasible for different retinal symptoms as
illustrated in Figure 3. The area under ROC curve is $0.979033$ and the area
under the Precision-Recall curve is $0.910335$.

Figure 3: This figure illustrates the qualitative results of contrast
enhancement algorithm from the original clinical images to the (b) histogram
equalization (c) contrast-limited adaptive histogram equalization.

3.2 Principal Component Analysis as Eigenface in the limit of Sparse Data

$\lambda_{k}=\frac{1}{M}\sum_{n=1}^{M}(u^{T}_{k}\Phi_{n})^{2}~{}~{}~{}~{}~{}(3)$

Eigenface [44] is classical and high-efficient image recognition technique
derived from the covariance matrix of the probability distribution over the
high-dimensional vector space of face images. Even with a single training
image, previous works [45, 46] of eigenface already established robust
automatic classification with confident accuracy ($85.6\%$) by combined
principal component analysis (PCA) and SVM classifiers. As a biological
feature, retinal images share similar properties with the human face for a
potential with eigenface recognition [47] included finite semantic layout
between facial features and ophthalmological features [48]. The eigenface
could be calculated [44] by maximizing the equivalent $(3)$, where $\Phi_{n}$
represent the face differ, $u_{k}$ is a chosen $k_{th}$ vector, $\lambda_{k}$
is the $k_{th}$ eigenvalue, and $M$ is a number of the dimension space. In our
experiment, we select the $k_{UNet}=40$ and $k_{RGB}=61$ to generate a
eigenface with highest accuracy for the U-Net-stream and RGB-stream
separately.

3.3 Support Vector Machine

Support Vector Machine is a machine learning technique for classification,
regression, and other learning tasks. Support vector classification (SVC) in
SVM, map data from an input space to a high-dimensional feature space, in
which an optimal separating hyperplane that maximizes the boundary margin
between the two classes is established. The hinge loss function is shown as:

$\
\frac{1}{n}\left[\sum_{i=1}^{n}max(0,1-y_{i}(\vec{w}\cdot\vec{x_{i}}-b))\right]+\lambda\left\|\vec{w}\right\|^{2}\hskip
14.22636pt(4)$

Where the parameter $\lambda$ determines the trade off between increasing the
margin-size and ensuring that the $\vec{x_{i}}$ lies on the right side of the
margin. We use radial basis function (RBF) and polynomial kernel for SVC,
which have been widely discussed [49] as a kernel-based fast SVC for images.

3.4 Contrast Enhancement

Contrast enhancement techniques play a vital role in image processing to bring
out the information that exists within a less dynamic range of that image. As
a major clinical feature, fundus [50, 51] structure is highly related [52] the
image contrast [53]. Here, we use histogram equalization for the contrast
enhancement in retinal images. Compare to the original images, images after
histogram equalization show the light color detail as lesions in Figure 3(b).
Images after contrast-limited adaptive histogram equalization (CLAHE) give
further features as areas of retinopathy in Figure 3(c).

##  4 Efforts on Retinal Dataset

Retina Image Bank (RIB) [54] is an international clinical project launched by
American Society of Retina Specialists in 2012, which allow retina specialists
and ophthalmic photographers around the world to share the existing clinical
cases online for medicine-educational proposes for patients and physicians in
developing countries lack training resource. Any researcher could join as a
contributor for dedicating the retinal images or as a visitor using the
medical images and label for non-commercial propose. With the recent success
on dataset collection, such as ImageNet [55], we believe that the effort of
sorting and mining the clinical labels from RIB is valuable. With a more
developer-friendly information pipeline, both Ophthalmology and Computer
Vision community could go further on the analytical researches on medical
informatics. Our proposed label collection, EyeNet is mainly based on the RIB
and following the RIB’s using guideline.

Figure 4: We use figure (i, j), where $i=a,b,c$ and $j=1,2,3,4$, to
demonstrate that our proposed method can capture the similar lesion areas as
the ophthalmologist’s manual annotations, i.e., the yellow sketches.

##  5 Experiments

In this section, we describe the implementation details and experiments we
conducted to validate our proposed method.

5.1 Label Collection

For experiments, the EyeNet is randomly divided into three parts: 70% for
training, 10% for validation and 20% for testing. All the training data have
to go through the PCA before SVM. All classification experiments are trained
and tested on the EyeNet.

5.2 Setup

The EyeNet has been processed to U-Net to generate a subset with a semantic
feature of the blood vessel. For the DNNs and Transfer Learning models, we
directly use the RGB images from the retinal dataset. EyeNet will be published
online after getting accepted. For the CLAHE processing, we use $adapthisteq$
function from the image toolbox in MATLAB.

5.3 Deep Convolutional Neural Networks

CNN has demonstrated extraordinary performance in visual recognition tasks
[55], and the state of the art is in a great many vision-related benchmarks
and challenges [56]. With little or no prior knowledge and human effort in
feature design, it yet provides a general and effective method solving variant
vision tasks in variant domains. This new development in computer vision has
also shown great potential for helping/replacing human judgment in vision
problems like medical imaging [16], which is the topic we try to address in
this paper. In this section, we introduce several baselines in multi-class
image recognition and compare their results on the EyeNet.

Baseline1-AlexNet

AlexNet [55] brought up a succinct network architecture, with 3 fully
connected layers, 5 convolutional layers, and the activation function is ReLU
[57].

Baseline2-VGG

The authors of VGG [41] exploit the filters (3x3) repeatedly to replace the
large filters (5x5,7x7) in traditional architectures. By increasing depths of
the network, it achieved better results on ImageNet with fewer parameters.

Baseline3-ResNet

Residual Networks [58], one of the most popular neural networks today, utilize
skip connections or short-cuts to jump over some layers. With skip
connections, the network essentially collapses into a shallower network in the
initial phase and this makes it easier to be trained, and then it expands its
layers as it learns more of the feature space.

Baseline4-SqueezeNet

In real world, medical imaging tasks usually require a small and effective
model to adapt to limited resources. As some deep neural networks cost several
hundred megabytes to store, SqueezeNet [59] adopting model compression
techniques has achieved the accuracy of AlexNet level with around 500 times
smaller models.

5.4 Transfer Learning

We use a transfer learning framework from the normalized ImageNet [55] to the
EyeNet for solving the small samples issue on the computational retinal visual
analytics. With sufficient and utilizable training classified model, Transfer
Learning resolves the challenge of Machine Learning in the limit of minimal
amount of training labels and it drastically reduces the data requirements.
The first few layers of DNNs learn features , similar to Gabor filters and
color blobs, and these features appear not to be specific to any particular
task or dataset and thus applicable to other datasets and tasks [60]. Our
experiments show the significant improvement after we apply the pretrained
parameters on our deep learning based models, referring to Table 3 and Table
4.

Hybrid-Ratio | RBF | Polyn.
---|---|---
0% : 100% |  0.8159 |  0.8391
40% : 60% |  0.8371 |  0.8381
50% : 50% |  0.9002 |  0.8632
60% : 40% |  0.8881 |  0.9040
100% : 0% |  0.8324 |  0.8241

Table 2: Accuracy comparison of the two-streams model with Radial basis
function (RBF) and polynomial kernel. We use the hybrid-ratio [61] of the
mixed weighted voting between two multi-SVCs trained from images over U-Net
and CLAHE.

Model | Pretrained | Random Init.
---|---|---
AlexNet |  0.7912 |  0.4837
VGG11 |  0.8802 |  0.7579
VGG13 |  0.8721 |  0.7123
ResNet18 |  0.8805 |  0.7250
SqueezeNet |  0.8239 |  0.5625

Table 3: Accuracy comparison of three DNNs baselines on EyeNet.

Model | Pretrained | Random Init.
---|---|---
AlexNet |  0.7952 |  0.4892
VGG11 |  0.8726 |  0.7583
VGG13 |  0.8885 |  0.7588
ResNet18 |  0.8834 |  0.6741
SqueezeNet |  0.8349 |  0.5721

Table 4: Accuracy comparison of three DNNs baselines on EyeNet.

5.5 Two-Streams Results

All SVM has implemented in Matlab with libsvm [62] module. We separate both
the original retinal dataset and the subset to three parts included 70%
training set, 20% test set, and 10% validation set. By training two multiple-
classes SVM models on both original EyeNet and the subset, we implement a
weighted voting method to identify the candidate of retina symptom. We have
testified different weight ratio as $Hybrid-Ratio$, SVM model with {Images
over CLAHE: Image over U-Net }, with different accuracy at Table 2. We have
verified the model without over-fitting by the validation set via
normalization on the accuracy with 2̃.03% difference.

5.6 Deep Neural Networks Results

All DNNs are implemented in PyTorch. We use identical hyperparameters for all
models. The training lasts 400 epochs. The first 200 epochs take a learning
rate of 1e-4 and the second 200 take 1e-5. Besides, we apply random data
augmentation during training. In every epoch, there is $70\%$ probability for
a training sample to be affinely transformed by one of the operations in
{flip, rotate, transpose}$\times${random crop}. Though ImageNet and our
Retinal Dataset are much different, using weights pretrained on ImageNet
rather than random ones has boosted test accuracy of any models with 5 to 15
percentages, referring to Table 3. Besides, pretrained models tend to converge
much faster than random initialized ones as suggested in Figure 4. The
performance of DNNs on our retinal dataset can greatly benefit from the
knowledge of other domains.

5.7 Neuron Visualization for Medical Images

Importantly, we verified the hypothesis that vessel-based segmentation and
contrast enhancement are two coherent features to decide the type of retinal
diseases. Using techniques of generating class activation maps introduced by
[63], we visualized feature maps of the final convolution layer of ResNet18
(which is one of our deep learning model baselines). We notice that the
features learned by deep learning models agree with our intuitions about
developing the two-stream machine learning model. In fact, in the clinical
diagnosis process, ”vessel patterns” and ”fundus structure” are also the two
most crucial features to identify the symptom of different diseases. These two
types of features actually cover more than 80% of retinal diseases [50, 51].

##  6 Conclusion and Future Work

In this work, we have designed a novel hybrid model for visual-assisted
diagnosis based on the SVM and U-Net. The performance of this model shows the
higher accuracy, 90.40%, over the other pre-trained DNNs models as an aid for
ophthalmologists. Also, we propose the EyeNet to benefit the medical
informatics research community. Finally, since our dataset not only contains
images but also text information of the images, image captioning and Visual
Question Answering [13, 14, 15] based on the retinal images are also the
interesting future directions. Our work may also help the remote rural area,
where there are no ophthalmologists locally, to screen retinal disease without
the help of ophthalmologists in the future.

## References

  * [1] Tan, O., Chopra, V., Lu, A.T.H., Schuman, J.S., Ishikawa, H., Wollstein, G., Varma, R., Huang, D.:  Detection of macular ganglion cell loss in glaucoma by fourier-domain optical coherence tomography.  Ophthalmology 116 (2009) 2305–2314
  * [2] Lalezary, M., Medeiros, F.A., Weinreb, R.N., Bowd, C., Sample, P.A., Tavares, I.M., Tafreshi, A., Zangwill, L.M.:  Baseline optical coherence tomography predicts the development of glaucomatous change in glaucoma suspects.  American journal of ophthalmology 142 (2006) 576–582
  * [3] Sharifi, M., Fathy, M., Mahmoudi, M.T.:  A classified and comparative study of edge detection algorithms.  In: Information Technology: Coding and Computing, 2002. Proceedings. International Conference on, IEEE (2002) 117–120
  * [4] Abràmoff, M.D., Garvin, M.K., Sonka, M.:  Retinal imaging and image analysis.  Reviews in biomedical engineering 3 (2010) 169–208
  * [5] Pizzarello, L., Abiose, A., Ffytche, T., Duerksen, R., Thulasiraj, R., Taylor, H., Faal, H., Rao, G., Kocur, I., Resnikoff, S.:  Vision 2020: The right to sight: a global initiative to eliminate avoidable blindness.  Archives of ophthalmology 122 (2004) 615–620
  * [6] Bhattacharya, S.:  Watermarking digital images using fuzzy matrix compositions and ($\alpha$, $\beta$)-cut of fuzzy set.  International Journal of Advanced Computing (2014)
  * [7] Lin, C.Y., Wu, M., Bloom, J.A., Cox, I.J., Miller, M.L., Lui, Y.M.:  Rotation-, scale-, and translation-resilient public watermarking for images.  In: Security and watermarking of multimedia contents II. Volume 3971., International Society for Optics and Photonics (2000) 90–99
  * [8] Cochocki, A., Unbehauen, R.:  Neural networks for optimization and signal processing.  John Wiley & Sons, Inc. (1993)
  * [9] Hannun, A., Case, C., Casper, J., Catanzaro, B., Diamos, G., Elsen, E., Prenger, R., Satheesh, S., Sengupta, S., Coates, A., et al.:  Deep speech: Scaling up end-to-end speech recognition.  arXiv preprint arXiv:1412.5567 (2014)
  * [10] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.:  Imagenet: A large-scale hierarchical image database.  In: Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, IEEE (2009) 248–255
  * [11] Rajpurkar, P., Zhang, J., Lopyrev, K., Liang, P.:  Squad: 100,000+ questions for machine comprehension of text.  arXiv preprint arXiv:1606.05250 (2016)
  * [12] Antol, S., Agrawal, A., Lu, J., Mitchell, M., Batra, D., Lawrence Zitnick, C., Parikh, D.:  Vqa: Visual question answering.  In: Proceedings of the ICCV. (2015) 2425–2433
  * [13] Huang, J.H., Alfadly, M., Ghanem, B.:  Vqabq: visual question answering by basic questions.  arXiv:1703.06492 (2017)
  * [14] Huang, J.H., Dao, C.D., Alfadly, M., Ghanem, B.:  A novel framework for robustness analysis of visual qa models.  arXiv:1711.06232 (2017)
  * [15] Huang, J.H., Alfadly, M., Ghanem, B.:  Robustness analysis of visual qa models by basic questions.  arXiv:1709.04625 (2017)
  * [16] Esteva, A., Kuprel, B., Novoa, R.A., Ko, J., Swetter, S.M., Blau, H.M., Thrun, S.:  Dermatologist-level classification of skin cancer with deep neural networks.  Nature 542 (2017) 115
  * [17] Gulshan, V., Peng, L., Coram, M., Stumpe, M.C., Wu, D., Narayanaswamy, A., Venugopalan, S., Widner, K., Madams, T., Cuadros, J., et al.:  Development and validation of a deep learning algorithm for detection of diabetic retinopathy in retinal fundus photographs.  Jama 316 (2016) 2402–2410
  * [18] Rajpurkar, P., Hannun, A.Y., Haghpanahi, M., Bourn, C., Ng, A.Y.:  Cardiologist-level arrhythmia detection with convolutional neural networks.  arXiv preprint arXiv:1707.01836 (2017)
  * [19] Rajpurkar, P., Irvin, J., Zhu, K., Yang, B., Mehta, H., Duan, T., Ding, D., Bagul, A., Langlotz, C., Shpanskaya, K., et al.:  Chexnet: Radiologist-level pneumonia detection on chest x-rays with deep learning.  arXiv preprint arXiv:1711.05225 (2017)
  * [20] Grewal, M., Srivastava, M.M., Kumar, P., Varadarajan, S.:  Radnet: Radiologist level accuracy using deep learning for hemorrhage detection in ct scans.  In: Biomedical Imaging (ISBI 2018), 2018 IEEE 15th International Symposium on, IEEE (2018) 281–284
  * [21] Bejnordi, B.E., Veta, M., van Diest, P.J., van Ginneken, B., Karssemeijer, N., Litjens, G., van der Laak, J.A., Hermsen, M., Manson, Q.F., Balkenhol, M., et al.:  Diagnostic assessment of deep learning algorithms for detection of lymph node metastases in women with breast cancer.  Jama 318 (2017) 2199–2210
  * [22] Gale, W., Oakden-Rayner, L., Carneiro, G., Bradley, A.P., Palmer, L.J.:  Detecting hip fractures with radiologist-level performance using deep neural networks.  arXiv preprint arXiv:1711.06504 (2017)
  * [23] Wang, X., Peng, Y., Lu, L., Lu, Z., Bagheri, M., Summers, R.M.:  Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In: 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), IEEE (2017) 3462–3471
  * [24] Staal, J., Abràmoff, M.D., Niemeijer, M., Viergever, M.A., Van Ginneken, B.:  Ridge-based vessel segmentation in color images of the retina.  TMI 23 (2004) 501–509
  * [25] Gertych, A., Zhang, A., Sayre, J., Pospiech-Kurkowska, S., Huang, H.:  Bone age assessment of children using a digital hand atlas.  Computerized Medical Imaging and Graphics 31 (2007) 322–331
  * [26] Rajpurkar, P., Irvin, J., Bagul, A., Ding, D., Duan, T., Mehta, H., Yang, B., Zhu, K., Laird, D., Ball, R.L., et al.:  Mura dataset: Towards radiologist-level abnormality detection in musculoskeletal radiographs.  arXiv preprint arXiv:1712.06957 (2017)
  * [27] Heath, M., Bowyer, K., Kopans, D., Moore, R., Kegelmeyer, P.:  The digital database for screening mammography.  Digital mammography (2000) 431–434
  * [28] Costa, J.A., Hero, A.:  Classification constrained dimensionality reduction.  In: Acoustics, Speech, and Signal Processing, 2005. Proceedings.(ICASSP’05). IEEE International Conference on. Volume 5., IEEE (2005) v–1077
  * [29] Fukunaga, K.:  Introduction to statistical pattern recognition.  Academic press (2013)
  * [30] Jimenez, L.O., Landgrebe, D.A.:  Supervised classification in high-dimensional space: geometrical, statistical, and asymptotical properties of multivariate data.  IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews) 28 (1998) 39–54
  * [31] Yang, Z., He, X., Gao, J., Deng, L., Smola, A.:  Stacked attention networks for image question answering.  In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. (2016) 21–29
  * [32] Trier, Ø.D., Jain, A.K., Taxt, T.:  Feature extraction methods for character recognition-a survey.  Pattern recognition 29 (1996) 641–662
  * [33] Srihari, R., Li, W.:  Information extraction supported question answering.  Technical report, CYMFONY NET INC WILLIAMSVILLE NY (1999)
  * [34] Somers, H.:  Example-based machine translation.  Machine Translation 14 (1999) 113–157
  * [35] Bahdanau, D., Cho, K., Bengio, Y.:  Neural machine translation by jointly learning to align and translate.  arXiv preprint arXiv:1409.0473 (2014)
  * [36] Ciresan, D., Giusti, A., Gambardella, L.M., Schmidhuber, J.:  Deep neural networks segment neuronal membranes in electron microscopy images.  In: Advances in neural information processing systems. (2012) 2843–2851
  * [37] Hariharan, B., Arbeláez, P., Girshick, R., Malik, J.:  Hypercolumns for object segmentation and fine-grained localization.  In: Proceedings of the IEEE conference on computer vision and pattern recognition. (2015) 447–456
  * [38] Seyedhosseini, M., Sajjadi, M., Tasdizen, T.:  Image segmentation with cascaded hierarchical models and logistic disjunctive normal networks.  In: Computer Vision (ICCV), 2013 IEEE International Conference on, IEEE (2013) 2168–2175
  * [39] Ronneberger, O., Fischer, P., Brox, T.:  U-net: Convolutional networks for biomedical image segmentation.  In: International Conference on MICCAI, Springer (2015) 234–241
  * [40] Long, J., Shelhamer, E., Darrell, T.:  Fully convolutional networks for semantic segmentation.  In: Proceedings of the IEEE conference on computer vision and pattern recognition. (2015) 3431–3440
  * [41] Simonyan, K., Zisserman, A.:  Very deep convolutional networks for large-scale image recognition.  arXiv:1409.1556 (2014)
  * [42] Khurana, A.:  Comprehensive ophthalmology.  New Age International Ltd (2007)
  * [43] Rezaee, K., Haddadnia, J., Tashk, A.:  Optimized clinical segmentation of retinal blood vessels by using combination of adaptive filtering, fuzzy entropy and skeletonization.  Applied Soft Computing 52 (2017) 937–951
  * [44] Turk, M.A., Pentland, A.P.:  Face recognition using eigenfaces.  In: Computer Vision and Pattern Recognition, 1991. Proceedings CVPR’91., IEEE Computer Society Conference on, IEEE (1991) 586–591
  * [45] Lyons, M.J., Budynek, J., Akamatsu, S.:  Automatic classification of single facial images.  IEEE Transactions on pattern analysis and machine intelligence 21 (1999) 1357–1362
  * [46] Wu, J., Zhou, Z.H.:  Face recognition with one training image per person.  Pattern Recognition Letters 23 (2002) 1711–1719
  * [47] Moghaddam, B., Wahid, W., Pentland, A.:  Beyond eigenfaces: Probabilistic matching for face recognition.  In: Automatic Face and Gesture Recognition, 1998. Proceedings. Third IEEE International Conference on, IEEE (1998) 30–35
  * [48] Akram, M.U., Tariq, A., Khan, S.A.:  Retinal recognition: Personal identification using blood vessels.  In: Internet Technology and Secured Transactions (ICITST), 2011 International Conference for, IEEE (2011) 180–184
  * [49] Kuo, B.C., Ho, H.H., Li, C.H., Hung, C.C., Taur, J.S.:  A kernel-based feature selection method for svm with rbf kernel for hyperspectral image classification.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 7 (2014) 317–326
  * [50] Crick, R.P., Khaw, P.T.:  A textbook of clinical ophthalmology: a practical guide to disorders of the eyes and their management.  (World Scientific)
  * [51] Akram, I., Rubinstein, A.:  Common retinal signs. an overview.  Optometry Today (2005)
  * [52] Tang, S., Huang, L., Wang, Y., Wang, Y.:  Contrast-enhanced ultrasonography diagnosis of fundal localized type of gallbladder adenomyomatosis.  BMC gastroenterology 15 (2015) 99
  * [53] Noyel, G., Thomas, R., Bhakta, G., Crowder, A., Owens, D., Boyle, P.:  Superimposition of eye fundus images for longitudinal analysis from large public health databases.  Biomedical Physics & Engineering Express 3 (2017) 045015
  * [54] :  Retina Image Bank: A project from the American Society of Retina Specialists.  (http://imagebank.asrs.org/about) Accessed: 2018-06-30.
  * [55] Krizhevsky, A., Sutskever, I., Hinton, G.E.:  Imagenet classification with deep convolutional neural networks.  In: Advances in NIPS. (2012) 1097–1105
  * [56] Xie, S., Girshick, R., Dollár, P., Tu, Z., He, K.:  Aggregated residual transformations for deep neural networks.  In: CVPR, IEEE (2017) 5987–5995
  * [57] Nair, V., Hinton, G.E.:  Rectified linear units improve restricted boltzmann machines.  In: Proceedings of the 27th international conference on machine learning (ICML-10). (2010) 807–814
  * [58] He, K., Zhang, X., Ren, S., Sun, J.:  (Deep residual learning for image recognition)
  * [59] Iandola, F.N., Han, S., Moskewicz, M.W., Ashraf, K., Dally, W.J., Keutzer, K.:  Squeezenet: Alexnet-level accuracy with 50x fewer parameters and¡ 0.5 mb model size.  arXiv:1602.07360 (2016)
  * [60] Yosinski, J., Clune, J., Bengio, Y., Lipson, H.:  How transferable are features in deep neural networks?  In: NIPS. (2014) 3320–3328
  * [61] Yang, C.H.H., Huang, J.H., Liu, F., Chiu, F.Y., Gao, M., Lyu, W., Tegner, J., et al.:  A novel hybrid machine learning model for auto-classification of retinal diseases.  arXiv preprint arXiv:1806.06423 (2018)
  * [62] Chang, C.C., Lin, C.J.:  Libsvm: a library for support vector machines.  TIST 2 (2011) 27
  * [63] Zhou, B., Khosla, A., A., L., Oliva, A., Torralba, A.:  Learning Deep Features for Discriminative Localization.  CVPR (2016)