body.tex

\jd{Can we drop “comprehensive” from the title? What about “modular construction of \ldots”?} \steve{agreed. done.}

\section{Introduction}\label{intro}
The COVID-19 pandemic has reemphasized the importance of compartmental epidemic models \citep{abou2020compartmental, massonis2021structural, adam2020special, currie2020simulation, lofgren2014mathematical, mcbryde2020role, enserink2020covid} and has resulted in a flood of new compartmental models (e.g., \cite{friston2020dynamic, fields2021age, chang2022stochastic, lavielle2020extension, balabdaoui2020age, leontitsis2021seahir, LeviEarn22}). 
This abundance of new model variants is expected given the number of public health modelers trying to use scientific understanding of emerging infectious diseases to contribute to public policy. Modelers must be able to build models rapidly to explore scenarios and generate high quality forecasts; public health recommendations have the biggest impact if they can be acted on promptly. However, the speed at which modelers can develop new models typically trades off with model quality. We therefore need tools that allow modelers to build models more quickly without sacrificing quality. 

One approach to this speed-quality trade-off is to build infectious disease models incrementally. Information is scarce early in an epidemic, and so early models should be simple to reflect ignorance. As epidemics progress, we learn more about the characteristics of the pathogen and its transmission; at the same time the public health landscape becomes clearer but more complex. Because policy choices require fast input from scientists, modelers need to add complexity to their models quickly if they are to be relevant to policy. 

Proceeding in this way eventually results in complex and fragile models, much of whose complexity is no longer relevant. This sort of complexity also makes it harder to add additional features to the model. Therefore, modelers need tools that make it easier to flexibly add and remove model structure.

\cite{savageau1988introduction} and \cite{voit1988recasting, voit1990s} made an early attempt to create such a toolbox by recasting the underlying differential equations of a model into a canonical form they call an ``S-model". Unfortunately this effort focused on the model's differential equations rather than its graphical structure, thus making it unsuitable for less mathematically inclined modelers. It does not seem to have been widely adopted.

\cite{friston2020dynamic} describe how the state space of a complex epidemiological model can be constructed from the product of different latent state dimensions (their Figure 1 shows an example with infection status, clinical status, testing status, and location), but the definition of which compartments are connected, and the rates of flow between them, is left up to the modeller.

