\section{Introduction} \label{sec:intro}

The rapid increase of observed exoplanets revolutionizes our understanding of planet formation \citep[e.g.,][for a review]{2015ARA&A..53..409W},
which has been accelerated thanks to the {\it Kepler} mission \citep{2010Sci...327..977B}.
One famous, astonishing result is that small-sized planets whose radius is larger than $1 R_{\oplus}$ and smaller than $ 4R_{\oplus}$,
are dominant in the observed population \citep[e.g.,][]{2010Sci...330..653H}.
This is consistent with the pioneering discovery of radial velocity observations \citep{2011arXiv1109.2497M}.

The ubiquity of observed small-sized planets has challenged the canonical picture of planet formation \citep[e.g.,][]{2004ApJ...604..388I}.
This is because there is no such analogue in the solar system,
and hence both the observational characterization of these planets and
the theoretical development of possible formation mechanisms have been required.
For the former, the follow-up observations of exoplanet-host stars reveal that
the radius range of small-sized planets is divided into two regimes \citep{2017AJ....154..109F}:
$1 R_{\oplus} \la R_p \la 1.8  R_{\oplus}$ and $1.8 R_{\oplus} \la R_p \la 4  R_{\oplus}$.
Based on the bulk density \citep[e.g.,][]{2014ApJ...783L...6W,2015ApJ...801...41R},
planets in the first regime can literally be called as ``super-Earth",
and those in the second regime may be named as ``sub-Neptune".
The presence of these bi-modal populations has been popularized as the "radius gap (or valley)", and triggered a number of theoretical studies.
Recent studies suggest that both the formation and evolution processes can generate the gap \citep[e.g.,][]{2018MNRAS.476..759G,2013ApJ...775..105O}.

The multiplicity is another important feature of observed small-sized planets \citep[e.g.,][]{2011Natur.470...53L},
which also sheds light on their formation mechanisms \citep[e.g.,][]{2013ApJ...775...53H,2016ApJ...822...54D}.
Intriguingly, \citet{2018AJ....155...48W} confirm the presence of the radius gap, focusing only on exoplanets in multi-planet systems.
This finding may suggest that physical processes inevitable for shaping multi-planet systems would contribute to the generation of the radius gap.

Here, we show through a simple theoretical argument and an observational analysis that
planetary migration produces different configurations of protoplanets with different masses
and that this difference and the subsequent giant impact lead to two distinguishable populations of small-sized multi-planet systems.
It is prominent that the dividing radius of these two populations broadly corresponds to the radius gap
when the gas surface density of the natal protoplanetary disk has a low value ($\sim 10^2$ g cm$^{-2}$) in the vicinity of the central star.
Thus, this work points out the importance of planetary migration and the inner disk properties
to better understand the observed properties of super-Earths and sub-Neptunes.

\section{Theoretical and Observational Analyses} \label{sec:mod}

\subsection{Exoplanet data} \label{sec:data}

We first introduce the exoplanetary data that will be used in the following analyses;
in these analyses, the mass ($M_p$), radius ($R_p$), semimajor axis ($a_p$) of planets and the mass ($M_s$) of the host stars are needed.

We obtain the {\it Kepler} data from \citet{keplercumulativetable}\footnote{
We accessed to https://exoplanetarchive.ipac.caltech.edu/cgi-bin/TblView/nph-tblView?app=ExoTbls\&config=cumulative
on 2021-01-21 at 16:17, and obtained the data with the size of 9564 $\times$ 49 columns.}
We focus only on multi-planetary systems orbiting around single stars,
which are comprised of both confirmed and candidate planets.
We filter out the systems that do not have certain quantities and hence our analyses cannot be applied to.
These are the surface gravity, effective temperature, and radius of the host stars, and the radius and orbital period of planets.
Note that the stellar parameters are needed both to identify the host stars as the main sequence ones, and to compute their masses;
for the former, we follow the approach of \citet{2018AJ....156...24M},
where dwarf stars are identified, based on the relationship between the surface gravity and effective temperature \citep{2016ApJS..224....2H}.
For the latter, we use the mass-luminosity and the mass-radius relations \citep{2018MNRAS.479.5491E}.\footnote{
Note that the stellar masses estimated by our approach are comparable to those obtained, using the Gaia data \citep{2020AJ....159..280B};
we find that the difference is small ($< 20$ \%) enough for most cases that our results do not change very much.}

We select planets whose radius is in the range of $0.1 R_{\oplus} \leq R_p  \leq 30 R_{\oplus}$.
Following \citet{2018ApJS..235...38T}, we adopt the disposition score of $> 0.6$,
where the trade-off between reliability and completeness is achieved properly (see the top panel of their figure 13).
We have confirmed that our results do not change very much even if planets with the disposition score of $> 0.5$ are chosen.
The semimajor axis of planets is computed, using the Kepler's third law.

In summary, we have the total of 348 planetary systems with 870 exoplanets observed by {\it Kepler}.\footnote{
Note that we have removed one planetary system, Kepler-1659, due to the unphysically high mass of one planet.}

\subsection{Mass-radius relation}

The planet mass is one fundamental quantity for our analyses and is computed from the mass-radius relation.

We utilize the pretrained model developed by \citet{2017ApJ...834...17C},
where an unbiased forecasting model was built upon for a probabilistic mass-radius relation with the Bayesian framework.
They found that the mass-radius relation is divided into two regimes for small-sized planets:
$ M_p \la 2  M_{\oplus}$, and $2 M_{\oplus} \la M_p \la 100  M_{\oplus}$,
and the corresponding regimes are referred to as the Terran worlds, and the Neptunian worlds, respectively.

Figure \ref{fig1} shows the computed masses for given radii.
We subsequently fit these data and confirm that
the resulting power-law indices for the two types of planets are similar to those of \citet{2017ApJ...834...17C}.
It is interesting that the Neptunian worlds reside in the radius range of $1.2 R_{\oplus} \la R_p \la 10  R_{\oplus}$,
and therefore they cover both (massive) super-Earths and the entire sub-Neptunes.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{fig1.pdf}
\caption{The mass-radius relation for exoplanets in our data set.
The cyan dots with the error bars represent planets in the Terran world (i.e., $M \la 2M_{\oplus}$),
and the light green dots with the error bars represent planets in the Neptunian worlds (i.e., $2M_{\oplus} \la  M \la100 M_{\oplus}$).
The corresponding fitting profiles are denoted by the blue dashed, and the green solid lines, respectively.}
\label{fig1}
\end{center}
\end{figure}

\subsection{Theoretical Prediction: the importance of planetary migration} \label{sec:theory}

The formation of small-sized, multi-planet systems is currently the target of active research
\citep[e.g., ][]{2010ApJ...719..810I,2013ApJ...775...53H,2013MNRAS.431.3444C,2017MNRAS.470.1750I},
and it is still far from a complete understanding.

Our theoretical argument assumes that type I migration and the subsequent giant impact play an important role in forming small-sized, multi-planet systems;
type I migration takes place due to the disk-planet interaction and becomes effective for low-mass planets \citep[$M_p \ga 1M_{\oplus}$, e.g.,][]{2012ARA&A..50..211K}.
Giant impact is the collision between protoplanets and serves as the final mass assembly for terrestrial planets in the solar system \citep[e.g.,][hereafter C01]{2001Icar..152..205C}.

Under this assumption, the importance of planetary migration to trigger the subsequent giant impact between neighbouring protoplanets may be determined
by the difference between the mutual spacing ($\Delta_{\rm mig}$) achieved by type I migration and that ($\Delta_{\rm esc}$) by the escape velocity;
$\Delta_{\rm esc}$ is the result of the pure gravitational interaction between the neighbouring protoplanets
and is the largest separation between them before they undergo giant impact \citep[e.g.,][]{2014ApJ...795L..15S}:
\begin{equation}
\Delta_{\rm esc} \simeq 2 a_p \frac{v_{\rm esc}}{v_{\rm Kep}},
\end{equation}
where $v_{\rm Kep}$ is the Keplerian velocity, and $v_{\rm esc}= \sqrt{2GM_p/R_p}$ is the escape velocity for the neighbouring protoplanets with comparable masses and radii.

Planetary migration affects the onset of the subsequent giant impact considerably if the following condition is met:
\begin{equation}
\Delta_{\rm mig} < \Delta_{\rm esc}.
\end{equation}
Mathematically, $\Delta_{\rm mig}$ can be computed from the consideration that
the speed ($v_{\rm mig}$) of (differential) type I migration, which reduces the mutual spacing ($b=|a_{p,i}-a_{p,j}|$) between the neighbouring protoplanets ($i,j$),
is compensated by the gravitational repulsion between them \citep{2010ApJ...719..810I}.
For $v_{\rm mig}$, it is written as \citep[e.g., ][]{2011MNRAS.417.1236H}
\begin{eqnarray}
v_{\rm mig}   \simeq 2 K_{\rm mig} \frac{\Sigma_{\rm g} a_p^2 M_p}{M_s^2} \left( \frac{H}{a_p} \right)^{-2} v_{\rm Kep},
\end{eqnarray}
where $\Sigma_{\rm g}$ and $H$ are the gas surface density and the pressure scale height of the natal protoplanetary disk, respectively,
and $K_{\rm mig}$ is a coefficient which accounts for the differential speed and the detail of disk properties.
For simplicity, we here adopt that $H/a_{p}=0.02$, $K_{\rm mig}=1$, and $M_{s} = 1 M_{\odot}$.
For the gravitational repulsion, the expansion ($\delta b$) of the mutual spacing is given by linear theory as \citep{1990A&A...227..619H}
\begin{equation}
\delta b \simeq 30 \left( \frac{b}{r_{\rm H}} \right)^{-5} r_{\rm H},
\end{equation}
where $r_{\rm H} = a (2M_{p}/3M_s)^{1/3}$ is the mutual Hill radius for the neighbouring protoplanets with comparable masses.
Since the expansion is driven by encounters which occur at every synodic period ($T_{\rm syn}$),
the expansion rate becomes
\begin{equation}
 \frac{d b}{dt} \simeq \frac{\delta b}{T_{\rm syn}},
\end{equation}
where $T_{\rm syn} \simeq 4 \pi a_p^2 / 3b v_{\rm Kep}$.
In summary, $\Delta_{\rm mig} $ can be computed, by equating $v_{\rm mig}$ with $\delta b/T_{\rm syn}$.

Figure \ref{fig2} shows the resulting behaviors of $\Delta_{\rm esc}$ and $\Delta_{\rm mig}$ as a function of planet mass.
In this plot, the mass-radius relation obtained in Figure \ref{fig1} is used to compute $\Delta_{\rm esc}$.
As an example, we adopt the value of $\Sigma_{\rm g}= 5 \times 10^{2}$ g cm$^{-2}$ at $a_p=0.1$ au to calculate $\Delta_{\rm mig}$.
In addition, the mass regime where type I migration becomes effective is denoted by the solid line of $\Delta_{\rm mig}$,
following \citet{2012ApJ...760..117H}.

We find that $\Delta_{\rm mig} > \Delta_{\rm esc}$ for the planet mass of  $M_p \la 1 M_{\oplus}$,
 $\Delta_{\rm mig} < \Delta_{\rm esc}$ for $1 M_{\oplus} \la M_p \la 10 M_{\oplus}$,
and $\Delta_{\rm mig} \simeq \Delta_{\rm esc}$ for $M_p \ga 10 M_{\oplus}$ in this particular setup.
Note that the mutual spacing of $20 r_{\rm H}$ corresponds to the period ratio of $\sim 1.5$,
which is consistent with the peak value of the observed exoplanets \citep[e.g.,][]{2020MNRAS.495.4192C}.
Our value is, however, larger than what the typical numerical simulations predict \citep[i.e., $\sim$ a few $r_{\rm H}$, e.g.,][]{2017MNRAS.470.1750I}.
This difference comes from a lower value of $\Sigma_{\rm g}$ (see Section \ref{sec:disc} for more discussion).

