Age-structured viral dynamics in a host with multiple compartments

1.   Introduction

Several mathematical models in the field of viral dynamics simply regard the host of an infection as a well-mixed compartment ^[1,2,3] and aim to theoretically determine the threshold, from an overall perspective for determining whether an infection is extinct or persistent. Under the assumption of spatial homogeneity, various functions of different tissues and organs may be overlooked. Given the complexity of the human body, it is essential to regard the host body as composed of multiple compartments and to explore the effect of connection between compartments on the infection ^[4,5,6].

Viral propagation in human immunodeficiency virus (HIV) infection occurs mainly in lymphoid tissues ^[5,7], which are composed of primary and secondary lymphoid organs. The primary lymphoid organs include the thymus and bone marrow, which are where lymphocytes are generated; the lymph nodes and the spleen constitute the secondary lymphoid organs, which maintain mature naive lymphocytes and initiate adaptive immune response ^[8]. Between these tissues, connections are constructed through afferent or efferent lymphatics (or both) ^[9,10,11], by which cluster or differentiation 4 (CD4)-positive or CD8-positive T lymphocytes and other cells can migrate to execute their functions. This gives rise to the need to explore how the connection of various tissues influences the viral dynamics and curative effect of treatment.

To analyze the spatiotemporal patterns of viral progression, the authors in ^[12] proposed a reaction-diffusion model by fitting the data to a radially symmetric propagation and numerically identified the initial propagation of the infection front followed by a state of stagnation due to an antiviral response. Infected CD4 T cells can actively migrate over different lymph nodes ^[13]. Accordingly, another scenario for studying the spatial heterogeneity of viral dynamics is using multi-compartment models to describe the occurrence of infection in various organs. In ^[14], a heterogeneous environment with various parameters for several grids was constructed, wherein the virus randomly spreading between the grids was modeled as local diffusion. The main result showed that the local dispersion of a virus can reduce the amplitude of viral oscillation. In a recent study ^[5], the spread of HIV infection throughout the lymphoid tissue was considered from a perspective of network structure. That study numerically determined that HIV infection may persist even under drug administration because a portion of lymphoid tissues fails to receive a sufficient amount of the dose. These studies clarify that considering the spatial heterogeneity is essential in studying viral dynamics.

Most within-host virus models assume transmission solely through virus-to-cell infection, although it is now appreciated that there are multiple infection routes resulting in heterogenous kinetics. The cell-to-cell transmission process of HIV infection has been explored in a compartment such as the lymph nodes and the brain, because CD4+ T lymphocytes are densely agammaegated and frequently interact in lymphoid tissues, and the virus can disseminate through direct transmission from infected cells to uninfected cells through the HIV virological synapse ^[15,16]. Among the studies that have considered the delayed intracellular reaction but without distinction on infection age, HIV infection was explored in ^[1] by applying ordinary and delayed differential models, and sustained oscillations of infection were revealed for certain culture parameters. In ^[17], the infection-age-dependent infectivity of individuals was considered and the dynamics of HIV transmission were studied in a population group. The age-dependent production rate of viral particles and death rate of productively infected cells were considered in ^[18], in which an in vivo age-structured model of HIV infection was constructed. In order to explore distinct intracellular and extracellular infection kinetics, it is necessary to extend previous virus models by incorporating both age-since-infection and a multi-compartment environment.

Global convergence in viral dynamics models is crucial for predicting the evolution of viral infection, and age-structure of infected cells tends to complicate matters. Our single-compartment model considered in this paper is related to that in ^[19], wherein a susceptible-infective (SI) epidemic model was incorporated with infection-age and an incidence function derived from the law of mass action. In that study, the system was addressed using an integrated semigroup to determine the persistence of disease; the results demonstrated global convergence to either disease-free or endemic equilibrium through the application of Lyapunov functionals. Notably, the probability of disease transmission was found to depend on the age of infection ^[19]. Another difficulty in the field arises because of the multigroup structure of the models. In the recent study ^[20], the authors investigated a multigroup (susceptible-infective-recovered) SIR model with age structure and elaborated on the convergence dynamics when the morality and removal rate were age-independent; because of the constant total host population, the model involved a constant boundary condition. Furthermore, in ^[21], an aged model in a network environment was considered, and both the age-dependent transmission rate and degree of the network connectivity were found to affect the spread of disease; on the basis of the constant total population, the authors demonstrated that global convergence to the endemic equilibrium occurs when it exists. In ^[22], a two-group model was adopted to investigate the globally asymptotic behavior of a SI epidemic model incorporating the age of infection. The Lyapunov functional employed in that study was skillfully constructed for the two-group structure. The age-dependent infectivity of infected cells renders it impossible to rewrite the viral infection model as the classic (without aged-structure) SIS or SIR epidemiological models or variant predator-prey models, and it increases the difficulty of exploring the global convergence dynamics, especially in deriving convergence to the endemic equilibrium state.

In this paper, we construct and analyze a virus model consisting of multiple compartments with distinct cell infection-age structured infectivity kernels describing infection spread through various organs. In section 2, we detail the model formulation, along with establishing boundedness and compactness of the generated semiflow. In sections 3 and 4, we derive the basic reproduction number, and show that it is a global threshold determining the viral persistence or extinction. Section 5 concerns extending analysis to non-strongly connected networks through a sequence of threshold values controlling the infection pattern in the whole system. Finally in Section 6, we numerically investigate the viral dynamics for different connection topologies, along with examples of distinct cell infectivity kernels based on data for cell-free and cell-to-cell infection modes.

2.   Model formulation and preliminaries

An infected host with $n$ compartments is described by

for $k\in\mathbf{N}_0: = \{1, 2, \cdots, n\}$, with the boundary condition

and the initial condition

where $L^1_+((0, +\infty), \mathbb{R})$ is the nonnegative cone of $L^1((0, +\infty), \mathbb{R})$. In the $k$th compartment, at time $t$, $T_k(t)$ denotes the concentration of healthy cells, and $i_k(t, a)$ is that of infected cells with age $a$ since infection. The parameters $\lambda_k$ and $d_k$ respectively represent the production and death rates of healthy cells, and $\alpha_k$ accounts for the rate of both virus-to-cell and cell-to-cell transmission. Here, we neglect the virus compartment due to its fast kinetics. The function $p_k(a)$ accounts for the net infectivity of infected cells from both infection routes, variant in the infection age, and $\delta_k(a)$ represents the death rate of infected cells. All these parameters describe the viral environment in the $k$th compartment. They may be different from compartment to compartment because of individual function. For instance, the infectivity kernel $p_k(a)$ may change according to primary transmission mode (cell-to-cell versus viral), or infection induced death rate, $\delta_k(a)$, may be different in distinct compartments based on predominant target cell type (e.g. T-cell in blood vs. macrophage in brain). The parameter $m_{jk}$, $j\neq k$, denotes the migration rate of the cells from the $k$th compartment to the $j$th compartment. Assuming no loss in the process of migration, it satisfies that $\sum^n_{j = 1}m_{jk} = 0$ and then $m_{kk} = -\sum^n_{j\neq k, j = 1}m_{jk}$. As mentioned in Section 1, the migration may be asymmetric between any two compartments. In addition, there is no isolated compartments to be considered, however this assumption is relaxed later in Section 5.

Here, we assume that the migration matrix $M = (m_{kj})_{n\times n}$ is irreducible and, for $k\in\mathbf{N}_0$,

$(A_1)$ $\lambda_k, d_k, \alpha_k > 0$;

$(A_2)$ $\delta_k(\cdot)\in \mathcal{C}([0, +\infty), \mathbb{R}_+)$, $d_k\leq \delta_k(a)\leq \delta_{k, {\rm max}} < \infty$, for $a\geq 0$;

$(A_3)$ $p_k(\cdot)\in \mathcal{C}([0, +\infty), \mathbb{R}_+)$, $0 < p_k(a)\leq p_{k, {\rm max}} < \infty$, $a\in I_0$, for some finite interval $I_0\subset[0, +\infty)$, and $p_k(a) = 0$ for $a\notin I_0$.

Denoting $T(t) = (T_1(t), \cdots, T_n(t))^{\tau}$ and $i(t, a) = (i_1(t, a), \cdots, i_n(t, a))^{\tau}$, where $(\cdot)^{\tau}$ means the transpose of a vector, the system (2.1) is equivalent to the vector form:

with the boundary condition

