projects:household.tex

The generated PDF

Sections break-down

The bibtex file household.bib

The main latex file

<project-file type=“source”/> <content> \documentclass{article}

\usepackage{amsmath} \usepackage{amsthm} \usepackage{booktabs} \usepackage{graphicx} \usepackage[displaymath,mathlines]{lineno}

\newcommand{\R}{\mathcal R} \newcommand{\diag}{\text{diag}} \newcommand{\avg}[1]{\langle#1\rangle}

% the long and ugly lines below to \begin{document} are copied from % http://phaseportrait.blogspot.com/2007/08/lineno-and-amsmath-compatibility.html % for a fix around the incompatibility of lineno and amsmath % To deal with amsmath, must redefine math environments later \AtBeginDocument{% % \let\oldequation\equation% \let\endoldequation\endequation% \renewenvironment{equation}%

{\linenomath\oldequation}{\endoldequation\endlinenomath}%

% \expandafter\let\expandafter\oldequationstar\csname equation*\endcsname% \expandafter\let\expandafter\endoldequationstar\csname endequation*\endcsname% \renewenvironment{equation*}%

{\linenomath\oldequationstar}{\endoldequationstar\endlinenomath}%

% \let\oldalign\align% \let\endoldalign\endalign% \renewenvironment{align}%

{\linenomath\oldalign}{\endoldalign\endlinenomath}%

% \expandafter\let\expandafter\oldalignstar\csname align*\endcsname% \expandafter\let\expandafter\endoldalignstar\csname endalign*\endcsname% \renewenvironment{align*}%

{\linenomath\oldalignstar}{\endoldalignstar\endlinenomath}%

% \let\oldflalign\flalign% \let\endoldflalign\endflalign% \renewenvironment{flalign}%

{\linenomath\oldflalign}{\endoldflalign\endlinenomath}%

% \expandafter\let\expandafter\oldflalignstar\csname flalign*\endcsname% \expandafter\let\expandafter\endoldflalignstar\csname endflalign*\endcsname% \renewenvironment{flalign*}%

{\linenomath\oldflalignstar}{\endoldflalignstar\endlinenomath}%

% \let\oldalignat\alignat% \let\endoldalignat\endalignat% \renewenvironment{alignat}%

{\linenomath\oldalignat}{\endoldalignat\endlinenomath}%

% \expandafter\let\expandafter\oldalignatstar\csname alignat*\endcsname% \expandafter\let\expandafter\endoldalignatstar\csname endalignat*\endcsname% \renewenvironment{alignat*}%

{\linenomath\oldalignatstar}{\endoldalignatstar\endlinenomath}%

% \let\oldgather\gather% \let\endoldgather\endgather% \renewenvironment{gather}%

{\linenomath\oldgather}{\endoldgather\endlinenomath}%

% \expandafter\let\expandafter\oldgatherstar\csname gather*\endcsname% \expandafter\let\expandafter\endoldgatherstar\csname endgather*\endcsname% \renewenvironment{gather*}%

{\linenomath\oldgatherstar}{\endoldgatherstar\endlinenomath}%

% \let\oldmultline\multline% \let\endoldmultline\endmultline% \renewenvironment{multline}%

{\linenomath\oldmultline}{\endoldmultline\endlinenomath}%

% \expandafter\let\expandafter\oldmultlinestar\csname multline*\endcsname% \expandafter\let\expandafter\endoldmultlinestar\csname endmultline*\endcsname% \renewenvironment{multline*}%

{\linenomath\oldmultlinestar}{\endoldmultlinestar\endlinenomath}%

% }

\begin{document}

\begin{linenumbers}

\title{Effective Degree Household \\Network Disease Model} \author{Junling Ma$^{1}$, P. van den Driessche$^{2}$, Frederick H. Willeboordse$^{3}$
Department of Mathematics and Statistics, University of Victoria,
Victoria BC V8W 3R4, Canada
\textbf{1} jma@math.uvic.ca
\textbf{2} pvdd@math.uvic.ca
\textbf{3} willeboordse@yahoo.com } \date{\today} \maketitle

\begin{abstract} An ordinary differential equation (ODE) epidemiological model for the spread of a disease that confers immunity, such as influenza, is introduced incorporating both network topology and households. Since most individuals of a susceptible population are members of a household, including the household structure as an aspect of the contact network in the population is of significant interest. Epidemic curves derived from the model are compared with those from stochastic simulations, and shown to be in excellent agreement. An expression for the basic reproduction number of the ODE model is derived analytically and interpreted in terms of the household structure. \end{abstract}