Given that the mutual spacing is one key parameter for giant impact,
it is natural to consider that
protoplanets with $\Delta_{\rm mig} \la \Delta_{\rm esc}$ undergo giant impact more efficiently/rapidly than those in $\Delta_{\rm mig} > \Delta_{\rm esc}$.
This consideration leads to the expectation that (proto)planets in the red shaded region may surely experience giant impact and obtain higher masses.
Note that the size of the region depends on the power-law index of the mass-radius relation;
assuming that $R_p \propto M_p^{\beta}$, $\Delta_{\rm mig} \propto M^{-1/12}$ and $\Delta_{\rm esc} \propto M^{(1/3-\beta)/2}$ in the unit of $r_{\rm H}$.
The region expands if $\beta > 1/3$ and shrinks if  $\beta < 1/3$ in the Neptunian worlds, compared with Figure \ref{fig2}.
However, the region surely exists if $\Sigma_{\rm g}$ is relatively small (see Section \ref{sec:disc}).

In the following, we conduct the observational analysis to verify our theoretical prediction.

\begin{figure
\begin{center}
\includegraphics[width=8cm]{fig2.pdf}
\caption{The mutual spacing as a function of planet mass for the case that $\Sigma_{\rm g}= 5 \times 10^{2}$ g cm$^{-2}$ at $a_p=0.1$ au.
Planetary migration achieves narrower mutual spacings for the mass range of $1 M_{\oplus} \la M_p \la 10 M_{\oplus}$ (see the red shaded region).
The mass-radius relation is used to plot $\Delta_{\rm esc}$ (see the black, dashed line).
For $\Delta_{\rm mig}$, the solid line denotes  the mass regime where type I migration becomes effective (see the magenta line).
The background cyan and light green regions represent the Terran and Neptunian worlds, respectively.}
\label{fig2}
\end{center}
\end{figure}

\subsection{Statistical measures}

The direct observabales (e.g., the period ratio of neighboring planets) can be used to characterize multi-planet systems \citep[e.g.,][]{2018AJ....155...48W,2018ApJ...860..101Z}.
However, as the number of planets in the systems increases, the complexity in systematically characterizing multi-planet systems as a whole increase \citep[e.g.,][]{2020AJ....159..281G}.

In order to resolve this issue, we adopt the so-called "statistical measures" that are devised by C01.
We here focus on two quantities that are computed directly from $M_p$, $a_p$, and $M_s$.

The first quantity is the mass concentration ($S_c$) that measures the degree to which the planet mass is concentrated in a certain location of the system
and is calculated as
\begin{equation}
    S_c = \mbox{max} \left( \frac{\sum_j  M_{p,j}}{\sum_j M_{p,j} [ \log_{10}{(a/a_{p,j})]^2}} \right),
\label{mc}
\end{equation}
where the summation is done for the index $j$.
By changing $a$, the maximum value is searched.
Based on $N$-body simulations,
the value of $S_c$ reflects the initial spatial distribution of protoplanets that undergo giant impact eventually
\citep[hereafter, H09 and HM12, respectively]{2009ApJ...703.1131H,2012ApJ...751..158H}.

The second quantity is the orbital spacing ($S_s$) which is somewhat similar to the averaged mutual spacing normalized by the mutual Hill radius.
The main difference is that $S_s$ is normalized by $M_p^{1/4}$ (not $M_p^{1/3}$)
and is motivated by the results of $N$-body simulations that explore the stability of multi-planet systems \citep{1996Icar..119..261C}:
\begin{equation}
    S_s = \frac{6}{N-1} \Bigg( \frac{a_{p,max} - a_{p,min}}{a_{p,max}+a_{p,min}}\Bigg)\Bigg(\frac{3 M_s}{2\ \overline{M}_p}\Bigg)^{1/4},
\label{oss}
\end{equation}
where $N$,  $a_{p,max}$, $a_{p,min}$, and $\overline{M}_p$ are the number, the maximum and the minimum values of the semi-major axis,
and the mean mass of the planets in a system, respectively.

\subsection{Observational Analysis: KS tests}

\begin{figure*}
\begin{minipage}{17cm}
\begin{center}
\includegraphics[width=8.5cm]{fig3a.pdf}
\includegraphics[width=8cm]{fig3b.pdf}
\caption{The results of our KS tests and the best case of the cumulative distributions of $S_c$ and $S_s$ (i.e., $M_{p,\rm th}=3.5M_\oplus$).
The left panel shows that $S_c$ can be used to divide the full sample of small-sized, multi-planet systems into two populations
in the mass range of $2.5 M_{\oplus} \la M_{p, \rm th} \la 6.3 M_{\oplus}$ (see the red shaded region).
In this region, the KS statistic (the red dashed line) is greater than the critical value (the black solid line),
and hence the null hypothesis is rejected with the 99.5\% confidence;
equivalently, low-mass ($M_{p,\rm ave} \leq M_{p, \rm th}$) and high-mass ($M_{p,\rm ave} > M_{p,\rm th}$) systems are not drawn from the same distribution.
On the contrary, $S_s$ does not show the similar feature.
The central and right panels show the best case,
where low-mass and high-mass systems are denoted by the blue dashed and the green solid lines, respectively.}
\label{fig3}
\end{center}
\end{minipage}
\end{figure*}

We now perform the Kolmogorov–Smirnov test (KS test) for the cumulative distributions of $S_c$ and $S_s$
to determine whether the population of small-sized, multi-planet systems can be divided into sub-groups.

The KS test is a non-parametric test and used to discover whether two samples are drawn from the same distribution, based on the null hypothesis.
The significance level (also referred to as $\alpha$) is the probability of rejecting the null hypothesis when it is actually true.
The KS test for the two samples returns a KS statistic ($\mathcal{K}$) and a p-value ($\mathcal{P}$).
Based on the value of $\alpha$, we calculate a critical value ($D_{\rm crit} $),
which is given as
\begin{equation}
    D_{\rm crit} = c(\alpha) \cdot \sqrt{\Bigg(\frac{1}{n_1} + \frac{1}{n_2}\Bigg)},
\label{dcrit}
\end{equation}
where $c(\alpha) = \sqrt{- \ln(\alpha/2)/2}$, and $n_1$ and $n_2$ are the two sample sizes we employ.
If $\mathcal{K}>D_{\rm crit}$ or $\mathcal{P}<\alpha$, then we can reject the null hypothesis at a level of $\alpha$.

We choose the average mass of planets in the systems ($M_{p, \rm ave}$) as a parameter to divide the full samples into two groups;
the current {\it Kepler} data suggest that at least the planet radius difference in the systems may be insignificant \citep{2018AJ....155...48W}.
This does not necessarily mean that the planet mass difference is small as well (see Figure \ref{fig1}).
However, the mass-radius relation is one best tool currently available in the literature to estimate the planet mass,
and hence we use the planet mass for our KS tests.
Practically, we bin the mass range of $1 M_{\oplus} \leq M_{p, \rm th} \leq 10 M_{\oplus}$ with the bin size of $0.5 M_{\oplus}$,
where $M_{p,\rm th}$ is a parameterized threshold mass to divide the full samples into two groups,
and compare these two samples: planetary systems with $M_{p,\rm ave} \leq M_{p, \rm th}$ and those with $M_{p,\rm ave} > M_{p,\rm th}$.

Figure \ref{fig3} summarizes our results;
the left panel shows the results of our KS tests, and the central and right ones depict the cumulative distributions of $S_c$ and $S_s$ for the best case (i.e., $M_{p, \rm th}=3.5 M_{\oplus}$).
We find that for the mass concentration,
the null hypothesis can be rejected for $2.48 M_{\oplus} < M_{p, \rm th} < 6.3 M_{\oplus}$ with the 99.5\% confidence (i.e., $\alpha = 0.005$);
this is evident from the KS statistic which is greater than the value of $D_{\rm crit}$ in this range.
Our analysis therefore suggests that the full samples can be divided into two groups at $M_{p, \rm th}=3.5 M_{\oplus}$.
On the other hand, similar distinction is not possible for the orbital spacing.

We discuss below how our observational analysis is consistent with our theoretical prediction,
and why the mass concentration and orbital spacing behave differently.

\section{Discussion} \label{sec:disc}

Our theoretical argument has suggested
that planetary migration should have played a more important role for planets with $M_{p} \ga 2 M_{\oplus}$ (Figure \ref{fig2}),
and our observational analysis has shown
that planetary systems with $M_{p, \rm ave} \la 3 M_{\oplus}$ have a higher value of $S_c$ than those with $M_{p, \rm ave} \ga 3 M_{\oplus}$ (Figure \ref{fig3}).
Based on these results, we here discuss how small-sized, multi-planet systems form.

\begin{figure
\begin{center}
\includegraphics[width=8cm]{fig4.pdf}
\caption{The cumulative distributions of $S_c$ and $S_s$ for the simulation results made by H09 and HM12.
The cumulative distribution of $S_c$ shows a clearer difference,
which comes from different distributions of protoplanets that eventually undergo giant impact.}
\label{fig4}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{fig5.pdf}
\caption{The optimal value of $\Sigma_{\rm g}$ as a function of the distance from the central star.
The maximum and minimum values of $\Sigma_{\rm g}$ are computed (see the black solid and dashed lines, respectively),
using the results of Figure \ref{fig3} (see the vertical dashed lines in the left panel).
For comparison, the MMSN model is denoted by the blue dotted line.}
\label{fig5}
\end{center}
\end{figure}

First, we point out that both of our analyses lead to a consistent picture that
super-Earth (i.e., $M_{p, \rm ave} \la 3 M_{\oplus}$) systems form out of narrow ring-like distributions of protoplanets
while sub-Neptune (i.e., $M_{p, \rm ave} \ga 3 M_{\oplus}$) systems form from radially wide distributions of protoplanets.
The difference in the spatial distribution of protoplanets is produced by planetary migration;
a chain of protoplanets tends to form readily with a wider radial extent for massive protoplanets due to faster migration.
Importantly, the existing $N$-body simulations already show that
different (narrow vs wide) spatial distributions of protoplanets leads to different (high vs low) values of $S_c$.
Figure \ref{fig4} shows the simulation results obtained by H09 and HM12;
the former targets the formation of terrestrial planets in the solar system and hence starts from narrow rings whose radial extent is comparable to the resulting value of $S_s$.
On the other hand, the latter does the formation of small-sized planets observed by {\it Kepler} and
begins with wide distributions whose radial extent is a few times larger than the resulting value of $S_s$.

It is obvious that quantitative comparison between the observations (Figure \ref{fig3}) and the simulations (Figure \ref{fig4}) cannot be made
due to the observational bias and idealization adopted in simulations.
For instance, the observational bias may partially wash out the effect of planetary migration;
this may be a reason of why there is no clear difference in the cumulative distribution of $S_s$ between low-mass and high-mass systems (Figure \ref{fig3}).
However, qualitative comparison may be useful;
it may be reasonable to consider that the difference in the cumulative distribution of $S_c$ is significant enough to infer the spatial distribution of protoplanets.

Second, we consider the value of $\Sigma_{\rm g}$ in the vicinity of the central star.
As discussed in Section \ref{sec:theory},
a lower value of $\Sigma_{\rm g}$ may be preferred,
in order to limit the effect of planetary migration (Figure \ref{fig2}).
Under the assumption that the population of low-mass and high-mass systems is divided at $M_{p, \rm th}$,
one can compute the optimal value of $\Sigma_{\rm g}$.
Figure \ref{fig5} shows the results.
It is interesting that the resulting values are much lower than the minimum-mass solar nebula (MMSN) model \citep{1981PThPS..70...35H}.
This difference is readily understood by the fact that there is no close-in planet in the solar system,
and hence the MMSN model is not sensitive to such a region.
Furthermore, the recent studies support the lower $\Sigma_{\rm g}$ in the vicinity of the central star \citep[e.g.,][]{2020MNRAS.495.4192C},
which may be caused by stellar magnetic fields \citep{2019A&A...629L...1H} and/or disk winds \citep{2015A&A...584L...1O}.

Finally, we comment on caveats involved in this work.
The most important one would be the observational bias as discussed above.
Since it would be complicated to reliably determine the observational bias for multi-planet systems,
it may be straightforward to run simulations and impose the bias on the simulation results as done by \citet{2020ApJ...897...72M,2020AJ....160..276H}.
This approach allows one to compare simulations and observations directly and quantitatively.
Detailed simulations are desired to verify the results of this work,
which is the target of our future work.

Thus, planetary migration and the properties of protoplanetary disks in the vicinity of the central stars are the key
to better understanding the origin and properties of observed super-Earths and sub-Neptunes.

\begin{acknowledgments}

The authors thank an anonymous referee for useful comments on our manuscript.
This research was carried out in part at JPL/Caltech, under a contract with NASA.
Y.H. is supported by JPL/Caltech.

\end{acknowledgments}

\bibliographystyle{aasjournal}


\section{Introduction}
Viral pandemics, antibiotic resistant bacteria, or fungal infections such as \textit{candida auris} are fundamental threats to human health \cite{aslam2018antibiotic,jeffery2018candida}. SARS-CoV-2 (COVID-19) shocked the world with its first appearance estimated in the Fall/Winter of 2019 to becoming a global crisis by March, 2020 when it was declared a global pandemic by the World Health Organization. We believe the ML community can help in an effort for global preparedness by developing computational infrastructure to scale computational drug discovery efforts. Extremely fast and high throughput computational techniques can be used at the beginning of pandemics or even as surveillance systems are screening billions of compounds against entire protemes to find the most promising leads.

In response to the COVID-19 pandemic, scientists across the globe began a massive drug discovery effort spanning traditional targeted combinatorial library screening \cite{sepay2021anti, kong2020covid, clyde2021high}, drug repurposing screens \cite{abo2020molecular}, and crowd-sourced community screening \cite{achdout2020covid}. In common to all these efforts was the ability to leverage off-the-shelf molecular docking programs rapidly. Molecular docking programs are essential preliminary drug discovery efforts as they predict the structure drug candidates form in complex with the protein targets of focus.

The first phase of computational drug discovery often starts with identifying regions of a protein (receptors) that are reasonable targets for small molecule binding, followed by searching small molecules for their ability to bind the receptor. These receptor-binding computations are performed using standard docking software, such as DOCK , AutoDock , UCSF Dock , and many others \cite{wang2016comprehensive}. These molecular docking programs perform a search in which the conformation of the ligand in relation to the receptor (pose) is evaluated recursively until the convergence to the minimum energy is reached. An affinity score, $\Delta$G (U total in kcal/mol), is then used to rank candidate poses. Protein-ligand docking software has two main outputs: a pose (ligand conformer in a particular complex with the protein and an associated score. Scores are approximates to the binding free energy, though the units and interpretation of scores depend on the exact docking protocol used. The inputs are an ensemble of 3D ligand conformations and a prepared protein receptor. We call a protein receptor a protein structure that has been annotated with a particular binding pocket coordinates (i.e. the binding site box).

Several open databases of commercially available compounds, not yet commercially available but with synthesizable properties, and compounds that have an unknown synthesizability are available to the scientific community. Some of the largest include ZINC , Enamine Real , GDB-13  and SAVI  with each containing on order $10^9$ compounds. Recently, \textit{Babuji} released an aggregate collection of over $10^9$ compounds in representations suitable for deep learning \cite{babuji2020targeting}. These representations included drug descriptors, fingerprints, images and canonical smiles. Searching collections of this size using traditional docking tools is not practical.

In this work, we demonstrate deep learning to accelerate docking, thereby expanding our ability search space of small molecules. We do so by first illustrating our large-scale effort of training surrogate models across the SARS-CoV-2 proteome. We illustrate the speed of ML-based surrogate docking models by performing inference on over 1 billion compounds in hours on supercomputer GPU nodes. ML-based surrogate docking models can be combined with classical docking through a prefilter regime where molecules are first screened with the fast surrogate model and the top scoring products are passed to traditional docking. We call this technique surrogate model prefilter then dock (SPFD). We perform a detailed analysis of SPFD by combining model accuracy performance, computational throughput, and virtual screening detection sensitive into a single analysis framework. Our analysis framework indicates our workflow is 10x as fast as the traditional only docking method with less than 0.1\% error of detecting top hits. The framework indicates our ability to improve this speedup is mainly limited by surrogate model accuracy rather than computing power. Our hope is this benchmark allows the community to make progress on the modeling aspect, leading to tangible 100x or even 1000x improvement in throughput for virtual ligand screening campaigns.

\section{Dataset Overview}
As part of our drug discovery campaign for SARS-CoV-2 \cite{clyde2021high}, we developed a database of docked protein-ligand across 15 protein targets and 12M compounds as well as the complexes' associated scores (see SI.\ref{si:docking} for docking pipeline details).

The database contains two related tasks. The first is predicting a ligand's docking score to a receptor based on 2D structural information from the ligand, and the second is a pan-receptor model which encodes the protein target in order to use a single model across different ligands and different targets.  These tasks are distinct from other drug discovery datasets as this benchmark is focused directly on surrogate model performance over the baseline computational drug discovery method of docking. A different approach to applying machine learning to docking is the use ML models a scoring function rather than the result of the optimization of the ligand conformation/position relative to the scoring function \cite{ragoza2017protein}. Other benchmarks are available to address to gap between docking and experimental binding free energy calculations such as DUD-E \cite{mysinger2012directory}.

The data has two modes of representation, either 2D or 3D, where the 3D data a ligand conformation in an SDF file. 2D ligand data is available in a CSV file containing the molecule's purchasable name, a SMILES string, and it's associated docking score in a particular complex.

The ligands available for each dataset are sorted into three categories ORD (orderable compounds from Mcule \cite{mcule}), ORZ (orderable compounds from Zinc \cite{sterling2015zinc}), and an aggregate collection which contains all the available compounds plus others (Drug Bank \cite{wishart2018drugbank}, and Enamine Hit Locator Library \cite{EDB-web}). Docking failures were treated as omissions in the data, which may be important consideration though typically the number of omissions accounts for 1-2\% at most of each sample.

The data is available here, \url{https://doi.org/10.26311/BFKY-EX6P}, and more information regarding persistence and usage is available on the data website and SI \cite{clydedata}.

\section{Method}\label{sec:method}
At a high level, the use of surrogate models for protein-docking ligand aims to accelerate virtual ligand screening campaigns. A surrogate model aims to replace the CPU-bound docking program with a trained model to filter incoming ligands such that the number of actual docking calculations is minimized and the number of missed compounds due to model accuracy is minimized. In other words, a surrogate model is trained and a cut-off is specific, say 1\%. The model is run over the proposed library to screen, and the top 1\% of ligands are then docked utilizing the program in order to have the exact scores and pose information as if typically docking. In this way, we do not see current surrogate models as a replacement for docking, but rather as a means of expanding it's use over large virtual libraries. We call this model SPFD, surrogate model prefilter then dock. This model has a single hyperparameter, $\sigma$ which determines after running the surrogate model over the library, what percentage of most promising predicted compounds do we actually dock utilizing traditional docking techniques.

\subsection{Data frame construction}
We used the protein-ligand docking results between the prepared receptors and compounds in the OZD library to build machine learning (ML) data-frames for each binding site. The raw docking scores (the minimum Chemgauss4 score over the ensemble of conformers in a ligand-receptor docking simulation) were processed. Because we were interested in determining strong binding molecules (low scores), we clipped all positive values to zero. Then, since we used the ReLu activation function at the output layer of the deep neural network, we transformed the values to positive by taking the absolute value of the scores.  The processed docking scores for each compound to each binding site then served as the prediction target. The code for model training can be found here: \url{https://github.com/2019-ncovgroup/ML-Code}.

The features used to train the models were computed molecular descriptors. The molecular descriptors were computed as described by Babuji et al [2020].  The full set of molecular features contains 2-D and 3-D descriptors, resulting in a total of 1,826 descriptors. The approximately 6 million docking scores per receptor and 1,826 descriptors were then joined into a data frame for each receptor.
\subsection{Learning curves}
We performed learning curve analysis with the 3CLPro receptor. A subset of 2M samples were obtained from the full set of 6M samples. The 2M sample dataset was split into train (0.8), validation (0.1), and test (0.1) sets. We trained the deep neural network on subsets of training samples, starting from 100K and linearly increasing to 1.6M samples (i.e., 80\% of the full 2M set). Each model was trained from scratch and we used the validation set to trigger early stopping and the test to calculate measures of generalization performance such as the mean absolute error of predictions.
\subsection{Model Details}
The model was a fully connected deep neural network with four hidden layers each followed by a dropout layer. The dropout rate was set to 0.1. Layer activation was done using the rectified linear unit activation function. The number of samples per gradient update (batch size) was set to 32. The model was compiled using mean squared error as the loss function and stochastic gradient descent with an initial learning rate of 0.0001 and momentum set to 0.9 as the optimizer. The implementation was python using Keras.

The model was trained by setting the initial number of epochs to 400. A learning rate scheduler was used that monitored the validation loss and reduced the learning rate when learning stagnated. The number of epochs with no improvement after which the learning rate was reduced (patience) was set to 20. The factor by which the learning rate will be reduced was set to 0.75 and the minimum allowable learning rate was set to 0.000000001. Early stopping was used to terminate training If after 100 epochs the validation loss did not improve.

Features were standardized by removing the mean and scaling to unit variance before the onset of training using the entire data frame (before the data frame was split into train and test partitions). The train and test partitions were based on a random 80:20 split of the input data frame. Hyperparameter optimization was performed (see SI.\ref{hyperopt}).

Inferencing was performed on Summit. The input was converted to Feather files using the python package feather, which is a wrapper around pyarrow.feather (see the Apache Arrow project at apache.org). Feather formatted files as input in our experience are read faster from disk than parquet, pickle, and comma separated value formats.

\section{Results}

\subsection{Identification of protein targets and binding receptors}
A total of thirty one receptors representing 9  SARS-CoV-2 protein conformations were prepared for docking. These are illustrated in Figure 2 and listed in Table 1. The quality of the receptors reflect what was available at the time the receptor was prepared. For example, whereas the NSP13 (helicase) structure in Table 1was based on homology modeling, today there exists x-ray diffraction models.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.49\textwidth]{images/fig3.png}
    \includegraphics[width=0.45\textwidth]{images/fig5.png}
    \caption{(left) \textbf{Histogram of protein-ligand docking of transformed docking scores for 3CL-M\textsubscript{pro}}. The distribution is from the ORZ dataset based on the transformed 2D scores. (right) \textbf{Learning curve between dataset size and MAE between random and flattened datasets.}}
    \label{fig:3}
\end{figure}

\subsection{Generation of training data}
Compounds in the OZD library were docked as described in SI~\ref{si:libraries}.
The results for the 3CLPro receptor demonstrate a normal distribution (Fig.~\ref{fig:3}). The best docking scores would be in the range of 12 to 18. The distribution of docking scores for the 3CLPro receptor is illustrative of the distributions for all the other receptors. As shown in the figure, there are very few samples with good docking scores relative to the entire set of samples.

\subsection{Sampling comparisons}
We constructed a set of data frames to investigate the impact of the number of samples, sampling approach, and the choice of drug descriptors as features. The number of samples was further investigated using learning curves. Because we are interested in predicting docking scores in the tail of the distribution where the best docking scores exist, we explored two sampling approaches. Lastly, we investigated the impact of using the Mordred 3-D descriptors as features of the compounds.

\begin{table}[h]
\centering
\begin{tabular}{@{}llll@{}}
\toprule
Dataset      & Count (samples) & Sampling method & Distribution (approximate) \\ \midrule
100K-random  & 100,000         & Random          & Normal                     \\
100K-flatten & 100,000         & Flatten         & Uniform                    \\
1M-random    & 1,000,000       & Random          & Normal                     \\
1M-flatten   & 1,000,000       & Flatten         & Uniform                    \\ \bottomrule
\end{tabular}
\caption{The four sampling approaches used to subset the approx. 6M docking scores for OZD.}
\label{tab:sampling}
\end{table}

We generated a dataset subset by sampling the approximately 6M samples in the OZD data complete data-frames. We examined four sampling approaches, differing by two parameters, as listed in Table~\ref{tab:sampling}: (1) the total number of samples drawn from the entire dataset (i.e., the count), and (2) the algorithm used to draw the samples (i.e., the sampling method).

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{images/fig4.png}
    \caption{Docking score histograms for each of the four sampling a) 100K-random, b) 100K-flatten, c) 1M-random and d) 1M-flatten approaches used to generate a subset by sampling the full dataset of available scores (approximately six million samples).}
    \label{fig:sampling}