and the initial condition $T(0) = T_{0} > 0, \; i(0, \cdot) = i_{0}(\cdot)\in L^1_+((0, +\infty), \mathbb{R}^n)$, where $\Lambda = (\lambda_1, \cdots, \lambda_n)^{\tau}$, $D = {\rm diag}(d_1, \cdots, d_n)$, $\Gamma = {\rm diag}(\alpha_1, \cdots, \alpha_n)$, ${\rm diag}(T(t)) = {\rm diag}(T_1(t), \cdots, T_n(t))$, $\tilde{D}(a) = {\rm diag}(\delta_1(a), \cdots, \delta_n(a))$, $P(a) = {\rm diag}(p_1(a), \cdots, p_n(a))$ and

Following the approach proposed in ^[23], we reformulate (2.4) as a semilinear Cauchy problem. To this end, we first consider the extended state spaces

and equip the space $\mathcal{X}$ with the norm

for $\mathbf{u} = (u_1, \cdots, u_n, v_1, \cdots, v_n, w_1(\cdot), \cdots, w_n(\cdot))^{\tau}\in \mathcal{X}$. Define the linear operator $A:Dom(A)\subset \mathcal{X}\rightarrow \mathcal{X}$, $Dom(A) = \mathbb{R}^n\times\{0_{\mathbb{R}^n}\}\times W^{1, 1}((0, +\infty), \mathbb{R}^n)$, by

where $W^{1, 1}$ is a Sobolev space, and the operator $F:\mathcal{X}_0\rightarrow \mathcal{X}$ by

By denoting $\mathbf{u}(t) = (T(t), 0_{\mathbb{R}^n}, i(t))^{\tau}$, we regard system (2.4) as the abstract Cauchy problem:

and assert, by [23,Theorems 2 and 3], that there exists a unique solution semiflow $\Psi(t):\mathcal{X}_{0+}\rightarrow \mathcal{X}_0$ related to system (2.7).

Next, by employing the boundary condition and the initial condition, we write (2.4) into the Volterra-type equation,

where $Q(t) = i(t, 0)$ and satisfies

Denote

First, we show a priority about the matrix $\Omega(a_0, a)$.

Lemma 2.1. For fixed $a_0\geq 0$, there exists a constant $\gamma > 0$ such that

for all $k, j\in\mathbf{N}_0$ and $a\geq a_0$, where $\underline{d} = \min_{k\in\mathbf{N}_0}\{d_k/2\} > 0$.

Proof. Denote $\tilde{d} = \max_{k\in\mathbf{N}_0}\{\delta_{k, {\rm max}}\} + \max_{k\in\mathbf{N}_0}\{-m_{kk}\}$ and $I$ as the $n\times n$ identity matrix. Then, for $a\geq a_0$, the matrix $-\int^a_{a_0}\tilde{D}(s)ds+\tilde{d}(a-a_0)I+(a-a_0)M$ is nonnegative and

In addition,

where the last inequality holds since the fact that the quasi-positive matrices $A = (a_{kj})\leq B = (b_{kj})$ with $0\leq a_{kj} = b_{kj}$ for all $k\neq j$ implies $\exp(A)\leq\exp(B)$. Since the matrix $-\underline{d}I+M$ is strictly diagonally dominant and each diagonal entry is negative, all its eigenvalues, say $\mu_j$, $j = 1, 2, \cdots, m$, $m\leq n$, have negative real parts [24,Theorem 6.1.10]. We write $-\underline{d}I+M$ in its Jordan canonical form, $J$, with $-\underline{d}I+M = PJP^{-1}$, $J = {\rm diag}(J_1, \cdots, J_{m})$ and each $J_j$ is a standard Jordan block related to the eigenvalue $\mu_j$. Obviously,

Since the real part of $\mu_j$ is negative,

as $a-a_0\rightarrow +\infty$. Thus, there exists a positive constant $\gamma$ such that $\exp((a-a_0)(-\underline{d}I+M))\leq \gamma \mathcal{I}$ for all $a\geq a_0$, where $\mathcal{I}$ is the matrix with a value of 1 for all entries, and the assertion is true.

Denote

and equip the space $\mathcal{Y}$ with the norm

Then the semiflow generated by (2.1) is point dissipative, as demonstrated in the following lemma.

Lemma 2.2. The solution semiflow generated by (2.1) is point dissipative. Explicitly, the subset of phase space

is positively invariant and attracts all nonnegative solutions, where $d_{{\rm min}} = \min_{k\in\mathbf{N}_0}\{d_k\}$.

Proof. It is easy to show that all solutions of (2.1) are nonnegative. Denote the total infected cells by $I(t) = \int^{\infty}_0i(t, a)da$ and the total cells by $N(t) = T(t)+I(t)$. Integrating the $i$-equation in (2.4) with respect to $a$ leads to

Obviously,

and then

By a comparison principle, we derive that

Since each solution of (2.1) remains nonnegative, as previously mentioned, the assertion holds true.

Deriving the uniform persistence and global properties of the solution dynamics is essential in studying viral infection and it is necessary to show that the semiflow $\Phi(t)$ generated by (2.1) is asymptotically smooth. Denote

By integrating along the characteristic lines and incorporating the boundary condition and the initial condition of the model, we write equation (2.1) in the following equivalent Volterra integral equation,

Few studies have investigated asymptotically smooth semiflows in a system incorporating both the age structure and the migration of a population. Here, we show the property in the semiflow $\Phi(t)$ by using Lemma 1 and a result reported in [25,Lemma 3.2.3].

Lemma 2.3. The semiflow $\Phi(t)$ generated by (2.10) ((2.1)) is asymptotically smooth.

Proof. It is easy to check that each function $\tilde{\mathcal{P}}_k(\cdot)$ is bounded, say with the bound $\bar{P}_k$, and Lipschitz continuous on $\mathbb{R}_+$, say with the Lipschitz constant $l_k$ (see ^[26]). Define $\Phi(t) = \Lambda_1(t)+\Lambda_2(t)$, where

for $\mathbf{x}\in \mathcal{Y}$. From (2.10) and Lemma 2.1,

Define $\Delta(r, t) = nr\gamma e^{-\underline{d}t}$. Then $\Delta(r, t)\rightarrow 0$ as $t\rightarrow +\infty$ and $\|\Lambda_1(t)\mathbf{x}\|\leq \Delta(r, t)$ for $\|\mathbf{x}\|\leq r$.

Let $B\subset \mathcal{Y}$ be a bounded subset such that $\Phi(t)B\subset B$. Choose $r_0 > 0$ such that $\|\mathbf{x}\|\leq r_0$ for all $\mathbf{x}\in B$. Note, from Lemma 2.2, that $\cup_{(T_0, i_0(\cdot))\in B}\{T(t)\}$ is bounded in $\mathbb{R}^n$ and then is precompact in $\mathbb{R}^n$. Hence, in order to show that $\Lambda_2(t)B$ is precompact, it suffices to verify that the set $\tilde{\Lambda}_2(t)B$ is precompact for

From ^[27], it is sufficient to verify the following conditions:

(ⅰ) $\lim_{h\rightarrow 0}\sum_{k\in\mathbf{N}_0}\int^{\infty}_0\{|i_k(t, a)-i_k(t, a+h)|\}da = 0$ uniformly for $i(t, a)\in\tilde{\Lambda}_2(t)B$.

(ⅱ) $\lim_{h\rightarrow +\infty}\sum_{k\in\mathbf{N}_0}\int^{\infty}_h\{|i_k(t, a)|\}da = 0$ uniformly for $i(t, a)\in\tilde{\Lambda}_2(t)B$.

From (2.11), $\sum_{k\in\mathbf{N}_0}\int^{\infty}_h\{|i_k(t, a)|\}da = 0$ for $h\geq t$, and then the criterion (ⅱ) holds for the set $\tilde{\Lambda}_2(t)B$. To show the criterion (ⅰ), we directly calculate for $k\in\mathbf{N}_0$ and $h\leq t$

Note that

and

According to the assumption $(A_2)$, there exists a positive constant $h_0$ such that $\frac{\Omega_{kj}(0, a+h)}{\Omega_{kj}(0, a)} < 2$ for $|h| < h_0$ since that $\Omega_{kj}(0, a)$ is continuous in $a$. Hence,

Solve the equation (2.10) to obtain, for $j\in\mathbf{N}_0$,

Accordingly, we have that

From (2.15), (2.16) and the Lipschitz constant, $l_k$, of $\tilde{\mathcal{P}}_k$, it holds that

By (2.12)-(2.14) and (2.17), it concludes that the criterion (ⅱ) holds for the set $\tilde{\Lambda}_2(t)B$, and then $\tilde{\Lambda}_2(t)B$ is compact. From [25,Lemma 3.2.3], the assertion is true.