% % % \section{Introduction} The development of realistic models is a key objective of mathematical epidemiology. The natural question that immediately arises then is which parts of a population's contact network are essential and which ones can effectively be averaged. A complicating factor is that the details of the contact network of the human population are relatively poorly understood. However, from real-world infection data as provided by for example the CDC, the actual epidemic curve can serve as an important benchmark for judging the validity of a model. Recent work has shown that homogeneous mixing in general is less reliable than models that explicitly incorporate a contact network structure for the population. In the homogeneous mixing model, see e.g. \cite{Brauer_2008}, any member of the population can infect or be infected by any other member with equal probability. Clearly this is unrealistic, but the model has the advantage that it can be analyzed relatively simply and therefore serve as a reference for other work. Models that have been successful in taking into account the population's network have been developed by, for example, Newman who used percolation theory to obtain a threshold condition for and final size of an epidemic \cite{Newman_PRE2002}, Pastor-Satorras and Vespignani who tracked the number of susceptible and infectious individuals with a certain number of contacts (i.e.~the individual's contact degree) \cite{PSV_2001a,PSV_2001b}, and Ball and Neal who used a branching process approximation \cite{Ball_1983,BallNeal_2008,Neal_2007}. Recently, Miller \cite{Miller_note} showed that the model developed by Volz \cite{Volz_JMB08} with the help of a generating function for the network degree distribution can be simplified to only three differential equations. Together with Lindquist, we showed that keeping track of the number of infectious and susceptible neighbors is a powerful approach in the \textit{effective degree} model \cite{LMDW_JMB2010}. All the network models mentioned above share in common that they are based on a (global) network topology that does not specifically deal with local structures. For example, an Erd\H{o}s-R\'enyi type random network \cite{Erdos-Renyi_1959} is often used to represent the randomness of the contacts between members of the population. While this is likely realistic to some degree for the contacts between friends and acquaintances, it appears that this may not reflect local structures. Local structures can for example be families, schools and work places where not only are contacts closer but the network is also homogeneous. For simplicity, homogeneous local structures with possibly a disease transmission rate different from that on the rest of the network are all termed households here. Indeed, it was recently shown by Nichols et al.~\cite{Nichols_2011} that household size is critical to varicella-zoster virus transmission in the tropics (Guinea Bissau), in that a larger household size compensates for a lower viral infectivity when compared to temperate climates. Recently the local structure has been represented by incorporating clusters on a contact network. For example, Newman \cite{Newman_PRL2009} developed techniques using percolation theory to analyze a random network that includes clusters. Miller \cite{Miller_PRE2009} used these techniques to find the disease threshold and epidemic size of an SIR model. This model assumes triangular clusters. However, in realistic contact networks, clusters are commonly households and work places, which are usually not limited to triangles, and the inter-cluster transmission rate is likely different from the intra-cluster transmission rate. Even more recently, Volz et al \cite{VMGM_PCB2011} introduced an SIR model that describes the disease dynamics on a contact network that has clusters of arbitrary size. They show numerically that their model agrees well with stochastic simulations, and clustering always slows the epidemic. We extend our effective degree SIR model \cite{LMDW_JMB2010} to include clusters of arbitrary size to represent households and workplaces. Even though we here restrict to SIR models, our effective degree SIS model \cite{LMDW_JMB2010} can be so extended to include household structure. Our goal is to analyze the extended SIR model and analytically derive disease threshold conditions. In order to investigate the effects of households on the spread of a disease, in Section \ref{sec:model}, a model is developed that combines households with a random inter-household network. The results are compared with stochastic numerical simulations, and excellent agreement is found. Section \ref{sec:R0} contains the derivation of analytic expressions for disease thresholds. The dependency of the threshold on parameters is discussed and compared with numerical solutions of our model. The results and implications are discussed in Section \ref{sec:discussion}. % % % \section{Model formulation} \label{sec:model}

We consider a susceptible-infectious-recovered (SIR) type of disease, where a susceptible individual is infected upon contact with an infectious individual at a rate $\beta_h$ if the infection is due to a contact with a household member and a rate $\beta_e$ if the infection is due to an external contact. An infected individual recovers with a rate $\gamma$. It is assumed that once recovered, the individual gains immunity and will not be infected again.This is a reasonable assumption for many diseases such as influenza and childhood diseases. For the purpose of modeling, immune individuals do not contribute to the dynamics and can hence be removed from the system.

Each individual of the population is part of a household, and the household sizes are uniformly distributed over a certain range. For simplicity, a lone individual is considered to be a member of a household with size 1. Each individual of a household is linked to every other individual of the same household. In other words, internally, households are homogeneously mixed. Individuals can also have external links to randomly selected individuals in other households. That is to say, externally, the links between individuals form an Erd\H{o}s–R\'enyi type random network.

The key idea of the model is to keep track of the internal links through household classes $H_{nj}$ where $n$ is the size of the household and $j$ the number of infectious individuals in the household, and to keep track of the external links with the help of the effective degree model developed earlier \cite{LMDW_JMB2010}. Thus, from the individual's point of view, each member of the population is grouped into both a household class as well as an effective degree class and the question is how transitions between classes in one group are distributed to classes in the other group. We now briefly elaborate on the household classes, the effective degree model and how they can be coupled before introducing our new model in Subsection \ref{subsec:model}.

\subsection{Intra-household transmission}

Households are classified by the household size $n$ and the number of infectious individuals $j$ in the household, denoted as $H_{nj}$. Note that we also use $H_{nj}$ to denote the number of households in the class. The household size $n$ decreases with time as infectious individuals in the household are removed upon recovery, since recovered individuals have no influence on the disease dynamics in the SIR scenario. Let $m$ be the maximum household size, then $1\leq n\leq m$, and $0\leq j\leq n$.

We assume that the individuals in a household are homogeneously mixed, and that the intra-household transmission rate is $\beta_{h}$. For a household in class $H_{nj}$, there are $n-j$ susceptible individuals, and these susceptible individuals are infected with a rate $\beta_{h}j$ by infectious household members. When a member of a household in class $H_{nj}$ is infected, the number of households in class $H_{nj}$ is reduced by one while the number of households in class $H_{n,j+1}$ is increased by one. Thus classes $H_{nj}$ move to $H_{n,j+1}$ with rate $\beta_h (n-j)jH_{nj}$ due to intra-household transmission. Similarly, classes $H_{nj}$ enter $H_{n-1,j-1}$ with a rate $\gamma jH_{nj}$ due to recovery of infected household members.

Taking into account both entry into and exit out of a class due to infection and recovery, we obtain the following expression for the household class dynamics % % H_nj-only % \begin{eqnarray} H_{nj}' & = & \beta_h \left[ (n-j+1)(j-1)H_{n,j-1}-(n-j)jH_{nj} \right] \nonumber \\ & & \!\!\!\! + \gamma\left[ (j+1)H_{n+1,j+1}-jH_{nj} \right] . \label{eq:Hnj-only} \end{eqnarray}

\subsection{Inter-household transmission}

Contacts with non-household members are tracked with the help of the effective degree model \cite{LMDW_JMB2010}, which we briefly describe here. In the effective degree model, depending on their disease status, individuals are part of a susceptible class $S_{si}$ or infectious class $I_{si}$ where the subscript $s$ indicates the number of susceptible neighbors and the subscript $i$ the number of infectious neighbors. Again, $S_{si}$ and $I_{si}$ also denote the number of individuals in each class. While in the original effective degree model all neighbors are considered, here, only \textit{external} neighbors are taken into account since the effective degree model is exclusively used for infections due to external contacts (the internal infections are tracked by the household classes $H_{nj}$). Let $M$ be the maximum external degree, then the external degree of an individual $k=s+i \leq M$, and $0\leq s,i\leq k$.

We assume that the transmission rate due to external infectious contacts is $\beta_{e}$. Consequently, a susceptible individual in class $S_{si}$ is infected by an external infectious neighbor with rate $\beta_{e}iS_{si}$. Due to the recovery of their external infectious neighbors, individuals in class $S_{si}$ move to class $S_{s,i-1}$ with rate $\gamma iS_{si}$ while due to the external infection of external susceptible neighbors, they move to class $S_{s-1,i+1}$ with a rate \[ \sum_{k=0}^{M}\sum_{u+v=k}\beta_e u v S_{uv} \frac{s S_{si}}{\sum_{k=0}^{M}\sum_{u+v=k} u S_{uv}} \] where the double sum in the numerator is the total number of external infections in the population, $s S_{si}$ is the total number of susceptible neighbors of class $S_{si}$, and the denominator is the total number of susceptible neighbors of susceptible individuals in the population. Infectious individuals in class $I_{si}$ recover with rate $\gamma I_{si}$ while they enter class $I_{s,i-1}$ with rate $\gamma i I_{si}$ due to the recovery of their external infectious neighbors. Due to infection of their external susceptible neighbors, infectious individuals in class $I_{si}$ enter class $I_{s-1,i+1}$ with rate \[ \sum_{k=0}^{M}\sum_{u+v=k}\beta_e v^{2}S_{uv} \frac{sI_{si}}{\sum_{k=0}^{M}\sum_{u+v=k}uI_{uv}} \] where the double sum in the numerator is the total number of infections amongst \textit{neighbors} of a susceptible individual, $sI_{si}$ is the total number of susceptible neighbors of class $I_{si}$, and the denominator is the total number of susceptible neighbors of infectious individuals in the population.

Thus the following set of equations for the effective degree dynamics is obtained

\begin{eqnarray} S_{si}' & = & -\beta_e iS_{si}+\gamma\left[(i+1)S_{s,i+1}-iS_{si}\right] \nonumber \\ & & +\frac{\sum_{k=0}^{M}\sum_{u+v=k}\beta_e uvS_{uv}} {\sum_{k=0}^{M}\sum_{u+v=k}uS_{uv}}\left[(s+1)S_{s+1,i-1}-sS_{si}\right] \label{eq:Ssi-only}\\ I_{si}' & = & \beta_e iS_{si}-\gamma I_{si}+\gamma\left[(i+1)I_{s,i+1}-iI_{si}\right] \nonumber \\ & & +\frac{\sum_{k=0}^{M}\sum_{u+v=k}\beta_e v^2S_{uv}} {\sum_{k=0}^{M}\sum_{u+v=k}uI_{uv}}\left[(s+1)I_{s+1,i-1}-sI_{si}\right]. \label{eq:Isi-only} \end{eqnarray} This is the effective degree SIR model presented in \cite{LMDW_JMB2010} which, with the help of stochastic simulations, was shown to be an accurate models for investigating the spread of a disease on a random network.

\subsection{Coupling terms \label{subsec:coupling}}

All the individuals in the population are members of a household (and hence part of a household class $H_{nj}$) and also have zero or more susceptible or infectious neighbors (making them part of an effective degree class $S_{si}$ or $I_{si}$). Consequently, changes due to the dynamics of $H_{nj}$ need to be reflected in $S_{si}$ and $I_{si}$, and vice versa through coupling terms.

\subsubsection{Coupling of the effective degree classes to the household classes}

As mentioned above, intra-household infections lead to decreases in the number of households in classes $H_{nj}$ and increases in the number of households in classes $H_{n,j+1}$. Similarly, inter-household infections also lead to decreases in the number of households in classes $H_{nj}$ and increases in the number of households in classes $H_{n,j+1}$, although in general at a different rate.

From the effective degree model external infections occur at a rate \[ \sum_{k=0}^{M}\sum_{s+i=k}\beta iS_{si}. \] The probability $p_{nj}$ that a newly infected individual is in household class $H_{nj}$ depends on the structure of the network. In the cases considered here, the links between \textit{individuals} form a random network with no degree correlation, no clustering, no self loops or multiple edges. Then, this probability is proportional to the number of susceptible individuals in the household, giving \begin{equation} \label{eq:pnj:node} p_{nj} = \frac{(n-j)H_{nj}}{\sum_{n=1}^m\sum_{j=0}^n (n-j)H_{nj}} \,. \end{equation} The term coupling the external infections to the household classes is then \begin{equation} \sum_{k=0}^{M}\sum_{s+i=k}\beta_e iS_{si}\frac{(n-j+1)H_{n,j-1}-(n-j)H_{nj}}{\sum_{n=1}^m\sum_{j=0}^n (n-j)H_{nj}} \,. \label{eq:Fnj} \end{equation}

\subsubsection{Coupling of the household classes to the effective degree classes}

A susceptible individual in the $S_{si}$ class can be infected by a household member and enter the $I_{si}$ class. From the household Equation (\ref{eq:Hnj-only}), such infections occur with rate \begin{equation} \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj}, \nonumber \end{equation} and these new infections now need to be distributed over the effective degree classes $S_{si}$ and $I_{si}$.

