Dynamics of Virus and Immune Response in Multi-Epitope Network

TL;DR

This paper models the complex interactions between viruses and immune responses using differential equations, revealing how immune dominance influences persistent viral variants and network structure during infection.

Contribution

It introduces a mathematical framework for virus-immune dynamics, analyzing stability and persistence, especially in multi-epitope scenarios with binary resistance states.

Findings

Equal viral fitness costs lead to a nested network of persistent strains.

The number of persistent strains is at most n+1 for n epitopes.

Immunodominance determines the structure of persistent viral networks.

Abstract

The host immune response can often efficiently suppress a virus infection, which may lead to selection for immune-resistant viral variants within the host. For example, during HIV infection, an array of CTL immune response populations recognize specific epitopes (viral proteins) presented on the surface of infected cells to effectively mediate their killing. However HIV can rapidly evolve resistance to CTL attack at different epitopes, inducing a dynamic network of interacting viral and immune response variants. We consider models for the network of virus and immune response populations, consisting of Lotka-Volterra-like systems of ordinary differential equations. Stability of feasible equilibria and corresponding uniform persistence of distinct variants are characterized via a Lyapunov function. We specialize the model to a "binary sequence" setting, where for $n$ epitopes there can be $2^n$ distinct viral variants mapped on a hypercube graph. The dynamics in several cases are analyzed and sharp polychotomies are derived characterizing persistent variants. In particular, we prove that if the viral fitness costs for gaining resistance to each epitope are equal, then the system of $2^n$ virus strains converges to a "perfectly nested network" with less than or equal to $n+1$ persistent virus strains. Overall, our results suggest that immunodominance, i.e. relative strength of immune response to an epitope, is the most important factor determining the persistent network structure.

Figures

Tables

Table 1. Table 1: Stable equilibria values for healthy and infected cells in Theorem 4.4 where six regimes sharply characterize distinct viral strain persistence scenarios ( 𝟎 0 \bf{0} components go extinct) when both immune responses z 1 subscript 𝑧 1 z_{1} & z 2 subscript 𝑧 2 z_{2} persist (when ℛ 2 > 𝒬 2 subscript ℛ 2 subscript 𝒬 2 \mathcal{R}_{2}>\mathcal{Q}_{2} ). Note that if ℛ 2 ≤ Q 2 subscript ℛ 2 subscript 𝑄 2 \mathcal{R}_{2}\leq Q_{2} , then only z 1 subscript 𝑧 1 z_{1} and y 1 subscript 𝑦 1 y_{1} , y 2 subscript 𝑦 2 y_{2} can persist and the dynamics are detailed in Theorem 4.3 .

Strains	,		$ℛ={ℛ}_4$		${ℛ}_4\le {𝒫}_2,$
	${ℛ}_4\le {𝒬}_2$
$x^{\ast }$	$\frac{1}{{𝒬}_2}$	$\frac{1}{{ℛ}_4}$	$\frac{1}{ℛ}$	$\frac{1}{ℛ}$	$\frac{1}{{𝒫}_2}$	$\frac{1}{{ℛ}_4}$
$y_1^{\ast }\left(\mathrm{𝟎𝟎}\right)$	$s_1$	$s_1$	$\frac{{𝒫}_2}{ℛ}-1+y_4^{\ast }$	$\frac{{𝒫}_2}{ℛ}-1$	0	0
$y_2^{\ast }\left(\mathrm{𝟏𝟎}\right)$	$s_2-s_1$	$s_2-s_1$	$1-\frac{{𝒬}_2}{ℛ}+s_2-s_1-y_4^{\ast }$	$1-\frac{{𝒬}_2}{ℛ}+s_2-s_1$	$s_2$	$s_2$
$y_3^{\ast }\left(\mathrm{𝟎𝟏}\right)$	0	0	$1-\frac{{𝒬}_2}{ℛ}-y_4^{\ast }$	$1-\frac{{𝒬}_2}{ℛ}$	$s_1$	$s_1$
$y_4^{\ast }\left(\mathrm{𝟏𝟏}\right)$	0	$1-\frac{{𝒬}_2}{{ℛ}_4}$	$y_4^{\ast }$	0	0	$1-\frac{{𝒫}_2}{{ℛ}_4}$
Stable Equilib.	${\overline{ℰ}}_2$	${\widetilde{ℰ}}_2$	$ℰ\left(y_4^{\ast }\right),0\le y_4^{\ast }\le 1-\frac{{𝒬}_2}{ℛ}$	${\widehat{ℰ}}_2$	${ℰ}_2^{†}$	${ℰ}_2^{‡}$

Equations

$$\frac{dX}{dt}$$

$$\frac{d{Y}_i}{dt}$$

$$\frac{d{Z}_j}{dt}$$

$$x$$

$$a_{ij}$$

$${\mathcal{R}}_i$$

$$\dot{x}$$

$${\dot{y}}_i$$

$${\dot{Z}}_j$$

$${\Lambda }_i:=\left\{j\in [1,n]:a_{ij}>0\right\},$$

$$\forall i\in [1,m]:a_{ij}=a_j>0\text{if}j\in {\Lambda }_i,$$

$${\mathcal{I}}_j:=\frac{a_j}{{\rho }_j}.$$

$$\dot{x}$$

$$\text{If}{y}_i\leftrightarrow y_k\text{and}d(y_i,y_w)<d(y_k,y_w),\text{then}{\mathcal{R}}_i>{\mathcal{R}}_k.$$

$${\Omega }_y=\left\{i\in [1,m]:y_i^{\ast }>0\right\}\text{and}{\Omega }_z=\left\{j\in [1,n]:Z_j^{\ast }>0\right\}.$$

$$\Omega$$

$$\sum _{i\in {\Omega }_y}{a}_{ij}{y}_i^{\ast }$$

$$\left(1+\sum _{i\in {\Omega }_y}{\mathcal{R}}_i{y}_i^{\ast }\right)$$

$$\sum _{j\in {\Omega }_z}{a}_{ij}{Z}_j^{\ast }$$

$${\mathcal{R}}_i{x}^{\ast }-1-\sum _{j\in {\Omega }_z}{a}_{ij}{Z}_j^{\ast }\le 0,\forall i\notin {\Omega }_y,\sum _{i\in {\Omega }_y}{a}_{ij}{y}_i^{\ast }-{\rho }_j\le 0,\forall j\notin {\Omega }_z.$$

$$y^{\prime }{A}^{\prime }={\rho }^{\prime },A^{\prime }{Z}^{\prime }=\frac{{\vec{\mathcal{R}}}^{\prime }}{1+y^{\prime }{\vec{\mathcal{R}}}^{\prime }}-\vec{1},$$

$$(y^{\prime }-y^{\prime \prime })A^{\prime }$$

$$\Rightarrow 0$$

$$\Rightarrow x^{\prime }$$

$$\vec{r}{}^{\prime }$$

$$y^{\prime }-y^{\prime \prime }\in Ker(A^{\prime }{)}^T\cap {\vec{\mathcal{R}}}^{\prime \perp }and{Z}^{\prime }-Z^{\prime \prime }\in Ker{A}^{\prime }.$$

$$\mathrm{Ker}(A^{\prime }{)}^T\cap {\vec{\mathcal{R}}}^{\prime \perp }=\left\{0\right\}and\mathrm{Ker}(A^{\prime })=\left\{0\right\}$$

$$0$$

$$\exists ϵ,M>0\text{and}T({\vec{w}}_0)\text{such that}M>y_i(t),Z_j(t)>ϵ,i\in {\Omega }_y,j\in {\Omega }_z,\forall t>T({\vec{w}}_0),\text{and}$$

$$\underset{t\to \infty }{\lim }{y}_i(t),Z_j(t)=0,i\notin {\Omega }_y,j\notin {\Omega }_z,\text{for every solution with initial condition}{\vec{w}}_0\in \Omega .$$

$$\sum _{i\in {\Omega }_y}{\mathcal{R}}_i{y}_i$$

$${\dot{y}}_i$$

$${\dot{Z}}_j$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}\underset{0}{\overset{t}{\int }}{y}_i(s)ds=y_i^{\ast },\underset{t\to \infty }{\lim }\frac{1}{t}\underset{0}{\overset{t}{\int }}{Z}_j(s)ds=Z_j^{\ast },$$

$$W(x,y,Z)$$

$$:=W_1+W_2+W_3,$$

$${\dot{W}}_1$$

Full text

1 Introduction

The dynamics of virus and immune response within a host can be viewed as a complex ecological system. Both predator-prey and competitive interactions are especially important during a host infection. The immune response predates on the pathogen, and distinct viral strains compete for a target cell population, while immune response populations compete for the virus since their proliferation occurs upon pathogen recognition. The immune response can cause significant mortality of the virus, which may lead to selection for immune-resistant viral variants within the host. For example during HIV infection, an extensive repertoire of CTL immune effectors recognize specific epitopes (viral proteins) presented on the surface of infected cells to effectively mediate their killing, however HIV can rapidly evolve resistance to CTL attack at different epitopes. The ensuing battle precipitates a dynamic network of interacting viral strains and immune response variants, analogous to an ecosystem of rapidly evolving prey countering attack from a diverse collection of predators.

While the virus-immune interactions may be quite complex, patterns and structure can emerge. For the cellular immune response, a consistent and reproducible hierarchy of T cell populations organize in response to multiple epitopes of a pathogen, according to their (vertical) immunodominance, i.e. relative expansion levels of the responding immune populations within the host Kloverpris . Vertical T cell immunodominance patterns are highly variable among HIV infected individuals and change over time, largely due to sequence variability in the viral “quasispecies” liu2013vertical . Rapidly evolving pathogens, such as HIV and HCV, can evade the immune response via mutations at multiple epitopes. The pattern of epitope mutations, called the escape pathway, is of significant interest, and there is some evidence that the viral evolution is predictable Barton . The fitness of an emerging viral mutant strain, along with the strength of the CTL response, certainly affect the selection pressure for a single epitope mutation Vitaly3 . However the concurrent interaction of diverse virus and immune response populations necessitate considering the whole system together in order to understand viral escape of multiple epitopes liu2013vertical . In this paper, we introduce and analyze mathematical models for the dynamics of virus and immune response in a network determined by interaction at multiple epitopes.

A large amount of work on modeling within-host virus dynamics has been based on the “standard” virus model; an ordinary differential equation system describing the coupled changes in target cells, infected cells, and free virus particles through time in an infected individual Perelson2 . The CTL immune response has been included in variations of the standard virus model by considering an immune effector population which kills and is activated by infected cells according to a mass-action (bilinear) rate, although other functional forms for the activation rate have been utilized browne2016 ; deBoer . Nowak et al. nowak1996population , along with other subsequent works bobko2015singularly ; iwasa2004some ; souza2011global , have considered the dynamics of multiple virus strains which are attacked by strain-specific CTL immune response populations (“one-to-one” virus-immune network). However the assumption of strain-specific immune response does not correspond to the biological reality that CTLs are specific to epitopes and, in general, multiple epitopes will be shared among virus strains.

Multi-epitope models have been utilized with different datasets cataloging several epitope specific CTL response and viral escape mutations, in order to quantify escape rates and patterns liu2013vertical ; Leviyang ; pandit2014reliable . Earlier work often considered escape dynamics at epitopes separately, but some recent work has emphasized the concurrent interaction of distinct CTLs with the virus at multiple epitopes. Ganusov et al. Vitaly1 ; Vitaly2 explicitly include multiple CTL clones specific to different epitopes in the standard virus model, and utilize statistical approaches on a linearized version of the model to estimate rates of escape. Althaus et al. and van Deutekom et al. have also considered multiple epitopes and viral strains in the standard virus model, although results were mostly based on stochastic simulations Althaus ; vanDeutekom .

Browne browne2016global recently analyzed the stability and uniform persistence of a multi-epitope virus-immune model with a perfectly nested interaction network. The model setup mirrors a tri-trophic chemostat ecosystem with a single resource (healthy cells), and a network of consumers (viral strains) and their predators (immune variants). The perfectly nested network constrains the viral escape pathway so that resistance to multiple epitopes is built sequentially in the order of the immunodominance hierarchy. The successive rise of more broadly resistant prey (coming with a fitness cost) and weaker but more generalist predators, in a perfectly nested fashion, is the route to persistence of nested bacteria-phage communities argued in jover2013mechanisms ; korytowski2015nested . The analysis of more complex interaction networks, which allow arbitrary viral escape pathways for building multi-epitope resistance, will be addressed in this paper.

Here we extend the previous work by analyzing a within-host virus model with a general interaction network of multiple variants of virus and immune response. Our results suggest that a diverse viral quasispecies is constructed by resistance mutations at multiple epitopes and the immunodominance hierarchy is the main factor shaping the viral escape pathway. These notions are supported by observations in HIV infection. Indeed, the efficacy and breadth of cognate CTL immune responses increase within-host HIV diversity, driving viral evolution so that different combinations of multiple epitope escapes become prevalent in the viral population pandit2014reliable . Also, recent studies have shown that immunodominance hierarchies in HIV are major determinants of viral escape from multiple epitopes Barton ; liu2013vertical , in particular immunonodominance was found to play a substantially larger role than the viral fitness costs and other factors liu2013vertical . Understanding the main factors shaping viral escape pathways and immune dynamics is important for design of effective vaccines and immunotherapies.

