Structure and Hierarchy of Influenza Virus Models Revealed by Reaction Network Analysis

Abstract

Influenza A virus is recognized today as one of the most challenging viruses that threatens both human and animal health worldwide. Understanding the control mechanisms of influenza infection and dynamics is crucial and could result in effective future treatment strategies. Many kinetic models based on differential equations have been developed in recent decades to capture viral dynamics within a host. These models differ in their complexity in terms of number of species elements and number of reactions. Here, we present a new approach to understanding the overall structure of twelve influenza A virus infection models and their relationship to each other. To this end, we apply chemical organization theory to obtain a hierarchical decomposition of the models into chemical organizations. The decomposition is based on the model structure (reaction rules) but is independent of kinetic details such as rate constants. We found different types of model structures ranging from two to eight organizations. Furthermore, the model’s organizations imply a partial order among models entailing a hierarchy of model, revealing a high model diversity with respect to their long-term behavior. Our methods and results can be helpful in model development and model integration, also beyond the influenza area.

1. Introduction

Influenza is an infectious respiratory disease, annually infecting 5–15% of the human population and causing epidemics that result in 3–5 million severe cases with 300,000–500,000 deaths each year [1]. The annual recurrence of epidemics is caused by the continuous alteration of seasonal influenza viruses, which enables them to efficiently escape the immune system even due to previous infections or vaccinations [2]. The major burden of disease in humans is caused by seasonal influenza A (IAV) and influenza B viruses, causing symptoms varying from mild respiratory disease characterized by fever, sore throat, headache and muscle pain to severe and in some cases lethal pneumonia and secondary bacterial infections [3].

The long-term spread of influenza viruses in the human population and the acute nature of influenza virus epidemics is driven by the global movement of these viruses. Differences in seasonal epidemics caused by influenza viruses are mainly driven by differences in the rates of virus evolution. The single-stranded RNA segments of influenza viruses, which are located inside the virus particle (or virion), evolve rapidly and thus can escape the host’s immune response very efficiently.

Several ordinary differential equations (ODEs) models have been developed to provide insight into within-host dynamics of influenza A virus infections (for reviews, see [4,5,6,7,8]). These models work in a time scale of days and describe the concentration dynamics of target cells, immune system components, viral load, and sometimes co-infecting pathogens. The models differ in terms of complexity and state space dimensions, which are between 3 and 15 for the models examined here. While the low-dimensional models can be analyzed completely and in a straightforward way (e.g., by calculating their fixed points and stability analysis), the characterization of the entire behavioral spectrum of complex models is more difficult (see for example [9]).

Here, we present an approach to understand the overall structure of these models that allows them to be related to each other in a simple way. To this end, we apply chemical organization theory [10,11] to obtain a hierarchical decomposition of each model into chemical organizations. A chemical organization is a sub-set of species (i.e., dimensions or model components, like, uninfected cells or viruses) that cannot generate any other species (property of closure) and that can self-maintain its own species, i.e., any species consumed by a process within the organization can be regenerated by a process within the organization. The organizations of an ODE model are rigorously related to its long-term dynamics in the following way: Given a stationary state of the ODE model, the set of species with strictly positive concentrations must be an organization [10]. The same is true for all practically relevant periodic and chaotic attractors [12]. Note that the advantage of this approach is that decomposition into organizations is based solely on the model structure (i.e., reaction rules) and thus is independent of kinetic details, like rate constants. The relation between measured data, ODE model, and organizations is depicted in Figure 1.

By applying the method to twelve models of influenza A virus infection, we found different types of model structures ranging from two to eight organizations. Furthermore, the models’ organizations imply a partial order among models. The resulting hierarchy of models can help to select a suitable model for certain data or serve as a framework for further model development.