A promising recent project to formalize the construction of compartmental models employs the language of category theory \citep{fong2018seven, Libkind2022an, libkind2021operadic, baez2023categorical, baez2017compositional}. This powerful approach addresses many of the concepts we discuss here; however, at its current stage of development it requires considerable knowledge of advanced mathematics to use effectively.  An ongoing project to implement the category theoretic approach in the Julia language can be found at \url{https://github.com/AlgebraicJulia/AlgebraicPetri.jl} \citep{algebraicjulia}. 

\cite{worden2017products} use the relatively simple language of graph theory to describe common methods of “multiplying” a set of compartmental models (\term{factor models}) into a new \term{product model} that incorporates the structure of all of its factors. The current paper is a result of the our efforts to build software implementations of the products described by Worden and Porco. 

\paragraph{Procedure for Model Multiplication}\label{genproc}

We view model multiplication as a four-step procedure:

\begin{enumerate}
    \item \textbf{Generate vertices} of the product model by combining the vertices of the factor models using a Cartesian product (i.e., generate new vertices by choosing one vertex from each factor model, e.g., old susceptible).
    \item \textbf{Generate edges} of the product model, typically by taking the Cartesian product of edges in each factor model with the vertices in the other (e.g., we would have an edge from old-susceptible to old-infected). In some cases, we may add additional edges, to allow state changes in multiple strata to occur simultaneously (e.g., a direct flow from medium susceptible to old infected).
	\item \textbf{Simplify} by first removing vertices we wish to disallow (e.g., old people with maternal antibody protection), with associated edges, and then removing any additional edges that we wish to disallow.
    \item \textbf{Adjust flows} by specifying how flow functions should accommodate the presence of additional strata (e.g., spatial structure constrains infection processes).
\end{enumerate}

The graph generated by using the default Cartesian products in steps 1 and 2 is sometimes called the “box product” (apparently due to Harary 1969, needs to be chased down).

Worden and Porco pointed out that flows may need to be adjusted (step 4); for example, when combining a standard SIR model with some other form of structure, we need to decide whether the susceptible population of any particular stratum can be infected by the infectious populations of other strata.
In other cases we might prohibit cross-stratum infection, for example by allowing infection within but not between geographic regions. 

\section{What is a Compartmental Model?}\label{dcm}

\FloatBarrier %% Why?
\begin{figure}
    \centering
    \includesvg[inkscapelatex = true, width=\textwidth]{images_redux/SIR_and_Age.svg}
    \caption{A standard SIR model and a simple age stratification model. Orange shading denotes the infectious compartment, while blue shading denotes non-infectious compartments.}
    \label{fig:sir_age_strat}
\end{figure}

\begin{figure}
    \centering
    \includesvg[inkscapelatex = true, width=\textwidth]{images_redux/simple_SIR_age.svg}
    \caption{The product of the two models shown in Figure~\ref{fig:sir_age_strat}. \jd{These two panels should be combined into one Figure for readability, I think.}}
    \label{fig:simple_sir_age}
\end{figure}

Figure~\ref{fig:sir_age_strat} shows two factor models: the ``Epi" model represents individuals in a population being infected and then recovering from some infectious disease; the ``Age" model represents individuals aging. Our goal is to design an algorithm to generate an appropriate product of these two factors (Figure~\ref{fig:simple_sir_age}). Each compartment in the product model inherits labels from one compartment in each of the factor models (e.g. SY = ``susceptible, young''); the flows also typically have two labels (e.g. old infection flow). 

\paragraph{Infection flows and infectious compartments}

All of the flows in the Age model are constant \pc flows; that is, the rate of flow between any two compartments is strictly proportional to the number of individuals in the ``from'' compartment. 
The flow from ``infected" to ``recovered" in the Epi model is also a constant \pc flow. In contrast, the \pc flow from Susceptible to Infected is typically proportional to the number of Infected -- thus, the \pc rate depends on occupancy of another compartment, indicated by the dashed arrow in Figure~\ref{fig:sir_age_strat}. In the specific context of compartmental models for epidemiological dynamics, \pc flows that are state-dependent will generally relate to the infection (transmission) process; we will call these flows ``infection'' flows, and call the compartments that determine their rates ``infectious'' compartments.

In the SIR model the infectious compartment is the ``to" compartment of the infection flow; this is not true for all models.  Models may also have multiple infectious compartments --- for example the infected population maybe separated on the basis of symptom severity. Figure \ref{fig:simple_sir_age} does not include dashed arrows connecting infection flows with infectious compartments; these relationships will be discussed in Section \ref{worden}.

\paragraph{Sources and sinks}

Consider the long-term outcome of the ``Age'' model in Figure \ref{fig:sir_age_strat}. Since the population is closed and people only flow from younger compartments to older ones, the entire population will eventually accumulate in the ``old" compartment. If we want to simulate the system over time scales comparable to the host lifespan, we need to add and remove people from the model (``vital dynamics''). This can be done with ``sources" (flows with no ``from'' compartment'') and ``sinks" (flows with no ``to" compartment). We could add a source flow into the ``young" compartment and a sink flow out of the ``old" compartment (or out of every compartment), representing the effects of birth and death. 

Sources are also useful in models of viral presence in wastewater. The rate at which new viral particles enter the wastewater compartment is proportional to the number of infected people in the population. But infected people don't \emph{become} contaminated wastewater; instead, we add a source whose flow rate depends on the infected population.

Multiplying factor models that have sources and sinks may require extra decisions. In most models incorporating vital dynamics, for example, newborns entering the youngest age class are assumed to be susceptible \citep{Earn2008}, but in cases where vertical transmission or maternal immunity are possible, we might need to consider other possibilities. 

