\section{Introduction}
\IEEEPARstart{A}{} rare-earth barium copper oxide (ReBCO) YBa$_2$Cu$_3$O$_{7-\delta}$ coated conductor on round core (CORC) cable was proposed a decade ago \cite{laan2009}
as a new concept for practical superconducting (SC) coated conductor cables.
This has since been developed into a much thinner CORC {\it wire} with a round core diameter of a few millimeters \cite{laan2016,weiss2017}
by using an alloy substrate with a thickness of as little as 30~$\mu$m.
Flexible deformation along the former  has ever been demonstrated, although such deformation degrades the critical current to some extent \cite{weiss2017,souc2010,laan2011}.
Furthermore, the current technology for fabricating CORC wires achieves the shortest conductor pitch (less than 10~mm \cite{weiss2017}) among all high-temperature superconductor (HTS) cable concepts.
Therefore, CORC wires show promise as practical SC wires for applications such as the high-field magnets used in particle accelerators and magnetic resonance imaging (MRI) machines.

To date, there have been relatively few reports on magnetization and transport alternating current (AC) loss measurements of CORC cables or wires without striation \cite{souc2010,majoros2014,terzioglu2017}.
A helically wound tape structure (i.e., a continuous inclination of the tape) has been reported to decrease the magnetization loss
because of the reduced average perpendicular field component penetrating the tape over a full pitch \cite{souc2010}.
However, other measurements have been reported showing higher AC loss in a helically wound tape in the fully penetrated state than in
a flat tape.
This fact is attributed to the anisotropy of the critical current density $J_{\rm c}$ with respect to the applied field \cite{majoros2014};
a helically wound tape has a higher $J_{\rm c}$ than that of a flat tape in a certain portion of the tape oriented parallel to the field.
This also results in a sizable AC loss reduction relative to the loss in a flat tape at low field strength \cite{majoros2014}.
For a single-layer CORC cable in the fully penetrated state, using more
tapes in a single layer decreases the AC loss proportionally \cite{majoros2014}.
Meanwhile, in the case of a multi-layered CORC cable in the fully penetrated state,
the AC loss increases in proportion to the number of layers in the cable \cite{majoros2014}.
An effect of a copper tube former on the magnetization and transport AC losses has been reported:
eddy current loss makes a large contribution to the total transport AC loss,
and the magnetic shielding reduces the magnetization AC loss at low field strength \cite{terzioglu2017}.

Regarding simulation studies,
numerical simulations of AC losses based on the thin-sheet approximation have been carried out, the
focus being the spiraled multi-layer coated conductors used for SC power transmission cables \cite{takeuchi2011,amemiya2013}.
Another approach to determining the electromagnetic properties of a power transmission cable composed of helically wound SC tapes is based on electric-circuit theory using the Bean model \cite{tominaka2009}.
A computational method was developed by transforming three-dimensional helicoidal structures into two-dimensional ones \cite{stenvall2013}.
To simulate a CORC cable, the
AC losses were simulated at first by a two-dimensional model neglecting the winding of the tapes \cite{solovyov2014}.
Later, an improved simulation study of the magnetization losses in CORC cables was carried out using three-dimensional models based on the H formulation \cite{terzioglu2017,sheng2017}.
Simulations with three-dimensional models using the T-A formulation have also been conducted to investigate the AC losses in CORC cables, twisted stacked-tape conductor cables, and double coaxial cables with various transport currents and background magnetic fields \cite{fu2018}. Comprehensive
simulations of the magnetization losses of CORC cables wound from the SC tapes with and without striation have been reported recently based on
a three-dimensional model using the T-A formulation \cite{wang2019}.

However, an analytical approach to modeling the AC losses of a CORC wire is currently lacking, partly because of the rather complex wire structure,
and therefore the theoretical limit of the AC loss reduction in a CORC wire remains elusive.
Herein, we develop an analytical method for modeling the magnetization loss of a helically wound SC tape. The
magnetization loss is simulated to obtain properties such as its dependences on $R$ and $L_{\rm p}$ (see Fig.~\ref{fig1}), and
the simulation results are verified by the analytical formula in the limiting case of $w_0/R \rightarrow 2\pi$ or $0$.
We consider the situation in which ReBCO tapes wound round an MRI magnet are exposed to a magnetic field that is ramped up to high strength \cite{yokoyama2017,yachida2017}.
The fully penetrated state is realized, resulting in the steady state in which the magnetic field profile on the tape is unchanged.
Although the effects of tape stacking and $J_{\rm c}$ anisotropy on the AC losses are important issue \cite{majoros2014} for the quantitative evaluation,
we ignore such effects and focus instead on the simplest model.
The presented models and results represent a first step toward further advanced research and development of a high-field magnet wound from CORC wires.

\section{Model of a helically wound superconducting tape on a hollow cylinder}
\begin{figure}[tb]
\centering
\includegraphics[width=3.0in]{fig1.eps}
\caption{Schematic of a helically wound superconducting (SC) tape of width $w_0$ and pitch $L_{\rm p}$ on a hollow cylinder of radius $R$.
The axes $(\xi,\zeta)$ are on the tape surface.}
\label{fig1}
\end{figure}