The outline of this paper is as follows. In Section 2, we formulate our general virus-immune network model. We also introduce the “binary sequence” example motivated by HIV and CTL immune response dynamics, where a viral strain is either completely susceptible $(0)$ or has evolved complete resistance $(1)$ to immune attack at a specific epitope, in which for $n$ epitopes, there can be $2^{n}$ distinct viral variants distinguished by their immune resistance profile. In Section 3, we characterize the structure of feasible equilibria in the general model, along with finding a Lyapunov function for stability and corresponding uniform persistence of distinct variants. Next, in Section 4 we analyze the binary sequence example, deriving some graph-theoretic properties of feasible equilibria, and classifying dynamics in several special cases. In particular, if we constrain the virus-immune response network to be “perfectly nested” (Sec. 4.1), “strain-specific” (Sec. 4.2), or have $n=2$ epitopes (Sec. 4.3), sharp polychotomies characterize persistent variants. In Section 4.4 we prove that if the viral fitness costs for gaining resistance to each epitope are equal, then the system of $2^{n}$ virus strains converges to a perfectly nested network with less than or equal to $n+1$ persistent virus strains.

2 Mathematical model

We consider the following general virus-immune dynamics model, as in Browne browne2016global , which includes a population of target cells ($X$), $m$ competing virus strains ($Y_{i}$ denotes strain $i$ infected cells), and $n$ variants of immune response ($Z_{j}$):

[TABLE]

The function $f(X)=b-cX$ represents the net growth rate of the uninfected cell population. The parameter $\beta_{i}$ is the infection rate and $\delta_{i}$ is the decay rate for infected cells infected with virus strain $i$. The parameter $\mu_{j}$ denotes the decay rate of the immune response population $j$. We assume immune killing and activation rates are mass-action, representative of these events occurring as immune response cells recognize epitopes on the surface of infected cells. The parameter $r_{ij}$ describes the killing/interaction rate of immune population $Z_{j}$ on a strain-$i$ infected cell, whereas $q_{j}r_{ij}$ describes the corresponding activation rate for $Z_{j}$ (proportional to interaction rate $r_{ij}$). In the present paper, we assume that virus load (the abundance of virions) is proportional to the amount of (productively) infected cells. This assumption has frequently been made for HIV since the dynamic of free virions occurs on a much faster time scale than the other variables.

The model can be rescaled by introducing the following quantities:

[TABLE]

The model then becomes:

[TABLE]

Here $\mathcal{R}_{i}$ represents the basic reproduction number of virus strain $i$. Note that $\rho_{j}$ represents the reciprocal of the immune response fitness excluding the (rescaled) avidity to each strain $j$. The $m\times n$ nonnegative matrix $A=\left(a_{ij}\right)$ describes the virus-immune interaction network, which determines each immune effector population’s avidity to the distinct viral strains. Each virus strain $i$ (cells infected with strain $i$), $y_{i}$, has a set of CTLs, $z_{j}$, that recognize and attack $y_{i}$. We call this set the epitope set of $y_{i}$, denoted by $\Lambda_{i}$, where

[TABLE]

i.e. $j\in\Lambda_{i}$, if $y_{i}$ is not completely resistant to CTL $Z_{j}$.

It is not hard to show that solutions to (2) remain non-negative and bounded for all time $t$ browne2016global . In Section 3, we analyze feasible equilibria and their stability in the model (2). Note that the general interaction network in (2) allows for cross-reactivity, i.e. targeting of multiple epitopes by an immune population $Z_{j}$. However in Section 4, we focus on a special case of the interaction network for epitope specific CTL immune responses, which we introduce below.

2.1 Example: “Binary Sequence Case”

Here we specialize system (2) in order to model viral escape from multiple epitopes targeted by specific CTL immune responses, as occurs during HIV infection Barton ; liu2013vertical . Suppose there are $n$ viral epitopes, each one recognized by their corresponding specific CTL variant population from the set $Z_{1},\dots Z_{n}$. To model the viral escape pathway, we consider two possible alleles for each epitope: the wild type (0) and the mutated type (1) which has escaped recognition from the cognate immune response. For each virus strain $y_{i}$, we associate a binary sequence of length $n$, $y_{i}\sim\mathbf{i}=\left(i_{1},i_{2},\dots,i_{n}\right)\in\left\{0,1\right\}^{n}$, coding the allele type at each epitope. We assume that each immune response ($Z_{j}$) targets its specific epitope at the specific rate $a_{j}$ for virus strains containing the wild-type (allele 0) epitope $j$, whereas $Z_{j}$ completely loses ability to recognize strains with the mutant (allele 1) epitope $j$, i.e.

[TABLE]

where $\Lambda_{i}$ is the strain $i$ (susceptible) epitope set defined earlier for model (2) by (3) (see Fig. 1). For example, the wild-type (founder) virus strain, denoted here by $y_{w}$, is represented by the sequence of all zeroes and epitope set $\Lambda_{w}=\left\{1,\dots,n\right\}$ since it is susceptible to attack by all immune responses. With assumption (4), we can define an immune reproduction number corresponding to each $Z_{j}$:

[TABLE]

Furthermore there are $m=2^{n}$ possible viral mutant strains distinguished by reproduction number $\mathcal{R}_{i}$ and epitope set $\Lambda_{i}$ (or equivalently the binary string $\mathbf{i}\in\left\{0,1\right\}^{n}$). System (2) can be rescaled as:

[TABLE]

where $z_{j}=a_{j}Z_{j}$, $s_{j}=1/\mathcal{I}_{j}$, $i=1,\dots,m\leq 2^{n}$ and $j=1,\dots,n$.

Although mutation rates are not explicitly included in model (6), the $2^{n}$ potential virus strains can be viewed in a mutational pathway network; each strain is a vertex in an $n$-dimensional hypercube graph, shown in Figs. 1 and 1 in the case of $n=3$ and $n=2$ epitopes. Viral strains $y_{i}\sim\mathbf{i}=(i_{1},\dots,i_{n})$ and $y_{k}\sim\mathbf{k}=(k_{1},\dots k_{n})$ are connected by an edge, which we denote $y_{i}\leftrightarrow y_{k}$, if the sequences $\mathbf{i}$ and $\mathbf{k}$ differ in exactly one bit, i.e. their Hamming distance – denoted by $d(y_{i},y_{k})$ – is one. In this way, $y_{i}$ can mutate into $y_{k}$ (or vice-versa) through a single epitope mutation if $y_{i}\leftrightarrow y_{k}$. Since each mutation of an epitope comes with a fitness cost, we assume that

[TABLE]

We say that an immune response $z_{j}$ is immunodominant over another immune response $z_{k}$ if $\mathcal{I}_{j}>\mathcal{I}_{k}$ and assume without loss of generality the ordered immunodominance hierarchy; $\mathcal{I}_{1}\geq\mathcal{I}_{2}\geq\dots\geq\mathcal{I}_{n}$.

System (6) generalizes many previous model structures in the sense that they can be seen as subgraphs of our “hypercube network”. For instance, the “strain-specific” (virus-immune response) network nowak1996population (also called “one-to-one network” in phage-bacteria models jover2013mechanisms ; korytowski2015nested ) is equivalent to restricting (6) to the $m=n$ viral strains which have mutated $n-1$ epitopes (Figure 1). The “perfectly nested network” restricts (6) to the $m=n+1$ viral strains which have sequential epitope mutations in the order of the immunodominance hierarchy (Figure 1). Nested networks were considered in HIV models browne2016global , along with phage-bacteria models korytowski2015nested , and may be a common persistent structure in ecological communities gurney2017network . The “full hypercube network” has been considered for modeling CTL escape patterns in HIV infected individuals Althaus ; vanDeutekom . The dynamics of system (6) with the full network on $n$ epitopes (consisting of $m=2^{n}$ viral strains), along with one-to-one and perfectly nested subgraphs (consisting of $m=n$ or $m=n+1$ strains), will be analyzed in Section  4.

3 Equilibria and Lyapunov function

A general non-negative equilibrium point, $\mathcal{E}^{*}=\left(x^{*},y^{*},Z^{*}\right)\in\mathbb{R}_{+}^{1+m+n}$, of system (2) will be characterized in terms of the positive virus and immune variant components. Define the “persistent variant sets” associated with $\mathcal{E}^{*}$ as:

[TABLE]

In addition, define the following subsets of $\mathbb{R}_{+}^{1+m+n}$:

[TABLE]

Here $\Gamma_{\Omega}$, consisting of only those state vectors having the same set of positive and zero components as equilibrium $\mathcal{E}^{*}$, is called the positivity class of $\mathcal{E}^{*}$. Notice that the dimension of the subset $\Gamma_{\Omega}$ is $1+|\Omega_{y}|+|\Omega_{z}|$, where the notation $|\Omega_{y}|$ ($|\Omega_{z}|$) denotes the cardinality of the set $\Omega_{y}$ ($\Omega_{z}$). The equilibrium $\mathcal{E}^{*}$ must satisfy the following equations:

[TABLE]

We note that $\mathcal{R}_{i}>1,\ i\in\Omega_{y}$ must hold, even in the absence of CTL response.

Following Hofbauer and Sigmund hofbauer1998evolutionary , we call an equilibrium $\mathcal{E}^{*}=(x^{*},y^{*},Z^{*})$ of (2) saturated if the following holds when $\mathcal{E}^{*}$ has zero components:

[TABLE]

Note that if $\mathcal{E}^{*}$ has all positive components, i.e. $\Omega_{y}=[1,m]$ and $\Omega_{z}=[1,n]$, the inequalities (11) trivially hold. Note also that each term in (11) is an “invasion” eigenvalue of the Jacobian matrix evaluated at $\mathcal{E}^{*}$ and thus a saturated equilibrium enjoys a weak stability against invasion by missing species $i\notin\Omega_{y},\ j\notin\Omega_{z}$. It immediately follows that a stable equilibrium must be saturated since the Jacobian cannot have a positive eigenvalue. As part of Theorem 3.1 later in this section, we will conversely show that every saturated equilibrium is stable. First, the following proposition states that there exists at least one saturated equilibrium. For completeness we provide its proof in Appendix A.1, however we also note that this proposition follows directly from Theorem 2 in hofbauer1990index .

Proposition 1

There exists a saturated equilibrium of system (2).

More properties of relevant equilibria can be ascertained. The following proposition states that any equilibria in the same positivity class must share the same value at $x^{*}$.

Proposition 2

If $\mathcal{E}^{\prime}=\left(x^{\prime},y^{\prime},Z^{\prime}\right)$ and $\mathcal{E}^{\prime\prime}=\left(x^{\prime\prime},y^{\prime\prime},Z^{\prime\prime}\right)$ are both equilibria in the same positivity class, $\Gamma_{\Omega}$, then $x^{\prime}=x^{\prime\prime}$.

Proof

Let $A^{\prime}$ denote the submatrix of $A$ which contains only the rows in $\Omega_{y}$ and columns in $\Omega_{z}$. Then the equilibrium conditions for $\mathcal{E}^{\prime}$ can be rewritten as:

[TABLE]

where $\rho^{\prime}$ is the row vector with components $\rho_{j}$ where $j\in\Omega_{z}$, $\vec{\mathcal{R}}^{\prime}$ is the column vector with components $\mathcal{R}_{i}$ where $i\in\Omega_{y}$ and $\vec{1}$ is the vector with all components one. Since $\mathcal{E}^{\prime\prime}$ satisfies the same conditions, we obtain:

[TABLE]

The previous proposition implies that if an equilibrium $\mathcal{E}^{*}=\left(x^{*},y^{*},Z^{*}\right)$ exists in positivity class $\Gamma_{\Omega}$, then any equilibrium $\mathcal{E}^{\prime}=\left(x^{*},y^{\prime},Z^{\prime}\right)$ belonging to $\Gamma_{\Omega}$ will satisfy the following Lotka-Volterra equilibria conditions within $\Gamma_{\Omega}$:

[TABLE]

where $A^{\prime}$ is the submatrix of $A$ with rows in $\Omega_{y}$ and columns in $\Omega_{z}$, $\rho^{\prime}$ is the row vector with components $\rho_{j}$ where $j\in\Omega_{z}$, $\vec{\mathcal{R}}^{\prime}$ is the column vector with components $\mathcal{R}_{i}$ where $i\in\Omega_{y}$, and $\vec{v}\,^{\prime}=\left(y^{\prime},Z^{\prime}\right)^{T}$. Note that if the cardinality of $\Omega_{y}$ and $\Omega_{z}$ are equal ($|\Omega_{y}|=|\Omega_{z}|$) and $A^{\prime}$ is non-singular, then clearly $\mathcal{E}^{*}$ is unique in its positivity class, $\Gamma_{\Omega}$. More generally, if $B^{\prime}$ is non-singular, then $\mathcal{E}^{*}$ is unique in $\Gamma_{\Omega}$. The following proposition sharpens the condition for uniqueness of an equilibrium within a positivity class, and shows that in such equilibria the number of virus strains either is equal to or exactly one more than the number of immune responses.

