\section{Introduction}
Weyl orbits \cite{Potter2014,Zhang2016} in topological semimetals are unique magnetic orbits that involve two Fermi arcs on opposite surfaces connected via the chiral Landau levels running through the bulk. Experiments studying the quantum oscillation and quantum Hall effect (QHE) related to such orbits have been carried out in recent years \cite{Zhang2019a,Zhang2017,Nakazawa2019,Nishihaya2019,Goyal2018,Galletti2019,Schumann2018,Lin2019,Galletti2018a,Zheng2017,Moll2016} and have sparked heated debates among research groups. For instance, as most of the experiments were conducted with the topological Dirac semimetal (DSM) candidate Cd$_3$As$_2$ \cite{Wang2013,Borisenko2014,Liu2014,He2014,Neupane2014}, some works attribute the QHE in those Cd$_3$As$_2$ films to the surface Dirac cones instead of the Weyl orbits \cite{Galletti2018a,Kealhofer2020}. This is because the surface states of Cd$_3$As$_2$ are shown to be two-dimensional (2D) Dirac cones by the angle-resolved photoemission measurements \cite{Yi2014a,Neupane2015,Roth2018a} instead of the open Fermi arcs as theoretically predicted \cite{Wang2013}. In order to distinguish the two mechanisms of QHE in Cd$_3$As$_2$, it is proposed as a key that the QHE based on Weyl orbits depends on the film thickness \cite{Zhang2019a,Goyal2018,Galletti2019,Zhang2016}, which was observed in a wedge-shaped film of Cd$_3$As$_2$ \cite{Zhang2019a}. However, a recent study suggests that the QHE in such wedge-shaped films may come from a completely different mechanism \cite{Cheng2020b}, and thus the existence of Weyl orbits remains elusive. As another example for the debates, some experiments \cite{Zhang2019a,Zhang2017} find that the energy levels of Weyl orbits are related to the quantum confinement subbands, which contradicts with some other works \cite{Nishihaya2019,Wang2017,Li2020b} distinguishing Weyl orbits with the confinement effect. Hence, further studies are still needed in order to understand the nature of both Weyl orbits and the surface states of Cd$_3$As$_2$.

Motivated by the connection between Weyl orbits and the quantum confinement effect, and by the progress in fabricating high-quality nanostructures of Cd$_3$As$_2$ \cite{Kealhofer2019,Galletti2018}, we investigate the QHE induced by confinement subbands of the chiral Landau levels in topological semimetal films. We study the QHE in Weyl semimetals (WSMs) in two cases of the confinement direction: with and without the surface states. While in the former, the energy levels agree with the results obtained from the semiclassical analysis as expected, we show that the chiral Landau subbands behave differently in the latter. By examining the difference between those two cases, we find a more general expression for the energy levels of Weyl orbits in a WSM. Furthermore, we find that in DSM films confined in their rotational symmetry axis, the quantum Hall and quantum spin Hall states coexist due to the spin conservation.