\end{figure}

Drawing samples at random preserves the original normal-shaped distribution (thus, the name Random). Alternatively, for a more balanced dataset, we sample scores with an alternative algorithm to create a roughly flattened, uniform-like distribution. To include the highly significant, top score samples, we retain the top ten thousand binding ligands. Figure~\ref{fig:sampling} shows the histograms of the docking scores subset with each of the four sampling scenarios for 3CL-M\textsubscript{pro}. The top ten thousand binding ligands are indicated in red. Note that the distribution of the full dataset can be roughly modeled as a normal distribution, as shown in Figure~\ref{fig:3}.

When examining the impact of including the Mordred 3-D descriptors in the feature set, we average the validation loss, validation MAE, and validation $r^2$ across the 31 models as we are interested in the aggregate performance of the models across the 31 receptors.  Our analysis of the inclusion of the Mordred 3D descriptors is presented in SI Table~\ref{tab:samplingresults}. Our results show no significant advantage to including the 3D descriptors. The results  show small improvements in the validation loss across all training data frames when using only the 2D descriptors, and the results are mixed when considering validation $r^2$, with two smaller data frames performing slightly better, and the two larger data frames performing slightly worse. While we do not consider the differences in most cases to be significant, we do demonstrate that adding the extra training parameters in the form of 3D descriptors does not improve the training performance of the model.