From Lemma 2.2, the semiflow $\Phi(t)$ is point dissipative and then there exists a positively invariant absorbing set under the semiflow $\Phi(t)$, and then $\Phi(t)$ maps any bounded subset of $\mathcal{Y}_+$ to a precompact set in $\mathcal{Y}_+$ and is then compact for any $t > 0$. From Lemma 2.3, the semiflow is asymptotically smooth. Combining these properties and the result on the existence of global attractors in [25,Theorem 3.4.6] (or see [28,Theorem 2.6]), the following lemma is implied.

Lemma 2.4. The solution semiflow $\Phi(t)$ generated by system (2.1) in $\mathcal{Y}_+$ admits a compact global attractor $\mathcal{A}\subset \mathcal{Y}_+$.

3.   Basic reproduction number

First, we demonstrate the existence of infection-free equilibrium by examining the convergence dynamics of the no-infection model.

Lemma 3.1. ([29,Sec. 2]) Consider the no-infection model

There exists a positive equilibrium $(\bar{T}_1, \cdots, \bar{T}_n)$ that is globally asymptotically stable with respect to $\mathbb{R}^n_+$.

Hence, the model (2.1) admits an infection-free equilibrium $\bar{E} = (\bar{T}, 0_{\mathbb{R}^n})^{\tau} = (\bar{T}_1, \cdots, \bar{T}_n, 0, \cdots, 0)^{\tau}$. The linearization of (2.4) at the equilibrium $\bar{E}$ is the system

To determine the basic reproduction number of (2.1), we consider the decoupled $i$-equation in (3.2):

Define the operators $\hat{A}: L_+\rightarrow L_+$ by

and

Let $\tilde{\Phi} = \mathcal{A}+\mathcal{B}$ and $\tilde{\Phi}_0$ be the restriction of $\tilde{\Phi}$ in $\mathcal{Y}_0$; that is, $Dom(\tilde{\Phi}_0) = \{(0, f)^{\tau}:\tilde{\Phi} f\in \mathcal{Y}_0\}$. It is easy to see that $(0, f)^{\tau}\in Dom(\tilde{\Phi}_0)$ if it holds that $f(0) = \int^{\infty}_0P(a)f(a)da$. From the argument in [30,Sec. 6], the basic reproduction number of system (2.1) is

where $\rho(\cdot)$ denotes the spectral radius of an operator (also for that of a matrix). In fact, for $(-\mathcal{A}^{-1})(x, g)^{\tau} = (0, f)^{\tau}$, it holds that

and then

Note that the operator $-\mathcal{B}\mathcal{A}^{-1}$ has the same spectral radius on $\mathcal{Y}$ and $\mathbb{R}^n\times \{0_{L}\}$. Hence,

By the result in ^[30] (Theorems 3.16, 3.17), $\mathcal{R}_0$ plays the threshold value to determine local stability of the infection-free equilibrium, as stated in the following theorem.

Theorem 3.1. The infection-free equilibrium $\bar{E}$ is locally asymptotically stable (unstable) in (2.1) if $\mathcal{R}_0 < 1$ ($\mathcal{R}_0 > 1$).

Now, as a concrete example we consider a two-compartment system, $n = 2$, with $\delta_k(\cdot) = \delta(\cdot)$ for $k = 1, 2$. Then $\Theta = \Theta_2$, where

Through a direct calculation, we derive that

where

Since the matrix $\Theta_2$ is nonnegative and irreducible, $\rho(\Theta_2)$ is real and positive ^[24]. Denote $\rho(\Theta_2) = e^{\xi}$ for $\xi\in\mathbb{R}$. Then

To explore the effect of the migration rate on the dynamics of (2.1), we choose $p_k(a) = e^{-pa}$, $\delta_k(a) = \delta$, $m_k = m$, $k = 1, 2$, for positive $p$, $\delta$ and nonnegative $m$. Then (3.6) is equivalent to

Referring to ^[31,32,33], we choose $\lambda_1 = \lambda_2 = 1$, $d_1 = d_2 = 1$, $\alpha_1 = 0.00065$, $\alpha_2 = 0.00035$, $p = 0.15$, $\delta = 0.4$ and adjust $m$ to determine the root of the function $G_m$. Figure 1 shows that (3.7) admits a positive root, which indicates that $\rho(\Theta_2) > 1$, when $m = 0$, and there is no nonnegative root to (3.7), which indicates that $\rho(\Theta_2) < 1$, when $m = 2$. Hence, different values of migration rate $m$ may affect the extinction or persistence of virus in (2.1). Explicitly, in this case, strong circulation facilitates clearing the virus. In Section 6, we shall discuss the effect of the migration rate on the viral dynamics in more cases.

4.   Threshold dynamics

4.1.   Extinction of virus

When the reproduction number $\mathcal{R}_0 < 1$, we will demonstrate that the virus becomes extinct. First, the following lemma claims that infection can not occur under non-viable initial infection.

Lemma 4.1. Consider the virus equation in an environment with a constant vector of healthy cells, $T_0$,

with $i_0(\cdot)\in L_+$ and $\int^{\infty}_0P(a)i_0(a)da = 0$. Then $i(t, \cdot) = 0(\cdot)$ for all $t\geq 0$.

Proof. The well know result in ^[34] of linear age-structured models reveals that

where $\tilde{Q}$ satisfies

The initial condition, with $\int^{\infty}_0P(a)i_0(a)da = 0$, leads to

Since $\tilde{Q}(0) = 0$, we have that $\tilde{Q}(t) = 0$ for all $t\geq 0$, and then $i(t, \cdot) = 0(\cdot)$ for all $t\geq 0$.

With both the multi-compartmental structure and infection-age, it is a challenge to demonstrate the global dynamics in (2.1) by constructing a Lyapunov function ^[20]. However, the following theorem shows that the virus goes to extinction in system (2.1) when the basic reproduction number is less than or equal to unity.

Theorem 4.1. When $\mathcal{R}_0\leq 1$, the infection-free equilibrium $\bar{E}$ is globally asymptotically stable in (2.1).

Proof. We first show the assertion for the case with $\mathcal{R}_0 < 1$. By the standard comparison theorem, the solution of (2.1) satisfies that

where $\tilde{T}_k(t)$ is the solution to (3.1). From the assumption $\mathcal{R}_0 < 1$, there exists a sufficiently small $\varepsilon > 0$ such that $\rho(\Theta_{\varepsilon}) < 1$, where

From Lemma 3.1 and (4.2), there is a $t_1 > 0$ such that $T_k(t)\leq \bar{T}_k+\varepsilon$, $k\in\mathbf{N}_0$, for $t\geq t_1$. Since the matrix $\Theta_{\varepsilon}$ is nonnegative and irreducible, $\rho\left(\Theta_{\varepsilon}\right)$ is a positive eigenvalue which corresponds to a positive eigenvector, say $v = (v_1, \cdots, v_n)$, with $v^{\tau}\Theta_{\varepsilon} = \rho(\Theta_{\varepsilon})v^{\tau}$. Define

We calculate the derivative of $W(t)$ along the solution of (2.4) and derive, for $t\geq t_1$,

and the equality holds if and only if $\int^{\infty}_0P(a)i(t, a)da = 0$ for $t\geq t_1$. We claim that the largest positively invariant subset $M_0$ of $\{\frac{dW(t)}{dt} = 0\}$ is $\{ (T, 0(\cdot)) | T\in\mathbb{R}^n \}$. For $\tilde{X}(t) = (\tilde{T}(t), \tilde{i}(t, a))\in M_0$, a solution of (2.4), there is a $\tilde{t}_1 > 0$ such that $\int^{\infty}_0P(a)\tilde{i}(t, a)da = 0$ for $t\geq \tilde{t}_1$. Denote $\tilde{q}_0(\cdot) = \tilde{i}(\tilde{t}_1, \cdot)$ and consider the time-rescaled solution $\tilde{Y}(t) = \tilde{X}(t+\tilde{t}_1)$, but remain the symbol $\tilde{Y}(t) = (\tilde{T}(t), \tilde{i}(t, a))$. Then

where $\hat{i}(t, a)$ satisfies (4.1) with $T_0 = \bar{T}+\varepsilon I$ and $\hat{i}_0(\cdot) = \tilde{q}_0(\cdot)$. Note that $\int^{\infty}_0P(a)\tilde{q}_0(\cdot)da = 0$. From Lemma 4.1, $\hat{i}(t, \cdot) = 0(\cdot)$ for all $t\geq 0$. From (4.3), $M_0$ is in fact the set$\{ (T, 0(\cdot)) | T\in\mathbb{R}^n \}$. In addition, by using the theory of asymptotically autonomous systems [35,Theorem 2.3] and ^[36], we show that $\bar{E}$ is globally asymptotically stable in system (2.1), when $\mathcal{R}_0 < 1$.

