\section{Introduction}
The discovery of ferromagnetism in two-dimensional (2D) monolayer CrI\textsubscript{3} and other 2D van der Waals (vdW) materials opened a new window for exploring low dimensional magnetism and its applications in spintronics\cite{gibertini2019magnetic, gong2019two, burch2018magnetism, huang2017layer, gong2017discovery, deng2018gate, bonilla2018strong, o2018room, soriano2019interplay}. The properties of 2D materials can be controlled by external parameters\cite{cristoloveanu2019concept,ye2012superconducting,cao2018unconventional,huang2018electrical,song2019switching} and are highly sensitive to stacking and twisting between the layers\cite{ponomarenko2013cloning,dean2013hofstadter,gorbachev2014detecting,cao2018correlated,cao2018unconventional}. In particular, with the discovery of superconductivity in twisted bilayer graphene\cite{cao2018correlated,cao2018unconventional}, there has been tremendous progress on exploring moir\'e superlattices both experimentally and theoretically~\cite{zhang2017interlayer,moon2014electronic,chen2019evidence,chen2019signatures,wallbank2015moire,ni2015plasmons}. In terms of magnetism, stacking order and twisting can significantly alter the interlayer exchange as the exchange is highly sensitive to atomic registries\cite{Chen_Science2019, Sivadas_NanoLett2018, Jiang_PRB2019, du2016weak,klein2018probing, tong2018skyrmions}.