We provide reaction network files for all models and a software tool for computing their organizations (https://github.com/stephanpeter/orgsflu).

2. Materials and Methods: Procedure for the Organizational Analysis

To illustrate our method, we follow a basic ODE model of influenza dynamics, namely the target cell limited model by Baccam et al. [13], called Baccam Model in the following. Note that we refer throughout this paper to a model by the first author’s name of the respective publication.

The Baccam Model is based on in vivo experimental data and contains three variables: the number of susceptible and uninfected target (epithelial) cells T, the number of infected cells I, and the number of infectious-viral titer V. The dynamical behavior of the variables is given by the ODEs shown in Figure 2a. That is, target cells become infected and thus converted to infected cells at a rate $\beta TV$, infected cells die spontaneously at rate $\delta I$, virus proliferates at a rate $pI$ and dies at a rate $cV$. The parameters $\beta$, $\delta$, p and c are, as usual, positive real numbers (cf. [13] for actual values).

Throughout this paper, the following coloring scheme for particular classes of species is used to improve readability:

• Uninfected (target) cells or those resistant/refractory to infection are marked in blue, e.g., T.
• Infected cells, partially or latently infected cells, and viruses are marked in magenta, e.g., I and V.
• Species belonging to the active immune response are marked in green. It is worth noting that the first two models analyzed in this paper [13,14] do not explicitly have immune system components.
• Bacterial co-infection species are marked in orange. These species are only occurring in Smith’s model [15].
• Text referring to any other species is marked in black, e.g., transient target cell states, passive immune system, or dead cells.

For simplicity, the models’ variable names also denote the related species. For example, V denominates not only the number of viruses in the ODE model (Figure 2a), but also the virus itself (e.g., Figure 2b).

2.1. Deriving the Reaction Network from the ODE System

In a first step, we need to obtain the reaction network underlying the ODE model. A reaction represents, for example, a cell infection by a virus, the generation of new viruses from an infected cell or the spontaneous death of a cell. The reaction rules can be derived from the ODEs in a straightforward way [16]. This step can also be performed by an online tool presented by Soliman and colleagues [16]. Note that in modeling one first creates the network and then derives the ODEs. For our analysis, we have to take the other direction.

For this purpose, we have to investigate the kinetic terms (kinetic laws) of the ODE (Figure 2a):

• The term $\beta TV$ represents the reaction $R_1:T+V\to I+V$, which in turn denotes the transformation of an uninfected target cell T to an infected cell I catalysed by the virus V.
• The terms $-\delta I$ and $-cV$ represent reactions $R_2:I\to \varnothing$ and $R_4:V\to \varnothing$ which are the outflow of infected cells I resp. virus V.
• The term $pI$ represents the reaction $R_3:I\to I+V$ which is the production of viruses V catalysed by infected cells I.

The set of species $\mathcal{S}=\{T,I,V\}$ together with the set of reactions $\mathcal{R}=\{R_1,R_2,R_3,R_4\}$ constitute the so-called reaction network $(\mathcal{S},\mathcal{R})$ associated with the model. The set of reactions together with their kinetic parameters are depicted in Figure 2c. Note that for clarity we use different types of underlining to highlight certain recurring kinetic terms in the ODE:

• Single underline for the transformation of uninfected cells into infected ones by the action of viruses.
• $\underset{=}{\text{Double underline}}$ for kinetic terms involving interferon.

We call the species on the left-hand side (LHS) of a reaction R support of R and write $supp\left(R\right)$, e.g., $supp\left(R_1\right)=\{T,V\}$ (see Figure 2d). Analogously, we call the set of species occurring on the right-hand side (RHS) of a reaction R products of R and denominate this set by $prod\left(R\right)$, e.g., $prod\left(R_1\right)=\{I,V\}$.

Furthermore, let us recall the stoichiometric matrix $\mathbf{N}=\left(n_{ij}\right)$ of a reaction network [17]. The element $n_{ij}$ in the i-th line and j-th column of $\mathbf{N}$ denotes the net-production of the i-th species in reaction $R_j$. The net-production $n_{ij}$ is the difference between the number of occurrences (stoichiometric coefficient) of species i in the RHS of reaction $R_j$ minus the number of occurrences of species i in the LHS of reaction $R_j$. For example, $n_{21}=1-0=1$ because the second species (I) appears in the first $R_1$ once as a reactant in the support of $R_1$ (LHS) but does not appear in $R_1$ as a product (RHS). For our example in Figure 2, the stoichiometric matrix becomes:

$$\begin{array}{c}\hfill \mathbf{N}=\left(\begin{array}{cccc}-1& 0& 0& 0\\ 1& -1& 0& 0\\ 0& 0& 1& -1\end{array}\right).\end{array}$$

(1)

2.2. Computing the Organizations from the Reaction Network

From the reaction network, we can compute the (chemical) organizations of the model. Each organization is a subset of species that is closed and self-maintaining [10,18]. In the following, let $S\subseteq \mathcal{S}$ be a subset of species and n be the total number of reactions of the reaction network ($n=4$ in our example).

We call S closed if and only if all reactions $R\in \mathcal{R}$ with $supp\left(R\right)\subseteq S$ fulfill $prod\left(R\right)\subseteq S$ too [10,18]. This means that the products of a reaction R with support in S are also in S. In other words, no species outside of S can be produced by the reactions “running on” S. As an example, we assume $S=\{T,I\}$. The reactions with support in S are $R_2$ and $R_3$. However, $R_3$ produces species V, which is not in S. Thus, S is not closed.

We call a vector $\mathbf{v}\in {\mathbb{R}}^n$ flux vector for S if and only if it fulfills

$$\begin{array}{c}\hfill {\mathbf{v}}_j\left\{\begin{array}{cc}>0,\hfill & \mathrm{if}supp\left(R_j\right)\subseteq S,\hfill \\ =0,\hfill & \mathrm{if}supp\left(R_j\right)⊈S.\hfill \end{array}\right.\end{array}$$

(2)

Thus, all flux vectors for S have in common that those components are strictly positive which correspond to reactions that can run on S, while all other entries are zero. As an example, consider $S_{Example}=\{T,I\}$ again. We know that the reactions $R_2$ and $R_3$ can "run on" it, i.e., they have support in $S_{Example}$. Thus, ${\mathbf{v}}_1={(0,1,3,0)}^T$ or ${\mathbf{v}}_2={(0,5,2,0)}^T$ are example flux vectors for $S_{Example}$.

We call S self-maintaining if and only if there exists (at least one) flux vector $\mathbf{v}\in {\mathbb{R}}^n$ for S that fulfills

$$\begin{array}{c}\hfill \mathbf{N}·\mathbf{v}\ge \mathbf{0},\end{array}$$

(3)

i.e., ${(\mathbf{N}·\mathbf{v})}_i\ge 0$ for all $i=1,\dots ,n$, where n is again the total number of reactions [10,19,20,21]. Roughly speaking, if S is self-maintaining, it has the potential to sustain the amount of its species above a certain level. Our example set $S_{Example}$ is not self-maintaining because ${(\mathbf{N}·\mathbf{v})}_2<0$ for all flux vectors for $S_{Example}$. That is, species number 2 (the infected cells I) can not be maintained by this set.

As mentioned in the beginning of this section, we call S an organization if and only if it is both closed and self-maintaining. Clearly $S_{Example}$ is not an organization, as it has neither of these properties. In Figure 2d, the so-called Hasse diagram of organizations of this model can be seen. In it, two organizations are linked by a line if the lower one is a subset of the upper one and there is no organization in between. The Hasse diagram for the Baccam Model contains only the two organizations $O_1^B=\varnothing$ and $O_2^B=\left\{T\right\}$ (see Figure 2d). The superscript “B” stands for Baccam Model and with the subscripts we refer to the organizations within a model. The organization $O_2^B$ represents an organism without any influenza A virus infection. Note that there is no organization with all species, i.e., representing the infected body.

2.3. The Role Organizations Play in the Dynamics

Given a fixed point x of an ODE describing a reaction system, the set of species with strictly positive concentrations in x is an organization. This is shown in [10]. This in turn means that, if a subset of species is not an organization, then the system does not have a fixed point with exactly the chosen subset of species. This is not only true for fixed points but also for other attractors [12]. Attractors are those states that a system approaches in the long-run and once reached never leaves anymore. Besides converging towards a fixed point, the long-term behavior can be also periodic oscillations or chaotic trajectories. In particular, it was shown that the long-term behavior of the system always tends at least to one organization [12]. Thus, organizations rule the long-term behavior of such dynamical reaction system. Note that the case of a system tending towards a fixed point is included in this statement as a special case. For example, the simulation results of the Baccam model (Figure 2a and Ref. [13]) suggest that after about two to three days the species begin to decay to finally arrive in an organization namely the empty set.

3. Results and Discussion

In the literature, there exist several mathematical models of IAV dynamics that are derived from experimental data, reviewed in Refs. [4,5,6,7]). These models differ in their complexity, e.g., the number of reactions and the number of species, depending on the available experimental data used for parameter fitting and questions to be answered. For example, models can include eclipse phases, an innate immune response, or an adaptive immune response.

After having exemplified our method in Section 2 by an analysis of the basic target cell limited model by Baccam et al. [13], we present now the full analysis for eleven additional more detailed influenza models of IAV dynamics, with up to 15 variables (species) and 45 reactions (cf.

Table 1. Overview of all models and organization types contained. An organization type like $\mathbf{X}\mathbf{X}$ denotes the type of species contained in an organization, according to our coloring scheme. The set of organization types of a model is called its signature.

Model	Number of Variables	Number of Reactions	Number of Organizations	Organizations & Signature
Baccam [13] 2006	3	4	2	$O_1^{Ba1}=\varnothing$
				$O_2^{Ba1}$	X
Miao [14] 2010	3	5	3	$O_1^M=\varnothing$
				$O_2^M$	X
				$O_3^M$	X	X
Baccam II [13] 2006	4	5	2	$O_1^{Ba2}=\varnothing$
				$O_2^{Ba2}$	X
Pawelek [23] 2012	5	9	2	$O_1^P=\varnothing$
				$O_2^P$	X
Smith [15] 2016	5	12	4	$O_1^{Sm}=\varnothing$
				$O_2^{Sm}$	X
				$O_3^{Sm}$	X
				$O_4^{Sm}$	X	X
Handel [24] 2010	7	12	5	$O_1^{Ha}=\varnothing$
				$O_2^{Ha}$	X
				$O_3^{Ha}$	X
				$O_4^{Ha}$	X	X
				$O_5^{Ha}=\left\{all\right\}$	X	X	X
Hernandez [26] 2012	7	16	2	$O_1^{He}$	X	X
				$O_2^{He}=\left\{all\right\}$	X	X	X
Cao [27] 2015	7	26	3	$O_1^C=\varnothing$
				$O_2^C$	X
				$O_3^C=\left\{all\right\}$	X	X	X
Saenz [28] 2010	8	12	4	$O_1^{Sa}=\varnothing$
				$O_2^{Sa}$, $O_3^{Sa}$, $O_4^{Sa}$	X
Hancioglu [29] 2007	10	44	3	$O_1^{Hcg}$	X
				$O_2^{Hcg}$	X	X
				$O_3^{Hcg}=\left\{all\right\}$	X	X	X
Bocharov [30] 1994	10	45	2	$O_1^{Bo}$	X
				$O_2^{Bo}$	X	X
Lee [9] 2009	15	37	8	$O_1^L$	X
				$O_2^L$, $O_3^L$, $O_4^L$, $O_5^L$, $O_6^L$	X	X
				$O_7^L$, $O_8^L=\left\{all\right\}$	X	X	X

for an overview at the end). Note that for our analysis we abstract from kinetic details, that is, the organizations are independent of particular settings of parameter values.

3.1. Target Cell Limited Model by Miao et al. (Miao Model, M, 2010)

The models by Miao et al. [14] are designed to fit experimental in vivo data from mice [6,7]. The first one ([14], Equation (1)) depends on measured time-series. The second one ([14], Equation (2)) is a simplified version of the first one, neglecting the terms depending on those time-series and still leading to a good fit within the first five days after infection [14]. Thus, we analyze this second model (Miao Model).

Compared to the basic Baccam Model, the Miao Model (Figure 3a) has the same three variables (named differently) and one additional reaction, $E_P\to 2E_P$. This reaction represents the self-replication of target cells $E_P$ taking place at a rate ${\rho }_E{E}_P$. The full set of reactions can be found in the Appendix A (Figure A2).

In the Hasse diagram of organizations (Figure 3b), a new “full” organization $O_3^M$ appears, which contains all three species. Thus, organization $O_3^M$ reflects the slight difference between the two models: in the Baccam Model, uninfected target cells T are only the susceptible ones and can not increase in number, but in the Miao Model uninfected cells $E_P$ are reproduced repeatedly by the organism. Thus, in the Baccam Model, infection is limited inherently by the limited number of uninfected target cells, while in the Miao Model the limitation of an infection in time and number of infected cells and viruses depends on other mechanisms:

3.2. Target Cell Limited Model with Delayed Virus Production (Baccam II Model, Ba2, 2006)

The Baccam II Model [13,22] contains one more species than the Baccam Model presented in the methods section above. That is, there are now two types of infected cells: those which do not yet produce viruses $I_1$ and those which actively produce viruses $I_2$. In addition, there is only one new reaction, which transforms infected cells of type $I_1$ into type $I_2$ at rate $k{I}_1$ (Figure 4a). However, the Hasse diagram of organizations remains the same when compared with the basic Baccam Model [13].

3.3. Innate and (Simple) Adaptive Immune Response (Pawelek Model, P, 2012)

The Pawelek Model [23] contains 11 parameters and was designed to fit in vivo experimental data of horses [6,7]. The model has five variables and nine reactions. Like the basic Baccam Model, it contains uninfected target cells (T), infected cells (I), and viruses (V). Furthermore, there are two new species: interferon (F) and uninfected cells that are refractory to infections (R) because of the antiviral effect induced by interferons.

Investigating the reaction network (Figure A4 in the Appendix A) derived from the differential equations (Figure 5a), we can see that, like in the basic Baccam Model, self-replication of uninfected cells T is missing. However, due to the two new species R and F, we have five new reactions, which are neither included in the Baccam Model nor in the Miao Model. One of these five reactions is the spontaneous decay of interferon F at a rate $dF$. The other four new reactions describe interactions between different species:

• The rate term $\varphi FT$ represents the transformation of uninfected target cells to refractory cells catalysed by interferon.
• The reverse shift back from refractory to simple uninfected cells is represented by the term $\rho R$.
• Furthermore, infected cells are deleted by the action of interferon at a rate $\kappa IF$.
• Interferon is produced in the presence of infected cells at a rate $qI$.

Even though we have more species and more reactions, we get the same small pattern of organizations as in the basic Baccam Model (Figure 5b). Both models have in common that there is no self-replication of target cells. This might be one reason for the missing of other and bigger organizations which could contain species related to infection and/or immune response. This in turn means that, like the Baccam Model, this model implicitly treats infections and immune responses as phenomena that can only appear in a limited (transient) time span. The Hasse diagram of organizations (Figure 5b) tells us that the system necessarily tends towards a state of healthiness, which is represented by the organizations $O_1^P=\left\{\right\}$ and $O_2^P=\left\{T\right\}$, showing absence of infection and immune response.

3.4. A Model Including Bacterial Co-Infection (Smith Model, Sm, 2016)

The Smith Model [15] contains 15 parameters and compared to experimental in vivo data from mice. It has five variables and 12 reactions. Like in the previous models, we have susceptible target cells (T) and viruses (V). Note that T is only consumed in this model but not produced. Contrarily to previous models, we have two kinds of infected cells ($I_1$ and $I_2$) here as well as bacteria (P). Bacteria P represent bacterial co-infection during or after virus infection. Infection is modeled as a transformation of susceptible target cells T into infected cells $I_1$ catalyzed by viruses V (see underlined terms in Figure 6a). Infected cells $I_1$ in turn spontaneously transform into $I_2$ at rate $k{I}_1$. Only infected cells $I_2$ produce viruses V at a rate $p{I}_2$. Furthermore, infected cells $I_2$ produce viruses V together with bacteria P. Bacteria P are self-replicating (rate term $rP$). Viruses V are the only species influencing bacteria P.

Figure 6b shows the Hasse diagram of organizations. The smallest one is the empty set. The biggest one is $O_4^{Sm}$, which contains susceptible target cells T and bacteria P. It represents an organism without viral but with bacterial infection. Between those two extreme organizations, we find $O_2^{Sm}=\left\{\begin{array}{c}\hfill \mathbf{T}\end{array}\right\}$ and $O_3^{Sm}=\left\{\begin{array}{c}\hfill \mathbf{P}\end{array}\right\}$. Thus, $O_2^{Sm}$ represents the healthy organism without any infection. $O_3^{Sm}=\left\{\begin{array}{c}\hfill \mathbf{P}\end{array}\right\}$ could mark the state after a viral-bacterial co-infection: after viral infection and because of the death of all target cells as well as all viruses only bacteria remain.

3.5. Innate and Adaptive Immune Response (Handel Model, Ha, 2009)

The Handel Model [24] contains eight parameters and was designed to fit experimental in vivo data from mice [6,7,25]. It has seven variables (see Figure 7) and only 12 reactions (see Figure 7a):

• Infection is catalyzed by viruses V and transforms uninfected cells U to latently infected cells E and viruses V are consumed thereby. Latently infected cells E transform into infected cells I autonomously, which in turn transform into dead cells D autonomously too. Finally, the transformation of dead cells D into non-infected cells U closes the circle.
• The remaining three species V, F and X form an almost totally separate subsystem since the only interaction with the four species from the "circle" mentioned above is the catalysis of the infection by viruses V.
• The interactions within the subsystem $\{V,F,X\}$ consisting of viruses V and immune responses F and X are as follows:

  –

  Viruses V catalyze the proliferation of F and X. In the Hernandez model, proliferation of interferon F is catalyzed by infected cells instead of viruses.

  –

  There is no direct interaction between innate immune response F and adaptive immune response X.

  –

  The adaptive immune response X deletes viruses directly. Innate immune response F inhibits the self-replication of the viruses which is represented by the denominator of the fraction ${\displaystyle \frac{p\mathbf{I}}{1+\kappa \mathbf{F}}}$. We ignore the inhibition because whether the rate is zero or not is independent of F.

The Hasse diagram Figure 7b shows five organizations. For the first time, it contains the empty set as well as the set of all species at the same time. Between these extremes, we find $O_2^{Ha}=\left\{\mathbf{U}\right\}$ representing the healthy organism. The Baccam, Miao, and Pawelek models exhibit the same organization. Structurally, the Hasse diagram of the Handel Model is very similar to that of the Smith Model (Figure 6b). The first reveals an autonomy of the adaptive immune response X, whereas the latter does this same for bacteria P.

Temporal Dynamics

For the Handel Model, we perform dynamical simulation in order to show how the organizational hierarchy helps also to understand transient short-term behavior. We start at $t=-20d$ with an uninfected state, i.e., $7\times {10}^9$ uninfected cells. After 20 days, at $t=0$, we add ${10}^4$ virus particles. The resulting seven-dimensional trajectory in state space is shown in Figure 8. Projecting this trajectory to organizations results in a more abstract view of the dynamics, shown as a dashed curved arrow in Figure 7b. The system starts in organization $O_2^{Ha}$ (uninfected organization), moves after adding virus particles at $t=0$ into organization $O_5^{Ha}$ (infected organization with immune response), and drops after 37 days into organization $O_4^{Ha}$ (immune response active, virus absent).

The projection of a state x to an organization O follows the procedure suggested by Dittrich and Speroni d.F. [10]: First, we generate a set S of those species whose concentration is greater than a particular threshold (here: ${10}^0=1$). Then, we generate the closure of this set by adding all species that can be produced from the set. Finally, we take the largest organization O that is a subset of that closure. For example: At $t=0$ by adding viruses to the system, we have $S=\{U,V\}$, whose closure is the set of all species, which is also an organization, here; thus, the state at $t=0$ is projected to organization $O_5^{Ha}$. At $t=60d$, we have $S=\{U,X,D\}$, whose closure is again $\{U,X,D\}$ and the largest organization contained is $O_4^{Ha}=\{U,X\}$. Thus, as can be seen in Figure 8, the system state is projected to organization $O_4^{Ha}$, in which it stays for $t>37d$.

3.6. Innate Immune Response and Resistance to Infection (Hernandez Model, He, 2012)

The 13 parameters of the Hernandez Model [26] were fitted to data from many different sources. The model contains seven variables and 16 reactions (see Figure 9). The species refer to viruses (V), interferon (F) and natural killer cells (K) as well as four types of respiratory tract epithelial cells: healthy/uninfected ($U_H$), partially infected ($U_E$), infected ($U_I$) and resistant to infection ($U_R$). Compared to the Pawelek Model, there are two qualitatively new species: partially infected cells $U_E$ and killers K.

Next, we state some remarks about the reactions:

• There is an infection reaction catalyzed by viruses like in all previous models but with one difference: during infection, healthy cells $U_H$ first transform to partially infected cells $U_E$ and only after that they transform spontaneously to infected cells $U_I$ at a rate $k_e{U}_E$.
• Interferon catalyzes the transformation of healthy cells to resistant cells $U_R$, like in the Pawelek Model. However, in the Pawelek Model, interferon removes infected cells. Here, interferon’s production is catalyzed by infected cells $U_I$ at a rate $a_I{U}_I$. There is no further influence of interferon on any other species.
• Infected cells are removed by natural killers K, which also delete partially infected cells in this model. The production of killers K is catalyzed by infected cells $U_I$ at a rate ${\mathsf{\Phi }}_K{U}_I$.
• Note that here we have an constant inflow of healthy cells $U_H$ at a rate $S_H$ (first differential equation). Thus, healthy cells cannot converge to zero.

The Hasse diagram of organizations consists only of two organizations (Figure 9b). For the first time, the empty set is not an organization because the empty set is not closed due to a constant inflow of healthy cells $U_H$ and killers K, represented by the reaction $\varnothing \to U_H$ and $\varnothing \to K$, respectively. The smallest organization $O_1^{He}=\{U_H,K\}$ can be regarded as a state of healthiness. Contrarily, the second organization $O_2^{He}$ contains all species (as in the Miao Model) and therefore can be interpreted as the infected organism exhibiting immune response to infection.

3.7. A Model with More Detailed Immune Response (Cao Model, C, 2015)

The Cao Model [27] consists of 20 parameters and has been derived by referring to experimental data from ferrets. The model has seven variables and 26 reactions. Like in most of the previous models, we have (uninfected) target cells (T), infected cells (I), viruses (V), resistant cells (R), and interferon (F). Furthermore, there are B cells B and antibodies A.

According to the ODE shown in Figure 10a, B cells are only influenced by viruses: viruses support the production of B cells (rate term $m_{1}V$), but the more B cells that are present, the more of them are destroyed by viruses again (term: $m_{1}VB$). B cells influence only one other species, namely antibodies B, which they produce. Antibodies A influence only one other species, namely viruses, which are destroyed by this reaction at rate $\mu VA$. Antibodies in turn are influenced by B cells positively and by viruses negatively.

There are three organizations in this model (Figure 10b): the empty set, the healthy organism without any infection and without any immune response ($O_2^C=\left\{\mathbf{T}\right\}$) and the organization $O_3^C$ containing all species.

3.8. Innate Immune Response and Eclipse Phase (Saenz Model, Sa, 2010)

The Saenz Model [28] requires 12 parameters and was designed to fit experimental in vivo data from horses [6,7]. The model contains eight variables and 12 reactions (Figure 11a). There are no adaptive immune response, no dead cells, and no natural killer cells. However, the model contains viruses V and interferon F. There is an eclipse phase ($E_1$ and $E_2$) here as well as prerefractory and refractory cells. In particular, epithelial cells are represented by six species: susceptible (T), eclipse phases ($E_1$ and $E_2$), infectious (I), prerefractory (W), and refractory (R). Thus, the new features are the inclusion of two eclipse phases and three steps for the transformation of uninfected cells to refractory cells.

The Hasse diagram of organizations is composed by four organizations: $O_1^{Sa}=\left\{\right\}$: the empty set; $O_2^{Sa}=\left\{\mathbf{T}\right\}$: representing a healthy organism; $O_3^{Sa}=\left\{\mathbf{R}\right\}$: there is no consuming reaction for refractory cells R; $O_4^{Sa}=\{\mathbf{T},\mathbf{R}\}$: also representing a healthy organism that contains refractory cells maybe as the remains of a previous infection.

The Hasse diagram is very similar to that from of the Handel Model. There are only two differences:

• The “full” organization is missing here. For sure, one of the reasons is that there is no reaction producing susceptible cells T. Thus, when viruses V or interferon F are present susceptible cells T can not survive and the “full” organization neither.
• Adaptive immune response is replaced by refractory cells in the organizations here.

3.9. Focusing on Innate and Adaptive Immunity (Hancioglu Model, Hcg, 2007)

The Hancioglu Model [29] contains 28 parameters, 10 variables, and 44 reactions. It has not been mathematically fitted to data but has been designed to meet specific general criteria [6,7]. The ODEs (Figure 12a) describe the dynamics of the following 10 species: viruses (V), healthy cells (H), infected cells (I), interferon (F) and resistant cells (R). The remaining species are new: antigen presenting cells (M), effector cells (E), plasma cells (P) antibodies (A) and antigenic distance (S). There are no species for an eclipse phase in this model.

Looking at the reaction network (Figure A8), we can see again a reaction for infection, i.e., the transformation of healthy cells H into infected cells I catalyzed by viruses V at a rate ${\gamma }_{HV}VH$ (single underlined in Figure 12a). Furthermore, interferon F is produced catalytically by antigen presenting M and infected cells I, decays spontaneously at a rate $a_{F}F$, and is additionally removed when converting healthy cells H into resistant cells R by the reaction $H+F\to R$ at rate $b_{HF}FH$ (double underlined in Figure 12a).

The Hancioglu Model has three organizations (Figure 12b):

• The smallest one is $O_1^{Hcg}$, which contains all the species responsible for the immune response.
• $O_1^{Hcg}$ is a subset of $O_2^{Hcg}$, which additionally contains species H and R, representing the healthy organism without infection but with the immune response turned on.
• $O_3^{Hcg}$ is the “full” organization containing all the species of the models and thus representing the organism with infection and immune response.

Thus, all the organizations represent meaningful states of the organism. However, there is no organization that only consists of healthy cells without any infection and immune response. Note that almost all the previous models except for Hernandez have such an organization.

3.10. Model with Delay Differential Equations (Bocharov Model, Bo, 1994)

The Bocharov Model [30] contains 49 parameters and was designed to fit experimental in vivo data from humans [6,7]. It includes 10 variables and 51 reactions (Figure A9). Only here and in the Lee Model (below) we have differential equations with delay, i.e., some rates depend on variable values from the past (Figure 13a). Because the delay does not matter in a steady-state, we can also neglect the delay when analyzing the chemical organizations of a delay differential equation model.

Note that this is by far the oldest model analyzed here. The names of the variables are a bit particular when compared to those of the previously analyzed models. As in all the other models, we have viruses $V_f$ and infected cells C. Furthermore, there are destroyed epithelial cells m as in the Handel Model. All other species belong to the immune response. Note that only in this model there is no state variable for uninfected, healthy cells. Bocharov et al. represent these healthy cells implicitly by subtracting the amounts of infected-cells C and destroyed epithelial cells m from the initial total amount of target epithelial cells $C^{*}$. Since all the other models analyzed here have a related variable, we inserted the variable $U=C^{*}-C-m$ for uninfected cells together with its ODE to make the model comparable to the others.

Due to the fact that the majority of the species belongs to the immune response, this is the case for most of the reactions too. These species of the immune response form exactly the organization $O_1^{Bo}$, only macrophages $M_V$ are missing.

Similarly to the Hancioglus model, the smallest organization $O_1^{Bo}$ is an organization with immune response but without infection (C, $V_f$). There is only one further organization that contains only one more species than $O_1^{Bo}$, namely uninfected cells U. This organization we already found in three of the previous models. However, for the first time, there is no bigger organization in this model. Thus, virus infection is necessarily transient.

3.11. Complex Dual-Compartment Model (Lee Model, L, 2009)

The Lee Model [9] is the most complex model considered here. It contains 48 parameters and was designed with respect to experimental in vivo data from mice [6,7]. It has 15 variables and 37 reactions (Figure A10). Like in the Bocharov Model, Lee et al. apply delay differential equations.

Note that this model is the only one analyzed here that distinguishes between lung compartment and lymphatic compartment. There is one species representing uninfected (healthy) cells ${\mathbf{E}}_{\mathbf{p}}$ and three species for modelling infection: ${\mathbf{E}}_{\mathbf{P}}^{*}$, ${\mathbf{D}}^{*}$ and viruses $\mathbf{V}$. The remaining species belong to the immune response, colored black when naive to infection, while colored green when activated for infection.

Note that we write a species in the organizations in Figure 14b in bold text, if it is “new”, that is, not contained in neither of its subset organizations. The Hasse diagram contains eight organizations. The smallest one is $O_1^L=\{{\mathbf{E}}_{\mathbf{P}},\mathbf{D},H_N,T_N,B_N\}$ and contains exactly the uninfected cells as well as the naive part of the immune response. The biggest organization contains all species. Between these two “extreme” organizations are six further organizations containing different parts of the activated part of the immune response.

3.12. Hierarchy of Influenza A Virus Models

In order to construct a hierarchical map of all investigated models, we characterize a model by a signature of organizations, which is a set of organization types. For example, the signature of the Handel Model (Figure 7b) is the set $\{\varnothing ,\mathbf{X},\mathbf{X},\mathbf{X}\mathbf{X},\mathbf{X}\mathbf{X}\mathbf{X}\}$. An organization type like $\mathbf{X}\mathbf{X}$ means that there is at least one organization that contains uninfected (target) cells ($\mathbf{X}$) and species of the active immune response $\left(\mathbf{X}\right)$. The deviation of the signatures for all models is shown in Table 1. Note that we ignore species colored black. We include the empty set ∅ because this distinguishes models without any inflow from those that possess an inflow of some species.

Now, we can obtain a partial order among models by defining that a model A is smaller or equal to another model B (A ≤ B), if the signature of A is a subset of the signature of model B. For example, the Hernandez Model is smaller than the Lee Model because $\{\mathbf{X}\mathbf{X},\mathbf{X}\mathbf{X}\mathbf{X}\}\subseteq \{\mathbf{X},\mathbf{X}\mathbf{X},\mathbf{X}\mathbf{X}\mathbf{X}\}$. This partial order among models leads to a hierarchical map of models, which is visualized by a Hasse-diagram in Figure 15. Note that a model A that is smaller than a model B according to this partial order can possess more species and reactions than B.

In Figure 15, we can see that all models have organizations with uninfected, healthy cells ($\mathbf{X}$). There are models that furthermore have infection ($\mathbf{X}$) and/or immune response ($\mathbf{X}$) in their organizations. There are exactly two models (Hancioglu and Hernandez Model) with immune response in all their organizations which means that these models implicitly assume immune response to be active all the time. Among the models neglecting immune response are those which have infection (Miao Model) or bacteria ($\mathbf{X}$) (Smith Model) in their organizations and also those that do not (Baccam and Baccam II Model). For models involving immune response, the situation is more complex. There are those that only have healthy cells in their organizations (Pawelek and Saenz Models). This means that these models implicitly exclude infection and immune response from the long run and thus treat them as transient phenomena a priori. The Bocharov Model is the only one that exhibits only healthy cells and immune response in its organizations but no infection. The remaining five models include all kinds of species (except for bacteria of course) in their organizations.

By looking at the hierarchy of models, it becomes evident that there is space for more models. Above the Smith and Handel Model, there could be one in which virus infection as well as bacterial coinfection can be simultaneously persistent (“fully persistent models” denote such hypothetical models in Figure 15). Another extreme case would be a “fully-transient model” in which we have only transient dynamics and all species would finally tend to zero. Such a model would be the smallest one in our partial order of models (Figure 15).

The derived hierarchical map of models might be used to choose the most appropriate model for a particular domain and data set: The model should contain at least one organization for each set of species that were experimentally observed to survive in the long run. If there are several models with such organizations, the one with the smallest organizations might be chosen to provide maximum efficiency in modeling. Table 1 as well as Figure 15 might guide the selection process, complementing established quantitative selection methods, such as those using the area under the viral load curve [31].

4. Conclusions

By analyzing twelve published IAV immune system models, we have shown that we can compute, independently of quantitative kinetic data like rate constants or kinetic laws, all chemical organizations of a typical IAV model, which provides a hierarchical decomposition of the model and an overview of its potential long-term behavior.

It turned out that the derived organizations are meaningful with respect to the model’s domain. That is, the composition of species inside an organization can be related to a particular state of the organism, like “healthy”, “infected”, or “virus controlled by active immune response”, and it is possible to annotate organizations accordingly.

By deriving an organizational signature from a model’s organizations, we obtained a novel classification scheme and a hierarchy of models with respect to their qualitative long-term behavior. Although the classification via organizational signature is quite coarse-grained, the analysis revealed still a high diversity of models. That is, the models have different potentiality with respect to which variables persist in the long turn and which vanish. Furthermore, the hierarchy map of models contains various empty territories, suggesting space for potential future models.

We envision as a practical use that our method and results can help to select the right model for a particular situation, to relate other models to the present ones, to obtain an overview of the potential long-term dynamics of complex models, and to support model development, for example, by providing a rapid consistency check. Note that measured long-term as well as transient data can be explained with respect to organizations, by defining a projection from a system state to an organization, as demonstrated for the Handel model (Figure 7 and Figure 8).

In addition, our approach is not limited to IAV models, but can be directly applied to other viruses in the same way since their dynamics are similarly modeled by ODEs [32]. Furthermore, our approach is open to include other dynamically relevant components like treatment and vaccination strategies. Another task for future work would be to study the transient virus immune system in more detail, for example, by mapping the basin of attraction of each organization or by systematically analyzing the transition dynamics in-between organizations [33,34].