When $\mathcal{R}_0 = 1$, we show it by contradiction. Suppose there is a $k_0\in\mathbf{N}_0$ such that $\limsup_{t\rightarrow +\infty}\int^{\infty}_0 i_{k_0}(t, a)da > 0$, then by Lemma 4.1 we have

From $T$-equation in (2.1),

We claim that $\limsup_{t\rightarrow +\infty} T_{k_0} < \bar{T}_{k_0}$. Denoting $\tilde{\mathcal{P}}(t) = {\rm diag}(\tilde{\mathcal{P}}_1(t), \cdots, \tilde{\mathcal{P}}_n(t))$, we rewrite the $T$-equation in system (2.1) into the vector form

From the variation of constant formula, the solution to (4.4) satisfies

From Lemma 3.1, we have $\lim_{t\rightarrow +\infty} (e^{-t(D-M)}T(0) + \int^t_0e^{-(t-s)(D-M)}\Lambda ds) = \bar{T}$. We next estimate the third term in (4.5). Since the assumption $\epsilon_0 > 0$, there exist sequences of positive numbers $\{t^{(k_0)}_l\}$ and $\{\tau_l\}$ with $\lim_{l\rightarrow \infty} t^{(k_0)}_l = +\infty$ such that $\tilde{\mathcal{P}}_{k_0}(t)\geq \frac{\epsilon_0}{2}$ for $t\in [t^{(k_0)}_l-\tau_l, t^{(k_0)}_l+\tau_l]$.

We claim that $\tilde{\mathcal{P}}'_{k_0}(t)$ is a bounded function. From the assumption $(A_3)$, it holds true whenever $\int^{\infty}_0 \frac{\partial i_{k_0}(t, a)}{\partial t} da$ is a bounded function of $t$. In fact, (2.1) and (2.2) imply that

It shows that $\int^{\infty}_0 \frac{\partial i_{k_0}(t, a)}{\partial t} da$ is a bounded function of $t$ due to the assumptions $(A_2)$ and $(A_3)$ and Lemma 2.2. Hence $\tilde{\mathcal{P}}'_{k_0}(t)$ is a bounded function and then the sequence $\{\tau_l\}$ can be chosen such that $\lim_{l\rightarrow \infty} \tau_l > 0$. It is easy to show that $T(t)$ has a positive lower bound. In addition, since the matrix $e^{-(D-M)}$ is positive and irreducible, the $k_0$-th component of $\int^t_0e^{-(t-s)(D-M)} \Gamma\tilde{\mathcal{P}}(s) T(s) ds$ has a positive lower bound. Therefore, there exist $\varepsilon_1 > 0$ and $t_1 > 0$ such that $T_{k_0}(t) < \bar{T}_{k_0} - \varepsilon_1$ for $t > t_1$.

From the monotonicity of the spectral of nonnegative matrices ([37,Corollary (1.5) in Ch2]), and the assumption $\mathcal{R}_0 = 1$, we have $\rho(\Theta_{\varepsilon_1}) < 1$, where

In addition, since the spectral of a matrix is continuous on each component, there exists a sufficiently small $\varepsilon_2 > 0$ such that $\rho(\Theta_{\varepsilon_1, \varepsilon_2}) < 1$, where

By using the method in case of $\mathcal{R}_0 < 1$, we can show that $i_k(t, \cdot)\rightarrow 0(\cdot)$ as $t\rightarrow +\infty$ for each $k\in\mathbf{N}_0$. This contradicts to the assumption on $i_{k_0}(t, \cdot)$. Therefore, we conclude that $\bar{E}$ is globally asymptotically stable in system (2.1), when $\mathcal{R}_0 = 1$.

4.2.   Uniform persistence

In this subsection, we study the uniform persistence for the model (2.4) when $\mathcal{R}_0 > 1$, by using the persistence theory developed in ^[38] (or see ^[19]). Recall that the semiflow $\Psi(t)$ is generated by the Cauchy problem (2.7) which is equivalent to (2.4). Define

Lemma 4.2. The subsets $U$ and $\partial U$ are positively invariant under the semiflow $\Psi(t)$. Furthermore, $\lim_{t\rightarrow +\infty}\Psi(t)\mathbf{u} = \bar{\mathbf{u}} = (\bar{T}, 0_{\mathbb{R}^n}, 0_{L})^{\tau}$ for each $\mathbf{u}\in \partial U$.

Proof. First, we show the positive invariance of the set $U$. Let $\mathbf{u}_0 = (T_0, 0, i_0(\cdot))^{\tau}\in U$. Denote

and define

Since $\mathbf{u}_0\in U$, $\Upsilon(0) > 0$. Through a direct calculation, we have

where $\delta_{{\rm max}} = \max_{1\leq k\leq n}\{\delta_{k, {\rm max}}\}$. Hence, we obtain that

for $t\geq 0$, and then $\Psi(t)U\subset U$.

Next, in order to show the positive invariance of the set $\partial U$, we consider $\mathbf{u}_0\in \partial U$. It is clear from Lemma 2.2 that $T(t)\leq\hat{T}$ for some $\hat{T} > 0$. Then the comparison theorem implies that

where $\hat{i}(t, a)$ is the solution of the following system

Since $\int^{\infty}_0i_0(a)da = 0$ and the assumption $(A_3)$, it holds that $\int^{\infty}_0P(a)\hat{i}(0, a)da = 0$. From Lemma 4.1, we obtain that $\hat{i}(t, a) = 0(\cdot)$ for all $t\geq 0$. The comparison in (4.6) implies that $i(t, a) = 0(\cdot)$ for all $t\geq 0$ and then $\partial U$ is positively invariant under the semifolw $\Psi(t)$. In addition, from Lemma 3.1, it is clear that for the solution remaining in $\partial U$ we have $T(t)\rightarrow \bar{T}$. Hence, $\lim_{t\rightarrow +\infty}\Psi(t)\mathbf{u} = \bar{\mathbf{u}}$ for each $\mathbf{u}\in \partial U$.

Next, we demonstrate uniform persistence of system (2.7).

Theorem 4.2. When $\mathcal{R}_0 > 1$, the semiflow $\Psi(t)$ generated by (2.7) is uniformly persistent with respect to $(U, \partial U)$; that is, there exists a constant $\varsigma > 0$ such that for each $\mathbf{u}\in U$,

where $d(\cdot, \cdot)$ is the distance associated to the norm in (2.6). Furthermore, there exists a compact subset $\mathcal{A}^0\subset U$ which is a global attractor for $\{\Psi(t)\}_{t\geq 0}$ in $U$.

Proof. In the following, we will prove that $W^s(\{\bar{\mathbf{u}}\})\cap U = \emptyset, {\rm where}\; W^s(\{\bar{\mathbf{u}}\}) = \{\mathbf{u}\in \mathcal{X}_{0+} | \lim_{t\rightarrow +\infty}\Psi(t)\mathbf{u} = \bar{\mathbf{u}}\}$. Recall from Lemma 4.2 that the infection-free equilibrium $\bar{\mathbf{u}}$ is GAS in $\partial U$. It is sufficient to show that there exists $\sigma > 0$ satisfying for each $\mathbf{u}\in\{\mathbf{v}\in U | \|\mathbf{v}-\bar{\mathbf{u}}\|\leq\sigma\}$ there exists $\tilde{t}_0\geq 0$ such that $\|\Psi(\tilde{t}_0)\mathbf{u}-\bar{\mathbf{u}}\| > \sigma$. Suppose by contradiction that for each integer $m\geq 0$ there is a $\mathbf{u}_m\in\{ \mathbf{v}\in U | \|\mathbf{v}-\bar{\mathbf{u}}\|\leq\frac{1}{m+1} \}$ such that $\|\Psi(t)\mathbf{u}_m-\bar{\mathbf{u}}\|\leq\frac{1}{m+1}$ for $t\geq 0$. Denote $\Psi(t)\mathbf{u}_m = (T^m(t), 0_{\mathbb{R}^n}, i^m(t, a))^{\tau} = (T^m_1(t), \cdots, T^m_n(t), 0_{\mathbb{R}^n}, i^m_1(t, a), \cdots, i^m_n(t))^{\tau}$ then $|T^m_j(t)-\bar{T}_j|\leq\frac{1}{m+1}, \; j\in\mathbf{N}_0$, for all $t\geq 0$. Consider

with $T^m_{0} = (T^m_{10}, \cdots, T^m_{n0})^{\tau}\geq 0$, $i^m_0(\cdot) = (i^m_{10}(\cdot), \cdots, i^m_{n0}(\cdot))^{\tau}\geq 0(\cdot)$ and $\sum_{k\in\mathbf{N}_0}\int^{\infty}_0i^m_{k0}(a)da > 0$. From the comparison theorem,