Proposition 3

Suppose the equilibrium $\mathcal{E}^{*}=\left(x^{*},y^{*},Z^{*}\right)$ exists in positivity class $\Gamma_{\Omega}$, where $\left(y^{*},Z^{*}\right)$ satisfy the linear system of equations (12) and the cardinality of $\Omega_{y}$ and $\Omega_{z}$ are $|\Omega_{y}|=m^{\prime}$ and $|\Omega_{z}|=n^{\prime}$. Then $\mathcal{E}^{*}$ is the unique equilibrium in $\Gamma_{\Omega}$, i.e. $\vec{v}=\left(y^{*},Z^{*}\right)^{T}$ is the unique solution to (12), if and only if ${\rm Ker}(A^{\prime})^{T}\cap\vec{\mathcal{R}}^{\prime\perp}=\left\{0\right\}$ and ${\rm Ker}(A^{\prime})=\left\{0\right\}$.

Moreover, if $\mathcal{E}^{*}$ is the unique equilibrium in $\Gamma_{\Omega}$, then one of the following holds:

• (i)

$m^{\prime}=n^{\prime}$, and $x^{*}=1/\left(1+(\vec{\rho}\,^{\prime})^{T}(A^{\prime})^{-1}\vec{\mathcal{R}}^{\prime}\right)$.

• (ii)

$m^{\prime}=n^{\prime}+1$, and $x^{*}=\vec{1}^{\,T}C^{-1}_{(n^{\prime}+1)}$, where $C^{-1}_{(n^{\prime}+1)}$ is the last column in the $(n^{\prime}+1)\times(n^{\prime}+1)$ matrix inverse of $C=\begin{pmatrix}A^{\prime}&\vec{\mathcal{R}^{\prime}}\end{pmatrix}^{T}$.

Proof

Suppose there exists equilibrium $\mathcal{E}^{*}$ in positivity class $\Gamma_{\Omega}$ and consider the corresponding $m^{\prime}\times n^{\prime}$ submatrix $A^{\prime}$ from linear system (12), where $m^{\prime}=|\Omega_{y}|$ and $n^{\prime}=|\Omega_{z}|$. Proposition 2 shows that any two equilibria $\mathcal{E}^{\prime}$ and $\mathcal{E}^{\prime\prime}$, contained in the same positivity class $\Gamma_{\Omega}$, have equal $x^{*}$ components, i.e. $x^{\prime}=x^{\prime\prime}$. Thus if these equilibria are distinct, then $y^{\prime}\neq y^{\prime\prime}$ or $Z^{\prime}\neq Z^{\prime\prime}$. In the proof of Proposition 2 it is shown that

[TABLE]

Therefore, the condition

[TABLE]

is equivalent to uniqueness of an equilibrium $\mathcal{E}^{*}$ in its positivity class $\Gamma_{\Omega}$.

Moreover $\hbox{Ker}(A^{\prime})^{T}\cap\vec{\mathcal{R}}^{\prime\perp}=\{0\}$ if and only if the augmented $(n^{\prime}+1)\times m^{\prime}$ matrix $C$ consisting of adding the final row $\vec{\mathcal{R}}^{\prime T}$ to $(A^{\prime})^{T}$ has trivial kernel. By the rank-nullity theorem, we obtain that $C$ has rank equal to $m^{\prime}$. Since rank cannot exceed the number of rows, $m^{\prime}\leq n^{\prime}+1$. Applying the rank-nullity theorem to $A^{\prime}$, gives rank equal to $n^{\prime}$ and $n^{\prime}\leq m^{\prime}$. Thus $n^{\prime}\leq m^{\prime}\leq n^{\prime}+1$.

In the case (i), $m^{\prime}=n^{\prime}$, the matrix $A$ is invertible and so from the equilibria equations (10), we obtain $x^{*}=1/\left(1+(\vec{\rho}\,^{\prime})^{T}(A^{\prime})^{-1}\vec{\mathcal{R}}^{\prime}\right)$. Finally, consider case (ii), $m^{\prime}=n^{\prime}+1$. Since $n^{\prime}={\rm rank}(A^{\prime})={\rm rank}\left((A^{\prime})^{T}\right)$, by the rank-nullity theorem we obtain that ${\rm null}\left((A^{\prime})^{T}\right)=1$. We claim that $\hbox{Ker}(A^{\prime})^{T}$ contains a vector $\vec{w}$ such that $\vec{w}^{T}\vec{\mathcal{R}}^{\prime}=1$ and $\sum_{i}w_{i}=x^{*}$. Since the matrix $C$ (defined in previous paragraph) has trivial Kernel, it is invertible. Let $\vec{w}$ be the last column of $C^{-1}$. Then it is not hard to see that $(A^{\prime})^{T}\vec{w}=\vec{0}$ and $(\vec{\mathcal{R}}^{\prime})^{T}\vec{w}=1$. Furthermore

[TABLE]

Thus $x^{*}=\sum_{i}w_{i}=\vec{1}^{T}C^{-1}_{(n^{\prime}+1)}$.

Notice from the above proof that if an equilibrium $\mathcal{E}^{*}$ is not unique in its positivity class $\Gamma_{\Omega}$, then $\Gamma_{\Omega}$ contains an infinite number (a continuum) of equilibria. Conversely, if $\mathcal{E}^{*}$ is unique in a positivity class $\Gamma_{\Omega}$ containing $n^{\prime}$ (persistent) immune responses, then there are either (i) $n^{\prime}$ virus strains or (ii) $n^{\prime}+1$ virus strains in $\Gamma_{\Omega}$.

Additionally some results on existence of saturated equilibria in Lotka-Volterra systems can be recast in our setting. For example, a sufficient condition for $\mathcal{E}^{*}$ to be saturated (and unique equilibrium in $\Gamma_{\Omega}$) is if the matrix $-B^{\prime}$ in (12) is a $P$-matrix, i.e. all principal minors of $-B^{\prime}$ are positive, by Theorem 15.4.5 in hofbauer1998evolutionary .

In what follows, we will be interested in the global behavior of solutions to system (2). In doing so, we will determine which viral strains and immune responses uniformly persist thieme1993persistence and which go extinct. Define the system to be $\Omega_{yz}$ permanent if

[TABLE]

We will sometimes use the terminology that $y_{i},Z_{j}\ i\in\Omega_{y},j\in\Omega_{z}$ are uniformly persistent and $y_{i},Z_{j}\rightarrow 0\ i\notin\Omega_{y},j\notin\Omega_{z}$ to signify the system being $\Omega_{yz}$ permanent.

In the spirit of permanence as a sufficient condition for existence of a unique interior rest point in Lotka-Volterra systems hofbauer1998evolutionary , we find the following proposition.

Proposition 4

If the system is $\Omega_{yz}$ permanent, then there is a unique equilibrium $\mathcal{E}$ in the positivity class $\Gamma_{\Omega}$.

Proof

By Theorem 6.2 in smith2011dynamical , $\Omega_{yz}$ permanence implies that there exists an equilibrium in the positivity class $\Gamma_{\Omega}$. Suppose by way of contradiction that there are two equilibria, $\mathcal{E}^{\prime}$ and $\mathcal{E}^{\prime\prime}$, in $\Gamma_{\Omega}$. By Proposition 2, $x^{\prime}=x^{\prime\prime}$. Since the remaining equilibria equations (10) are linear, it can be shown that the line through $\mathcal{E}^{\prime}$ and $\mathcal{E}^{\prime\prime}$ consist entirely of equilibria. Then, we can find equilibria arbitrarily close to the boundary of $\Gamma_{\Omega}$. This contradicts the fact that the system is $\Omega_{yz}$ permanent.

Now we state the main theorem of this section concerning the persistence of viral and immune variants of model (2). It builds on results by Browne browne2016global concerning a special case of (2), namely the case of perfectly nested subnetwork (see Fig. 1) in the “binary sequence” system (6). In particular, the notion of saturated equilibria allows us to significantly extend persistence and stability results to the arbitrary interaction networks in general model (2).

Theorem 3.1

Suppose that $\mathcal{E}^{*}=\left(x^{*},y^{*},Z^{*}\right)$ is a non-negative equilibrium of system (2) with positivity class $\Gamma_{\Omega}$. Suppose further that $\mathcal{E}^{*}$ is saturated, i.e. the inequalities (11) hold. Then $\mathcal{E}^{*}$ is locally stable and $x(t)\rightarrow x^{*}$ as $t\rightarrow\infty$.

Furthermore, if $\mathcal{E}^{*}$ is the unique equilibrium in its positivity class $\Gamma_{\Omega}$ and the inequalities (11) are strict, then $y_{i},Z_{j}\rightarrow 0$ for all $i\notin\Omega_{y},j\notin\Omega_{z}$. If $i\in\Omega_{y}$ and $a_{ij}=0\ \forall j\in\Omega_{z}$, i.e. $\Lambda_{i}\cap\Omega_{z}=\emptyset$, then $y_{i}\rightarrow y_{i}^{*}$ and $x^{*}=1/\mathcal{R}_{i}$. In addition, omega limit sets corresponding to positive initial conditions are contained in invariant orbits satisfying

[TABLE]

and for each $i\in\Omega_{y},j\in\Omega_{z}$, $y_{i}$ and $Z_{j}$ persist (the system is $\Omega_{yz}$ permanent) with asymptotic averages converging to equilibria values, i.e.

[TABLE]

In the case that there are less than or equal to two persistent viral strains with non-empty epitope sets (restricted to $\Omega_{z}$), i.e. $|\left\{i\in\Omega_{y}:\Lambda_{i}\cap\Omega_{z}\neq\emptyset\right\}|\leq 2$, then $\mathcal{E}^{*}$ is globally asymptotically stable.

Proof

Consider the following Lyapunov function:

[TABLE]

where the term with logarithm should be omitted if the corresponding coordinate in the particular equilibrium is zero. Then taking the time derivatives, we obtain

[TABLE]

Thus

[TABLE]

where we make use of equilibrium conditions (10) in the final line. By the assumed inequalities (11), we obtain that $\dot{W}\leq 0$, and thus $\dot{W}$ is a Lyapunov function at the equilibrium $\mathcal{E}^{*}$. Noting that $\mathcal{E}^{*}$ is the unique minimizer of $W$, we obtain that $\mathcal{E}^{*}$ is (locally) stable. Additionally, since $\dot{W}\leq 0$ and $W\rightarrow\infty$ as $y_{i}$, $z_{j}$ goes to [math] or $\infty$ for $i\in\Omega_{y},j\in\Omega_{z}$, we find that for any solution there exists $p,P>0$ such that $p\leq y_{i},z_{j}\leq P$ for $i\in\Omega_{y},j\in\Omega_{z}$. Applying La Salle’s Invariance principle (See e.g. Thm 2.6.1 in hofbauer1998evolutionary ), the $\omega$-limit set corresponding to any solution of (2) with positive initial conditions is contained in the largest invariant set, $\mathcal{L}$, where $\dot{W}=0$. Clearly $\dot{W}=0\Rightarrow x=x^{*}$, thus $x=x^{*}$ in $\mathcal{L}$. This implies that $x(t)\rightarrow x^{*}$ as $t\rightarrow\infty$ for all solutions with positive initial conditions. Furthermore if the inequalities (11) are strict, then $y_{i}=0,z_{j}=0$ in $\mathcal{L}$ for $i\notin\Omega_{y},j\notin\Omega_{z}$. This implies that omega limit sets corresponding to positive initial conditions are contained in invariant orbits satisfying (14) and (13). The differential equations in (14) can be integrated to obtain asymptotic averages as follows:

[TABLE]

Since $\mathcal{E}^{*}$ is the unique equilibrium in $\Gamma_{\Omega}$, utilizing a similar argument as hofbauer1998evolutionary , we find that for each $i\in\Omega_{y},j\in\Omega_{z}$, solutions $y_{i},Z_{j}$ in the invariant set $\mathcal{L}$ satisfy:

[TABLE]

For any solution:

[TABLE]

Therefore, $y_{i}$ and $Z_{j}$, $i\in\Omega_{y},j\in\Omega_{z}$, are uniformly weakly persistent. The key compactness hypotheses of Corollary 4.8 from smith2011dynamical are satisfied, and thus weak uniform persistence implies strong uniform persistence for these variants. In addition, note that if $i\in\Omega_{y}$ and $a_{ij}=0\ \forall j\in\Omega_{z}$, i.e. $\Lambda_{i}\cap\Omega_{z}=\emptyset$, then $\dot{y}_{i}=0$ and therefore $y_{i}=y_{i}^{*}$ for this component on the invariant set $\mathcal{L}$ by (17).