\paragraph{Absolute and \pc flows}

Most of the flows discussed above are \pc flows, where the total flow between compartments is given by a \pc rate multiplied by the population in the ``from" compartment
(for infection flows, the \pc \emph{rates} depend on the infectious compartments).
In contrast, \term{absolute flows} are specified in terms of the total flow. 
Source flows, which have no from compartment, can only be implemented as absolute flows (this is how births are typically implemented \citep{Earn2008}). If population-level vaccination rates are recorded in public health data, it may be simplest to implement these as absolute flows (although specifying an absolute flow from a non-source population may sometimes lead to a mathematically ill-posed model \citep{gharouni2022testing}).

\paragraph{Implementation}

To define a compartmental model, we:

\begin{enumerate}
\item specify compartments (nodes) and flows (edges). A flow is specified by defining: from/to compartments (one of which will be empty for a source or sink flow); a rate type (\pc or total); and a rate function that describes how the rate depends on parameters, time and other flows. 

\item compute the rate vector \ratevector, based on parameters, time and the state vector \statevector. The length of the rate vector is the number of flows in the model.

\item compute the flow vector \flowvector\ by multiplying each element of the rate vector by a value from the state vector (for \pc rates) or by 1.

\end{enumerate}

It will sometimes be convenient to have different elements of the flow vector have the same to and from compartments; for example if individuals in one compartment may be subject to infection by individuals in several different compartments. Such flows will typically be added together downstream.

\jd{Above are the changes that I already wanted to make... I need to try now to read the rest before changing anything. I switched the definitions of \flowvector\ and \ratevector\ above only, so watch out for that below.}

\jd{Need to add a \P\ about discrete vs.~continuous and Markovian vs.~deterministic downstream of the flow vectors.}

\subsection{Parameterizing product models}

In addition to understanding how to construct compartments and flows of a product model from the compartments and flows of its factors, we would like to be able to compute parameter values for the product model using the known values from the factor models.
Product model parameters are often related to parameters in the original model factors in simple mechanistic ways. However, there is an enormous range of possible relationships between the parameters of the factor models and the parameters of their product.
Some parameters, such as those relating to aging, may be constant across all strata of a product model. Others, such as recovery time, may be constant with respect to some strata (e.g., location) but variable with respect to others (e.g., age). In some cases, we may want to simplify relationships due to data constraints --- for example, we may choose to model recovery time as constant across ages for a particular purpose. Ultimately, the question of how factor-model parameters should be generalized to the product model depends on the goals of the model. 

If a parameter describes processes that are purely host-dependent it will be constant across strata for different strains of the pathogen. The converse could also occur, if there are parameters that are purely pathogen-dependent. 
In other cases, the value of a parameter at each stratum in the product model may be derived from factor-model parameters by a simple rule. For example, it might be useful to parameterize the recovery time of hosts in different strata as proportional to some baseline value (e.g., the value in healthy adults). If $\alpha$ is the parameter value in the factor model and $\betavec$ denotes the values of the derived parameters at each different strata in the product model then $\betavec = \alpha\vec{w}$ where $\vec{w}$ is a vector of weights.

Relationships among factor-model parameters can be constrained in other ways. When multiple flows emanate from a single compartment, we may want to parameterize them as a partition of the total outflow. For example, an exposed host may have multiple possible fates such as asymptomatic infection, mild infection, or severe symptoms. In that case the factor model in question will have three parameters ($\alpha_i \in (0,1), i = 1 \ldots 3$) that sum to one ($\sum \alpha_i = 1)$. Every stratum of the corresponding product model will have three parameters derived from the original $\alpha$ values,  but each stratum may have different partitions. For example people in different age groups may be more or less likely to experience severe, mild, or no symptoms. 