We consider a helically wound SC tape on a hollow cylinder as shown in Fig.~\ref{fig1}.
The tape width and thickness and the radius of the hollow cylinder are $w_0$, $d_0$, and $R$, respectively.
The helically wound SC tape in the present simple model is assumed to be sufficiently thin ($d_0 \ll w_0$) that
it can be approximated as a helical SC tape surface with infinitesimal thickness as depicted in Fig.~\ref{fig1}.
The SC tape surface is described by the coordinates (see Appendix \ref{appendix_reduced-faraday})
\begin{equation}
\label{helicoid-coordinate}
\setlength{\nulldelimiterspace}{0pt}
\left\{
\begin{IEEEeqnarraybox}[\relax][c]{l's}
x=R \cos[ (k\zeta-\xi/R )/\sqrt{1+(kR)^2} ],\\
y=R \sin[(k\zeta-\xi/R)/\sqrt{1+(kR)^2}],\\
z=(\zeta+kR\xi)/\sqrt{1+(kR)^2},
\end{IEEEeqnarraybox}
\right.
\end{equation}
where $k=2\pi/L_{\rm p}$ and $L_{\rm p}$ is the pitch length of the helical conductor.
The helically wound SC tape surface with a full pitch $L_{\rm p}$ corresponds to $-w_0/2\le\xi\le w_0/2$ and $0\le\zeta \le L_{\rm tape}$,
with $L_{\rm tape}=L_{\rm p} \sqrt{1+(k R)^2}$ being the tape length corresponding to $L_{\rm p}$.
The $\xi$ and $\zeta$ axes are orthogonal to each other on the helical SC tape surface as shown in Fig.~\ref{fig1}.

\section{Response to an applied field on a helical superconducting tape}
\subsection{Reduced Faraday--Maxwell equation}
\begin{figure}[tb]
\centering
\includegraphics[width=1.3in]{fig2.eps}
\caption{Current streamlines and profile of the scalar function $g(\xi,\zeta)$ on a helically wound SC tape surface with $w_0=2$~mm and $L_{\rm p}=12$~mm
for $\beta=3$~mT/s and $R=1.5$~mm.
The tape length for one pitch of the helical conductor is $L_{\rm tape} \approx 15.26$~mm.
}
\label{fig2}
\end{figure}

We assume that the superconductor is in the steady state, which is plausible at high applied fields stronger than the full-flux-penetration field.
We therefore neglect the contribution from the transport current because the magnetization losses are dominant over the transport losses at high field strength \cite{kajikawa2016}. We also ignore
the magnetic field due to the screening current because that field is much weaker than the applied field.
In an applied field $\bm{B}_{\rm a}(t)=B_{\rm a}(t)\hat{\bm{x}}$ that is ramped constantly at a rate $\beta={\rm d}B_{\rm a}(t)/{\rm d}t$,
the Faraday--Maxwell equation in the steady state is
\begin{equation}
\bm{\nabla}\times\bm{E}=-\frac{\partial \bm{B}}{\partial t} \approx -\frac{{\rm d} \bm{B}_{\rm a}(t)}{{\rm d} t} = -\beta \hat{\bm{x}},
\label{Faraday-maxwell-equation-steady}
\end{equation}
where $\bm{E}$ (resp.\ $\bm{B}$) is the electric (resp.\ magnetic) field and
$\hat{\bm{x}}$ is the unit vector in the $x$ direction.
With the thin-sheet approximation \cite{zhang2017},
we take account of only the response to the applied field component perpendicular to the tape surface.
Consequently, the Faraday--Maxwell equation reduces to the following two-dimensional equation on the tape surface \cite{higashi2019,higashi2018} (see Appendix \ref{appendix_reduced-faraday}):
\begin{equation}
\frac{\partial}{\partial \xi}\left(\rho_{\rm sc}\frac{\partial g}{\partial \xi} \right)+ \frac{\partial}{\partial \zeta}\left( \rho_{\rm sc}\frac{\partial g}{\partial \zeta} \right)=\beta \cos\left( \frac{k\zeta-\xi/R}{\sqrt{1+(kR)^2}}\right),
\label{reduced-FM-equation}
\end{equation}
where $g(\xi,\zeta)$ a scalar function defined on the SC tape surface.
The contour lines of $g(\xi,\zeta)$ describe the current streamlines on the tape surface, as depicted by the solid lines in Fig.~\ref{fig2}.
The SC property is incorporated through the electric field ($\bm{E}$)--current density ($\bm{J}$) characteristics, namely $\bm{E}=\rho_{\rm sc}\bm{J}$, with
\begin{equation}
\rho_{\rm sc}(|\bm{J}|)=\frac{E_{\rm c}}{J_{\rm c}}\left( \frac{|\bm{J}|}{J_{\rm c}} \right)^{n-1}
\end{equation}
being the isotropic SC nonlinear resistivity. Herein,
the electric field criterion is fixed to $E_{\rm c}=1$~$\mu$V/cm.

\subsection{Boundary conditions}
The specifications of the helically wound SC tape are summarized in Table~\ref{table1}.
We simulate the current streamlines and the magnetization losses by solving Eq.~(\ref{reduced-FM-equation}) numerically using the commercial software COMSOL Multiphysics\textsuperscript{\textregistered} \cite{comsol}.
The Dirichlet boundary condition $g(\xi=\pm w_0/2,\zeta)=0$ is imposed on $g(\xi,\zeta)$ along the long edges of the SC tape,
whereas the periodic boundary condition $g(\xi,\zeta)=g(\xi,\zeta+L_{\rm tape})$ is imposed at the terminals of the SC tape with a full pitch. Herein,
the field sweep rate is set to $\beta=3$~mT/s because the typical field sweep rate of an MRI magnet is of the order of milliteslas per second \cite{yokoyama2017,yachida2017}.

\section{Results and discussion}

\begin{table}
\caption{Specifications of helically wound SC tape}
\label{table}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{|p{25pt}|p{106pt}|p{100pt}|}
\hline
Symbol&
Quantity&
Dimensions of helically wound SC tape \\
\hline
$w_0 $& tape width& 2~mm \\
$d_0 $&  tape thickness &2 $\mu$m\\
$L_{\rm p}$ & pitch length of helical conductor& 2--200~mm \\
$R$ & radius of hollow cylinder & 1/$\pi$--20~mm\\
\hline
&&Superconducting characteristics\\
\hline
$J_{\rm c}$& critical current density & 5$\times$10$^{10}$~A/m$^2$\\
$n$ & $n$-value & 29 \\
\hline
\end{tabular}
\label{table1}
\end{table}

\subsection{Current streamlines}
By solving Eq.~(\ref{reduced-FM-equation}) numerically,
the current density vector on the tape surface can be calculated via
\begin{equation}
\bm{J}(\xi,\zeta)=J_\xi \hat{\bm{\xi}}+J_\zeta\hat{\bm{\zeta}},~~  J_\xi=-\frac{\partial g}{\partial \zeta},~~J_\zeta=\frac{\partial g}{\partial \xi},
\label{current-density}
\end{equation}
where $\hat{\bm{\xi}}$ (resp.\ $\hat{\bm{\zeta}}$) is the unit vector in the $\xi$ (resp.\ $\zeta$) direction.
The color scale in Fig.~\ref{fig2} shows $g(\xi,\zeta)/(J_{\rm c}w_0/2)$. The solid lines are the current streamlines, which
are divided into two parts within a full pitch because of the helical winding,
thereby reducing the effective tape length to $L_{\rm tape}/2$.

\subsection{Magnetization loss}
\begin{figure}[tb]
\centering
\includegraphics[width=3.5in]{fig3.eps}
\caption{Dependence of loss power per unit tape length on cylinder radius (plotted double logarithmically).
The circles indicate where the SC tape overlaps.
}
\label{fig3}
\end{figure}

\begin{figure}[tb]
\centering
\includegraphics[width=3.5in]{fig4.eps}
\caption{Dependence of loss power per unit tape length on the pitch length of the helical conductor (plotted double logarithmically).
The circles indicate where the SC tape overlaps.
}
\label{fig4}
\end{figure}

\begin{figure}[tb]
\centering
\includegraphics[width=3.5in]{current-density-profile_corc_Lp10w0_R1w0_w0-2mm.eps}
\caption{Current density profile on a helical tape in the two and three dimensional models.
The pitch length and the radius of the hollow cylinder are set to $L_{\rm p}=10 w_0$ and $R=w_0$, respectively.
The other parameters are common.
}
\label{fig:current}
\end{figure}

The loss power $P$ per unit tape length is evaluated by (see Appendix \ref{appendix_coordinate-transformation})
\begin{equation}
P=\frac{d_0}{L_{\rm tape}}\int {\rm d}\xi\int{\rm d}\zeta \rho_{\rm sc}(\xi,\zeta)\left[ \left( \frac{\partial g}{\partial \xi}\right)^2+\left( \frac{\partial g}{\partial \zeta}\right)^2\right].
\end{equation}

Figure~\ref{fig3} shows the double logarithmic plot of the loss power per unit tape length versus $R$ for different values of $L_{\rm p}$.
At $R/w_0=1/2\pi$, the helical SC tape overlaps with itself irrespective of $L_{\rm p}$,
and for $L_{\rm p}/w_0=1$, the SC tape overlaps irrespective of $R$.
The circles in Fig.~\ref{fig3} indicate
where the tape overlaps.

For $R\gtrsim w_0$,
the loss power saturates close to the analytically evaluated value (the dashed line in Fig.~\ref{fig3}) in the thin-filament limit of $w_0/R\rightarrow 0$ via Eq.~(\ref{loss_thin-filament}).
For $n=29$, Eq.~(\ref{loss_thin-filament}) yields $P(w_0/R\rightarrow 0)\approx0.165$~mW/m,
which is close to $(2/\pi)P_{\rm flat}$ (where $P_{\rm flat}$ is the loss power of the flat tape) in the Bean limit of $n\rightarrow \infty$
as evidenced by a recent magnetization measurement \cite{myers2019} [see also Eq.~(\ref{loss_thin-filament})].

Meanwhile, for $R\lesssim w_0$, with decreasing $R$ the loss power becomes much
lower than that in the limit of $w_0/R\rightarrow 0$ (the dashed line in Fig.~\ref{fig3}).
Upon decreasing $R/w_0$, the helically wound tape behaves more like a tubular wire, as discussed in Ref.~\cite{majoros2014}.
Indeed, the loss power for $L_{\rm p}/w_0=100$ at $R/w_0=1/2\pi$ shows quantitatively good agreement with the analytically evaluated loss power (the dashed-dotted line in Fig.~\ref{fig3})
for a tubular wire via Eq.~(\ref{loss_cylinder}),
although there is some deviation from the theoretical value for a tubular wire for $L_{\rm p}/w_0=1$ and 2 because of the short pitch.
The slope of the dashed-dotted line is two in the Bean limit
because the loss power for a tubular superconductor increases with $R^2$ as in Eq.~(\ref{loss_flat-tape_tube}).
The initial $R$ slope for each value of $L_{\rm p}/w_0$ is less than two, although it approaches two with increasing $L_{\rm p}/w_0$.
We ascribe the initial slope being less than two to the finite pitch effect specific to the helically wound tape.

Figure~\ref{fig4} shows the dependence of the loss power per unit tape length on $L_{\rm p}$ in double logarithm form for different values of $R$.
The loss power shows weak $L_{\rm p}$ dependence for $L_{\rm p}/w_0 \lesssim 10$.
At $R/w_0=4$, the loss power coincides with that evaluated via Eq.~(\ref{loss_thin-filament}) in the thin-filament limit (the dashed line in Fig.~\ref{fig4}) independent of $L_{\rm p}$.
Meanwhile, for $R/w_0 <4$,
the loss power becomes lower than the theoretical value in the thin-filament limit
and becomes even lower with decreasing $R$ as already mentioned in Fig.~\ref{fig3}.
The circles in Fig.~\ref{fig4} indicate where the tape overlaps.

Figure~\ref{fig:current} shows the profiles of the current density $J_{\zeta}(v,\zeta=0.5L_{\rm p})$ in the direction of the long dimension of a helical tape in two and three dimensional models.
The profile of $J_{\zeta}(v,\zeta=0.5L_{\rm p})$ on a helical tape shows a quantitative agreement between two and three dimensional models.
Here, the axes $v$ and $\zeta$ are defined via $(x,y,z)=(R\cos k(\zeta-v),R\sin k(\zeta-v),\zeta)$.

Next, we explain
physically why the loss power decreases for $R \lesssim w_0$.
The length scale for the penetration of magnetic flux vortices decreases as $R$ decreases, thereby
reducing the AC losses \cite{mawatari2011}.
In the present study, we consider the steady state wherein the applied field should far exceed the full-flux-penetration field.
We therefore consider that it is
the reduced flux-penetration length scale, rather than the screening of the applied magnetic field, that reduces the loss power.

Next, we evaluate approximately the theoretically expected maximum reduction rate of the loss power upon reducing $R$.
We evaluate the ratio of the loss power for $R/w_0=1/2\pi$ to that for $R/w_0\gg 1$ (i.e., $w_0/R\rightarrow 0$) by using Eqs.~(\ref{loss_thin-filament}) and (\ref{loss_cylinder}):
\begin{equation}
\frac{P(w_0/R\rightarrow 2\pi)}{P(w_0/R\rightarrow 0)}=\frac{2+1/n}{\pi^{1+1/n}}\rightarrow \frac{2}{\pi}~~(n\rightarrow \infty).
\label{reduction-rate}
\end{equation}
This suggests that by thinning the wire radius down to $R\ll w_0$, the loss power of the helical tape becomes smaller than $P_{\rm flat}$ by a factor of $(2/\pi)^2\approx 0.405$
in accordance with a previous study \cite{mawatari2011}. Our numerical results
confirm the reduction rate of the loss power.
At $L_{\rm p}/w_0=100$, the ratio of the loss power for $R/w_0=1/2\pi$ to that for $R/w_0=4$ is
$P(R/w_0=1/2\pi)/P(R/w_0=4)\approx 0.619$, which agrees well with $(2+1/n)/\pi^{1+1/n}\approx 0.623$ [Eq.~(\ref{reduction-rate})] for $n=29$.

\begin{figure}[tb]
\centering
\includegraphics[width=3.1in]{fig5.eps}
\caption{Critical line on which the helically wound SC tape begins to overlap.
}
\label{fig5}
\end{figure}

Next, we consider the condition that must be met for the helically wound SC tape not to overlap.
This is determined by $L_{\rm p}>w_0/\cos \alpha$ (see Appendix \ref{appendix_reduced-faraday}.1).
By solving the condition with respect to $L_{\rm p}$, we obtain the condition that $L_{\rm p}>w_0/\sqrt{1-(w_0/2\pi R)^2}$ (see Fig.~\ref{fig5}).

\subsection{CORC wire with multiple tapes}
We discuss the effect of multiple tapes on magnetization losses in a single layer CORC wire.
The response to an applied field is sufficiently large at high fields far above a full flux penetration field $B_{\rm p}$ that the interaction between SC tapes is relatively small and negligible,
although the magnetic field shielding by neighboring tapes should be important to model a CORC wire wound by several tapes at low fields \cite{majoros2014,wang2019}.
The magnetic field due to the shielding current is of the order of $\mu_0 J_{\rm c}d_0$ ($\approx 0.126$ T for the parameters adopted in this paper).
It is sufficiently small comparing to the background field.
Thus, we can ignore the contribution from a shielding current and treat total losses as a simple sum of the magnetization loss on an individual helical tape.

Our assumption should be plausible when the background field strength is as high as a few Tesla or at least above $B_{\rm p}$.
$B_{\rm p}$ for a single helical tape is approximately evaluated by the formula for a flat tape, $B_{\rm p}=(\mu_0 J_{\rm c}d_0/\pi)[1 + \ln(w_0/d_0)] \approx 0.316$ T \cite{brandt1996},
where $\mu_0$ is the magnetic permeability of a vacuum.

\section{Summary}
We investigated theoretically the magnetization loss of a helically wound SC tape on a hollow cylinder in the steady state under a ramped field.
For $R\gtrsim w_0$, the loss power saturates close to the value expected for a loosely twisted SC tape.
Meanwhile, for $R\lesssim w_0$, it becomes lower than that value.
Compared with the loss power $P_{\rm flat}$ for a flat SC tape,
the loss power for a helical SC tape is roughly that for a loosely twisted tape ($\sim 2/\pi$ times lower than $P_{\rm flat}$) for $R\gtrsim w_0$.
Further reduction of the loss power value is expected by thinning the round core diameter down to $R\lesssim w_0$.
For a thin wire with a radius of a few millimeters,
the loss power per unit tape length is almost independent of $L_{\rm p}$.

\section*{Acknowledgment}
We would like to thank Y.~Yoshida (AIST) for providing us with the properties and dimensions of commercial superconducting CORC\textsuperscript{\textregistered} wires.
We also thank R.~Toyomoto (Kyoto Univ.) and N.~Amemiya (Kyoto Univ.) for sharing unpublished AC-loss measurement and simulation data with us.
The present work is based on the results obtained from a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO).

\appendices
\renewcommand{\thesectiondis}[2]{\Alph{section}:}
\section{Derivation of reduced Faraday--Maxwell equation}
\label{appendix_reduced-faraday}
\renewcommand{\theequation}{A.\arabic{equation}}
\setcounter{equation}{0}

\renewcommand{\thesubsectiondis}[0]{\Alph{section}.\arabic{subsection}.}
\setcounter{subsection}{0}
\subsection{Model of helical tape}
\label{appendix_model-helical}
We use the following helical coordinates to model the helically wound SC tape:
\begin{equation}
\label{helix-coordinates}
\setlength{\nulldelimiterspace}{0pt}
\left\{
\begin{IEEEeqnarraybox}[\relax][c]{l's}
x=u \cos\theta , \\
y=u \sin\theta,\\
z=\theta/k+v,
\end{IEEEeqnarraybox}
\right.
\end{equation}
where $k=2\pi/L_{\rm p}$ and $L_{\rm p}$ is the pitch length of the helical conductor.
The tape surface for one pitch length of the helical conductor corresponds to $u=R$, $0\le v \le w_0/\cos\alpha$, and $0\le\theta\le2\pi$, where
$R$ is the radius of the hollow cylinder, $w_0$ is the tape width, and $\alpha$ is the lay angle defined by $\cos\alpha=kR/\sqrt{1+(kR)^2}$.

In general,
the current density vector is described in terms of the current vector potential $\bm{T}$ as $\bm{J}=\bm{\nabla}\times \bm{T}$ so as to satisfy $\bm{\nabla}\cdot\bm{J}=0$.
In the thin-sheet approximation, we regard the SC tape as a strip with infinitesimal thickness,
and therefore the electric current flows on the tape surface at $u=R$.
In this case, the current vector potential can be expressed in terms of the tape-normal unit vector $\bm{\nabla}u$ as $\bm{T}=\bm{\nabla}\times[g(\theta,v)\bm{\nabla}u]$.
Here, the scalar function $g(\theta,v)$ describes the current streamlines on the tape surface.
The current density on the tape surface at $u=R$ is calculated as
\begin{equation}
\bm{J}=\bm{\nabla}g(\theta,v)\times\bm{\nabla}u=J_v\hat{\bm{v}}+J_\theta\hat{\bm{\theta}},
\label{current-density-u-v-theta}
\end{equation}
where $\hat{\bm{v}}$ and $\hat{\bm{\theta}}$ are the unit vectors in the $v$ and $\theta$ directions, respectively, and
\begin{equation}
J_v=-\frac{1}{R}\frac{\partial g}{\partial \theta},~~~J_\theta=\frac{\sqrt{1+(kR)^2}}{kR}\frac{\partial g}{\partial v}.
\end{equation}

\subsection{Reduced Faraday--Maxwell equation on a helical tape}
The electric field in the helical coordinates (\ref{helix-coordinates}) is $\bm{E}=E_u\hat{\bm{u}}+E_v\hat{\bm{v}}+E_\theta\hat{\bm{\theta}}$ with
\begin{align}
\hat{\bm{u}}&=\bm{\nabla}u=\hat{\bm{x}}\cos\theta+\hat{\bm{y}}\sin\theta,\label{u-hat}\\
\hat{\bm{v}}&=\bm{\nabla}\theta/k+\bm{\nabla}v=\hat{\bm{z}},\label{v-hat}\\
\hat{\bm{\theta}}&=(\sqrt{1+(ku)^2}/k)\bm{\nabla}\theta+\bm{\nabla}v/\sqrt{1+(ku)^2}\nonumber \\
&=ku(-\hat{\bm{x}}\sin\theta+\hat{\bm{y}}\cos\theta+\hat{\bm{z}}/ku)/\sqrt{1+(ku)^2} \label{theta-hat}.
\end{align}
Then, the electric field is expressed as
$\bm{E}=E_u\bm{\nabla}u+\tilde{E}_v\bm{\nabla}v+\tilde{E}_\theta\bm{\nabla}\theta$,
where $\tilde{E}_v=E_v+E_\theta/\sqrt{1+(ku)^2}$
and $\tilde{E}_\theta =E_v/k+E_\theta\sqrt{1+(ku)^2}/k$.
Using the relations $\bm{\nabla}u\times\bm{\nabla}v=-(\sqrt{1+(ku)^2}/ku)\hat{\bm{\theta}}$, $\bm{\nabla}v\times\bm{\nabla}\theta=-\hat{\bm{u}}/u$, and
$\bm{\nabla}\theta \times \bm{\nabla}u=-\hat{\bm{v}}/u$,
we obtain
\begin{align}
\bm{\nabla}\times\bm{E}
&=-\left( \frac{\partial \tilde{E}_v}{\partial u}-\frac{\partial E_u}{\partial v} \right) \frac{\sqrt{1+(ku)^2}}{ku}\hat{\bm{\theta}} \nonumber\\
&-\left( \frac{\partial\tilde{E}_\theta}{\partial v}-\frac{\partial \tilde{E}_v}{\partial \theta}\right)\frac{\hat{\bm{u}}}{u}-\left(\frac{\partial E_u}{\partial \theta}-\frac{\partial \tilde{E}_\theta}{\partial u} \right)\frac{\hat{\bm{v}}}{u}.
\label{faraday-induction-polar}
\end{align}
Because we adopt the thin-sheet approximation, we may consider the response to the perpendicular field component alone.
Therefore, projecting Eq.~(\ref{Faraday-maxwell-equation-steady}) with eq.~(\ref{faraday-induction-polar}) onto the unit vector normal to the tape surface $\hat{\bm{u}}$ gives
\begin{align}
\frac{\partial}{\partial v}&\left( E_v/ku+\sqrt{1+(ku)^2}E_\theta/ku\right)\nonumber\\
&-\frac{1}{u}\frac{\partial}{\partial \theta}\left( E_v+E_\theta/\sqrt{1+(ku)^2}\right)=\beta \cos\theta
\label{Faraday-E-helical}
\end{align}
by using $\hat{\bm{u}}\cdot \hat{\bm{\theta}}=0$ and $\hat{\bm{u}}\cdot \hat{\bm{v}}=0$.
Combining Eq.~(\ref{Faraday-E-helical}) with $(E_\theta,E_v)=\rho_{\rm sc} (J_\theta,J_v)$, the Faraday--Maxwell equation on the tape surface $u=R$ is reduced to
\begin{align}
\frac{\partial }{\partial v}&\left[ \rho_{\rm sc} \frac{1+(kR)^2}{(kR)^2}\frac{\partial g}{\partial v}-\frac{\rho_{\rm sc}}{kR}\frac{1}{R}\frac{\partial g}{\partial \theta} \right]\nonumber \\
&+\frac{1}{R}\frac{\partial}{\partial \theta}\left[\frac{\rho_{\rm sc}}{R}\frac{\partial g}{\partial \theta}-\frac{\rho_{\rm sc}}{kR}\frac{\partial g}{\partial v}  \right]=\beta\cos\theta.
\label{reduced-Faraday}
\end{align}

Although the derivation of the reduced Faraday--Maxwell equation (\ref{reduced-Faraday}) is clear in the coordinates of Eq.~(\ref{helix-coordinates}),
Eq.~(\ref{reduced-Faraday}) includes the divergent terms in the loosely wound limit of $kR\rightarrow 0$ as is obvious from geometric considerations. We
therefore work in the coordinates in which the spatial axes $(\xi,\zeta)$ are orthogonal on the tape surface.
An arbitrary position on the tape surface in the $(\xi,\zeta)$ coordinates is expressed in terms of $(v,\theta)$ as
\begin{equation}
\xi=v  \cos\alpha,~~\zeta=\sqrt{(R\theta)^2+(\theta/k)^2}+v\sin\alpha,
\end{equation}
from which we have that
\begin{equation}
v=\frac{\sqrt{1+(kR)^2}}{kR}\xi,~~\theta=\frac{k\zeta-\xi/R}{\sqrt{1+(kR)^2}}.
\label{v-theta}
\end{equation}
Thus, on the tape surface $u=R$, from Eqs.~(\ref{helix-coordinates}) and (\ref{v-theta}), we obtain Eq.~(\ref{helicoid-coordinate}).

For $\bm{r}=x\hat{\bm{x}}+y\hat{\bm{y}}+z\hat{\bm{z}}$, we have that
\begin{align}
\hat{\bm{\xi}}&=\partial \bm{r}/\partial \xi=-\hat{\bm{\theta}}/ku+(\sqrt{1+(ku)^2}/ku)\hat{\bm{v}},\label{xi-hat}\\
\hat{\bm{\zeta}}&=\partial \bm{r}/\partial \zeta=\hat{\bm{\theta}}\label{zeta-hat}.
\end{align}
Thus, we obtain
\begin{eqnarray}
\hat{\bm{v}}=(\hat{\bm{\zeta}}+ku\hat{\bm{\xi}})/\sqrt{1+(ku)^2},~~\hat{\bm{\theta}}=\hat{\bm{\zeta}}.
\label{v-theta-hat}
\end{eqnarray}
As for the reduced Faraday--Maxwell equation, substituting the expressions
\begin{eqnarray}
\frac{1}{R}\frac{\partial g}{\partial \theta}=\frac{1}{R}\left( \frac{\partial g}{\partial \xi}\frac{\partial \xi}{\partial \theta}+\frac{\partial g}{\partial \zeta}\frac{\partial \zeta}{\partial \theta}\right)=\frac{\partial g}{\partial \zeta}\frac{1}{\cos\alpha},\\
\frac{\partial g}{\partial v}=\frac{\partial g}{\partial \xi}\frac{\partial \xi}{\partial v}+\frac{\partial g}{\partial \zeta}\frac{\partial \zeta}{\partial v}=\frac{\partial g}{\partial \xi}\cos\alpha+\frac{\partial g}{\partial \zeta}\sin\alpha
\label{g-derivative}
\end{eqnarray}
into Eq.~(\ref{reduced-Faraday}) yields Eq.~(\ref{reduced-FM-equation}).
Furthermore, by substituting Eqs.~(\ref{v-theta-hat})--(\ref{g-derivative}) into Eq.~(\ref{current-density-u-v-theta}),
we have the screening current density Eq.~(\ref{current-density}) on the tape surface $u=R$.

\section{Vectors in $(\xi,u,\zeta)$ coordinates}
\renewcommand{\theequation}{B.\arabic{equation}}
\setcounter{equation}{0}
\label{appendix_coordinate-transformation}
First, by using Eq.~(\ref{v-theta-hat}),
the coordinate transformation from the $(\theta,u,v)$ coordinates to the $(\xi,u,\zeta)$ ones can be achieved for an arbitrary vector $\bm{O}=(O_x,O_y,O_z)$.
Second, in the same way, $\bm{O}$ in the $(\xi,u,\zeta)$ coordinates is transformed to that in the $(x,y,z)$ ones by using Eqs.~(\ref{u-hat})--(\ref{theta-hat}), (\ref{xi-hat}), and (\ref{zeta-hat}):
$
\bm{O}=O_u\hat{\bm{u}}+O_v\hat{\bm{v}}+O_\theta\hat{\bm{\theta}}=O_u\hat{\bm{u}}+O_\xi\hat{\bm{\xi}}+O_\zeta\hat{\bm{\zeta}}=O_x\hat{\bm{x}}+O_y\hat{\bm{y}}+O_z\hat{\bm{z}},
$
where $O_\xi\equiv kuO_v/\sqrt{1+(ku)^2}$ and $O_\zeta\equiv O_\theta+O_v/\sqrt{1+(ku)^2}$.
In the case of the thin-sheet approximation, we have on the tape surface $u=R$ that
\begin{equation}
\label{vector-thin-film-approx}
\setlength{\nulldelimiterspace}{0pt}
\left\{
\begin{IEEEeqnarraybox}[\relax][c]{l's}
O_x\approx(O_\xi \sin\alpha-O_\zeta \cos\alpha)\sin\theta, \\
O_y\approx-(O_\xi \sin \alpha-O_\zeta \cos \alpha)\cos\theta, \\
O_z\approx (O_\xi \cos\alpha +O_\zeta \sin\alpha).
\end{IEEEeqnarraybox}
\right.
\end{equation}
Herein, note that $O_u$ vanishes in the thin-sheet approximation when $\bm{O}$ refers to the current density or the electric field, neither of which has an out-of-plane component.
Thus we have from Eq.~(\ref{vector-thin-film-approx}) that $O^2_x+O^2_y+O^2_z=O^2_\xi+O^2_\zeta$.
Note also that the Jacobian $\partial(x,y,z)/\partial(\xi,u,\zeta)=1$.

\section{Magnetization loss in thin-filament or cylinder limit}
\renewcommand{\theequation}{C.\arabic{equation}}
\setcounter{equation}{0}
\label{appendix_gapped-cylinder}
\subsection{Thin-filament limit of $w_0/R \rightarrow 0$}
\label{appendix_thin-filament-limit}
\begin{figure}[tb]
\centering
\includegraphics[width=3.5in]{fig6.eps}
\caption{Schematic of a helical tape in the thin-filament limit of $w_0/R \rightarrow 0$.}
\label{fig6}
\end{figure}
In the thin-filament limit, we may take the limit of $\xi/R\rightarrow 0$ while keeping $kR$ finite (see Fig.~\ref{fig6}).
In this limit, the scalar function $g(\xi,\zeta)$ is almost uniform in the $\zeta$ direction,
and therefore the second term on the left-hand side (l.h.s.) of Eq.~(\ref{reduced-FM-equation}) can be neglected because $\partial g/\partial \zeta \approx 0$.
Thus, we obtain the equation for the scalar function:
\begin{equation}
\frac{\partial}{\partial \xi}\left( \rho_{\rm sc}\frac{\partial g}{\partial \xi} \right) \approx \beta\cos\phi(\zeta),
\label{eq_thin-filament}
\end{equation}
with $\phi(\zeta)=k\zeta/\sqrt{1+(kR)^2}$.
Equation~(\ref{eq_thin-filament}) has been solved previously with the Dirichlet boundary condition $g(\xi=\pm w_0/2,\zeta)=0$ \cite{higashi2019}.
The loss power per unit tape length coincides with that for a twisted tape \cite{higashi2019} and is readily evaluated as
\begin{align}
P\left({w_0\over R}\rightarrow 0\right)
&=\frac{B(\frac{2n+1}{2n},\frac{1}{2})}{\pi}\left( \frac{\beta w_0}{2E_{\rm c}} \right)^{1\over n}\frac{J_{\rm c} d_0 w^2_0 \beta}{2(2+1/n)}
\label{loss_thin-filament}\\
&\to (2/\pi)J_{\rm c}d_0w^2_0 \beta/4,~~~(n\to \infty), \nonumber
\end{align}
where $B(p,q)=B(q,p)=2\int_0^{\pi/2}{\rm d}\theta \cos^{2p-1} \theta \sin^{2q-1} \theta$ is the beta function and $p$ and $q$ are positive real numbers.

\subsection{Cylinder limit of $w_0/R \rightarrow 2\pi$}
\begin{figure}[tb]
\centering
\includegraphics[width=3.5in]{fig7.eps}
\caption{Schematic of a helical tape in the limit of $w_0/R \rightarrow 2\pi$. By taking the limit of $kR \rightarrow 0$ as well, a helical tape approaches to a gapped cylindrical superconductor.}
\label{fig7}
\end{figure}
In the loosely wound limit of $L_{\rm p} \rightarrow \infty$ (i.e., $kR\rightarrow 0$),
$g(\xi,\zeta)$ can be viewed as uniform in the $\zeta$ direction,
and thus the second term on the l.h.s.\ of Eq.~(\ref{reduced-FM-equation}) may be neglected.
We note that $\theta=(k\zeta-\xi/R)/\sqrt{1+(kR)^2}$ on the right-hand side of Eq.~(\ref{reduced-FM-equation}) should be kept over a full pitch to take account of the helical winding effect of the tape.
However, we encounter severe difficulty in evaluating the loss power analytically in this case.
Instead, we take the cylinder limit of $w_0/R\rightarrow 2\pi$ in addition to the limit of $kR \rightarrow 0$ to understand the dependence of the loss power on $R$ (see Fig.~\ref{fig7}).
By taking the limit of $kR\rightarrow 0$, we obtain
\begin{equation}
\frac{\partial}{ \partial \xi} \left( \rho_{\rm sc} \frac{\partial g}{\partial \xi}\right) \approx \beta \cos(-\xi/R).
\label{eq_gapped-cylinder}
\end{equation}
In this case, we can neglect the $kR\zeta$ term because it is not necessary to keep $\theta=(k\zeta-\xi/R)/\sqrt{1+(kR)^2}$ over a full pitch $L_{\rm p}$ in the limit of $w_0/R \rightarrow 2\pi$.
Equation~(\ref{eq_gapped-cylinder}) lacks the information about the helical structure of the tape in the $\zeta$ direction
and thus describes Faraday's law of induction in the steady state for an infinitely long gapped cylindrical superconductor.
The loss power of a tubular superconductor per unit tape length is readily evaluated in the limit of $w_0/R \rightarrow 2\pi$ as
\begin{equation}
P\left({w_0\over R}\rightarrow 2\pi\right)
=
2B\left({1\over2},{2n+1 \over 2n}\right)  J_{\rm c}d_0R^2\beta \left( {\beta R\over E_{\rm c}} \right)^{1/n}.
\label{loss_cylinder}
\end{equation}
Meanwhile, in the Bean limit of $n\rightarrow \infty$, the loss power of a gapped cylindrical superconductor per unit tape length is reduced to
\begin{align}
P_{\rm gapped}&
= 2J_{\rm c}d_0 \beta R^2\left[ 1-\cos(w_0/2R) \right]  \nonumber \\
&
\left\{
\begin{IEEEeqnarraybox}[\relax][c]{l's}
\approx J_{\rm c}d_0 w^2_0 \beta/4~~(w_0/R\ll 1),\\
\rightarrow 4J_{\rm c}d_0R^2\beta~~(w_0/R \rightarrow 2\pi),
\end{IEEEeqnarraybox}
\right.
\label{loss_flat-tape_tube}
\end{align}
which coincide with the loss power for a flat tape in a swept field \cite{brandt1996,higashi2019} and a tubular superconductor (the $R$ dependence is in accord with Ref.~\cite{mawatari2011}), respectively.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section*{Author summary}

\nolinenumbers

\section*{Abstract}

System identification techniques---projection pursuit regression models (PPR\reviseold{s}) and convolutional neural networks (CNNs)---provide state-of-the-art performance in predicting visual cortical neurons' responses to arbitrary input stimuli.
\reviseold{However, the constituent kernels recovered by these methods are often noisy and lack coherent structure, making it difficult to understand the underlying component features of a neuron's receptive field.}
In this paper, we show that using \reviseold{a dictionary of diverse kernels with complex shapes} learned from natural scenes based on efficient coding theory, as the front-end for PPR\reviseold{s} and CNNs can improve their performance in neuronal response prediction \reviseold{ as well as algorithmic data efficiency and convergence speed}.  Extensive experimental results also indicate that these sparse-code kernels provide important information on the component features of a neuron's receptive field. In addition, we find that models with the \reviseold{complex-shaped} sparse code front-end are significantly better than models with a standard orientation-selective Gabor filter front-end for modeling V1 neurons that have been found to exhibit complex pattern selectivity. \reviseold{We show that the relative performance difference due to these two front-ends can be used to produce a sensitive metric for detecting  complex selectivity in V1 neurons.}

\section*{Author Summary}
Convolution neural networks and projection pursuit regression models are two state-of-the-art approaches to characterizing the neural codes or the receptive fields of neurons in the visual system. However, the constituent kernels recovered by these methods are often noisy and difficult to interpret. Here, we propose an improvement of these standard methods by using a set of neural codes learned from natural scene images based on the convolutional sparse coding theory as priors or the front-end for these methods. \reviseold{We found that this approach improves the model performance in predicting neural responses with less data and with faster convergence for fitting}, \reviseold{and allows a possible interpretation of the constituents of the receptive fields in terms of the dictionary learned from natural scenes.}  \reviseold{The relative performance difference due to these two front-ends has been shown to produce an effective metric for detecting  complex selectivity in V1 neurons.}

\section*{Introduction}

The neural code of neurons in the primary visual cortex has been investigated for decades~\cite{10, 11, 12}. A number of quantitative approaches have been developed to characterize and model the receptive fields (RFs) of V1 neurons, notably energy models~\cite{1}, spike-triggered average~\cite{2, Ringach2001} and spike-triggered covariance methods~\cite{3, 4}, linear-nonlinear cascades~\cite{5}, sub-unit models~\cite{6} and general linear models based on handcrafted nonlinear feature spaces~\cite{7, 8, 9}. The typical neuronal RFs assumed or recovered by the previous methods are Gabor filter models followed by a threshold or quadratic nonlinearity~\cite{10,11,12}, for simple and complex cells respectively. However, a recent study~\cite{tang} discovered that a large proportion of neurons in the superficial layer of V1 of awake macaque monkeys are highly selective to specific complex features~\cite{4, 16}, suggesting that many V1 neurons act also as complex pattern detectors rather than just Gabor-based oriented edge detectors. \reviseold{This finding underscores the traditional models' limitations in capturing the full characteristics of neural tuning}.

\smallskip

Sparse coding or efficient coding theory~\cite{18, 28, olshausenovercomplete, CSC, gregor2010learning} has provided great insight into the computational principles underlying the structure of the receptive fields of V1 neurons.
\reviseold{In the highly overcomplete regime~\cite{olshausen2013highly}, particularly when coupled with a convolutional mechanism~\cite{CSC}}, sparse coding yields a diverse set of features, including center-surround filters, corner, curvature, and cross detectors.
In this paper, we investigate the idea that a set of basis filters learned from natural scene images based on convolutional sparse coding could be used to improve PPR\reviseold{s} and CNNs in neuronal RF recovery tasks. \reviseold{In short}, the basic idea is to use these learned codes as the front-end for PPR\reviseold{s} and CNNs to predict the neuronal responses. Our experiments show that CNNs with a diverse \reviseold{complex-shape} code front-end not only achieve the state-of-the-art performance in neuronal response prediction, but are also more data-efficient to train and \reviseold{faster in convergence}. This result suggests that  \reviseold{this more diverse set of } basis functions learned from natural scenes might provide a better approximation of the underlying components of the neurons' RFs than the \reviseold{noisy} filters \reviseold{, characterized by white-noise like power spectrum and lacking in coherent structures, often} learned by CNN trained from scratch.

Earlier studies showed that sparse coding principle leads to the development of Gabor-like receptive fields~\cite{16, 18}, whereas recent work suggested that receptive fields or kernels with more complex feature selectivity would arise when there are more neurons (basis functions) available than necessary, i.e. overcomplete representation, \revise{obtained} either with more features or denser spatial sampling (\revise{due to} convolution) ~\cite{18, olshausenovercomplete, 28, CSC, 27}. These complex features and codes are consistent with the diversity and the complexity of the feature tuning observed in V1~\cite{tang, Victor2006, Livingstone2017, Dobbins1987, Sillito1995, Ringach2002}. Therefore, we conjecture that neurons selective to higher-order complex features, such as corners, curvatures, and junctions, might be better modeled by the \reviseold{complex-shaped kernels provided by convolutional sparse coding} \revise{than by a less diverse set of kernels.} To test this conjecture, we compare the performance of CNN models with two different front-end types (\reviseold{complex-shaped} codes versus Gabor wavelets). Indeed, we found that models with the complex front-end are much better than models with the classical Gabor and/or Laplacian front-end for modeling and predicting the activities of V1 neurons, and this improvement is larger for neurons previously determined to have complex pattern selectivity. This observation
further supports \reviseold{ a plausible role} of the convolutional sparse coding theory for understanding neural codes in V1.

\smallskip
\noindent\underline{\textbf{Contributions:}}
\reviseold{We introduce two novel classes of neural response prediction models \reviseold{for receptive field characterization}:(1) convolution matching pursuit regression models (CMPRs), a variant of PPRs, which select the kernel from the sparse coding dictionaries instead of learning from scratch;
(2) fixed-kernel CNNs (FKCNNs), which are CNNs with fixed first-layer kernels.  FKCNNs  achieves state-of-the-art performance with improved data-efficiency \reviseold{and faster convergence.}}
\reviseold{Both methods allow possible interpretation of the constituents of the receptive fields in terms of the dictionary learned from natural scenes, and provide  useful metrics for characterizing and discriminating the complexity of tuning in V1 neurons.}

\section*{Materials and Methods}

\subsection*{Ethics statement}

\revise{All procedures involving animals for generating the data of this paper were in accordance with the Guide of Institutional Animal Care and Use Committee (IACUC) of Peking University Animals, and approved by the Peking University Animal Care and Use Committee (LSC-TangSM-5). }

\subsection*{Dataset: stimuli and neuronal responses}

The neuronal data studied in this paper were published in a previous study~\cite{tang}. They are the responses of neurons in the V1 superficial layer of two awake (Macaca mulatta) monkeys A and E, with respect to a set of stimuli. They are obtained using large scale two-photon imaging with the calcium indicator GCaMP5s as reported in~\cite{tang}. The calcium signals in response to visual stimuli of 1142 neurons from monkey A and 979 neurons from monkey E are imaged during a fixation task. The response of a cell is computed as the standard $\Delta F/F_0$, based on the averaged activity within an ROI (region of interest) during stimulus presentation in each trial.

The stimulus set is designed to test the hypothesis that many
neurons in V1 are not only orientation-tuned (OT) but also selective
to specific high-order (HO) complex features. The stimulus set contains 9500 binary (black and white) images generated from 138 basic prototypes, shown in Fig~\ref{fig2} by rotating and scaling. These prototypes are grouped into five major categories:
orientation stimuli (bars and gratings), curvature
stimuli (curves, solid disks, and concentric rings),
corner stimuli (lines or solid corners), cross stimuli
(lines crossing one another), and composition
stimuli (patterns created by combining multiple elements
from the first four categories).
The entire stimulus set consists of 1,600 OT patterns and 7,900 more complex HO patterns. The full set was shown to monkey A while a half-set (with half of the rotations) was shown to monkey E. Note that even for the half set, oriented bars at 48 orientations with rotation increment of $7.5^{o}$ were tested.

Each stimulus is shown in a $3^{o} \times 3^{o}$  (90 x 90 pixels) aperture. The neurons possess overlapping receptive fields, all within a radius of $0.5^{o}$ in visual angle. All recorded neurons have classical RFs of diameters
well below $1^{o}$ in visual angle around the
stimulus center~\cite{tang}. Hence, to reduce the number of parameters in the models, we crop the central 40 $\times$ 40 pixels of all stimuli and downscale them to 20 $\times$ 20 pixels as model input. The input size is at least twice the size of any neuronal RFs and thus covers all of them. In this study, we discard neurons without a significant preference for stimulus patterns. Specifically, we discard neurons that have less than 50 stimuli with response $\Delta F/F_0$ greater than 0.2, where 0.1 is one standard deviation of the measurement. This measure reduces our samples to 781 neurons for monkey A and 632 neurons for monkey E. \reviseold{We used this selection criteria to ensure there is a sufficient number of stimuli with responses two standard deviations above noise level. We assumed that is necessary for the approaches to work well. Most of the neurons we eliminated based on this criteria have a low response to all pattern stimuli, rather than being very sharply tuned \revise{(see Supplementary Information)}.}

\begin{figure}[htp]
  \centering
  \includegraphics[width=8.52cm, height=6.0cm]{images/new_fig_1.png}
  \caption{{\bf Image stimuli:} The 138 prototype stimuli, including simple orientation patterns as well as complex ones, which generate all 9,500 stimuli used through rotation and translation.}
  \label{fig2}
\end{figure}

\subsection*{Gabor code and \reviseold{complex-shape} sparse code front-end }

The key contribution of this paper is the use of neural sparse codes derived from natural scenes as a front-end for receptive field models trained to characterize and predict neuronal responses. According to sparse coding theory~\cite{olshausenovercomplete}, the number of neurons used for natural scenes encoding can determine the complexity and diversity of the neural sparse codes. For example, for an input patch of 12 x 12 pixels, each \revise{gray-level} image is a vector in 144-dimensional space, with each pixel representing one dimension. When the number of neurons allocated to this patch is equal to 144, the representation is called a complete representation; otherwise, the representation will be either undercomplete or overcomplete. It has been shown that when the resource constraint is complete \revise{or} undercomplete, the basis functions learned are mostly Gabor filters or Gabor wavelets~\cite{18}. But in the overcomplete scenario, the basis learned tends to exhibit greater diversity and \reviseold{complexity~\cite{18, olshausenovercomplete, CSC, 27, 28, olshausen2013highly}}.

\reviseold{In this study, we use the complex-shaped kernels obtained using the convolutional sparse coding (CSC) method~\cite{CSC}.} These are feedforward kernels that compute the coefficients of a set of complementary basis by linear projection for reconstructing the input images with the coefficients subject to sparsity constraint.
The use of convolution in sparse coding is shown to increase efficiency of patch-based sparse coding~\cite{CSC} that otherwise tends to produce Gabor filters or Gabor wavelets. CSC yields a variety of complex tunings reminiscent of the feature selectivity observed in~\cite{tang}. We train this unsupervised CSC model with a kernel size of 9 $\times$ 9 pixels on 1,000,000 25 $\times$ 25px image patches sampled from the 5000 128 $\times$ 128px natural scene images used in the original paper~\cite{CSC}. The model is trained with dictionaries of either 16, 32, 64, 128, or 256 channels. Some of the learned dictionaries are shown in Fig~\ref{fig1}. We observed that as we continue increasing the number of channels of the dictionary past a certain point, filters appear to have more redundancy instead of showing more diversified shapes. This model is complete in the sense that there are more than 625 neurons in the first hidden layer for an input image patch of 625 pixels. \revise{In early sparse coding works [16, 17], when convolution is not used, there is a greater amount of diversity of feature detectors in each hypercolumn. With convolution, overcompleteness can be achieved by a small set of feature detectors sampling the visual space more densely.  We adopt the convolution strategy to minimize the number of parameters in the model.}

\begin{figure}[!h]
\centering
  \includegraphics[width=5cm, height=5cm]{images/64_filter.png}
  \includegraphics[width=3cm, height=5cm]{images/128_9_filter.png}
  \includegraphics[width=5cm, height=5cm]{images/256_9_filter-1.png}
\caption{{\bf \reviseold{Complex-shaped convolutional} sparse code dictionary:} Learned feature dictionaries from the CSC algorithm \revise{applied to grayscale natural images} with 64, 128, and 256 elements, each with a kernel size of 9px.
}
\label{fig1}
\end{figure}

In our study, we compare the use of these \reviseold{complex-shaped} sparse code front-ends with Gabor wavelet front-ends \cite{daugman1988, lee96}, which have long been the mainstream models for V1 neurons. We also compare them against the Laplacian of Gaussian (LOG) filter front-end, which has been used to model retinal ganglion cells and V1 cells with center-surround RFs. Fig~\ref{fig5} shows an example of a set of Gabor~\cite{lee96} and LOG wavelets in contrast to a set of \reviseold{complex-shaped} sparse codes selected from the learned dictionary of the CSC algorithm. Note that the Gabor wavelet codes and the \reviseold{complex-shaped} sparse codes have similar coverage in the spatial frequency domain, as shown in the aggregated power spectrum in Fig~\ref{fig5}. \revise{Their spatial frequency spectra highly overlap with each other (see Supplementary Information)}, suggesting that their difference in performance is not likely due to different coverage in the frequency domain.

\begin{figure}[htp]
  \centering
  \includegraphics[width=13cm]{images/fig3.png}
  \caption{\reviseold{{\bf Gabor, Laplacian of Gaussian versus \reviseold{complex-shaped} sparse code front-ends:} From left to right we have: \textbf{A}: Gabor wavelets with 3 scales and 8 orientations. \textbf{B}. Laplacian of Gaussian kernels with 3 scales. \textbf{C}. The aggregated power spectrum of the Gabor wavelets. \textbf{D}: 24 \reviseold{complex-shaped} sparse kernels selected from the 64-element CSC model for use in our Fixed Kernel CNN and Projection Pursuit models. \textbf{E}: Aggregated power spectrum of the \reviseold{complex-shaped} sparse codes. }}
  \label{fig5}
\end{figure}

\subsection*{Receptive field models for neuronal response prediction}

In this section, we study the effect of incorporating the \reviseold{complex-shaped} sparse code front-end into two classes of supervised learning models for modeling V1 neurons. The first class is the pursuit regression models, which perform layer-by-layer non-parametric regression. We extend the standard PPR approach~\cite{21, 23, 22} in two ways: first, by incorporating the convolution operation, and second, by incorporating the \reviseold{complex-shaped} sparse code as a front-end, i.e. using these codes as the feature kernels in the first layer of the model.
The second class is CNNs. We compare baseline CNN models with kernels learned from scratch by a data-driven approach~\cite{zhang}, against models using fixed sparse code kernels in their first layer. One model of each type is learned for each neuron in monkey A and E. For fair comparison, we constrain the number of the parameters of the proposed models to be the same or lower than the baseline models.

\smallskip
\noindent\underline{\textbf{Projection pursuit regression (PPR)}}~\cite{21} has been previously used to model neural responses, achieving remarkable performance~\cite{liu2016spatial} \reviseold{and allowing analysis of receptive field structure~\cite{20}}. This method does not require stimuli with specific statistical properties and is thus suitable for analyzing the neuronal responses to our complex pattern stimuli. The key idea of PPR is to optimize one kernel of the RF model at a time, allowing the first learned kernel to provide important information about the neuronal RF while the consecutive kernels only fit the residuals of the previous ones. This means that instead of the structure of preferred inputs being distributed evenly across many kernels, as in a standard CNN, it is concentrated as much as possible into a single, more \reviseold{structurally coherent} kernel.

The detailed PPR algorithm is shown in Algorithm~\ref{alg1}, where $\cdot$ denotes the dot product and \reviseold{$g(x)$ stands for a quadratic function $a*x^2 + b*x + c$ with learned parameters $a, b, c$}. This algorithm will sequentially learn layers of kernels $F$, which are the same size as the input image $x$. The regression of the neuronal response $y$ on the input stimulus $x$ is defined as the dot product between $x$ and the learned kernel, followed by a quadratic nonlinearity. After all layers have been learned, the kernels will be reordered based on pursuit contribution $I$, and redundant kernels will be eliminated. Finally, the function $\psi(S, F) = \sum_{m=1}^{|F|} g({F}_{m} \cdot S)$ will be used to predict the neuronal response with respect to input stimulus $S$.

\begin{algorithm}\
\small
\caption{PPR Algorithm}\label{alg1}
\begin{algorithmic}[1]
\Procedure{PPR}{$X,Y,M$}\newline \Comment{input set of images $X$, corresponding neuronal responses $y$, and number of kernels to be learned $M$.}
\State $\text{Define quadratic function } \reviseold{g(x) = a*x^2 + b*x + c}$
\State ${r}\gets {Y}, m\gets 1$
\While{$m\leq M$}
\State $I(F'_m) = 1 - \sum_{i=1}^{n} ({r}_{i} - g(F'_{m} \cdot X, 2))^2 / \sum_{i=1}^{n} {r}_{i}^2$
\State ${F}_{m}\gets argmax(I(F'_m))$
\State ${r}\gets {r} - \reviseold{g({F}_{m} \cdot X)}$
\State $m\gets m + 1$
\EndWhile\label{euclidendwhile}
\State $F \gets \{F_1, \cdots, F_M\}$ \\
\Return $\text{reduce\_redundant(F)}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\smallskip
\noindent Next, we introduce two variants of PPR that improve the neuronal prediction results and kernel interpretability. The first variant is
{\em Convolutional projection pursuit regression}, which introduces the convolution operation~\cite{lecun1995convolutional} to the kernel layers. The second variant is {\em Convolutional matching pursuit regression}~\cite{22}, which limits the feature kernels selected by the convolutional projection pursuit regression to a fixed dictionary of filters.

\smallskip
\noindent\underline{\textbf{Convolutional projection pursuit regression (CPPR)}}:
The original PPR model requires the kernels $F$ to be the same size as the input images, leading to learning redundant parameters. More importantly, monkey V1 neurons' receptive field sizes are usually much smaller than image stimuli, and thus the kernels learned by PPR can not accurately describe the size of the neuronal RF. Here, we propose a modified version of PPR, CPPR, which allows using a smaller kernel of size $d\times d$ and performs convolution over the input images of size $w\times w$ at each layer. Then the result after convolution (of shape $(w-d+1) \times (w-d+1)$) goes through an \reviseold{absolute value layer} and a max-pooling layer to get the output, as shown in Equation~\ref{equ1}. All other procedures are the same as in Algorithm~\ref{alg1} except that we take convolution between kernels and the image instead of the simple dot product. Ideally, CPPR can learn kernels that are the same size as the receptive field, and the receptive field location can be chosen by the max-pooling layer.

\begin{equation}\label{equ1}
    \psi_i(X, F) = \sum_{m=1}^{|F|} \reviseold{g \ ( \ max\ (\ abs({F}_{m} \ast X)))}
\end{equation}
In Equation~\ref{equ1}, $\ast$ denotes the convolution operation, $abs$ denotes the $absolute value$ layer and $max$ denotes the max pooling layer.

\smallskip
\noindent\underline{\textbf{Convolutional matching pursuit regression (CMPR)}}:
We find that most of the kernels learned from PPR or CPPR algorithm tend to be rather noisy and lack coherent forms. \reviseold{We define the term ``noisy'' as having high power at high spatial frequencies.} Inspired by the matching pursuit regression algorithm~\cite{23}, we design CMPR methods to incorporate the \reviseold{complex-shaped} sparse code front-end into the CPPR models. Specifically, we modify the CPPR approach by selecting the feature kernels $F$ exclusively from an existing dictionary, instead of learning them from scratch. Here, we take the \reviseold{complex-shaped} sparse codes with 64 kernels (Fig~\ref{fig2}) and expand them eight-fold by rotating each component at $45^{o}$ step by step to allow better orientation steering of the kernels. The resulting dictionary $D$ has 512 kernels. If the dictionary $D$ spans the space which we intend to learn from CPPR, then CMPR should have similar performance to CPPR. More importantly, the kernels selected using the CMPR model will be much \reviseold{less noisy} than those learned by CPPR. We also tried to use \reviseold{complex-shaped} sparse codes without rotational expansions as the dictionary, but the performance was substantially worse.

\smallskip
\noindent\underline{\textbf{Convolutional neural networks (CNNs)}} are the state-of-the-art models for predicting neuronal responses~\cite{zhang, 24, 29, 31, Yamins2014, Yamins2016}. They outperform all other baseline models, such as Gabor-based
standard models for V1 cells and the various variants of generalized linear models~\cite{7, 8, 24, 29}.

We use the structure of an earlier paper~\cite{zhang}, which is optimized for this set of pattern data. A one-layer CNN model passes the input images through a series of linear-nonlinear operations, each of which consists of convolution, ReLU nonlinearity, and max pooling (Fig~\ref{fig4}). Finally, the outputs of the above operations are linearly combined to form the predicted neuronal response. Here we experiment with two specific CNN structures also tested in~\cite{zhang} for comparison: one with 4 kernels in the convolutional layer ($CNN_4$) and another with 9 kernels ($CNN_9$). The earlier study~\cite{zhang} explored different CNN structures and established that adding more layers to the model does not increase the fitting ability of the model for this dataset, \revise{so we do not experiment with deeper models.} $CNN_9$ has been shown to be the state-of-the-art model for this dataset, superior to a variant of the projection pursuit methods.

\begin{figure}[htp]
  \centering
  \includegraphics[width=10cm]{images/new_CNN_str.png}
  \caption{{\bf Structure of CNN:}
  The input image is of size $20\times 20$ pixels, downsampled from the $40\times 40$ image of the experimental stimulus centered at the neurons' aggregated receptive fields. The model has three stages: convolution and thresholding, max-pooling, and a final fully-connected linear combination to generate the predicted neuronal response. \\
  }
  \label{fig4}
\end{figure}

\smallskip
\noindent\underline{\textbf{Fixed kernel convolutional neural networks (FKCNNs)}}:
Even though the standard CNN models have the best predicting accuracy, the recovered features (i.e. the 4 or 9 kernels in the first layer) in the hidden layer often do not have \reviseold{coherent} forms, making it difficult to understand the constituent components of the neuronal RF. Here, \reviseold{ we explore whether  the use of the sparse code front-end can improve its performance by reducing overfitting and whether the constituent components of a neuron's receptive field can be related to the statistical priors of natural scenes.}

FKCNNs use the baseline CNN model's structure (Fig~\ref{fig4}), but replace the first layer with either the sparse code dictionaries learned from the sparse coding algorithm (Fig~\ref{fig1}) or Gabor wavelets, which are considered standard receptive field models of V1. During training, the substituted first layer is kept fixed. In order to make a fair comparison with the standard CNNs, we fix the number of parameters learned by the FKCNNs to be roughly the same. Specifically, we select 24 distinct complex or Gabor kernels of size $9 \times 9$ as the first convolutional layer's kernels. This results in the same number of parameters to learn as the $CNN_4$ and half the number of parameters as in $CNN_9$, since the kernels themselves are not learned. The 24 kernels are selected in the following manner: we start with a set of 64 kernels learned with CSC (Fig~\ref{fig2}) and discard the kernels that are too noisy or redundant. Then we randomly subsample 20 sets of 24 kernels as the fixed first layer from the remaining 35 kernels. We construct FKCNN models from these sets of 24 kernels as well as the entire set of 35 kernels and tested their performance on the pattern dataset. The results have an average performance in terms of Pearson correlation of 0.445 with a standard deviation of 0.0079, comparable to the normal CNN models, suggesting that a random subset with size 24 from these 35 kernels is complete enough to model this dataset.

In addition to the \reviseold{complex-shaped} sparse code dictionary, we also test the Gabor wavelet ( Gabor-CNN or GCNN) as shown in Fig~\ref{fig5}. Since some V1 neurons are known to have center-surround receptive fields, we also add three Laplacian of Gaussian feature channels of different scales to supplement the Gabor wavelet dictionary (we denote it as GLCNN). The purpose of introducing GCNN and GLCNN here is to assess whether the \reviseold{complex-shaped} sparse codes provide a better dictionary for modeling V1 neurons than the standard models.

\noindent\underline{\textbf{\reviseold{Comparison of model architectures of the two approaches:}}}

\reviseold{The PPR-based methods and the CNN-based methods we explored for neural response prediction  have the following three notable differences in their designs: a) for training, PPR-based methods train one kernel at a time, in contrast to CNN, which train all kernels and parameters simultaneously; b) for CNNs, the nonlinearity is ReLU, while PPR methods use the absolute value or quadratic non-linearity; c) for CNNs, the max pooling retains some spatial information, and PPR methods pool everything.}

 \reviseold{These are the defining and intrinsic features of the two approaches. PPR emphasizes finding the primary components and then fitting the successive residuals while CNNs seek to find potentially many different ways to compose the receptive field.  As for the choice of nonlinearity, a quadratic nonlinearity is used in the original PPR method, giving PPR-based models three tunable parameters and thus extra degrees of freedom over a ReLU or sigmoid function. This is an intrinsic advantage of the PPR-based methods. In CMPR, the filters are fixed, and only the prediction scale and bias need to be adjusted; using ReLU or sigmond will deprive them of learnable parameters. As for the pooling choice, the key idea of the projection pursuit is to find the best filter at the best location, which max pooling addresses.  Retaining some spatial information in CNN-based methods confers more flexibility and expressiveness, allowing the CNN to synthesize the receptive fields by combining information across both feature channels and space. This is an intrinsic advantage of the CNN approach. Thus, even though in theory, the two approaches could use different pooling and nonlinearity strategies, our design choices are typical and suit the relative advantages of their respective methods.  }

\section*{Results}
In this section, we first evaluate the performance of different models quantitatively and then discuss the model's interpretability from various aspects. \reviseold{In particular, our quantitative performance evaluation is based on neuronal response prediction accuracy, convergence rate, and data efficiency}. We use the Pearson correlation between predicted and actual neuronal responses as an evaluation metric to measure prediction accuracy, which is widely used in many papers~\cite{4, 24, zhang}. An example can be found in Fig~\ref{fitting_curve}.

\begin{figure}[htp]
  \centering

  \includegraphics[width=8cm]{images/neu_stderr_fitting.png}
  \caption{ \reviseold{{\bf Predicted tuning curve versus actual tuning curve:} An example of a real neuronal tuning curve to \revise{the 1900 stimuli from} the testing dataset (not used during training), ordered by the response magnitude (black curve, \revise{with $\pm$ 1 standard error estimated from the trial-by-trial variance}). The red curve shows the model's predicted response to these stimuli. The Pearson correlation (0.61 in this example) between these two curves is calculated as a quantitative evaluation metric.}}
  \label{fitting_curve}
\end{figure}

\subsection*{Quantitative Evaluation of the Methods}

\noindent\underline{\textbf{Synthetic Data:}} We evaluate these novel methods by first testing them using synthetic neurons. The synthetic neurons are modeled as the nonlinear transform of one or the sum of multiple linear filter outputs. We choose the ReLU activation function as our nonlinear transform because previous work found that all activation functions behave essentially the same~\cite{zhang}. Then, they are stimulated by a set of Gaussian white noise images (20 x 20 pixels), and the response is given by the following linear-nonlinear model, as described by Equation~\ref{equ2}. This is a common ground truth evaluation method, also used by~\cite{4, 24}.

\begin{equation}\label{equ2}
    y_i = \sum_{s \in S} ReLU(I_i \cdot s) + \eta
\end{equation}
where $\cdot$ denotes the inner dot product between the input image and the artificial RF $s$; $ReLU$ is the rectified linear unit function that maps negative values to 0; $\eta$ represents Gaussian noise with zero mean and variance equal to $\frac{\sqrt{I_i \cdot s}}{10}$; the set of artificial RFs $S$ is a subset of $ 9 \times 9$ filters padded with zeros in the surrounding to make them the same size as the input stimulus ($20 \times 20$). We select $S$ either from our learned \reviseold{complex-shaped} sparse codes or cropped patches of our pattern stimuli.
We train the models with 10000 samples of noise inputs $I$ until convergence, and then test them with another set of 2000 noise input. As shown in Fig~\ref{fig6_2}, for these simple neurons, all the approaches described in the previous section are able to achieve near perfect testing performance when we only have one RF in $S$. However as we gradually increase the number of filters in $S$, the pursuit regression models' performance seems to suffer whereas the CNN models roughly retained their good performance.
Fig~\ref{fig6_2} shows one example of the recovered receptive fields using different pursuit regression models and CNN's visualization~\cite{26}. Even though the visualization of the CNN model of the neuron is consistent with $S$, the recovered kernels in the hidden layer are very noisy and impossible to interpret, whereas the FKCNN's component kernels  \reviseold{are statistical priors of natural scenes with coherent structured patterns}. Specifically, panel D of Figure~\ref{fig6_2} demonstrates the top 4 ranked kernels of FKCNN method. \reviseold{We can see that the first kernel (the most important one) matches with the artificial RF. In fact, the importance value of the first kernel, is much higher than those of the rest kernels. Therefore in this case, we have reasons to believe that the first kernel resembles the RF.
For the second top-ranked kernel, it matches with the vertical bar in the RF. Therefore, it can partially response to testing stimulus and derive neuronal responses that has a relatively high correlation with the actual neuronal responses.}

\begin{figure}[htp]
  \centering
  \includegraphics[width=9cm]{images/GT_performance.png}
  \includegraphics[width=11cm]{images/GT_ALL.png}
  \includegraphics[width=5cm]{images/CNN_kernel_GT.png}
  \includegraphics[width=5cm]{images/FKCNN_GT.png}
  \caption{\reviseold{{\bf Model performance on artificial data:} This figure shows the performance (A) and the recovered receptive fields (B-D) of the different methods. On the second row (B) from left to right, the seven images are the receptive field $s$ recovered from each model type: PPR's filter, CPPR's filter, CMPR's selection, CNN's visualization, FKCNN's visualization, and GaborCNN's visualization. On the last row, the left images (C) are the learned kernels of the CNN. The right image (D) shows the top 4 highest-contribution kernels of the FKCNN model, \reviseold{whose importance scores are 0.782, 0.503, 0.469, and 0.391 respectively.}}}
  \label{fig6_2}
\end{figure}

\smallskip
\noindent\underline{\textbf{Real Neuronal Data:}} We test all seven methods on the calcium imaging neuronal data.
For the three projection pursuit based approaches, all models can select up to 10 kernels. The CMPR model was tested using a 512 feature dictionary derived from 64 \reviseold{complex-shaped} sparse code filters learned from CSC (Fig~\ref{fig1}), with 8 rotations each. The CMPR and CPPR models have roughly the same number of parameters to learn, which is about 25\% of that used in PPR. It is worth noting that CMPR takes longer to train because of the selection process. Fig~\ref{fig8} shows that CPPR and CMPR achieve prediction accuracy better than PPR with fewer parameters, which underscores the efficacy and efficiency of the convolution approach~\cite{zhang}.

\begin{figure}[htp]
  \centering
  \includegraphics[width=6.5cm]{images/performance_ALL_A.pdf}
  \includegraphics[width=6.5cm]{images/performance_ALL_E.pdf}
  \caption{{\bf Model performance on neuronal data:} \revise{Box-plots of response prediction performance based the Pearson correlation  between true and predicted responses of all modeled neurons. The performance of all models on Monkey A is shown on the left and performance on Monkey B is shown on the right.} Overall, the performance is significantly better for Monkey E because of higher-quality neuronal recordings for this monkey.}
  \label{fig8}
\end{figure}

\begin{table}[!ht]
    \centering
    \scalebox{0.78}{
    \begin{tabular}{c|c|c|c|c|c|c|c|c}
    \hline
        Model & PPR & CPPR & CMPR & $CNN_4$ & $CNN_9$ & GCNN & GLCNN& FKCNN  \\ \hline

        PPR & [0.278] &-32.7\%$^{***}$ & -30.9\%$^{***}$ &-48.9\%$^{***}$ & -61.2\%$^{***}$ & -35.6\%$^{***}$ & -42.1\%$^{***}$ & -62.9\%$^{***}$ \\ \hline

        CPPR & 32.7\%$^{***}$ & [0.369] & 1.3\% & -12.2\%$^{*}$ & -21.4\%$^{**}$ & -2.2\% & -7.0\%$^{*}$ & -22.8\%$^{**}$ \\ \hline

        CMPR & 30.9\%$^{***}$ & -1.3\% & [0.364] & -13.7\%$^{*}$ & -23.0\%$^{**}$ & -3.6\% & -8.5\%$^{*}$ & -24.5\%$^{**}$ \\ \hline

        $CNN_4$ & 48.9\%$^{***}$ & 12.2\%$^{*}$ & 13.7\%$^{*}$ & [0.414] & -8.2\%$^{*}$ & 9.8\%$^{*}$ & 4.8\%$^{*}$ & -9.4\%$^{*}$ \\ \hline

        $CNN_9$ & 61.1\%$^{***}$ & 21.4\%$^{**}$ & 23.0\%$^{**}$ & 8.2\%$^{*}$ & [0.448] & 18.8\%$^{**}$ & 13.4\%$^{**}$ & -1.1\% \\ \hline

        GCNN & 35.6\%$^{***}$ & 2.1\% & 3.5\% & -9.8\%$^{*}$ & -18.8\%$^{**}$ & [0.377] & -4.8\%$^{*}$ & -20.2\%$^{**}$ \\ \hline

        GLCNN & 42.1\%$^{***}$ & 7.0\%$^{*}$ & 8.5\%$^{*}$ & -4.8\%$^{*}$ & -13.4\%$^{*}$ & 4.8\%$^{*}$ & [0.395] & -14.7\%$^{**}$ \\ \hline

        FKCNN & 62.9\%$^{***}$ & 22.8\%$^{**}$ & 24.5\%$^{**}$ & 9.4\%$^{**}$ & 1.1\% & 20.2\%$^{**}$ & 14.7\%$^{***}$ & [0.453]  \\ \hline

    \end{tabular}
    }
    \caption{
    \revise{
{\bf Relative performance difference between the different methods on Monkey A:} The average performance of all models are put on the diagonal of this table.
The relative performance difference is calculated as (A-B)/min(A, B), where A is the average performance of the method on the row and B is that of the method on the column.\\
* refers to a $p$-value between 1e-10 and 0.05 \\
** refers to a $p$-value between 1e-40 and 1e-10 \\
*** refers to a $p$-value between 0 and 1e-40.\\
Otherwise, it means that the difference between two methods is not significant.}
}
    \label{tab1}
\end{table}

Fig~\ref{fig8} and Table~\ref{tab1} show that our FKCNN model performs significantly better than the $CNN_4$ model which has the same number of parameters and comparable to $CNN_9$ model (state-of-the-art for this set of data~\cite{zhang}) for both monkeys. The fact that the FKCNN's \reviseold{complex-shaped} sparse code dictionary serves as an effective front-end suggests that the \reviseold{complex-shaped} sparse code might indeed span a similar or even better feature space than that spanned by the noisy constituent components of CNN for characterizing the real neurons.

\smallskip
\noindent\underline{\textbf{Interpretation of Recovered Receptive Field Components:}}
The CMPR and CPPR models have comparable performance, but the receptive fields recovered by the CPPR method often do not have an \reviseold{coherent structure}. \reviseold{Fig~\ref{fig9} compares the kernels recovered by the standard PPR method (Images (B1) and (B2)), those recovered by CPPR (Images (C1) and (C2)), and those recovered by CMPR (Images (D1) and (D2)) for two example neurons. The PPR and CPPR methods yield filters that are noisy, while CMPR reveals the understandable key component features preferred by the neurons.}

\reviseold{In Fig~\ref{fig9}, the kernels from the CNN model's convolution layer (Images (E1) and (E2)) appear to be noisy and lack coherent structure, whereas the top FKCNN's component kernels are
curves and edges(Images (F1) and (F2))}. We rank the FKCNN's 24 kernels by their importance scores and showed the top four kernels in Fig~\ref{fig9}(fourth row). \revise{The importance score of each kernel component in Fig~\ref{fig9} is the response prediction performance on the testing data performance in terms of Pearson correlation when we retain only that filter and delete all other filters in the FKCNN after the training process}.
\revise{We can see that the top-ranked FKCNN kernels (shown in F1 and F2 of Fig~\ref{fig9}) are reasonable components for constructing the top preferred stimuli (shown in A1 and A2 of Fig~\ref{fig9}).} \revise{ While the receptive fields of the neurons are difficult to visualize because of the non-linearity,  the most preferred stimulus extracted from the entire stimulus set by each model is a reasonable reflection of the receptive fields characterized by the model. We show the most preferred stimulus of each model to the left of the top kernels of each of the models. We observe that the top stimulus extracted by CMPR and FKCNN are indeed among the top 3 preferred stimuli of these two neurons, while models based on PPR, CPPR and CNN seem to have identified less relevant features.}

\begin{figure}[htp]
  \centering
  \includegraphics[width=10cm]{images/fig8_new.pdf}
  \caption{{\bf The recovery of two example neurons' receptive field components:} Two neurons are shown in the figure, separated by the double lines. Images (A1) and (A2) show each neuron's top 20 preferred stimuli (size 40x40). \revise{The most preferred stimulus and the two kernels (size 20x20) of PPR are shown in images (B1) and (B2) for two example neurons respectively. The most preferred stimulus and the two kernels of CPPR are shown in images (C1) and (C2), and those of CMPR are shown in images (D1) and (D2). The most preferred stimulus and four learned convolutional kernels from $CNN_4$ (size 9x9) for each neuron are shown in images (E1) and (E2). The most preferred stimulus and four highest-importance kernels in the FKCNN model (size 9x9) and their importance scores underneath the kernels are shown in images (F1) and (F2).}}
  \label{fig9}
\end{figure}

\subsection*{\reviseold{Data efficiency and convergence rate}}

\reviseold{In this section, we compare the FKCNN, GCNN, CPPR, and CMPR with the baseline CNN and PPR models on two additional aspects beyond the already-discussed test performance: data efficiency and convergence rate. The results are shown in Fig~\ref{fig11} and the statistical significance (tested with paired t-test) of performance differences is reported in Table~\ref{tab-fig9}}.

\smallskip
\noindent\underline{\textbf{Data efficiency:}} \reviseold{In order to test the data efficiency of different models, we randomly subsample a fraction (100$\%$, 50$\%$, 25$\%$, or 12.5$\%$) of the training set to train the model while keeping the validation and test sets fixed to evaluate the model performance. Then we measure the drop of model performance with regard to the performance trained with $100\%$ and report this drop and statistical significance of Monkey A only in Table~\ref{tab-fig9}, because we observe the same pattern for Monkey E.}
\reviseold{From Plot A and B, we can see that FKCNN and GCNN require less data than the baseline $CNN_9$ to achieve the same or better predictive performance.}
\reviseold{From Plot C and D, we can see that CMPR requires less data than CPPR to achieve the same or better predictive performance.}

\smallskip
\noindent\underline{\textbf{Convergence rate:}} FKCNN and GCNN converge much faster than $CNN_9$ during training, i.e. the number of epochs needed for the model to have stabilized test performance is fewer than that of $CNN_9$. The models with the fixed front-end converge faster and require less data during training because the gradients do not need to propagate to the first convolutional layer as opposed to the baseline CNN models.

\begin{figure}[htp]
  \centering
  \includegraphics[width=4cm]{images/data_MA.png}
  \includegraphics[width=4cm]{images/Data_ME.png}
  \includegraphics[width=4cm]{images/monkeyA_dataeff_ppr.png}
  \includegraphics[width=4cm]{images/monkeyE_dataeff_ppr.png}
  \includegraphics[width=4cm]{images/conv_monkeyA_PPR.png}
  \includegraphics[width=4cm]{images/conv_monkeyE_PPR.png}
  \caption{\reviseold{{\bf Model data-efficiency and convergence test:} Plots A, B, C and D show model performance as a function of the proportion of data used for each monkey and model type. Plots E and F summarize the distributions of training times (measured in epochs) across neurons for each model type and monkey in box-and-whisker plots. Please note that we did not compare the PPR, CPPR, and CMPR's convergence rate because they each have different training procedures, which don't allow fair comparison.}}
  \label{fig11}
\end{figure}

\begin{table}[h]
    \centering
    \scalebox{0.9}{
    \reviseold{
   \begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
{\bf Model} & \multicolumn{3}{|l|}{\bf Performance Drop} & \multicolumn{3}{|l|}{\bf T-test of averaged difference}\\ \hline
 & $50\%$ Data & $25\%$ Data & $12.5\%$ Data & FKCNN & GCNN & CNN$_9$ \\ \hline
FKCNN & $3.3 \%$ & $17.9 \%$ & $33.1 \%$ & / & - & +++ \\ \hline
GCNN & $4.7 \%$ & $18.8 \%$ & $28.9 \%$ & + & / & +++ \\ \hline
CNN$_9$ & $8.1 \%$ & $25.2 \%$ & $38.6 \%$ & - - - & - - - & / \\ \hline
\end{tabular}}
}

\scalebox{0.9}{
\reviseold{
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
{\bf Model} & \multicolumn{3}{|l|}{\bf Performance Drop} & \multicolumn{3}{|l|}{\bf T-test of averaged difference}\\ \hline
 & $50\%$ Data & $25\%$ Data & $12.5\%$ Data & PPR & CPPR & CMPR \\ \hline
PPR \  & $13.6 \%$ & $17.6 \%$ & $35.5 \%$ & / & +++ & - \\ \hline
CPPR \  & $8.8 \%$ & $34.0 \%$ & $41.4 \%$ & - - - & / & - - - \\ \hline
CMPR \  & $5.1 \%$ & $27.4 \%$ & $35.8 \%$ & + & +++ & / \\ \hline
\end{tabular}}}
    \caption{\reviseold{\textbf{The performance drop when using less data on Monkey A.} For statistical significance, we adopt the same notation from Table~\ref{tab1}.}}
    \label{tab-fig9}
\end{table}

\subsection*{Assessment of Higher Order Selectivity}

\noindent\underline{\textbf{HO versus OT selectivity:}} The standard models of V1 receptive fields are Gabor filters~\cite{1, 2}. Yet, we find that the FKCNN with a \reviseold{complex-shaped} sparse code front-end outperforms that with a Gabor front-end, even with Laplacian of Gaussian filters added (Fig~\ref{fig8}). In an earlier paper~\cite{tang}, we classified these neurons into simple orientation tuned (OT) and \reviseold{higher-order tuned} (HO) groups using a very stringent criterion. Specifically, a neuron is categorized as OT if there exists one or more single oriented bar or edge stimuli among the set of stimuli that elicit a response above 50\% of the peak response of that neuron to the entire stimuli set. Thus, the neuron can be considered to exhibit higher-order (HO) selectivity only when all the preferred stimuli are HO. About 40\% of the neurons in the superficial layer of V1 are classified as HO in this manner. The most preferred patterns of these neurons highly resemble the dictionary filters learned from the sparse coding algorithm \reviseold{(some examples of this can be found in section 3 of the Supplementary Information)}. Thus it is intuitive to conjecture that the \reviseold{complex-shaped} sparse code FKCNN would outperform the Gabor-based one (GCNN) by a larger margin for cells classified as HO than those classified as OT. Table 3 shows that to be the case on average, with the relative gain from the complex FKCNN over the GCNN being 19\% higher for HO cells in Monkey A, and 18\% higher for Monkey E. This correspondence provides supporting evidence for the previous work's claim that HO neurons' receptive fields are composed of more complex feature elements.

\begin{table}[h]
    \centering
    \scalebox{0.9}{
    \reviseold{
   \begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
\multicolumn{1}{|l|}{\bf } & \multicolumn{2}{|l|}{\bf Monkey A} & \multicolumn{2}{|l|}{\bf Monkey E}\\ \hline
$Model$ & FKCNN & GCNN & FKCNN & GCNN\\ \hline
OT+HO neurons & $0.418 \pm 0.0037 $& $0.362 \pm 0.0036$ & $0.550 \pm 0.0043$ & $0.483 \pm 0.0041$\\ \hline
OT neurons & $0.424 \pm 0.0038$ & $0.389 \pm 0.0032$ & $0.588 \pm 0.0045$ & $0.544 \pm 0.0037$\\ \hline
HO neurons & $0.408 \pm 0.0035$ & $0.319 \pm 0.0040$ & $0.488 \pm 0.0041$ & $0.387 \pm 0.0047$ \\ \hline
\end{tabular}}}
    \caption{\reviseold{{\bf Average and standard error model performance on neurons with different types.}}}
    \label{tab2}
\end{table}

To further probe the greater improvement for the complex front-end for HO neurons, beyond the average difference, we propose the following complexity score to test whether a neuron's higher-order selectivity can be predicted by the level of that improvement. This score measures the superiority of the complex basis over the Gabor-based basis for a given neuron.

\begin{equation}\label{equ3}
    score(neuron) = \frac{Corr_{FKCNN}(neuron) - Corr_{GCNN}(neuron)}{Corr_{FKCNN}(neuron) + Corr_{GCNN}(neuron)}
\end{equation}

We find that the HO and OT neurons indeed form visibly distinct distributions (Fig~\ref{fig12}) under this measure (equ~\ref{equ3}). A linear classifier trained to distinguish them based only on the measure achieves higher-than-chance, but not perfect accuracy (73\% for Monkey A, 75\% for Monkey E, compared to baselines of 62\% and 61\% respectively). These results demonstrate there is consistency between the complexity score with our earlier classification of HO and OT groups ~\cite{tang}. It is unlikely that either method perfectly corresponds to the ground truth of neuronal tuning, but their consistency despite very different methodologies suggests that some such complex tuning does exist. The FKCNN-based metric introduced here has the additional benefit of efficiency: it does not require extensive testing with a specific set of thousands of complex pattern stimuli, rather just enough data to fit the two types of model.

\begin{figure}[htp]
  \centering
  \includegraphics[width=6cm]{images/final_a_disc_hist.png}
  \includegraphics[width=6cm]{images/E_gabor_diffs.png}
  \caption{{\bf Distributions of HO and OT neurons' scores:} These figures show histograms of the HO and OT neurons' scores, measured by the difference in performance between the sparse FKCNN and Gabor models, for each monkey. Note that since the classes have different numbers of neurons, the histograms have been normalized by the total number of neurons in each group. The two cell types are fairly separable by this metric\reviseold{, with a linear separability of 73\% for Monkey A and 75\% for Monkey E, compared to baselines of 62\% and 61\% respectively.}}
  \label{fig12}
\end{figure}

\begin{figure}[htp]
  \centering
  \includegraphics[width=13cm]{images/OT_HO.png}
  \caption{{\bf Classification of OT and HO neurons:} The figure shows the RF recovered by the CMPR algorithm for all 387 OT neurons and 232 HO neurons. The order is ranked by score, as described in equation~\ref{equ3}, in increasing order. Along the left margin, the proportion of each row with simple orientation kernels selected by CMPR is shown, clearly decreasing as the score increases, showing a correspondence between these metrics of HO tuning.}
  \label{fig10}
\end{figure}

\reviseold{We further explored the consistency of the original HO/OT classification with the top kernel extracted by the CMPR method.}
We split the OT and HO neurons up, rank each group by its FKCNN/GCNN score in increasing order, and show each neuron's primary CMPR kernel  (Fig~\ref{fig10}). The CMPR kernels can be classified into simple oriented and complex groups, and the proportion of oriented filters in each row is shown on the left side, clearly decreasing as the FKCNN models increasingly outperform the GCNN models, and lower on average for the HO neurons. This metric, an estimate of the probability of CMPR selecting a complex filter rather than a simple oriented one, is again consistent with the other two, and is qualitatively matched by the appearance of the filters themselves as well. The findings in \cite{tang} thus are doubly reinforced by our basis-selecting methods, both in the comparison between a \reviseold{complex-shape basis}  and a Gabor-based one for CNN modeling, and the specific nature of the filters chosen from the complex basis for CMPR modeling.

\section*{\reviseold{Discussion}}

\reviseold{In this paper, we explore the novel idea of incorporating a sparse code dictionary as a front-end to the projection pursuit regression and convolutional neural network approaches for neural response prediction and receptive field characterization of neurons in the primary visual cortex of macaque monkeys.
We show that the diverse \reviseold{complex-shaped} sparse codes learned from natural scenes based on the convolutional sparse coding principle can serve as an effective front-end to improve the predictive performance of both pursuit regression and CNN models.
Imposing these dictionary front-ends on RF models not only improves predictive performance, but also leads to faster convergence and higher data efficiency.} \reviseold{More importantly, it allows us to interpret the constituent components of a neuron's receptive field in terms of the dictionary or statistical priors learned from natural scenes using sparse coding theory.}

\reviseold{It is important to note that the set of dictionary elements that can yield good performance is not unique. Different subsets of kernels drawn from the learned sparse codes yield comparable results as long as there is sufficient diversity in the kernel elements. Indeed, adding three center-surround (Laplacian of Gaussian) kernels to the mix of Gabor-wavelet kernels is found to yield a significant improvement in performance. We also found subsets of diverse kernels in the first-layer of the classical CNN AlexNet \cite{alexnet}  trained for image classification performs worse than  the complex-shaped sparse codes, but better than standard Gabor wavelets (see Supplementary Information). This suggests that convolutional methods learned with either a sparse coding objective or an image classification objective yield a diverse set of kernels that  span well the manifold of natural scene priors used by V1.
Given that the set of fixed kernels that can yield good performance is not unique, the kernel features \reviseold{selected} should be more appropriately considered as the basis for spanning the space of subcomponents of the neurons' receptive fields, rather than the actual subcomponents themselves. Nevertheless, these subcomponents give us a glimpse of the potential complexity and diversity of the cells' subcomponents.}

Compared to standard CNNs for neural response predictions~\cite{zhang, 29}, \reviseold{our approach provides information on the relative importance of more structured, coherent intermediate-layer kernels for predicting a neuron's response. This allows us to relate the constituent components of the neurons' receptive fields to the sparse codes learned from natural scenes, aiding our understanding over the noisy components learned by a standard CNN.}
 \reviseold{
 Efficient coding theory  \cite{18, olshausenovercomplete, 28, CSC, gregor2010learning}  predicts the emergence of receptive fields with coherent  features.
It is an open question whether these coherently structured sparse codes provide intrinsic advantages for neural computation in the downstream visual areas.
Here, we show that models with these sparse code kernels do generalize better, as evident by their data efficiency and faster convergence in learning.}   \reviseold{Perhaps,} having \reviseold{sparse code} kernels that  are meaningful visual features might also facilitate learning and inference in the nervous system \reviseold{as well}.

\reviseold{ The pattern stimuli we used in this study bear strong similarity to the complex-shaped sparse codes used in our front-end.  This stimulus set was designed to evaluate the hypothesis that  higher-order features known to be encoded in V2 and beyond are also encoded by V1 neurons~\cite{tang}. While it provides a great variety of higher-order features in different orientations and positions, it also features abundant orientation stimuli, over 1600 oriented edges, lines and gratings, in different lengths, widths,  end-point positions, and many orientations in step as fine as 7.5 degrees, projected to different positions within the 1-degree diameter visual receptive fields of the neurons located within a 500 um $\times$ 500 um cortical patch. Thus, we believe there is sufficient representation of orientation stimuli for our regression approach to select the oriented kernels, if they are appropriate.  The observation that the top kernels selected by the CMPR method for OT neurons are more likely to be oriented kernels than those selected for HO neurons (Figure~\ref{fig10}) reassures us of this assumption.}

\reviseold{The high-order features such as corners, curvatures, junctions and crosses are given "equal opportunity" in our study \revise{ in that they are}  represented at a much higher rate than they are normally be observed in random natural scene patches. These higher-order features however were considered to be highly informative for vision \cite{attneave, Biederman91} and are widely assumed to be encoded by neurons in the extrastriate cortex~\cite{connorVanEssen, pasupathy2006}.  What is the consequence of this bias in the pattern stimuli tested on our findings regarding receptive field structure?  In the Supplementary Information, we provide results of these methods on \revise{the neurons of one monkey tested on} natural image stimuli. \revise{Natural images}  tend to have strong lower frequency components and feature predominantly uniform luminance regions or oriented edges in their random local patches.
Both the metrics based on FKCNN and the CMPR methods (Fig~\ref{fig12} and Fig~\ref{fig10} respectively), when trained on responses to natural images, fail to discriminate between the HO and OT neurons. These findings can be interpreted in two ways. One interpretation is that the pattern stimuli are important for detecting and discriminating the HO and OT neurons. \revise{ Natural images, when presented in 4X receptive field aperture as in our case, tended to evoke very sparse responses. When the amount of natural images is limited as in this case, the regression-based CMPR and CNN-based methods are not strong enough to detect the true nature of higher-order tuning in natural image responses. The second interpretation is that the particular higher-order tuning exhibited on the pattern stimuli, the preference for specific complex shapes like corners and crosses, is biased by the set of stimuli chosen, and may be different from the ground-truth higher-order tuning that produces the responses to natural images.
}}

\reviseold{Given these observations, our findings here alone do not constitute evidence for the diversity and complexity of V1 neurons. The arguments and evidence in support of selectivity and tuning of V1 neurons to higher-order features were put forth in \cite{tang}. Thus, our contribution here is primarily methodological. If we accept the premise that V1 HO neurons do exist, our approaches not only demonstrate some improvement in the CNN-based and PPR-based approaches in neural response modeling, but also provide useful metrics to detect  the HO neurons identified in the original study. We bolster that HO vs OT classification with two additional quantitative metrics, based on model properties instead of solely the neural response, and made possible by the use of fixed-kernel dictionaries in our models.  Recent iterative optimization approaches, such as the most exciting image (MEI) method~\cite{walker2019inception}, have also revealed complexity and diversity in V1 receptive fields in mice from natural images, even though earlier studies had suggested they have Gabor-like receptive fields.} \revise{This trend towards recognition of receptive field complexity is also present in work on the auditory cortex~\cite{auditory}.}

\reviseold{Our study shows that \reviseold{the diverse set of}   complex-shaped kernels provided by convolution sparse coding of natural scenes form an effective front-end for PPR-based and CNN-based approaches for characterizing V1 neurons' receptive fields in a data efficient manner.
\reviseold{This finding shows that an overcomplete sparse code dictionary learned from natural scenes can potentially be related to the constituent elements of V1 receptive fields. This is not to say that they are the same, just that the sparse codes span the space of those elements better than the other alternatives such as Gabor filters, Alexnet filters or data-driven filters learned from scratch.}
By making pairwise comparisons of the different methods, we gain some insight into the incremental benefit of each specific feature in the different methods.  The observation that CPPR is better than PPR indicates that max-pooling over convolution might provide an additional degree of freedom of matching flexibility, reminiscent of the classical deformable part models for object recognition~\cite{girshick}; that the FKCNN is better than CNN (and CMPR better than CPPR) suggests the features encoding natural scene priors provide a better embedding space; that the FKCNN is better than CMPR suggests that  compositional learning of multiple filters simultaneously in a neural network approach might further enhance matching flexibility to create a better embedding space; that the FKCNN is better than the GCNN suggests that a more diverse and complex set of kernels is important for modeling the responses of these V1 neurons, particularly those with higher order feature selectivity. \reviseold{Finally, the fact \revise{that} convolutional mechanisms are useful for learning the sparse codes and for providing flexibility in matching during model fitting in both the PPR-based and CNN-based methods suggests the potential relevance of convolution mechanisms for cortical information processing.}}

\section*{Acknowledgements}

This work was supported by an NSF grant CISE RI  1816568, and its Undergraduate Research Experience (URE) supplement, and the NIH grant R01 EY030226-01A1 to Lee, and National Natural Science Foundation of China (Grant No.31730109 and Grant No.U1909205) to Tang.
Ziniu Wu and Harold Rockwell were also supported by an interdisciplinary training fellowship in computational neuroscience at Carnegie Mellon sponsored by NIH NIDA grant 5T90 DA023426.

\nolinenumbers