Finally, we show that if $|\left\{i\in\Omega_{y}:\Lambda_{i}\cap\Omega_{z}\neq\emptyset\right\}|\leq 2$, then $\mathcal{E}^{*}$ is globally asymptotically stable. Without loss of generality suppose that $\left\{i\in\Omega_{y}:\Lambda_{i}\cap\Omega_{z}\neq\emptyset\right\}=\left\{1,2\right\}$. Then $\mathcal{R}_{1}y_{1}+\mathcal{R}_{2}y_{2}=\mathcal{R}_{1}y_{1}^{*}+\mathcal{R}_{2}y_{2}^{*}$ by (13) since $y_{i}=y_{i}^{*}\ i>2$ on the invariant set $\mathcal{L}$ for any other strains where $\Lambda_{i}\cap\Omega_{z}=\emptyset$. Now clearly there exists time $t_{1}$ such that $y_{1}(t_{1})=y_{1}^{*}$ by (17), and this implies that $y_{2}(t_{1})=y_{2}^{*}$ by the previous sentence. Without loss of generality assume $t_{1}=0$. Differentiating (13) twice and evaluating at time $t=0$, we obtain:

[TABLE]

where the last equality comes from the fact that $y_{k}(0)=y_{k}^{*}$ for $k=1,2$, and $a_{kj}=0$ for $k>2$. The only way the above equation can be satisfied is if $Z_{j}(0)=Z_{j}^{*}$ for all $j$. Then $y_{i}(0)=y_{i}^{*},\ Z_{j}(0)=Z_{j}^{*}$, $i\in\Omega_{y},j\in\Omega_{z}$. Thus $\mathcal{E}^{*}$ is globally asymptotically stable in this case.

More general results can be obtained in special cases, in particular, the inequalities (11) need not be strict for persistence results in certain cases discussed in Section 4. However, the global convergence to persistent variants is still an open question when there are more than two persistent immune responses. Note that the proof for global stability of “strictly saturated” equilibria $\mathcal{E}^{*}$ with $|\left\{i\in\Omega_{y}:\Lambda_{i}\cap\Omega_{z}\neq\emptyset\right\}|\leq 2$ does not extend to higher dimensional equilibria. Our numerical simulations that we have conducted support the global stability of $\mathcal{E}^{*}$ in general. We conjecture the following:

Conjecture 1

An equilibrium $\mathcal{E}^{*}$ of (2), which is unique in its positivity class $\Gamma_{\Omega}$ and satisfies inequalities (10) strictly, is globally asymptotically stable for positive initial conditions (regardless of the dimension of $\Gamma_{\Omega}$).

We briefly discuss Theorem 3.1 in the context of Lotka-Volterra (L-V) systems. L-V systems take the general form

[TABLE]

The multi-trophic version of model (18) with $m$ prey and $n$ predators may be set up with similar prey-predator interaction rates as our model (see (12)). The main difference between our chemostat-type system (2) and the L-V model is the explicit inclusion of the resource, $x$, modulating competition in the “prey” species (virus strains) in (2), as opposed to possible logistic competition terms, $b_{ij}<0,\ 1\leq i,j\leq m$ in L-V models. Theorem 3.1 reduces the dynamics of our model (2) (in the case of a “$\Gamma_{\Omega}$-unique” strictly saturated equilibrium) to the invariant L-V system (14) with $|\Omega_{y}|$ prey and $|\Omega_{z}|$ predators subject to the constraint (13) on prey. Global stability of equilibria in L-V systems have been studied extensively (see e.g. takeuchi1996global ). For instance, global stability of distinct equilibria was proved in the case of $m=n$ ($N=2n$ in (18)) and uniform competition coefficient for prey ($b_{ij}=b\ 1\leq i,j\leq n$) for a perfectly nested network korytowski2017persistence . For general LV systems, Goh proved that a strictly saturated equilibrium $\mathbf{v}^{*}$ is globally asymptotically stable for system (18) if there exists a positive diagonal matrix $D$ such that $DB+B^{T}D$ is negative definite (except possibly at equilibrium $\mathbf{v}^{*}$) Goh . Unfortunately none of these global stability results on L-V models apply to the invariant system (14), in part due to the lack of intraspecific competition coefficients in (14). Indeed the simple one predator-one prey L-V planar system of the form (14) displays oscillatory solutions in absence of the additional constraint (13). The constraint (13) seems to induce saturated equilibria to be attractors, leading to our conjecture (Conjecture 1) stating that in Theorem 3.1, variants $y_{i},z_{j}\ i\in\Omega_{y},j\in\Omega_{z}$ are not only uniformly persistent but also globally converge to their equilibrium values.

4 Special Cases of Multi-Epitope Model

In this section we consider the “binary sequence” case of model (2), which was introduced in Section 2.1. Recall that the assumption (4) leads to the rescaled and simplified system (6). Biologically, we are considering the situation where $n$ immune response populations $z_{j}$ each target the corresponding epitope $j$ in the virus strains at a rate solely dependent on the allele type of epitope $j$; (0) wild-type or (1) mutated form conferring full resistance to $z_{j}$. The avidity of immune response $z_{j}$ and (wild-type) epitope $j$ is described by the immune reproduction number $\mathcal{I}_{j}$ given by (5). The resulting model (6) has $m=2^{n}$ potential viral strains distinguished by their (basic) reproduction number, $\mathcal{R}_{i}$, and epitope set, $\Lambda_{i}$. Recall that in model (6) the viral strains can be viewed in a mutational pathway network, where each strain, $y_{i}\ i=1,\dots,2^{n}$, is a vertex in the hypercube graph, $Q_{n}$, corresponding to the strain’s epitope set, $\Lambda_{i}$, represented as a binary sequence of length $n$, $\mathbf{i}=(i_{1},\dots,i_{n})$. See Section 2.1 for relevant explanation and definitions.

We can establish restrictions on the positivity class of feasible equilibria in model (6) based on graph-theoretic considerations of the viral strains viewed as vertices in the hypercube graph. Indeed, we will show that equilibria with persistent viral strains forming a cycle of order $2^{j}$ can only occur in the degenerate cases. A (simple) cycle is defined as a closed path in the hypercube graph from containing no other repetition of vertices other than the starting and ending vertex. For $2\leq j\leq n$, there are $2^{n-j}$ disjoint $2^{j}$-cycles which cover the vertices of the hypercube graph $Q_{n}$. In particular, it is well-known that there is a Hamiltonian cycle covering $Q_{n}$, i.e. a cycle that includes every edge in the graph. The following proposition concerns feasibility of equilibria with cycles in the viral hypercube graph (here the important aspect is the set of vertices on the cycle).

Proposition 5

For $n\geq 2$, consider model (6). Let $2\leq j\leq n$ and $\ell=2^{j}$. Suppose there is a simple $\ell$-cycle in the representative hypercube graph $Q_{n}$: $y_{k_{1}}\leftrightarrow y_{k_{2}}\leftrightarrow\dots\leftrightarrow y_{k_{\ell}}\leftrightarrow y_{k_{1}}$. If $\sum_{i=1}^{\ell}(-1)^{i}\mathcal{R}_{k_{i}}=\sum_{i=1}^{\ell}(-1)^{d(y_{k_{i}},y_{w})}\mathcal{R}_{k_{i}}\neq 0$, then there does not exist an equilibrium with $y^{*}_{k_{1}},y^{*}_{k_{2}},\dots,y^{*}_{k_{\ell}}>0$.

Proof

Fix $2\leq j\leq n$ and let $\ell=2^{j}$. Without loss of generality, consider an $\ell$-cycle with strains labeled as $y_{1}\leftrightarrow\dots\leftrightarrow y_{\ell}\leftrightarrow y_{1}$. The simple cycle of length $\ell$ can be seen as an embedded hypercube graph on $2^{j}$ vertices, $Q_{j}$, corresponding to binary strings where $n-j$ slots are fixed to be all zeros or all ones. Without loss of generality suppose that these fixed slots are indices $[j+1,n]$ corresponding to $z_{j+1},\dots z_{n}$, i.e. either $[j+1,n]\subseteq\Lambda_{i}$ or $[j+1,n]\cap\Lambda_{i}=\emptyset$ for all $i=1,\dots,\ell$. The viral strains can be grouped into the classes $\mathcal{C}_{k}$, $k=0,\dots j$, where $k$ is the number of mutated alleles (1) appearing in the epitope set $[1,j]$, i.e. $y_{i}\in\mathcal{C}_{k}$ if $k=j-|\Lambda_{i}\cap[1,j]|$. Notice that in $\mathcal{C}_{k}$ there are exactly $\binom{j-1}{k}$ viral strains, $y_{i}$, with epitope 1 being wild-type allele (0), i.e. $1\in\Lambda_{i}$. The analogous statement holds for any epitope in $[1,j]$. Thus, for equilibria of (6), we find that

[TABLE]

Note here that $\sum_{i=1}^{\ell}(-1)^{i}\mathcal{R}_{i}=\sum_{i=1}^{\ell}(-1)^{d(y_{i},y_{w})}\mathcal{R}_{i}$ since we must traverse the simple cycle through “distance-one” single epitope mutations.

A couple remarks about the above proposition are in order. First, we highlight the of case $j=2$ of Proposition 5, which states the generic non-existence of equilibria with persistent viral strains forming a 4-cycle in the associated hypercube graph. In particular, if viral strains $y_{k_{1}},\dots,y_{k_{4}}$ form a cycle, then there exists an equilibrium $\mathcal{E}^{*}$ with $y_{k_{1}}^{*},\dots,y_{k_{4}}^{*}>0$ only if $\sum_{i=1}^{4}(-1)^{p_{i}}\mathcal{R}_{k_{i}}=0$ where $p_{i}=d(y_{k_{i}},y_{w})$. In Section 4.3, the degeneracy of “four-cycle equilibria” will be detailed further for the case of $n=2$ epitopes. In general, for $n\geq 2$, there are $2^{n-2}$ disjoint 4-cycles which cover the vertices of $Q_{n}$, and their unions form larger cycles with order as powers of two. Each of these cycles of order $2^{j},\ j\geq 2,$ can be seen as an embedded $j$-dimensional hypercube. Thus the proposition establishes the degeneracy of any equilibrium with positive viral components forming an embedded hypercube subgraph in the associated hypercube $Q_{n}$, which suggests that these combinations of viral strains generally do not persist together in model (6).

In following subsections, we analyze a few special cases of the multi-epitope model (6) where stable equilibria can be sharply characterized by quantities derived from the parameters using Theorem 3.1. We assume without loss of generality that immune responses $z_{1},\dots,z_{n}$ are ordered according to an immunodominance hierarchy:

[TABLE]

Section 4.1 summarizes previously published analysis of model (6) when the network is constrained to be perfectly nested, and Section 4.2 contains prior results, along with a new proposition, about the strain-specific (one-to-one) network. These cases provide nice applications of our general Theorem 3.1, and also are important for the analysis of the full network. We consider the full network with $m=2^{n}$ virus strains in section 4.3 for $n=2$ epitopes, and in section 4.4 for arbitrary $n$ in the special case of equal fitness costs for each epitope mutation.

4.1 Nested Network

While the virus-immune epitope interaction network generally can be quite complex, patterns of viral escape and dynamic immunodominance hierarchies often emerge. In observations of HIV infection, the initial CTL response occurs at a few immunodominant epitopes and is followed by viral mutations at these epitopes conferring resistance, along with a fall in these specific CTLs and rise in subdominant CTLs liu2013vertical . This pattern continues, albeit at diminishing rates as time proceeds, resulting in viral strains with resistance at multiple epitopes and corresponding fitness costs, along with subdominant CTLs of increasing breadth. An idealized description of this process is a perfectly nested network, where resistance to multiple epitopes is built sequentially according to the immunodominance hierarchy. Nested networks have been of recent interest in explaining the biodiversity and structure of bacteria-phage communities jover2013mechanisms ; korytowski2015nested ; weitz2013phage , and there is some evidence that nestedness is a feature of HIV-CTL dynamics kessinger2015inferring ; liu2013vertical ; vanDeutekom . Below we summarize recent work browne2016global , characterizing the stability of equilibria and uniform persistence or extinction of the populations for system (6) in the case of a perfectly nested network.

The perfectly nested network consists of $n$ epitope specific CTLs, $z_{1},\dots,z_{n}$, and $m=n+1$ virus strains $y_{1},\dots,y_{n+1}$ where the epitope set of $y_{i}$ is $\Lambda_{i}=\left\{i,\dots,n\right\}$ (having escaped immune responses $z_{1},\dots,z_{i-1}$). The equations are

[TABLE]

As before we assume a fitness cost for each mutation, and here we also assume a strict immunodominance:

[TABLE]

Out of a multitude of non-negative equilibria ($>2^{n}$), there are $2n+2$ which can be stable, and their stability depends upon quantities derived from parameters as follows. For $k\geq 1$ define:

[TABLE]

Then, for each $k\in[0,n]$, define the following equilibria:

[TABLE]

Equilibrium $\widetilde{\mathcal{E}}_{k+1}$ represents the appearance of escape mutant $y_{k+1}$ from equilibrium $\bar{\mathcal{E}}_{k}$.

The main result of browne2016global is as follows:

Theorem 4.1 (browne2016global )