for $t\geq 0$, where $\tilde{i}^m(t, \cdot)$ is the solution of the following system

or equivalently

for $t\geq 0$, $\mathbf{x}(0)\in\overline{ \hat{\mathcal{D}}(\hat{A}) }$, the closure of $\hat{\mathcal{D}}(\hat{A}) = \{0_{\mathbb{R}^n}\}\times (W^{1, 1})^n$, where $\mathbf{x}(t) = (0_{\mathbb{R}^n}, \tilde{i}^m(t, \cdot))^{\tau}$ and the operators $\hat{A}, \hat{L}_m$ are defined as

Since $\mathcal{R}_0 > 1$, there is a $m_0 > 0$ large enough such that, for $m\geq m_0$,

The dominant eigenvalue $\lambda^*_m$ of system (4.9) satisfies the characteristic equation $\Delta_m(\lambda) = 0$, where

Note that $\mathcal{R}^m_0 > 1$ implies the existence of $\lambda^*_m > 0$ such that $\Delta_m(\lambda^*_m) = 0$. Denote the solution semiflow of (4.9) by $\{\hat{\Psi}_m(t)\}_{t\geq 0}$ and let $\hat{\Pi}_m(\hat{\mathbf{x}}):\hat{\mathcal{X}}\rightarrow\hat{\mathcal{X}}$, $\hat{\mathcal{X}}: = \mathbb{R}^n\times L$, be the projection of given $\hat{\mathbf{x}}\in \hat{\mathcal{X}}$ on the eigenspace associated to the dominant eigenvalue $\lambda^*_m$. By the result in [19,Section 2.3], we see an explicit expression for the following projection

for $\hat{\mathbf{x}}_0 = (0_{\mathbb{R}^n}, i^m_0(\cdot), )^{\tau}$ and

By using (4.11) and directly calculating $(\lambda I-\hat{A})^{-1}$, $(I-\hat{L}_m(\lambda I-\hat{A})^{-1})^{-1}$ respectively, we obtain that

if and only if

and

Since $\Delta_m(\lambda^*_m) = 0$, we have

By using ${\rm Re}\; \lambda^*_m > 0$, (4.10), (4.11) and (4.12), it deduces that $\lim_{t\rightarrow +\infty}\|\tilde{i}^m(t)\| = +\infty$. The comparison (4.8) implies $\lim_{t\rightarrow +\infty}\|v^m(t)\| = +\infty$, which is a contradiction to the boundedness of the solution. Hence, $W^s(\{\bar{\mathbf{u}}\})\cap U = \emptyset$. Accordingly, we derive from [38,Theorem 4.1,Theorem 4.2] that the semiflow $\{\Psi(t)\}_{t\geq 0}$ is uniform persistence with respect to $(U, \partial U)$. In addition, the result in [28,Theorem 3.7] implies that there exists a compact global attractor $\mathcal{A}^0\subset U$ for $\{\Psi(t)\}_{t\geq 0}$ in $U$.

Remark 1. According to Theorem 4.2, the uniform persistence of (2.4) indicates that there exists a constant $\zeta > 0$ such that for each solution in (2.4) satisfies

We further note that there exists a constant $\zeta_0 > 0$ such that

for each $k\in\mathbf{N}_0$. Otherwise, there is a $k_0\in\mathbf{N}_0$ such that $lim_{t\rightarrow\infty} \int^{\infty}_0i_{k_0}(t, a)da = 0$, which contradicts to the $i_{k_0}$-equation in (2.1) and the fact (4.13).

5.   Non-strongly connected within-host structure

The dichotomy of viral persistence or extinction within a host is based on the assumption of irreducible migration matrix. However, the migration matrix may be reducible due to functions of afferent or efferent lymphatics. In this case, the theory of direct graphs developed in ^[39] provides a method to describe directional connections between all within-host compartments (see also ^[40]). By building $n$ vertices and assigning a directed edge from vertex $j$ to vertex $k$ when a flow from the $j$-th compartment to $k$-th compartment is available, we can associate the matrix $M$ with a direct graph. Then the graph of $M$ reflects the connection structure within a host, and the matrix $M$ is irreducible if and only if for any pair of two compartments there is a path from one compartment to the other. In such a case, we say that the viral environment assumes a strongly connected within-host structure, whereas the case with a reducible matrix $M$ is referred to non-strongly connected structure.

When the within-host structure in non-strongly connected, that is, the matrix $M$ is reducible, one can use a permutation operator to reach a triangular block form, still denoted by matrix $M$,

where $p\leq n$ and each $M_{ll}$, $1\leq l \leq p$, is a $r_l\times r_l$ irreducible square matrix. We denote $\mathbf{P}_0 = \{1, 2, \cdots, p\}$ and write $\mathbf{N}_0 = \cup_{l\in\mathbf{P}_0} \mathbf{N}_0^{(l)}$, where $\mathbf{N}_0^{(l)} = \{\Sigma^{l-1}_{s = 1} r_s+1, \Sigma^{l-1}_{s = 1} r_s+2, \cdots, \Sigma^{l}_{s = 1} r_s \}$. By this way, we divide the whole within-host environment into $p$ parts of which each part consists of strongly connected compartments. We indicate, by block-$l$, the part consisting of $k$-th compartment for $k\in\mathbf{N}_0^{(l)}$. Explicitly, given $l\in \mathbf{P}_0$ with $r_l\geq 2$, the compartments in block-$l$ are strongly connected, whereas $r_l = 1$ means that the block-$l$ consists of a single compartment. For later convenience, we denote $\mathbf{N}_0^{(l)} = \{\xi_l, \cdots, \xi^l \}$, that is, $\xi_l = \Sigma^{l-1}_{s = 1} r_s+1$ and $\xi^l = \Sigma^{l}_{s = 1} r_s$. Notably, under the form of (5.1), there is a connection from the block-$l_1$ into the block-$l_2$ only for $1\leq l_1 < l_2\leq p$. Denote $\underline{\mathbf{P}}_0^{(l)} = \{\tilde{l} | 1\leq\tilde{l} < l, \; M_{l\tilde{l}} \neq0 \}$ for $1 < l\leq p$ and $\bar{\mathbf{P}}_0^{(l)} = \{\tilde{l} | l < \tilde{l}\leq p, \; M_{\tilde{l}l} \neq0 \}$ for $1\leq l < p$. $\underline{\mathbf{P}}_0^{(l)}$ is the set of indexes of blocks where each one is connected to block-$l$ through a directional pathway, whereas $\bar{\mathbf{P}}_0^{(l)}$ is that of blocks where each one is connected by a directional pathway from block-$l$.

Now, for $l\in \mathbf{P}_0$, we say that the block-$l$ is infection free if

for all solutions, and is infected if the infected population is uniformly persistent, that is, from Remark 1, there exists a constant $\zeta_0 > 0$ such that

We set the order of the blocks according to (5.1) and employ mathematical induction to analyze the virus dynamics in (2.1). We first introductorily define a sequence of basic reproduction numbers. In the system with only block-1, the basic reproduction number can be formulated as in (3.5) via replacing the matrix $M$ by $M_{11}$. However, when there are more than one block, it is necessary to reformulate the threshold value to determine the viral dynamics. Explicitly, for $k\in\mathbf{N}_0^{(1)}$, we have

Denote $X(t) = (T_1(t), \cdots, T_{\xi^1}(t))^{\tau}$, $y(t, a) = (i_1(t, a), \cdots, i_{\xi^1}(t, a))^{\tau}$, $\mathcal{L}^{(1)} = (\lambda_1, \cdots, \lambda_{\xi^1})^{\tau}$, $\mathcal{D}^{(1)} = {\rm diag}(d_1, \cdots, d_{\xi^1})$, $\mathcal{G}^{(1)} = {\rm diag}(\alpha_1, \cdots, \alpha_{\xi^1})$, ${\rm diag}(X(t)) = {\rm diag}(T_1(t), \cdots, T_{\xi^1}(t))$, $\mathcal{\tilde{D}}^{(1)}(a) = {\rm diag}(\delta_1(a), \cdots, \delta_{\xi^1}(a))$, $\mathcal{P}^{(1)}(a) = {\rm diag}(p_1(a), \cdots, p_{\xi^1}(a))$ and

with

We rewrite (5.2) into

with

In order to determine the threshold value to conclude viral persistence or extinction in (5.2), we denote

and $\bar{T}^{(1)}$ satisfies $\mathcal{L}^{(1)}-\mathcal{D}^{(1)}\bar{T}^{(1)}+\mathcal{M}^{(1)}\bar{T}^{(1)} = 0$. Note that $\mathcal{R}_0^{(1)}$ is not the basic reproductive number for isolated block-1 except that $\bar{\mathbf{P}}_0^{(1)} = \emptyset$.