When examining the impact of both the training set size (1M or 100K) and sample selection from either a random distribution or flattened distribution, we average the validation loss, validation MAE, and validation $r^2$ for each trained receptor model that represents one of the thirty one different protein pockets. SI Table~\ref{tab:samplingresults} shows the differences between the means. A negative value for the validation loss and validation MAE differences would indicate 1M samples achieved a higher quality model, and a positive value for the validation $r^2$ difference would indicate 1M samples achieved a higher quality model. The results indicate that 1M samples from a flattened distribution perform better than 100K samples for all three metrics, whereas 1M samples from a random distribution achieved better metrics for the validation loss and MAE, however the 1M samples from a random distribution had a lower validation $r^2$.

To better understand the differences between the 1M data sets, the Pearson correlation coefficient was calculated between predicted and the observed values from the validation set for each pocket model. In the case of the v5.1-1M samples, the validation set had 200,000 samples. The mean of the PCC across the set of pocket models was calculated for each 1M data set and the V5.1-1M-random is 0.853 and the V5.1-1M-flatten is 0.914.

\subsection{Learning curve analysis}
To further explore the optimal sample size, we generated learning curves for the 3CLPro receptor model and assume 3CLPro will be indicative of other receptors. Using the entire dataset, which contains approximately 6M samples, imposes a significant computational burden for training a deep neural network model for each receptor and performing HP tuning.  Regardless of the learning algorithm, supervised learning models are expected to improve generalization performance with increasing amounts of high-quality labeled data. At a certain sample size, however, the rate of model improvement achieved by adding more samples significantly diminishes. The trajectory of generalization performance as a function training set size can be estimated using empirical learning curves.

The range at which the learning curve starts to converge indicates the sample size at which the model starts to exhaust its learning capacity. Figure~\ref{fig:3} shows the learning curve where the mean absolute error of predictions is plotted versus the training set size. The curve starts to converge at approximately 1M samples, implying that increasing sample size beyond this range is not expected to improve predictions.

\subsection{Model Accuracy}
The observations of the data frame comparisons (sec. 4.2.3) and learning curves (sec. 4.2.4), specifically:
\begin{itemize}
    \item The 1613 MOrdred descriptors performed better than the 1826 features in most cases
    \item The 1M data frames performed better than the 100K data frames in most cases
    \item The mean Pearson correlation coefficient of the 1M-flatten was higher than that of the 1M-random data frame
\end{itemize}
led to the selection of training data that was used to create the models listed in SI table 4. These models were then used in inferencing to help guide the selection of compounds for use in whole cell viral inhibition assays.
\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{images/fig6.png}
    \caption{Comparison of the 31 receptor models with the 2000 best scoring compounds from ORD and ORZ.}
    \label{fig:topmeans}
\end{figure}
\subsection{Inference across 3.8 billion compounds}
We divided the 4 billion compounds into 4 input data sets to enable better utilization of resources. ENA, G13, ZIN, OTH. We also constructed a set of compounds from the MCULE data set that could be easily purchased (organic synthesis already done). The MCULE subset was named ORD. The inferencing rate was approximately 50,000 samples per second per GPU, and all 6 GPUs per summit node were used.

We analyzed the results for each receptor by selecting the top 2000 scoring compounds and computing mean, standard deviation, maximum, and minimum values. We present two examples of these results in Figure~\ref{fig:topmeans}. It is interesting to note that the range represented by the maximum and minimum predicted scores for the best 2000 scoring compounds is remarkably different between these two, and in fact ZIN was representative of the others (G13, ENA, and OTH). One working hypothesis is that the compounds in ORD have been shown to be synthesizable, whereas compounds in the other sets are not necessarily synthesizable as these are virtual libraries.

\subsection{Deep Learning Inference Speed}
\begin{figure}[t]
    \centering
    \includegraphics[width=0.49\textwidth]{images/RESplotNeurips.png}
    \includegraphics[width=0.49\textwidth]{images/speedup.png}

    \caption{(left) Regression enrichment surface $(n=200,000)$ based on the surrogate model for 7BQY \cite{clyde2020regression}. The $x$-axis represents $\sigma$ which determines the level of filtering the model is used for (i.e., after predicting over the whole library, what percentage of compounds then used in the next stage docking). The $y$-axis is the threshold for determining if a compound is a hit or not. The point $(10^{-1},10^{-3})$ is shaded with 100\% detection. This implies the model over a test set can filter out 90\% of compounds without ever missing a compound with a score in the $10^{-3}$ percentile. In concrete numbers, we can screen 200,000 compounds with the model, take the top 20,000 based on those inference scores, and dock them. The result is running only 20,000 docking calculations, but those would contain near 100\% of the top 200 compounds (as if one docked the entire dataset). (right) Based on equation (1) we compute the relative speedup using surrogate models over traditional workflows with fixed parameters library size (1 billion compounds) and $T_\text{D}=1.37\frac{\text{samples}}{\text{seconds node}}$. The horizontal line indicates where current GPU, surrogate model, throughput is, $T_\text{SPF}$, and the vertical lines correspond to the RES plot values for hit threshold equal to $10^{-3}$. The right-most vertical line implies a VLS campaign with surrogate models where the surrogate GPU-based model can with accuracy $>99\%$ detect the top 10\% from the bottom 90\% implying a 10x speedup over traditional methods. By adding surrogate models as a pre-filter to docking, scientists can dock 10x more in the same amount of time with little detectable loss.}
    \label{fig:my_label}
\end{figure}

We discuss the relative speedup of utilizing a pre-filter surrogate model for
docking campaigns against a traditional docking campaign. We define two
workflows for performing protein-ligand docking over a library of compounds:
D (traditional docking, no surrogate-prefilter) and SPFD (surrogate-prefilter
then dock). We construct timing models of both of these workflows to
understand the relationship between computational accuracy, computational
performance (time and throughput), and pre-filter hyperparameter ($\sigma$). We distinguish between the surrogate model's accuracy, which pertains to how well the surrogate model fits the data, and workflow accuracy, which pertains to how well the results of the whole SPFD workflow compares to the results of just traditional docking workflow.

As discussed in Section \ref{sec:method}, the choice of pre-filter
hyperparameter
is a decision about workflow accuracy for
detecting top leads. Model accuracy influences the workflow accuracy, but the
workflow accuracy can be adjusted with respect to a fixed model accuracy (see
figure~\ref{fig:my_label} where the vertical lines each correspond to the same
underlying model with fixed accuracy but differentiate the overall workflow
accuracy with respect to traditional docking). Therefore we can interpret
$\sigma$ as a trade-off between the workflow throughput and workflow accuracy.
For example, $\sigma=1$ is always 100\% workflow accurate since traditional
docking is run on the whole library when $\sigma=1$, but $\sigma=1$ is even slower than traditional docking as it implies docking the whole library as well as utilizing the surrogate model. We can determine the
overall workflow accuracy with the model by looking at the RES plot, which is
of course has as a factor the performance characteristics of the surrogate
model. Given a particular model's accuracy versus
performance characteristics, different levels of pre-filtering ($\sigma$),
correlate to different tolerances to detecting top-scoring compounds.

For the following analysis, we fix node types for simplicity. Assume the traditional protein-ligand docking software has a throughput $T_{\text{D}}$ in units $\frac{\text{samples}}{\text{(seconds)(node)}}$, and the surrogate models have a throughput $T_{\text{SPF}}$. Let $t_\text{D}$ and $t_\text{SPFD}$ be the wall-clock time of the two workflows. The time of the traditional workflow, $t_\text{D}$, and the time of the surrogate prefilter then dock workflow, $t_\text{SPFD}$, are
\begin{equation}
t_\text{D}=\frac{L}{T_\text{D}}\quad\text{and}\quad t_\text{SPFD} = \frac{\sigma L}{T_\text{D}} + \frac{L}{T_\text{SPF}}.
\end{equation}
Notice, that $t_\text{SPFD}$ is simply the sum of the time of running the surrogate model over the library, $L/T_\text{SPF}$, and the time of traditional docking the highest scoring $\sigma L$ compounds.

\begin{equation}
    \text{Speedup} = \frac{T_\text{SPF}}{T_\text{D} + \sigma T_\text{SPF}}
\end{equation}

This implies that the ideal speedup of our workflow is directly dependent on the throughput of both the docking calculation, surrogate model, and the parameter $\sigma$. $\sigma$ is indirectly dependent on the model accuracy. If the surrogate model was completely inaccurate, even though $\sigma=10^{-3}$ implies a 1000x speed up, no hits would be detected. If one wants to maximize workflow accuracy, that is not miss any high scoring compounds compared to traditional docking, then they must supply a threshold for hits (corresponding to the $y$-axis of RES). Suppose this threshold is $y_\text{thres}=10^{-3}$. If they wanted to maximize not missing any compounds, they should set $\sigma$ to $10^{-1}$ based on this model's RES plot since that is the smallest value of $\sigma$ such that the detection accuracy of the surrogate model is near 100\%. But, it is not always the case downstream tasks require 100\% detection---hence $\sigma$ is a true hyperparameter.

 We infer $T_{\text{SPF}}$ on a Summit (Oak Ridge Leadership Computing Facility) and on an A100 ThetaGPU node (Argonne Leader Computing Facility). Both tests were using 64 nodes, 6 GPUs per node, but the throughput was computed per GPU. We found the V100 summit node was capable of 258.0K$\frac{\text{samples}}{\text{(s)(node)}}$ while the A100 nodes were 713.4K $\frac{\text{samples}}{\text{(s)(node)}}$. We infer $T_\text{D}$ as 1.37 $\frac{\text{samples}}{\text{(s)(node)}}$ based on a CPU docking run over 4,000 summit nodes with 90\% CPU utilization from \cite{clyde2021high}. Thus, we can compute the speedup based on a Summit node head-to-head comparing setting $T_\text{SPF}$ to 258.0K and $T_\text{D}$ to 1.37 in Equation 1, resulting in a speedup of 10X for $\sigma=0.1$. Based on the RES analysis in Figure~\ref{fig:my_label}, $\sigma$ of 0.1 corresponds to a model accuracy of near 100\% for filtering high scoring leads ($>1$\% of library). If one is willing to trade-off some loss of detection, say 70\% detection of high scoring leads, then the speedup is 100X. The extreme case, a choice of $\sigma=10^{-3}$, implies a speedup 1000X but means only roughly 20\% of the top scoring leads may appear at the end.