Different strata of a product model may interact.
A standard formulation for the force of infection of a pathogen (i.e., the per-susceptible rate of infection) is
$\Lambda = \frac{\beta I}{N}$. Suppose we now stratify this model to represent a scenario where each person lives in one of three different locations but may come in contact with anyone in any location. Our model would then have three infectious proportions ($\vec{\chi} = (I_1/N_1, I_2/N_2, I_3/N_3)$) and three susceptible compartments ($\vec{S} = (S_1, S_2, S_3)$.
Using the standard formulation, we could then calculate the force of infection as:
$$
    \vec{\Lambda} =  B \vec{\chi}, 
$$
where $B$ is a $3 \times 3$ matrix of transmission parameters $B = \left(\beta_{ij}\right)$.

The generic version of this approach makes the number of required parameters increase quadratically with the number of locations. Epidemiological modelers have devised many ways to add structure $B$ to reduce the required number of parameters \cite{andemay85, andemaybook, grenande85}. 
\jd{Need recent cites here.}
For example, the likelihood of a person residing in one location coming into contact with a person somewhere else may be assumed to depend only on the distance between the two locations, allowing the construction of large $B$ matrices with a relatively small number of parameters. 

There are, of course, other ways to handle this kind of parameter simplification. Many situations will allow for a mapping that relates the default parameter space generated by model products to a smaller parameter space based on straightforward assumptions. \steve{Added the following to this paragraph because it is how I think about $\beta$ in product models.} One useful approach is to decompose $\beta_{ij}$ into susceptibility, $\sigma_i$, and infectivity, $\tau_j$, parameters.
$$
    \beta_{ij} = \sigma_ip_{ij}\tau_j,
$$
where $p_{ij}$ is the rate of contact between individuals in $i$ and $j$ \steve{Should $p_{ij}$ be considered part of the parameter decomposition or, as I prefer, as additional `data' that help put the parameters on a more interesting scale?}. If $i$ and $j$ represent age groups, these parameterizations could distinguish between the effects of age on non-contact processes that impact infectivity (e.g., masking \steve{something better here?}) versus those that impact susceptibility (e.g., age-dependent immunity). \steve{Balancing contact rates stuff}

\section{Cartesian Model Products}\label{worden}

Returning to Figure \ref{fig:sir_age_strat}, recall that the dashed arrow in the SIR model indicates that the flow from $S$ to $I$ is an infected flow and that its magnitude is partly determined by the number of infected people. When we combined the SIR and age models in Figure \ref{fig:simple_sir_age} we omitted dashed arrows. \cite{worden2017products} describe two separate products, each with the dashed lines drawn differently. In their \term{naive product}, each susceptible group can only be infected by infectious people in the same group (Figure \ref{fig:naive_product}). In their \term{modified  product}, each susceptible group can be infected by infectious people in any group (see Figure \ref{fig:modified_product}). 

These two models have the same graph, given by the box product of the factor models. They differ only in the form of the infection flows. The graph underlying a product model is not always the box product (or a subset), though: Worden and Porco's \term{strong product} adds additional edges. In the model illustrated in Figure \ref{fig:simple_sir_age}, individuals cannot move directly from the young susceptible compartment to the middle infected one; to make that transition they must take one step first and then the other. If we have strata corresponding to different pathogens, we might use a strong product to allow for the possibility of infection with multiple pathogens in a single event. Strong products might also be preferred for implementations with discrete time steps, where the probability of two independent events in the same time step cannot be neglected. 

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Age_stratified_SIR_Naive.svg}
    \caption{The naive product of the two models from Figure \ref{fig:sir_age_strat}. Blue denotes non-infectious compartments, yellow/orange/red denote infectious compartments. The force of infection is only influenced by the infected population within the same age stratum. In this example, people of different age groups have no contact (or very limited contact) with each other.}
    \label{fig:naive_product}
\end{figure}

\begin{figure}
    \centering
    \includesvg[    width=\textwidth]{images_redux/Age_stratified_SIR_Modified.svg}
    \caption{The modified product of the two models from Figure \ref{fig:sir_age_strat}. Unlike in Figure \ref{fig:naive_product}, individuals make epidemiological contacts across age strata, so the force of infection for each age stratum is influenced by the infected population in all age strata.}
    \label{fig:modified_product}
\end{figure}

\FloatBarrier

Age-stratification will typically use a modified product, because people of all ages commonly interact with each other. Spatial stratification can use different approaches. When movement of individuals is simulated explicitly, \cite[e.g.][]{mohammadi2023importation}, it makes sense to use a naive product, where we assume that people in different locations do not physically interact. Models without explicit movement, \cite[e.g.][]{dietz1995structured}, are more suited to a modified product, which would allow people based in different locations to have some amount of contact. 

The naive product restricts interaction to people in the same stratum; the modified product allows people in any strata to interact. In many cases, it will be computationally efficient to consider an intermediate “generalized” approach, that allows for people in each stratum to interact with people in “neighboring” strata.

Below we show three different ways an SI model could be stratified with location. Figures \ref{fig:spat_n} and \ref{fig:spat_m} show the naive and modified products respectively. Figure \ref{fig:spat_g} shows a generalized product where interactions can only occur within a single geographic region or between neighboring regions. For example, an infected person in the Toronto region could infect a susceptible person in Toronto or Ottawa but not one in Montreal, while the neighbourhood of an infected person in Ottawa includes all three regions.
This model is equivalent to, but more efficient than, the modified product with $\beta_{31}$ and $\beta_{13}$ set to zero. 

\FloatBarrier

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Spatial_stratified_SI_Naive.svg}
    \caption{The naive product of an SI model with location model including Toronto, Ottawa, and Montreal. The force of infection at any given location depends only on the infectious population at that location.}
    \label{fig:spat_n}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Spatial_stratified_SI_Modified.svg}
    \caption{The modified product of an SI model with location model including Toronto, Ottawa, and Montreal; the forces of infection depend on infectious populations in all three regions.}
    \label{fig:spat_m}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Spatial_stratified_SI_General.svg}
    \caption{A generalized product of an SI model with location model including Toronto, Ottawa, and Montreal. The forces of infection depend on the infectious population at the same location and in neighboring locations.}
    \label{fig:spat_g}
\end{figure}

\FloatBarrier

%In this spatially stratified example, if we let the total population of the model be $N$, the transmission matrix for the naive product would be $\textrm{Diag}(\beta_i/N)$
%% \[
%%     T = \begin{bmatrix}
%%         \frac{\beta_1}{N} & 0 & 0 \\
%%         0 & \frac{\beta_2}{N} & 0 \\
%%         0 & 0 & \frac{\beta_3}{N}
%%     \end{bmatrix}
%% \]
%While the modified product transmission matrix is given by
%% \[
%%     T = \begin{bmatrix}
%%         \frac{\beta_{11}}{N} & \frac{\beta_{12}}{N} & \frac{\beta_{13}}{N} \\
%%         \frac{\beta_{21}}{N} & \frac{\beta_{22}}{N} & \frac{\beta_{23}}{N} \\
%%         \frac{\beta_{31}}{N} & \frac{\beta_{32}}{N} & \frac{\beta_{33}}{N}
%%     \end{bmatrix}
%% \]

\section{Examples}\label{unco}

While Cartesian products allow us to construct a wide range of compartmental models, they cannot account for many useful stratified models. In this section we discuss a number of examples where the model products we have defined so far are insufficient.

\subsection{Linear Chains}\label{chains}

\subsection{Models with Testing}\label{testing}

One such example involves modeling the effects of testing for infection, inspired by the dynamics of testing during the COVID-19 pandemic \citep{gharouni2022testing}. Consider the epidemiological model in Figure \ref{fig:testify_epi} and the testing process depicted in Figure \ref{fig:testify_states}. The modified product of these two models includes a compartment for untested individuals at the hospital. However, this product is not what we want (Figure \ref{fig:testify_desired}). 
The key difference is that untested individuals entering the hospital are typically tested (i.e., moved from ``untested'' to ``awaiting results''); in our model, we will assume that they are always tested.
Therefore, the ``untested hospitalized" compartment in product model is always empty; the flow that goes to that compartment should instead be directed to the ``hospitalized/awaiting test result" compartment. Constructing the desired model would thus require an extra step to redirect this flow and remove the superfluous compartment.

\ben{can we condense/combine Figs 8-10 without making them impossible to read? (Fig 10 is pretty challenging in any case ...)}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/SIR_w_Hospital.svg}
    \caption{A simple epidemiological model that we will expand to include testing. In this model, some exposed individuals will develop asymptomatic or mild illness, in which case they stay in the community during their infectious period (and potentially transmit to others); those who instead develop severe illness will be hospitalized. (This model allows neither for within-hospital transmission nor for disease-induced mortality either inside or outside the hospital.)}
    \label{fig:testify_epi}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Testing.svg}
    \caption{A simple testing model. Individuals who test negative will, over time, revert back to the ``untested" status. This is not the case for those that test positive; at least during the early stages of the COVID-19 pandemic, someone who had tested positive for COVID-19 would assume that they were immune and would not be re-tested even if they developed COVID-like symptoms.}
    \label{fig:testify_states}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/SIR_w_Hospital_and_Testing.svg}
    \caption{The desired result of combining Figure \ref{fig:testify_epi} with Figure \ref{fig:testify_states}. Note the missing grey ``untested'' box associated with the hospital location; exposed individuals going into the hospital (enlarged, grey downward arrow starting at $E$) flow into the purple ``awaiting results'' subcompartment.}
    \label{fig:testify_desired}
\end{figure}

\FloatBarrier


\subsection{Multistrain Models and a Weak Product}\label{wp}

Many epidemics involve multiple co-circulating strains of the same pathogen \citep{gog2002dynamics, williams2021localization}. In the case of COVID-19 such variants have significant implications for the efficacy of vaccines \citep{abu2021effectiveness, koyama2020emergence} and diagnostic tests \citep{vasireddy2021review}. In more complex models, including multiple strains rapidly inflates the size of both the state space and the parameter space \citep{kryazhimskiy2007state}. One way to limit the size of these unwieldy models while continuing to include the effects of multiple strains in our model is to disallow the possibility of superinfection (i.e., an individual being infected with multiple strains at the same time). \ben{should we worry about distinguishing superinfection (infection A then B) vs coinfection (simultaneous inf with A and B)?} It would therefore be useful to define a weak product similar to the operations proposed by \cite{worden2017products} but which excludes all states corresponding to a superinfected status. One way to do this, which works well for two-strain models, is to use the standard Cartesian product but include only flows that emanate from compartments with no inflow (\term{sources}) or enter compartments with no outflow (\term{sinks}). Flows in this category typically represent initial infection (i.e., flows out of a susceptible class that is a source) or final recovery (i.e., flows into a recovered class that is a sink). In making this restriction, we exclude all the flows within the set of compartments that represent infected states, such as from an exposed to an infected compartment).
For convenience we denote this operation by $\boxplus$ and call it the \term{weak product}.

Figure~\ref{fig:ms_ns_2d} depicts a two-strain SIR model without superinfection, corresponding to the weak product of two SIR models.
Unfortunately, the weak product is not associative!
Figures \ref{fig:wp_act} and \ref{fig:wp_nonasoc} depict two different results for the weak product of three SIR models with different grouping. If the factor models for the three strains are $M_1$, $M_2$, $M_3$ then Figure \ref{fig:wp_act} depicts $(M_1\boxplus M_2)\boxplus M_3$ and Figure \ref{fig:wp_nonasoc} depicts $M_1 \boxplus (M_2 \boxplus M_3)$. Figure \ref{fig:ms_ns_3d} depicts the desired result for a three-strain SIR model with no super-infection.

\ben{I got a little bit confused by the flow here. The three-strain model, and why the weak product doesn't work in this case, is presented (are there other applications of the weak product other than multi-strain models? Seems a little limited if the \emph{only} thing it is good for is two-strain models without superinfection). Then we say ``but we can't do that with these products''. Then we talk about labeled partitions. Then we present a \emph{different} two-strain model, which I guess represents another problematic case for the weak product (but which can also be solved with labeled partitions)? The only thing we say about this model (in the last few lines of the body of the paper) is that it's problematic [and could be solved by labeled partitions?]}

It is possible to create a version of the weak product defined above that will produce the model shown in Figure \ref{fig:ms_ns_3d}. However, it requires us to distinguish between compartments that are \emph{global} sources or sinks and compartments that are sources or sinks with respect to one of the three strains specifically. That is to say, while a global sink must have no outflows, a weaker condition says that a compartment is a sink with respect to a specific pathogen if every compartment that can be reached via the outflow has the same infection status with respect to that pathogen as the original compartment. Programmatically we achieve this by introducing a concept of `labeled partitions' which separates the vertices of the model into disjoint sets corresponding to the vertices' status with respect to a specific pathogen. Each dimension of stratification in the model corresponds to a different labeled partition with each stratum corresponding to a different disjoint set. In this way we can define sources and sinks with respect to a specific set of labels rather than globally. For example, we can say a compartment $A$ is a sink with respect to a specific labeled partition if every compartment that can be reached after being in $A$ is in the same set as $A$. Figure \ref{fig:msms} outlines a compartmental model with one source compartment but two sink compartments and Figure \ref{fig:msms2} shows the weak product of two such models.
An unfortunate aspect of this construction is that several of the compartments can only be reached by individuals after they are already dead (!).
If there are relatively few such compartments a modeler may choose simply to leave them in the model and treat them all as a single compartment. But if there are many such ``zombie compartments", or if computational efficiency is a pressing concern, they could be removed  from the model.

\FloatBarrier
\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Two_strain_SIR_no_super.svg}
    \caption{A two-strain SIR model admitting no superinfection. Red Compartments indicate an infectious population whereas the population in blue compartments are not infectious}
    \label{fig:ms_ns_2d}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Three_strain_SIR_no_super_1.svg}
    \caption{A model corresponding to the product $(M_1\boxplus M_2)\boxplus M_3$. The ``SRS" and ``RSS" compartments are not sinks or sources in $M_1\boxplus M_2$; hence they have no paths to ``SRR" and ``RSR" respectively.}
    \label{fig:wp_act}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Three_strain_SIR_no_super_2.svg}
    \caption{A model corresponding to the product $M_1\boxplus (M_2\boxplus M_3)$. The ``SRS" and ``SSR" compartments are not sources or sinks in $(M_2 \boxplus M_3)$; hence they have no paths to ``RRS" and ``RSR" respectively.}
    \label{fig:wp_nonasoc}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Three_strain_SIR_no_super_actual.svg}
    \caption{A three-strain SIR model admitting no superinfection. This model cannot be constructed using only the products defined in this article.}
    \label{fig:ms_ns_3d}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/SIR_w_Death.svg}
    \caption{A single strain model with two sinks and one source.}
    \label{fig:msms}