Next, we start the process to determine the viral dynamics from the first block and the following is a straightforward result from Theorem 4.1 and Theorem 4.2.

Proposition 5.1. The first block is infection free when $\mathcal{R}_0^{(1)}\leq 1$, whereas it is infected when $\mathcal{R}_0^{(1)} > 1$.

For given $l\in \mathbf{P}_0$, we suppose that from block-1 to block-($l-1$) are all determined to be infected or infection free. In an infection free block, it satisfies $\lim_{t\rightarrow}T_k(t) = \bar{T}_k$ for constants $\bar{T}_k > 0$ derived as in Lemma 3.1 and the virus goes extinction. Next, we check on the connection from these $l-1$ blocks to the block-$l$. Explicitly, we say that the block-$l$ is susceptible from an infected block if there exists $\tilde{l}\in\{1, 2, \cdots, l-1\}$ such that the block-$\tilde{l}$ is infected and $M_{l\tilde{l}}\neq \mathbf{0}$.

Here, we denote by an index to determine whether the block-$l$ is susceptible from an infected block. First, we denote

Then there is a (no) connection from the block-$\tilde{l}$ to the block-$l$ when $\mathcal{S}(M_{l\tilde{l}}) = 1$ ($\mathcal{S}(M_{l\tilde{l}}) = 0$). We further denote

Then the block-$l$ is susceptible from an infected $\tilde{l}$ block with $\tilde{l} < l$ if $\chi^{(l)} = 1$, whereas it is not susceptible from any infected block-$\tilde{l}$ with $\tilde{l} < l$ if $\chi^{(l)} = 0$.

Suppose that the block-$l$ is not susceptible from any infected block-$\tilde{l}$ with $\tilde{l} < l$, i.e. $\chi^{(l)} = 0$. Then the subsystem with $k\in\mathbf{N}_0^{(l)}$ has the limit system

Denote $X(t) = (T_{\xi_l}(t), \cdots, T_{\xi^l}(t))^{\tau}$, $y(t, a) = (i_{\xi_l}(t, a), \cdots, i_{\xi^l}(t, a))^{\tau}$, $\mathcal{L}^{(l)} = (\bar{\lambda}_{\xi_l}, \cdots, \bar{\lambda}_{\xi^l})^{\tau}$,

$\mathcal{D}^{(l)} = {\rm diag}(d_{\xi_l}, \cdots, d_{\xi^l})$, $\mathcal{G}^{(l)} = {\rm diag}(\alpha_{\xi_l}, \cdots, \alpha_{\xi^l})$, ${\rm diag}(X(t)) = {\rm diag}(T_{\xi_l}(t), \cdots, T_{\xi^l}(t))$, $\mathcal{\tilde{D}}^{(l)}(a) = {\rm diag}(\delta_{\xi_l}(a), \cdots, \delta_{\xi^l}(a))$, $\mathcal{P}^{(l)}(a) = {\rm diag}(p_{\xi_l}(a), \cdots, p_{\xi^l}(a))$ and

with

We rewrite (5.6) into

with

In order to determine the threshold value to conclude viral persistence or extinction in (5.8), we denote

and $\bar{T}^{(l)}$ satisfies $\mathcal{L}^{(l)}-\mathcal{D}^{(l)}\bar{T}^{(l)}+\mathcal{M}^{(l)}\bar{T}^{(l)} = 0$. Note that $\mathcal{R}_0^{(l)}$ is not the basic reproductive number for isolated block-$l$ except that $\underline{\mathbf{P}}_0^{(l)} \cup \bar{\mathbf{P}}_0^{(l)} = \emptyset$.

Accordingly, the virus dynamics in the block-$l$ can be determined as the following.

Theorem 5.2. For $1 < l\leq p$, if $\chi^{(l)} = 1$, that is, the block-$l$ is susceptible from an infected block-$\tilde{l}$ with $\tilde{l} < l$, then the block-$l$ is infected; if $\chi^{(l)} = 0$, that is, the block-$l$ is not susceptible from any infected block-$\tilde{l}$ with $\tilde{l} < l$, then the block-$l$ is infection free when $\mathcal{R}_0^{(l)}\leq 1$, whereas it is infected when $\mathcal{R}_0^{(l)} > 1$.

Proof. Consider the block-$l$ with $\mathcal{X}^{(l)} = 1$. It holds that, for $k\in\mathbf{N}_0^{(l)}$,

Also note that there exists $t_0$ and $\check{T}_k$ such that $T_k(t)\geq\check{T}_k$ for all $k\in\mathbf{N}_0^{(l)}$ for $t\geq t_0$. Hence, we consider an auxiliary system

with $\tilde{i}_k(t, 0) = \alpha_k\check{T}_k\int^{\infty}_0p_k(a)\tilde{i}_k(t, a)da$ and $\tilde{i}_k(0, a) = i_{k0}(a)$. By a comparison theory, it reveals that $i_k(t, a)\geq \tilde{i}_k(t, a)$ for $a\geq 0$ and $t\geq t_0$. Denote $\tilde{I}_k(t) = \int^{\infty}_0\tilde{i}_k(t, a)da$ for $k\in\mathbf{N}_0^{(l)}$, then it satisfies that

where $\tilde{\zeta}_k: = \sum_{j\in\cup_{\tilde{l} = 1}^{l-1}\mathbf{N}_0^{(\tilde{l})}}m_{kj}\zeta_0^{(\tilde{l})}$ and $\zeta_0^{(\tilde{l})}$ is the lower bound in (4.14). Note that the assumption $\mathcal{X}^{(l)} = 1$ implies $\tilde{\zeta}_k > 0$. We further consider the auxiliary system

Since the matrix $M_{ll}$ is irreducible, there exists a positive equilibrium which is globally asymptotically stable as in Lemma 3.1. Thus we reveal that

for some positive constants $\check{\zeta}_k$. In addition, since $i_k(t, a)\geq \tilde{i}_k(t, a)$ for $a\geq 0$ and $t\geq t_0$, we see that $\int^{\infty}_0i_k(t, a)da \geq \int^{\infty}_0\tilde{i}_k(t, a)da$ and then

which concludes that the virus population uniformly persists in block-$l$.

Next, we suppose that the block-$l$ is not susceptible from any infected block-$\tilde{l}$ with $\tilde{l} < l$, that is $\mathcal{X}^{(l)} = 0$. Then in the block-$l$ we obtain a limit system as in (5.4) for $k\in\mathbf{N}_0^{(l)}$. Therefore, the value of $\mathcal{R}_0^{(l)}$ will determine the viral dynamics and it completes the proof.

Remark 2. Note that, for $k\in\mathbf{N}_0^{(l)}$, the value of $\bar{\lambda}_k$ in (5.7) may depend on the connection from block-$\tilde{l}$ with $\tilde{l} < l$. From Theorem 5.2, we see that even if the block-$l$ is not susceptible from any infected block-$\tilde{l}$ with $\tilde{l} < l$, the involved viral dynamics may be affected by the connection from an infection free block.

6.   Effects of connection structure and distinct compartments

Considering each compartment with or without afferent and efferent lymphatics, there are several different connection structures to build the whole system. In this section, we investigate the viral extinction and persistence in the model (2.1) in cases of strongly and non-strongly connected structures, along with compartments containing distinct infection characteristics.

6.1.   A strongly connected circular chain

Deep lymph nodes of the head and neck are arranged in a vertical chain along the internal jugular vein and a circular chain consisting of occipital nodes, submental nodes, submandibular nodes, buccal or facial nodes and parotid nodes ^[41]. Hence, it is necessary to consider the group as a circular chain and to study the influence of removing one connection between lymphoid tissues on the viral infection. In this subsection, we consider three-compartment models with two types of connection matrices; one with a complete connection and the other without the path from compartment 3 to compartment 2, see Figure 2. To explore the effect of the connection topology on the viral infection, we assume identical compartments, $\lambda_k = \lambda$, $d_k = d$, $\alpha_k = \alpha$, $p_k(\cdot) = p(\cdot)$, and $\delta_k(\cdot) = \delta(\cdot)$, $k = 1, 2, 3$, and an identical migration rate, $m_{kj} = m$ for $k\neq j$. In the case with the complete connection, it yields that

A direct calculation gives

where $\phi(a) = \exp\left(-\int^a_0 \delta(s)ds\right)$. Denoting $\mathcal{Q}_1(m, a) = \frac{1}{3}+\frac{2}{3}e^{-3ma}$ and $\mathcal{Q}_2(m, a) = \frac{1}{3}-\frac{1}{3}e^{-3ma}$, we derive that