This paper is organized as follows. In Sec. \ref{sec2}, we introduce the models describing our WSM and DSM. In Sec. \ref{sec3}, we revisit the QHE based on Weyl orbits by considering a WSM film with chiral surface states on their boundaries. We then carry out the same calculation but in a film confined in a direction so that it has no surface states. In Sec. \ref{sec4}, we explain the results obtained in the preceding section and examine the QHE in a WSM confined in an arbitrary direction. We also study the QHE in a DSM confined in its rotational symmetry axis and consider a specific case of $(001)$ Cd$_3$As$_2$ film to show that the effect is experimentally observable. Finally, our results will be summarized in Sec. \ref{sec5}.
\section{Models}\label{sec2}
In this work, we describe a WSM by the minimal model \cite{Yang2011,Imura2011}
\begin{align}
	H_{\text{W}}(\mathbf{k}) = & t(\sin k_x\sigma_x + \sin k_y\sigma_y) \nonumber\\
	&+ [M + t'(\cos k_x + \cos k_y + \cos k_z)]\sigma_z\label{Eq: Weyl}
\end{align}
with $\sigma_i$ being the Pauli matrices. Notice that the dimensionless $k_i$ denote the product $q_ia_i$, where $q_i$ are the usual crystal momenta and $a_i$ lattice constants. This cubic lattice model is centrosymmetric and breaks the time-reversal symmetry. Hereafter, we use $t=t'=1$ and $M=-2.5$ so that the Hamiltonian gives a pair of Weyl nodes at $\mathbf{k}_W=(0,0,\pm\pi/3)$, as depicted in Fig.~\hyperref[fig:1]{1(a)}. If we confine this WSM in the $x$ direction, the topological surface states called Fermi arcs will appear. In our minimal model, the arcs are straight lines connecting the projections of the two bulk Weyl cones.

When the Weyl points are located close to the $\Gamma$-point, we can expand the Hamiltonian to obtain a continuum model
\begin{equation}
	H_{\text{W}}(\mathbf{k}) = t(k_x\sigma_x + k_y\sigma_y) + [m - t'(k_x^2+k_y^2+k_z^2)/2]\sigma_z,\label{Eq: Weylconti}
\end{equation}
where $m=M+3t'$.
We will use the lattice Hamiltonian for numerical calculations and its continuum counterpart for obtaining the analytical expressions.

\begin{figure}
	\includegraphics[width=\linewidth]{fig1
	\caption{\label{fig:1} (a) WSM/DSM with two Weyl/Dirac nodes (orange points) aligned along the $z$ axis. (b)/(c) A slab confined in the $x$/$z$ direction in a perpendicular magnetic field $\mathbf{B}$.}
\end{figure}

On the other hand, the DSM is described by an effective Hamiltonian of Cd$_3$As$_2$ and Na$_3$Bi \cite{Wang2013,Wang2012a}
\begin{equation}
	H_{\text{D}}(\mathbf{q}) = \tilde{\mathcal{E}}_0(\mathbf{q}) + \begin{pmatrix}
		\tilde{\mathcal{M}}(\mathbf{q}) & Aq_+ & 0 & 0 \\
		Aq_- & -\tilde{\mathcal{M}}(\mathbf{q}) & 0 & 0\\
		0 & 0 & \tilde{\mathcal{M}}(\mathbf{q}) & -Aq_-\\
		0 & 0 & -Ak_+ & -\tilde{\mathcal{M}}(\mathbf{q})
	\end{pmatrix} \label{Eq: Diracconti}
\end{equation}
with $q_\pm = q_x\pm iq_y$, $\tilde{\mathcal{E}}(\mathbf{q}) = C_0 + C_1q_z^2 + C_2(q_x^2+q_y^2)$, and $\tilde{\mathcal{M}}(\mathbf{q}) = M_0 + M_1q_z^2 + M_2(q_x^2+q_y^2)$. We transform this into a tetragonal lattice model using the substitution
\begin{equation}
	q_i\rightarrow \frac{1}{a_i}\sin k_i,\qquad q_i^2\rightarrow \frac{2}{a_i^2}(1-\cos k_i).\label{Eq: transform}
\end{equation}
Here, $a_x=a_y=a$ and $a_z=c$ are the tetragonal lattice constants. The Hamiltonian then becomes
\begin{align}
	H_{\text{D}}(\mathbf{k}) =& c_0 + c_1\cos k_z + c_2(\cos k_x + \cos k_y) \nonumber\\
	&+ t(\sin k_x\sigma_z\tau_x - \sin k_y\tau_y) \nonumber\\
	&+ [m_0 + m_1\cos k_z + m_2(\cos k_x + \cos k_y)]\tau_z,\label{Eq: Dirac}
\end{align}
where the Pauli matrices $\sigma_i$ and $\tau_i$ represent the spin and orbital degrees of freedom, respectively. These models of DSMs preserve both time-reversal and inversion symmetries, and have a rotational symmetry axis parallel to the $z$ axis. If we confine the DSM in any direction not parallel to $z$, each of its surfaces has two disconnected Fermi arcs with opposite spin and chirality. From now on, unless stated otherwise, we choose the parameters $c_0=c_1=c_2=0$, $t = m_1 = m_2 = 1$ and $m_0=-2.5$ so that the material has a pair of Dirac points at $\mathbf{k}_D=(0,0,\pm\pi/3)$ [Fig.~\hyperref[fig:1]{1(a)}].

We will examine the WSM and DSM in a slab geometry confined in the $i$ ($i=x$ or $z$) direction, as shown in Figs.~\hyperref[fig:1]{1(b)} and \hyperref[fig:1]{1(c)}, with periodic boundary conditions (PBCs) in the other two directions. The slabs are subjected to a uniform magnetic field $\mathbf{B}\|\hat{i}$, whose magnitude can be written in terms of the magnetic flux $\Phi$ threading a unit cell as $B=\Phi/(a_ja_k)=\Phi_0\phi/(a_ja_k)$. Here, $i\neq j\neq k$, $\Phi_0=h/e$ is the magnetic flux quantum, and $\phi=\Phi/\Phi_0$.
\section{Quantum Hall effect in Weyl semimetal films}\label{sec3}
\subsection{Films confined in the $x$ direction}
First, we consider the WSM confined in the $x$ direction [Fig.~\hyperref[fig:1]{1(b)}] with a magnetic field $\mathbf{B}\|\hat{x}$ represented by vector potential $\mathbf{A}=(0,0,Bay)$, where $y$ is dimensionless. The Hamiltonian of our WSM is then modified in accordance with the Peierls substitution as (Appendix \ref{App:Lattice})
\begin{align}
	\mathcal{H}_{\text{W}} = \sum_{x,y,k}\frac{1}{2}\Big\{&d^{\dagger}_{xyk}[M + t'\cos(k+2\pi\phi y)]\sigma_zd_{xyk}\nonumber\\
	&+ d^{\dagger}_{xyk}(t'\sigma_z-it\sigma_x)d_{(x+1)yk} \nonumber\\
	&+ d^{\dagger}_{xyk}(t'\sigma_z-it\sigma_y)d_{x(y+1)k}\Big\} + h.c.\label{Eq: wslabx}
\end{align}
with $d^{\dagger}_{xyk}$ ($d_{xyk}$) being creation (annihilation) operator, $k$ momentum in the $z$ direction. This Hamiltonian gives an energy spectrum [Fig.~\hyperref[fig:2]{2(a)}] that was shown in Ref. \cite{Igarashi2017} with a similar model to have the energy levels in agreement with those of Weyl orbits \cite{Potter2014,Zhang2016}
\begin{equation}
	\epsilon_n = \frac{\pi t}{k_a/(2\pi\phi) + N_x+1}(n+\gamma),\quad n=0,1,\ldots\label{Eq: Weyl orbit}
\end{equation}
Here, $k_a$ is the length of the surface Fermi arcs, $N_x$ the number of lattice sites in the $x$ direction, and $\gamma$ the phase offset that depends on the distance between the Weyl points \cite{Borchmann2017}.
Additionally, to confirm the quantum Hall state of this system, we compute the Chern numbers of some energy gaps formed by the discrete levels of Weyl orbits using the Streda formula \cite{Streda1982} (Appendix \ref{App:Chern}). The nonzero Chern numbers indicate the existence of a quantized Hall conductance and are equal to the number of chiral modes at the edges of our system. Moreover, such a quantum Hall state depends on the film thickness as the energy levels do [Eq.~(\ref{Eq: Weyl orbit})], which is regarded as a signature of the quantum Hall effect based on Weyl orbits \cite{Zhang2019a,Goyal2018,Galletti2019}.
\begin{figure*}
	\includegraphics[width=\linewidth]{fig2
	\caption{\label{fig:2} Energy spectra against the magnetic flux $\phi$ for (a) a WSM slab confined in the $x$ direction [Eq.~(\ref{Eq: wslabx})], and (b) a WSM bulk [Eq.~(\ref{Eq: wbulkx})] with $k_x=\pi \zeta/(N_x+1)$. The green lines represent semiclassical result determined by Eq.~(\ref{Eq: Weyl orbit}) for (a) $k_a = 2\pi/3$ and $\gamma=0.7$, (b) $k_a=0$ and $\gamma=1$. (c) The bulk Landau bands of a WSM confined in the $x$ direction and subjected to a magnetic flux $\phi=0.01$. The green dots indicate confinement subbands formed from the chiral Landau bands. The Chern numbers of some energy gaps are shown in (a) and (b). All the results are computed with $N_x=20$.}
\end{figure*}

As mentioned in some experiments \cite{Zhang2017,Zhang2019a}, the Weyl orbit levels are related to the quantum confinement effect in topological semimetal films. In order to illustrate this idea, we consider those levels [Eq.~(\ref{Eq: Weyl orbit})] in the limit where the Fermi arcs vanish, $k_a=0$, and consider only the tunneling process of the quasi-particles via the bulk chiral Landau levels. In this case, we see that the energy levels take the form $\epsilon_n'=\pi t(n+1)/(N_x+1)$ and do not change with the magnetic field. On the other hand, we know that in the bulk WSM the linear dispersion of the chiral Landau levels at the Weyl points is given by
\begin{equation}
	E_c=\pm tk_x.\label{Eq: linear chiral}
\end{equation}
A comparison between $\epsilon_n'$ and $E_c$ indicates that the momentum $k_x$ is quantized in a similar way to an infinite quantum well problem, i.e., $k_x=(n+1)\pi/(N_x+1)$. This indicates that the energy levels of Weyl orbits in the limit $k_a=0$ are confinement subbands stemming from the chiral Landau levels. To further substantiate this argument, we carry out some numerical calculations as follows. We impose an additional PBC on the Hamiltonian~(\ref{Eq: wslabx}) in the $x$ direction to eliminate the surface states and preserve only the bulk spectrum. The Hamiltonian is then given by
\begin{align}
	\mathcal{H}_{\text{W}} = &\sum_{k_x,y,k_z}\frac{1}{2}\boldsymbol{\Big(}d^{\dagger}_{k_xyk_z}\big\{t\sin k_x\sigma_x + [M + t'\cos k_x \nonumber\\
	&+ t'\cos(k_z+2\pi\phi y)]\sigma_z\big\}d_{k_xyk_z} \nonumber\\
	&+ d^{\dagger}_{k_xyk_z}(t'\sigma_z-it\sigma_y)d_{k_x(y+1)k_z}\boldsymbol{\Big)} + h.c..\label{Eq: wbulkx}
\end{align}
In order to obtain only the bulk spectrum of our WSM slab, i.e., to add the effect of quantum confinement, we apply the particle-in-a-box method (PiBM) (Appendix \ref{App:confinement}) by diagonalizing this Hamiltonian only at the momenta
\begin{equation}
	k_x = \frac{\pi \zeta}{N_x+1},\qquad\zeta=1,2\ldots,N_x.\label{Eq: quantization x}
\end{equation}
The spectrum is shown in Fig.~\hyperref[fig:2]{2(b)}, where we see that distinct energy gaps still exist even in the absence of surface states. The energy levels also agree well with the semiclassical result given by Eq.~(\ref{Eq: Weyl orbit}) for $k_a=0$, as denoted by the green lines. To find their origin, we show the Landau bands obtained from Hamiltonian~(\ref{Eq: wbulkx}) for a fixed magnetic flux [Fig.~\hyperref[fig:2]{2(c)}]. Then, by adding the confinement effect [Eq.~(\ref{Eq: quantization x})], each continuous Landau band becomes a set of discrete confinement subbands. We see that the energy levels in Fig.~\hyperref[fig:2]{2(b)} are actually the confinement subbands formed from the chiral Landau levels, as we have predicted. In the presence of boundaries, these chiral Landau subbands hybridize with the surface Fermi arcs and bend towards zero energy in the low field regime, giving rise to the Weyl orbit levels shown in Fig.~\hyperref[fig:2]{2(a)}. As a result, we can conclude that the quantization of Weyl orbits gives energy levels that are chiral Landau subbands hybridizing with the surface Fermi arcs.

From this interpretation of the Weyl orbit levels, we gain two new perspectives about the QHE in topological semimetal films. First, the QHE induced by Weyl orbits is intrinsically two-dimensional instead of three-dimensional as being claimed before \cite{Nishihaya2019,Wang2017,Li2020b} since the energy levels originate from the quantum confinement effect. Further evidence for this 2D nature is that the Hall resistance of our WSM film is a factor of the Klitzing constant $R_K=h/e^2$, in agreement with the experiments, whereas in a so-claimed 3D QHE induced by the charge-density wave \cite{Tang2019a}, the Hall resistance is much smaller than $R_K$. Second, if our material somehow has the bulk chiral Landau levels but with no open Fermi arcs on its surfaces, e.g. the arcs are combined into a closed Fermi loop \cite{Potter2014}, the gaps between the chiral Landau subbands still remain. Hence, a thickness-dependent QHE or quantum oscillation in relatively thick topological semimetal films is not conclusive evidence for observing either Weyl orbits or surface Fermi arcs, in contrast to the usual expectation \cite{Zhang2019a,Goyal2018,Galletti2019,Moll2016}. Moreover, since such a QHE can take place even without the surface Fermi arcs, a question then naturally arises; is the QHE observable in our WSM when it is confined in the $z$ direction, which has no nontrivial states on its boundaries?
\subsection{Films confined in the $z$ direction}\label{sec: 3B}
We now consider a WSM confined in the $z$ direction [Fig.~\hyperref[fig:1]{1(c)}] with a magnetic field $\mathbf{B}\|\hat{z}$ given by vector potential $\mathbf{A}=(0,Bax,0)$. Similar to the previous case, the bulk spectrum of our WSM decomposes into 1D Landau bands dispersing along the $z$ direction, as shown in Fig.~\hyperref[fig:3]{3(b)}, and the chiral level still evolves differently from other bands. If we introduce confinement in $z$, we expect that the chiral Landau subbands also form energy gaps distinct from others, and the QHE in relatively thick films will thus be observable.
\begin{figure}
	\includegraphics[width=\linewidth]{fig3}
	\caption{\label{fig:3} (a) Energy spectrum against the magnetic flux $\phi$ for a WSM slab confined in the $z$ direction. The green lines represent energy levels determined by Eq.~(\ref{Eq: level z}). The Chern numbers of some energy gaps are shown. (b) The bulk Landau bands of a WSM confined in the $z$ axis and subjected to a magnetic flux $\phi=0.01$. The green dots indicate confinement subbands formed from the chiral Landau band. All the results are computed with $N_z=20$.}
\end{figure}
The Hamiltonian of our system now reads
\begin{align}
	\mathcal{H}_{\text{W}} =&\sum_{x,k,z}\frac{1}{2}\boldsymbol{\Big(}d^{\dagger}_{xkz}\big\{[M + t'\cos(k+2\pi\phi x)]\sigma_z \nonumber\\
	&+ t\sin(k+2\pi\phi x)\sigma_y\big\}d_{xkz} + d^{\dagger}_{xkz}t'\sigma_zd_{xk(z+1)} \nonumber\\
	&+ d^{\dagger}_{xkz}(t'\sigma_z-it\sigma_x)d_{(x+1)kz}\boldsymbol{\Big)} + h.c.\label{Eq: HWz}
\end{align}
with $k$ being momentum in the $y$ direction. The energy spectrum of this Hamiltonian is shown in Fig.~\hyperref[fig:3]{3(a)}. Since this slab does not have topological surface states, we can also obtain its spectrum analytically by using the PiBM. We start with adding the effect of magnetic field $\mathbf{B}$ on the continuum model [Eq.~(\ref{Eq: Weylconti})], which is done by replacing the momenta $k_x$ and $k_y$ with ladder operators $l$, $l^{\dagger}$ as
\begin{equation}
	k_x\rightarrow \sqrt{\pi\phi}(l^{\dagger}+l),\quad k_y\rightarrow -i\sqrt{\pi\phi}(l^{\dagger}-l).\label{eq: ladder}
\end{equation}
The Hamiltonian then becomes
\begin{equation}
	H_{\text{W}}(k_z) = \begin{pmatrix}
		\mathcal{M}(k_z) & 2\sqrt{\pi\phi}tl\\ 2\sqrt{\pi\phi}tl^{\dagger} & -\mathcal{M}(k_z)
	\end{pmatrix},
\end{equation}
with $\mathcal{M}(k_z) = m - t'k_z^2/2 - 2\pi\phi t'\left(l^{\dagger}l + \dfrac{1}{2}\right)$. Then, using the trial wavefunctions $(\alpha_1\ket{\nu-1}, \alpha_2\ket{\nu})^T$ for $\nu=1,2,\ldots$ and $(0,\ket{0})^T$ for $\nu=0$, where $\nu$ is the band index, we can obtain the spectrum of $H_{\text{W}}(k_z)$ from the secular equations as
\begin{equation}
	\det\begin{vmatrix}
		K_{\nu k_z} + \pi t'\phi - E & 2\sqrt{\pi\phi}t\sqrt{\nu}\\
		2\sqrt{\pi\phi}t\sqrt{\nu} & -K_{\nu k_z} + \pi t'\phi - E
	\end{vmatrix} = 0,
\end{equation}
for $\nu=1,2,\ldots$, and
\begin{equation}
	-K_{(\nu=0)k_z} + \pi t'\phi - E = 0,\qquad\text{for }\nu=0.
\end{equation}
Here, $K_{\nu k_z} = m - \frac{t'}{2}k_z^2 - 2\pi t'\phi \nu$. The 1D Landau bands of the WSM are then given by
\begin{align}
	E_0(k_z) &= -m + \frac{t'}{2}k_z^2 + \pi t'\phi,\\
	E_{\nu}(k_z) &= \pm\sqrt{K^2_{\nu k_z} + 4\pi t^2\phi \nu} + \pi t'\phi.
\end{align}
Now, we can transform the dispersion of the zeroth level into a lattice version by substituting $k_z^2\rightarrow2(1-\cos k_z)$, which yields
\begin{equation}
	E_0(k_z) = -M - 2t' - t'\cos k_z + \pi t'\phi.\label{Eq: quadratic chiral}
\end{equation}
Finally, we take into account the effect of quantum confinement by employing the PiBM, i.e., replacing $k_z=\pi\zeta/(N_z+1)$ with $\zeta=1,2,\ldots,N_z$. The subbands of chiral Landau level read
\begin{equation}
	\epsilon_0(\zeta) = -M - 2t' - t'\cos\frac{\pi\zeta}{N_z+1} + \pi t'\phi,\label{Eq: level z}
\end{equation}
and they evolve linearly with respect to the magnetic flux $\phi$, as shown by the green lines in Fig.~\hyperref[fig:3]{3(a)}. On the other hand, the Landau bands with $\nu>0$ move away from $0$ and make the gaps between the chiral Landau subbands observable. These gaps also have nonzero Chern numbers, indicating the existence of the QHE.
\section{Discussion}\label{sec4}
Based on the spectrum in Fig.~\hyperref[fig:3]{3(a)}, we make two inferences:

First, we see that when our WSM is confined in the $z$ direction, the dependence of its chiral Landau subbands on $\phi$ deviates considerably from that of the slab perpendicular to the $x$ axis. In particular, if confinement is in the $x$ direction, and we neglect the surface Fermi arcs, the subbands stay constant as the field strength increases [Fig.~\hyperref[fig:2]{2(b)}], which agrees with the semiclassical equation. On the other hand, if the WSM is confined in $z$, the subbands depend linearly on the flux and are unevenly spaced, and thus behave differently from Eq.~(\ref{Eq: Weyl orbit}). From Eq.~(\ref{Eq: level z}), we know that such a difference stems from the second-order terms in $k$ of Hamiltonian~(\ref{Eq: Weylconti}) while Eq.~(\ref{Eq: Weyl orbit}) was obtained by using a linear dispersion of the chiral Landau levels \cite{Potter2014}. However, another problem then comes up: those quadratic terms in $k$ contribute substantially to the spectrum in Fig.~\hyperref[fig:3]{3(a)} but does not affect the one in Fig.~\hyperref[fig:2]{2(b)}. We investigate this problem and revisit the semiclassical Weyl orbits in subsection \ref{subsec: 1}.

Second, since a DSM can be regarded as a combination of two WSMs with opposite spin and chirality, we expect that a $(001)$ film of DSM will have a spectrum with nonzero spin Hall conductance. We demonstrate this idea in subsection \ref{subsec: 2}, and show that such a quantum spin Hall effect (QSHE) is observable in the DSM candidate Cd$_3$As$_2$.
\subsection{Dependence of the chiral Landau levels on the $k^2$-terms in $H_{\text{W}}(\mathbf{k})$}\label{subsec: 1}
We study how the $k^2$-terms in Hamiltonian~(\ref{Eq: Weylconti}) affect the dispersion of the zeroth Landau bands by finding their analytical expressions when the magnetic field $\mathbf{B}$ is applied along an arbitrary direction. We transform the vectors $(k_x,k_y,k_z)$ of the crystal frame into the $(k_1,k_2,k_3)$ of the magnetic frame ($\hat{k}_3\|\mathbf{B}$) using a 3D rotation matrix. For simplicity and without loss of generality, we assume that $\mathbf{B}\perp\hat{k}_y$ and rotate the vectors about the $y$ axis ($k_2\equiv k_y$) as
\begin{equation}
	\begin{pmatrix} k_1 \\ k_2 \\ k_3 \end{pmatrix} =
	\begin{pmatrix}
		\cos\theta & 0 & -\sin\theta\\
		0 & 1 & 0\\
		\sin\theta & 0 & \cos\theta
	\end{pmatrix}\begin{pmatrix}
		k_x \\ k_y \\ k_z
	\end{pmatrix}
\end{equation}
with $\theta$ being the angle between $k_3$ and $k_z$. The continuum Hamiltonian of our WSM then reads
\begin{equation}
	H_{\text{W}}(\mathbf{k}) = t(\cos\theta k_1 + \sin\theta k_3)\sigma_x + k_2\sigma_y + \mathcal{M}(\mathbf{k})\sigma_z ,\label{Eq: H rotate}
\end{equation}
where $\mathcal{M}(\mathbf{k}) = m - t'(k_1^2+k_2^2+k_3^2)/2$. The eigenvalues of this Hamiltonian are given by $E = \pm\sqrt{\mathcal{M}^2 + t^2[(\cos\theta k_1+\sin\theta k_3)^2+k_2^2]}$, and the Weyl points are located at $(-k_W\sin\theta,0,k_W\cos\theta)$ and $(k_W\sin\theta,0,-k_W\cos\theta)$, $k_W=\sqrt{2m/t'}$.
\subsubsection{Bulk chiral Landau levels}
To find the chiral Landau bands analytically, we consider the parameters $\{t,t',m\}$ that satisfy the constraint $t'k_W/t=\eta=\pm1$, which keeps the velocity at the Weyl points isotropic. In the vicinity of the first Weyl point $(-k_W\sin\theta,0,k_W\cos\theta)$, the Hamiltonian reads
\begin{align}
	&H_{\text{W}}(\mathbf{k}') = tk_1'(\cos\theta\sigma_x + \eta\sin\theta\sigma_z) + tk_2\sigma_y \nonumber\\
	&+ tk_3'(\sin\theta\sigma_x - \eta\cos\theta\sigma_z) - \frac{t'}{2}(k_1'^2+k_2^2+k_3'^2)\sigma_z,
\end{align}
where $k_1'=k_1+k_W\sin\theta$ and $k_3'=k_3-k_W\cos\theta$. When a magnetic field $\mathbf{B}$ is applied, the momenta $k_1$ and $k_2$ are quantized in terms of the ladder operators similar to Eq.~(\ref{eq: ladder}). We then have $H_{\text{W}}(k_z) = H_1(k_z) + H_2(k_z)$, where the first-order term is
\begin{align}
	&H_1(k_z) = tk_3'\begin{pmatrix} -\eta\cos\theta & \sin\theta \\ \sin\theta & \eta\cos\theta\end{pmatrix} \nonumber\\
	&+ t\sqrt{\pi\phi}\begin{pmatrix}
		\eta\sin\theta(l^{\dagger}+l) & \cos\theta(l^{\dagger}+l) - (l^{\dagger}-l) \\
		\cos\theta(l^{\dagger}+l) + (l^{\dagger}-l) & -\eta\sin\theta(l^{\dagger}+l)
	\end{pmatrix},
\end{align}
and the second-order one reads
\begin{equation}
	H_2(k_z) = -\frac{t'}{2}\left[2\pi\phi(2l^{\dagger}l + 1) + k_3'^2\right]\sigma_z.
\end{equation}
We now make an approximation by solving the chiral Landau level from $H_1$ and treating $H_2$ as a perturbation. With the eigenvector $\ket{L_{c+}} = [1/\sqrt{2(1+\cos\theta)}](\eta\sin\theta, 1+\cos\theta)^T\ket{0}$, the chiral Landau level at the first Weyl point is given by $E_+^{(0)}(k_3') = \eta tk_3'$. The first-order correction from the perturbation $H_2$ is
\begin{equation}
	E_+^{(1)}(k_3') = \braket{L_{c1}|H_2|L_{c1}} = t'\left(\pi\phi + \frac{k_3'^2}{2}\right)\cos\theta .
\end{equation}
The dispersion of the zeroth Landau level is then expressed as
\begin{eqnarray}
	E_+(k_3) &=& \eta t(k_3-k_W\cos\theta) \nonumber\\
	&&+ t'\left[\pi\phi + \frac{(k_3-k_W\cos\theta)^2}{2}\right]\cos\theta .\label{Eq: lls +}
\end{eqnarray}
The first-order correction from $H_2$ allows us to reproduce both Eq.~(\ref{Eq: linear chiral}) and Eq.~(\ref{Eq: quadratic chiral}), and hence we neglect all higher-order corrections

A similar calculation yields the chiral Landau level crossing the second Weyl point $(-\sin k_W,0,\cos k_W)$ as
\begin{eqnarray}
	E_-(k_3) &=& -\eta t(k_3+k_W\cos\theta) \nonumber\\
	&&+ t'\left[\pi\phi + \frac{(k_3+k_W\cos\theta)^2}{2}\right]\cos\theta .\label{Eq: lls -}
\end{eqnarray}
These expressions demonstrate that the contribution from $H_2$ to the chiral Landau levels vanishes when the magnetic field is perpendicular to the line connecting the two Weyl points, i.e., $\mathbf{B}\|\hat{x}$. Besides, the effect of $H_2$ becomes more significant when the magnetic field direction approaches the $z$ axis. To assess the reliability of these results, we compare Eqs.~(\ref{Eq: lls +}) and (\ref{Eq: lls -}) with the Landau bands obtained from the lattice model of Eq.~(\ref{Eq: H rotate}), as shown in Fig.~\hyperref[fig:4]{4}.
\begin{figure*}
	\includegraphics[width=\linewidth]{fig4}
	\caption{\label{fig:4} Landau bands of our WSM bulk in a magnetic field with (a) $\theta=60^{\circ}$, (b) $\theta=45^{\circ}$, and (c) $\theta=30^{\circ}$. The Landau bands (magenta lines) are obtained from the lattice model of Hamiltonian (\ref{Eq: H rotate}). The green lines represent the analytical results given by Eqs.~(\ref{Eq: lls +}) and (\ref{Eq: lls -}). The blue lines are also those equations but being modified to a lattice version by Eq.~(\ref{Eq: transform}). The green and blue lines approach each other as the Weyl points come close to the $\Gamma$-point.}
\end{figure*}
\subsubsection{QHE in a WSM slab confined in an arbitrary direction}
After getting a more general expression for the chiral Landau bands to explain the difference between Figs.~\hyperref[fig:2]{2(b)} and \hyperref[fig:3]{3(a)}, it is interesting to seek an expression for the Weyl orbit levels taking into account the effect of the $k^2$-terms. The Onsager - Bohr - Sommerfeld quantization for a classical orbit reads
\begin{equation}
	\oint\mathbf{p}\cdot d\mathbf{r} = 2\pi\hbar(n+\gamma).
\end{equation}
After some calculations, we get
\begin{align}
	\epsilon_n = \frac{B}{t'\cos\theta}\left[A + B - \sqrt{2AB + B^2 + t^2 - 2\pi\phi t'^2\cos^2\theta}\right]\label{Eq: minh Weyl orbits}
\end{align}
with
\begin{equation}
	A = t'\left[\frac{\pi(n+\gamma)}{L_z'}\cos\theta + k_W\sin^2\theta\right]\text{ and } B = 2\pi\phi t\frac{L_z'}{k_a}.
\end{equation}
Here, $L_z'=(N_z+1)/\cos(\Theta-\theta)$ is dimensionless, $k_a = 2k_W\sin\Theta$, $\theta$ and $\Theta$ determine the directions of magnetic field and confinement, respectively. In the limits $\Theta=\theta\rightarrow 90^{\circ}$ and $\Theta=\theta\rightarrow 0^{\circ}$, this equation reproduces Eq.~(\ref{Eq: Weyl orbit}) and Eq.~(\ref{Eq: level z}).
\begin{figure*}
	\includegraphics[width=\linewidth]{fig5}
	\caption{\label{fig:5} Energy spectra against the magnetic flux $\phi$ for a WSM slab confined in the $k_3$ direction determined by (a) $\Theta=60^{\circ}$, (b) $\Theta=45^{\circ}$, and (c) $\Theta=30^{\circ}$. The green lines represent semiclassical result determined by Eq.~(\ref{Eq: minh Weyl orbits}) with $\gamma=0.7$. The Chern numbers of some energy gaps are shown. All the results are computed with a thickness corresponding to $N_3=20$.}
\end{figure*}

We now compare Eq.~(\ref{Eq: minh Weyl orbits}) with the results obtained from the lattice model of Hamiltonian (\ref{Eq: H rotate}). A WSM confined in the $k_3$ direction and subjected to a perpendicular magnetic field is described by
\begin{align}
	\mathcal{H}_{\text{W}} = &\sum_{x_1,k,x_3}\frac{1}{2}\boldsymbol{\Big(}a^{\dagger}_{x_1kx_3}\Big\{ t\sin (k+2\pi\phi x_1)\sigma_y \nonumber\\
	&+ \big[M + t'\cos(k+2\pi\phi x_1)\big]\sigma_z\Big\} a_{x_1kx_3}\nonumber\\
	&+ a^{\dagger}_{x_1kx_3}\left(t'\sigma_z - it\cos\theta\sigma_x\right)a_{(x_1+1)kx_3}\nonumber\\
	&+ a^{\dagger}_{x_1kx_3}\left(t'\sigma_z - it\sin\theta\sigma_x\right)a_{x_1k(x_3+1)}\boldsymbol{\Big)} + h.c. .
\end{align}
The energy spectra of this WSM are shown in Fig.~\hyperref[fig:5]{5} for different tilting angles, which shows a transition from Fig.~\hyperref[fig:2]{2(a)} to Fig.~\hyperref[fig:3]{3(a)}, and agree well with Eq.~(\ref{Eq: minh Weyl orbits}). When the confinement direction deviates from the $x$ axis, both the spacing and the dependence on magnetic flux $\phi$ of the chiral Landau subbands change. Interestingly, we can see how the subbands evolve by computing the Chern numbers of the gaps between them.

In DSMs, we assume that the physics is somewhat similar. As the experiments about Weyl orbits are often conducted in $(112)$ films of Cd$_3$As$_2$, this result may give better explanations to those magneto-transport studies, e.g. the high values of the Landau indices in the quantum Hall measurements.
\subsection{Quantum spin Hall effect in topological Dirac semimetal thin films}\label{subsec: 2}
Finally, we show how the QHE and QSHE induced by chiral Landau subbands take place in the DSM films. For our DSM [Eq.~(\ref{Eq: Dirac})], a magnetic field $\mathbf{B}\|\hat{z}$ preserves the spin-$z$ component ($S_z$) and lifts the spin degeneracy even in the absence of Zeeman coupling. A DSM film grown along its rotational symmetry axis ($z$) and subjected to $\mathbf{B}$ is described by
\begin{align}
	\mathcal{H}_{\text{D}} = \sum_{x,k,z}\frac{1}{2}\Big\{&d^{\dagger}_{xkz}\big[m_0\tau_z + m_2\cos(k+2\pi\phi x)\tau_z \nonumber\\
	&- t\sin(k+2\pi\phi x)\tau_y + \lambda_z\phi\sigma_z\tau_g\big]d_{xkz} \nonumber\\
	&+ d^{\dagger}_{xkz}(m_2\tau_z - it\sigma_z\tau_x)d_{(x+1)kz}\nonumber\\
	& + d^{\dagger}_{xkz}m_1\tau_zd_{xk(z+1)}\Big\} + h.c.,
\end{align}
where $\lambda_z=\dfrac{h^2}{8\pi m_ea^2}$ and
\begin{equation}
	\tau_g = \begin{pmatrix}
		g_s & 0 \\ 0 & g_p
	\end{pmatrix},
\end{equation}
$g_s$ and $g_p$ are the effective $g$-factors for $s$ and $p$ orbitals, respectively. For simplicity, we choose $\lambda_zg_s=2$ and $\lambda_zg_p=1$. Solving the eigenvalues of $\mathcal{H}_D$ numerically gives an energy spectrum as shown in Fig.~\hyperref[fig:6]{6(a)}, which is composed of two sets of chiral Landau subbands with opposite chirality and spin polarization. They evolve in different directions with respect to the magnetic field and thus cross each other to form the energy gaps that can be seen as an overlap of two separate gaps. This spectrum is also explained by the quantum confinement picture as shown in Fig.~\hyperref[fig:6]{6(b)}. According to this figure, we should notice that the two sets of subbands cross if $\pi/(N_z+1)<|\mathbf{k}_D|$, i.e.,
\begin{equation}
	L_z>\frac{\pi}{|\mathbf{q}_D|}.\label{Eq: thickness}
\end{equation}
To find how the chiral Landau subbands depend on the magnetic field, we also follow the same calculation as presented in Sec.~\ref{sec: 3B} and get an analytical expression of those subbands
\begin{subequations}
	\begin{eqnarray}
		\epsilon_{0\uparrow}(q_z) = C_0 &+& M_0 + (C_1+M_1)q_z^2 \nonumber\\
		&+& (C_2+M_2)\frac{eB}{\hbar} + \frac{\mu_Bg_sB}{2},
	\end{eqnarray}
	\begin{eqnarray}
		\epsilon_{0\downarrow}(q_z) = C_0 &-& M_0 + (C_1-M_1)q_z^2 \nonumber\\
		&+& (C_2-M_2)\frac{eB}{\hbar} - \frac{\mu_Bg_pB}{2}.
	\end{eqnarray}
\end{subequations}
To make a comparison with the lattice model spectrum [Fig.~\hyperref[fig:6]{6(a)}], we let $C_i=0$ and transform the equations using Eq.~(\ref{Eq: transform}), which gives
\begin{subequations}
	\label{eq:whole}
	\begin{equation}
		\epsilon_{0\uparrow}(\zeta) = m_0 + 2m_2 + m_1\cos\frac{\pi\zeta}{N_z+1} - \pi m_2\phi + \lambda_zg_s\phi,\label{Eq: Dirac subbands up}
	\end{equation}
	\begin{equation}
		\epsilon_{0\downarrow}(\zeta) = -m_0 - 2m_2 - m_1\cos\frac{\pi\zeta}{N_z+1} + \pi m_2\phi - \lambda_zg_p\phi.\label{Eq: Dirac subbands down}
	\end{equation}
\end{subequations}
These expressions show that the presence of Zeeman interaction still keeps the chiral Landau subbands evolving linearly with the magnetic flux but changes their slopes. Specifically, the $g_s$-factor modifies the slope of the spin-up chiral Landau subbands whereas $g_p$ affects the spin-down ones, which reflects the band inversion of DSMs.
\begin{figure}
	\includegraphics[width=\linewidth]{fig6}
	\caption{\label{fig:6} (a) Energy spectrum against the magnetic flux $\phi$ for a DSM slab. The red lines denote spin-up chiral Landau subbands given by Eq.~(\ref{Eq: Dirac subbands up}). The blue lines denote spin-down chiral Landau subbands given by Eq.~(\ref{Eq: Dirac subbands down}). The Chern numbers of some energy gaps are shown. (b) Bulk Landau bands of the DSM confined in the $z$ axis and subjected to a magnetic flux $\phi=0.01$. The red and blue dots indicate confinement subbands formed from the chiral Landau bands. All the results are computed with $N_z=20$.}
\end{figure}

Due to  our choice of parameters, the spin-up chiral Landau subbands are hole-like, and the spin-down ones are electron-like. As a result, when our material has an additional boundary, the spin-up subbands are bent downward and give right-handed edge states while those with spin-down disperse upward giving left-handed edge states. The combination of these edge modes may give rise to the coexistence of chiral and helical edge states. For example, we consider an energy gap of Chern number $\mathcal{C}=1$, which is a combination of $\mathcal{C}_{\downarrow}=7$ and $\mathcal{C}_{\uparrow}=-6$ gaps. When the Fermi level lies within this gap, the spin Hall conductance is quantized as \cite{Murakami2006,Yang2006}
\begin{equation}
	\sigma_s = -C_s\frac{e}{2\pi}
\end{equation}
with $C_s=(C_{\uparrow}-C_{\downarrow})/2=-6.5$. At the boundary, energy levels are bent and form $7$ left-handed spin-down and $6$ right-handed spin-up edge states, which would result in a total of $1$ chiral edge state and $6$ helical edge states. Hence, both quantum Hall and quantum spin Hall phases coexist in our system whose time-reversal symmetry is broken by the magnetic field. In this case, the helical edge states exist because $S_z$ is approximately conserved \cite{Zhang2014}. This scenario is well-known in graphene \cite{Abanin2006,Fertig2006,Young2014}, where the helical edge states are protected by an additional symmetry instead of the time-reversal symmetry as in topological insulators.
\subsubsection*{Estimation for Cd$_3$As$_2$}
Since a quantized Hall conductance has recently been observed in $(001)$ films of Cd$_3$As$_2$ \cite{Kealhofer2020} under a strong magnetic field, we roughly estimate whether the QHE and QSHE induced by the chiral Landau subbands can be observed in such films. Depending on the growth condition, Cd$_3$As$_2$ can be either centrosymmetric or noncentrosymmetric \cite{Sankar2015b}. Here, we consider a centrosymmetric Cd$_3$As$_2$, which is more popular and can be described by Hamiltonian~(\ref{Eq: Diracconti}) with parameters obtained by fitting with the \textit{ab initio} calculation \cite{Cano2017}. Additionally, as the distance between the Dirac points of Cd$_3$As$_2$ obtained by the Landau level spectroscopy measurements \cite{Akrap2016,Jeon2014,Desrat2018,Krizman2019} is one order smaller than the \textit{ab initio} calculations \cite{Wang2013,Ali2014a,Conte2017a}, we also consider the hyperbolic model proposed by Ref. \cite{Jeon2014}. A detailed procedure is presented in Appendix \ref{App:quantization}. The spectra computed for these two models are shown in Fig.~\hyperref[fig:7]{7}, where, in both cases, we use $g_s=18.6$ \cite{Jeon2014} and choose $g_p=10$ for simplicity.
\begin{figure}
	\includegraphics[width=\linewidth]{fig7}
	\caption{\label{fig:7}Energy spectra against the magnetic flux $\phi$ for a $(001)$ film of Cd$_3$As$_2$ described by two different models: (a) Hamiltonian (\ref{Eq: Diracconti}) with parameters given by Ref. \cite{Cano2017}, the film thickness is $50$ nm, and (b) the hyperbolic model proposed by Ref. \cite{Jeon2014}, the film thickness is $120$ nm. Some Chern numbers are shown in both figures by comparing with Fig.~\hyperref[fig:6]{6(a)}.}
\end{figure}

In general, both spectra have three types of gaps formed by the chiral Landau subbands: those with nonzero spin Chern number, those with nonzero Chern number and zero spin Chern number, and one with vanishing Chern and spin Chern numbers. We are interested in the gaps with nonzero spin Chern number, which do not exist if the film is too thin but is not observable if the film is too thick. If we assume an energy gap $\Delta\epsilon\ge10$~K $\sim 1$ meV to be observable, the majority of gaps in Fig.~\hyperref[fig:7]{7(a)} and all gaps in Fig.~\hyperref[fig:7]{7(b)} satisfy the condition. Hence, the QHE induced by the chiral Landau subbands can be observed in relatively thick films compared to the conventional quantum wells ($5$-$20$ nm).

If we use this simple model [Fig.~\hyperref[fig:7]{7(a)}] to explain the QHE observed in Ref. \cite{Kealhofer2020}, it may capture two features of the experiment, i.e., the absence of $\nu=1$ plateau and the high resistance at strong magnetic field. Nevertheless, just like the explanation given in that reference, it is also unable to clarify the difference in activation energy between the states at even and odd filling factors. On the other hand, if the actual distance between the Dirac nodes agrees with the Landau level spectroscopy measurements, the QHE induced by chiral Landau subbands may not be observable in such a quantum Hall experiment since the film thickness does not satisfy Eq.~(\ref{Eq: thickness}).
\section{Conclusion}\label{sec5}
Our study gives a generic and simplified picture of the QHE induced by confinement subbands stemming from the chiral Landau levels in topological semimetal films. Using a minimal model of WSM, we have demonstrated that the energy levels of Weyl orbits originate from the confinement subbands of the chiral Landau levels that hybridize with the surface Fermi arcs. We have then studied how the $k^2$-terms in a WSM Hamiltonian affect the evolution of the chiral Landau bands with respect to the magnetic field. We have found a general expression for the Weyl orbit levels and shown how the QHE takes place in a WSM film confined in an arbitrary direction. Furthermore, when examining a DSM confined in its rotational symmetry axis, we may have not only explained the QHE recently observed but also predicted the coexistence of both quantum Hall and quantum spin Hall states in such a system.

\begin{acknowledgements}
	This work was supported by the Japan Society for the Promotion of Science KAKENHI (Grant Nos.
	JP19K14607 and JP20H01830) and by CREST, Japan Science and Technology Agency (Grant No. JPMJCR18T2).
\end{acknowledgements}


\section{Introduction}

The Zika virus is a flavivirus that upon infection in humans causes an illness, known as Zika fever,
identified commonly with macular or papular rash, mild fever and arthritis
\cite{Brasil2016,WHOurl2016}. It is mainly a vector-borne disease carried by
the genus \textit{Aedes} of mosquitoes \cite{WHOurl2016,Fernandes2016}, while in a lesser amount
it is also transmitted via sexual interaction \cite{Petersen2016,Coelho2016},
and contamination by blood transfusion is under investigation \cite{Motta2016}.
The Zika virus was first isolated in primates from the Zika forest in Uganda in $1947$ \cite{Dick1952}.
Evidences of the virus in humans were found in Nigeria in $1968$ \cite{Moore1975}.
An epidemic occurred in $2007$ on Micronesia \cite{Duffy2009}, followed
by multiple outbreaks on several Pacific Islands between $2013$
and $2014$ \cite{Cao2014,Tognarelli2016}. The first Zika virus autochthonous case
in Brazil was reported around April, $2015$ \cite{SVS2017}, and nearly $30{,}000$ cases of infection
were already notified by January $30$, $2016$ \cite{Faria2016},
along with the Pan American Health Organization being informed in the same month about
locally-transmitted cases on numerous continental and island territories of America \cite{Hennessey2016}.
The Brazilian Ministry of Health registered $215{,}319$ probable cases of Zika fever
($130{,}701$ of which were confirmed) until the $52$th epidemiological
week (EW) of $2016$ \cite{SVS2017}. The Zika epidemic has been causing concern in
the international medical community, health authorities and population, specially due to an association
between the Zika virus and other diseases, such as newborn microcephaly \cite{Petersen2016,Valentine2016}
and  Guillain-Barré syndrome \cite{Santos2016}, whose correlation to the
Zika virus was considered by the World Health Organization
a \textit{``scientific consensus''} \cite{WHOApril2016}.

In this epidemic scenario, the development of control and prevention strategies for the mass
infection is a critical issue. A mathematical model capable of \fanswer{providing} a description
of the infected people throughout an outbreak is \fanswer{a useful} tool that can be employed to identify
effective and vulnerable aspects on disease control programs \cite{Naheed2014,Huppert2013,Lizarralde2017}.
Furthermore, for an epidemic model to be truly useful it must undergo a judicious process of \emph{validation}
\cite{Oberkampf2010,Cunhajr2017}, which consists in comparing model predictions with real
data in order to evaluate if they are realistic. In general, the first predictions of a model
do not agree with the observations, possibly due to inadequacies in the model hypotheses or
because of a poor choice for the model nominal parameters. The first case invalidates the model,
but the second can be amended through a procedure known as \emph{model calibration}, where a
set of parameters that promote a good agreement between predictions and observations is sought.

This work is one of the results in a rigorous ongoing process of identification and
validation of representative models to describe Zika virus outbreaks in Brazil
\cite{Cunhajr_ccis2016,Cunhajr_cnmac2017_1}.
For this purpose, a SEIR-SEI mathematical model is adapted to the Brazilian scenario.
This specific SEIR-SEI description has been successfully used before for the outbreaks in
Micronesia \cite{Funk2016} and French Polynesia \cite{Kucharski2016}. Some
assumptions were also based on similar studies performed over SEIR dynamical systems \cite{Wang2014,Wanga2017,Safi2013,Cai2017}. The nominal values of the model parameters
belong to characteristics of the Zika infection
and its vector, quantitatively estimated in the literature or published by health
organizations. Predictions are obtained from numerical simulation and further heuristic
manipulation, followed by a comparison to real data of the outbreak as
an initial effort to validate the model. In sequence, a rigorous process of model calibration
is performed through the formulation and solution of an inverse problem.

The rest of this paper is organized as follows. In section~\ref{EpidemicModel},
the mathematical model is described and estimation of nominal values for
the model parameters is discussed. In section~\ref{model_calibration},
the forward and inverse problems are formulated and solved, where results are detailed
and a subsequent comparison between model predictions
and experimental data is made to calibrate the model.
Finally, in section~\ref{concl_remaks}, the main contributions of this work
are emphasized and a path for future works is suggested.

\section{Epidemic model for Zika virus dynamics}
\label{EpidemicModel}

\subsection{Model hypotheses}
\label{ModHyp}

This work utilizes a variant of the Ross-Macdonald model \cite{Smith2012}
for epidemic predictions, separating the populations into a SEIR-SEI framework
(susceptible, exposed, infectious, recovered) \cite{Brauer2012,Brauer2008,Martcheva2015}.
Each category represents the health condition of an individual inside such group at time $t$,
with respect to the considered infection. The susceptible group, denoted by $S(t)$, represents
those who are uncontaminated and are able to become infected. The exposed portion of the population,
$E(t)$, comprehends anyone that is carrying the pathogen but is still incapable of transmitting the disease.
While the infectious individuals, $I(t)$, can spread the pathogen and may display symptoms associated
with the illness. Finally, the recovered group, $R(t)$, contains those who are no longer infected. The
populations of humans and vectors are segmented into the aforesaid classes (excepting the group of
recovered vectors), as Figure~\ref{model_schematic} depicts
schematically with the accordingly subscripts. $S_h$, $E_h$, $I_h$ and $R_h$ amass the number of
people at each stage of the model description, and $S_v$, $E_v$, $I_v$ signifies proportion of vectors
($0 \leq S_v, \, E_v, \, I_v \leq 1$ and $ S_v + E_v + I_v = 1$).

\begin{figure}[h!]
\centering
\includegraphics[height=5.7cm]{Diagram_SEIRmodel.pdf}
\vspace{-2mm}
\caption{Schematic representation of the SEIR--SEI model for the Zika virus outbreak description.}
\label{model_schematic}
\end{figure}

\fanswer{Demographic} changes in the number of humans are not considered because the timescale of infection
is much faster than the timescale of birth and deaths for the latter to significantly alter the
development of the disease in the Brazilian context (supplementary material B provides tools to test this
assertion). The total vector population is maintained constant during the analysis,
although variations of the proportion of vectors on the particular
compartments are introduced via birth and death rates. The vector in
question is regarded as a hypothetical mosquito apt to being infected
or infectious throughout all its lifetime --- which means the model
accounts only for the adult stage of their life cycle --- and also unable
to recover.

The time elapsed while an individual is on the aforementioned
exposed group is known as the latent period of an organism and, in this work,
is adopted as equivalent to the commonly called incubation period (the time elapsed
between being infected and exhibiting symptoms), since data is extremely
sparse on the latter for humans \cite{Chan2012} (namely, the intrinsic incubation period).
Both terms are used interchangeably hereafter, and the concepts do not differ
on the mosquitoes case (the extrinsic incubation period) \cite{Lessler2016}.
In addition, all the members of the susceptible group are treated as equally capable
of being infected and the recovered ones as completely immunized.

\subsection{Model equations}
\label{ModEq}

The evolution of individuals through the SEIR-SEI groups is governed by the following
(nonlinear) autonomous system of ordinary differential equations
\begin{equation}
\begin{aligned}
\dod{S_h}{t} &= - \beta_h \, S_h \, \frac{I_v}{N_v} \:,
			& \dod{S_v}{t} &= \delta\,N_v - \beta_v \, S_v \, \frac{I_h}{N_h} - \delta \, S_v \:, \\[1em]
\dod{E_h}{t} &= {\beta_h} \, S_h \, \frac{I_v}{N_v} - \alpha_h \, E_h \:,
			& \dod{E_v}{t} &= \beta_v \, S_v \, \frac{I_h}{N_h} -(\alpha_v + \delta) \, E_v \:,\\[1em]
\dod{I_h}{t} &= \alpha_h \, E_h - \gamma \, I_h \:,
			& \dod{I_v}{t} &= \alpha_v \, E_v - \delta \, I_v \:, \\[1em]
\dod{R_h}{t} &= \gamma \, I_h \:,
			& \dod{C}{t} &= \alpha_h \, E_h \:,
\end{aligned}
\label{model_eqs}
\end{equation}
\noindent
where $N$ represents the total population and $1/\!\alpha$ the disease's incubation period
(each with the corresponding subscript of $h$ for human's and $v$ for vector's),
$1/\!\delta$ means the vector lifespan, $1/\!\gamma$ is the human infection period --- which is defined
in this work as the interval of time that a human is infectious --- and $\beta$ identifies the transmission rate,
specifically $\beta_h$ is the mosquito-to-human rate and $\beta_v$ the human-to-mosquito rate.

The transmission terms $\beta_h\,S_h\,I_v/N_v$ and $\beta_v\,S_v\,I_h/N_h$ are composed by a
number of susceptible individuals ($S_h,\,S_v$, respectively), a transmission rate
($\beta_h,\,\beta_v$) and the probability of the contact being made with an infectious member
of the other population ($I_v/N_v,\,I_h/N_h$). Both transmissions terms come from the assumption
that the rate of contacts is constant, which characterizes a frequency-dependent transmission
\cite{Begon2002}.

The compartmentalization hypothesis requires setting the variables at the initial time of the
analysis, $t_0$, such that their sum equals the total population in each case, e.g, $S_v(t_0)+E_v(t_0)+
I_v(t_0) = N_v = 1$. \fanswer{By so doing}, for all $t>t_0$, $S_h(t)$, $E_h(t)$, $I_h(t)$ and $R_h(t)$ will
always add up to $N_h$; and $N_v(t)$ will be the equilibrium solution of the initial value problem
\begin{equation}
\begin{aligned}
\dod{N_v}{t} &= \delta\,(1 - N_v(t)) \:, & N_v(t_0) = 1 \:,
\end{aligned}
\end{equation}
\noindent
namely, $N_v(t) = 1$. This consideration allows the simplification of treating the total vector population
as a parameter $N_v$, instead of a variable, since it stays constant throughout the analysis.
If one wishes to treat the vector population as a variable, a recruitment parameter would need to
be added in place of $\delta\, N_v$ as well as another differential equation to account for the changes in
$N_v(t)$.

The $\dif C/\!\dif t$ equation allows evaluation of the cumulative number $C(t)$ of infectious people
until the time $t$, that is, the amount of humans so far that contracted the disease and have
passed through or are in the infectious group at the given time.

Additionally, a set of $M=52$ points to represent the number of new infectious cases of Zika fever
at each week is defined as follows:
\begin{equation}
\mathcal{N}_{w} = C_{w} - C_{w-1} \:, \quad \mathcal{N}_{1} = C_1, ~~w=2, \cdots, 52 \:,
\label{NC_eq}
\end{equation}
\noindent
where $C_{w}$ is the cumulative number of infectious humans
in the \emph{w}-th EW.

Figure~\ref{data} organizes the data of cumulative number of infectious and new cases per
week provided by the Brazilian Ministry of Health \cite{Data2016} (supplementary material A) for 2016,
where \sanswer{the evolution of the infection can be seen}.
\begin{figure}[h]
\centering
\subfigure[Cumulative number]{\includegraphics[scale=0.32]{data_C.pdf}}
\subfigure[New cases]{\includegraphics[scale=0.32]{data_NewCases.pdf}}\\
\vspace{-3mm}
\caption{Outbreak data.}
\label{data}
\end{figure}
The expected behavior for the model is that it generates a seemingly year long outbreak for 2016
that matches $C(t)$ and $\mathcal{N}_{w}$ to the data of Figure~\ref{data}. The end of this outbreak is
marked by the exposed and infectious portions of vectors reaching very low values, effectively killing off
the contamination process to the point that $S_h(t)$ stabilizes.

\subsection{Nominal parameters}
\label{NomParam}

The preliminary values for the parameters of the set of Eqs.~(\ref{model_eqs}) come from the
related literature concerning the Zika infection,  vector-borne epidemic models,
the \textit{Aedes aegypti} mosquito (which is the main vector for
Zika, Dengue and Yellow fever in Brazil) and publications provided by
health organizations and government agencies.
The Brazilian Institute of Geography and Statistics
reports that Brazil had approximately $N_h = 206 \times 10^6$ people by July, $2016$ \cite{DOU},
and $N_v$ is set $1$ to entail an entire vector population.
The adopted extrinsic incubation period is  $1/\!\alpha_v = 9.1$ days \cite{Funk2016}.
This value agrees with other statistical confidence intervals (CI)
that are presented for the parameter in another works ($95$\% CI: $7.3$--$9.3$ days \cite{Ferguson2016})
and is close to the numbers suggested by experimental studies for the time necessary for the virus
to reach the mosquito's saliva after an infectious blood meal ($5$ \cite{Wong2013,Li2012}
and $7$ days \cite{Chouin-Carneiro2016}).
A systematic review and pooled analysis of the literature and
case studies available in \cite{Ott2016} estimates that the median intrinsic incubation period is
$5.9$ days ($95$\% CI: $4.4$–-$7.6$). This values is selected for $1/\!\alpha_h$ in this work and is
compatible with the range of $3$--$12$ days recommended by multiple sources
\cite{Valentine2016,Ioos2014,ECDC2015},
also formerly used in previous studies \cite{Villela2017}.
The aforementioned literature
analysis in \cite{Ott2016} also concludes that $9.9$ ($95$\% CI: $6.9$-–$21.4$) days is the mean time
until an infected has no detectable virus in blood. Considering the assumption that the infectiousness
in Zika infection ends $1.5$--$2$ days before the virus becomes
undetectable \cite{Funk2016,Ferguson2016}, the chosen value for the human infectious
period is $1/\!\gamma = 9.9-2 = 7.9$ days.
As for the vector lifespan $1/\!\delta$,
``\textit{the adult stage of the mosquito is considered to last an average of eleven days in
the urban environment}" according to \cite{Otero2006}. This is the assumed value for the parameter
in this work, which is also consistent with the usual life expectancy for
the mosquito in Rio de Janeiro, Brazil
\cite{Freitas2007}, and comes close to the average of $2$--$3$ weeks considered in
biological studies about the species \cite{Nelson1986} and by the Centers for Disease Control and
Prevention \cite{CDCurl}.
Lastly, the time between a
mosquito being infected and it infecting a human, $1/\!\beta_h$, and the time between a human
infection and a mosquito taking an infectious blood meal, $1/\!\beta_v$, is estimated
in \cite{Ferguson2016} as an average of $11.3$ days ($95$\% CI: $8.0$--$16.3$)
and $8.6$ days ($95$\% CI: $6.2$-$11.6$), respectively.

\section{Calibration of the epidemic model}
\label{model_calibration}

\subsection{Forward Problem}
\label{forward_problem}

The epidemic model of section~\ref{EpidemicModel}, supplemented by an
appropriate set of initial conditions, is a continuous-time dynamical system
\begin{equation}
	\dot{\mathbf{x}}(t) = \mathbf{f} \left( \mathbf{x}(t), \mathbf{p}\right) \:, \quad
	\mathbf{x}(t_0) = \mathbf{x}_{0}
	\label{ivp_def}
\end{equation}
\noindent
where $\mathbf{x}(t) = \left( S_h(t), E_h(t), I_h(t), R_h(t), S_v(t), E_v(t), I_v(t), C(t) \right) \in \R^8$
is the vector of states at time $t$,
$\mathbf{x}_{0} = \left( S_h^i, E_h^i, I_h^i, R_h^i, S_v^i, E_v^i, I_v^i, C^i \right) \in \R^8$
is a prescribed initial condition vector referring to the initial time $t_0$ of
the analysis, the vector
$\mathbf{p} = \left( N_h, \beta_h, \alpha_h, \gamma, N_v, \beta_v, \alpha_v, \delta \right) \in \R^8$
lumps the model parameters and
$\mathbf{f}: U \subset \R^8 \times \R^8 \to \R^8$ is a nonlinear map which
gives the evolution law of this dynamics, defined (for fixed $t$) on the open set
\begin{equation}
	U = \left\lbrace (\mathbf{x}(t),\mathbf{p}) \in \R^8 \times \R^8 ~\big |~ x_n(t) > 0 ~ \mbox{and} ~ p_n > 0, ~\mbox{for} ~ n=1,\cdots, 8 \right\rbrace .
\end{equation}

The forward problem consists in providing initial conditions (IC) and a set of parameters,
represented by the pair $\boldsymbol \alpha = (\mathbf{x}_{0},\mathbf{p})$,
and compute by means of numerical integration
the model response $\mathbf{x}(t)$ from which a scalar
observable $\boldsymbol \phi (\boldsymbol \alpha, t)$ is obtained. In the forward problem,
$\boldsymbol \alpha$ represents all IC and system parameters
from Eq.(\ref{model_eqs}), while $\boldsymbol \phi (\boldsymbol \alpha, t)$ is the new cases
$\mathcal{N}_{w}$ system response from Eq.(\ref{NC_eq}).

Since the map $\mathbf{f}$ has a polynomial nature in $\mathbf{x}$, it is a continuously
differentiable function in $\mathbf{x}$. Thus, the existence and uniqueness theorem for ordinary
differential equations guarantees that the initial value problem of (\ref{ivp_def}) has an
unique solution. Besides that, one can also show that this solution \sanswer{is} continuously dependent on
$\boldsymbol \alpha$, as well as the \emph{forward map} $\boldsymbol \phi$ \cite{Hirsch2013,Perko2006}.

The evaluation of the system response in the \emph{forward problem} is performed
numerically in this work via a Runge-Kutta (4,5) method and the scalar observable
of interest $\mathcal{N}_{w}$ is  used to assess the simulation when compared with real
data of the 2016 outbreak made available by the Brazilian Ministry of Health
\cite{Data2016} (supplementary material A). The referred data consists of probable cases of infected people per EW,
registered by sanitary outposts and health institutions throughout the country
when the common symptoms of Zika fever were exhibited by an individual.
In accordance with the hypothesis that one only displays symptoms when inside
$I_h(t)$, the $C(t)$ variable models this discrete accumulating data on a continuous sense and
$\mathcal{N}_{w}$ provides the corresponding influx per EW. The $C(t)$ time series
is also observed as a criteria for good fitting, since it signifies the general impact
of the epidemic and because a reasonable $\mathcal{N}_{w}$ result does not necessarily implies
an acceptable cumulative number of infectious individuals for all $t$ in comparison to the data.

The initial time of the analysis was established as the first EW of 2016.
The remaining individuals in both populations are assumed susceptible at first, meaning
$S_{h}^{i} = N_h - E_{h}^{i} - I_{h}^{i} - R_{h}^{i}$ and $S_{v}^{i} = N_v - E_{v}^{i} - I_{v}^{i}$.
The initial values for the exposed and infectious groups are set equal,
$E_{h}^{i} = I_{h}^{i}$ and $E_{v}^{i} = I_{v}^{i}$.
Likewise, the number of infected humans at the initial time must be $I_{h}^{i} = C^{i}$,
given its definition. The value of $C^i$ is taken as the number of confirmed Zika cases in Brazil
on the first EW of $2016$ \cite{Data2016}, $8{,}201$ individuals, and the recovered
 group is assumed equal to the suspected number of infected in $2015$, according
to the data available \cite{Faria2016}, $R_{h}^{i} = 29{,}639$ individuals.
As for the proportion of infectious vectors in the first week, to work around the lack of
data for this initial condition, repetitive manual
estimations were tried until the resulted time series of $\mathcal{N}_w$ presented reasonable
values compared to the real data.
It became clear that the system response is very
sensible to $I_{v}^{i}$, as slight variations in
its value are required to achieve feasible results. In the process of choosing its value,
the matching of the $\mathcal{N}_w$ curve's peak to the amplitude of infection is also a priority,
since this is the main interest region for evaluation of the outbreak.
The nominal values of the parameters exhibited viable $\mathcal{N}_w$ curves around
$I_{v}^{i}= 2.2 \times 10^{-4}$.

Figure~\ref{graph_forward_problem} presents the SEIR-SEI model response for the nominal set of parameters
from section~\ref{NomParam}, supplied with the above IC, on \sanswer{an}
epidemiological week temporal domain consisting of one to fifty-two weeks ($7$ to $365$ days), compared
with the data of the outbreak (red dots).

\begin{figure}[h]
\centering
\subfigure[Cumulative number]{\includegraphics[scale=0.32]{forward_C.pdf} \label{graph_forward_C}}
\subfigure[New cases]{\includegraphics[scale=0.32]{forward_NewCases.pdf}\label{graph_forward_NewCases}}\\
\vspace{-3mm}
\caption{Outbreak data (red) and model response using the nominal parameters (blue).}
\label{graph_forward_problem}
\end{figure}

The general shape of the curves in Figure~\ref{graph_forward_problem} do provide qualitative information
regarding the evolution of the infection, even though the portrayed descriptions
are not quantitatively realistic. This inherent pattern agreement and numerical
mismatch \sanswer{suggests} that the model response may differ from the real data due to the use of
unsuitable values for the parameters or incorrect IC assumptions.
The search for parameter and initial condition values that make the simulations fit well to
the observed data defines model calibration (or system identification),
being the object of interest of the next section.

\subsection{Inverse Problem}
\label{inverse_problem}

The model calibration problem seeks to find a set of parameters
that, to a certain degree, makes the model response as close
as possible to the empirical observations (reference data), once,
due to the erros on model conception and reference data acquisition,
it is (practically) impossible for the  forward map to reproduce the
outbreak observations.

The mathematical setting for this case considers the
\emph{parameters vector} $\boldsymbol \alpha$ defined
in the \emph{parameter space} $E = \R^{12}$, since here
$\boldsymbol \alpha$ comprises all IC and system parameters
from Eq.(\ref{model_eqs}) , excepting $N_h$, $N_v$, $R_h^i$ and $C^i$,
which are kept fixed in their nominal values.
For the purpose of comparison between observations and predictions,
a discrete set with $M$ time-instants is considered, so that
scalar observations and predictions are respectively lumped into
$\vec{y} = (y_1, y_2, \cdots, y_M)$ and
$\bm{\phi} \left( \bm{\alpha} \right) = \left(\phi_1, \phi_2, \cdots , \phi_M \right)$,
both defined in the \emph{data space} $F=\R^{M}$.
Note that the \emph{forward map} $\boldsymbol \phi: E \to F$ associates to
each \sanswer{parameters} vector $\bm{\alpha}$ an observable vector
$\bm{\phi} \left( \bm{\alpha} \right)$ where the component
represents the number of new cases in each week of the year,
i.e., $\bm{\phi}_{w} = \mathcal{N}_{w}$.
In practice, the parameters vector is restricted to be on the convex
\emph{set of admissible values}
$C = \left\lbrace \boldsymbol \alpha \in E ~ \big | ~ \textbf{lb}\, \leq
\boldsymbol \alpha \leq \,\textbf{ub} \right\rbrace$, in which
$\textbf{lb}$ and $\textbf{ub}$ are lower and upper vector bounds
for $\boldsymbol \alpha$, respectively.

In formal terms, given an \emph{observation vector} $\vec{y} \in F$ and a
\emph{prediction vector} $\bm{\phi} \left(\boldsymbol \alpha \right) \in F$,
the calibration aims at finding a vector of parameters $\boldsymbol \alpha^{\ast}$ such that
\begin{equation}
\boldsymbol \alpha^{\ast} = \argmin_{\boldsymbol \alpha \in C} \, J(\boldsymbol \alpha)  \:,
\label{inv_prob_def}
\end{equation}
\noindent
for a misfit function
\begin{equation}
J(\boldsymbol \alpha) = ||\vec{y} - \bm{\phi} \left(\boldsymbol \alpha \right) ||^2
= \left\lbrace \sum_{m=1}^{M} \Big| y_m - \phi_m \left(\boldsymbol \alpha  \right) \Big|^{2} \right \rbrace.
\label{misfit_func}
\end{equation}

This is the \emph{inverse problem} associated
to the epidemic model. In general this type of problem
is extremely nonlinear, with none or low regularity,
multiple solutions (or even none),
being much more complicated to attack in comparison
with the forward problem
\cite{Aster2012,Yaman2013}.
A schematic representation of the forward and the inverse problem
associated to the epidemic model is shown in Figure~\ref{InvProb_schematic}.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.40]{direct_inv_problem.pdf}
\vspace{-3mm}
\caption{Schematic representation of the forward and inverse problems associated to the epidemic model.}
\label{InvProb_schematic}
\end{figure}

This inverse problem attempts to estimate a finite number of parameters on a
finite dimension space, being defined in terms of a typical nonlinear misfit function.
\fanswer{Therefore, Theorem 4.5.1} of \cite{Chavent2010} can be invoked to guarantee
a proper sense of well-posedness (existence, uniqueness, unimodality and local stability)
for the \fanswer{inverse} problem.

The Trust-Region-Reflective method \fanswer{(TRR)} is employed here to numerically approximate a
solution for the inverse problem (\ref{inv_prob_def}). The main idea of the method is to minimize a simpler function
that reflects the behavior of $J(\boldsymbol \alpha)$
in a neighborhood (trust-region) around $\boldsymbol \alpha$. The simpler function is defined as dependent on
the trial step $\mathbf{s}$, characterizing the Trust-Region subproblem, and its computation is optimized
by restricting the subproblem to a two-dimensional subspace.
The subspace is linear spanned by a multiple of the gradient $\mathbf{g}$ and
(in the bounded case) a vector obtained in a scaled modified Newton step, used for the
convergence condition ${D(\boldsymbol \alpha)}^{-2}\,\mathbf{g}(\boldsymbol \alpha) = 0$,
where $D$ is a diagonal matrix that depends on $\boldsymbol \alpha$, $\mathbf{g}$, $\mathbf{lb}$,
and $\mathbf{ub}$ \cite{Coleman1996}. Finally, the trial step is found through the subproblem as
\begin{equation}
\mathbf{s^*} =
\argmin_\mathbf{s} \left\lbrace
\frac{1}{2} \mathbf{s}^\mathrm{\!T} Q \mathbf{s} + \mathbf{g}^\mathrm{\!T} \mathbf{s} \: \mid \:
{||D\,\mathbf{s}||}_2 \leq \Delta \right\rbrace \:,
\label{subproblem_eq}
\end{equation}
where $\Delta$ is a scalar associated with the trust region size; $Q$ is a matrix involving
$D$, a Jacobian matrix (also dependent on $\boldsymbol \alpha$, $\mathbf{g}$, $\mathbf{lb}$,
and $\mathbf{ub}$) and an approximation of the Hessian matrix \cite{Coleman1996}.
The quadratic approximation in Eq.(\ref{subproblem_eq}) has well-behaved solutions \cite{Conn2000}
and if $J(\boldsymbol \alpha + \mathbf{s}) < J(\boldsymbol \alpha) $
then $\boldsymbol \alpha$ is updated to $ \boldsymbol \alpha + \mathbf{s} $ and the process iterates,
otherwise $\Delta$ is decreased. In addition, a reflection step also occurs if a given step
intersects a bound: the reflected step is equivalent to the original step except in the
intersecting dimension, where it assumes the opposite value after reflection.

The TRR algorithm also requires an initial guess for each parameter, identified in the
next section as ``TRR input''. The stopping criteria are the norm of the step and the
change in the value of the objective function, with a tolerance of $10^{-7}$.
Supplementary material B provides details on the \fanswer{software} used for implementation.

\subsection{Numerical experiments for calibration}
\label{NumExpCalib}

Figure~\ref{graph_inverse_problem_N} presents the best result for the $\mathcal{N}_{w}$ system response
fitting problem using the nominal parameters and IC from sections~\ref{NomParam} and \ref{forward_problem}
as initial guesses for the TRR algorithm. The upper and lower bounds used for the parameters
were set compatible with the literature suggested intervals and are presented in Table~\ref{tab_inverse_N},
along with the parameters and IC values resulted from the calibration (``TRR output''). The
\sanswer{$\mathbf{ub}$} for $\delta$ was assumed \fanswer{lower} than the \sanswer{$\mathbf{lb}$}
for $\alpha_v$ to maintain consistency with the model
interpretation. The minima for $S_h^i$ and $S_v^i$ were set to
$0.9\,N_h$ and $0.99$, respectively, to establish a high number
of susceptible individuals as is expected for the beginning of an outbreak. Also, the
\sanswer{lower and upper bound} ($0.999$) for $S_v^i$ were
motivated by noticing how variations in $E_v^i$ and $I_v^i$ of order $10^{-3}$ already
bring significant changes in the system response. The lack of available data
for the exposed and infectious groups at the onset of the epidemic was circumvented by
appointing its minimum and maximum possible values as $\mathbf{lb}$ and \fanswer{$\mathbf{ub}$}, i.e.,
$E_v^i$ and $I_v^i$ were restricted between zero and one, while
$E_h^i$ and $I_h^i$ were bounded by zero and $N_h$.

Additionally, to ensure the model hypotheses of compartmentalization
and constant population, two additional fitting points were defined, $ \Sigma_h = S_h^i + E_h^i + I_h^i$ and
$\Sigma_v = S_v^i + E_v^i + I_v^i$, which were set to match $N_h - R_h^i$ and
$N_v=1$ on Eq.(\ref{inv_prob_def}), respectively. However, the algorithm is only capable
of \textbf{approximating} $\Sigma_h$ and $\Sigma_v$ to their intended values. So to
account for these minor differences, the resulting values of  $(N_h - R_h^i) - \Sigma_h$ and
$1 - \Sigma_v$ were added to the TRR output of $S_h^i$ and $S_v^i$. These corrections did not
impact the calibration, since the scale of the differences would always be, correspondingly,
$10^{-4}$ and $10^{-2}$ (at most), which are below the sensibility of  $S_h^i$ and $S_v^i$; thus,
they were only exacted to keep the hypotheses rigorously sustained, otherwise the sum of the
compartments in each population would quickly tend to $N_h$ and $N_v$ in an asymptotic fashion.

\begin{figure}[h]
\centering
\subfigure[Cumulative number]{\includegraphics[scale=0.32]{inverse_C_N.pdf} \label{graph_inverse_C_N}}
\subfigure[New cases]{\includegraphics[scale=0.32]{inverse_NewCases_N.pdf}\label{graph_inverse_NewCases_N}}\\
\vspace{-3mm}
\caption{Outbreak data (red) and calibrated model response (blue) from Table~\ref{tab_inverse_N}.}
\label{graph_inverse_problem_N}
\end{figure}

Clearly, Figure~\ref{graph_inverse_problem_N} is a reasonable
description of the outbreak: the general shape of the infection evolution is attained,
all parameters are within realistic possibilities, the $C(t)$ curve overshoots the data by merely $6.00\%$,
the peak value of $\mathcal{N}_w$ differs from the empirical data maximum by $7.87\%$ and is
only one week off. However, taking into consideration the order of
magnitude of the first data point ($C^i = 8{,}201$) and the scale of the infection
($215{,}319$ probable cases until the $52$th EW \cite{SVS2017}), the TRR output for
$I_h^i$ ($253{,}360$) is probably too high, even though there is no reference value
to compare with the number of infectious individuals at the beginning of 2016,
making it difficult to ascertain on a deterministic manner what is a feasible value for $I_h^i$.

\begin{table}[h!]
\vspace{2mm}
\centering
\begin{tabular}{cllllll}
\toprule
$\boldsymbol \alpha$    & TRR input            & $\mathbf{lb}$ & $\mathbf{ub}$ & TRR output & Reference \\
\midrule
$\alpha_{h}$ & $1/5.9$              & $1/12$        & $1/3$      & $1/12$  &
							\cite{Ott2016,Ioos2014,ECDC2015,Valentine2016,Villela2017,Ferguson2016} \\
$\alpha_{v}$ & $1/9.1$              & $1/10$        & $1/5$      & $1/10$  &
							\cite{Wong2013,Chouin-Carneiro2016,Ferguson2016} \\
$\gamma$     & $1/7.9$              & $1/8.8$       & $1/3$      & $1/8.8$ &
							\cite{Villela2017,Funk2016,Ferguson2016} \\
$\delta$     & $1/11$               & $1/21$        & $1/11$     & $1/16.86$  &
							\cite{CDCurl,Nelson1986,Villela2017,Freitas2007} \\
$\beta_{h}$  & $1/11.3$             & $1/16.3$      & $1/8$      & $1/16.3$    &
							\cite{Ferguson2016} \\
$\beta_{v}$  & $1/8.6$              & $1/11.6$      & $1/6.2$    & $1/11.6$    &
							\cite{Ferguson2016} \\
$S_{h}^{i}$  & $205{,}953{,}959$    & $0.9\times N_h$ & $N_h$    & $205{,}700{,}000$  &
							 ----------- \\
$E_{h}^{i}$  & $8{,}201$            & $0$           & $N_h$      & $15{,}089$ &
							 ----------- \\
$I_{h}^{i}$  & $8{,}201$            & $0$           & $N_h$      & $253{,}360$  &
							 ----------- \\
$S_{v}^{i}$  & $0.99956$             & $0.99$        & $0.999$   & $1$        &
							 ----------- \\
$E_{v}^{i}$  & $2.2 \times 10^{-4}$ & $0$           & $1$        & $0$       &
							 ----------- \\
$I_{v}^{i}$  & $2.2 \times 10^{-4}$ & $0$           & $1$        & $0$     &
							 ----------- \\
\bottomrule
\end{tabular}
\caption{TRR setup for the Figure~\ref{graph_inverse_problem_N} calibrated response. The values referring to
parameters are in $\text{days}^{-1}$, human IC are expressed in individuals, and vector quantities in proportion. }
\label{tab_inverse_N}
\end{table}

Figure~\ref{graph_inverse_problem_10k} allows examination of the system behavior when the $I_h^i$
value is around $C^i$, by depiction of another result to the inverse problem when the upper
bound of the initial number of infectious individuals is set to $10{,}000$. The same restriction
was made over $E_h^i$ merely to simplify the analysis. Table~\ref{tab_inverse_10k}
displays the resulting parameters and IC.

The model response for a $10{,}000$ $I_h^i$ restriction also presents acceptable
predictions of the general shape and numbers of the outbreak, even though it is less
accurate than Figure~\ref{graph_inverse_problem_N} on a fitting criteria for $\mathcal{N}_w$.
For comparison, the $\mathcal{N}_w$ peak and data maximum difference
increased to $10.57\%$ and two weeks, while the overshoot on the
$C(t)$ time series actually reduced to $5.74\%$.

\begin{figure}[h!]
\centering
\subfigure[Cumulative number]{\includegraphics[scale=0.319]{inverse_C_10k.pdf}
\label{graph_inverse_C_10k}}
\subfigure[New cases]{\includegraphics[scale=0.32]{inverse_NewCases_10k.pdf}
\label{graph_inverse_NewCases_10k}}\\
\vspace{-3mm}
\caption{Outbreak data (red) and calibrated model response (blue) from Table~\ref{tab_inverse_10k}.}
\label{graph_inverse_problem_10k}
\end{figure}

\begin{table}[h!]
\vspace{2mm}
\centering
\begin{tabular}{cllllll}
\toprule
$\boldsymbol \alpha$    & TRR input            & $\mathbf{lb}$ & $\mathbf{ub}$ & TRR output & Reference \\
\midrule
$\alpha_{h}$ & $1/5.9$              & $1/12$        & $1/3$      & $1/12$  &
							\cite{Ott2016,Ioos2014,ECDC2015,Valentine2016,Villela2017,Ferguson2016} \\
$\alpha_{v}$ & $1/9.1$              & $1/10$        & $1/5$      & $1/10$  &
							\cite{Wong2013,Chouin-Carneiro2016,Ferguson2016} \\
$\gamma$     & $1/7.9$              & $1/8.8$       & $1/3$      & $1/3$ &
							\cite{Villela2017,Funk2016,Ferguson2016} \\
$\delta$     & $1/11$               & $1/21$        & $1/11$     & $1/21$  &
							\cite{CDCurl,Nelson1986,Villela2017,Freitas2007} \\
$\beta_{h}$  & $1/11.3$             & $1/16.3$      & $1/8$      & $1/10.40$    &
							\cite{Ferguson2016} \\
$\beta_{v}$  & $1/8.6$              & $1/11.6$      & $1/6.2$    & $1/7.77$    &
							\cite{Ferguson2016} \\
$S_{h}^{i}$  & $205{,}953{,}959$    & $0.9\times N_h$ & $N_h$ & $205{,}953{,}534$  &
							 ----------- \\
$E_{h}^{i}$  & $8{,}201$            & $0$           & $10{,}000$ & $6{,}827$ &
							 ----------- \\
$I_{h}^{i}$  & $8{,}201$            & $0$           & $10{,}000$ & $10{,}000$  &
							 ----------- \\
$S_{v}^{i}$  & $0.99956$            & $0.99$        & $0.999$    & $0.999586$     &
							 ----------- \\
$E_{v}^{i}$  & $2.2 \times 10^{-4}$ & $0$           & $1$        & $4.14 \times 10^{-4}$       &
							 ----------- \\
$I_{v}^{i}$  & $2.2 \times 10^{-4}$ & $0$           & $1$        & $0$     &
							 ----------- \\
\bottomrule
\end{tabular}
\caption{TRR setup for the Figure~\ref{graph_inverse_problem_10k} calibrated response. The values referring to
parameters are in $\text{days}^{-1}$, human IC are expressed in individuals, and vector quantities in proportion. }
\label{tab_inverse_10k}
\end{table}

To compare the two systems defined by the parameters and IC from each table,
Figure~\ref{graph_inverse_problem_IH_comparison} portrays their $I_h$ response. The system
from Table~\ref{tab_inverse_N} has an almost monotonically decreasing $I_h$ curve, except for a
slight local maximum around the $6$th EW, while the system from Table~\ref{tab_inverse_10k}
presents a significant increase in the number of infectious individuals by the same time, which
correspond to the weeks right before the peak infection. As stated, the lack of empirical data
for the current number of infectious at each EW makes it impossible to determine
what is a quantitatively reasonable prediction for $I_h$ values.
But this work assumes that a more possible scenario involves a
$I_h$ time series that also follows the general shape of an epidemic curve around the weeks of maximum
infection, especially considering that Figure~\ref{graph_inverse_IH_N}
implicates the notion that most people were infected strictly before $2016$, which does not seems
the case suggested by the Brazilian outbreak data of probable cases of infected per EW \cite{SVS2017}
when compared to numbers available for $2015$ \cite{Faria2016}. Thus, with this qualitative criteria
in mind, the system from Table~\ref{tab_inverse_10k} is selected for a further analysis over its
behavior dependency to the $I_h^i$ initial condition.

\begin{figure}[h!]
\centering
\subfigure[$I_h(t)$ response from Table~\ref{tab_inverse_N}]
{\includegraphics[scale=0.32]{inverse_IH_N.pdf} \label{graph_inverse_IH_N}}
\subfigure[$I_h(t)$ response from Table~\ref{tab_inverse_10k}]
{\includegraphics[scale=0.32]{inverse_IH_10k.pdf}\label{graph_inverse_IH_10k}}\\
\vspace{-3mm}
\caption{Comparison between $I_h(t)$ responses from Table~\ref{tab_inverse_N} (left) and Table~\ref{tab_inverse_10k} (rigth).}
\label{graph_inverse_problem_IH_comparison}
\end{figure}

Figure~\ref{graph_inverse_problem_range} displays the $C(t)$, $\mathcal{N}_w$ and $I_h(t)$ responses per EW
for various $I_h^i$ on the system with parameters from Table~\ref{tab_inverse_10k}. To simplify the analysis,
$E_h^i$ is considered equal to $I_h^i$ in each case. The remaining IC are the same from the Table.
The pattern suggests that a $I_h^i$ increase on the system with this given set of parameters
continuously \sanswer{escalates} the $C(t)$ and $\mathcal{N}_w$ curves, eventually overshooting
the data by far, and \sanswer{reduces} the variations of $I_h(t)$ curve around its local maximum.
Figure~\ref{graph_inverse_problem_range} allows one to make better predictions about the outbreak by
analyzing the multiple possible scenarios of epidemic evolution over different values for the IC missing
empirical information. Clearly, the system response in all cases is qualitative reliable in simulating
the outbreak ($C(t)$ and $\mathcal{N}_w$ shape) and can quantitatively approximate the real data values
for some values of $I_h^i$.

\begin{figure}[h!]
\centering
\subfigure[Cumulative number]{\includegraphics[scale=0.32]{multi_C.pdf}
\label{graph_inverse_multi_C}}
\subfigure[New cases]{\includegraphics[scale=0.319]{multi_NewCases.pdf}
\label{graph_inverse_multi_NewCases}}\\
\subfigure[$I_h$ response]{\includegraphics[scale=0.32]{multi_IH.pdf} \label{graph_inverse_multi_IH}}
\subfigure[$I_h$ response magnified]{\includegraphics[scale=0.32]{multi_IH_zoom.pdf}\label{graph_inverse_multi_IH_zoom}}\\
\caption{Multiple $C(t)$, $\mathcal{N}_w$ and $I_h(t)$ responses per EW using the parameters from Table~\ref{tab_inverse_10k}. The values on the legend correspond to the used $I_h^i$ and $E_h^i$ to generate the curves. The other IC are also from Table~\ref{tab_inverse_10k}. The bottom right graph is the $I_h$ comparison magnified around the local maximum region. The red circles are the previously used real data.}
\label{graph_inverse_problem_range}
\end{figure}

\section{Concluding remarks}
\label{concl_remaks}

A SEIR-SEI epidemic model to describe the dynamics of the 2016 Zika virus outbreak in Brazl is developed
and calibrated in this work. Nominal quantities for the parameters are selected from the related literature
concerning the Zika infection, the \emph{Aedes aegypti} genus of mosquitoes, vector-born epidemic models and
information provided by health organizations. The calibration process is done through the solution of an
inverse problem with the aid of a Trust-Region-Reflective method, used to pick the best parameter values
that would fit the model response for the number of new infectious cases per week into the disease's
empirical data. Results within realistic values for the parameters are presented, stating reasonable
descriptions with the curve shape similar to the outbreak evolution and
proximity between the estimated peak value and data for maximum number of infected during 2016. Further
analysis of the results about the lack of data for an initial condition is performed, exhibiting a range
of values over which the system response keeps its quantitative reliability to a certain degree.

This work is part of a long  project of  modeling and prediction of epidemics related to the Zika virus in the  Brazilian context \cite{Cunhajr_ccis2016,Cunhajr_cnmac2017_1}. In upcoming studies the authors intend
to take into account the uncertainties underlying the model parameters via Bayesian updating
and employ an Active Subspace approach \cite{Constantine2015,Constantine2018,Constantine2014}
to explore relevant scenarios in parametric studies.

\section*{Acknowledgments}

The authors are indebted to the Brazilian agencies
CNPq (National Council for Scientific and Technological Development),
CAPES (Coordination for the Improvement of Higher Education Personnel)
and FAPERJ (Research Support Foundation of the State of Rio de Janeiro)
for the financial support given to this research. They are also grateful to
Jo\~{a}o Peterson and Vin\'{\i}cius Lopes, both engineering students at UERJ,
for the collaboration in the initial stages of this work.
The anonymous reviewers made a series of comments
and suggestions that greatly enriched this paper, for which the authors
are very grateful.

\section*{References}

\section{Forward Problem}

The default values for the parameters and initial conditions (IC) are the TRR output from Table 2,
section 3.3 of the main document.

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,caption = Forward problem - main file]
Nh = 206e6;
bH = 1/10.40;
aH = 1/12;
yH = 1/3;

Nv = 1;
bV = 1/7.77;
aV = 1/10;
dV = 1/21;

param = [Nh bH aH yH Nv bV aV dV];

SH0 = 205953534;
EH0 = 6827;
IH0 = 10000;
RH0 = C2015;
SV0 = 0.999586;
EV0 = 4.14e-4;
IV0 = 0;
C0  = NewCasesData(1);

IC = [SH0 EH0 IH0 RH0 SV0 EV0 IV0 C0];

dt = 0.1;
tspan = 7:dt:365;

opt = odeset('RelTol',1.0e-6,'AbsTol',1.0e-9);
[time y] = ode45(@(t,y)rhs_SEIR_SEI(t,y,param),tspan,IC,opt);

SH = y(:,1); EH = y(:,2); IH = y(:,3); RH = y(:,4);
SV = y(:,5); EV = y(:,6); IV = y(:,7);

C  = y(:,8);

NewCases = C(1);
for i = 2:52
    NewCases(i) = C( (7*i-7)/dt +1 ) - C( (7*(i-1)-7)/dt +1 );
end

figure(1)
plot(time/7,C)
hold on
scatter([1:52],CData)

figure(2)
stem(NewCasesData)
hold on
stem(NewCases)
\end{lstlisting}

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,caption = Forward problem - right hand side]
function ydot = rhs_SEIR_SEI(t,y,param)

Nh  = param(1); bH  = param(2); aH  = param(3); yH  = param(4);
Nv  = param(5); bV  = param(6); aV  = param(7); dV  = param(8);

ydot = zeros(size(y));
ydot(1) = - bH*y(1).*y(7)/Nv;
ydot(2) = bH*y(1).*y(7)/Nv - aH*y(2);
ydot(3) = aH*y(2) - yH*y(3);
ydot(4) = yH*y(3);
ydot(5) = dV - bV.*y(5).*y(3)/Nh - dV*y(5);
ydot(6) = bV*y(5).*y(3)/Nh - aV*y(6) - dV*y(6);
ydot(7) = aV*y(6) - dV*y(7);
ydot(8) = aH*y(2);

end
\end{lstlisting}

\section{Inverse Problem}

See section 2.3 and 2.4 of the main document for the nominal parameters and initial conditions used as
initial guesses here. The default results are those from Table 2, section 3.3 of the main document.

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,
caption = Inverse problem main file]
Nh = 206e6;
bH = 1/11.3;
aH = 1/5.9;
yH = 1/7.9;

Nv = 1;
bV = 1/8.6;
aV = 1/9.1;
dV = 1/11;

IH0 = NewCasesData(1);
EH0 = IH0;
RH0 = C2015;
SH0 = Nh - (IH0 + EH0 + RH0);

IV0 = 0.00022;
EV0 = IV0;
SV0 = Nv - (EV0 + IV0);

C0  = NewCasesData(1);

inital_guess = [bH aH yH bV aV dV SH0 EH0 IH0 SV0 EV0 IV0];
fixed_TRRparam = [Nh Nv C0];

data = [NewCasesData Nh-RH0 Nv];

data_pos = [7:7:365 400 450];

lb = [1/16.3, 1/12, 1/8.8, 1/11.6, 1/10, 1/21, 0.9*Nh, ...
	0,     0,  0.99,  0,  0];
up = [   1/8,  1/3,  1/3,   1/6.2,  1/5,  1/11,    Nh, ...
    10000, 10000, 0.999, Nv, Nv];

dt = 1;
tspan = 7:dt:365;

F = @(variable_TRRparam,data_pos)TRR_FunctionOutput(...
	variable_TRRparam,data_pos,fixed_TRRparam,tspan,dt);

options = optimoptions('lsqcurvefit','Display','iter', ...
	  'Algorithm','trust-region-reflective', ...
	  'MaxIterations',100000,'MaxFunctionEvaluations',  ...
	  500000,'FunctionTolerance',1e-007,'StepTolerance', ...
	  1e-007,'FiniteDifferenceType','central');

result_TRRparam = lsqcurvefit(F,inital_guess,data_pos,data,lb, ...
		  up,options);

new_param = [Nh result_TRRparam(1:3) Nv result_TRRparam(4:6)];
new_IC = [result_TRRparam(7:9) RH0 result_TRRparam(10) ...
	  result_TRRparam(11:12) C0];

Xh = Nh - RH0 - sum(result_TRRparam(7:9));
Xv = Nv - sum(result_TRRparam(10:12));
new_IC(1) = result_TRRparam(7) + Xh;
new_IC(5) = result_TRRparam(10) + Xv;

[time y] = ode45(@(t,y)rhs_SEIR_SEI(t,y,new_param),tspan,new_IC);
opt = odeset('RelTol',1.0e-6,'AbsTol',1.0e-9);

SH = y(:,1); EH = y(:,2); IH = y(:,3); RH = y(:,4);
SV = y(:,5); EV = y(:,6); IV = y(:,7);

C  = y(:,8);

NewCases = C(1);
for i = 2:52
    NewCases(i) = C( (7*i-7)/dt +1 ) - C( (7*(i-1)-7)/dt +1 );
end

figure(1)
plot(time/7,C)
hold on
scatter([1:52],CData)

figure(2)
stem(NewCasesData)
hold on
stem(NewCases)
\end{lstlisting}

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,
caption = Inverse problem output function]
function F = TRR_FunctionOutput(variable_TRRparam,data_pos,fixed_TRRparam,tspan,dt)

IC = [variable_TRRparam(7:12) fixed_TRRparam(3)];

opt = odeset('RelTol',1.0e-6,'AbsTol',1.0e-9);
[time y] = ode45(@(t,y)TRR_rhs_SEIR_SEI(t,y,variable_TRRparam, ...
		fixed_TRRparam),tspan,IC,opt);

SumH = y(1,1) + y(1,2) + y(1,3);

SumV = y(1,4) + y(1,5) + y(1,6);

C  = y(:,7);

NewCases = C(1);
for i = 2:52
    NewCases(i) = C( (7*i-7)/dt +1 ) - C( (7*(i-1)-7)/dt +1 );
end

F = [NewCases SumH SumV];

end
\end{lstlisting}

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,
caption = Inverse problem right hand side]
function ydot = TRR_rhs_SEIR_SEI(t,y,variable_TRRparam, ...
				fixed_TRRparam)

Nh   = fixed_TRRparam(1);
bH  = variable_TRRparam(1);
aH  = variable_TRRparam(2);
yH  = variable_TRRparam(3);

Nv = fixed_TRRparam(2);
bV  = variable_TRRparam(4);
aV  = variable_TRRparam(5);
dV  = variable_TRRparam(6);

ydot = zeros(size(y));

ydot(1) = -bH*y(1).*y(6)/Nv;
ydot(2) = bH*y(1).*y(6)/Nv - aH*y(2);
ydot(3) = aH*y(2) - yH*y(3);
ydot(4) = dV - bV.*y(4).*(y(3)/Nh) - dV*y(4);
ydot(5) = bV*y(4).*(y(3)/Nh) - aV*y(5) - dV*y(5);
ydot(6) = aV*y(5) - dV*y(6);
ydot(7) = aH*y(2);

end
\end{lstlisting}

\section{Changing human population}

A recruitment term ($\eta\,N_h$) and a mortality rate ($\mu$) are needed to implement a
changing human population. The human variables in Eqs.~(1) in section 2.2 of the main document become

\begin{equation}
\begin{aligned}
\dod{S_h}{t} &= \eta\,N_h -\beta_h \, {S_h} \, \frac{I_v}{N_v} - \mu\,S_h \:,\\
\dod{E_h}{t} &= \beta_h \, S_h \, \frac{I_v}{N_v} - \alpha_h \, E_h - \mu\,E_h \:,\\
\dod{I_h}{t} &= \alpha_h \, E_h - \gamma \, I_h - \mu\,I_h \:,\\
\dod{R_h}{t} &= \gamma \, I_h  - \mu\,R_h \:.
\end{aligned}
\end{equation}

\noindent
Also, a new differential equation must be added to the model to track the varying population

\begin{equation}
\dod{N_h}{t} = \eta\,N_h - \mu\,(S_h+E_h+I_h+R_h) \:,
\end{equation}

\noindent
and, consequently, its respective initial condition, $N_h(0)$, is required for the forward problem.
Code~\ref{vital} provides the right hand side of this variation of the model.

\begin{lstlisting}[style=Matlab-editor,basicstyle=\footnotesize,
caption = Vital Dynamics - rigth hand side,label = vital]
function ydot = rhs_SEIR_SEI(t,y,phys_param)

nH  = phys_param(1);
bH  = phys_param(2);
aH  = phys_param(3);
yH  = phys_param(4);
bV  = phys_param(5);
aV  = phys_param(6);
dV  = phys_param(7);
mH  = phys_param(8);

ydot = zeros(size(y));

ydot(1) = nH*y(9) - bH*y(1).*y(7) - mH*y(1);
ydot(2) = bH*y(1).*y(7) - aH*y(2) - mH*y(2);
ydot(3) = aH*y(2) - yH*y(3) - mH*y(3);
ydot(4) = yH*y(3) - mH*y(4);
ydot(5) = dV - bV.*y(5).*(y(3)/y(9)) - dV*y(5);
ydot(6) = bV*y(5).*(y(3)/y(9)) - aV*y(6) - dV*y(6);
ydot(7) = aV*y(6) - dV*y(7);
ydot(8) = aH*y(2);
ydot(9) = nH*y(9) - mH*(y(1)+y(2)+y(3)+y(4));

end
\end{lstlisting}

An estimation for the values of $\eta$ and $\mu$ (per thousand individuals in a year) can be found in
\cite{IBGE} for 2015 (the most recent officially published mortality and birth rates for Brazil).
The solution of the forward problem with this code can be used to observe how the human vital dynamics do
not affect $\mathcal{N}_{w}$ (neither $C(t)$) in a one year outbreak for this model. Of course, $S_h(t)$
will behave different, as expected.

\section*{References}