The analysis of SPFD implies that speedups are essentially determined by $\sigma$ while $T_\text{SPF}$ does not have as large of an effect (this is based on how fast ML inference currently is). As a hyperparameter, $\sigma$ is dependent on the workflow's context and, in particular, what the researchers are after for that virtual ligand screening campaign. We can say, though, at least informally, model accuracy and $\sigma$ are highly related. The more accurate the models are, the better the RES plot gets as one is willing to trust the ML model for filtering the best compounds from the rest. In figure~\ref{fig:my_label}, the $x$-axis of both plots are similar. The accuracy of a particular $\sigma$ is found by setting one's level of desired detection, the $y$-axis of the RES plot, and then checking the $(\sigma, y)$ point to see how accurate the model is there. The choice of $\sigma$ is subjective based on how accurate one needs the model for their $y$-axis threshold for hits. We focus mainly on the case of no loss of detection, which means $\sigma=0.1$ for our particular trained models. In order to focus on the theoretical model of relating computational accuracy, confidence (again, in a colloquial sense), and computational performance, we simplify over a richer model of performance analysis assuming uniformity of task timing and perfect scaling. Furthermore, we chose a head-to-head comparison of a particular node type's CPU performance to GPU performance. At the same time, we could have compared the best non-surrogate model workflow times to the best surrogate model workflow times.

\section{Discussion}
We demonstrate an accelerated protein-ligand docking workflow called surrogate model prefiltering then dock (SPFD) which is at least 10x faster than traditional docking with nearly zero loss of detection power. We utilize neural network models to learn a surrogate to the CPU-bound protein-ligand docking code. The surrogate model has a throughput over six orders of magnitude than the standard docking code. By combining these workflows, utilizing the surrogate model as a prefilter, we can gain a 10x speedup in the traditional docking software without losing any detection ability (for hits defined as the best scoring 0.1\% of a compound library). We utilize regression enrichment surfaces to perform this analysis. The regression enrichment surface plot is more illustrative than the typical accuracy metrics reported from deep learning practices. Figure~\ref{fig:my_label} showcases our initial models at this benchmark show a 10x speedup without loss of detection (or 100x speedup with 70\% detection). This 10x speedup means if a current campaign takes one day to run on library size $L$, in the same amount of time without missing leads, one can screen ten times as many compounds. Given the potential for 100x or even 1000x speedup for docking campaigns, we hope to advance the ability of surrogate models to filter at finer levels of discrimination accurately.

Utilizing SPFD we report a 10x speedup to traditional docking with little to no loss of accuracy or methodical changes needed besides just utilizing ML-models as a prefilter. On the other hand, we see these results as conceptually tying together the model application (hit threshold and adversity to detection loss, choice of $\sigma$) with more traditional analysis such as performance characteristics and model performance evaluation. Our results of the timing analysis also imply a need for more accurate models that are more sensitive to screening large libraries. Our analysis, outlined in Figure~\ref{fig:my_label}, highlight the choice of prefilter threshold as the limiting factor for seeing orders of magnitude speedup. In particular, focusing on speedups which show \textit{no loss} of detection, model accuracy must be pressed forward as there is no path to accelerating traditional docking workflows without more accurate surrogate models. Given out of box modeling technique can speed up virtual screening 10x with no loss of detection power for a reasonable hit labeling strategy (top $0.1$\%), we believe the community is not far from 100x or even 1000x. The way to get there is to boost our model accuracies or develop techniques to recover hits in lossy SPFD regimes (such as not improving model performance but decreasing $\sigma$ to $10^{-2}$ and applying another technique to recover the 30\% loss of detection power). We believe this benchmark is important, as successful early drug discovery efforts are essential to rapidly finding drugs to emerging novel targets. Improvements in this benchmark will lead to orders of magnitude improvement in drug discovery throughput.

\section*{Acknowledgements}
This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.

We acknowledge other members of the National Virtual Biotechnology Laboratory (NVBL) Medical Therapeutics group. We acknowledge computing time via the COVID19 HPC Consortium.

Research was supported by the DOE Office of Science through the National Virtual Biotechnology Laboratory, a consortium of DOE national laboratories focused on response to COVID-19, with funding provided by the Coronavirus CARES Act.  SJ also acknowledges support from ASCR DE-SC0021352.

\textit{Scientific and Technical Information Only For All Information.}\\
Unless otherwise indicated, this information has been authored by an employee or employees of UChicago Argonne, LLC., operator of Argonne National Laboratory, with the U.S. Department of Energy. The U.S. Government has rights to use, reproduce, and distribute this information. The public may copy and use this information without charge, provided that this Notice and any statement of authorship are reproduced on all copies.
While every effort has been made to produce valid data, by using this data, User acknowledges that neither the Government nor operating contractors of the above national laboratories makes any warranty, express or implied, of either the accuracy or completeness of this information or assumes any liability or responsibility for the use of this information. Additionally, this information is provided solely for research purposes and is not provided for purposes of offering medical advice. Accordingly, the U.S. Government and operating contractors of the above national laboratories are not to be liable to any user for any loss or damage, whether in contract, tort (including negligence), breach of statutory duty, or otherwise, even if foreseeable, arising under or in connection with use of or reliance on the content displayed in this report.

\printbibliography

\section{Dataset Details and Use}

\subsection{Persistence, Ethics, and License}