\end{figure}

\begin{figure}
    \centering
    \includesvg[width=\textwidth]{images_redux/Two_strain_SIR_w_Death_b.svg}
    \caption{The weak product of two of the single strain models depicted in Figure \ref{fig:msms}. The grey compartments are superfluous as they correspond to changes in infection status occurring after death.}
    \label{fig:msms2}
\end{figure}

\FloatBarrier

\subsection{Levels of immunity}\label{immunity}

\subsection{STD models}\label{std}

\steve{
    These are some ideas from Ben about STD models. I've deleted the section that they originally appeared in.
    
    STD models with varying mixing rates?
    
    Activity-weighted mixing, as in STD models: the relevant value is an activity-weighted average, $\sum_j \beta_{ij} c_j I_j/\sum c_j N_j$ (or something like that)
}

\section{Conclusion}\label{conc}

Adding new strata to simple epidemiological models is closely related to taking the Cartesian product of digraphs. Modellers who want to combine sets of simple models into a single large stratified model would benefit from a toolkit based on well-defined mathematical operations. This toolkit must contain a variety of operations representing a useful subset of the numerous ways that separate strata in a model can interact. Ideally, it should also be specified in such a way that it can act \emph{automatically}, once users have clearly specified which set of rules should be used for constructing the model product. We have developed a mathematical formalism for defining such operations and used it to restate two previously proposed model operations, the naive and modified products, which represent extremes of a spectrum of interactions between strata. The naive product corresponds to the case where different strata never interact, while the modified product corresponds to scenarios where any stratum can interact with any other stratum. We generalize these previously proposed operations to a third operation that allows any level of interaction between model strata, for example to construct geographically stratified models where interactions can occur within a single location and its neighbours but not more distantly.