It is reasonable to assume that an infection has the same probability of occurring on any susceptible individual. Thus, the probability that it occurs in class $S_{si}$ is given by \begin{equation} \frac{S_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}}, \nonumber \end{equation} where the denominator is the total number of susceptible individuals in the population. Hence, the number of individuals that leave the $S_{si}$ class due to internal infection is given by \begin{equation} A_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{S_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \, . \end{equation}

An individual in class $S_{si}$ can also move to class $S_{s-1,i+1}$ due to an external susceptible neighbor being internally infected. On average, each newly infected individual has \begin{equation} \frac{\sum_{k=0}^{M}\sum_{s+i=k}sS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \nonumber \end{equation} susceptible neighbors. The probability that such a susceptible neighbor is in class $S_{si}$ class is given by \begin{equation} \frac{sS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}sS_{si}} \nonumber \end{equation} so that the transition rate from $S_{si}$ to $S_{s-1,i+1}$ due to internal infections of an external neighbor becomes \begin{equation} B_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{sS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}}, \nonumber \end{equation} with an analogous term $C_{si}$ for the transition from $S_{s+1,i-1}$ to $S_{si}$ given in Table \ref{table:edCouplingTerms}.

Similarly, each newly infected individual has \begin{equation} \frac{\sum_{k=0}^{M}\sum_{s+i=k}iS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \nonumber \end{equation} infectious neighbors on average and the probability that such a neighbor is in the $I_{si}$ class is given by \begin{equation} \frac{sI_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}sI_{si}}. \nonumber \end{equation} Therefore, individuals in class $I_{si}$ move to class $I_{s-1,i+1}$ due to internal infections of an external neighbor with rate \begin{equation} D_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj}{\displaystyle \frac{\sum_{k=0}^{M}\sum_{s+i=k}iS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}}}{\displaystyle \frac{sI_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}sI_{si}}}, \end{equation} with an analogous term $E_{si}$ for the transition from $I_{s+1,i-1}$ to $I_{si}$ as given in Table \ref{table:edCouplingTerms}.