If $\mathcal{R}_{1}>\mathcal{Q}_{1}$, let $k$ be the largest integer in $[1,n]$ such that $\mathcal{R}_{k}>\mathcal{Q}_{k}$, otherwise let $k=0$. If $\mathcal{R}_{k+1}\leq\mathcal{Q}_{k}$, then for $1\leq i\leq k$, $y_{i},z_{i}$ are uniformly persistent (if $\mathcal{R}_{k+1}>\mathcal{Q}_{k}$, $y_{k+1}$ persists), and the other variants globally converge to zero, $x(t)\rightarrow 1/\mathcal{Q}_{k}$ (if $\mathcal{R}_{k+1}>\mathcal{Q}_{k}$, $x(t)\rightarrow 1/\mathcal{R}_{k+1}$ and $y_{k+1}(t)\rightarrow 1-\frac{\mathcal{Q}_{k}}{\mathcal{R}_{k+1}}$) as $t\rightarrow\infty$. Additionally, the corresponding equilibria ($\widetilde{\mathcal{E}}_{k+1}$ or $\overline{\mathcal{E}}_{k}$) are locally stable (globally asymptotically stable when $k=0,1,2$), asymptotic averages converge to equilibria, and the global attractor satisfies (13) and (14).

Note that the proof can be obtained through application of Theorem 3.1.

Theorem 4.1 suggests a stable diverse set of viral strains and immune response which can be built up by the nested accumulation of epitope resistance and rise of subdominant CTLs. The diversity achieved depends upon potential breadth $n$ and “immune invasion” number at epitope $k\leq n$, $\mathcal{R}_{k}/\mathcal{Q}_{k}$, which depends upon the strengths of CTL directed at the $k$ epitopes in immunodominance hierarchy and the viral fitness costs of $k$ sequential mutations, along with initial fitness $\mathcal{R}_{1}$. Observe that since $\mathcal{R}_{k}$ decreases and $\mathcal{Q}_{k}$ increases with breadth $k$, Theorem 4.1 implies exclusion of $y_{k+1}$ is more likely as the breadth increases. Additionally, it is shown in browne2016global that the rate of $y_{k+1}$ invasion decreases as the breadth $k$ increases, which is consistent with several studies showing rate of HIV viral escape from CTL responses slows down after acute infection, along with relatively few escapes asquith2006inefficient ; Vitaly1 .

Overall, the analysis in the case of the nested network confirms some patterns of multi-epitope viral escape and reinforces the importance of strong immune responses directed at conserved epitopes (high fitness cost for resistance) in order to control HIV with CTL response. However, constraining multi-epitope resistance to be built in a nested fashion leaves out other potential mutational pathways.

The question remains, with $n$ epitopes targeted by distinct immune responses, what are all of the potential escape patterns and stable equilibria? We will answer this question for the simplest case of $n=2$ epitopes in Section 4.3, but first we consider the “strain-specific” (or one-to-one) subnetwork.

4.2 Strain-Specific network

Consider the case where $m=n+1$ with $\Lambda_{i}=\left\{i\right\}$ and $\Lambda_{n+1}=\emptyset$ in model (6), i.e. in system (2), $A$ is a $n+1\times n$ matrix comprised of the diagonal matrix ${\rm diag}\left(a_{1},\dots,a_{n}\right)$ and a row of zeros. This particular assumption of a “one-to-one” interaction network, where each immune response population attacks a unique specific viral strain, has been considered in wolkowicz1989successful ; korytowski2015nested ; bobko2015singularly . In this case, model (6) reduces to the following $n+1$-strain and $n$-immune variant model:

[TABLE]

Note that virus strain $y_{n+1}$ is immune to all immune cells. First we summarize previous results on system (24), which all can be derived from our general analysis in Section 3. For these results, we impose the additional assumption that $\mathcal{R}_{i}$ are decreasing with $i$, along with our immundominance hierarchy (19), to avoid degeneracy,

[TABLE]

For each $k\in[1,n]$, an equilibrium $\mathcal{E}^{*}$ with persistent variant sets $\Omega_{y}=\Omega_{z}=[1,k]$, must have component $x^{*}=1/\mathcal{P}_{k}$ by equation (10) (or by Proposition 3), where:

[TABLE]

More generally we can define $\mathcal{P}_{\mathcal{J}}=1+\sum_{i\in\mathcal{J}}\mathcal{R}_{i}s_{i}$ for any subset $\mathcal{J}\subseteq[1,n]$. Then for each $k\in[0,n]$, the following equilibria are found:

[TABLE]

Let $k\in[1,n+1]$ be maximal such that

[TABLE]

where the sum on the right vanishes if $k=1$. There are two cases depending on whether (for $k<n+1$)

[TABLE]

or not. If the inequality holds (for $k<n+1$), then there is a unique saturated equilibrium, $\mathcal{E}^{\dagger}_{k}$, with $\Omega_{z}=\Omega_{y}=[1,k]$ and if it does not or $k=n+1$, then there is a unique saturated equilibrium, $\mathcal{E}^{\ddagger}_{k}$, with $\Omega_{z}=[1,k-1],\Omega_{y}=[1,k]$. The inequalities (11) are strict in either case so the system is $\Omega_{yz}$ permanent and $y_{i},z_{j}\to 0$ for all $i\notin\Omega_{y},j\notin\Omega_{z}$, by Theorem 3.1.

In fact, much more can be said about the asymptotic behavior of solutions. In particular, the other conclusions of Theorem 3.1 imply that the components $y_{i},z_{j}$ where $i\notin\Omega_{y},j\notin\Omega_{z}$ converge to zero, the asymptotic means of the persistent variants converge to equilibria values, and the global attractor satisfies (14) which consists of $k$ or $k-1$ distinct planar Lotka-Volterra predator-prey differential equations for $(y_{i},z_{i})$ constrained by the relation $\sum_{i}\mathcal{R}_{i}y_{i}=\sum_{i}\mathcal{R}_{i}s_{i}$. This fact was utilized in a chemostat model by Wolkowicz wolkowicz1989successful to prove that a unique saturated equilibrium is globally asymptotically stable for $k=1,2$ (including when $\Omega_{z}=\left\{1,2\right\},\Omega_{y}=\left\{1,2,3\right\}$). These results are summarized in the following theorem.

Theorem 4.2 (wolkowicz1989successful )

If $\mathcal{R}_{1}>\mathcal{P}_{1}$, let $k$ be the largest integer in $[1,n]$ such that $\mathcal{R}_{k}>\mathcal{P}_{k}$, otherwise let $k=0$. If $\mathcal{R}_{k+1}\leq\mathcal{P}_{k}$, then for $1\leq i\leq k$, $y_{i},z_{i}$ are uniformly persistent (if $\mathcal{R}_{k+1}>\mathcal{P}_{k}$, $y_{k+1}$ also persists), and the other variants globally converge to zero, $x(t)\rightarrow 1/\mathcal{P}_{k}$ (if $\mathcal{R}_{k+1}>\mathcal{P}_{k}$, $x(t)\rightarrow 1/\mathcal{P}_{k+1}$ and $y_{k+1}(t)\rightarrow 1-\frac{\mathcal{P}_{k}}{\mathcal{R}_{k+1}}$) as $t\rightarrow\infty$. Additionally, the corresponding equilibria ($\mathcal{E}^{\dagger}_{k}$ or $\mathcal{E}^{\ddagger}_{k+1}$) are locally stable (globally asymptotically stable when $k=0,1,2$), asymptotic averages converge to equilibria, and the global attractor satisfies (13) and (14).

If the assumption of strictly decreasing reproduction numbers are relaxed, then the strain-specific system (24) can have multiple degenerate saturated equilibria. However the full hypercube network for $n$ epitopes containing $2^{n}$ virus strains (model (6)) allows us to relax this particular assumption on reproduction numbers, avoid any degeneracy, and explicitly identify two possible saturated “strain-specific equilibria”. Essentially, $\mathcal{E}^{\dagger}_{n}$ and $\mathcal{E}^{\ddagger}_{n+1}$ are the only strain-specific equilibria which can be saturated in the full hypercube network. The proposition below (proved in Appendix A.2) contains this detailed result on strain-specific equilibria.

Proposition 6

Consider system (6) on the full network with $n$ epitopes ($m=2^{n}$ virus strains) under the assumption of mutational fitness costs (7) and the viral strains, $y_{i}$, ordered so that $\Lambda_{i}=\left\{i\right\}$ for $i=1,\dots,n$ and $\Lambda_{n+1}=\emptyset$. Let $\mathcal{P}_{n}$ be defined by (26). Then an equilibrium $\mathcal{E}^{*}$ with a strain-specific subgraph, i.e. $\Omega_{y}\subseteq[1,n+1]$, is saturated if and only if one of the following holds:

• i.

$\mathcal{R}_{n+1}\leq\mathcal{P}_{n}$* and $\left(|\Lambda_{\ell}|-1\right)\mathcal{P}_{n}+\mathcal{R}_{\ell}\leq\sum\limits_{i\in\Lambda_{\ell}}\mathcal{R}_{i}\quad\forall\ell\in[n+2,2^{n}]$, in which case $\Omega_{y}=\Omega_{z}=[1,n]$.*

• ii.

$\mathcal{R}_{n+1}>\mathcal{P}_{n}$* and $\left(|\Lambda_{\ell}|-1\right)\mathcal{R}_{n+1}+\mathcal{R}_{\ell}\leq\sum\limits_{i\in\Lambda_{\ell}}\mathcal{R}_{i}\quad\forall\ell\in[n+2,2^{n}]$, in which case $\Omega_{y}=[1,n+1]$ and $\Omega_{z}=[1,n]$.*

Therefore by applying Theorem 3.1 to the above Proposition 6, we obtain conditions for the stability of either (i) $\mathcal{E}^{\dagger}_{n}$ (or (ii) $\mathcal{E}^{\ddagger}_{n+1}$), and the corresponding uniform persistence of $y_{1},\dots,y_{n}$ (and $y_{n+1}$ in case of (ii)). Furthermore these are the only potential stable equilibria with persistent strain set contained in the strain-specific subnetwork ($\Omega_{y}\subseteq[1,n+1]$) for system (6).

4.3 Dynamics for full network on $n=2$ epitopes

If two epitopes are concurrently targeted by two distinct specific immune responses, which escape pathway will the virus follow and what mutant strains persist? In the nested network, we assumed that the virus escaped the most immunodominant response. However, in general, both CTL pressure and virus fitness cost determine selective advantage of a resistant mutant. For a single epitope, an escape mutant $y_{2}$ invades the wild-type $y_{1}$, if its reproductive number is large enough given the CTL pressure on $y_{1}$. In particular, $y_{2}$ persists if $(f-s_{1})\mathcal{R}_{1}>1$, where $\mathcal{R}_{2}=f\mathcal{R}_{1}$ and $s_{1}=1/\mathcal{I}_{1}$ with $f$ the fitness proportion of the wild-type reproductive number $\mathcal{R}_{1}$ and $\mathcal{I}_{1}$ the immune response reproduction number (see pruss2008global or the case $n=1$ in Sections 4.1 or 4.2). For the general case of $n=2$ epitopes, although the situation is fundamentally more complex, we will sharply characterize the dynamics in this section. The results suggest that immunodominance may play a larger role than viral fitness in determining the structure of the persistent virus-immune network.

The full network for $n=2$ epitopes (shown in Fig. 1) consists of 2 CTL populations, $z_{1}$ and $z_{2}$, and $m=4$ virus strains, $y_{i}\ i=1,\dots,4$, each with an associated binary string describing their resistance profile; $y_{1}\sim 00$, $y_{2}\sim 10$, $y_{3}\sim 01$, $y_{4}\sim 11$. Recall that a 0 in the $j^{th}$ bit of the binary string signifies susceptible to $z_{j}$, whereas 1 signifies resistance; for example $y_{2}\ (10)$ is resistant to $z_{1}$ but susceptible to $z_{2}$. We assume that $z_{1}$ is strictly immunodominant, i.e. $\mathcal{I}_{1}>\mathcal{I}_{2}$ or $s_{1}<s_{2}$. Each epitope escape comes with a fitness cost as before, therefore the viral reproduction numbers satisfy $\mathcal{R}_{1}>\max(\mathcal{R}_{2},\mathcal{R}_{3})\geq\min(\mathcal{R}_{2},\mathcal{R}_{3})>\mathcal{R}_{4}$. We note that the fitness cost for resistance to $z_{1}$ may be greater than or equal to the fitness cost to $z_{2}$, in other words the fitness of mutant $y_{2}$ is less than or equal to $y_{3}$ ($\mathcal{R}_{2}\leq\mathcal{R}_{3}$). Conversely, it may be the case that resistance to the dominant immune response, $z_{1}$, comes at less cost than resistance to the weaker response ($\mathcal{R}_{2}>\mathcal{R}_{3}$). Our underlying assumptions are summarized below:

[TABLE]

For clarity, we write the 7 equations in model (6) for this case $n=2$ with the chosen index notation:

[TABLE]