Note that $\mathcal{R}_0 > 1$ if and only if there exists a positive constant $\xi$ such that $\det(I-e^{-\xi}\Theta) = 0$, which is equivalent to that the determinant

where

Equation (6.1) holds if and only if

Hence, there exists a positive $\xi$ satisfying $\det(I-e^{-\xi}\Theta) = 0$ if and only if

that is

Thus, the migration of cells between three nodes has no influence on the viral dynamics.

However, in the case without the path from compartment 3 to compartment 2,

where the infection-free equilibrium $(\bar{T}_1, 0, \bar{T}_2, 0, \bar{T}_3, 0)$ is obtained as in Lemma 3.1. A direct calculation reveals that

and then

where

In this circular chain, we choose parameters $d = 0.09$, $\lambda = 10$, $p(a) = e^{-0.15a}$, $\delta = 0.4$ ^[42] and change the value of $\alpha$ to see possible effect of migration rate, $m$, on the viral persistence. When $\alpha = 0.0046$, there is a $m^* > 0$ such that $\mathcal{R}_0 < 1$ for $0\leq m < m^*$ and $\mathcal{R}_0 > 1$ for $m > m^*$. Hence sufficiently slow migration of cells will drives the virus to extinction, whereas faster migration of cells will lead to viral persistence, which is different from the non-influence of migration in the completely connected circular chain. Moreover, when $\alpha = 0.0042$, there exist positive constants $m^{**} < m^{***}$ such that $\mathcal{R}_0 < 1$ for $0\leq m < m^{**}$ or $m > m^{***}$ and $\mathcal{R}_0 > 1$ for $m^{**} < m < m^{***}$. It means that the influence of migration on the viral dynamics is complicated, from eradicating virus in the whole system with sufficiently small migration rate to initiating viral infection by medium value of migration rate, and eradicating virus again with sufficiently large migration rate.

6.2.   Two non-strongly connected blocks

Consider a host consisting of four compartments as in Figure 4, where the 1st and 2nd compartments (involved in block-1) are strongly connected and also for the 3rd and 4th compartments (involved in block-2). The only connection between the two blocks is from the 2nd compartment to the 3rd compartment, that is, the block-2 may be susceptible from the block-1 when the latter one is infected. Denote $m = m_{32}$. Hence, $\mathcal{R}_0^{(1)} = \rho(\Theta^{(1)})$, where

$\mathcal{D}^{(1)} = {\rm diag}(d_1, d_2)$, $\mathcal{G}^{(1)} = {\rm diag}(\alpha_1, \alpha_2)$, $\mathcal{\tilde{D}}^{(1)}(a) = {\rm diag}(\delta_1(a), \delta_2(a))$, $\mathcal{P}^{(1)}(a) = {\rm diag}(p_1(a), p_2(a))$,

and $\bar{T}^{(1)} = (\bar{T}^{(1)}_1, \bar{T}^{(1)}_2)^{\tau} = (\mathcal{D}^{(1)} - \mathcal{M}^{(1)})^{-1}\mathcal{L}^{(1)}$. By direct calculation, it yields that

The inner integral in $\Theta^{(1)}$ is a positive matrix because of the irreducibility of matrix $\mathcal{M}^{(1)}$, and then $\Theta^{(1)}$ is positive for all $m\geq 0$. We see how the value of $m$ affects that of $\mathcal{R}_0^{(1)}$ in the following. All proofs in this subsection will be postponed to Appendix.

Proposition 6.1. The value of $\mathcal{R}_0^{(1)}$ is decreasing in $m$.

On the other hand, $\mathcal{R}_0^{(2)} = \rho(\Theta^{(2)})$, where

$\mathcal{D}^{(2)} = {\rm diag}(d_3, d_4)$, $\mathcal{G}^{(2)} = {\rm diag}(\alpha_3, \alpha_4)$, $\mathcal{\tilde{D}}^{(2)}(a) = {\rm diag}(\delta_3(a), \delta_4(a))$, $\mathcal{P}^{(2)}(a) = {\rm diag}(p_3(a), p_4(a))$,

and $\bar{T}^{(2)} = (\bar{T}^{(2)}_1, \bar{T}^{(2)}_2)^{\tau} = (\mathcal{D}^{(2)} - \mathcal{M}^{(2)})^{-1}\mathcal{L}^{(2)}$. Accordingly, we can also see how the value of $m$ affects that of $\mathcal{R}_0^{(2)}$.

Proposition 6.2. The value of $\mathcal{R}_0^{(2)}$ is increasing in $m$.

Biologically, a unidirectional connection between two blocks benefits the one with exportation in viral eradication but damages the other one with importation.

In these two non-strongly connected blocks, we fix the viral environment in each block, including migration rates within the blocks, and change the migration rate between two blocks to explore its influence on the viral dynamics. Explicitly, we choose parameters relevant to block-1 by $\alpha_1 = \alpha_2 = 0.0085$, $d_1 = d_2 = 0.09$, $\lambda_1 = \lambda_2 = 10$, $m_{21} = m_{12} = 1$, $\delta = 0.4$, $p(a) = e^{-0.15a}$ and those relevant to block-2 by $\alpha_3 = \alpha_4 = 0.00275$, $d_3 = d_4 = 0.09$, $\lambda_3 = \lambda_4 = 10$, $m_{43} = m_{34} = 1$, $\delta = 0.4$, $p(a) = e^{-0.15a}$. As depicted in Figure 5, the value of $\mathcal{R}_0^{(1)}$ is decreasing in the migration rate, $m$, whereas the value of $\mathcal{R}_0^{(2)}$ is increasing in $m$. There exist positive constants $m_* < m_{**}$ with $\mathcal{R}_0^{(1)}(m_*) = 1$ and $\mathcal{R}_0^{(2)}(m_{**}) = 1$. Hence, according to Theorem 5.1, we observe the change of viral dynamics from viral persistence in the whole system with small value of $m$ ($0\leq m < m_*$) to virus extinction in the whole system with medium value of $m$ ($m_*\leq m \leq m_{**}$), and to virus persistence only in the block-2 when the value of $m$ is sufficiently large ($m > m_{**}$).

6.3.   Distinct gamma distributed infectivity kernels with delay

Recent studies show that for HIV the viral production and cell infection kernels are gamma distributed ^[43,44] after an intracellular delay for integration of viral genome. The predominant infected cell type or infection mode may determine the shape characteristics of the kernels. In particular, cell-to-cell transmission results in faster infection kinetics ^[44], and certain tissue compartments may have cells arranged closer together facilitating this infection model ^[45]. In order to model, compartments with distinct infection kinetics, we consider infected cell death rates $\delta_k(a)$ and infectivity kernels $p_k(a)$ of the following piecewise form:

where $\tau_k$ is the intracellular delay, $\nu_k$ is the cell death rate after integration of viral genome, and $g_k(s)$ is the gamma-distributed ($g_k \sim \ \rm{Gamma}(\kappa_k, \theta_k)$) infectivity kernel, for compartment $k\in\mathbf{N}_0^{(l)}$. Note that the mean of the gamma distribution $g_k(\kappa_k, \theta_k)$ is $\mu_k = \kappa_k/\theta_k$, representing the mean time of (secondary) cell infection for an infected cell in the absence of death. In the special case that each compartment has identical delays $\tau_k = \tau$ and gamma p.d.f.s, $g_k = g(\kappa, \theta)$, with shape parameter $\kappa = n$ as a positive integer, then the linear chain trick can yield the following equivalent system of delay differential equations:

where $Y_{m}(t): = \int_0^{\infty} \frac{\theta^m a^{m-1}}{(m-1)!} e^{-\theta a} i(t, a+\tau) da, \quad m = 1, \dots n$, $\tilde D = {\rm diag}(\nu_k)$, and $T(t-\tau)\circ Y_n(t-\tau)$ refers to entry-wise product of the vectors. Note that previous examples fall under this scenario with $n = 1$ and $\tau = 0$. Thus for the previous results on how migration rate affects $\mathcal R_0$ in Figure 3 and Figure 5, we can verify extinction versus persistence by transforming to the ODE to obtain numerical solutions. For the general case of patches with distinct intracellular delays and/or gamma distribution parameters in (6.4), however the system will be infinite-dimensional as our original PDE model (2.1) and we utilize a finite-difference method for example simulations below.