Magnetic skyrmions~\cite{roessler2006spontaneous} are nanoscale  vortex-like spin textures that were first observed in non-cetrosymmetic bulk magnetic materials such as MnSi~\cite{muhlbauer2009skyrmion,yu2010real}, (FeCo)Si~\cite{munzer2010skyrmion} and FeGe~\cite{yu2011near}. Skyrmions are topological defects, and they exhibit novel transport phenomena such as topological Hall effect and topological Nerst effect~\cite{hamamoto2015quantized, shiomi2013topological}. In recent years skyrmions received ample attention due to their potential for  spintronics applications and memory storage devices~\cite{fert2013skyrmions}. In most cases, skyrmions are stabilized by interplay of Dzyaloshinskii-Moriya (DM) interaction and external magnetic field\cite{muhlbauer2009skyrmion,yu2010real,banerjee2014enhanced}. In this article, we explore the possibility of stabilizing magnetic skyrmions in the absence of an external magnetic field in moir\'e superlattices. We consider a ferromagnetic (FM) monolayer on an antiferromagnetic (AFM) substrate with N\'eel order. Twisting the FM layer by an angle $\theta$ produces moir\'e patterns as shown in Fig. 1(a), (c). Ferromagnetic coupling between the substrate and the FM monolayer leads to an alternating exchange field for the moir\'e superlattice as shown in Fig. 1(b). Our setup is motivated from Ref. \onlinecite{tong2018skyrmions}. However unlike Ref. \onlinecite{tong2018skyrmions} which includes dipole-dipole interaction to stabilize magnetic skyrmions, we consider DM interaction which is the primary interaction for magnetic skyrmions in chiral magnets\cite{muhlbauer2009skyrmion, yu2010real, johnson1996magnetic, dubowik1996shape, Banerjee_NatPhys2013}.
\begin{figure}[t]
\includegraphics[width=\columnwidth]{Fig1.png}
\caption{(a) Schematic of our setup of a FM monolayer twisted on top of an AFM substrate with N\'eel order. We consider honeycomb lattice for both the monolayer and the substrate. (b) Exchange field due to interlayer coupling. (c) Moir\'e pattern due to small twisting angle between FM monolayer (green) and top most layer of AFM substrate (red and black). Three local atomic registries $\textbf{R}_0$, $\textbf{R}_1$ and $\textbf{R}_2$ are zoomed which matches the atomic registries of different interlayer translations.}
\label{Layer+substrate}
\end{figure}

Our main results are summarized in Fig.~\ref{L_vs_B_{Max}} and~\ref{A_vs_B_{Max}}. We show that (i) skyrmion crystal (SkX) is stabilized as a function of exchange coupling between the layers ($J_{\rm ex}$) and moir\'e periodicity. (ii) Even though SkX can be stabilized for a wide range of twisting angle, we find the optimal moir\'e periodicity to be about $L=9L_D$, where $L_D=(J/D)a$ is the intrinsic length scale for skyrmions. (iii) Unlike chiral magnets under magnetic field, we find an extended region of SkX for easy-axis anisotropy. (iv) We show that a large fraction of the topological charge $q=\frac{1}{4\pi}\int d^2r \widehat{\textbf{m}}.(\partial_{x}\widehat{\textbf{m}}\times\partial_{y}\widehat{\textbf{m}})$ of the magnetic skyrmions is concentrated at the edges and splits into three parts for large moir\'e periodicity and large easy axis anisotropy. This effect arises due to the anisotropic shape of the skyrmion.\\
\begin{figure*}[t]
\includegraphics[width=\textwidth]{Fig2.png}
\caption{$A=0$ phase diagram at different moir\'e periods and exchange field. (a) Moir\'e period $L$ vs  maximum interlayer exchange field $B_{Max}$ phase diagram with spiral (Sp) and skyrmion crystal (SkX) phases. Here $L_D=(J/D)a$, where $a$ is the lattice constant and the purple diamonds correspond to error bars. The dashed lines at $L=5L_D$ and $L=8L_D$ represent the cross section of Fig.~\ref{A_vs_B_{Max}}(a) and Fig.~\ref{A_vs_B(L=5L_D)} respectively. (b,c,d) Magnetization texture and (e,f,g) topological charge density in $1\times 2$ moir\'e super cells for $L/L_D = (5, 8, 11)$ and $B_{Max}J/D^2=1.73$ as marked by $\bigstar, \blacktriangle, \blacksquare$ symbols in (a). The colors represent out of plane component ($m_z$) of magnetization and topological charge density ($\chi$) in (b,c,d) and (e,f,g) respectively. The arrows show the in-plane component of magnetization in (b,c,d) and the dotted lines represent $1\times 2$ moir\'e super cells in (b,c,d,e,f,g).
}
\label{L_vs_B_{Max}}
\end{figure*}

\section{Model}
Before we delve into the analysis of the effective magnetic Hamiltonian, we first describe our setup. As mentioned above, we follow the procedure of Ref. \onlinecite{tong2018skyrmions} to derive the effective interlayer exchange field. We consider a FM monolayer twisted on top of an AFM substrate,  both on a honeycomb lattice with the same lattice constant, $a$. For twisting angle $\theta$, the moir\'e period is given by $L= a/2\sin(\theta/2)$. For small angle $\theta$ and/or lattice mismatch $\delta$, $L\thickapprox a/\sqrt{\theta^2+\delta^2}$ (large period), the local atomic registries on length scale smaller than $L$ but larger than $a$ matches the atomic stacking of different interlayer translation $\textbf{r}$ as shown in Fig.~\ref{Layer+substrate}(c). Hence moir\'e superlattice can be described by interlayer translation vector $\textbf{r(R)}$ that gives atomic registry at position $\textbf{R}$. The interlayer exchange coupling between AFM substrate and FM layer is different at different positions due to different atomic stacking of monolayer and substrate and this leads to spatially dependent exchange field $\textbf{B}\textbf{(R)}$ as shown in Fig.~\ref{Layer+substrate}(b). For example, at position $\textbf{R}_1$ the coupling aligns the spins of 2D layer (green) in positive z-direction when it sits on top of AFM sublattice with spins up (black) and the spins align in opposite direction  at $\textbf{R}_2$ when it is on top of AFM sublattice with spins down (red). The interlayer exchange field at interlayer translation $\textbf{r}$ is give by~\cite{tong2018skyrmions}
\begin{figure*}[t]
\includegraphics[width=\textwidth]{Fig3.png}
\caption{$L=8L_D$ phase diagram. (a) Anisotropy $A$ versus maximum interlayer exchange field $B_{Max}$ phase diagram with ferromagnetic (FM), spiral (Sp), skyrmion crystal (SkX) and mixed state (Sp+FM). Here $L_D=(J/D)a$, where $a$ is the lattice constant. (b,c,d) Magnetization texture in $1\times 2$ moir\'e supercells and (e,f,g) spin structure factor for ($\bigstar, \blacktriangle, \blacksquare$) symbols marked in (a). The parameters corresponding to these symbols are the following: $\bigstar$ = $\{ A=-0.6D^2/J$, $B_{Max}=1.15D^2/J \}$, $ \blacktriangle$ = $\{ A=0.3D^2/J$, $B_{Max}=0.38D^2/J \}$ and $\blacksquare$ = $\{ A=1.14D^2/J$, $B_{Max}=0.58D^2/J \}$.
}
\label{A_vs_B_{Max}}
\end{figure*}
\begin{eqnarray}
\textbf{B}^{\tau}(\textbf{r})=\sum_{j,\tau^\prime}J^{\tau^\prime}_{ex}(\textbf{r}-\bm{\xi}_{\tau}+\textbf{R}_j)\textbf{m}_{j,\tau^\prime}
\end{eqnarray}
where $J_{ex}$ is the interlayer coupling coefficient, $\textbf{m}$ is the magnetic moment of top most layer of AFM substrate,  $\left\lbrace\tau,\tau^\prime\right\rbrace=\left\lbrace A,B \right\rbrace$ represents the two in-equivalent sites in unit cell, $\bm{\xi}_A=0$, $\bm{\xi}_B=\left\lbrace 0,a \right\rbrace$ and the summation $j$ is over the Bravais lattice. The total interlayer exchange field per unit cell is given by summing the fields of site $A$ and $B$
\begin{eqnarray}
\textbf{B}(\textbf{r})=\textbf{B}^{A}(\textbf{r})+\textbf{B}^{B}(\textbf{r}).
\end{eqnarray}
This approximation holds when the interlayer coupling, $J_{ex}$ is small as compared to intralayer coupling, $J$. In general,  $J_{ex}< J$ holds for vdW magnets as the interplane exchange is expected to be much smaller than intraplane exchange.
We used the following coupling form that decays exponentially at long distances
\begin{eqnarray}
J^{\tau^\prime}_{ex}(\textbf{r}-\bm{\xi}_{\tau}+\textbf{R}_j)=J^{0}_{ex}e^{-\frac{\sqrt{(\textbf{r}-\bm{\xi}_{\tau}+\bm{\xi}_{\tau^\prime}+\textbf{R}_j)^2+d^2}}{r_0}}
\end{eqnarray}
where $d$ is the interlayer separation and $r_{0}$ is the decay length. In our calculations we used $d=a$, $r_0=a$ and this leads to
\begin{eqnarray}
B_{Max}\thickapprox 1.9244\times 10^{-2}J^{0}_{ex}.
\end{eqnarray}

Next we describe our model for the monolayer. We consider a magnetic model for 2D honeycomb lattice which is relevant to 2D vdW magnets such as trihallides~\cite{huang2017layer}

\begin{eqnarray}
H&=& -J\sum_{{r,\mu}} {\bf S}_{r}. ( {\bf S}_{r+\hat{\delta}_{\mu}} )
-D \sum_{{r,\mu}}[
\hat{d}_{\mu}. ({\bf S}_{r} \times {\bf S}_{r+\hat{\delta}_{\mu}} )]
\nonumber\\
&& -A_{c} \sum_{{r,\mu}}[
({\bf S}_{r}.  \hat{d}_{\mu} ) ({\bf S}_{r+\hat{\delta}_{\mu}}.  \hat{d}_{\mu} )
]+A_{s} \sum_{{r}}(S_{r}^z)^2\nonumber\\
&&-\sum_r {\bf B}({\bf R}_r)\cdot {\bf S}_r
\label{Eq:Gnrl_Hmlt_1}
\end{eqnarray}
where $\vec{S}_{r}$ is local moment at site $r$ and $\hat{\delta}_{\mu}$ are the three nearest neighbors on the honeycomb lattice.
$J$ is the ferromagnetic Heisenberg exchange coupling, $D$ is the DM coupling~\cite{Moriya_PR1960}. $A_{c}$ and $A_s $ are the compass and single-ion anisotropies respectively. The DM vector $\hat{d}_{i}=\hat{z}\times\hat{\delta}_{i}$ is set by the symmetry and originates due to the inversion symmetry breaking on the surface. However, the microscopic mechanism for the DM interaction in insulating vdW magnets is different than metallic multilayers since the former is due to the superexchange interaction in the presence of an electric field as originally discussed by Moriya\cite{Moriya_PR1960} whereas the latter can be ascribed to asymmetric interaction paths proposed by Fert and Levy\cite{Levy_PRB1981, Fert_PRL1980}.
$\bf B(R)$ is the interlayer exchange field with the twisted substrate.
To explore the phase diagram of $H$, we consider the free energy functional in the continuum $\mathcal{F}[{\bf m}] = \int d^2r F({\bf m})$ where ${\bf m}({\bf r})$ is the local magnetization.  We set the lattice constant, $a=1$. The free energy density has the following four components
\begin{eqnarray}
F({\bf m}) =F_{iso}+F_{DM}+F_{aniso}+F_{moire}
\end{eqnarray}
where
\begin{eqnarray}
F_{iso}&=&F_{0}(\textbf{m})+\frac{3}{2}(J/2)\sum_{{\alpha}}(\nabla m^{\alpha})^{2}
\end{eqnarray}
\begin{eqnarray}
F_{DM}&=&-\frac{3}{2}D(m^{z}\partial_{x}m^{x}-m^{x}\partial_{x}m^{z})\nonumber\\
&& +\frac{3}{2}D(m^{y}\partial_{y}m^{z}-m^{z}\partial_{y}m^{y})
\end{eqnarray}
\begin{eqnarray}
F_{aniso}&=&-\frac{3}{2} A_{c} [(m^{x})^{2}+(m^{y})^{2}]+A_{s}(m^{z})^{2}\nonumber\\
&&+\frac{3A_{c}}{4} [ m^{y} (\partial_{x}m^{y}+\partial_{y}m^{x})\nonumber\\
&&-m^{x} (\partial_{x}m^{x}-\partial_{y}m^{y}) ]\\
F_{moire}&=&-{\bf B(R)} \cdot {\bf m}.
\end{eqnarray}
We absorb factor $\frac{3}{2}$ in $J$, $D$ and $A_{c}$ and define the effective anisotropy $A=A_{c}+A_{s}$ which can be positive (easy plane) or negative (easy axis or Ising). The last two terms in $F_{aniso}$, given as $\sim[ m^{y} (\partial_{x}m^{y}+\partial_{y}m^{x}) -m^{x} (\partial_{x}m^{x}-\partial_{y}m^{y}) ]$ have no contribution to the free energy for systems with periodic boundary conditions since they are total derivatives. To obtain the ground state spin configuration $\textbf{m}$, we solve the coupled Landau-Lifshitz-Gilbert (LLG) equations~\cite{1353448}
\begin{eqnarray}
\frac{d\textbf{m}}{dt}=-\gamma \textbf{m}\times \textbf{B}^{eff}+\alpha \textbf{m} \times \frac{d\textbf{m}}{dt},
\end{eqnarray}
where  $\textbf{B}^{eff}=-\delta H /\delta \textbf{m}$, $\gamma$ is gyromagnetic ratio and $\alpha$ is Gilbert damping coefficient. We start from different initial states\footnote{We consider multiple random configurations as well as ferromagnetic states pointing along (0,0,1), (1,1,5), (1,1,4), (1,1,3), (1,1,2), (1,1,1), (1,1,0), (2,2,1), (3,3,1), (4,4,1), (-1,-1,0), (-1,-1,-1), (-1,-1,-2), (-1,-1,-3), (0,0,-1) directions.} and compare the energies of final states to get the actual ground state. To solve LLG equations  numerically we used mid point method~\cite{d2005numerical} by discretizing the effective magnetic field $\textbf{B}^{eff}=-\delta H /\delta \textbf{m}$ on a $1\times2$ moir\'e supercells. We used $L_D=(J/D)a=10a$ to construct the phase diagrams and the results were also verified at various points for larger values of $L_D$. The magnitude of the magnetization was kept constant at each grid point after each time step enforcing the hard spin constraint, $|{\bf m}|^2=1$ and periodic boundary conditions were imposed at the boundaries. There is unavoidable error due to solving the LLG equations on a discrete lattice. The error is small as long as the variation of the magnetization is small within neighboring sites. We checked that our results are robust as a function of increasing system.

Our method can only capture magnetic orders that are commensurate with the 1x2 moir\'e supercell. However, in the limit of $J_{ex} \rightarrow 0$, the ground state is a spiral with the wave vector $q \simeq D/J$ for weak anisotropy\cite{banerjee2014enhanced} which is in general not commensurate with the moir\'e supercell. Therefore, our method is not suitable to explore the weak $J_{ex}$ limit. Such incommensurate phases and incommensurate to commensurate phase transitions have recently been analyzed in Ref. \citenum{Hejazi_arXiv2020}. The parameters used in this article lie outside the regime of incommensurate phases\cite{Hejazi_arXiv2020}, yet it is not possible to extrapolate our results to $J_{ex} \rightarrow 0$ limit due to this reason.

\section{Results}
We start by exploring the interplay between the moir\'e periodicity $L$ and $B_{Max}$ as shown in Fig.~\ref{L_vs_B_{Max}} for $A=0$. $B_{Max}$ is the maximum value of interlayer exchange magnetic field. As shown in Fig 1(b), the moir\'e supercell splits into two triangles of opposite alternating effective magnetic field which is maximum at the centers and vanishingly small at the corners of the triangles. Therefore, in the limit of large $B_{Max}$, the magnetization aligns with the effective magnetic field at the center of the triangles, creating two ferromagnetic domains separated with a domain wall whose chirality is set by D. The only degree of freedom left that determines the ground state is the magnetization at the corners of the triangles which creates an effective triangular lattice. In particular, if the magnetization at the corners are ferromagnetic, the state is a skyrmion whereas a stripe order give rise to a spiral phase. We find that the ground state is a spiral for low exchange field. As we increase the field, the SkX phase starts at the moir\'e period $L/L_{D}\approx9$ which corresponds to the optimum angle between the layer and the  substrate. As we further increase the field, we get SkX for a range of values around optimum period. This range increases with increasing the exchange field. Unlike skyrmions in chiral magnets where the size of the skyrmion is set by $L_D$, here we find that the size of skyrmions is determined by the moir\'e period. On the other hand, $L_D$ determines the boundary length between the interior and the exterior of the skyrmions. For small moir\'e period, skyrmions are small and their shape is nearly circular as shown in Fig.~\ref{L_vs_B_{Max}}(b). Fig~.\ref{L_vs_B_{Max}}(c) shows that as the period increases, the size of skyrmion also increases and it takes the triangular shape of the exchange field. The corners of skyrmions get sharper with increasing $L$. Unlike the skyrmions in chiral magnets, we find that a large fraction of the topological charge is concentrated at the edges of the skyrmion. This fraction increases with increasing $L$. There is also a small fraction of opposite charge between the skyrmions which decreases with increasing $L$. This charge arises due to the anti-vortices between the skyrmions\cite{lin2015skyrmion}. For large $L$, the topological charge further splits into three parts due to the triangular shape of the skyrmion as shown in Fig.~\ref{L_vs_B_{Max}}(g).\\

\begin{figure}
\includegraphics[width=\columnwidth]{Fig4}
\caption{Evolution of magnetization texture, topological charge density $\chi$ and spin structure factor $I(Q)$ for $AJ/D^2=\{-1.2, -0.6, 0\}$ at $L=8L_D$.}
\label{SkX_vs_A}
\end{figure}
Next, we explore the effects of anisotropy, $A$ at around optimal angle $L=8L_{D}$ as shown in Fig.~\ref{A_vs_B_{Max}}. At low exchange field, we obtain spiral phase for a wide range of $A$, a small ferromagnetic phase at lowest negative values of $A$ as well as a mixed (FM+Sp) phase at largest positive values of $A$. As the field increases, initially we get SkX near the lowest negative values of $A$ that corresponds to easy axis anisotropy. By further increasing the field, the range of SkX gradually increases and eventually, occupies the whole phase diagram. Fig.~\ref{A_vs_B_{Max}}(b,c,d) show the local magnetization and Fig.~\ref{A_vs_B_{Max}}(e,f,g) show the spin structure factor $I(Q)\propto|\langle\textbf{m}_{Q}\rangle|^2$ for the three phases. $\textbf{m}_Q$ is the Fourier transform of the
magnetization. Unlike an isotropic SkX which has six peaks on the circle in the spin structure factor, we find four peaks lie on circle and
two peaks lie inside the circle. This is due to the anisotropic triangular shape of the SkX. The spiral has two peaks at $Q = \pm Q_0$ and the mixed state (FM+Sp) has three peaks including the $Q=0$ from the FM and $Q=\pm Q_0$ from the spiral phase. The intensity of $Q=0$
peak increases with increasing $A$ and decreases with increasing exchange field $B_{Max}$. On the other hand, the intensities corresponding to the spiral wave vectors have the opposite behavior of $Q=0$ with $A$ and $B_{Max}$. \\
The properties of SkX also depends on the anisotropy $A$. Fig.~\ref{SkX_vs_A} shows the magnetization, spin structure factor and topological charge density as a function of $A$ for moir\'e period $L=8L_D$. For $A=-1.2D^2/J$ the skyrmion has a sharp boundary wall
where magnetization changes abruptly and then it changes slowly inside the skyrmion. The sharpness of boundary wall decreases with
increasing $A$ and the change in magnetization inside the skyrmion increases with increasing $A$. For $A=-1.2D^2/J$, a large fraction of topological charge is concentrated at the boundary wall of skyrmion and a small fraction lies inside the skyrmion. There is also a small fraction of opposite charge between the skyrmions. The concentration of charge at boundary decreases with increasing $A$ and the central charge increases with increasing $A$. The fraction of opposite charge between the skyrmions also increases with increasing $A$\cite{lin2015skyrmion}.

We also studied the effects of anisotropy for a non-optimal angle at $L = 5L_D$. As shown in Fig.~\ref{A_vs_B(L=5L_D)}, SkX is highly suppressed in this case but still persists for large exchange field and easy axis anisotropy. Suppression of SkX is due to the fact that the moir\'e supercell is too small with respect to the optimal size of the skyrmions.

\begin{figure}[h]
\includegraphics[width=\columnwidth]{Fig5.png}
\caption{Anisotropy $A$ vs  maximum interlayer exchange field $B_{Max}$ phase diagram with ferromagnetic (FM), spiral (Sp), skyrmion crystal (SkX) and mixed state (Sp+FM) phases at $L=5L_D$. Here $L_D=(J/D)a$, where $a$ is the lattice constant.}
\label{A_vs_B(L=5L_D)}
\end{figure}

Our results apply to a wide range of vdW magnets since the primary constraint in our model is the requirement that $J_{ex}<J$ which is in general satisfied for vdW magnets. Our phase diagrams span a wide range of parameters and for a fixed twisting angle, our results scale with $AJ/D^2$.

\section{Conclusion}
We have shown that skyrmion crystal can be stabilized in moir\'e superlattices in the absence of external magnetic field. We found a large SkX phase for easy axis anisotropy which can be essential to stabilize skyrmions in vdW magnets such as CrI$_3$\cite{huang2017layer}. In particular, for optimal moir\'e periodicity SkX occupies the majority of the phase diagram. We find that the properties of the skyrmion can be tuned with the moir\'e periodicity and anisotropy. Unlike skyrmions in chiral magnets, the topological charge density depends on the size of the skyrmions sand it is concentrated at the edges for skyrmions with large L.

\section{Acknowledgements}
We thank Sumilan Banerjee for useful discussions. MA is supported by Fulbright scholarship and ASU startup grant and OE is supported by NSF-DMR-1904716. We acknowledge the ASU Research Computing Center for HPC resources.


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

Ordinary differential equations (ODE) models are widely used in the sciences \citep[e.g., see][]{Strogatz2014, Murdoch2003, BeuterGlassMackeyTitcombe2003,Diekmann2000, AndersonMay1992}, but are often limited by the difficulty of incorporating time delays. Such delays are more easily incorporated into models using other mathematical frameworks. For example, delays can be modeled using integral equations or integro-differential equations to incorporate distributed delays into dynamic models, or using delay differential equations (DDEs) to incorporate fixed delays \citep{Wearing2005,Feng2016,Feng2007,Ruan2006,Lin2018,Roussel1996}. In an ODE framework, the linear chain trick has long been used to incorporate exponential and Erlang distributed delays \citep{Metz1991,Metz1986,Nisbet1989,Smith2010,Diekmann2017}, and recently this approach has been generalized to a much broader family of delay or dwell time distributions \citep{Hurtado2019}.

This generalized linear chain trick \citep[GLCT;][]{Hurtado2019} allows modelers to incorporate delay and dwell time distributions that include  (but are not limited to) the phase-type family of univariate probability distributions \citep{Bladt2017ch3,Horvath2016,Reinecke2012a,BuTools2}, which include hyperexponential, hypoexponential, Coxian, and the previously mentioned exponential and Erlang distributions. The phase-type family of distributions is the set of all possible absorption time distributions for continuous time Markov chains (CTMCs) with one or more transient states and a single absorbing state. The GLCT also permits the use of similar time-varying versions of such distributions \citep{Hurtado2019}. Together, the GLCT and the tools and techniques associated with phase-type distributions enable modelers to draw from a richer set of ODE model assumptions when constructing new models, and provide a framework for more clearly seeing how underlying stochastic model assumptions are reflected in the structure of mean field ODE models.

In this paper, we illustrate how to use the GLCT to formulate new ODE models that explicitly incorporate phase-type distributed delays. We do this by generalizing multiple different models from the literature. These include ODE models with no explicit delays, DDE models, and distributed delay models in the form of integro-differential equations. In section \ref{sec:glct} we review the GLCT framework for phase-type distributions, as well as the standard Linear Chain Trick (LCT). We then generalize multiple models starting in section \ref{sec:TGI}, where we generalize a tumor growth inhibition model by \citet{Simeoni2004}. In section \ref{sec:opioid} we then generalize a prescription opioid epidemic model by \citet{Battista2019}, then a within-host immune-pathogen model by \citet{Hurtado2012} in section \ref{sec:inhost}, and finally a model of cell-to-cell spread of HIV by \citet{Culshaw2003} in section \ref{sec:HIV}.

\subsection{Generalized Linear Chain Trick}\label{sec:glct}
For our purposes below, we here provide a statement of the Generalized Linear Chain Trick (GLCT) for phase-type distributions. More generally, the GLCT in \citet{Hurtado2019} extends the version below to also include time-varying parameters (analogous to extending homogeneous Poisson process rates to the time-varying rates of inhomogeneous Poisson processes).

The (continuous) phase-type distributions are a family of matrix exponential distributions that can otherwise be thought of as the absorption time distributions for CTMCs with a single absorbing state. They are parameterized in terms of a vector \textbf{v}, which is the initial distribution vector across the set of transient states\footnote{The full initial distribution vector would be [\textbf{v}, $v_0$].}, and the transient state block \textbf{M} of the transition rate matrix, which together define the corresponding CTMC. The general form for the density function, cumulative distribution function, and moments of a phase-type distribution are  \begin{subequations} \begin{align}
	f(t) =&\; \text{\textbf{v}}\,e^{\mathbf{M}t}\,(-\mathbf{M}\mathbf{1}) \\
	F(t) =&\; 1 - \text{\textbf{v}}\,e^{\mathbf{M}t}\,\mathbf{1} \\
	E(T^j)=&\; j!\,\text{\textbf{v}}\,(-\mathbf{M})^{-j}\mathbf{1}
\end{align} \end{subequations} where $\mathbf{1}$ is an appropriately long column vector of ones.

For more on phase-type distributions, see \citet{HurtadoRichards2020, Hurtado2019, Bladt2017ch3, BuTools2, Horvath2016, Reinecke2012a}.

\begin{theorem}[\textbf{GLCT for phase-type distributions}] \label{th:glct}

	Assume individuals enter a state (call it state X) at rate $\mathcal{I}(t)\in\mathbb{R}$ and that the distribution of time spent in state X follows a continuous phase-type distribution given by the length $k$ initial probability vector $\textnormal{\textbf{v}}$ and the $k\times k$ matrix $\mathbf{M}$. Then partitioning X into $k$ sub-states X$_i$, and denoting the corresponding amount of individuals in state X$_i$ at time $t$ by $x_i(t)$, then the mean field equations for these sub-states $x_i$ are given by
	\begin{equation} \frac{d}{dt}\mathbf{x}(t)=\textnormal{\textbf{v}}\,\mathcal{I}(t) + \mathbf{M}^\text{T}\,\mathbf{x}(t) \label{eq:GLCT}\end{equation}
	where the rate of individuals leaving each of these sub-states of X is given by the vector $(-\mathbf{M\,1})\circ\mathbf{x}$, where $\circ$ is the Hadamard (element-wise) product of the two vectors, and thus the total rate of individuals leaving state X is given by the sum of those terms, i.e., $(-\mathbf{M\,1})^\text{T}\mathbf{x}=-\mathbf{1}^\text{T}\mathbf{M}^\text{T}\mathbf{x}$.
\end{theorem}

The standard Linear Chain Trick (LCT) is well known (see \citet{Hurtado2019} and references therein) and is a special case of Theorem \ref{th:glct} above. However, it is usually stated without the above matrix-vector notation. The following is a formal statement of the standard LCT, but here we have slightly generalized it to include generalized Erlang distributions (i.e., the sum of $k$ independent exponentially distributed random variables, each with potentially different rates $r_i$) as this only changes a few subscripts in the mean field equations. See \citet{Smith2010} for a similar statement of the standard LCT (for Erlang distributions).

\begin{corollary}[\textbf{Linear Chain Trick for Erlang and Hypoexponential Distributions}] \label{th:lct}
	~\\
	Consider the GLCT above. Assume that the dwell-time distribution is a generalized Erlang (hypoexponential) distribution with rates $\mathbf{r}=[r_1,\;r_2,\;\ldots\;,r_k]^\text{T}$, where $r_i>0$, or an Erlang distribution with rate $r$ (all rates $r_i=r$) and shape $k$ (or if written in terms of shape $k$ and mean $\tau=k/r$, use $r=k/\tau$). Then the corresponding mean field equations are
		\begin{equation} \label{eq:LCT}\begin{split}
		\frac{dx_1}{dt} =&\; \mathcal{I}(t) - r_1\,x_1 \\
		\frac{dx_2}{dt} =&\;  r_1\,x_1 - r_2\,x_2 \\
		&\vdots \\
		\frac{dx_k}{dt} =&\; r_{k-1}\,x_{k-1} - r_{k}\,x_{k}.
		\end{split}\end{equation}
\end{corollary}
\begin{proof}
	The phase-type distribution formulation of the generalized Erlang distribution described above is given by \begin{equation} \label{eq:vM} \textnormal{\textbf{v}} = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0  \end{bmatrix}  \qquad \text{and} \qquad \mathbf{M} = \begin{bmatrix} -r_1 & r_1 & 0 & \cdots & 0 \\
	0 & -r_2 & r_2  & \ddots & 0 \\
	\vdots & \ddots & \ddots  & \ddots & \ddots \\
	0 & 0 & \ddots &  -r_{k-1} & r_{k-1} \\
	0 & 0 & \cdots &  0 & -r_k
	\end{bmatrix}. \end{equation} Substituting these into eq. \eqref{eq:GLCT} yields the desired result. If all $r_i=r$ then the phase-type distribution is an Erlang distribution with rate $r$ and shape $k$ (mean $k/r$ and coefficient of variation $1/\sqrt{k}$).
\end{proof}

\subsubsection{Using The GLCT To Derive New ODEs With Phase-Type Distributed Delays}

As we will show below, it is relatively straightforward to derive, from an existing ODE model with exponential or Erlang delays, a more general ODE model with phase-type distributed delays. Such models can also be similarly derived from existing integral equations \citep{Hurtado2019}. In this section, we describe how modelers can also use the GLCT to derive ODEs with phase-type distributed delays from existing ODEs and DDEs.

One important application of the LCT is that it can be used to approximate delay differential equations (DDEs) with ODEs, as discussed in \citet{Smith2010}. To do this, DDEs can be thought of in the context of distributed delay models as the result of assuming a delay distribution with a point mass at $\tau$. This distribution can be approximated by an Erlang distribution with mean $\tau$ and a very small coefficient of variation (i.e., a large shape parameter), which yields an ODE via the LCT. That Erlang distribution assumption can then be replaced with a more general phase-type distribution assumption to yield ODEs via the GLCT, as illustrated in the following example.

\textbf{Example:} Consider a simple birth-death model where the recruitment rate ($f$) of new adults at time $t$ is a function of the number of adults $\tau$ time units in the past. Suppose $x$ tracks the number of adults in a population, that there is an assumed maturation time of $\tau$ time units, and that adults die with per-capita rate $g(x)$. Then the dynamics of $x$ could be modeled by the DDE

\begin{equation}
\frac{dx}{dt} =\; f(x(t-\tau)) - g(x(t))\,x(t).
\end{equation}

The LCT can be used to approximate the above DDE with an ODE, by replacing the assumption of a fixed time delay (or maturation time, in this example) with an Erlang distributed delay, which can have arbitrarily small variance. This is accomplished by using a shape parameter $k\gg1$ (recall the coefficient of variation is given by $1/\sqrt{k}$), and assuming a rate $r=k/\tau$ which yields the desired Erlang distribution with mean $r/k=\tau$. By the LCT, this introduces $k$ new state variables to the model, which track this delayed quantity, and in this context these new states can be thought of as a sequence of immature stages. Applying the LCT in this manner yields the following ODE approximation of the above DDE, where these immature stages are denoted $w_i$, $i=1,\ldots,k$. \begin{subequations} \label{eq:DDE2ODE}
\begin{align}
\frac{dw_1}{dt} =&\; f(x) - r\,w_1 \\
\frac{dw_j}{dt} =&\; r\,w_{j-1}  - r\,w_j, \quad\; j=2,\ldots,k \\
\frac{dx}{dt} =&\;  r\,w_k - g(x)\,x
\end{align}\end{subequations}

This approximation can then be used as an intermediate step to derive a phase-type distributed delay model. To do this, simply write the above system of ODEs in matrix form \textit{a la} the GLCT,  \begin{subequations} \label{eq:DDE2PTODE}
\begin{align}
\frac{d\mathbf{w}}{dt} =&\; f(x)\text{\textbf{v}} + \mathbf{M}^\text{T} \mathbf{w} \\
\frac{dx}{dt} =&\;  -\mathbf{1}^\text{T}\mathbf{M}^\text{T}\mathbf{w} - g(x)\,x,
\end{align}\end{subequations}

where the term $r\,w_k$ is instead written using the more generic expression given in Theorem \ref{th:glct} for computing the overall loss rate from the intermediate states $\mathbf{w}$, which in this case is $-\mathbf{1}^\text{T}\mathbf{M}^\text{T}\mathbf{w} = r\,w_k$, and the vector \textbf{v} and matrix \textbf{M} are of the form of eqs. \eqref{eq:vM}, with $r_i=k/\tau$

The above steps leading to the derivation of eqs. \eqref{eq:DDE2PTODE} show how a DDE can be approximated by changing the fixed delay assumption to an Erlang distributed delay, and then generalized to a phase-type distributed delay. In this more general form, \textbf{v} and \textbf{M} can correspond to any phase-type distribution, not just an Erlang distribution. Furthermore, this assumption could be relaxed even further by allowing entries in \textbf{v} or \textbf{M} to vary over time, or with one or more state variables, as described in \citet{Hurtado2019}.

\section{Results}\label{sec:results}

In the sections below, we illustrate this process of deriving new models using the GLCT by generalizing various biological models taken from the peer reviewed literature. We first highlight some of the key assumptions of these models related to delays and the time spent in different states, viewing each model as a mean field model corresponding to some unspecified stochastic model. We then derive from each model an ODE model that incorporates phase-type distributed delays or dwell times, thereby generalizing or approximating the original model.

\subsection{Model 1: Tumor Growth Inhibition (TGI) Model} \label{sec:TGI}

\begin{figure}[!tbh]
	\centering\includegraphics[trim = 0 450 430 0, clip, width=0.8\textwidth]{TGI.pdf}
	\caption{Schematic diagram of the Tumor Growth Inhibition model (TGI) by \citet{Simeoni2004}. See the main text for further details. \label{fig:TGI}}
\end{figure}

\citet{Simeoni2004} introduced a simple model of tumor growth inhibition that employs the standard Linear Chain Trick (LCT) to incorporate an Erlang distributed time to cell death following tumor cell damage from treatment. The model was subsequently analyzed using standard approaches from dynamical systems \citep{Magni2006,Magni2008}, and has been used elsewhere in the study of tumor growth and the development of cancer treatments \citep[e.g.,][]{Rocchetti2009,Simeoni2013}.  The Simeoni model has also been extended to a ``competing Poisson processes" like assumption (compare Fig. 6 in \citet{Hurtado2019} to Fig. 2 in \citet{Rocchetti2009} and Fig. 1 in \citet{Terranova2013}) in order to model tumor cell death arising from the combined effects of two drugs with no pharmacokinetic interaction.

The basic TGI model is given by eqs. \eqref{eq:GF} and \eqref{eq:TGI} below. In the absence of pharmacological treatment, the amount of cycling (replicating) tumor cells at time $t$ ($z_0(t)$) grows according to the overall growth rate\footnote{This growth rate function is an approximation of the piece-wise function that is equal to $\lambda_0 z_0$ when $w<\lambda_0/\lambda_1$, and $\lambda_1 z_0/w$ when $w\geq \lambda_0/\lambda_1$. See \citet{Magni2006} for details.} \begin{equation} \label{eq:GF}
GF(z_0,w) =\; \frac{\lambda_0\,z_0(t)}{\bigg[1 + \bigg(\frac{\lambda_0}{\lambda_1}\,w(t)\bigg)^\psi \bigg]^\frac{1}{\psi}}.
\end{equation}

Treatment begins at time $t_0>0$, and accordingly the effect of that treatment $c(t)=0$ for $0\leq t\leq t_0$. Once treatment begins, cells that are damaged by the treatment then progress through a series of states Z$_i$, $i=1,\ldots,n$, prior to cell death (see Fig. \ref{fig:TGI}). Together, the full model is given by \begin{subequations}\label{eq:TGI} \begin{align}
\frac{dz_0(t)}{dt} =&\; GF(z_0(t),w(t)) - k_0\,c(t)\,z_0(t) \\
\frac{dz_1(t)}{dt} =&\;  k_0\,c(t)\,z_0(t) - k_1\,z_1(t) \\
\frac{dz_i(t)}{dt} =&\;  k_1\,z_{i-1}(t) - k_1\,z_i(t), \quad\; i = 2,\ldots,n \\
w(t) =&\; \sum_{i=0}^n z_i; \;\; z_0(0)=w_0,\;\; z_i(0)=0, \; i=1,\ldots,n
\end{align}
\end{subequations}

where $w$ is the total amount of tumor cells, and $k_0$ and $c(t)\geq0$ determine the rate of initial tumor cell damage from the treatment. Alternatively, from a more mathematical perspective, $k_0$ and $c(t)$ determine the distribution of time spent in the base state Z$_0$, which follows the first event time distribution under a non-homogeneous Poisson process with rate $r(t)=k_0c(t)$ (see \citet{Hurtado2019} for details). Parameters $n$ and $k_1$ are the shape and rate parameters, respectively, for the Erlang distributed time until cell death for the cells damaged by the treatment. The treatment is assumed to have no effect on the time until cell death after the initial damage to the cell.

To extend this model to instead assume a more general phase-type distributed time to cell death, the equations for $z_i$, $i=1,\ldots,n$ in eqs. \eqref{eq:TGI} can be written in matrix form, using Theorem \ref{th:glct}, where \begin{equation}\textnormal{\textbf{v}} = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0  \end{bmatrix}  \qquad \text{and} \qquad \mathbf{M} = \begin{bmatrix} -k_1 & k_1 & 0 & \cdots & 0 \\
0 & -k_1 & k_1  & \ddots & 0 \\
\vdots & \ddots & \ddots  & \ddots & \ddots \\
0 & 0 & \ddots &  -k_{1} & k_{1} \\
0 & 0 & \cdots &  0 & -k_1
\end{bmatrix}. \end{equation}

This yields the more compact, and more general, set of equations below, where $\mathbf{x}=[z_1,z_2,\ldots,z_n]^\text{T}$. \begin{equation}\label{eq:TGIglct} \begin{split}
\frac{dz_0(t)}{dt} =&\; GF(z_0(t),w(t)) - k_0\,c(t)\,z_0(t) \\
\frac{d\mathbf{x}(t)}{dt} =&\;  k_0\,c(t)\,z_0(t)\,\text{\textbf{v}} + \mathbf{M}^\text{T}\mathbf{x} \\
w(t) =&\; \sum_{i=0}^n z_i \\
z_0(0)=&w_0,\; z_i(0)=0, \; i=1,\ldots,n.
\end{split}
\end{equation}

Note that eqs. \eqref{eq:TGIglct} generalize the TGI model in the sense that these equations accommodate any phase-type distribution assumption for the time to cell death following the initial effect of treatment, not just the Erlang distribution assumed in the original TGI model. Additionally, this matrix-vector form of the original TGI model (i.e., assuming an Erlang distribution) can still be used with some benefit for both computational and mathematical analyses of the TGI model given by eqs. \eqref{eq:TGI}, where those analyses can take advantage of the matrix-vector form of these more general equations \citep{HurtadoRichards2020}.

\subsection{Model 2: Perscription Opioid Epidemic Model}\label{sec:opioid}

The prescription opioid epidemic model by \citet{Battista2019} is a system of ordinary differential equations with no explicit time delays, and (implicit) exponentially distributed dwell times in multiple states. The model assumes individuals are in one of four different states: S, P, A, and R. Here $S$ is the size of the susceptible class. These individuals are not using opioids or recovering from addiction. $P$ is the number of prescribed users (those who are prescribed the drugs and using them but have no addiction). $A$ is the number of addicted individuals who can be using either prescribed or ilicit opioids, and $R$ is the number of individuals undergoing treatment to recover from addiction.

The model describing how individuals transition among these states is given by the following equations, where the dot over each state variable indicates a time derivative. \begin{subequations}\label{eq:opiod}\begin{align}
\dot{S}&= - \alpha S - \beta _A SA - \beta _P SP + \epsilon P + \delta R + \mu (P+R)+\mu^{*}A \\
\dot{P}&= \alpha S - (\epsilon +\gamma +\mu ) P \label{eq:opiodP}\\
\dot{A}&= \gamma P + \sigma R + \beta _A SA + \beta _P SP - (\zeta +\mu ^{*})A \\
\dot{R}&= \zeta A - (\delta +\sigma +\mu ) R.
\end{align}\end{subequations}

The term $\alpha S$ is the rate of individuals transitioning from the susceptible state to the prescribed state after being prescribed opioids per unit time, $\beta _A S A$ is the rate of those transitioning from state S to state A after interacting with addicted individuals, and similarly $\beta _P SP$ represents the rate of individuals who transition from the susceptible class to the addicted class after exposure to opioids via perscription opiod users who have extra or unsecured drugs. The terms $\epsilon P$ and $\delta R$ are the rate individuals leave the prescribed users class without becoming addicted and then reenter the susceptible class at per-capita rate $\epsilon$, and those who leave the rehabilitation state after treatment and reenter the susceptible state at per-capita rate $\delta$. The rates $\mu P$, $\mu R$ and $\mu^*A$ are the death rates for the prescribed, rehabilitated, and addicted classes (to ensure constant population size, deaths are replaced instantaneously by new susceptible individuals). The term  $\gamma P$ is the rate that individuals leave the prescribed class by becoming addicted to their prescription opioids, $\zeta A$ is the rate at which addicted individuals initiate treatment, and $\sigma R$ is the rate at which individuals who are undergoing treatment reenter the addiction class.

Note that from the model terms described above, we may assume that prescription users remain in the prescribed state for an exponentially distributed amount time \citep{Hurtado2019}. Thus, another way to interpret these model terms is that (focusing on eq. \eqref{eq:opiodP}, for example) the proportion of individuals which leave the prescribed user state and go to the susceptible state is $\frac{\epsilon}{\epsilon+\gamma+\mu}$, and thus the net rate of individual entering the susceptible state from the prescribed state is $\frac{\epsilon}{\epsilon+\gamma+\mu}\,(\epsilon+\gamma+\mu)P=\epsilon\,P$. Similarly, the proportion of individuals who go on to become addicted, and who die, are given by $\frac {\gamma}{\epsilon+\gamma+\mu}$ and $\frac {\mu}{\epsilon+\gamma+\mu}$, respectively.

To generalize this model, we aim to replace the implicit assumption of exponentially distributed dwell times in each state, and replace those with the more general phase-type distributions instead.

If we assume the dwell time distribution for the prescribed user state P is a continuous phase-type distribution parameterized by the $n\times 1$ parameter vector $\bf{v_P}$ and $n\times n$  matrix $\bf{M_P}$, then to total number of individuals in state P is given by the sum of the $n$ sub-states $P_i$, $i=1,\ldots, n$. Let ${\bf{x}}=\left[P_1,P_2, \dots, P_n\right]^{\text{T}}$. Then by the GLCT (Theorem \ref{th:glct}), the mean field equations for our prescribed user sub-states are  \begin{equation} \label{eq:dx} \dot{{\bf{x}}} = {\bf{v_P}}\alpha S + {\bf{M_P}}^{\text{T}} {\bf{x}}.\end{equation}

Observe that if we let $\bf{v_P}$ be a one dimensional row vector with its first and only entry being a 1 and let $\bf{M_P}$ be a $1\times 1$ matrix with the entry $-\left(\epsilon +\gamma +\mu \right)$, we recover eq. \eqref{eq:opiodP}.

Recall that individuals who leave the prescribed user state either transition to the addicted state, the susceptible state, or they die. We can denote these proportions as $F_{PA}$, $F_{PS},$ and $F_{PD}$, respectively, where $F_{ij}\in[0,1]$ and $F_{PA}+F_{PS}+F_{PD}=1$. Note that in the original model $F_{PA}=\frac{\gamma}{(\epsilon+ \gamma +\mu)}, F_{PS}=\frac{\epsilon}{(\epsilon+ \gamma +\mu)},$ and $F_{PD}=\frac{\mu}{(\epsilon+ \gamma +\mu)}$. Combining the above with eq. \eqref{eq:dx}, this yields \begin{subequations}\label{eq:opiod2}\begin{align}
\dot{S}=&\; - \alpha S - \beta _A SA - \beta _P SP + F_{PS} (-{\bf{M_P 1}})^{\text{T}} {\bf{x}} + \delta R +F_{PD} (-{\bf{M_P 1}})^{\text{T}} {\bf{x}}+\mu R +\mu ^{*}A \\
\dot{{\bf{x}}}=&\; {\bf{v_P}}\alpha S + {\bf{M_P}}^{\text{T}} {\bf{x}} \\
\dot{A}=&\; F_{PA} (-{\bf{M_P 1}})^{\text{T}} {\bf{x}} + \sigma R + \beta _A SA + \beta _P SP - (\zeta +\mu ^{*})A \\
\dot{R}=&\; \zeta A - (\delta +\sigma +\mu ) R.
\end{align}\end{subequations}

Similarly, we can generalize the addicted and rehabilitated states with phase-type dwell time distributions, assuming the respective phase-type distributions are parameterized by $\mathbf{v_A}$, $\mathbf{M_A}$, $\mathbf{v_R}$, and $\mathbf{M_R}$. Let ${\bf{y}}=\left[A_1,A_2, \dots, A_k\right]^{\text{T}}$ denote the $k$ sub-states of $A$, and ${\bf{z}}=\left[R_1,R_2, \dots, R_m\right]^{\text{T}}$ the $m$ sub-states of $R$. This yields the generalized model:  \begin{subequations} \label{eq:opiodglct} \begin{align} \begin{split}
\dot{S}=&\; - \alpha S - \beta _A SA - \beta _P SP + (F_{PS}+F_{PD}) (-{\bf{M_P 1}})^{\text{T}} {\bf{x}} \;+ \\
        &\; + (F_{RS}+F_{RD}) (-{\bf{M_R 1}})^{\text{T}}\,{\bf{z}}  + F_{AD} (-{\bf{M_A 1}})^{\text{T}}\,{\bf{y}} \end{split} \\
\dot{{\bf{x}}} =&\; {\bf{v_P}}\alpha S + {\bf{M_P}}^{\text{T}} {\bf{x}} \\
\dot{{\bf{y}}}=&\; {\bf{v_A}}\big(F_{PA} (-{\bf{M_P 1}})^{\text{T}} {\bf{x}} +  F_{RA} (-{\bf{M_R 1}})^{\text{T}}{\bf{z}} + \beta _A SA + \beta _P SP\big) + {\bf{M_A}}^{\text{T}}\,{\bf{y}}\\
\dot{{\bf{z}}}=&\; {\bf{v_R}}(F_{AR} (-{\bf{M_A 1}})^{\text{T}}\,{\bf{y}})+ {\bf{M_R}}^{\text{T}} {\bf{z}}.
\end{align} \end{subequations}

It is worth noting that the original model eqs. \eqref{eq:opiod} are a special case of eqs. \eqref{eq:opiodglct}, as are any intermediate extensions of the original model obtained by applying the standard Linear Chain Trick (LCT) to impose Erlang distributed dwell times on one or more of the four main states.

\subsection{Model 3: Within-Host Model of Immune-Pathogen Interactions} \label{sec:inhost}

In \citet{Hurtado2012}, a specific (adaptive) immune response was added to the innate immune response model introduced by \citet{Reynolds2006}. The scaled version of this within-host model, as stated in \citet{Hurtado2012}, is  \begin{subequations}\label{eq:inhost}\begin{align}
	\frac{dp}{dt} =&\; k_{pg}p(1-p) - \frac{k_m p}{\mu_p+p} -K(y)\,n\,p \label{eq:inhostp} \\
	\frac{dn}{dt} =&\; \frac{n+k_p \,p}{x_n + n+k_p \,p} - \mu_n \, p \\
	\frac{dy_0}{dt} =&\; \frac{(np)^\alpha}{x_y^\alpha + (np)^\alpha} - \mu_{y0}\,y_0 \\
	\frac{dy}{dt} =&\; \mu_{y0}\,y_0  - \mu_y\,y
	\end{align}
\end{subequations}

In this model, $p$ is the scaled pathogen (bacteria) population size, which follows a logistic growth model in the absence of an immune response. The second term in eq. \eqref{eq:inhostp} models the effect of some baseline local immune defenses capable of neutralizing a small population of pathogen, and mathematically introduces a strong Allee effect into the model. The level of innate immune activity $n$ increases in response to the presence of pathogen, as well as from a positive feedback loop, and the interaction of this innate immune activity and pathogen stimulates progenitor cells ($y_0$) that mature into active specific immune components ($y$), e.g., B-cells, which augment the pathogen-killing capacity of the innate immune components (i.e., which increase $K(y)$).  For further details on this model, see \citet{Hurtado2012} and \citet{Reynolds2006}.

In this model, the delay in activating the specific immune response can be thought of as an exponentially distributed maturation time (with mean $1/\mu_{y0}$) and the duration of the active immune response (i.e., the dwell time of mature specific immune components modeled by $y$) is also exponentially distributed (with mean $1/\mu_y$).

Both of these dwell time distribution assumptions can be replaced by phase-type distributions with respective parameters $\mathbf{v_{y0}}$, $\mathbf{M_{y0}}$, $\mathbf{v_{y}}$ and $\mathbf{M_y}$, respectively.  To do this, we first partition state Y$_0$ into sub-states X$_i$, $i=1,\ldots,m$, and the state Y into sub-states Z$_j$, $j=1,\ldots,n$, where $y_0=\sum_{i=1}^m x_i$ and $y=\sum_{j=1}^n z_j$, and we let  $\mathbf{x}=[x_1,\ldots,x_m]^\text{T}$ and $\mathbf{z}=[z_1,\ldots,z_n]^\text{T}$. The GLCT (Theorem \ref{th:glct}) then yields the more general model \begin{subequations}\label{eq:inhostglct}\begin{align}
	\frac{dp}{dt} =&\; k_{pg}p(1-p) - \frac{k_m p}{\mu_p+p} -K(y)\,n\,p  \\
	\frac{dn}{dt} =&\; \frac{n+k_p \,p}{x_n + n+k_p \,p} - \mu_n \, p \\
	\frac{d\mathbf{x}}{dt} =&\; \mathbf{v_{y0}}\frac{(np)^\alpha}{x_y^\alpha + (np)^\alpha} + \mathbf{M_{y0}}^\text{T}\,\mathbf{x} \\
	\frac{d\mathbf{z}}{dt} =&\; -\mathbf{1}^\text{T}\mathbf{M_y}^\text{T}\mathbf{v_y} + \mathbf{M_y}^\text{T}\mathbf{z}.
	\end{align}
\end{subequations}

Note that the above model could be similarly be extended to include a phase-type distributed time lag in the activation of the non-specific immune response (modeled by $n$). However, given the relatively fast activation time of this response, we have omitted this extension in the above model.

\subsection{Model 4: Cell-To-Cell Spread of HIV}\label{sec:HIV}

\citet{Culshaw2003} introduced an integro-differential model for the cell-to-cell spread of HIV, which incorporates a distributed time delay in the time between cells becoming infected and infectious. They then derive from this general model multiple other models which differ only in the specific assumptions on the form of this delay distribution.

In their most general model, state variable $C(t)$ represents the concentration of healthy cells at time $t$, and $I(t)$ is the concentration of infected cells. The model is as follows: \begin{subequations} \label{eq:CIide}\begin{align}
	\frac {dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
	\frac {dI}{dt} =&\;  k'_i\int_{-\infty}^{t} C(u)I(u)F(t-u)du-\mu_{I}I(t).
	\end{align} \end{subequations}

Parameter $r_C$ is the net growth rate of the healthy cell population, $C_M$ is an effective carrying capacity of the system, $k_I$ is an infection rate parameter, $k_I'/k_I$ is the fraction of cells surviving the incubation period, and $\mu_I$ is the per capita death rate of infected cells (implicitly, the infected cell lifetime is exponentially distributed with mean $1/\mu_I$). Initial values for $C$ and $I$ must be functions defined over all $s\in (-\infty,0]$ and are denoted $\phi(s)\geq 0$ and $\psi(s)\geq0$, respectively.

In \citet{Culshaw2003}, the delay kernel $F(u)$ is assumed to be of the form \begin{equation} \label{eq:CIF} F(u) =\; \frac{\alpha^{n+1}u^n}{n!}e^{-\alpha\,u} \end{equation} which is just the density function for an Erlang distribution with rate $\alpha$ and shape $n+1$ (and thus, mean $(n+1)/\alpha$ and coefficient of variation $1/\sqrt{n+1}$), and the \textit{weak} and \textit{strong} kernels are just the particular cases where the shape parameter is 1 (i.e., an exponential distribution with rate $\alpha$) or 2 (Erlang with rate $\alpha$ and shape 2), respectively.

Three models are then derived in \citet{Culshaw2003} from this more general integro-differential equation model, which we will summarize here then extend further using the LCT and GLCT (see \citet{Culshaw2003} for a comparison of the dynamics of these three models).

First, assuming $F(u)=\delta(u)$ is a Dirac delta function (point mass at zero) yields the model with no delay, given in \citet{Culshaw2003} by equations  \begin{equation}\begin{split}
\frac {dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac {dI}{dt} =&\;  k'_i C(t)I(t) - \mu_{I}I(t),
\end{split} \end{equation}

with initial conditions $C(0)=c_0\geq 0$ and $I(0)=I_0\geq 0$.

Second, assuming $F(u)=\delta(u-\tau)$ is a Dirac delta function at time $\tau>0$ yields the delay differential equation below (as written in \citet{Culshaw2003}) with the same initial conditions as eqs. \eqref{eq:CIide}. \begin{equation}\begin{split}
\frac {dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac {dI}{dt} =&\;  k'_i C(t-\tau)I(t-\tau) - \mu_{I}I(t).
\end{split} \end{equation}

Third, assuming a ``weak kernel" (i.e., exponentially distributed delay with rate $\alpha$) yields \begin{equation}\label{eq:CXI}\begin{split}
\frac{dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac{dX}{dt} =&\;  \alpha\,C(t)I(t) - \alpha\,X(t) \\
\frac{dI}{dt} =&\;  k'_I\,X(t) - \mu_{I}I(t).
\end{split} \end{equation}

Observe that, by making the simple substitution $Y(t)=\frac{k_I}{\alpha}X(t)$, we can write the following alternative model which is equivalent to eqs. \eqref{eq:CXI} by Culshaw \textit{et al}, but is more natural in terms of the units of $X$ and $Y$ and in the context of the LCT: \begin{equation} \label{eq:CYI} \begin{split}
\frac{dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac{dY}{dt} =&\;  k_I\,C(t)I(t) - \alpha\,Y(t) \\
\frac{dI}{dt} =&\;   \frac{k'_I}{k_I} \alpha Y(t) - \mu_{I}I(t).
\end{split} \end{equation}

From eqs. \eqref{eq:CYI}, it is straightforward to derive two additional models using the LCT and GLCT. Using the LCT, eqs. \eqref{eq:CIlct} below are a more general form of eqs. \eqref{eq:CYI} that correspond to any choice of a non-negative integer value of $n$ and $\alpha>0$ for the delay kernel $F$ in eq. \eqref{eq:CIF} (i.e., any Erlang distribution with shape $n+1$ and rate $\alpha$).  \begin{equation} \label{eq:CIlct} \begin{split}
\frac{dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac{dY_1}{dt} =&\;  k_I\,C(t)I(t) - \alpha\,Y_1(t) \\
\frac{dY_i}{dt} =&\;  \alpha\,Y_{i-1}(t) - \alpha\,Y_i(t), \quad\; i=2,\ldots,n+1 \\
&\vdots \\
\frac{dI}{dt} =&\;   \frac{k'_I}{k_I} \alpha Y_{n+1}(t) - \mu_{I}I(t).
\end{split} \end{equation}

Using the above equations as a guide, re-writing them in matrix form as suggested by the GLCT (Theorem \ref{th:glct}), yields the more general set of equations below, which are the desired set of model equations for which the Erlang distribution assumption (with parameters $n+1$ and $\alpha$) has been replaced by a phase-type distribution parameterized by the length $k$ vector \textbf{v} and $k\times k$ matrix \textbf{M}, where $\mathbf{y}=[Y_1,\ldots,Y_k]^\text{T}$. \begin{equation} \label{eq:CIglct} \begin{split}
\frac{dC}{dt} =&\;  r_C C(t) \left(1-\frac{C(t)+I(t)}{C_M}\right) -k_IC(t)I(t) \\
\frac{d\mathbf{y}}{dt} =&\;  k_I\,C(t)I(t)\,\mathbf{v} + \mathbf{M}^\text{T}\,\mathbf{y}\\
\frac{dI}{dt} =&\;  -\frac{k'_I}{k_I}\mathbf{1}^\text{T}\,\mathbf{M}^\text{T}\,\mathbf{y} - \mu_{I}I(t).
\end{split} \end{equation}

Further generalizations, e.g., to time-varying \textbf{v}, \textbf{M}, or survival fraction $f=k'_I/k_I$, are also possible \citep{Hurtado2019}.

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

Mean field state transition models, written as ordinary differential equations (ODEs), are widely used throughout the sciences, but too often they include overly simplistic assumptions regarding time delays and the duration of time individual entities spend in specific states. In this paper, we illustrate how those assumptions can be refined within the context of the GLCT, and used to derive new models. We have derived multiple new dynamical systems models based on existing models taken from the literature, to illustrate the relative ease of deriving such models using the GLCT. These examples span a range of biological application areas and different model types (DDEs, ODEs, and integro-differential equations), which reflect a mix of different implicit and explicit delay assumptions. Using the GLCT, those delay assumptions were replaced or generalized to yield new ODEs that incorporate phase-type distributed delays and dwell times.

These straightforward generalizations illustrate how modelers can incorporate phase-type distributed delays and dwell times into ODE models, and how those underlying (implicit) stochastic model assumptions are reflected in the corresponding mean field ODE model structure. Importantly, these alternative model formulations can also be used in the computational and mathematical analysis of models that only assume Erlang distributions and could otherwise be derived using the standard linear chain trick (LCT). For an example, see \citet{HurtadoRichards2020}, where we illustrate the potential computational benefits of using a GLCT formulation of models with Erlang dwell time assumptions when computing numerical solutions to such models.

These generalized models also lay the groundwork for incorporating Coxian, hyperexponential, hypoexponential (i.e, generalized Erlang) and other phase-type distributions into these and similar mean field ODE models. Statistical tools such as \texttt{BuTools} \citep{BuTools2,BuToolswww} allow modelers to fit phase-type distributions to data, thereby allowing modelers to build approximate empirical distributions into ODE models using the GLCT. However, it is important to note that there are certain limitations to approximating some delay or dwell time assumptions with phase-type distributions. For example, delay distributions with compact support (e.g., a continuous uniform distribution) may not be well approximated by phase-type distributions.

In closing, we hope these new techniques prove to be helpful to modelers in their efforts to build better models, to check the consequences of certain simplifying assumptions, and to gain better intuition for how underlying assumptions are reflected in the structure of ODE model equations.

\subsection*{Acknowledgements}
The authors thank Dr. Deena Schmidt for conversations and comments that improved this manuscript. The authors also thank the organizers and sponsors of the Second International Conference on Applications of Mathematics to Nonlinear Sciences (AMNS-2019), held June 27-30, 2019, in Pokhara, Nepal, and the editors of this Thematic Issue.

\subsection*{Funding}
This work was supported by a grant awarded to PJH by the Sloan Scholars Mentoring Network of the Social Science Research Council with funds provided by the Alfred P. Sloan Foundation; and this material is based upon work supported by the National Science Foundation under Grant No. DEB-1929522.

\subsection*{Disclosure statement}
The authors declare that they have no conflict of interest.

\printbibliography

\end{document}


\section{Introduction}
The goal of semantic segmentation is to assign a semantic category to each pixel in the image. It has been one of the most important tasks in computer vision that enjoys a wide range of applications such as image editing and scene understanding.
Recently, deep convolutional neural network (CNN) based methods \cite{fcn_pami,deeplab,dilated} have been developed for semantic segmentation and achieved significant progress.
However, such approaches rely on learning supervised models that require pixel-wise annotations, which take extensive effort and time.
To reduce the effort in annotating pixel-wise ground truth labels, numerous weakly-supervised methods are proposed using various types of labels such as image-level \cite{ahn2018learning,kolesnikov2016seed,pathak2015constrained,pinheiro2015image}, video-level \cite{Chen_IJCV_2020,Zhong_ACCV_2016,Tsai_ECCV_2016}, bounding box \cite{papandreou2015weakly,dai2015boxsup,khoreva_CVPR17}, point-level \cite{Bearman_ECCV16}, and scribble-based \cite{lin2016scribblesup,Vernaza_CVPR17} labels.
In this work, we focus on using image-level labels which can be obtained effortlessly, yet a more challenging case under the weakly-supervised setting.

Existing algorithms mainly consist of three sequential steps to perform weakly-supervised training on the image-level label: 1) predict an initial category-wise response map to localize the object, 2) refine the initial response as the pseudo ground truth, and 3) train the segmentation network based on pseudo labels.
Although promising results have been achieved by recent methods \cite{ahn2018learning,huang2018weakly,wang2018weakly,wei2018revisiting}, most of them focus on improving the second and the third steps. Therefore, these approaches may suffer from inaccurate predictions generated in the first step, i.e., initial response.
Here, we aim to improve the performance of initial predictions which will benefit succeeding steps.

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{teaser.pdf}\\
	\caption{Existing weakly-supervised semantic segmentation methods based on image-level supervisions usually apply the class activation map (CAM) to obtain the response map as the initial prediction.
	%
	However, this response map can only highlight the discriminative parts of the object (top).
	%
	We propose a self-supervised task via sub-category exploration to enforce the classification network learn better response maps (bottom).
	}
	\label{fig:teaser}
 	\vspace{-3mm}
\end{figure}

In order to predict the initial response map for each category, numerous approaches based on the
class activation map (CAM) model \cite{zhou2016learning}
have been developed.
Essentially, these methods train a classification network and use its learned weights in the classifier as the cues to compute weighted sums of feature maps, which can be treated as the response map.
However, such response maps may only focus on a portion of the object, instead of localizing the entire object (see top of Figure \ref{fig:teaser}).
One explanation is that the objective of the classifier does not need to ``see'' the entire object for optimizing the loss function.
This impairs the classifier's ability to locate the objects.

At the core of our technique is to impose a more challenging task to the network for learning better representations, while not jeopardizing the original objective.
To this end, we propose a simple yet effective method by introducing a self-supervised task that discovers sub-categories in an unsupervised manner, as illustrated at the bottom of Figure \ref{fig:teaser}.
Specifically, our task consists of two steps: 1) perform clustering on image features extracted from the classification network for each annotated parent class (e.g., 20 parent classes on the PASCAL VOC 2012 dataset \cite{PASCAL_VOC_2010_Data}), and 2) use the clustering assignment for each image as the pseudo label to optimize the sub-category objective.

On one hand, the parent classifier establishes a feature space through supervised training as the guidance for unsupervised sub-category clustering.
On the other hand, the sub-category objective provides additional gradients to enhance feature representations and leverage the sub-space of the original feature space to obtain better results.
As such, the classification model takes a more challenging task and is not limited to the easier objective of learning only the parent classifier.
Moreover, to ensure better convergence in practice, we iteratively alter the two steps of feature clustering and pseudo training the sub-category objective.

We conduct extensive experiments on the PASCAL VOC 2012 dataset \cite{PASCAL_VOC_2010_Data} to demonstrate the effectiveness of our method, with regard to generating better initial response maps to localize objects.
As a result, our approach leads to favorable performance for the final semantic segmentation results against state-of-the-art weakly-supervised approaches.
Furthermore, we provide extensive ablation studies and analysis to validate the robustness of our method.
Interestingly, we notice that the network is able to differentiate sub-categories with respect to their object size/type, context, and coexistence with other categories.
The main contributions of this work are summarized as follows:
\begin{itemize}
    \item We propose a simple yet effective method via a self-supervised task to enhance feature representations in the classification network. This  improves the initial class activation maps for weakly-supervised semantic segmentation as well.

    \item We explore the idea of sub-category discovery via iteratively performing unsupervised clustering and pseudo training on the sub-category objective in a self-supervised fashion.

    \item We present extensive study and analysis to show the efficacy of the proposed method, which significantly improves the quality of initial response maps and leads to better semantic segmentation results.

\end{itemize}

\begin{figure*}[t]
	\centering
	\includegraphics[width=0.95\linewidth]{framework.pdf}\\
	\caption{
		Proposed framework for generating the class activation map.
		%
		Given input images $I$, we first feed them into a feature extractor $E$ to obtain their features $f$.
		%
		Then, we adopt unsupervised clustering on $f$ and obtain sub-category pseudo labels $Y_s$ for each image.
		%
		Next, we train the classification network to jointly optimize the parent classifier $H_p$ with ground truth labels $Y_p$ for parent classes and the sub-category classifier $H_s$ using the sub-category pseudo labels obtained in the clustering stage.
		%
		By iteratively performing unsupervised clustering on image features and pseudo training the classification module, we use the jointly optimized classification network to produce the final activation map $M$.
	}
	\label{fig: framework}
	\vspace{-3mm}
\end{figure*}

\section{Related Work}
Within the context of this work, we discuss methods for weakly-supervised semantic segmentation (WSSS) using image-level labels, including approaches that focus on initial prediction and refinement for generating pseudo ground truths. In addition, algorithms that are relevant to unsupervised representation learning are discussed in this section.

\paragraph{Initial Prediction for WSSS.}

Initial cues are essential for segmentation task since it can provide reliable priors to generate segmentation maps. The class activation map  \cite{zhou2016learning} is a widely used technique for localizing the object. It can highlight class-specific regions that are served as the initial cues. However, since the CAM model is trained by a classification task, it tends to activate to the small discriminative part of the object, leading to incomplete initial masks.

Several methods have been developed to alleviate this problem. Numerous approaches \cite{singh2017hide, wei2017object} deliberately hide or erase the region of an object, forcing models to seek more diverse parts.
However, those methods either hide fixed-size patches randomly or require repetitive model training and response aggregation steps.
A number of variants \cite{zhang2018adversarial, li2018tell} have been proposed to extend the initial response via an adversarial erasing strategy in an end-to-end training manner, yet such strategies may gradually expand their attention to non-object regions, leading to inaccurate attention maps.
Recently, the SeeNet approach \cite{hou2018self} applies self-erasing strategies to encourage networks to use both object and background cues, which prevent the attention from including more background regions.
Instead of using the erasing scheme, the FickleNet method \cite{lee2019ficklenet} introduces stochastic feature selection to obtain diverse combinations of locations on feature maps. By aggregating the localization maps, they acquire the initial cue that contains a larger region of the object.

Different from the methods that mitigate the problem by discovering complementary regions via iterative erasing steps or consolidating attention maps, our proposed approach aims at enforcing the network to learn harder on a more challenging task via self-supervised sub-category exploration, thereby enhancing feature representations and improving the response map.

\paragraph{Response Refinement for WSSS.}
Numerous approaches \cite{ahn2018learning,fan2018cian, fan2018associating,huang2018weakly,kolesnikov2016seed,wang2018weakly,wei2018revisiting} are proposed to refine the initial cue via expanding the region of attention map.
The SEC method \cite {kolesnikov2016seed} proposes a loss function that constrains both global weighted rank pooling and low-level boundary to expand the localization map.
To improve the network training, the MCOF scheme \cite{wang2018weakly} uses a bottom-up and top-down framework which alternatively expands object regions and optimize the segmentation network, while the MDC method \cite{wei2018revisiting} expands the seeds by employing multiple branches of convolutional layers with different dilation rates.
Moreover, the DSRG approach \cite{huang2018weakly} refines initial localization maps by applying a seeded region growing method during the training of the segmentation network.
Other approaches are developed via affinity learning.
For instance, the AffinityNet \cite{ahn2018learning} considers pixel-wise affinity to propagate local responses to nearby areas, while \cite{fan2018cian, fan2018associating} explore cross-image relationships to obtain complementary information that can infer the predictions.

Nevertheless, initial seeds are still obtained from the CAM method. If these seeds only come from the discriminative parts of objects, it is difficult to expand regions into non-discriminative parts. Moreover, if the initial prediction produces wrong attention regions, applying the refinement step would cover even more inaccurate regions.
In this paper, we focus on improving the initial prediction, which leads to more accurate object localization and benefits the refinement step.

\paragraph{Unsupervised Representation Learning.}
Unsupervised learning has been widely studied in the computer vision community. One advantage is to learn better representations of images and apply learned features on any specific domain or dataset where annotations are not always available.
Self-supervised learning \cite{de1994learning} utilizes a pretext task to replace the labels annotated by humans with ``pseudo-labels'' directly computed from the raw input data.
A number of methods \cite{noroozi2016unsupervised, pathak2016context, paulin2015local} are developed but require expert knowledge to carefully design a pretext task that may lead to good transferable features.
To reduce the domain knowledge requirement, Coates and Ng \cite{coates2012learning} validate that feature-learning systems with K-means can be a scalable unsupervised learning module that can train a model of the unlabeled data for extracting meaningful features.
Furthermore, a recent approach \cite{caron2018deep} employs a clustering framework to extract useful visual features by
alternating between clustering the image descriptors and updating the weights of the CNN by predicting the cluster assignments, in order to learn deep representations specific to domains where annotations are scarce.
In this work, we propose to learn a self-supervised method that explores the sub-category in the classification network, i.e., using unsupervised signal to enhance feature representations while improving initial response maps for weakly-supervised semantic segmentation.

\section{Weakly-supervised Semantic Segmentation}
In this section, we describe our framework for weakly-supervised semantic segmentation, including details of how we explore sub-categories to improve initial response maps and generate final semantic segmentation results.

\subsection{Algorithm Overview}

To obtain the initial response, we follow the common practice of training a classification network and utilize the CAM method \cite{zhou2016learning} to obtain our baseline model.
The CAM method typically only activates on discriminative object parts, which are not sufficient for the image classification task. To address this issue, We propose to integrate a more challenging task into the objective: self-supervised sub-category discovery, in order to enforce the network to learn from more object parts.

Firstly, for each annotated parent class, we determine $K$ sub-categories by applying K-means clustering on image features.
With the clustering results, we then assign each image with a pseudo label, which is identified as the index of the sub-category. Finally, we construct a sub-category objective to jointly train the classification network.
By iteratively updating the feature extractor, two classifiers, and sub-category pseudo labels, the enhanced features representations lead to better classification, and thereby gradually produce response maps that attain to more complete regions of the objects.
The overall process is illustrated in Figure \ref{fig: framework}.
Then, we use the method in \cite{ahn2018learning} to expand response maps, which are used as pseudo ground truths to train the segmentation network.
Also note that, our method focuses on the initial prediction, so it is not limited to certain region expansion or segmentation training methods.

\paragraph{Preliminaries: Initial Response via CAM.}
We adopt the CAM to generate the initial response using a typical classification network, whose architecture consists of convolutional layers as the feature extractor $E$, followed by global average pooling (GAP) and one fully-connected layer $H_p$ as the output classifier.
Given an input image $I$, the network is trained with image-level labels $Y_p$ using a multi-label classification loss $\mathcal{L}_p$, following \cite{zhou2016learning}.
After training, the activation map $M$ for each category $c$ can be obtained via directly applying classifier $H_p$ on the feature maps $f = E(I)$:
\begin{equation}
    M^c(x,y) = \theta_p^{c \top} f(x,y),
    \label{eq:cam}
\end{equation}
where $\theta_p^c$ is the classifier weight for the category c, and $f(x,y)$ is the feature at pixel $(x,y)$. The response map is further normalized by the maximum value in $M^c$.

\subsection{Sub-category Exploration}
The activation map for each image using \eqref{eq:cam} provides typically highlights only the discriminative object parts.
However, from the perspective of a classifier, discovering the most discriminative part of the object is already sufficient for optimizing the loss function $\mathcal{L}_p$ in classification.
As the learning objective is based on the classification scores, it is inevitable for the CAM model to generate
incomplete attention maps.
To address this issue, we integrate a self-supervised scheme to enhance feature representations $f$ while improving the response maps via exploring the sub-category information, in which $f$ appears to be an important cue to compute the activation map via \eqref{eq:cam}.

\paragraph{Sub-Category Objective.}
To assign a more challenging problem to the classification model, we introduce a task to discover sub-categories in an unsupervised manner.
For each parent class $p_c$, we define $K$ sub-categories $s_c^k$, where $k = \{1,2, ... ,K\}$.
For each image $I$ with the parent label $Y_p^c$ in $\{0, 1\}^{c}$, the corresponding sub-category label for the category $c$ is denoted as $Y_s^{c,k}$ in $\{0, 1\}^{k}$.
We also note that, if the label of one parent class does not exist (i.e., $Y_p^c = 0$), the labels of all sub-categories would be also 0, i.e., $Y_s^{c,k} = 0, k = \{1,2,...,K\}$.
Our objective is to learn a sub-category classifier $H_s$ parameterized with $\theta_s$, while sharing the same feature extractor $E$ with $H_p$.
Similar to the parent classification loss $\mathcal{L}_p$, we adopt the standard multi-label classification loss $\mathcal{L}_s$ with a larger and fine-grained label space $Y_s$.

\vspace{-3mm}
\paragraph{Sub-category Discovery.}
As there is no ground truth label for sub-category to directly optimize the above sub-category objective $\mathcal{L}_s$, we generate pseudo labels via unsupervised clustering.
Specifically, we perform clustering for each parent class on image features extracted from the feature extractor $E$.
The clustering objective for each class $c$ can be written as:
\begin{equation}
    \min_{D \in \mathbb{R}^{d \times k}} \frac{1}{N^c} \sum_{i=1}^{N^c} \min_{Y_s^c} ||f - TY_s^c||_2^2, \;\;\text{s.t.}, Y_{s}^{c \top} 1_k = 1,
    \label{eq:clustering}
\end{equation}
where $T$ is a $D \times K$ centroid matrix, $N^c$ is the number of images containing the class $c$, and $f = E(I) \in \mathbb{R}^D$ is the extracted feature.
We use the clustering assignment $Y_s^c$ for each image as the sub-category pseudo label to optimize $\mathcal{L}_s$.

\vspace{-3mm}
\paragraph{Joint Training.}
After obtaining sub-category pseudo labels $Y_s$ from the above clustering process, we jointly optimize the feature representations $f = E(I)$ and two classifiers, i.e., $H_p$ and $H_s$:
\begin{equation}
    \min_{\theta_p, \theta_s} \frac{1}{N} \sum_{i=1}^{N} \mathcal{L}_p(H_p(f_i), Y_p) + \lambda \mathcal{L}_s(H_s(f_i), Y_s),
    \label{eq:joint}
\end{equation}
where $N$ is the total number of images and $\lambda$ is weight to balance two loss functions.
With this method, the parent classification learns a feature space through supervised training via $\mathcal{L}_p$, while the sub-category objective $\mathcal{L}_s$ explores the feature sub-space and provides additional gradients to enhance feature representations $f$, which is used to compute CAM via \eqref{eq:cam}.

\begin{algorithm}[!t]
\vspace{0.1in}
\caption{Learning Sub-category Discovery for CAM}\label{alg:wwws}
\begin{algorithmic}
\State \textbf{Input:} Image $I$; Parent Label $Y_p$; Category Number $C$; \\ Sub-category Number $K$
\State \textbf{Output:} Class Activation Map $M^c$\\

\State \textbf{Model:} Feature extractor $E$; Parent Classifier $(H_p; \theta_p)$; \\ Sub-category Classifier $(H_s; \theta_s)$ \\

\State Optimize $\{E, H_p\}$ with $Y_p$ via $\mathcal{L}_p$
\While{Training}

\State Extract features via $f = E(I)$

\For{$c \gets 1$ to $C$}
\State Generate pseudo labels $Y_s^c$ with $f$ via \eqref{eq:clustering}
\EndFor

\State Optimize $\{E, H_p, H_s\}$ with $\{Y_p, Y_s\}$ via \eqref{eq:joint}

\EndWhile
\\
Compute $M^c$ via \eqref{eq:cam}
\end{algorithmic}
\end{algorithm}
\vspace{-1mm}

\paragraph{Iterative Optimization.}
The proposed unsupervised clustering scheme in \eqref{eq:clustering} relies on the feature $f$ to discover sub-category pseudo labels.
As such, the learned features via only the objective $\mathcal{L}_p$ could be less discriminative for the clustering purpose.
To mitigate this issue, we adopt an iterative training method by alternatively updating \eqref{eq:clustering} and \eqref{eq:joint}.
Therefore, features $f$ are first enhanced through the sub-category objective, and in turn facilitate the clustering process to generate better pseudo ground truths, which are then used to learn better feature representations in network training.
The overall optimization for generating final class activation maps is summarized in Algorithm \ref{alg:wwws}.

\subsection{Implementation Details}
In this section, we describe implementation details of the proposed framework and the following procedures to produce final semantic segmentation results.
All the source code and trained models are available at \url{https://github.com/Juliachang/SC-CAM}.

\paragraph{Classification Network.}
In this work, the ResNet-38 architecture \cite{wu2019wider} is used for the CAM model, and the training procedure is similar to that in \cite{ahn2018learning}.
The network consists of 38 convolution layers with wide channels, followed by a $3 \times 3$ convolution layer with 512 channels for better adaptation to the classification task, a global average pooling layer for feature aggregation, and two fully-connected layers for image and sub-category classification, respectively.
The model is pre-trained on the ImageNet \cite{deng2009imagenet} and is then finetuned on the PASCAL VOC 2012 dataset.
We use the typical techniques based on the horizontal flip, random cropping, and color jittering operations to augment the training data set.
We also randomly scale input images to impose scale invariance in the network.

We implement the proposed framework with PyTorch and train on a single Titan X GPU with 12 GB memory.
To train the classification network, we use the Adam optimizer \cite{kingma2014adam} with initial learning rate of 1e-3 and the weight decay of 5e-4.
In practice, we use $\lambda = 5$ and $K=10$ in all the experiments unless specified otherwise.
For iterative training, we empirically find that the model converges after training for 3 rounds.
In the experimental section, we show studies for the choice of $K$ and iterative training results.

\begin{table}[!t]
	\caption{Performance comparison in mIoU (\%) for evaluating activation maps on the PASCAL VOC training and validation sets.}
	\vspace{-1mm}
	\small
	\centering
	\renewcommand{\arraystretch}{1.1}
	\setlength{\tabcolsep}{6pt}
	\begin{tabular}{lcc|cc}
		    \toprule
	        & \multicolumn{2}{c}{Training Set} & \multicolumn{2}{c}{Validation Set} \\
		    \midrule
		    Method & CAM & CAM+RW & CAM & CAM+RW \\
		    \midrule
		    AffinityNet \cite{ahn2018learning} & 48.0 & 58.1 & 46.8 & 57.0 \\
		    Ours & 50.9 & 63.4 & 49.6 & 61.2 \\
		\bottomrule
	\label{table: compare_cam}
	\end{tabular}
	\vspace{-8mm}
\end{table}

\paragraph{Semantic Segmentation Generation.}
Based on the response map generated by our method as in Algorithm \ref{alg:wwws}, we adopt the random walk method via affinity \cite{ahn2018learning} to refine the map as pixel-wise pseudo ground truths for semantic segmentation.
In addition, as a common practice, we use dense conditional random fields (CRF) \cite{crf} to further refine the response to obtain better object boundaries.
To train the segmentation network, we utilize the Deeplab-v2 framework \cite{deeplab} with the ResNet-101 architecture \cite{He_2016_CVPR} as the backbone model.

\begin{figure*}[t]
	\centering
	\includegraphics[width=1\linewidth]{response_compare.pdf}\\
	\caption{Sample results of initial responses.
	%
	Our method often generates the response map that covers larger region of the object (i.e., attention on the body of the animal), while the response map produced by CAM \cite{zhou2016learning} tends to highlight small discriminative parts.}
	\label{fig: cam_viz}
	\vspace{-1mm}
\end{figure*}

\begin{figure}[!t]
	\centering
	\includegraphics[width=0.85\linewidth]{k.pdf}\\
	\vspace{-1mm}
	\caption{Ablation study for $K$. We show that the proposed method performs robustly with respect to $K$ and is consistently better than the original CAM that did not apply clustering to discover sub-categories. We mark the value of mIoU of the original CAM at $K = 1$ and the improved mIoUs are presented.}
	\label{fig: k}
	\vspace{-3mm}
\end{figure}

\section{Experimental Results}
In this section, we first present the main results and analysis of the initial response generated by our method.
Second, we show the final semantic segmentation performance on the PASCAL VOC dataset \cite{PASCAL_VOC_2010_Data} against the state-of-the-art approaches.
More results can be found in the supplementary material.

\subsection{Evaluated Dataset and Metric}
We evaluate the proposed approach on the PASCAL VOC 2012 semantic segmentation benchmark \cite{PASCAL_VOC_2010_Data} which contains 21 categories, including one background class.
Each image contains one or multiple object classes. Following previous weakly-supervised semantic segmentation methods, we use augmented 10,528 training images present in \cite{hariharan2011semantic} along with their image-level labels to train the network.
To evaluate the training set, we use the set without augmentation which has 1,464 examples.
We adopt 1,449 images in the validation set and 1,456 images in the test set to compare our results with other methods.
For all experiments, the mean Intersection-over-Union (mIoU) ratio is used as the evaluation metric.
The results for the test set are obtained from the official PASCAL VOC evaluation website.

\subsection{Improvement on Initial Response}
In Table \ref{table: compare_cam}, we show the mean IoU of the segments computed using the CAM on both the training and validation sets.
We present results after applying the refinement step to the activation map, i.e., CAM + random walk (CAM + RW).
Table \ref{table: compare_cam} shows that our approach significantly improves the IoU over AffinityNet \cite{ahn2018learning} by almost 3$\%$ using CAM and more than 4$\%$ for CAM+RW.
The improved initial response maps facilitate the downstream task in generating pixel-wise pseudo ground truths for training the semantic segmentation model.

In Figure \ref{fig: cam_viz}, we show comparisons of generated CAMs by the conventional classification loss $\mathcal{L}_p$ \cite{zhou2016learning} and the proposed method via sub-category discovery summarized in Algorithm \ref{alg:wwws}.
Visual results show that our method is able to localize more complete object regions, while the original CAM only focuses on discriminative object parts.
We also note that this is essentially critical for the refinement stage that takes the response map as the input.

\begin{table}[!t]
	\caption{Segmentation quality of the initial response at different rounds of training on the PASCAL VOC 2012 validation set. We show there is a gradual improvement on both mIoU and F-Score metrics.}
	\vspace{-2mm}
	\small
	\centering
	\renewcommand{\arraystretch}{1.1}
	\setlength{\tabcolsep}{6pt}
	\begin{tabular}{lcc}
		\toprule
		Round & mIoU (\%) $\uparrow$ & F-Score $\uparrow$ \\
		\midrule

		\#0 (CAM) & 46.8 & 65.1 \\
		\#1 & 48.0 & 65.6 \\
		\#2 & 48.7 & 66.6 \\
		\#3 & \textbf{49.6} & \textbf{67.0} \\

		\bottomrule
	\label{table: iter_improve}
	\end{tabular}
	\vspace{-5mm}
\end{table}

\subsection{Ablation Study and Analysis}
To demonstrate how our method helps improve feature representations and allow the network pay more attention to other object parts via exploiting the sub-category information, we present extensive analysis in this section.
Here, all the experimental results are based on the PASCAL VOC validation set.

\vspace{-2mm}
\paragraph{Effect of Sub-Category Number $K$.}
We first study how the sub-category number $K$  affect the performance of the proposed method.
In Figure \ref{fig: k}, we use $K = \{5, 10, 20, 50\}$, and show that the proposed method performs robustly with respect to $K$ (within a wide range) and consistently better than the original CAM method (i.e., $K = 1$).
The results also validate the necessity and importance of using more sub-categories (i.e., $K>1$) to generating better response maps.
Considering the efficiency and accuracy, we use $K = 10$ for each parent class in all the experiments.
\textcolor{black}{As a future work, it is of great interest to develop an adaptive method to determine the sub-category number \cite{sarfraz2019efficient}, which can reduce the redundant sub-categories and make the approach more efficient.}

\vspace{-2mm}
\paragraph{Iterative Improvement.}
To demonstrate the effectiveness of our iterative training process, we show the gradual improvement on the segment quality in Table \ref{table: iter_improve}.
We present the results of mIoU and F-Score that accounts for both the recall and precision measurements, in which they are important cues to validate whether the activation map is able to cover object parts.
Compared to the results in round \#0, which is the original CAM, our method gradually improves both metrics as training more rounds.

\begin{figure*}[t]
	\centering
	\includegraphics[width=1\linewidth]{clustering_viz.pdf}\\
	\vspace{-2mm}
	\caption{Clustering results of the last round model (\#3). We show 3 clusters for each parent class and demonstrate that our learned features are able to cluster objects based on their size (\textit{Aeroplane, Bird, Cow}), context (\textit{Aeroplane, Bird, Person}), type (\textit{Boat, Bird}), pose (\textit{Cow}), and interaction with other categories (\textit{Person}).
	}
	\label{fig: cluster_viz}
\end{figure*}

\begin{figure*}[t]
	\centering
	\includegraphics[width=1\linewidth]{viz_weight.pdf}\\
	\caption{Visualizations of weights based on the
	t-SNE method that illustrates the relationships on semantic-level between parent classifier and the person sub-category classifier. We show that one person sub-category is usually close to one parent class, as they often co-appear in the same image, as shown in example images on two sides.}
	\label{fig: weight_viz}
\end{figure*}

\paragraph{Clustering Results.}
Since the ground truth labels are not available for sub-categories, we present visualizations of clustering results in Figure \ref{fig: cluster_viz} to measure the quality, in which each parent class shows 3 example clusters.
Our method is able to cluster objects based on their size (\textit{Aeroplane, Bird, Cow}), context (\textit{Aeroplane, Bird, Person}), type (\textit{Boat, Bird}), pose (\textit{Cow}), and interaction with other categories (\textit{Person}).
For instance, persons with different categories, e.g., horse, motobike, and boat, are clustered into different groups.
This visually validates that our learned feature representations are enhanced via the sub-category objective in an unsupervised manner.
More visual comparisons are presented in the supplementary material.

\paragraph{Weight Visualization.}
In order to understand how our learning mechanism improves the clustering quality, we visualize the distribution of the classifier weights, i.e., $\theta_p$ and $\theta_s$, via t-SNE \cite{maaten2008visualizing}.
As such, we are able to find the relationship between the parent classifier $H_p$ and the sub-category module $H_s$.
Figure \ref{fig: weight_viz} shows the visualization of weights, in which we take the sub-categories of person (denoted as yellow cross symbols) as the example, since the person category has more interactions with other parent classes (denoted as solid circles).
It illustrates that one person sub-category is often close to one parent class, e.g., sub-category \textit{person} and parent class \textit{bike}, which makes sense as those two categories usually co-appear in the same image (see example images in Figure \ref{fig: weight_viz} on two sides).

\begin{figure*}[!t]
	\centering
	\includegraphics[width=1\linewidth]{segmap_viz.pdf} \\
	\vspace{-2mm}
	\caption{Qualitative results on the PASCAL VOC 2012 validation set.
	(a) Input images. (b) Ground truth. (c) Our results.
	}
	\label{fig: segmap_viz}

\end{figure*}

\vspace{2mm}
\begin{table*} [t]
	\caption{Semantic segmentation performance on the PASCAL VOC 2012 validation set. Bottom group contains results with CRF refinement, while the top group is without CRF. Note that 11/20 classes obtain improvements using our approach w/ CRF. The best three results are in \textcolor{red}{red}, \textcolor{PineGreen}{green} and \textcolor{blue}{blue}, respectively.}

	\vspace{-1mm}
	\footnotesize
	\centering
	\renewcommand{\arraystretch}{1.3}
	\setlength{\tabcolsep}{2pt}
	\begin{tabular}{lccccccccccccccccccccc|c}
		\toprule

		Method & bkg & aero & bike & bird & boat & bottle & bus & car & cat & chair & cow & table & dog & horse & motor & person & plant & sheep & sofa & train & tv & mIoU \\
		\midrule

		AffinityNet  \cite{ahn2018learning} & 88.2 & 68.2 & \textcolor{blue}{30.6} & \textcolor{PineGreen}{81.1} & 49.6 & 61.0 & 77.8 & 66.1 & 75.1 & \textcolor{blue}{29.0} & 66.0 & 40.2 & \textcolor{PineGreen}{80.4} & 62.0 & \textcolor{PineGreen}{70.4} & \textcolor{red}{73.7} & \textcolor{blue}{42.5} & 70.7 & \textcolor{PineGreen}{42.6} & \textcolor{blue}{68.1} & 51.6 & 61.7 \\

		Ours (w/o CRF) & 88.1 & 49.6 & 30.0 & 79.8 & \textcolor{blue}{51.9} & \textcolor{PineGreen}{74.6} & \textcolor{PineGreen}{87.7} & 73.7 & \textcolor{PineGreen}{85.1} & \textcolor{PineGreen}{31.0} & \textcolor{blue}{77.6} & \textcolor{PineGreen}{53.2} & \textcolor{blue}{80.3} & \textcolor{PineGreen}{76.3} & \textcolor{blue}{69.6} & 69.7 & 40.7 & 75.7 & \textcolor{PineGreen}{42.6} & 66.1 & \textcolor{blue}{58.2} & \textcolor{blue}{64.8} \\

		\midrule

		MCOF \cite{wang2018weakly}  & 87.0 &  \textcolor{red}{78.4} & 29.4 & 68.0 & 44.0 & 67.3 & 80.3 & \textcolor{blue}{74.1} & 82.2 & 21.1 & 70.7 & 28.2 & 73.2 & 71.5 & 67.2 & 53.0 & \textcolor{PineGreen}{47.7} & 74.5 & 32.4 & \textcolor{PineGreen}{71.0} & 45.8 & 60.3 \\

		Zeng et al. \cite{zeng2019joint} &  \textcolor{red}{90.0} & \textcolor{PineGreen}{77.4} &  \textcolor{red}{37.5} & \textcolor{blue}{80.7} &  \textcolor{red}{61.6} & 67.9 & 81.8 & 69.0 & \textcolor{blue}{83.7} & 13.6 & \textcolor{PineGreen}{79.4} & 23.3 & 78.0 & \textcolor{blue}{75.3} &  \textcolor{red}{71.4} & 68.1 & 35.2 & \textcolor{blue}{78.2} & 32.5 & \textcolor{red}{75.5} & 48.0 & 63.3\\

		FickleNet \cite{lee2019ficklenet} &
		\textcolor{PineGreen}{89.5} & \textcolor{blue}{76.6} & \textcolor{PineGreen}{32.6} & 74.6 & 51.5 & \textcolor{blue}{71.1} & \textcolor{blue}{83.4} & \textcolor{PineGreen}{74.4} & 83.6 & 24.1 & 73.4 & \textcolor{blue}{47.4} & 78.2 & 74.0 & 68.8 &  \textcolor{PineGreen}{73.2} &  \textcolor{red}{47.8} &  \textcolor{red}{79.9} & \textcolor{blue}{37.0} & 57.3 & \ \textcolor{red}{64.6} & \textcolor{PineGreen}{64.9} \\

		Ours (w/ CRF) & \textcolor{blue}{88.8} & 51.6 & 30.3 & \textcolor{red}{82.9} & \textcolor{PineGreen}{53.0} & \textcolor{red}{75.8} &  \textcolor{red}{88.6} &  \textcolor{red}{74.8} &  \textcolor{red}{86.6} &  \textcolor{red}{32.4} &  \textcolor{red}{79.9} &  \textcolor{red}{53.8} &  \textcolor{red}{82.3} &
		\textcolor{red}{78.5} & \textcolor{PineGreen}{70.4} & \textcolor{blue}{71.2} & 40.2 & \textcolor{PineGreen}{78.3} & \textcolor{red}{42.9} & 66.8 & \textcolor{PineGreen}{58.8} &  \textcolor{red}{66.1} \\

		\bottomrule
	\label{table: compare_sota_detail}
	\end{tabular}
	\vspace{-5mm}
\end{table*}

\begin{table}[!t]
\caption{Comparison of weakly-supervised semantic segmentation methods on the PASCAL VOC 2012 val and test sets. In addition, we present methods that aim to improve the initial response with $\checkmark$ in the ``Init. Res.'' column.}
\label{tab:orbital_data}
\centering
\small
\renewcommand{\arraystretch}{1.1}
\setlength{\tabcolsep}{3.5pt}
\begin{tabular}{lcccc}
\toprule
\text{Method} & {Backbone} & {Init. Res.} & {Val} & {Test} \\
\midrule
    MCOF \textsubscript{CVPR'18} \cite{wang2018weakly} & ResNet-101 &  & 60.3 & 61.2 \\

    DCSP \textsubscript{BMVC'17} \cite{chaudhry2017discovering} & ResNet-101 &  & 60.8 & 61.9 \\

    DSRG \textsubscript{CVPR'18} \cite{huang2018weakly} & ResNet-101 &  & 61.4 & 63.2 \\

    AffinityNet \textsubscript{CVPR'18} \cite{ahn2018learning} & Wide ResNet-38 & & 61.7 & 63.7 \\

    SeeNet \textsubscript{NIPS'18} \cite{hou2018self} & ResNet-101 & \checkmark & 63.1 & 62.8 \\

    Zeng \textit{et al} \textsubscript{ICCV'19} \cite{zeng2019joint} &  DenseNet-169 &  & 63.3 & 64.3 \\

    BDSSW \textsubscript{ECCV'18} \cite{fan2018associating} & ResNet-101 & & 63.6 & 64.5 \\

    OAA \textsubscript{ICCV'19} \cite{jiang2019integral} & ResNet-101 & \checkmark & 63.9 & 65.6 \\

    CIAN \textsubscript{CVPR'19} \cite{fan2018cian} & ResNet-101 &  & 64.1 & 64.7 \\

    FickleNet \textsubscript{CVPR'19} \cite{lee2019ficklenet} & ResNet-101 & \checkmark & 64.9 & 65.3 \\

    Ours & ResNet101 & \checkmark & \textbf{66.1} & \textbf{65.9} \\

\bottomrule
\label{table: compare_sota}
\end{tabular}
\vspace{-8mm}
\end{table}

 \vspace{-2mm}
\subsection{Semantic Segmentation Performance}
After generating the pseudo ground truths as the results in Table \ref{table: compare_cam} (i.e., CAM + RW), we use them to train the semantic segmentation network.
We first compare our method with recent work using the ResNet-101 backbone or other similarly powerful ones in Table \ref{table: compare_sota}.
On both validation and testing sets, the proposed algorithm performs favorably against the  state-of-the-art approaches.
We also note that, most methods focus on improving the refinement stage or network training, while ours improves the initial step to generate better object response maps.

In Table \ref{table: compare_sota_detail}, we show detailed results for each category on the validation set.
We compare two groups of results with (bottom) or without (top) applying the CRF \cite{crf} refinement to the final segmentation outputs.
Compared to the recent FickleNet \cite{lee2019ficklenet} method that also focuses on improving the initial response map, the proposed algorithm performs favorably for the segmentation task in terms of the mean IoU.
We also note that, our results without applying CRF (mIoU as $64.8\%$) already achieves similar performance compared with the FickleNet (mIoU as $64.9\%$).
In Figure \ref{fig: segmap_viz}, we present some examples of the final semantic segmentation results, and show that our results are close to the ground truth segmentation.

\section{Conclusions}
In this paper, we propose a simple yet effective approach to improve the class activation maps by introducing a self-supervised task to discover sub-categories in an unsupervised manner. Without bells and whistles, our approach performs favorably against existing weakly-supervised semantic segmentation methods.
Specifically, we develop an iterative learning scheme by running clustering on image features for each parent class and train the classification network on sub-category objectives.
Unlike other existing schemes that aggregate multiple response maps, our approach generates better initial predictions without introducing extra complexity or inference time to the model.
We conduct extensive experimental analysis to demonstrate the effectiveness of our approach via exploiting the sub-category information.
Finally, we show that our algorithm produces better activation maps, thereby improving the final semantic segmentation performance.

\vspace{-3mm}
{\flushleft \bf{Acknowledgments.}}
This work is supported in part by the NSF CAREER Grant \#1149783, and gifts from eBay and Google.

{\small
\bibliographystyle{ieee_fullname}