The dynamics are rigorously characterized in the Theorems 4.3 and 4.4 stated below. The theorems together present a sharp polychotomy which delineates the stability of nine potentially “strictly saturated” distinct equilibria, along with a degenerate case where a continuum of equilibria exists, and the corresponding uniform persistence of variants in each case. First, we derive these nine potential equilibria with their corresponding component values. In this way, we can capture all of the possible stable, non-negative equilibria and avoid listing equilibria that are always unstable or have negative components. First, consider the equilibria types encountered in our analysis of the nested and strain-specific subnetworks adapted to this case of $n=2$ epitopes. Note that any equilibrium, $\mathcal{E}^{*}$, with $z_{2}^{*}>0$ and $z_{1}^{*}=0$ or $y_{3}^{*}>0$ and $y_{2}^{*}=0$ will either be not saturated or have negative components when $s_{1}<s_{2}$. Then there are four equilibria with one or less persistent immune response, ranging from the infection-free equilibrium $\bar{\mathcal{E}}_{0}$ to the “escaped single immune” equilibrium $\widetilde{\mathcal{E}}_{1}$. Also, there are the nested and “strain-specific” type equilibria $\overline{\mathcal{E}}_{2},\widetilde{\mathcal{E}}_{2}$ and $\mathcal{E}^{\dagger}_{2},\mathcal{E}^{\ddagger}_{2}$ (with 2 or 3 virus strains). Next notice that Propositions 3 and 5 restrict the possibility of having an equilibrium with all four viral strains persistent to a degenerate case when $\mathcal{R}_{1}-\mathcal{R}_{2}-\mathcal{R}_{3}+\mathcal{R}_{4}=0$.

We are left to check a new equilibrium type $\widehat{\mathcal{E}}_{2}$ where $y_{1}$, $y_{2}$, $y_{3}$ coexist, and by Proposition 3, its component $\widehat{x}$ can be derived as follows:

[TABLE]

where $\mathcal{R}=\mathcal{R}_{2}+\mathcal{R}_{3}-\mathcal{R}_{1}$. The other components of $\widehat{\mathcal{E}}_{2}$ can be derived utilizing equilibrium equations (10). Therefore the distinct regimes for the 10 relevant equilibria forms are determined by the values of the viral and immune reproductive numbers, along with the following quantities:

[TABLE]

The following equilibria have two immune responses present:

[TABLE]

The four equilibria with one or zero immune responses are:

[TABLE]

The theorems classifying dynamics for model (30) with strict immunodominance, $s_{1}<s_{2}$, are as follows (proofs are in Appendix A.3):

Theorem 4.3 (Stability of equilibria with one or zero immune response)

Consider the model with two epitopes (30) under the assumptions (29) and suppose positive initial conditions, i.e. $x(0),y_{i}(0),z_{j}(0)>0$ for all $i=1,\dots 4,j=1,2$. Then the following results hold:

• i.

If $\mathcal{R}_{1}\leq 1$, then $\bar{\mathcal{E}}_{0}$ is GAS (globally asymptotically stable).

• ii.

If $1<\mathcal{R}_{1}\leq\mathcal{Q}_{1}$, then $\widetilde{\mathcal{E}}_{0}$ is GAS.

• iii.

If $\mathcal{R}_{1}>\mathcal{Q}_{1}\geq\mathcal{R}_{2}$, then $\bar{\mathcal{E}}_{1}$ is GAS.

• iv.

If $\mathcal{Q}_{2}\geq\mathcal{R}_{2}>\mathcal{Q}_{1}$, then $\widetilde{\mathcal{E}}_{1}$ is GAS.

Theorem 4.4 (Stability of equilibria with two immune responses present)

Consider the model with two epitopes (30) under the assumptions (29) and suppose positive initial conditions, i.e. $x(0),y_{i}(0),z_{j}(0)>0$ for all $i=1,\dots 4,j=1,2$. Then the stability of equilibria (33) are characterized as follows:

if $\mathcal{R}<\mathcal{Q}_{2}$ and $\mathcal{R}_{4}\leq\mathcal{Q}_{2}$, then $\bar{\mathcal{E}}_{2}$ is GAS. 2. 2.

if $\mathcal{R}<\mathcal{R}_{4}$, then $\widetilde{\mathcal{E}}_{2}$ is GAS. 3. 3.

if $\mathcal{Q}_{2}<\mathcal{R}<\mathcal{P}_{2}$ and $\mathcal{R}_{4}<\mathcal{R}$, then $\widehat{\mathcal{E}}_{2}$ is (locally) stable. Additionally, $\lim_{t\rightarrow\infty}x(t)=\widehat{x}=\frac{1}{\mathcal{R}}$, $\lim_{t\rightarrow\infty}y_{4}(t)=0$, and

[TABLE]

Furthermore $y_{i},z_{j},\ \ i=1,2,3,j=1,2$ are uniformly persistent. 4. 4.

if $\mathcal{P}_{2}<\mathcal{R}$ and $\mathcal{R}_{4}\leq\mathcal{P}_{2}$, then $\mathcal{E}_{2}^{\dagger}$ is GAS. 5. 5.

if $\mathcal{Q}_{2}<\mathcal{P}_{2}<\mathcal{R}_{4}<\mathcal{R}$, then $\mathcal{E}_{2}^{\ddagger}$ is GAS. 6. 6.

If $\mathcal{R}=\mathcal{R}_{4}$, then there is a continuum of saturated equilibria which forms a line connecting $y_{4}=0$ and $y_{3}=0$ boundaries at $\widehat{y}_{3}$ and $\widetilde{y}_{4}$, respectively, in the $(y_{3},y_{4})$ plane.

We make the following observations concerning the above theorems. First, note that the inequalities in the hypotheses of Theorems 4.3 and 4.4 cover all possible parameter combinations under the conditions (29). Indeed, observe that $\mathcal{P}_{2}\leq\mathcal{Q}_{2}\Leftrightarrow\mathcal{R}\leq 0$ and of course $\mathcal{Q}_{2},\mathcal{P}_{2}>0$, therefore the case $\mathcal{P}_{2}\leq\mathcal{Q}_{2}$ falls under case 1 or case 2 of the theorem. Additionally, $\mathcal{R}>\mathcal{Q}_{2}\Rightarrow\mathcal{R}_{2}>\mathcal{Q}_{2}$, separating this case from the cases considered in Theorem 4.3. Note also that any equilibrium, $\mathcal{E}^{*}$, with $z_{2}^{*}>0$ and $z_{1}^{*}=0$ or $y_{3}^{*}>0$ and $y_{2}^{*}=0$ will not be saturated when $s_{1}<s_{2}$.

For the case $\mathcal{R}_{4}=\mathcal{R}$ and $\mathcal{Q}_{2}<\mathcal{R}<\mathcal{P}_{2}$ , both equilibria $\widetilde{\mathcal{E}}_{2}$ and $\widehat{\mathcal{E}}_{2}$ are saturated, but the inequalities (11) are not strict. In this case, there are a continuum of saturated (neutrally stable) non-negative equilibria with $x(t)\rightarrow\frac{1}{\mathcal{R}}$ and $y_{4}(t)\rightarrow 1-\frac{\mathcal{Q}_{2}}{\mathcal{R}}-y_{3}^{*}$, where $0\leq y_{3}^{*}\leq 1-\frac{\mathcal{Q}_{2}}{\mathcal{R}}$. The endpoints of the line of equilibria in the $(y_{3},y_{4})$ plane correspond to $\widetilde{\mathcal{E}}_{2}$ and $\widehat{\mathcal{E}}_{2}$, where $\widetilde{y}_{3}=0$ and $\widehat{y}_{4}=0$, respectively. The bifurcation as $\mathcal{R}_{4}$ increases above $\mathcal{R}$, where Case 3 transitions to Case 6 and then to Case 2 of Theorem 4.4, is illustrated in Figure A1 in Appendix A.3. When $\mathcal{R}_{4}\neq\mathcal{R}$, it is not possible for all $y_{i}$ to persist, as the dynamics will fall into one of the other cases in Theorems 4.3 and 4.4. These results illustrate Proposition 5, which precludes the possibility of equilibria with all components $y_{i}$ positive except in the case $\mathcal{R}_{4}\neq\mathcal{R}$ since the viral strains form a 4-cycle in the associated hypercube graph.

An interesting finding is that strain $y_{2}\ (10)$ is present in all equilibria (and uniformly persistent) with mutation, even if it has a much higher fitness cost than $y_{3}\ (01)$, i.e. $\mathcal{R}_{3}>>\mathcal{R}_{2}$. Thus, if escape occurs in the two epitope setting, the viral quasispecies always includes the mutant $y_{2}$ with resistance to the immunodominant epitope, $z_{1}$. We can characterize the (linearized) invasion rate of an escape mutant, say $y_{2}$, from the single epitope, $z_{1}$, by calculating $\frac{\dot{y}_{2}}{y_{2}}$ at equilibrium $\overline{\mathcal{E}}_{1}$, to obtain $\lambda_{2}=\gamma_{2}\left(\frac{\mathcal{R}_{2}}{\mathcal{Q}_{1}}-1\right)=\frac{\mathcal{R}_{2}}{1+\mathcal{R}_{1}s_{1}}-1$. The mutant $y_{2}$ then converges to $y_{2}^{*}=1-\frac{1+\mathcal{R}_{1}s_{1}}{\mathcal{R}_{2}}$, in equilibrium $\widetilde{\mathcal{E}}_{1}$. The comparable quantities for mutant $y_{3}$ in escaping the single epitope $z_{2}$ (in the absence of $z_{1}$) are $\lambda_{3}=\gamma_{3}\left(\frac{\mathcal{R}_{3}}{1+\mathcal{R}_{1}s_{2}}-1\right)$ and $y_{3}^{*}=1-\frac{1+\mathcal{R}_{1}s_{2}}{\mathcal{R}_{3}}$. Even when $y_{3}$ (escaping $z_{2}$) has larger invasion characteristics in the single epitope escape setting ($y_{3}^{*}>y_{2}^{*}$ and $\lambda_{3}>\lambda_{2}$), $y_{2}$ will still persist and may actually exclude $y_{3}$ in the two epitope model. Figure 2,2 depict simulations of this scenario of $y_{2}$ excluding $y_{3}$ despite the larger selective advantage of $y_{3}$ in the single epitope setting. Figure 2 also shows further simulations of (30), which are all consistent with assertions of Theorem 4.3 and Theorem 4.4.

There are two special cases of model (30) where analysis further suggests the superior role of immunodominance over viral reproductive fitness. First, if we relax the assumption of strict immunodominance hierarchy, i.e. we allow that $s_{1}=s_{2}$, then a “strictly saturated” non-negative equilibrium $\mathcal{E}^{*}$ has the following properties: $z_{1}^{*}>0\Leftrightarrow z_{2}^{*}>0$ and $y_{2}^{*}>0\Leftrightarrow y_{3}^{*}>0$. The calculations for this case are presented in the Appendix A.4. Essentially, non-trivial strictly saturated equilibria are of the form $\widehat{\mathcal{E}}_{2}$, $\mathcal{E}_{2}^{\dagger}$ or $\mathcal{E}_{2}^{\ddagger}$ with the corresponding parameter regimes as defined in Table 1. Therefore, in the case of equal immunodominance ($s_{1}=s_{2}$), no matter the fitness of the distinct mutant strains $y_{2}$ and $y_{3}$, both strains can only persist together. Next, consider the case where the viral fitness cost to each epitope is equal, say the cost is $c=1-f$, and the strict immunodominance holds, $s_{1}>s_{2}$. Then $\mathcal{R}_{2}=\mathcal{R}_{3}=f\mathcal{R}_{1}$ and $\mathcal{R}_{4}=f^{2}\mathcal{R}_{1}$. This scenario will be analyzed in further generality for $n$ epitopes in the next section. The result for this special case of model (30) is that one of the nested equilibria, $\overline{\mathcal{E}}_{i}$ or $\widetilde{\mathcal{E}}_{i}$ ($i=0,1,2$) will be globally asymptotically stable. Thus, in contrast to the previous special case ($s_{1}=s_{2}$) where viral fitness did not determine persistence, here (in this special case of equal viral fitness costs) we find that persistence is determined by immunodominance.

The above results indicate that the immunodominance hierarchy is the most important factor in directing epitope escape, more so than viral fitness cost. In addition, the persistent variants depend upon reproductive numbers and quantities defined in (32), implying that the escape pathway depends upon the entire interacting system, not just the parameters associated with the single epitopes and corresponding resistant viral mutant strains. Both the relative importance of immunodominance and the multi-epitope dependence align with findings in an in vivo study of HIV patients liu2013vertical . From a broader ecological point of view, our results suggest top-down control of food webs, where top predators have more influence than intermediate species, along with the interconnectedness of complex ecological networks.

4.4 Uniform fitness costs for escape on full $n$-epitope network

We consider the full network for $n$ epitopes, in the special case where each viral mutation of an epitope incurs a uniform fitness cost, $c=1-f$, where $f\in(0,1)$ is the ratio between reproduction number of mutant and descendent strain. In particular, indexing the wild-strain here as $y_{1}=y_{w}$ with fitness $\mathcal{R}_{1}$, then we make the following assumption on the $2^{n}$ viral strains in the full network:

[TABLE]

where $d(y_{i},y_{1})$ is the Hamming distance between the associated binary sequences of the viral strains as defined earlier. Note that since the wild-strain $y_{1}$ is susceptible at all epitopes ($y_{1}\sim 0\cdots 0$), a viral strain $y_{i}$ with $d(y_{i},y_{1})=p,\ p\in[0,n]$ has mutated $p$ epitopes and thus has a (susceptible) epitope set of cardinality $n-p$, i.e. $|\Lambda_{i}|=n-p$. Then, we can write the model (6) as follows:

[TABLE]

As before, we assume the immunodominance hierarcy (19), $s_{1}\leq\cdots\leq s_{n}$. We denote the viral strains associated with “nested” (sequential) escape, in addition to the wild strain $y_{1}$, as $y_{2},\dots,y_{n+1}$. This means that for $i=1,\dots,n$, the epitope set of $y_{i}$ is $\Lambda_{i}=\left\{i,\dots,n\right\}$ (having escaped immune responses $z_{1},\dots,z_{i-1}$) and $\Lambda_{n+1}=\emptyset$ so $y_{2}\thicksim(1,0,\cdots,0),y_{3}=(1,1,0,\cdots,0),\cdots,y_{n+1}=(1,1,\cdots,1)$.

Our main result of this section is, informally, that viral escape from $n$ epitopes follows a “nested pattern” when the immunodominance hierarchy (19) is strict. In other words, one of the nested equilibria, $\widetilde{\mathcal{E}}_{k+1}$ or $\overline{\mathcal{E}}_{k}$ (introduced in Section 4.1) is stable. Roughly speaking, if viral fitness costs of resistance are the same independent of CTL type, it is better for a mutant strain to escape the current dominant CTL type than a sub-dominant one. In the Appendix A.5, we show that two other classes of equilibria, namely “strain-specific” and “one-mutation” equilibria, are unstable in this case of equal fitness costs, even when inequalities (19) are not strict. Now the main theorem is stated below.

Theorem 4.5

Consider model (36), the full network on $n$ epitopes ($m=2^{n}$) with equal fitness costs (35) and strict immunodominance hierarchy (strict inequalities in (19)). Suppose $y_{i}$, $i\in[1,n+1]$, is indexed so that $\Lambda_{i}=\left\{i,\dots,n\right\}$ for $i=1,\dots,n$, $\Lambda_{n+1}=\emptyset$. Then $y_{i}(t)\rightarrow 0$ as $t\rightarrow\infty$ for all $i\in[n+2,2^{n}]$, and Theorem 4.1 holds in system (36). In particular, $y_{i},z_{i}$ are uniformly persistent for $1\leq i\leq k\leq n$ when $\mathcal{R}_{k}>\mathcal{Q}_{k}$ (and $y_{k+1}$ is also persistent if $\mathcal{R}_{k+1}>\mathcal{Q}_{k}$).

Proof

If $\mathcal{R}_{1}>\mathcal{Q}_{1}$, let $k$ be the largest integer in $[1,n]$ such that $f^{k-1}\mathcal{R}_{1}=\mathcal{R}_{k}>\mathcal{Q}_{k}$, otherwise let $k=0$. We will apply Theorem 3.1. It suffices to check the invasion rate for $y_{\ell}$ in (11) for $\ell\in[n+2,2^{n}]$ since Theorem 4.1 establishes the result in the perfectly nested submodel consisting of $y_{1},\dots,y_{n+1}$ and $z_{1},\dots,z_{n}$. First suppose that $f^{k}\mathcal{R}_{1}=\mathcal{R}_{k+1}\leq\mathcal{Q}_{k}$. Let $\ell\in[n+2,2^{n}]$ and it suffices to consider $\Lambda_{\ell}\cap[1,k]$ since the calculations will be considered at equilibrium $\overline{\mathcal{E}}_{k}$ where $z_{i}^{*}=0$ for all $i\geq k+1$. Suppose that $d(y_{\ell},y_{1})=p$ where $p\in[1,n-1]$. Then $\mathcal{R}_{\ell}=f^{p}\mathcal{R}_{1}$. Consider the invasion rate for $y_{\ell}$ at the equilibrium $\overline{\mathcal{E}}_{k}$:

[TABLE]

If $p>k$, then clearly

[TABLE]

by the definition of $k$. If $p\leq k$, $|\Lambda_{\ell}|=n-p\geq n-k$ and note that $\Lambda_{\ell}\neq\Lambda_{k+1}=\left\{k+1,\dots,n\right\}$ (where $\Lambda_{k+1}$ is the epitope set of $y_{k+1}$). Thus, $[1,k]\cap\Lambda_{\ell}\neq\emptyset$. By (23), for $2\leq i\leq k-1$,

[TABLE]

Similarly $\overline{z}_{k}-\overline{z}_{k-1}<0$, and thus at equilibrium $\overline{\mathcal{E}}_{k}$, we find

[TABLE]

If $p=k$, then

[TABLE]

If $p\leq k-1$, then $[1,k-1]\cap\Lambda_{\ell}\neq\emptyset$ so

[TABLE]

Therefore, the equilibrium $\overline{\mathcal{E}}_{k}$ is saturated with inequalities (11) strictly holding and $\overline{\mathcal{E}}_{k}$ is unique in its positivity class. Thus Theorem 3.1 can be applied to obtained the conclusions of Theorem 4.1. If $f^{k}\mathcal{R}_{1}=\mathcal{R}_{k+1}>\mathcal{Q}_{k}$, then a similar argument works with $x^{*}=\frac{1}{\mathcal{R}_{k+1}}$ instead of $x^{*}=\frac{1}{\mathcal{Q}_{k}}$, which shows that $\widetilde{\mathcal{E}}_{k+1}$ is stable in that case.

In Figure 3, we illustrate Theorem 4.5 by numerical solution of the model (36) in the case of $n=3$ epitopes. The simulations show that after some transient dynamics, only the viral strains associated with the nested network, $y_{1},y_{2},y_{3},y_{4}$, persist. In other words, the full network of $n=3$ epitopes (displayed in Figure 1) converges to the perfectly nested subnetwork (displayed in Figure 1).

5 Discussion

In this paper, we analyzed a virus model consisting of target cells, multiple virus strains and several immune response populations. The interaction of virus and immune response is described by a network reflecting the avidity of each distinct immune response in recognizing each particular virus strain. We find some general conditions on stability and feasibility of equilibria, along with uniformly persistent virus and immune response variants by utilizing Lyapunov function techniques.

We specialize the model to consider the scenario where the immune response populations are $n$ different CTL lines, each specific to a particular epitope, and there are $2^{n}$ virus strains containing all possible combinations of (resistance conferring) epitope mutations. In this case, the viral strains in the virus-immune network can be translated to an $n$-dimensional hypercube graph representative of the potential pathways of immune escape by the virus. The number of uniformly persistent viral strains and CTL populations can be built up in an ordered fashion dependent on derived invasion thresholds for “strain-specific” and “perfectly nested” subgraphs. For the full network including $2^{n}$ viral strains, in certain cases we sharply characterize dynamics and stability of distinct equilibria. In particular, we find that for

• (i)

$n=2$ epitopes: the escape pathway (persistent viral strain set) always includes the mutant resistant to the immunodominant epitope, even when it suffers a relatively high fitness cost,

• (ii)

equal fitness costs for each of $n$ epitopes: the network of $2^{n}$ viral strains always converges to a perfectly nested subgraph (ordered according to immunodominance) with less than or equal to $n+1$ strains.

Our theoretical results indicate that a diverse viral “quasispecies” can be built through resistance mutations at multiple epitopes and the immunodominance hierarchy is the most important factor determining the escape pathway. Prior research on HIV data has found that CTL immune responses drive within-host HIV diversity pandit2014reliable and immunodominance hierarchy is the major factor shaping viral escape from multiple epitopes Barton ; liu2013vertical . Thus our model framework and analysis can shed light on the dynamic network of interacting virus and immune response variants, which may be important to understand for development of effective immunotherapies.

Future research can build upon the results presented here in several ways. First, since rapid evolution of HIV due to CTL pressure is motivation for our model, mutation between the different viral strains can be explicitly included in the model. Preliminary simulations from stochastic versions of the model show that qualitative dynamics are preserved under mutation. Rigorous global perturbation arguments in the deterministic setting may be attempted to show the effect of small mutation rates, however for the general case, this would rely upon the conjecture of global stability for equilibria. Thus, another important theoretical question is proving global stability when uniform persistence occurs. Finally, further extensions of our work can be applicable to coevolution of virus and cross-reactive antibodies during HIV infection luo2015competitive , along with other complex ecological networks.

Acknowledgement

We would like to thank the editor and anonymous reviewers for comments which have led to an improved manuscript. We also acknowledge support by the Simons Foundation Grant 355819.

Appendix A Appendix

A.1 Existence of saturated equilibrium

Proof (Proof of Proposition 1)

We find it somewhat easier to work with the unscaled system (1) so we introduce some new parameters into it. Let $\epsilon>0$ be so small that $b-\epsilon(m+n)>0$ and $0\leq\lambda\leq 1$ be a homotopy parameter. Our perturbed system is given by:

[TABLE]

We refer to the vector field on the right side as $G(W,\lambda,\epsilon)$, where $W=(X,Y,Z)\in\mathbb{R}^{m+n+1}_{+}$. Then $G(W,1,0)$ is the vector field given in equations (1).

Straightforward calculation establishes that

[TABLE]

for $d=\min\{c,\delta_{i},\mu_{j}\}$. Fix $p>b/d$ and let

[TABLE]

We claim that there are no equilibria of $G(W,\lambda,\epsilon)$ on the boundary of $\mathcal{U}$. Notice that any equilibrium $E=(X,Y,Z)\in\mathcal{U}$ satisfies $Y_{i}>0$ and $Z_{j}>0$ because of the perturbation terms and because $b>\epsilon(m+n)$. If $E$ belongs to the boundary of $\mathcal{U}$, then either $X=0$ or $X+\sum_{i}Y_{i}+\sum_{j}Z_{j}/q_{j}=p$. Each of these is easy to rule out. Suppose, for example the latter holds. Then the left side of the differential inequality above vanishes at $E$ so we have $0\leq b-pb$, a contradiction to $p>b/d$.

Now we can employ degree theory, as in Hofbauer and Sigmund hofbauer1998evolutionary . By Homotopy invariance of degree we have that $\hbox{deg}(G(\bullet,1,\epsilon),\mathcal{U})=\hbox{deg}(G(\bullet,0,\epsilon),\mathcal{U})$ and the latter is easy to compute since it is linear and has a unique equilibrium $\tilde{E}\in\mathcal{U}$

[TABLE]

As the degree is nonzero, $G(\bullet,1,\epsilon)$ has at least one equilibrium in the interior of $\mathcal{U}$ for each small $\epsilon$. Now, the argument for the existence of a saturated equilibrium of $G(\bullet,1,0)$ follows as in Theorem 13.4.1 in Hofbauer & Sigmund’s text hofbauer1998evolutionary by taking limits as $\epsilon\to 0$.

A.2 Strain-specific subnetwork

Proof (Proof of Proposition 6)

Consider system (6) under the prescribed assumptions of Proposition 6. We claim that any equilibrium $\mathcal{E}^{*}$ with $\Omega_{z}\subseteq[1,n]$ and $\Omega_{y}\subseteq[1,n+1]$ is saturated only if $\Omega_{z}=[1,n]$ and $\Omega_{y}\supseteq[1,n]$. In other words $\mathcal{E}^{\dagger}_{n}$ and $\mathcal{E}^{\ddagger}_{n+1}$ are the only strain-specific equilibria which can be saturated in the full hypercube network. Suppose by way of contradiction that we can choose $i,j\in[1,n]$ such that $i\in\Omega_{z}$ and $j\notin\Omega_{z}$. Then consider strain $y_{\ell}$ such that $\Lambda_{\ell}=\left\{i,j\right\}$, i.e. $y_{\ell}\leftrightarrow y_{i}$ and $y_{\ell}\leftrightarrow y_{j}$. At equilibrium $\mathcal{E}^{*}$, its invasion rate (on of the “saturated” conditions (11)) is given by $\frac{\dot{y}_{\ell}}{\gamma_{\ell}y_{\ell}}=(\mathcal{R}_{\ell}-\mathcal{R}_{i})/\mathcal{P}$, where $\mathcal{P}=\mathcal{P}_{\mathcal{J}}$ or $\mathcal{R}_{n+1}$ with $\mathcal{J}\subset[1,n]$. Since $\mathcal{R}_{\ell}>\mathcal{R}_{i}$, equilibrium $\mathcal{E}^{*}$ can not be saturated. The conditions for equilibria $\mathcal{E}^{\dagger}_{n}$ to be saturated is as follows in the general model:

[TABLE]

An analogous statement holds for $\mathcal{E}^{\ddagger}_{n+1}$ in case (ii).

A.3 Two-epitope model: dynamics when $s_{1}<s_{2}$

Proof (Proof of Theorems 4.3 and 4.4)

We apply Theorem 3.1 for each equilibrium and case. Note that the Lyapunov function $W$ is the same as in the proof of Theorem 3.1 except that $s_{j}$ replaces $\rho_{j}$. First in Theorem 4.3:

Case i.:

[TABLE]

Case ii.:

[TABLE]

Case iii.:

[TABLE]

Case iv.:

[TABLE]