\ben{this needs to be restated in a more positive and useful way!}
Ultimately, all we have really done here is discuss ways to reduce the size and complexity of compartmental model products by eliminating unnecessary flows and compartments. If time and resources were unlimited we could treat every compartmental model as though it was a complete graph with every compartment connected to every other compartment. Doing this would certainly simplify the process of finding the product of two models, but simulating the result would involve spending a lot of time computing things like the probability of someone spontaneously aging thirty years and developing immunity to a pathogen they never actually encountered. Our goal however is not only to produce mathematically accurate models, but also to use those models to provide policy makers with useful and timely advice. This means that finding ways to reduce the time needed create relevant models and generate their output is a highly beneficial exercise.

Several challenges remain for anyone wishing to further develop a model construction toolkit. (1) While we have considered the structure of model products in detail, we have been less thorough in building a parallel scheme for the rules used to generalize from factor-model parameters to the parameters of the corresponding product model; a set of product types similar to the naive, modified, weak, etc. products described here would be useful. Many models (e.g., models with infection status testing) also have asymmetries in their structure that cannot be reproduced with Cartesian-like products. This suggests the need for additive operations to supplement the multiplicative operations presented here. In fact, such operations already exist in the category-theoretic approach to model operations, one reason why it could be a worthwhile project to unite the category theory and graph theory approaches. This would involve finding ``type-graphs'' that cause the category theoretic operations known as ``pull-backs'' and ``push-outs'' to  reproduce the results of graph theory operations \citep{fong2018seven, Libkind2022an, libkind2021operadic, baez2023categorical, baez2017compositional}.

We are heavily motivated by the desire to develop software to facilitate model construction. One insight of our investigations is the utility of a system of so-called ``labeled partitions'', which divide the compartments of a model into mutually exclusive groups. Each group in such a division will contain all compartments that are in the same level of some dimension of stratification and the groups can be labeled accordingly. By applying several such divisions to a model, one for each dimension of stratification, it becomes possible to specify important subsets of the model compartments. Using this system of labels and partitions provides an easy way to address issues like the non-commutativity of the weak product and the presence of ``zombie compartments'' discussed in Section \ref{wp}.

Although theoretical and practical challenges with the application of binary operations on model space remain, our approach forms the basis of a powerful toolkit for the construction of complex, stratified, compartmental models.

\ben{beyond further streamlining and clarification, the two things I might like to see in this paper (if they fit in) are (1) something about the decomposition of transmission matrices into (susceptibility $\times$ contact $\times$ infectivity, as in the \href{https://canmod.github.io/macpan2/articles/state_dependent_rates.html}{state-dependent rates vignette} \ldots and (2) more on labeled partitions, which we introduce at the beginning of the paper and only get back to at the very end, a little bit. Can we add stuff on labeled partitions without driving ourselves crazy or making this take much longer \ldots ?}