We explore the dynamics in this general setting by considering parameters representative of the study ^[44], where distinct modes of infection produce different infectivity kernels. Here we consider the case of three compartments connected in a circular chain lacking path from compartment-3 to compartment-2, as in Figure 2(b), with infectivity kernels as (6.4). In Figure 6, we display time-dependent solutions for 3 distinct scenarios; all compartments utilize cell-free infection, all compartments utilize cell-to-cell infection, and finally a mix where the first two compartments display cell-free while compartment-3 utilizes cell-to-cell infection. We also investigate how varying intracellular delay, $\tau_3$, and gamma distribution mean, $\mu_3$, of compartment-3 affect the overall $\mathcal R_0$. The fixed parameters utilized are $d = 0.09$, $\lambda = 100$, $\nu = 0.4$, $\alpha = 0.0011$, $m_{ij} = m = 1$, and in Figure 6(a)-(c), intracellular delay $\tau = 2/3$ and gamma distribution mean $\mu_1, \mu_2, \mu_3$ either equal to $4/5$ (cell-free infection) or $1/2$ (cell-to-cell infection) with shape parameter $\kappa$ fixed at 4. Shorter values of $\tau_k$ and $\mu_k$ result in larger $\mathcal R_0$, along with faster infection dynamics and larger number of infected cells. Thus, compartments which facilitate cell-to-cell infection may be important for allowing the virus to persist and establish reservoirs, which has been found in studies of HIV under drug treatment ^[46] or during initial bottleneck of host infection ^[16].

7.   Conclusion and discussion

The assumption of well-mixed viral environment may ignore diversity of viral dynamics, for example, the existence of viral sanctuary sites and variant viral clearances, and then underestimate the infection within a host. We propose a multi-compartment model incorporating age-since-infection to imitate the migration of cells through afferent and efferent vessels between diverse organs. A threshold, depending on the rate of cell migration, the structure of organic connection and infection kinetics, was found to determine the extinction or persistence of infection. Accordingly, in this paper, we theoretically discussed the effect of the connection structure on the extinction or persistence of infection. Comparing systems with complete connection (with bi-directional connection between any two compartments of the system) and only lacking one connection, we see the influence of the migration rate in latter case. In both strongly connected systems, the influence may non-monotonously depend on the migration rate. As for a non-strongly connected system, we constructed a sequence of threshold values to conclude the infection pattern within a host and revealed that increasing the migration rate from a strongly connected block to another strongly connected block does not always help to eradicate virus in the whole system. Furthermore our inclusion of cell infection-age allows for organ-specific infectivity kinetics based on predominance of cell-free versus cell-to-cell transmission or target cell type. In particular, we demonstrated that individual compartments dominated by cell-to-cell infection may precipitate viral persistence in a strongly connected system with distinct gamma-distributed infectivity kernels representative of experimental studies. These findings complement and extend research on models consisting of a well-mixed compartment.

In this study, the global convergence to the infection-free equilibrium was demonstrated in the case when the basic reproduction number is less than or equal to unity. On the other hand, the global convergence to the unique infection state is a challenging problem. In ^[20] for studying a multi-group SIR epidemic model with an age structure, an analogous problem was partially solved by using the method of Lyapunov functional and graph theory. However, the convergence dynamics to an infection state in this study remains as an unsolved problem since the model incorporates both the age structure and the migration behavior of cells. Moreover, the connection topology between a group significantly affects the evolution dynamics, and hence attracts researchers to contribute toward either theoretical analysis or application-oriented studies. A portion of a host body can be imitated by a strongly connected structure in which the viral dynamics holds clearance or uniformly persists in the complete part. However, biological evidences show that a functional system within a host, for example, the lymphatic system or the blood circulation system, consists of several such portions which are not strongly connected. Motivated by the study in ^[40] for infectious diseases via transportation networks, we systematically constructed a sequence of threshold values to find out the infection pattern in the whole system.

We explicitly explored the effect of connection topology or cell infection-age kinetics on the viral dynamics in two basic structures: (a) a strongly connected circular chain with or without a certain directional connection, and (b) a system with one pair of two-compartment blocks connected via a unidirectional pathway. In the first example, either with or without the indicated directional connection, it remains the type of strongly connected. Nevertheless, in the latter case, it may show multiple switches between viral eradication and viral infection in the whole system when increasing the identical migration rate of cells. The second example was devoted to study the connection topology between two strongly-connected blocks. In particular, we explored the effect of the unidirectional migration between two compartments, each is located in a strongly-connected block. An explicit formula shows that increasing the migration rate between two blocks may first change the viral persistence only in the first block to viral eradication in the whole system, and then, for even larger migration rate, trigger the infection in only the flow-in block. Biologically, this unidirectional connection between two blocks benefits the one with exportation in viral eradication but damages the other one with importation. Although these example structures may be limited in number of compartments, they reveal the necessity in distinguishing the connection topology within a host. Finally, we investigated how incorporating intracellular delays and distinct gamma distributed infectivity kernels based on cell-free and cell-to-cell infection data ^[44] can affect viral dynamics in the strongly connected circular chain. We find that the faster kinetics associated with cell-to-cell infection promote viral persistence and spread, even when only a single compartment in the chain has this transmission mode.

Variable time-dependent drug efficacy is another issue that requires consideration. For drug treatments, such as those administrated through targeted therapy, the treatment focuses on the targeted organ (the target compartment), and the drug efficacy is constant in each compartment. However, under several traditional therapies, drugs may spread throughout the entire body of a host via the circulatory system, in which the drug efficacy would no longer be constant in each compartment. Therefore, it is necessary to formulate the drug circulation among neighboring compartments along afferent and efferent vessels. These considerations will be explored in our next research.

Acknowledgments

C. J. Browne was partially supported by a U.S. National Science Foundation grant (DMS-1815095); C.-Y. Cheng was partially supported by the Ministry of Science and Technology of Taiwan, R.O.C. (MOST 108-2115-M-153-003).

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

To show the dependence of $\mathcal{R}_0^{(1)}$ and $\mathcal{R}_0^{(2)}$ on the value of $m$, we begin first from a property of the following exponential matrix $S_0(\beta)$.

Lemma A.1.

Given $c, d>0$, the value of each component of the matrix

is decreasing in $\beta$, for $\beta>0$.

Proof. The following matrix can be written into with its diagonal form

where

with

Note that $\psi_1(\beta)>0$, $(\psi_1(\beta))^2-(\psi_2(\beta))^2 = 4cd>0$, and then $\psi_1(\beta)-\psi_2(\beta)>0$ and $\psi_1(\beta)+\psi_2(\beta)>0$. In addition, direct calculations give

Hence, from (7.1), we have

where

By direct calculations and using (7.2), it yields that

Define

It is easy to show that $k(x)>0$ for all $x>0$ and $h(x, y) < 0$ for all $x>0$. Finally, the fact $\psi_1(\beta)>0$ implies $S'_{ij}(\beta) < 0$ for $i, j = 1, 2$, and it completes the assertion.

Based on this property, we see the dependence of $\mathcal{R}_0^{(1)}$ on $m$.

Proof of Proposition 6.1:

Proof. Since the spectral of a positive matrix positively depends on each component, it suffices to show that, for a fixed $a>0$, the value of each component of the matrix

is decreasing in $m$. For given $a>0$, if $\int^a_0\delta_2(s)da\geq \int^a_0\delta_1(s)da$, we rewrite $S(m)$ into

where

Note that $\beta(m)>0$ and $\beta'(m)>0$ for all $m>0$. From Lemma A.1, the value of each component of the matrix $\exp(-\int^a_0(\mathcal{\tilde{D}}^{(1)}(s)-\mathcal{M}^{(1)})ds)$ is decreasing in $m$. A similar argument can be proceeded by exchanging $\delta_1(\cdot)$ and $\delta_2(\cdot)$ in (7.3) when $\int^a_0\delta_1(s)da\geq \int^a_0\delta_2(s)da$. In addition, since the spectral of a positive matrix positively depends on each component, the value of $\mathcal{R}_0^{(1)}$ is decreasing in $m$.

Proof of Proposition 6.2:

Proof. In $\Theta^{(2)}$, $\bar{T}^{(2)}$ is the only term tdepending on the value of $m$ and, since the inner integral is a positive matrix, the spectral of $\Theta^{(2)}$ depends increasingly on $\bar{T}^{(2)}_2$. It is easy to see that the spectral bound of the matrix $-(\mathcal{D}^{(2)} - \mathcal{M}^{(2)})$ is negative and then by the Perron-Frobenious Theorem, the matrix $(\mathcal{D}^{(2)} - \mathcal{M}^{(2)})^{-1}$ is positive. Hence, the value of $\bar{T}^{(2)}_2$ positively depends on the value of $F(m): = m\bar{T}^{(1)}_2(m)$. Moreover, since

we conclude that the value of $\mathcal{R}_0^{(2)}$ is increasing in $m$.