The 2D version of the dataset is hosted publicly by a third party provider, FigShare for redundancy (for which they determine persistence), \url{https://doi.org/10.6084/m9.figshare.14745234}. The rest of the data, including all the 3D structures, is hosted by Argonne's Leadership Computing Center and accessible via a Globus endpoint with documentation hosted by GitHub. The authors are confident the data will be persistent across FigShare, GitHub, ALCF, and Globus.

The authors believe this data presents a minimal ethical risk to the community at large. The proposed dataset contains computer-generated protein-ligand structures and computed scores. Information of this sort, albeit at a smaller scale, is widely available on the web currently, and releasing this particular dataset would not set any new standards (in terms of a qualitative assessment of data type). The authors believe the biggest risk of releasing this dataset would be localized to one's interpretation of resulting models, theories, or endeavours based on inductive reasoning from the data alone---but, this risk is typical of any scientific dataset.

\subsection{Data Generation Methods}
The original data was generated by Open Eye Scientific's FRED docking protocols, and was aggregated, cleaned, and standardized with naming convetions. The code was generated with this software: \url{github.com/inspiremd/Model-generation}.

\subsection{Docking Protocol}\label{si:docking}
The training and testing datasets for these experiments were generated using 31 protein receptors, covering 9 diverse SARS-CoV-2 viral target protein conformations, that target  (1) 3CLPro (main protease, part of the non-structural protein/ NSP-3), (2) papain like protease (PLPro), (3) SARS macrodomain (also referred to as ADP-ribosyltransferase, ADRP), (4) helicase (NSP13), (5) NSP15 (endoribonuclease), (6) RNA dependent RNA polymerase (RDRP, NSP7-8-12 complex), and (7) methyltransferase (NSP10-16 complex). For each of these protein targets, we identified a diverse set of binding sites along the protein interfaces using two strategies: for proteins that had already available structures with bound ligands, we utilized the X-ray crystallography data to identify where ligand densities are found and defined a pocket bound by a rectangular box surrounding that area; and for proteins that did not have ligands bound to them, we used the FPocket toolkit  that allowed us to define a variety of potential binding regions (including protein interfaces) around which we could define a rectangular box. This process allowed us to expand the potential binding sites to include over 90 unique regions for these target proteins. We use the term target to refer to one binding site. The protocol code can be found here: \url{https://github.com/2019-ncovgroup/HTDockingDataInstructions}.

\subsection{Preparation of Ligand Libraries}\label{si:libraries}
Two ligand libraries were prepared. The first was the orderable subset of the Zinc15 database (we refer to this as OZD) and the second was the orderable subset of the MCULE  compound database (we refer to this as ORD). The generation of the orderable subsets was primarily a manual activity that involved finding all compounds that are either in stock or available to ship in three weeks across a range of suppliers. Consistent SMILE strings and drug descriptors for the orderable subsets of the Zinc15 and MCULE compound databases were generated as described by Babuji et al [2020]. Drug descriptors for the Zinc15 and MCULE compound databases  can be downloaded from the nCOV Group Data Repository at \url{https://2019-ncovgroup.github.io}.

\subsection{Docking Protocols}
We used OpenEye Toolkits  for docking six million (6M) small-molecules from the OZD database. For each ligand from the database, we calculate a single Chemgauss4 score as a surrogate for binding affinity. For each ligand in the database (provided as a SMILES string), we create an ensemble of structures, sampled from various proteinization states (tautomers) and 3-D conformers. Typically, this results in  approximately one thousand 3-D structures for each ligand. If the ligand’s stereochemistry is ambiguous from the provided SMILES, enantiomers are enumerated prior to conformer generation. Each conformer is then docked to the protein target using FRED or HYBRID depending on the availability of a bound ligand in the crystal structure of the specific target. Scores are calculated over the best pose over the ligand ensemble using the Chemgauss4 scoring function. The scoring function is unitless and aims to measure the fitness of a ligand pose in an active site via a numerical score. Poses with more negative scores are more likely to be correctly docked. In order to produce a single score for each database SMILES entry, the minimum over the ensemble is returned. The typical range of the scoring function is between -20 and 0, though the scoring function range is unbounded.

\subsection{Data Cleanliness}
Failure analysis of docking runs is available in the SI of \cite{clyde2021high}.

\subsection{Model Hyperparameter Optimization}
\label{hyperopt}
The CANDLE framework was used subsequently used to tune the deep neural network for future training and screening activities. The CANDLE compliant deep neural network was tuned in two phases. The first involved using two CANDLE hyperparameter optimization workflows - mrlMBO and GA. Each differs in the underlying ML techniques used to optimize the hyperparameters. The second phase involved implementing and testing new sample weighting strategies in an attempt to weight the samples at the good end of the distribution more heavily during training. Results of the GA and mlrMBO workflows produced a model architecture that had a 6.6\% decrease in the validation mean absolute error and a 2.8\% increase in the validation R-squared metrics.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.65\textwidth]{images/fig7.png}
    \caption{(left) Effects of sample weighting strategies on the default and optimized model. }
    \label{fig:weight}
\end{figure}

Efforts to decrease the error in the good tail of the distribution (where the docking scores are best) focused on adding sample weights to the model while training. We investigated linear and quadratic weighting strategies. We applied the weighting strategies to both the default model as well as the hyperparameter optimized model. The linear strategy weights the sample proportionally with the docking score, while the quadratic scales with the square of the docking score. These strategies generic in that they can be applied to basically any training target value. To analyze the impact of the weighting strategies, we computed the mean absolute error on bins of predicted scores with a bin interval of one. These results are presented in Figure~\ref{fig:weight}.

\subsection{Model Scores}\label{si:models}
See table~\ref{tab:sisamples}.
\begin{table}[h]
\centering
\begin{tabular}{@{}llllllll@{}}
\toprule
pocket   (model)            & epochs & loss  & mae   & r2    & val\_loss & val\_mae & val\_r2 \\ \midrule
3CLPro\_7BQY\_A\_1\_F       & 513    & 0.338 & 0.454 & 0.870 & 0.426     & 0.505    & 0.838   \\
ADRP\_6W02\_A\_1\_H         & 599    & 1.020 & 0.782 & 0.786 & 1.326     & 0.876    & 0.724   \\
NPRBD\_6VYO\_A\_1\_F        & 453    & 0.302 & 0.427 & 0.848 & 0.356     & 0.466    & 0.822   \\
NPRBD\_6VYO\_AB\_1\_F       & 599    & 0.482 & 0.540 & 0.800 & 0.601     & 0.597    & 0.752   \\
NPRBD\_6VYO\_BC\_1\_F       & 427    & 0.566 & 0.586 & 0.899 & 0.702     & 0.653    & 0.876   \\
NPRBD\_6VYO\_CD\_1\_F       & 523    & 0.474 & 0.534 & 0.816 & 0.602     & 0.597    & 0.769   \\
NPRBD\_6VYO\_DA\_1\_F       & 587    & 0.485 & 0.541 & 0.854 & 0.591     & 0.595    & 0.824   \\
NSP10-16\_6W61\_AB\_1\_F    & 283    & 0.490 & 0.546 & 0.902 & 0.615     & 0.613    & 0.878   \\
NSP10-16\_6W61\_AB\_2\_F    & 387    & 0.523 & 0.565 & 0.901 & 0.655     & 0.628    & 0.877   \\
NSP10\_6W61\_B\_1\_F        & 433    & 0.576 & 0.590 & 0.908 & 0.677     & 0.631    & 0.893   \\
Nsp13.helicase\_m1\_pocket2 & 338    & 0.553 & 0.577 & 0.867 & 0.663     & 0.633    & 0.843   \\
Nsp13.helicase\_m3\_pocket2 & 434    & 0.485 & 0.538 & 0.878 & 0.582     & 0.585    & 0.855   \\
NSP15\_6VWW\_A\_1\_F        & 406    & 0.526 & 0.563 & 0.837 & 0.640     & 0.621    & 0.804   \\
NSP15\_6VWW\_A\_2\_F        & 441    & 0.336 & 0.451 & 0.876 & 0.417     & 0.506    & 0.849   \\
NSP15\_6VWW\_AB\_1\_F       & 599    & 0.234 & 0.376 & 0.829 & 0.298     & 0.423    & 0.784   \\
NSP15\_6W01\_A\_1\_F        & 596    & 0.473 & 0.533 & 0.835 & 0.595     & 0.597    & 0.795   \\
NSP15\_6W01\_A\_2\_F        & 451    & 0.313 & 0.434 & 0.888 & 0.378     & 0.475    & 0.865   \\
NSP15\_6W01\_A\_3\_H        & 530    & 0.759 & 0.679 & 0.784 & 0.967     & 0.754    & 0.727   \\
NSP15\_6W01\_AB\_1\_F       & 470    & 0.261 & 0.396 & 0.829 & 0.316     & 0.434    & 0.796   \\
NSP16\_6W61\_A\_1\_H        & 583    & 1.044 & 0.795 & 0.787 & 1.339     & 0.888    & 0.728   \\
PLPro\_6W9C\_A\_2\_F        & 512    & 0.343 & 0.458 & 0.858 & 0.427     & 0.508    & 0.825   \\
RDRP\_6M71\_A\_2\_F         & 461    & 0.311 & 0.430 & 0.855 & 0.384     & 0.479    & 0.823   \\
RDRP\_6M71\_A\_3\_F         & 498    & 0.495 & 0.548 & 0.859 & 0.599     & 0.602    & 0.830   \\
RDRP\_6M71\_A\_4\_F         & 463    & 0.382 & 0.481 & 0.837 & 0.465     & 0.528    & 0.803   \\
RDRP\_7BV1\_A\_1\_F         & 394    & 0.312 & 0.433 & 0.853 & 0.378     & 0.477    & 0.823   \\
RDRP\_7BV1\_A\_2\_F         & 531    & 0.497 & 0.546 & 0.848 & 0.603     & 0.597    & 0.817   \\
RDRP\_7BV1\_A\_3\_F         & 451    & 0.453 & 0.524 & 0.849 & 0.550     & 0.583    & 0.818   \\
RDRP\_7BV1\_A\_4\_F         & 420    & 0.385 & 0.481 & 0.873 & 0.476     & 0.536    & 0.844   \\
RDRP\_7BV2\_A\_1\_F         & 589    & 0.304 & 0.428 & 0.821 & 0.369     & 0.469    & 0.784   \\
RDRP\_7BV2\_A\_2\_F         & 422    & 0.441 & 0.516 & 0.839 & 0.562     & 0.581    & 0.798   \\
RDRP\_7BV2\_A\_3\_F         & 510    & 0.466 & 0.531 & 0.830 & 0.579     & 0.590    & 0.791   \\ \bottomrule
\end{tabular}
\caption{Table of released model's training details and validation scores. Released model files and corresponding code is avilable from the project GitHub.}
\label{tab:sisamples}
\end{table}

\subsection{Modeling Feature Details}
\begin{table}[t]
\centering
\begin{tabular}{@{}lllll@{}}
\multicolumn{5}{l}{1613 Features}                                                             \\ \midrule
\textbf{Model}      & \textbf{epoch} & \textbf{val loss} & \textbf{val MAE} & \textbf{val $\bm{r^2}$} \\ \midrule
V5.1-100K-flatten-2 & 337            & 0.80              & 0.66             & 0.71            \\
V5.1-100K-random-2  & 336            & 0.80              & 0.66             & 0.71            \\
V5.1-1M-flatten-2   & 484            & 0.60              & 0.59             & 0.81            \\
V5.1-1M-random-2    & 455            & 0.49              & 0.52             & 0.68            \\ \midrule
\multicolumn{5}{l}{}                                                                          \\
\multicolumn{5}{l}{1826 Features}                                                             \\ \midrule
\textbf{Model}      & \textbf{epoch} & \textbf{val loss} & \textbf{val MAE} & \textbf{val $\bm{r^2}$} \\ \midrule
V5.1-100K-flatten-2 & 313            & 0.97              & 0.74             & 0.85            \\
V5.1-100K-random-2  & 330            & 0.81              & 0.67             & 0.71            \\
V5.1-1M-flatten-2   & 462            & 0.60              & 0.59             & 0.81            \\
V5.1-1M-random-2    & 456            & 0.52              & 0.54             & 0.67            \\ \bottomrule
\end{tabular}
\caption{Impact of including Mordred 3-D descriptors in the training data for the different sampling strategies.}
\label{tab:samplingresults}
\end{table}

\begin{table}[t]
\centering
\begin{tabular}{@{}lllll@{}}
\toprule
\textbf{Sample Selection Strategy}                      & \textbf{epoch}   & \textbf{val loss}   & \textbf{val mae}   & \textbf{val r2}   \\ \midrule
V5.1-100K-flatten    & 337    & 0.80    & 0.66    & 0.71    \\
V5.1-100K-random     & 336    & 0.80    & 0.66    & 0.71    \\
V5.1-1M-flatten      & 484    & 0.60    & 0.59    & 0.81    \\
V5.1-1M-random       & 455    & 0.49    & 0.52    & 0.68    \\
\multicolumn{5}{l}{}                                        \\ \midrule
\multicolumn{5}{l}{Difference of 1M Samples - 100K Samples} \\ \midrule
\multicolumn{1}{l|}{\textbf{Sample Selection Strategy}} & \textbf{$\bm{\Delta}$ epoch} & \textbf{$\bm{\Delta}$ val loss} & \textbf{$\bm{\Delta}$ val mae} & \textbf{$\bm{\Delta}$ val r2} \\
flatten              & 147    & -0.20   & -0.07   & 0.11    \\
random               & 119    & -0.31   & -0.14   & -0.03   \\ \bottomrule
\end{tabular}
\caption{Comparison of 1M samples to 100K samples.}
\label{tab:sisampling}
\end{table}

\end{document}


\section{Introduction}
Let $G$ be a simple graph and let $E(G), V(G)$ denote its edge and vertex set respectively. For
$v \in V(G)$, denote by $N_G(v)$ the set of neighbors of $v$ in $G$, and let $d_G(v) = |N_G(v)|$ be the degree
of $v$ in $G$. For two disjoint subgraphs $H_1, H_2$ of $G$, $N_{H_1}(H_2)$ denotes the set of neighbors of vertices of $H_2$ in $G$
that are contained in $V(H_1)$.
For any subset $S \subseteq V(G)$,  $G[S]$ is the subgraph of $G$ induced on $S$,
$G - S$ denotes the subgraph $G[V(G) \setminus S]$, and $N(S) = \cup_{v \in S} N_G(v)$.
For any disjoint $A, B \subseteq V(G)$,  $E_G(A, B)$ is  the set of all edges  with one end-vertex
in $A$ and the other end-vertex in $B$.  If $u$ and $v$ are adjacent in $G$, we write $u\sim v$.

We say that $G$ is hamiltonian if there exists a cycle which contains every vertex of $G$.
We say that $G$ is $H$-free if there does not exist an induced copy of $H$ in $G$. Denote by $c(G)$ the number of components of $G$.
Let $t \geq 0$ be a real number. We say a graph $G$ is $t$-tough if for each cut set $S$ of $G$ we have $t \cdot c(G-S) \leq |S|$. The toughness of a graph $G$, denoted $\tau(G)$, is the maximum value of $t$ for which $G$ is $t$-tough if $G$ is non-complete, and is defined to be $\infty$ if $G$ is complete. Chv\'{a}tal introduced the notion of toughness in his 1973 paper \cite{MR316301}, where he also conjectured the existence of a constant $t$ such that every $t$-tough graph on at least three vertices is hamiltonian. The conjecture has been verified for certain special classes of graphs, but remains open in general. Recent work has proven the conjecture for $2K_2$-free graphs \cite{broersma2k2, shan2k2, ota2021toughness}, $(P_2 \cup P_3)$-free graphs \cite{shanp2up3}, $(K_2 \cup 2K_1)$-free graphs \cite{1tough}, and planar chordal graphs. We refer the reader to \cite{bauersurvey2} for a survey on more related results.

In this paper, we support Chv\'{a}tal's conjecture by proving the following result:

\begin{THM}\label{theorem1}
If $G$ is a 3-tough $(K_2 \cup 3K_1)$-free graph on at least 3 vertices, then $G$ is hamiltonian.
\end{THM}

The remainder of this paper is organized as follows: in Section \ref{preliminaries}, we discuss results necessary for the proof of Theorem \ref{theorem1} and in Section \ref{section3} we prove Theorem \ref{theorem1}.

\section{Preliminaries}\label{preliminaries}

In this section, we give results necessary to complete the proof of Theorem \ref{theorem1}.

\begin{LEM}[Dirac \cite{hamconn1}, Ore \cite{hamconn2}]\label{1.5}
    Let $G$ be a graph on $n$ vertices such that $\delta(G) \geq \frac{n+1}{2}$. Then $G$ is hamiltonian-connected.
\end{LEM}

\begin{LEM}[Bauer et al. \cite{MR1336668}]\label{bauer}
 Let $G$ be a $t$-tough graph on $n \geq 3$ vertices with $\delta(G) > n/(t + 1) - 1$. Then $G$ is hamiltonian.
\end{LEM}

\begin{LEM}[Li et al. \cite{1tough}]\label{k2cup2k1}
    Let $R$ be an induced subgraph of $P_4$, $K_1 \cup P_3$ or $K_2 \cup 2K_1$. Then
every $R$-free 1-tough graph on at least three vertices is hamiltonian.
\end{LEM}

The following lemma is a consequence of Menger's theorem, which can be found in \cite{bondy}. For a positive integer $k$, define
$[1,k] = \{1,2,\cdots, k\}$.

\begin{LEM}\label{3.2.6}
Let $G$ be a $k$-connected graph and $X_1,X_2$ be distinct subsets of $V(G)$.
Then there exist $k$ internally disjoint paths $P_1, \ldots, P_k$ such that
\begin{enumerate}[(a)]
    \item $|V(P_i) \cap X_1| = |V(P_i) \cap X_2| = 1$, and $P_i$ is internally disjoint from each $X_1$ and $X_2$.
    \item if $|X_i| \geq k$ for some $i \in [1,2]$, then $V(P_j) \cap X_i \ne V(P_\ell) \cap X_i$ for all distinct $j, \ell \in [1,k]$.
    \item if $|X_i| < k$ for some $i \in [1,2]$, then every vertex of $X_i$ is an end-vertex of some path $P_j$ for $j \in [1,k]$.
\end{enumerate}
\end{LEM}
The following lemma provides some structural properties of $(K_2 \cup 3K_1)$-free graphs.

\begin{LEM}\label{lemma:cut-structure}
Let $G$ be a connected ($K_2 \cup 3K_1$)-free graph, and $S
\subseteq V(G)$ be a cut set such that $G - S$ has at least three components. Then we have
the following statements:
\begin{enumerate}[(a)]
	\item If $G - S$ has a nontrivial component, then $G - S$ has exactly three components. \label{parta}
	\item If $G - S$ has a nontrivial component, then the component is ($K_2 \cup K_1$)-free.
\end{enumerate}
\end{LEM}

\textbf{Proof}.\quad   For  part $(a)$,
let $D$ denote a nontrivial component of $G - S$. Assume for the sake of contradiction that $G - S$ has more than three  components.
Taking an edge from $D$ and a single vertex from three other components, respectively, gives an induced copy of $K_2 \cup 3K_1$.
This gives a contradiction to the ($K_2 \cup 3K_1$)-freeness of $G$.

For part $(b)$,      $G - S$ must have exactly three  components by part (a). Assume for the sake of contradiction
    that the nontrivial component  is not $(K_2 \cup K_1)$-free.
    Then taking an induced copy of $K_2 \cup K_1$ from this component and one vertex each from the other two components gives
     an induced copy of $K_2 \cup 3K_1$, which contradicts
the ($K_2 \cup 3K_1$)-freeness of $G$.
\qed
\\

Note that in any $(K_2 \cup 3K_1)$-free graph $G$, the components yielded by any cut set $S$ such that $c(G-S) \geq 3$ must be $(K_2 \cup K_1)$-free. The following lemmas deal with the structure of $(K_2 \cup K_1)$-free graphs. Lemma $\ref{lemma:trivialcomponents}$ is
used in the proof of Lemma \ref{mindegreeD}.

\begin{LEM}\label{lemma:trivialcomponents}
If $G$ is a ($K_2 \cup K_1$)-free graph and $S$ is a cut set of $G$, then
every component of $G - S$ is trivial.
\end{LEM}

\textbf{Proof}.\quad Assume there exists some nontrivial component of $G - S$. Since $S$ is a cut set,
it must disconnect $G$ into at least two components. Taking an edge from a nontrivial component
and a vertex from another component gives an induced copy of $K_2 \cup K_1$, contradicting the
$(K_2 \cup K_1)$-freeness of $G$.
\qed

The independence number of a graph $G$, denoted $\alpha (G)$, is the size of a largest independent set of $G$.

\begin{LEM}\label{mindegreeD}
    Let $t > 0$ be real and $G$ be a $(K_2 \cup K_1)$-free graph on $n$ vertices with $\alpha(G) \leq \frac{n}{t+1}$. Then $\delta(G) \geq n - \frac{n}{t+1}$.
\end{LEM}

\textbf{Proof}.\quad
Assume  $\delta(G) < n-\frac{n}{t+1}$.
Let $v\in V(G)$ with $d_G(v)=\delta(G)$, and let
$W=V(G)\setminus N_G(v)$.  Then $|W| =|V(G)|-\delta(G) > \frac{n}{t+1}$.
As $G$ is $(K_2\cup K_1)$-free, and $N_G(v)$
is a cut set of $G$, every component of $G-N_G(v)$ is trivial by Lemma~\ref{lemma:trivialcomponents}. Since $W=G-N_G(v)$,  $W$ is an independent set of $G$. However, $|W| > \frac{n}{t+1} \geq \alpha(G)$, giving a contradiction.
\qed

\section{Proof of Theorem \ref{theorem1}}\label{section3}
For any $u,v \in V(G)$, we call a path $P$ connecting $u$ and $v$ a $(u,v)$-path.
Let $P$ be an $(x,y)$-path and $Q$ be a $(y,z)$-path such that $V(P) \cap V(Q) = \{y\}$. Then $xPyQz$ denotes a path from $x$ to $z$.
If $P$ is an $(x,y)$-path and $Q$ is a disjoint $(z, w)$-path where $y \sim z$, we denote by $xPyzQw$ the concatenation
of $P$ and $Q$ through the edge $yz$.
Let $\overset{\rightharpoonup }{C}$
be a cycle with some fixed orientation. For any $u,v \in V(C)$, we denote by $u\overset{\rightharpoonup }{C} v$ the path from $u$ to $v$ following the orientation of $C$. Similarly, we denote by $u \overset{\leftharpoonup }{C} v$ the inverse path from $u$ to $v$.
The immediate successor of $u$ on $\overset{\rightharpoonup }{C}$ is denoted by $u^+$.

\proof[Proof of Theorem~\ref{theorem1}]
Let $G$ be a 3-tough $(K_2 \cup 3K_1)$-free graph. We may assume that $G$ is not complete, otherwise there exists a hamiltonian cycle. Therefore, $G$ is 6-connected. By Theorem \ref{bauer},
we may assume $\delta(G) \leq \frac{n}{4} -1$. Since $\delta(G) \geq 6$, we get $n \geq 28$.
Let $C$ be a longest cycle of $G$.

\begin{CLAIM}\label{claim1}
$|V(C)| \geq \frac{3n}{4}$.
\end{CLAIM}

\textbf{Proof}.\quad We assume first that there exist $u,v \in V(G)$ with $u \not \sim v$ such that $|N(u) \cup N(v)| \leq \frac{n}{4}$. Let $S = N(u) \cup N(v)$.

Note that the components of $G - S$ cannot all be trivial, as this would imply $\frac{|S|}{c(G-S)} \leq \frac{\frac{n}{4}}{\frac{3n}{4}} < 3$, but
$G$ is 3-tough. Thus, there must exist a nontrivial component $D$ in $G - S$. Since each of $u,v$ are components of $G - S$ and $G - S$ has a nontrivial component, it follows that $c(G - S) \geq 3$. Thus $G - S$ has exactly 3 components by Lemma \ref{lemma:cut-structure}.
Since $G$ is 3-tough, it follows that $\alpha(G) \leq \frac{n}{4}$.

Also, we know that $\delta(D) \geq |V(D)| - \frac{n}{4} \geq \frac{|V(D)| + 1}{2}$, as $n \geq 12$ and $|V(D)| \geq \frac{3n}{4}-2$ by Lemma \ref{mindegreeD}. Thus $D$ is hamiltonian-connected by Lemma \ref{1.5}. Since $G$ is 2-connected, by Lemma \ref{3.2.6} we can find in $G$ two disjoint paths $P_1$ from $u$ to some $x_1 \in V(D)$, and $P_2$ from $u$ to some $x_2 \in V(D)$ where $x_1 \ne x_2$, and each $P_i$ is internally disjoint from $D$. As $D$ is hamiltonian-connected, we can find a hamiltonian path $Q$ in $D$ from $x_1$ to $x_2$.
Define $C' = uP_1x_1Qx_2P_2u$. Then $|V(C')| \geq |V(D)| + 3 \geq \frac{3n}{4} - 2 + 3 = \frac{3n}{4}+1 > \frac{3n}{4}$.
As $C$ is a longest cycle of $G$, $|V(C)| \geq |V(C')| \geq \frac{3n}{4}$.

We then assume that for any $u,v \in V(G)$ with $u \not\sim v$, it holds that
$|N(u) \cup N(v)| > \frac{n}{4}$. Note that $\delta(G) \leq \frac{n}{4} - 1$ by our earlier assumption. Let $u \in V(G)$ with $d_G(u) = \delta(G)$.
Define $G_1 = G - (N_G(u) \cup \{u\})$.
We know that $G_1$ is $(K_2 \cup 2K_1)$-free, as there are no edges between $G_1$ and $u$ and the original graph is $(K_2 \cup 3K_1)$-free.
If $G_1$ is 1-tough, then it has a hamiltonian cycle by Lemma \ref{k2cup2k1} which has at least $\frac{3n}{4}$ vertices, thus
$|V(C)| \geq \frac{3n}{4}$. Thus we may assume $G_1$ is not 1-tough. Let $W$ be a tough set of $G_1$, i.e. $W$ is a cut set of $G_1$ such that $\frac{|W|}{c(G_1-W)} = \tau(G_1)$.

We claim below that $c(G_1-W) = 2$.
We first note that $G_1 - W$ has at least one nontrivial component. Otherwise, all components of $G_1 - W$ are trivial,
i.e. $c(G_1 - W) = |V(G_1)| - |W|$. Since $\tau(G_1) < 1$, we have $|W| < c(G_1 - W)$. Then $c(G_1 - W) > \frac{1}{2}|V(G_1)| \geq
\frac{1}{2}(n - \frac{n}{4}) = \frac{3n}{8}$. Let $S = N(u) \cup \{u\} \cup W$. Then $c(G - S) = c(G_1-W)$ and $\frac{|S|}{c(G-S)} = \frac{|S|}{c(G_1-W)} \leq \frac{5n/8}{3n/8} < 3$, which is a contradiction as $G$ is 3-tough. Thus $G_1 - W$ has a nontrivial component.

Next, assume that $G_1 - W$ has more than 2 components. Then by the argument above, at least one of the components is nontrivial. Taking an edge from a nontrivial component
and a vertex from each of the two others gives a copy of $K_2 \cup 2K_1$, which contradicts the $(K_2 \cup 2K_1)$-freeness of $G_1$. Thus, we must have
$c(G_1 - W) = 2$.

Let $D_1, D_2$ be the two components of $G_1 - W$.
Then as $G_1$ is not 1-tough, we get $|W| = 1$. Note that $|V(D_i)| \geq 2$. Otherwise, say $|V(D_2)| = 1$ and let $V(D_2) = \{v\}$, then $|N(u) \cup N(v)| \leq \frac{n}{4} - 1 + 1 = \frac{n}{4}$ which contradicts the assumption that
$|N(u) \cup N(v)| > \frac{n}{4}$ for any nonadjacent $u,v \in V(G)$.
Since both $D_1$ and $D_2$ are nontrivial and $G_1$ is $(K_2 \cup 2K_1)$-free, each $D_i$ is a complete graph.

Then by Lemma \ref{3.2.6}, we can find in $G$ two disjoint paths $P_1$ from some $x_1 \in V(D_1)$ to some $x_2 \in V(D_2)$, and $P_2$ from some $y_1 \in V(D_1)$ to some $y_2 \in V(D_2)$ where $x_1 \ne y_1, x_2 \ne y_2$. As each $D_i$ is complete, we can find a hamiltonian path $Q_i$ in $D_i$ from $x_i$ to $y_i$. Then the cycle
$$C' = x_1P_1x_2Q_2y_2P_2y_1Q_1x_1$$
satisfies $|V(C')| \geq |V(D_1)| + |V(D_2)| + 2 \geq \frac{3n}{4}$. As $C$ is a longest cycle, we have $|V(C)| \geq |V(C')| \geq \frac{3n}{4}$. $\hfill\blacksquare$
\\

Assume that $C$ is not hamiltonian, as otherwise we are done. Thus $G - V(C)$ has components. Orient $C$ in the clockwise direction and denote the orientation by $\overset{\rightharpoonup }{C}$.

\newcommand{\{x_i^+ \ | \ 1 \leq i \leq 6\}}{\{x_i^+ \ | \ 1 \leq i \leq 6\}}

\begin{CLAIM}\label{hproperties}
Let $H$ be any component of $G - V(C)$. Then we have the following statements:
\begin{enumerate}[(a)]
    \item $|N_C(H)| \geq 2\tau(G) \geq 6$.
    \item for any two $x, y \in N_C(H)$, $xy \not \in E(C)$.
    \item for any two $x, y \in N_C(H)$, $x^+y^+ \not \in E(G)$.
    \item $H$ is a trivial component.
\end{enumerate}
\end{CLAIM}

\textbf{Proof}.\quad
Let $H$ be a component of $G - V(C)$, and $x,y \in N_C(H)$.
Note that since $G$ is 3-tough, $2\tau(G) \geq 6$.

For part $(a)$, assume $|N_C(H)| < 2\tau(G)$.  Then $\frac{|N_C(H)|}{c(G - N_C(H))} < \frac{2\tau(G)}{2} = \tau(G)$, contradicting the toughness of $G$. Thus we have $|N_C(H)| \geq 2\tau(G) \geq 6$.

For part $(b)$, $|N_C(H)| \geq 6$ by part $(a)$. If there exist distinct $x,y \in N_C(H)$ such that $xy \in E(C)$, then let $h_1 \in N_H(x), h_2 \in N_H(y)$.
Assume without loss of generality that $y = x^+$.
As $H$ is connected, there exists some $(h_1, h_2)$-path $P$ in $H$. Then the cycle $C' = xh_1Ph_2y\overset{\rightharpoonup }{C} x$ is a cycle longer than $C$,
contradicting the maximality of $C$.

For part $(c)$, assume for the sake of contradiction that $x^+y^+ \in E(G)$. Let $h_1 \in N_H(x), h_2 \in N_H(y)$.
Assume without loss of generality that $x$ appears before $y$ on the cycle.
Again, as $H$ is connected, there exists some $(h_1, h_2)$-path $P$ in $H$. Then the cycle $C' = x h_1 P h_2 y \overset{\leftharpoonup }{C} x^+ y^+ \overset{\rightharpoonup }{C} x$ is a cycle
longer than $C$, contradicting the maximality of $C$.

For part $(d)$, note that $\{x^+ \ | \ x \in N_C(H)\}$ is an independent set of $G$ by Claim \ref{hproperties} $(c)$. We assume that $H$ is nontrivial, then taking an edge
from $H$ and three vertices from the independent set $\{x^+ \ | \ x \in N_C(H)\}$ gives an induced copy of $(K_2 \cup 3K_1)$.
$\hfill\blacksquare$

\begin{CLAIM}\label{components}
$c(G - V(C)) \leq 3$.
\end{CLAIM}

\textbf{Proof}.\quad
By Claim $\ref{hproperties}$ $(d)$, each component of $G - V(C)$ is trivial. For the sake of contradiction, assume $G - V(C)$ has at least 4 components, and let $x,y,z,w$ be 4 of them.
Then the set $S = V(C) \setminus N_C(\{x,y,z,w\})$ is independent. Otherwise, taking an edge
from the set $S$ and three of $\{x,y,z,w\}$ gives an induced copy of $K_2 \cup 3K_1$.
Then there must exist
3 consecutive vertices $u_1, u_2, u_3 \in V(C)$. Otherwise $|S| \geq \frac{n}{4}$, therefore the independent set $S \cup \{x,y,z,w\}$
has size at least $\frac{3n/4 + 4}{3} > \frac{n}{4}$, contradicting $\tau(G) \geq 3$. Note that none of $\{x,y,z,w\}$
can be adjacent to two consecutive $u_i$, as otherwise we may easily extend the cycle $C$.
Without loss of generality, assume that $u_2 = u_1^+, x \sim u_1, y \sim u_2, z \sim u_3$.
Then we must have $w \sim u_1$, as otherwise taking the edge $xu_1$ and the three vertices
$\{y,z,w\}$ gives an induced copy of $K_2 \cup 3K_1$. Similarly, we must have
$w \sim u_2$. Thus the cycle $C' = u_1 w u_2 \overset{\rightharpoonup }{C} u_1$ is longer than $C$, which is a contradiction.
Therefore $c(G - V(C)) \leq 3$. $\hfill\blacksquare$

\begin{CLAIM}\label{claimwz}
Let $H$ be any component of $G - V(C)$, and $w,z \in N_C(H)$ be any two distinct vertices. Then for any $w_1 \in N_C(w^+)$,
\begin{enumerate}[(a)]
    \item if $w, w_1, z$ appears in the order $w, z, w_1$ along $\overset{\rightharpoonup }{C}$, then
    $z^+ \not\sim w_1^+$.
    \item if $w, w_1, z$ appear in the order $w, w_1, z$ along $\overset{\rightharpoonup }{C}$, then
    $z^+ \not\sim w_1^-$.
\end{enumerate}
\end{CLAIM}

\textbf{Proof}.\quad Let $h_1 \in N_H(w), h_2 \in N_H(z)$. As $H$ is connected, we can find in $H$ an $(h_2, h_1)$-path $P$. For part (a),
assume $z^+ \sim w_1^+$. Then define $$C' = w \overset{\leftharpoonup }{C} w_1^+ z^+ \overset{\rightharpoonup }{C} w_1 w^+ \overset{\rightharpoonup }{C} z h_2 P h_1 w.$$ For part (b), assume $z^+ \sim w_1^-$. Then define $$C' = w \overset{\leftharpoonup }{C} z^+ w_1^- \overset{\leftharpoonup }{C} w^+ w_1 \overset{\rightharpoonup }{C} z h_2 P h_1 w.$$ In either case, $C'$ is a longer cycle than $C$, contradicting the maximality of $C$. $\hfill\blacksquare$

By Claim $\ref{hproperties}$ $(d)$, we let $x \in V(G) \setminus V(C)$ be a component of $G - V(C)$. By Claim
$\ref{hproperties}$ $(a)$, $|N_C(x)| \geq 2\tau(G) \geq 6$. We let $x_1, x_2, \cdots, x_6 \in N_C(x)$ be all distinct vertices
and assume they appear in the order $x_1, \cdots, x_6$ along $\overset{\rightharpoonup }{C}$. Note that $\{x_1^+, \cdots, x_4^+\}$ is an independent set in $G$ by Claim \ref{hproperties} $(c)$.
Assume without loss of generality that $|V(x_4 \overset{\rightharpoonup }{C} x_1)| \geq (|V(C)|-2)/2$. Label the
vertices on this segment by $x_4, x_4^+, y_1, y_2, \cdots, y_t, x_1$. Since $x_4^+ \sim y_1$,
we know that $|N(y_1) \cap \{x_1^+, \cdots x_4^+\}| \geq 2$, as otherwise taking the edge
$x_4^+y_1$ and three vertices from the independent set $\{x_1^+, \cdots, x_4^+\}$ gives an induced copy of
$K_2 \cup 3K_1$. We claim that $|N(y_2) \cap \{x_1^+, \cdots, x_4^+\}| = 0$. Otherwise,
by the same argument as above, we have $|N(y_2) \cap \{x_1^+, \cdots, x_4^+\}| \geq 2$.

If $x_1^+ \in (N(y_1) \cap N(y_2) \cap \{x_1^+, \cdots, x_4^+\})$, then there exists $x_j^+$
with $2 \leq j \leq 4$ such that $x_j^+ \sim y_2$, contradicting Claim $\ref{claimwz}$.

If $x_4^+ \in N(y_1) \cap N(y_2) \cap \{x_1^+, \cdots, x_4^+\})$, then similarly there exists $x_j^+$
with $1 \leq j \leq 3$ such that $x_j^+ \sim y_1$, contradicting Claim $\ref{claimwz}$.

If $(N(y_1) \cap N(y_2) \cap \{x_1^+, \cdots, x_4^+\}) \cap \{x_2^+, x_3^+\} \ne \emptyset$, we may assume that $N(y_1) \cap \{x_1^+, \cdots, x_4^+\} = \{x_4^+, x_j^+\}$
and $N(y_2) \cap \{x_1^+, \cdots, x_4^+\} = \{x_1^+, x_j^+\}$, otherwise we can apply Claim \ref{claimwz}. Then if $x \sim y_1$,
the cycle $C' = x_j x y_1 \overset{\leftharpoonup }{C} x_j^+ y_2 \overset{\rightharpoonup }{C} x_j$is longer than $C$. Thus $x \not \sim y_1$. Then taking the independent set
$\{x, x_1^+, x_2^+, x_3^+\} \setminus \{x_j^+\}$ and the edge $y_1x_j^+$ gives an induced copy of $K_2 \cup 3K_1$, which is a
contradiction.

Thus $N(y_1) \cap N(y_2) \cap \{x_1^+, \cdots, x_4^+\} = \emptyset$. We may assume that
$N(y_1) \cap \{x_1^+, \cdots, x_4^+\} = \{x_3^+, x_4^+\}$ and $N(y_2) \cap
\{x_1^+, \cdots, x_4^+\} = \{x_1^+, x_2^+\}$, else we can apply Claim \ref{claimwz} as in previous arguments.
Note that if $x \sim y_1$, the cycle $C' = x_1 x y_1 \overset{\leftharpoonup }{C} x_1^+ y_2 \overset{\rightharpoonup }{C} x_1$ is longer than $C$. Thus
$x \not \sim y_1$.
Then taking the set $\{x, x_1^+, x_2^+\}$ and the edge $y_1x_3^+$ gives an induced copy
of $K_2 \cup 3K_1$, which is a contradiction. Thus $N(y_2) \cap \{x_1^+, x_2^+, \cdots, x_4^+\} = \emptyset$.

We claim that $|N(y_3) \cap \{x_1^+, \cdots, x_4^+\}| \geq 2$. Otherwise, taking 3 vertices from the independent set
$\{x_1^+, \cdots, x_4^+\}$ and the edge
$y_2y_3$ gives an induced copy of $K_2 \cup 3K_1$. Recall that the last vertex before $x_1$ on the segment
$x_4 \overset{\rightharpoonup }{C} x_1$ is labelled by $y_t$. Define $Y$ as the set of all even-indexed $y_i$. Then $Y = \{y_2, y_4, \cdots, y_t\}$ in the
case that $t$ is even and $Y = \{y_2, y_4, \cdots, y_{t-1}\}$ in the case that $t$ is odd. By applying a similar
argument as in the $y_2$ case to all even-indexed $y_i$, we see that $N(Y) \cap \{x_1^+, \cdots, x_4^+\} = \emptyset$.

Note that $Y$ is an independent set of $G$, as otherwise
taking an edge with end-vertices in $Y$ and 3 vertices from the set $\{x_1^+, \cdots, x_4^+\}$ gives
an induced copy of $K_2 \cup 3K_1$. Thus, $Y \cup \{x_1^+, \cdots, x_4^+\}$ is also an independent set in $G$.

By Claim $\ref{hproperties}$ $(d)$ and Claim $\ref{components}$, we have $|V(C)| \geq n - 3$.
Thus $|Y| \geq \floor{\frac{1}{2}(\frac{|V(C)|-2}{2})} \geq \floor{\frac{n-5}{4}}$.
Therefore $|Y \cup \{x_1^+, \cdots, x_4^+\}| > \frac{n}{4}$, contradicting $\tau(G) \geq 3$. This completes the proof.
\qed

\section*{Acknowledgements}
The authors wish to thank Dr. Songling Shan for her guidance during the research process and her valuable
feedback.

\bibliographystyle{plain}