In case iii, on the attracting invariant set $\mathcal{L}$ (same $\mathcal{L}$ as in proof of Theorem 3.1), we have $x=1/\mathcal{Q}_{1}$, $y_{3}=y_{4}=0$ and $z_{2}=0$. If $\mathcal{R}_{1}<\mathcal{Q}_{1}$, then $\dot{W}=0\Rightarrow y_{2}=0$, otherwise $\mathcal{R}_{1}=\mathcal{Q}_{1}\Rightarrow\dot{y}_{2}=0$. Either way the equation $\dot{x}=0$ implies that $y_{1}=y_{1}^{*}$ and $y_{2}=0$. Thus $\dot{y}_{1}=0$, which implies $z_{1}=z_{1}^{*}$. So the invariant omega limit set, which is nonempty by compactness of orbits, is $\mathcal{L}=\{\bar{\mathcal{E}}_{1}\}$ and we get convergence.

In case iv, i.e. $\mathcal{Q}_{1}<\mathcal{R}_{2}\leq\mathcal{Q}_{2}$, then we prove $\widetilde{\mathcal{E}}_{1}$ is GAS. Here, $\dot{W}=0$ iff $x=x^{*}=\frac{1}{\mathcal{R}_{1}}$, $y_{3}=y_{4}=0$, and $z_{2}=0$. (Note if $\mathcal{Q}_{2}>\mathcal{R}_{2}$, $z_{2}=0$ is immediate. If not, we can still reason that the asymptotic average of $z_{2}$ must converge to $z_{2}^{*}=0$, hence $z_{2}=0$.) In addition by Theorem 3.1, we can obtain $y_{2}=y_{2}^{*}=1-\frac{\mathcal{Q}_{1}}{\mathcal{R}_{1}}$. The last relation combined with $x=x^{*}$ implies that $y_{1}=y_{1}^{*}=s_{1}$. Then $\dot{y}_{1}=0$ implies $z_{1}=z_{1}^{*}$. Thus $\mathcal{L}$ consists solely of the equilibrium $\widetilde{\mathcal{E}}_{1}$.

Next in Theorem 4.4:

Case 1: $(x^{*},y^{*},z^{*})=\bar{\mathcal{E}}_{2}$:

[TABLE]

Notice that $\dot{W}\leq 0$ when $\mathcal{R}\leq\mathcal{Q}_{2}$ and $\mathcal{R}_{4}\leq\mathcal{Q}_{2}$. Also, the equilibrium $\bar{\mathcal{E}}_{2}$ is non-negative when $\mathcal{R}_{2}>\mathcal{Q}_{2}$. If $\mathcal{R}<\mathcal{Q}_{2}$ and $\mathcal{R}_{4}<\mathcal{Q}_{2}$, applying Theorem 3.1, all inequalities (11) are strict and only two strains, $y_{1}$ and $y_{2}$ have non-empty epitope sets, and therefore $\bar{\mathcal{E}}_{2}$ is globally asymptotically stable. If $\mathcal{R}_{4}=\mathcal{Q}_{2}$, then the differential equations in (14) hold along with $\dot{y}_{4}=0$ and $\sum_{i=1,2,4}\mathcal{R}_{i}y_{i}=\mathcal{Q}_{2}-1$. Taking asymptotic averages as in the proof of Theorem 3.1 in browneNest, we obtain that indeed $y_{4}=0$, and we similarly obtain that $\bar{\mathcal{E}}_{2}$ is globally asymptotically stable.

Case 2: $(x^{*},y^{*},z^{*})=\widetilde{\mathcal{E}}_{2}$:

[TABLE]

Notice that $\dot{W}\leq 0$ when $\mathcal{R}<\mathcal{R}_{4}$. Also, the equilibrium $\widetilde{\mathcal{E}}_{2}$ is non-negative when $\mathcal{R}_{4}>\mathcal{Q}_{2}$. Applying Theorem 3.1, we obtain that $\widetilde{\mathcal{E}}_{2}$ is globally asymptotically stable.

Case 3: $(x^{*},y^{*},z^{*})=\widehat{\mathcal{E}}_{2}$:

[TABLE]

Notice that $\dot{W}\leq 0$ when $\mathcal{R}_{4}\leq\mathcal{R}$. Also, the equilibrium $\widehat{\mathcal{E}}_{2}$ is non-negative when $\mathcal{Q}_{2}<\mathcal{R}<\mathcal{P}_{2}$. The result is a direct consequence of Theorem 3.1.

Case 4: $(x^{*},y^{*},z^{*})=\mathcal{E}^{\dagger}_{2}$:

[TABLE]

Notice that $\dot{W}\leq 0$ when $\mathcal{R}\geq\mathcal{P}_{2}$ and $\mathcal{R}_{4}\leq\mathcal{P}_{2}$. The equilibrium $\mathcal{E}^{\dagger}_{2}$ is non-negative when $\min(\mathcal{R}_{2},\mathcal{R}_{3})>\mathcal{P}_{2}$. Global stability follows from Theorem 3.1.

Case 5: $(x^{*},y^{*},z^{*})=\mathcal{E}^{\ddagger}_{2}$:

[TABLE]

Notice that $\dot{W}\leq 0$ when $\mathcal{R}_{4}\leq\mathcal{R}$. Also, the equilibrium $\mathcal{E}^{\dagger}_{2}$ is non-negative when $\mathcal{R}_{4}>\mathcal{P}_{2}$. Global stability follows from Theorem 3.1.

Case 6: $(x^{*},y^{*},z^{*})=\mathcal{E}(y_{4}^{*})$ ($\mathcal{R}_{4}=\mathcal{R}$):

[TABLE]

Notice that $\dot{W}\leq 0$ and the equilibrium $\mathcal{E}(y_{4}^{*})$ is non-negative when $\mathcal{Q}_{2}\leq\mathcal{R}<\mathcal{P}_{2}$. Thus $\mathcal{E}(y_{4}^{*})$ is saturated, though not strictly. Local stability follows from Theorem 3.1.

In Figure A1, we illustrate the bifurcation at $\mathcal{R}_{4}=\mathcal{R}$.

A.4 Two-epitope model: dynamics when $s_{1}=s_{2}$

We analyze the feasible stable equilibria with immune response for the case of equal immunodominance of $z_{1}$ and $z_{2}$ ($s_{1}=s_{2}=s$) in model (30). First, consider equilibria with one immune response present. Since $s_{1}=s_{2}$, without loss of generality, we can take $z_{1}^{*}>0$. From the equilibrium equations, $y_{1}^{*},y_{3}^{*}>0\Leftrightarrow\mathcal{R}_{1}=\mathcal{R}_{3}$, which is not possible since $\mathcal{R}_{1}>\mathcal{R}_{3}$. Also, an equilibrium with $y_{3}^{*}>0$ (and $z_{2}^{*}=0$) will not be saturated since $\mathcal{R}_{1}>\mathcal{R}_{3}$, therefore take $y_{3}^{*}=0$. Similarly we can take $y_{4}^{*}=0$. Thus, consider equilibria $\widetilde{\mathcal{E}}_{1}$ and $\overline{\mathcal{E}}_{1}$. The equilibrium $\widetilde{\mathcal{E}}_{1}$ can not be stable since $s_{1}=s_{2}\Rightarrow\mathcal{Q}_{1}=\mathcal{Q}_{2}$, which does not permit conditions for case iv of Theorem 4.3 to be satisfied. The equilibrium $\overline{\mathcal{E}}_{1}$ will not be “strictly saturated” under conditions for case iii of Theorem 4.3 because there are a continuum of equilibria of the form $\overline{\mathcal{E}}_{1}(z_{2}^{*})$ where $\bar{z}_{1}=\frac{\mathcal{R}_{1}}{\mathcal{Q}_{1}}-1-z_{2}^{*},\bar{z}_{2}=z_{2}^{*}$.

Now consider the possibility of equilibria with $z_{1}^{*},z_{2}^{*}>0$. Observe that at least two of the viral strain components $y_{1}^{*},y_{2}^{*},y_{3}^{*}$ are positive and $y_{2}^{*}=y_{3}^{*}$ from the $\dot{z}_{1},\dot{z}_{2}$ equations. Therefore $y_{2}^{*},y_{3}^{*}>0$ and for the case $s_{1}=s_{2}$, non-trivial strictly saturated equilibria are of the form $\widehat{\mathcal{E}}_{2}$, $\mathcal{E}_{2}^{\dagger}$ or $\mathcal{E}_{2}^{\ddagger}$ with the corresponding parameter regimes as defined in Table 1.

A.5 Instability of “one-mutation” and “strain-specific” equilibria for (36)

We assume equal viral fitness costs for mutation from each of $n$ epitopes, yielding system (36), as described in Section 4.4. Consider the set $\mathcal{S}_{1}$ of $n$ viral strains which have exactly one mutation, i.e. $y_{i}\in\mathcal{S}_{1}\Rightarrow d(y_{i},y_{1})=1$. For clarity here, we label these strains as $\mathcal{S}_{1}=\left\{y_{i}^{1}\ |\ i=1,\dots,n\right\}$ where $y_{i}^{1}$ has escaped $z_{i}$ but is susceptible all other immune responses. Note that $y_{1}^{1}$ is strain $y_{2}$ with our “nested” indexing introduced in Section 4.1. The subsystem only containing these viral strains looks as follows:

[TABLE]

By Proposition 3 and (10), a positive equilibrium $\mathcal{E}^{*}_{1}=(x^{*},y^{*},z^{*})$ of system (37) satisfies

[TABLE]

where $I_{n}$ is the $n\times n$ identity matrix. Here we find that:

[TABLE]

Assuming our hierarchy $s_{i}\leq s_{i+1}$, then $y_{i}^{1*}>0$ if $s_{1}>\sum_{i>1}(s_{n}-s_{i})$ and $z_{i}^{*}>0$ if $f\mathcal{R}_{1}\left(n-1-\sum_{i}s_{i}\right)>n-1$. If these conditions are satisfied, then the equilibrium $\mathcal{E}^{*}_{1}$ is saturated in the subsystem (37) where $y_{i}\in\mathcal{S}_{1}$. However, if we consider a larger network of viral strains, then equilibrium $\mathcal{E}^{*}_{1}$ is always unstable in this case with equal fitness costs for mutation. Indeed, consider the wild strain $y_{1}\sim 0\cdots 0$ and strain $y_{3}\sim 110\cdots 0$, which correspond to backward and forward mutations from the strain $y_{1}^{1}\in\mathcal{S}_{1}$. We calculate their invasion rates at equilibrium $\mathcal{E}^{*}_{1}$. Note that $y_{1}$ has reproduction number $\mathcal{R}_{1}$ and $y_{3}$ has reproduction number $\mathcal{R}_{3}=f^{2}\mathcal{R}_{1}$. In order for $\mathcal{E}^{*}_{1}$ to be saturated (stable), we find that:

[TABLE]

However, $(n-1)(f-1)^{2}>0$ for all $f\neq 1,n>1$, giving a contradiction. Thus $\mathcal{E}_{1}$ can only be stable when restricted to the subsystem (37) consisting of strains in $\mathcal{S}_{1}$, and becomes unstable in the larger network in this case.

Similarly, we can show that the “strain-specific” equilibria, $\mathcal{E}^{\dagger}_{n}$ and $\mathcal{E}^{\ddagger}_{n+1}$ (introduced in Section (4.2)), are always unstable in this case of uniform fitness costs. Indeed, for clarity here, denote $y^{o}_{i}$ as the viral strain with epitope set $\Lambda^{o}_{i}=\left\{i\right\}$, and $y_{n+1}$ as the strain which has completely escaped all $z_{i}$, i.e. $\Lambda_{n+1}=\emptyset$. Then $\mathcal{E}^{\dagger}_{n}$ (and $\mathcal{E}^{\ddagger}_{n+1}$) consist of equilibria where $y^{o*}_{i}>0$ (and $y_{n+1}^{*}>0$). First, consider the case that $\mathcal{R}_{n+1}=f^{n}\mathcal{R}_{1}\leq\mathcal{P}_{n}^{o}=1+\sum_{i}s_{i}\mathcal{R}^{o}_{i}=1+f^{n-1}\mathcal{R}_{1}\sum_{i}s_{i}$ and $f^{n-1}\mathcal{R}_{1}>\mathcal{P}_{n}^{o}$, so that $\mathcal{E}^{\dagger}_{n}$ is positive, but $\mathcal{E}^{\ddagger}_{n}$ is not positive. Then consider the invasion rate of viral strain $y_{n-1}\sim 1\cdots 100$ (with $\Lambda_{n-1}=\left\{n-1,n\right\}$ and reproduction number $\mathcal{R}_{n-1}=f^{n-2}\mathcal{R}_{1}$). Utilizing Proposition 6, if $\mathcal{E}^{\dagger}_{n}$ is stable, then

[TABLE]

which is clearly a contradiction since $\left(f^{2}-2f+1\right)=(f-1)^{2}>0$ for $f<1$. A similar argument applies to $\mathcal{E}^{\ddagger}_{n+1}$ in the case $f^{n}\mathcal{R}_{1}>\mathcal{P}_{n}^{o}$. Thus the “strain-specific” equilibria are always unstable for system (36).