Note that \begin{equation} \sum_{k=0}^{M}\sum_{s+i=k}iS_{si}=\sum_{k=0}^{M}\sum_{s+i=k}sI_{si} \,. \label{eq:Balance:SI} \end{equation} This is because the left hand side is the number of infectious neighbors of all the susceptible nodes, the right hand side is the susceptible neighbors of infectious nodes, and both correspond to the number of edges between infectious and susceptible nodes. In fact, it can be proved analogously to a computation done for the effective degree model \cite{LMDW_JMB2010}, by verifying that \begin{equation} \frac{d}{dt}\left[\sum_{k=0}^{M}\sum_{s+i=k}iS_{si}-\sum_{k=0}^{M}\sum_{s+i=k}sI_{si}\right] =-\gamma\left[\sum_{k=0}^{M}\sum_{s+i=k}iS_{si}-\sum_{k=0}^{M}\sum_{s+i=k}sI_{si}\right]\,. \label{cond:1} \end{equation} Thus, starting with the same initial condition, the two quantities must be equal for all $t\geq0$. With the help of Equation \eqref{cond:1}, the coupling term $D_{si}$ can be simplified to \vspace{-5mm} \begin{equation} D_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj}{\displaystyle \frac{sI_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}}}\label{eq:Isi} \,, \end{equation} and the term $E_{si}$ similarly simplified. The terms that couple household classes to the effective degree classes are summarized in Table \ref{table:edCouplingTerms}.

\begin{table}[!ht] \begin{tabular}{ll} \toprule

   \parbox{4cm}{Classes that leave $S_{si}$ class and enter $I_{si}$ due to internal infections.} & 
      $\displaystyle{ A_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{S_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} }$  \\ \midrule  
   \parbox{4cm}{Classes that leave $S_{si}$ due to the internal infection of a neighbor.} & 
      $\displaystyle B_{si} = sA_{si}$ \\ \midrule
   \parbox{4cm}{Classes $S_{s+1,i-1}$ that move to $S_{si}$ due to internal infection of a neighbor.} & 
      $\displaystyle{ C_{si} =  B_{s+1,i-1} }$ \\   \midrule   
   \parbox{4cm}{Classes that leave  $I_{si}$ due to the internal infection of a neighbor.} & 
      $ \displaystyle{ D_{si} = \sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj}
                     \frac{ sI_{s,i} }{ \sum_{k=0}^{M} \sum_{s+i=k} S_{si} } } $ \\ \midrule
   \parbox{4cm}{Classes $I_{s+1,i-1}$ that move to $I_{si}$ due to internal infection of a neighbor.} & 
          $ E_{si} = D_{s+1,i-1} $ \\
\bottomrule

\end{tabular} \caption{Coupling terms between the effective degree and the household classes. } \label{table:edCouplingTerms} \end{table}

\subsection{Effective Degree Household Model \label{subsec:model}}

Adding the coupling terms, we therefore obtain: \begin{align*} H_{nj}' & = H_{nj}'[\mathrm{uncoupled}] + \sum_{k=0}^{M}\sum_{s+i=k}\beta_e iS_{si}\frac{(n-j+1)H_{n,j-1}-(n-j)H_{nj}}{\sum_{n=1}^m\sum_{j=0}^n (n-j)H_{nj}} \\ S_{si}' & = S_{si}'[\mathrm{uncoupled}] - A_{si} - B_{si} + C_{si} \\ I_{si}' & = I_{si}'[\mathrm{uncoupled}] + A_{si} - D_{si} + E_{si}, \end{align*} giving the complete model as: \begin{eqnarray} H_{nj}' & = & \beta_h\left[(n-j+1)(j-1)H_{n,j-1}-(n-j)jH_{nj}\right] \nonumber \\ & & +\gamma\left[(j+1)H_{n+1,j+1}-jH_{nj}\right] \nonumber \\ & & + \sum_{k=0}^{M}\sum_{s+i=k}\beta_e iS_{si}\frac{(n-j+1)H_{n,j-1}-(n-j)H_{nj}}{\sum_{n=1}^m\sum_{j=0}^n (n-j)H_{nj}} \label{hhModel:H} \\ S_{si}' & = & -\beta_e iS_{si}+\gamma\left[(i+1)S_{s,i+1}-iS_{si}\right] \nonumber \\ & & +\frac{\sum_{k=0}^{M}\sum_{u+v=k}\beta_e uvS_{uv}}{\sum_{k=0}^{M}\sum_{u+v=k}uS_{uv}} \left[(s+1)S_{s+1,i-1}-sS_{si}\right]\nonumber \\ & & -\sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{S_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}}\nonumber \\ & & +\sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{(s+1)S_{s+1,i-1}-sS_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \label{hhModel:S}\\ I_{si}' & = & \beta_e iS_{si}-\gamma I_{si}+\gamma\left[(i+1)I_{s,i+1}-iI_{si}\right] \nonumber \\ & & +\frac{\sum_{k=0}^{M}\sum_{u+v=k}\beta_e v^2S_{uv}}{\sum_{k=0}^{M}\sum_{u+v=k}uI_{uv}} \left[(s+1)I_{s+1,i-1}-sI_{si}\right] \nonumber \\ & & +\sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{S_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \nonumber \\ & & +\sum_{n=1}^{m}\sum_{j=0}^{n}\beta_{h}(n-j)jH_{nj} \frac{(s+1)I_{s+1,i-1}-sI_{si}}{\sum_{k=0}^{M}\sum_{s+i=k}S_{si}} \label{hhModel:I} \end{eqnarray} with nonnegative initial conditions. We term this model the \textit{effective degree household model}.

We now check that our model is consistent. Since all individuals are in both a household as well as an effective degree class, the total number of susceptible (infectious) individuals, infections and recoveries must be identical in both types of classes. Indeed, a quick computation verifies that \[ \frac{d}{dt}\left[\sum_{n=1}^{m}\sum_{j=0}^n(n-j)H_{nj}-\sum_{k=0}^{M}\sum_{s+i=k}S_{si}\right]=0, \] where the first term in the brackets is the total number of susceptible individuals in the household classes, and the second term is the total number of susceptible individuals in the effective degree classes. Similarly, \[ \frac{d}{dt}\left[\sum_{n=1}^{m}\sum_{j=0}^njH_{nj}-\sum_{k=0}^{M}\sum_{s+i=k}I_{si}\right] =-\gamma\left[\sum_{n=1}^{m}\sum_{j=0}^njH_{nj}-\sum_{k=0}^{M}\sum_{s+i=k}I_{si}\right], \] where the left hand side is the change in the sum of the total number of infectious individuals in the population in both the household and effective degree classes while the right hand term is the total number of recoveries in both types of classes. The above two equations guarantee that \[ \sum_{n=1}^{m}\sum_{j=0}^n(n-j)H_{nj} = \sum_{k=0}^{M}\sum_{s+i=k}S_{si} \] and \[ \sum_{n=1}^{m}\sum_{j=0}^njH_{nj}=\sum_{k=0}^{M}\sum_{s+i=k}I_{si} \] must hold for all time $t\geq0$, given that they are true initially at $t=0$.

% % % \subsection{Comparison with stochastic simulations} \label{sec:simulations}

In order to check the accuracy of the effective degree household model, Equations (\ref{hhModel:H})–(\ref{hhModel:I}) were numerically evaluated and the results compared with stochastic simulations.

\subsubsection{The stochastic model} For the stochastic simulations, an array of nodes was constructed that keeps track of whether a node is susceptible, infectious or recovered, and its links to the other nodes. The network between the nodes is created in two steps. First, either a random household size $n_H$ with a uniform random distribution between 1 and the maximum household size $m$, or a fixed household size $n_H$ is chosen, and then $n_H$ consecutive nodes are assigned to be in the same household and fully linked. Next, a inter-household network is created using the configuration model \cite{Volz_JMB08}. This is done by repeatedly choosing two randomly selected nodes of different households and bi-directionally linking them until the desired total average degree $\langle k_T \rangle$ or the average external degree $\langle k \rangle$ is reached to form the external links ($k_T=k+n_H-1$). Thus the external links form an Erd\H{o}s–R\'enyi type network between individuals. The simulation of the disease dynamics on the constructed network is the same as used for our effective degree model \cite{LMDW_JMB2010}.

For all the simulations, the number of nodes is $N=100,000$; to reduce demographic stochasticity, we start from a large initial condition (100 initially infectious individuals); the average infectious period is taken as $5$ days ($\gamma=0.2$). For comparison, we numerically solve the ODE model using the same parameter values and initial condition as the stochastic stimulations.

\subsubsection{The agreement between ODE and stochastic models}

Figures \ref{fig:individuals} and \ref{figLfixedhh} show the comparison of the epidemic curve (the number of infectious individuals as a function of time) of the ODE model, and the average of 20 stochastically simulated epidemic curves. In Figure \ref{fig:individuals}, the network structure is held constant, while the inter- and intra-household transmission rates are varied. In Figure \ref{fig:fixedhh}, the transmission rates and the external degree distribution are held constant, while the household size distribution changes. As can be seen, the agreement between the ODE model and the stochastic simulations is excellent for all parameter values and household structures in the simulations.

% % % \section{The basic reproduction number} \label{sec:R0}

\subsection{Derivation of formula $\R_0$}

One of the most important disease parameter is the basic reproduction number $\mathcal R_0$. We apply the second generation matrix method \cite{DHM_JMB1990, DW_MB2002} to compute $\mathcal R_0$ for this model. This method splits the Jacobian matrix of the linearized system at the disease free equilibrium into two matrices: a new infection matrix $F$, whose $ij$ entry is the rate of new infections coming from fully susceptible classes into infected class $i$ caused by individuals in infected class $j$; and a transition matrix $V$, whose $ij$ entry is the rate of individuals from infected class $i$ entering infected class $j$. Then the basic reproduction number $\R_0=\rho(FV^{-1})$ where $\rho(\cdot)$ denotes the spectral radius. To apply this method, we need to identify the infected classes and the fully susceptible classes. Following Lindquist et al.~\cite{LMDW_JMB2010}, we choose $S_{k0}$, $k=0,1,2,\ldots,M$ to be the fully susceptible classes, while $S_{si}$ for all $i\geq1$ are deemed ``infected'' classes as they behave somewhat like exposed classes, because they have at least one infectious neighbor. All $I_{si}$ are infected classes. For the household classes, $H_{n0}$ are fully susceptible classes, while $H_{nj}$ for all $j\geq1$ are infected classes.

At the disease free equilibrium (DFE), $S_{k0}=N_k$, $S_{si}=0$ for all $i\geq1$, $H_{n0}=H_n$ being the household size distribution, and $H_{nj}=0$ for all $j\geq1$. Individuals from the fully susceptible classes $S_{k0}$ come into infected classes $S_{k-1,1}$ from external infection of a neighbor with rate \[ \sum_{k=0}^{M}\sum_{s+i=k}\beta_e si S_{si} \frac{kS_{k0}}{\sum_{k=0}^{M}\sum_{s+i=k} s S_{si}}\,, \] and from within-household infection of a neighbor with rate \[ \sum_{n=1}^m\sum_{j=1}^n\beta_h (n-j)jH_{nj} \frac{kS_{k0}}{\sum_{k=0}^{M}\sum_{s+i=k} S_{si}}\,, \] Individuals from the classes $S_{k0}$ come into an infected classes $I_{k0}$ from within-household infection with rate \[ \sum_{n=1}^m\sum_{j=1}^n\beta_h (n-j)jH_{nj} \frac{S_{k-1,1}}{\sum_{k=0}^{M}\sum_{s+i=k} S_{si}}\,. \] Individuals from $H_{n0}$ come into infected classes $H_{n1}$ with rate \[ \sum_{k=0}^{M}\sum_{s+i=k}\beta_e iS_{si}\frac{nH_{n0}}{\sum_{n=1}^m\sum_{j=1}^n(n-j)H_{nj}} \,. \] These terms are the only terms that cause fully susceptible individuals to move to infected classes. We order the variables of the infected classes in the following order: $S_{01}$, $S_{11}, S_{02}$, …, $I_{00}$, $I_{01}$, $I_{11}$, $I_{02}$, …, $H_{11}$, $H_{21}$, $H_{22}$, …, i.e., order the variables as $S$, $I$, $H$, and then by degree for $S$ and $I$ and household size for $H$, then by the number of infected neighbors for $S$ and $I$, and the number of infected household members for $H$. Then, at DFE, the matrix of new infections $F$ has the following block form, \begin{equation} F=\left[ \begin{array}{ccc} F_{SS} & 0 & F_{SH} \\ 0 & 0 & 0 \\ F_{HS} & 0 & 0 \end{array} \right]\,. \label{eq:F} \end{equation} The entry in row $S_{k-1,1}$, column $S_{uv}$ of $F_{SS}$ is \begin{equation} F_{SS}(S_{k-1,1}, S_{uv}) = \beta_e uv \frac{kS_{k0}}{\sum_{u=1}^{M}uN_u}\,, \label{eq:F:SS:1} \end{equation} and $0$ for all other rows (i.e., rows $S_{si}$ with $i\geq2$); the entry in row $S_{k-1,1}$, column $H_{nj}$ of $F_{SH}$ is \begin{equation} F_{SH}(S_{k-1,1}, H_{nj}) = \beta_h (n-j)j\frac{kS_{k0}}{\sum_{u=0}^{M}N_u}, \label{eq:F:SH:1} \end{equation} and $0$ for all other rows; the entry in row $H_{n1}$, column $S_{si}$ of $F_{HS}$ is \begin{equation} F_{HS}(H_{n1}, S_{si}) = \beta_e i\frac{nH_{n0}}{\sum_{n=1}^mnH_{n0}} \,, \label{eq:F:HS:1} \end{equation} and $0$ for all other rows.

Since the $F$ matrix has rows and columns of zeros for $I_{si}$ for all $s$ and $i$, and $H_{nn}$ for all $1\leq n\leq m$, these variables do not contribute to $\R_0$ \cite{DHR_JRSI2009}. Thus they can be removed when computing $\R_0$.

For the remainder of the state variables, the transition matrix $V$ is block diagonal $V=V_{SS}\bigoplus V_{HH}$ with $V_{SS}$, the block that corresponds to $S_{si}$, being block bidiagonal. Each diagonal block $V^S_{kk}$ corresponds to a degree $k$ with $V^S_{kk}=\diag\{i(\beta_e+\gamma)\}$. The off-diagonal block $V^S_{k,k+1}$ has a column of zeros followed by a diagonal matrix $\diag\{-i\gamma\}$. Also, $V_{HH}$, the diagonal block of $V$ that corresponds to $H_{nj}$, is block super-bidiagonal, with the diagonal blocks $V^H_{nn}$ corresponding to a household size $n$. For $n\geq 2$, the block $V^H_{nn}$ is sub-bidiagonal, with the $j$-th diagonal entry equal to $\beta_h(n-j)j+\gamma j$, and the $(j+1,j)$ entry below it equal to $-\beta_h(n-j)j$, for $1\leq j\leq n-1$. The off-diagonal block $V^H_{n,n+1}$ has a zero column followed by $\diag\{-2\gamma, \ldots, -n\gamma\}$ for the entry corresponding to $H_{nj}$.

To further simplify the computation of $\R_0$, we need the following lemma; see the proof in Appendix A.

\noindent{\bf Lemma} If the new infection matrix $F$ has a zero row, say, row $i$, and the inverse of the transition matrix $V^{-1}$ has all zeros in row $i$ except for the main diagonal entry, then the $i^{\text{th}}$ variable does not contribute to $\rho(FV^{-1})$.

For all $i\geq2$, variables $S_{M-i,i}$ correspond to rows of zeros in $F$, and to rows of zeros in $V^{-1}$ except for the main diagonal. Thus, from the above lemma, these variables do not contribute to $\R_0$ and so can be eliminated when computing $\R_0$. Variables $S_{si}$ can thus be iteratively removed for all degrees $k=s+i$ and $i\geq2$. We now use $F$ and $V$ to denote the reduced matrices with variables $I_{si}$, $H_{nn}$, and $S_{si}$ for $i\geq2$ eliminated. The elimination leaves a diagonal $V_{SS}=(\beta_e+\gamma)I_d$, where $I_d$ is the identity matrix or order $M$; and $F_{SS}$ has rank $1$, \begin{equation} F_{SS}=\frac{\beta_e}{\sum_{u=1}^{M}uN_u}a_1b_1^T \label{eq:F:SS} \end{equation} where $a_1=[N_1, \ldots, kN_k, \ldots, MN_M]^T$, and $b_1^T=[0, 1, 2, \ldots, k-1, \ldots, M-1]$. Also $F_{SH}$ has the form \begin{equation} F_{SH}=\frac{\beta_h}{\sum_{u=0}^{M}N_u}a_2b_2^T \label{eq:F:SH} \end{equation} where $a_2=a_1$, and $b_2^T=[(n-j)j]$ for $2\leq n\leq m$ , $1\leq j\leq n - 1$ (note that $H_{nn}$ are eliminated). The block $F_{HS}$ has the form \begin{equation} F_{HS}=\frac{\beta_e}{\sum_{n=2}^{m}nH_{n0}}a_3b_3^T \label{eq:F:HS} \end{equation} where $a_3=[c_2,c_3, \ldots, c_m]^T$, $c_n$ is a row vector of order $n-1$ (again, because $H_{nn}$ are eliminated) with first term $nH_{n0}$ followed by zeros, and $b_3^T$ is the row vector of ones of order $M$. Hence, $F$ has rank $2$, and so does $FV^{-1}$.

To find the spectral radius $\rho(FV^{-1})$, because $FV^{-1}$ has rank $2$, its characteristic equation is \begin{equation} \label{eq:characteristic} \chi(x) = x^2-\text{Tr}(FV^{-1})x+\vert FV^{-1}\vert_2 =0\,, \end{equation} where $\vert\cdot\vert_2$ is the sum of all the $2\times2$ minors. Note that \[ \text{Tr}(FV^{-1}) = \text{Tr}\left(F_{SS}V_{SS}^{-1}\right)=\frac{\beta_e}{\beta_e+\gamma} \frac{\sum_{k=1}^M k(k-1)N_k}{\sum_{k=1}^M kN_k}=\frac{\beta_e}{\beta_e+\gamma} \left(\frac{\avg{k^2}}{\avg{k}}-1\right) \,, \] where $\avg{k}$ is the average degree of the contact network, and $\avg{k^2}$ is the second moment of the degree distribution. Because $F_{SS}V_{SS}^{-1}$ has rank $1$, $\vert F_{SS}V_{SS}^{-1}\vert_2=0$. Thus, \[ \vert FV^{-1}\vert_2 = -\sum_{k=1}^M\sum_{n=2}^m\sum_{j=1}^{n-1} A_{njk}B_{knj} \] where $A_{njk}$ is the row $H_{nj}$, column $S_{k-1,1}$ entry of the matrix $F_{HS}V_{SS}^{-1}$, and $B_{knj}$ is the row $S_{k-1,1}$, column $H_{nj}$ entry of the matrix $F_{SH}V_{HH}^{-1}$. Note that, from \eqref{eq:F:SH} and \eqref{eq:F:HS}, both matrices have rank $1$ \[ A_{njk}=\begin{cases} \frac{\beta_e}{\beta_e+\gamma}\frac{nH_{n0}}{\sum_{n=1}^{m}nH_{n0}}\,, & j=1 \,,\\ 0\,, & j\geq2\,. \end{cases} \] and \[ B_{knj}=\frac{\beta_h}{\sum_{k=1}^{M}N_k}kN_kb_2^TV_{HH,nj}^{-1}, \] where $V_{HH,nj}^{-1}$ is the $H_{nj}$ column of the matrix $V_{HH}^{-1}$. Since $V$ is a nonsingular M-matrix, $V^{-1}_{HH,nj}$ has non-negative entries. Note that the total population $N=\sum_{n=1}^mnH_{n0}=\sum_{k=0}^MN_k$, giving \[ \vert FV^{-1}\vert_2 = -\frac{\beta_e}{\beta_e+\gamma}\frac{\avg{k}}{N}\sum_{n=2}^m \beta_hnH_{n0}b_2^TV_{HH,n1}^{-1} \] Thus, the characteristic equation \eqref{eq:characteristic} has a unique positive root, $\rho(FV^{-1})=\R_0$, which can be calculated explicitly as the largest root of \eqref{eq:characteristic}. From \cite{DW_MB2002, DHM_JMB1990}, $\R_0=1$ is a threshold between local stability of the disease free equilibrium ($\R_0<1$) and instability ($\R_0>1$).

To determine an equivalent, simpler form of this threshold, note that \eqref{eq:characteristic} has a unique positive root, which is less than $1$ if and only if \[ \chi(1) = 1-\text{Tr}(FV^{-1})+\vert FV^{-1}\vert_2 > 0\,, \] that is, if and only if \begin{equation} \label{eq:R0} \R = \frac{\beta_e}{\beta_e+\gamma} \left(\frac{\avg{k^2}}{\avg{k}}-1\right) + \frac{\beta_e}{\beta_e+\gamma}\frac{\avg{k}}{N}\sum_{n=2}^m \beta_hnH_{n0}b_2^TV_{HH,n1}^{-1} < 1 \,. \end{equation} The first term in \eqref{eq:R0} is the basic reproduction number of the network SIR model with no household structure \cite{LMDW_JMB2010, Newman_PRE2002}; the second term is the reproduction number of the following household model: \begin{linenomath} \begin{align} H_{nj}' = & \frac{\beta_e}{\beta_e+\gamma}\avg{k}\sum_{p=1}^m\sum_{q=1}^p\beta_h(p-q)qH_{pq} \frac{(n-j+1)H_{n,j-1}-(n-j)H_{nj}}{\sum_{p=1}^m\sum_{j=0}^n(p-j)H_{pj}} \nonumber \\ & +\beta_h[(n-j+1)(j-1)H_{n,j-1}-(n-j)jH_{nj}] \nonumber \\ & +\gamma\left[(j+1)H_{n+1,j+1}-jH_{nj}\right]\,. \label{eq:R0:household:equiv} \end{align} \end{linenomath} Each household infection on average brings $\displaystyle \frac{\beta_e}{\beta_e+\gamma}\avg{k}$ infections in other households, randomly distributed among the households. Thus, the threshold $\R$ is the sum of the number of secondary external infections and the number of secondary household infections (caused by itself and all secondary external infections) caused by a typical infectious individual in a completely susceptible population.

As $\beta_h\to\infty$, $V_{HH}$ is block diagonal, with each diagonal entry $(V^H_{nn})^{-1}$ being lower triangular. Specifically, the nonzero entries of row $H_{nj}$ of $(V^H_{nn})^{-1}$ are all $1/[\beta_h(n-j)j]$ with $2\leq n\leq m$ and $1\leq j\leq n-1$. Together with the definition of $b_2$, \[ b_2^TV^{-1}_{HH,n1}=\frac{n-1}{\beta_h} \,, \] for all $2\leq n \leq m$. Thus, as $\beta_h\to\infty$ \begin{equation} \label{eq:R0:limit} \R = \frac{\beta_e}{\beta_e+\gamma} \left(\frac{\avg{k^2}}{\avg{k}}-1\right) + \frac{\beta_e}{\beta_e+\gamma}\frac{\avg{k}}{N}\sum_{n=2}^m n(n-1)H_{n0} \,. \end{equation} As shown in Appendix B, this is the percolation threshold \cite{LMDW_JMB2010, Newman_PRE2002} of the household network in which nodes are households and links are the inter-household links.

\subsection{Dependency of the disease thresholds on parameters}

From \eqref{eq:characteristic} and \eqref{eq:R0}, $\R$ and $\R_0$ both are increasing function of $\beta_e$, and it can be shown that they are also increasing functions of $\beta_h$, and decreasing functions of $\gamma$. The monotonicity with $\beta_e$ and $\beta_h$ is illustrated in Figure \ref{fig:individuals}, as $\R_0$ (and hence $\R$) is an increasing function of the initial growth rate.

Volz et al. \cite{VMGM_PCB2011} found that, while keeping the total degree distribution constant, clusters (and thus household structure) always slow down the epidemic, and reduce the epidemic size. Thus, the increase in epidemic size and growth rate that we found may be the effect of the change in the total degree distribution. To further check the effect of households on disease dynamics, stochastic simulations were carried out to compare the spread of the disease on a network with households to the spread on a network without households having a degree distribution identical to that of the household network. The results of the simulations are depicted in Figure \ref{fig:config}, where it can again clearly be seen that the household structure slows down epidemics, and reduce the peak sizes.

The dependence of the epidemic curve on the household size is further illustrated in Figure \ref{fig:randhh}. For this figure, rather than the average total degree $\langle k_T \rangle$, the \textit{external} degree $\langle k \rangle$ was used. This has no impact on the network construction but does imply a higher degree average for nodes (i.e., for household sizes greater than one, $\langle k_T \rangle > \langle k \rangle$). The inset shows the dependence of the final size on the average household size. As expected, due to the full intra-household connectivity, the final size increases for increasing household sizes. For comparison, a household size of one is included giving a curve representing the spread of a disease on an Erd\H{o}s–R\'enyi type random network.

Intuitively, it seems reasonable that given the same transmission rate of a disease, larger households will lead to an epidemic curve characteristic of higher transmission rates. Indeed, this is the explanation given for the (on average) early onset of symptoms of varicella-zoster infection in Guinea Bissau where the infectivity of the virus is lower than in temperate climates yet the disease spreads to children at a similar age. As can be seen in Figure \ref{fig:bebh}, this is also the case in the household model. By increasing the average household size from 3 to 12, a drop in the transmission rate from $\beta_e = \beta_h = 0.15$ to $\beta_e = \beta_h=0.08$ yields nearly the same epidemic curve.

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

It is well known that the contact network in the human population is neither homogeneous nor random. While a random network is clearly closer to the actual situation than a homogeneous network, many other important structures, such as households, workplaces and schools, can not be represented by such random networks.

When modeling diseases, it is not obvious that such differences are of great importance. However, it has long been realized that certain diseases mainly spread locally, for example hand-foot-mouth disease often spreads through kindergartens. In order to understand the disease dynamics, it is therefore imperative to investigate the significance of local structures.

Households, workplaces and schools all share in common a substantially different contact network structure from the contact network in the general population in that they are either completely homogeneous (households) or to varying degrees highly homogeneous (schools, workplaces). However, even in schools and workplaces, there are many homogeneously linked groups with contacts similar to those of households. It is therefore a good (first) approximation to consider a hybrid contact network consisting of households that have external Erd\H{o}s–R\'enyi type random links.

In this paper, we have developed an ODE model for the spread of an SIR type disease on such a network and demonstrated its validity with the help of detailed stochastic simulations. It is found that under generic circumstances, households have a drastic impact on disease dynamics (it is notable that recent field work on the varicella-zoster virus in Guinea Bissau confirmed this rather dramatically).

Our model is more complex that the Volz et al model \cite{VMGM_PCB2011}, which also yields excellent agreement with stochastic simulations as our model does. However, with our model, we have been able to find analytic disease threshold conditions. Not surprisingly, increasing inter- or intra-household transmission rates, and decreasing the recovery rate (increasing the mean infectious period) all have the effect of making the disease transmit faster and invade more easily,

In addition, because of the nature of effective degree models, we can easily extend our model to Susceptible-Infectious-Susceptible (SIS) type diseases (where the acquired immunity is negligible, which is true for many bacteria infectious); see \cite{LMDW_JMB2010} for details. The disease threshold conditions for an SIS model will be different from the SIR model.

Thus our household model can be the solid foundation for computational and analytical studies for disease control in realistic networks with household or workplace structures.

\vspace{5mm} \noindent \textbf{Acknowledgments}

\noindent This work is partially supported by NSERC Discovery Grants (JM, PvdD) and MITACS (PvdD).

\appendix

\section{Proof of the Lemma} Let $F=[f_{ij}]$, and $V^{-1}=[w_{ij}]$. Then, the $uv$ entry of $FV^{-1}$ is \begin{equation} \label{eq:FV-1:ij} \left(FV^{-1}\right)_{uv}=\sum_k f_{uk}w_{kv} \,. \end{equation} If $u\neq i$ and $v\neq i$, then $f_{ik}$, $f_{ki}$, $w_{ik}$, and $w_{ki}$ only appear in the terms $f_{ui}w_{iv}$, which are zero because $w_{iv}=0$. Thus, $\left(FV^{-1}\right)_{uv}$ does not contain $f_{ik}$, $f_{ki}$, $w_{ik}$, and $w_{ki}$ for all $k$, and $u\neq i$ and $v\neq i$. If $f_{ij}=0$ for all $j$, then $FV^{-1}$ has all zeros in the $i^{\text{th}}$ row, and so the $i^{\text{th}}$ row and column of $FV^{-1}$ does not contribute to $\rho(FV^{-1})$. \qed

\section{The household network} As $\beta_h\to\infty$, any infected household member causes the entire household to be immediately infected. Thus all the members of a household are effectively coupled as a single node. In this limit, each household can be treated as a single node (individual), with its degree equal to the sum of the external degrees of all its members. We call this network the \emph{household network}. The disease threshold on the household network thus depends on the first and second moments of its degree distribution, denoted as $\avg{\bar k}$ and $\avg{\bar k^2}$. to be distinguished from the moments of the external degree distribution $\avg{k}$ and $\avg{k^2}$. At the disease free equilibrium, the first and second moments of the household size distribution are \[ \avg{n}=\sum_{n=1}^mn\frac{H_{n0}}{H} = \frac{N}{H}\,, \] and \[ \avg{n^2}=\sum_{n=1}^mn^2\frac{H_{n0}}{H} \,, \] where $H=\sum_{n=1}^mH_{n0}$ is the total number of households.

Let $k_i$ be the external degree of an individual in a household size of $n$, and let $K_n=k_1+k_2+\cdots+k_n$ be the total external degree of individuals in the household, which is a random variable of the household size $n$. Then, \[ \avg{\bar k} = E[K_n] = E\{E[k_1+k_2+\ldots k_n|n]\}=E[n\avg{k}]=\avg{k}\avg{n} \,. \] Also \begin{align*} \avg{\bar k^2}& = E[K_n^2] \\ & = E\{E[(k_1+k_2+\ldots+ k_n)^2|n]\} \\ & = E[n\avg{k^2}+(n^2-n)\avg{k}^2] \\ & = \avg{n}\avg{k^2}+(\avg{n^2}-\avg{n})\avg{k}^2\,. \end{align*} Thus, on the household network with inter-household transmission rate $\beta_e$, the percolation threshold \cite{LMDW_JMB2010, Newman_PRE2002} is \[ \R_0 = \frac{\beta_e}{\beta_e+\gamma}\frac{\avg{\bar k^2} -\avg{\bar k}}{\avg{\bar k}} \,. \]

On the other hand, \eqref{eq:R0:limit} becomes \begin{align*} \R & = \frac{\beta_e}{\beta_e+\gamma} \left(\frac{\avg{k^2}-\avg{k}}{\avg{k}} + \frac{\avg{k}}{N}\sum_{n=2}^m n(n-1)H_{n0} \right)\\ & = \frac{\beta_e}{\beta_e+\gamma} \left(\frac{\avg{k^2}-\avg{k}}{\avg{k}} + \frac{\avg{k}(\avg{n^2}-\avg{n})}{\avg{n}} \right)\\ & = \frac{\beta_e}{\beta_e+\gamma} \frac{\avg{\bar k^2}-\avg{\bar k}}{\avg{\bar k}} \end{align*} by the above formula for $\avg{\bar k}$ and $\avg{\bar k^2}$. Thus, as $\beta_h\to\infty$, \eqref{eq:R0:limit} is the same as the percolation threshold $\R_0$ for the household network. \bibliographystyle{plain} \bibliography{household}

\clearpage