Investigation of a Two-Patch Within-Host Model of Hepatitis B Viral Infection

Keoni Castellano; Omar Saucedo; Stanca M. Ciupe

doi:10.58301/gxvv8f17

Home Journals Meridian Journal of Mathematical SciencesVolume 10, Issue 410.58301/gxvv8f17

Meridian Journal of Mathematical Sciences

Article10 November 2024

Investigation of a Two-Patch Within-Host Model of Hepatitis B Viral Infection

Keoni Castellano^1† (unauthenticated), Omar Saucedo^1† (unauthenticated), Stanca M. Ciupe^1†* (unauthenticated)

¹https://ror.org/02smfhw86grid.438526.e0000 0001 0694 4940Department of Mathematics, Virginia Tech, 225 Stanger St., Blacksburg, 24060 VA United States; https://ror.org/02smfhw86grid.438526.e0000 0001 0694 4940Virginia Tech Center for the Mathematics of Biosystems, Virginia Tech, 765 West Campus Drive, Blacksburg, 24060 VA United States

† These authors contributed equally to this work.

* Corresponding author.

Meridian J. Math. Sci. 2024, 10(4), 181; https://doi.org/10.58301/gxvv8f17

Abstract

Chronic infection with hepatitis B virus (HBV) can lead to formation of abnormal nodular structures within the liver. To address how changes in liver anatomy affect overall virus-host dynamics, we developed within-host ordinary differential equation models of two-patch hepatitis B infection, one that assumes irreversible and one that assumes reversible movement between nodular structures. We investigated the models analytically and numerically, and determined the contribution of patch susceptibility, immune responses, and virus movement on within-patch and whole-liver virus dynamics. We explored the structural and practical identifiability of the models by implementing a differential algebra approach and the Monte Carlo approach for a specific HBV data set. We determined conditions for viral clearance, viral localization, and systemic viral infection. Our study suggests that cell susceptibility to infection within modular structures, the movement rate between patches, and the immune-mediated infected cell killing have the most influence on HBV dynamics. Our results can help inform intervention strategies.

Keywords:Hepatitis B virusAsymptotic dynamicsIdentifiability analysisTwo-patch modelOrdinary differential equations

Introduction

Chronic hepatitis B, caused by the hepatitis B virus (HBV), remains a major health problem, affecting more than 257 million carriers globally (Schweitzer et al. 2015; Trépo et al. 2014). HBV is transmitted vertically from mother to child, through sexual contact, and after exposure with infected blood or body fluids (Guidotti and Chisari 2006). The risk of progression from acute to chronic disease is age dependent, with adult immunocompetent patients clearing the virus and infants and children developing chronic infection, followed by liver disease such as cirrhosis and hepatocellular carcinoma (Lee 1997; Locarnini et al. 2015). Chronic hepatitis B is rarely cured with available therapeutics (Dusheiko 2013; Pawlotsky 2023), hence the need for interdisciplinary studies where mathematical modeling and empirical data can merge to better describe the disease evolution following hepatitis B viral infection.

Current mathematical models describing HBV dynamics from acute to chronic infection are based on classical within-host models (Nowak et al. 1996) (with variations to account for hepatitis B virus characteristics), and are mainly validated with HBV DNA titers in the serum of infected patients and chimpanzees (Ciupe, Ribeiro, Nelson, Dusheiko, and Perelson 2007; Ciupe et al. 2014; Ciupe, Ribeiro, Nelson, and Perelson 2007; Murray et al. 2005; Dahari et al. 2009; Neumann et al. 2010; Goyal et al. 2017; Tsiang et al. 1999; Ciupe et al. 2024, 2021). Recently, multiscale mathematical models have been developed, with the aim of integrating intracellular aspects of hepatitis B replication into the cellular dynamics (Reinharz et al. 2021; Ciupe et al. 2024; Hailegiorgis et al. 2023; Shekhtman et al. 2024; Kadelka et al. 2021). They are validated using in-vitro and humanized chimeric or transgenic mice data, and have been instrumental in addressing the role of immune system (both cytolytic killing and non-cytolytic cure) on viral resolution, and the role of viral markers (such as the s- and e-antigens, and the covalently closed circular DNA) on viral persistence (Shekhtman et al. 2024; Kadelka et al. 2021; Ciupe et al. 2024).

One aspect of chronic hepatitis B that is less studied, using mathematical models, concerns the effect that hepatitis B virus has on the reshaping of the liver as the disease progresses from healthy, to fibrotic, to cirrhotic (Fig. 1A). The progression to fibrotic and cirrhotic disease is driven by repeated cycles of inflammatory responses, liver cells (hepatocytes) death, hepatocytes regeneration, and accumulation of fibrous tissue (Roehlen et al. 2020; Yuan et al. 2019). All these processes lead to the activation of collagen-producing cells, resulting in excessive accumulation of extracellular matrix, the formation of abnormal nodular structures, overall liver dysfunction, and often the need for liver transplantation (Bataller et al. 2005).

**Fig. 1** Proposed model of liver disease evolution following hepatitis B viral infection: (A) Liver transition from healthy to fibrosis to cirrhosis; (B) Hepatitis B virus localization during acute and chronic stages of hepatitis B viral infection. During chronic disease, liver is restructured into abnormal nodular structures separated by collagen (yellow streaks, right panel). Figure created in https://BioRender.com

The death of hepatocytes, due to chronic hepatitis B infection, has been shown to trigger compensatory proliferation of the remaining hepatocytes (Michalopoulos and Bhushan 2021). Fibrosis and cirrhosis affect viral replication, with increased liver stiffness that promotes the progression of hepatitis B disease (Cast et al. 2019). In this paper, we investigated whether the HBV persistence is influenced by differences in hepatocyte turnover (replication and death) within the nodular structures formed during fibrosis and cirrhosis (Fig. 1B, right panel). To address this, we developed two-patch mathematical models of heterogeneous hepatitis B virus replication and used them to determine how patch-specific hepatocyte susceptibility, immune response-induced hepatocyte death, and HBV travel rates between patches influence viral persistence. The results can help elucidate how the host-virus interactions in a heterogeneous architecture influence the transition between viral resolution and viral persistence.

Mathematical Model Development

Single-patch model of hepatitis B virus infection. The mathematical model of homogeneous hepatitis B virus infection considers the interaction between uninfected liver cells T, infected liver cells I, and hepatitis B virus V. It assumes that uninfected liver cells are produced at rate s, die at natural death rate d, and are infected at rate

β

. Infected liver cells die at immune-response mediated rate

δ

(

\geq d

) and produce new virus at rate p. Virus is cleared at rate c. The mathematical model is given by:

\begin{matrix} \begin{matrix} \frac{dT}{dt} & = s - d T - β T V, \\ \frac{dI}{dt} & = β T V - δ I, \\ \frac{d V_{1}}{dt} & = p I - c V, \end{matrix} \end{matrix}

1

with initial conditions

T_{1} (0) = s / d

,

I (0) = I_{0}

,

V (0) = V_{0}

. This model was first used (in the context of hepatitis B virus infection) by Nowak et al. (1996) and utilized further in Ciupe, Ribeiro, Nelson, and Perelson (2007); Murray et al. (2005); Dahari et al. (2009); Tsiang et al. (1999); Neumann et al. (2010); Ciupe et al. (2024); K. Wang et al. (2010); Yousfi et al. (2011); Meena et al. (2023); Hews et al. (2010); Lewin et al. (2001). Its corresponding basic reproduction number (representing the average number of new viruses or infected cells generated by a single virus in a naive host) is given by

\begin{matrix} R_{0} = \frac{β p}{c δ} \frac{s}{d} . \end{matrix}

It is easy to show that virus is cleared when

R_{0} < 1

and persists when

R_{0} > 1

(Ciupe and Heffernan 2017). In this study, we expanded model Eq. 1 by taking into account the anatomy of the liver during HBV-induced liver disease. It has been shown that, following hepatitis B infection, liver is restructured into abnormal nodular structures separated by collagen (Fig. 1A and B, right panels). To address differences in virus-host interactions in nodular structures, we developed a two-patch mathematical model and used it to investigate how heterogeneity in infection changes the conditions for viral clearance and viral persistence. Special attention was given to parameters describing liver injury,

δ

, and liver regeneration, s. This two-patch model is a first step in including heterogeneity, with multi-patch models being the obvious next step.

Two-patch model of hepatitis B virus infection. We expanded model Eq. 1 to account for virus dynamics in two liver nodular structures, obtaining a two-patch model of HBV infection. In each nodular structure, we modeled the interaction between uninfected liver cells

T_{j}

, infected liver cells

I_{j}

, and hepatitis B virus

V_{j}

, with

j = {1, 2}

. As before, we assumed that uninfected liver cells are produced at rates

s_{j}

(nodular structure-dependent), die at natural death rate d, and get infected at rate

β

(nodular structure-independent). Infected liver cells are killed at rate

δ

and produce new virus at rate p (nodular structure-independent). Virus is cleared at rate c (nodular structure-independent). Lastly, we assumed that HBV moves between the nodular structures at rates

ϕ_{ij}

(nodular structure-dependent), for

i, j = {1, 2}

and

i \neq j

. The mathematical model becomes:

\begin{matrix} \begin{matrix} \frac{d T_{1}}{dt} & = s_{1} - d T_{1} - β T_{1} V_{1}, \\ \frac{d I_{1}}{dt} & = β T_{1} V_{1} - δ I_{1}, \\ \frac{d V_{1}}{dt} & = p I_{1} - c V_{1} + ϕ_{21} V_{2} - ϕ_{12} V_{1}, \\ \frac{d T_{2}}{dt} & = s_{2} - d T_{2} - β T_{2} V_{2}, \\ \frac{d I_{2}}{dt} & = β T_{2} V_{2} - δ I_{2}, \\ \frac{d V_{2}}{dt} & = p I_{2} - c V_{2} - ϕ_{21} V_{2} + ϕ_{12} V_{1}, \end{matrix} \end{matrix}

2

with initial conditions

T_{1} (0) = s_{1} / d

,

T_{2} (0) = s_{2} / d

,

V_{1} (0) = V_{0}

,

I_{1} (0) = I_{0}

,

I_{2} (0) = V_{2} (0) = 0

. This means that the HBV infection starts in nodular structure 1.

We considered two special cases: (i) the one-directional two-patch model, which assumes that the HBV is seeded in nodular structure 1, moves into nodular structure 2, and never returns into nodular structure 1, i.e. $ϕ_{12} = ϕ$ and $ϕ_{21} = 0$ ; and (ii) the two-directional two-patch model, which assumes that the HBV is seeded in nodular structure 1, and moves at the same rates between the two patches, i.e. $ϕ_{12} = ϕ_{21} = ϕ$ .

Results for the One-Directional Two-Patch Model

The one-directional two-patch model, which assumes that the HBV is seeded in nodular structure 1, moves into nodular structure 2, and never returns into nodular structure 1, is described by:

\begin{matrix} \begin{matrix} \frac{d T_{1}}{dt} & = s_{1} - d T_{1} - β T_{1} V_{1}, \\ \frac{d I_{1}}{dt} & = β T_{1} V_{1} - δ I_{1}, \\ \frac{d V_{1}}{dt} & = p I_{1} - c V_{1} - ϕ V_{1}, \\ \frac{d T_{2}}{dt} & = s_{2} - d T_{2} - β T_{2} V_{2}, \\ \frac{d I_{2}}{dt} & = β T_{2} V_{2} - δ I_{2}, \\ \frac{d V_{2}}{dt} & = p I_{2} - c V_{2} + ϕ V_{1}, \end{matrix} \end{matrix}

3

Asymptotic Analysis Results for the One-Directional Two-Patch Model Eq. 3

We investigated the long-term behavior of the one-directional two-patch model Eq. 3, where HBV is seeded in nodular structure 1, moves into nodular structure 2, and never returns into nodular structure 1. Model Eq. 3 has three non-negative equilibrium solutions. The infection-free equilibrium, representing viral clearance in both patches, is given by:

\begin{matrix} E_{0} = (T_{1}^{0}, I_{1}^{0}, V_{1}^{0}, T_{2}^{0}, I_{2}^{0}, V_{2}^{0}) = (\frac{s_{1}}{d}, 0, 0, \frac{s_{2}}{d}, 0, 0), \end{matrix}

and is always biologically realistic. The one-patch chronic equilibrium, representing virus clearance in nodular structure 1 and persistence in nodular structure 2, is given by:

\begin{matrix} E_{1} = (T_{1}^{1}, I_{1}^{1}, V_{1}^{1}, T_{2}^{1}, I_{2}^{1}, V_{2}^{1}) = (\frac{s_{1}}{d}, 0, 0, \frac{c δ}{β p}, \frac{cd}{β p} (R_{0}^{1D} - 1), \frac{d}{β} (R_{0}^{1D} - 1)), \end{matrix}

where:

\begin{matrix} R_{0}^{1D} = \frac{β p s_{2}}{c d δ} . \end{matrix}

4

Equilibrium

E_{1}

is biologically realistic if and only if

R_{0}^{1D} > 1

. Lastly, model Eq. 3 has two two-patch chronic equilibria,

\begin{matrix} \begin{matrix} E_{2} & = (T_{1}^{2}, I_{1}^{2}, V_{1}^{2}, T_{2}^{2}, I_{2}^{2}, V_{2}^{2}) \\ = (\frac{δ (c + ϕ)}{p β}, \frac{d (c + ϕ)}{p β} (R_{eff}^{1D} - 1), \frac{d}{β} (R_{eff}^{1D} - 1), \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 β d p},) \\ (\frac{c d δ (R_{0}^{1D} - 1) - ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 β p δ}, \frac{c d δ (R_{0}^{1D} - 1) + ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 β c δ}), \end{matrix} \end{matrix}

5

and

\begin{matrix} \begin{matrix} E_{3} & = (T_{1}^{3}, I_{1}^{3}, V_{1}^{3}, T_{2}^{3}, I_{2}^{3}, V_{2}^{3}) \\ = (\frac{δ (c + ϕ)}{p β}, \frac{d (c + ϕ)}{p β} (R_{eff}^{1D} - 1), \frac{d}{β} (R_{eff}^{1D} - 1), \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 β d p}), \\ (\frac{c d δ (R_{0}^{1D} - 1) - ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 β p δ}, \frac{c d δ (R_{0}^{1D} - 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 β c δ}) . \end{matrix} \end{matrix}

6

where

\begin{matrix} R_{eff}^{1D} = \frac{β p s_{1}}{d δ (c + ϕ)}, \end{matrix}

7

and

\begin{matrix} \begin{matrix} D & = d^{2} δ^{2} ({[c (R_{0}^{1 D} - 1) + ϕ (R_{eff}^{1 D} - 1)]}^{2} + 4 c ϕ (R_{eff}^{1 D} - 1)) . \end{matrix} \end{matrix}

Equilibrium

E_{2}

is biologically realistic if and only if

R_{eff}^{1D} > 1

. By contrast,

E_{3}

is not biologically realistic since

I_{2}^{3}

is always negative. Next, we studied the asymptotic stability of the equilibrium solutions

E_{0}

,

E_{1}

and

E_{2}

.

Proposition 1

The infection-free equilibrium

E_{0}

is locally asymptotically stable if

\begin{matrix} max {R_{0}^{1D}, R_{eff}^{1D}} < 1, \end{matrix}

and is unstable otherwise.

Proof

We linearized Eq. 3 at the disease-free equilibrium

E_{0}

,

\begin{matrix} \frac{d}{dt} [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{0}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & 0 \\ 0 & 0 & 0 & - d & 0 & - \frac{β s_{2}}{d} \\ 0 & 0 & 0 & 0 & - δ & \frac{β s_{2}}{d} \\ 0 & 0 & ϕ & 0 & p & - c \end{matrix}] [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{0}) . \end{matrix}

The corresponding Jacobian matrix at the disease-free equilibrium

E_{0}

,

\begin{matrix} J (E_{0}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & 0 \\ 0 & 0 & 0 & - d & 0 & - \frac{β s_{2}}{d} \\ 0 & 0 & 0 & 0 & - δ & \frac{β s_{2}}{d} \\ 0 & 0 & ϕ & 0 & p & - c \end{matrix}], \end{matrix}

can be written as:

\begin{matrix} J (E_{0}) = [\begin{matrix} A_{01} & 0 \\ Φ & A_{02} \end{matrix}], \end{matrix}

where

0

is the

3 \times 3

zero matrix,

\begin{matrix} A_{01} & = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} \\ 0 & - δ & \frac{β s_{1}}{d} \\ 0 & p & - c - ϕ \end{matrix}], A_{02} = [\begin{matrix} - d & 0 & - \frac{β s_{2}}{d} \\ 0 & - δ & \frac{β s_{2}}{d} \\ 0 & p & - c \end{matrix}], Φ = [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & ϕ \end{matrix}] . \end{matrix}

The characteristic equation of model Eq. 3 at

E_{0}

is:

\begin{matrix} 0 = P_{0} (λ) & = det (J (E_{0}) - λ 1_{6}) = det (A_{01} - λ 1_{3}) det (A_{02} - λ 1_{3}) \\ = {(λ + d)}^{2} P_{01} (λ) P_{02} (λ) . \end{matrix}

where

\begin{matrix} P_{01} (λ) & = |\begin{matrix} - δ - λ & \frac{β s_{1}}{d} \\ p & - c - ϕ - λ \end{matrix}| = λ^{2} + (c + δ + ϕ) λ + δ (c + ϕ) (1 - R_{eff}^{1D}), \\ P_{02} (λ) & = |\begin{matrix} - δ - λ & \frac{β s_{2}}{d} \\ p & - c - λ \end{matrix}| = λ^{2} + (c + δ) λ + c δ (1 - R_{0}^{1D}), \end{matrix}

and

1_{k}

is the

k \times k

identity matrix. Two eigenvalues,

λ_{1, 2} = - d

are always negative. It is easy to see that the solutions of

P_{01} (λ)

are negative or have negative real parts if

R_{eff}^{1D} < 1

and the solutions of

P_{02} (λ)

are negative or have negative real parts if

R_{0}^{1D} < 1

. Thus, when

\begin{matrix} max {R_{0}^{1D}, R_{eff}^{1D}} < 1, \end{matrix}

the infection-free equilibrium

E_{0}

is locally asymptotically stable. If, by contrast, either

R_{0}^{1D} > 1

or

R_{eff}^{1D} > 1

, the infection-free equilibrium

E_{0}

is unstable. This completes the proof.

□

Proposition 2

The one-patch chronic equilibrium

E_{1}

exists when

R_{0}^{1D} > 1

, is locally asymptotically stable if

\begin{matrix} R_{eff}^{1D} < 1 < R_{0}^{1D}, \end{matrix}

and is unstable otherwise.

Proof

We linearize model Eq. 3 at equilibrium

E_{1}

,

\begin{matrix} \frac{d}{dt} [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{1}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{β p s_{2}}{c δ} & 0 & - \frac{c δ}{p} \\ 0 & 0 & 0 & - d + \frac{β p s_{2}}{c δ} & - δ & \frac{c δ}{p} \\ 0 & 0 & ϕ & 0 & p & - c \end{matrix}] [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{1}) . \end{matrix}

The corresponding Jacobian matrix at equilibrium

E_{1}

,

\begin{matrix} J (E_{1}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{β p s_{2}}{c δ} & 0 & - \frac{c δ}{p} \\ 0 & 0 & 0 & - d + \frac{β p s_{2}}{c δ} & - δ & \frac{c δ}{p} \\ 0 & 0 & ϕ & 0 & p & - c \end{matrix}], \end{matrix}

can be written as:

\begin{matrix} J (E_{1}) = [\begin{matrix} A_{11} & 0 \\ Φ & A_{12} \end{matrix}], \end{matrix}

where

\begin{matrix} A_{11} = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} \\ 0 & - δ & \frac{β s_{1}}{d} \\ 0 & p & - c - ϕ \end{matrix}], A_{12} = [\begin{matrix} - \frac{β p s_{2}}{c δ} & 0 & - \frac{c δ}{p} \\ - d + \frac{β p s_{2}}{c δ} & - δ & \frac{c δ}{p} \\ 0 & p & - c \end{matrix}], \end{matrix}

and

0

,

Φ

are as before (see Proposition 1.). The characteristic equation of model Eq. 3 at

E_{1}

is:

\begin{matrix} 0 & = P_{1} (λ) = det (J (E_{1}) - λ 1_{6}) = det (A_{11} - λ 1_{3}) det (A_{12} - λ 1_{3}) \\ = (λ + d) P_{11} (λ) P_{12} (λ), \end{matrix}

where

\begin{matrix} P_{11} (λ) & = |\begin{matrix} - δ - λ & \frac{β s_{1}}{d} \\ p & - c - ϕ - λ \end{matrix}| = λ^{2} + (c + δ) λ + c δ (1 - R_{eff}^{1D}), \\ P_{12} (λ) & = |\begin{matrix} - \frac{β p s_{2}}{c δ} - λ & 0 & - \frac{c δ}{p} \\ - d + \frac{β p s_{2}}{c δ} & - δ - λ & \frac{c δ}{p} \\ 0 & p & - c - λ \end{matrix}| \\ = λ^{3} + (c + δ + \frac{β p s_{2}}{c δ}) λ^{2} + (\frac{β p s_{2}}{δ} + \frac{β p s_{2}}{c}) λ \\ + c d δ (R_{0}^{1D} - 1) . \end{matrix}

One eigenvalue

λ = - d

is always negative. It is easy to show that the eigenvalues of

P_{11} (λ)

are negative or have negative real parts if

R_{eff}^{1D} < 1

. For

P_{12} (λ)

, we apply the Routh-Hurwitz conditions which state that the solutions are negative or have negative real part if

b_{i} > 0

for

i = {1, 2, 3}

and

b_{1} b_{2} > b_{3}

, where:

\begin{matrix} P_{12} (λ) = λ^{3} + b_{1} λ^{2} + b_{2} λ + b_{3}, \end{matrix}

and

\begin{matrix} b_{1} & = c + δ + \frac{β p s_{2}}{c δ}, \\ b_{2} & = \frac{β p s_{2}}{δ} + \frac{β p s_{2}}{c}, \\ b_{3} & = c d δ (R_{0}^{1D} - 1) . \end{matrix}

Observe that

b_{1}, b_{2} > 0

always, and

b_{3} > 0

when

R_{0}^{1D} > 1

. Moreover,

\begin{matrix} \begin{matrix} b_{1} b_{2} - b_{3} & = (c + δ + \frac{β p s_{2}}{c δ}) (\frac{β p s_{2}}{δ} + \frac{β p s_{2}}{c}) \\ - c d δ (R_{0}^{1D} - 1) \\ (c + \frac{β p s_{2}}{c δ}) (\frac{β p s_{2}}{δ} + \frac{β p s_{2}}{c}) \\ + \frac{δ β p s_{2}}{c} + c d δ > 0 . \end{matrix} \end{matrix}

8

Therefore, when

\begin{matrix} R_{eff}^{1D} < 1 < R_{0}^{1D}, \end{matrix}

the one-patch chronic equilibrium

E_{1}

is locally asymptotically stable. Conversely, if

R_{eff}^{1D} > 1

or

R_{0}^{1D} < 1

, the one-patch chronic equilibrium

E_{1}

is unstable. This completes the proof.

□

Proposition 3

The biologically realistic two-patch chronic equilibrium

E_{2}

exists and is locally asymptotically stable if

\begin{matrix} R_{eff}^{1D} > 1 . \end{matrix}

Proof

We linearize model Eq. 3 at equilibrium

E_{2}

,

\begin{matrix} \begin{matrix} \frac{d}{dt} [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{2}) \\ = [\begin{matrix} - d - d (R_{eff}^{1D} - 1) & 0 & - \frac{δ (c + ϕ)}{p} & 0 & 0 & 0 \\ d (R_{eff}^{1D} - 1) & - δ & \frac{δ (c + ϕ)}{p} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{d (R_{0}^{1D} + 1)}{2} - \frac{ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 c δ} & 0 & - \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 d p} \\ 0 & 0 & 0 & \frac{d (R_{0}^{1D} - 1)}{2} + \frac{ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 c δ} & - δ & \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 d p} \\ 0 & 0 & ϕ & 0 & p & - c \end{matrix}] \\ \times [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{2}) . \end{matrix} \end{matrix}

The corresponding Jacobian matrix at

E_{2}

is,

\begin{matrix} \begin{matrix} J (E_{2}) & = [\begin{matrix} A_{11} & 0 \\ Φ & A_{22} \end{matrix}] . \end{matrix} \end{matrix}

where

\begin{matrix} A_{11} & = [\begin{matrix} - d - d (R_{eff}^{1D} - 1) & 0 & - \frac{δ (c + ϕ)}{p} \\ d (R_{eff}^{1D} - 1) & - δ & \frac{δ (c + ϕ)}{p} \\ 0 & p & - c - ϕ \end{matrix}], \\ A_{22} & = [\begin{matrix} - \frac{d (R_{0}^{1D} + 1)}{2} - \frac{ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 c δ} & 0 & - \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 d p} \\ \frac{d (R_{0}^{1D} - 1)}{2} + \frac{ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 c δ} & - δ & \frac{c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}}{2 d p} \\ 0 & p & - c \end{matrix}], \end{matrix}

and

0

and

Φ

are as before (see Proposition 1). The characteristic equation of model Eq. 3 at

E_{2}

is:

\begin{matrix} 0 = P_{2} (λ) = det (J (E_{2}) - λ 1_{6}) = det (A_{11} - λ 1_{3}) det (A_{22} - λ 1_{3}) = P_{21} (λ) P_{22} (λ) . \end{matrix}

For both polynomials, we apply the Routh-Hurwitz condition. Let

\begin{matrix} P_{21} (λ) = λ^{3} + a_{1} λ^{2} + a_{2} λ + a_{3}, \end{matrix}

where

\begin{matrix} a_{1} & = c + d + δ + ϕ + d (R_{eff}^{1D} - 1), \\ a_{2} & = d (c + δ + ϕ) + d (c + δ + ϕ) (R_{eff}^{1D} - 1), \\ a_{3} & = δ d (c + ϕ) (R_{eff}^{1D} - 1) . \end{matrix}

Polynomial

P_{21} (λ)

has negative roots or roots with negative real parts when

a_{i} > 0

for all i and

a_{1} a_{2} > a_{3}

. Note that

a_{i} > 0

when

R_{eff}^{1D} > 1

. Moreover, it is easy to show that

a_{1} a_{2} - a_{3} > 0

when

R_{eff}^{1D} > 1

.

Similarly, let

\begin{matrix} P_{22} (λ) = λ^{3} + b_{1} λ^{2} + b_{2} λ + b_{3}, \end{matrix}

where

\begin{matrix} b_{1} & = c + δ + d + β (V_{2}^{2}) \\ = c + δ + \frac{d (R_{0}^{1D} + 1)}{2} + \frac{ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}}{2 c δ}, \\ b_{2} & = (T_{2}^{2}) β p + (V_{2}^{2}) β c + (V_{2}^{2}) β δ + c d + c δ + d δ \\ = \frac{c + δ}{2 c δ} (c d δ (R_{0}^{1D} - 1) + ϕ δ d (R_{eff}^{1D} - 1) + \sqrt{D}) + d δ + c d + c δ \\ + \frac{1}{2 d} (c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1) - \sqrt{D}), \\ b_{3} & = (T_{2}^{2}) β d p + (V_{2}^{2}) β c δ + c d δ \\ = c d δ (R_{0}^{1D} + 1) + ϕ δ d (R_{eff}^{1D} - 1), \end{matrix}

and

V_{2}^{2}

and

T_{2}^{2}

are the target cell and virus equilibrium values for the nodular structure 2 as given by Eq. 5. Polynomial

P_{22} (λ)

has negative roots or roots with negative real parts when

b_{i} > 0

for all i and

b_{1} b_{2} > b_{3}

. Note that when

R_{eff}^{1D} > 1

,

b_{i} > 0

for all i. Moreover,

\begin{matrix} \begin{matrix} b_{1} b_{2} - b_{3} & = β^{2} (c + δ) V_{2}^{2} + T_{2} V_{2} β^{2} p \\ + V_{2} β (c + δ) (c + δ + 2 d) \\ + T_{2} β p (c + δ) \\ + (c + δ) (c + d) (δ + d), \end{matrix} \end{matrix}

is positive when

R_{eff}^{1D} > 1

. Then, by the Routh-Hurwitz condition,

P_{22} (λ)

has negative roots or roots with negative real part when

R_{eff}^{1D} > 1

. This guarantees the local stability analysis of the chronic equilibrium

E_{2}

. Conversely, when

R_{eff}^{1D} < 1

, equilibrum

E_{2}

does not exist. This concludes the proof.

□

Structural Identifiability Results for Model Eq. 3

We investigated the structural identifiability of model Eq. 3 using the differential algebra approach and the DAISY platform (Bellu et al. 2007) (see the Materials and Methods section for further details) under two assumptions: (i) viral load is measured in both patches and the initial conditions of model Eq. 3 are either known or unknown; and (ii) susceptible liver cells are measured in both patches and the initial conditions of model Eq. 3 are either known or unknown.

Proposition 4

When measurements for the $V_{1} (t)$ and $V_{2} (t)$ variables are given, system Eq. 3 is unidentifiable under unknown initial conditions. In particular, parameters $π_{1} = {β, c, d, δ, ϕ}$ are globally identifiable and parameters $π_{2} = {s_{1}, s_{2}, p}$ are unidentifiable. If, additionally, initial conditions are known, then all parameters are globally structurally identifiable.

Proof

The input-output equation for model Eq. 3, given data for

y_{1} = V_{1} (t)

and

y_{2} = V_{2} (t)

, is:

\begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) = \frac{d^{3} y_{1}}{d t^{3}} y_{1} - \frac{d^{2} y_{1}}{d t^{2}} \frac{d y_{1}}{dt} + \frac{d^{2} y_{1}}{d t^{2}} y_{1}^{2} β \\ + \frac{d^{2} y_{1}}{d t^{2}} y_{1} (c + d + δ + ϕ) - {(\frac{d y_{1}}{dt})}^{2} (c + δ + ϕ) \\ + \frac{d y_{1}}{dt} y_{1}^{2} β (c + δ + ϕ) + \frac{d y_{1}}{dt} y_{1} d (c + δ + ϕ) \\ + y_{1}^{3} β δ (c + ϕ) + y_{1}^{2} (- β p s_{1} + c d δ + d δ ϕ) = 0, \\ g (y_{1}, y_{2}, π_{1}, π_{2}) = \frac{d^{2} y_{1}}{d t^{2}} y_{2} ϕ + \frac{d y_{1}}{dt} \frac{d y_{2}}{dt} ϕ - \frac{d y_{1}}{dt} y_{2}^{2} β ϕ \\ - \frac{d y_{1}}{dt} y_{2} ϕ (d + δ) + \frac{d^{3} y_{2}}{d t^{3}} y_{2} - \frac{d^{2} y_{2}}{d t^{2}} \frac{d y_{2}}{dt} \\ + \frac{d^{2} y_{2}}{d t^{2}} y_{2}^{2} β + \frac{d^{2} y_{2}}{d t^{2}} y_{2} (c + d + δ) - {(\frac{d y_{2}}{dt})}^{2} (c + δ) \\ + \frac{d y_{2}}{dt} y_{1} δ ϕ + \frac{d y_{2}}{dt} y_{2}^{2} β (c + δ) \\ + \frac{d y_{2}}{dt} y_{2} d (c + δ) - y_{1} y_{2}^{2} β δ ϕ - y_{1} y_{2} d δ ϕ \\ + y_{2}^{3} β c δ + y_{2}^{2} (- β p s_{2} + c d δ) = 0 . \end{matrix}

Parameters

{π_{1}, π_{2}}

are globally structurally identifiable if:

\begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) = f (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) implies π_{1} = {\hat{π}}_{1} and π_{2} = {\hat{π}}_{2} . \end{matrix}

It is easy to see that:

\begin{matrix} {β = \hat{β}}, {d = \hat{d}}, {δ = \hat{δ}}, {c = \hat{c}}, {ϕ = \hat{ϕ}}, {p \times s_{1} = \hat{p} \times {\hat{s}}_{1}}, {p \times s_{2} = \hat{p} \times {\hat{s}}_{2}} . \end{matrix}

Therefore

{β, c, d, δ, ϕ}

are globally structurally identifiable and

{s_{1}, s_{2}, p}

are unidentifiable. When initial conditions are known, we have:

\begin{matrix} \frac{d V_{1}}{dt} (0) = p I_{1} (0) - c V_{1} (0) - ϕ V_{1} (0) + ϕ V_{2} (0) . \end{matrix}

As long as

I_{1} (0) \neq 0

, p (and by default

s_{1}

and

s_{2}

) are globally structurally identifiable. Hence, model Eq. 3 is globally structurally identifiable for known initial conditions. A summary of the results is given in Table 1.

□

Table 1 Identifiability analysis for model Eq. 3, performed using the DAISY software (Bellu et al. 2007). For the simulations with initial conditions, all initial conditions are known

Observe States	Initial Conditions Known	Initial Conditions Unknown
Model with $V_{1} (t)$ and $V_{2} (t)$ data	Globally structurally identifiable { $β$ , c, d, $δ$ , p, $ϕ$ , $s_{1}$ , $s_{2}$ }	Globally structurally identifiable {c, d, $β$ , $ϕ$ , $δ$ }; Unidentifiable { $s_{1}$ , $s_{2}$ , p}; Correlations { $s_{1} \times p = {\hat{s}}_{1} \times \hat{p}}$ , { $s_{2} \times p = {\hat{s}}_{2} \times \hat{p}$ }
Model with $T_{1} (t)$ and $T_{2} (t)$ data	Globally structurally identifiable { $β$ , c, d, $δ$ , p, $ϕ$ , $s_{1}$ , $s_{2}$ }	Globally structurally identifiable {c, d, $δ$ , $ϕ$ , $s_{1}$ , $s_{2}$ }; Unidentifiable { $β$ , p}; Correlations { $β \times p = \hat{β} \times \hat{p}$ }

Proposition 5

When measurements for the $T_{1} (t)$ and $T_{2} (t)$ variables are given, system Eq. 3 is unidentifiable when initial conditions are unknown. In particular, parameters $π_{1} = {c, d, δ, ϕ, s_{1}, s_{2}}$ are globally structurally identifiable, and parameters $π_{2} = {β, p}$ are unidentifiable. If, additionally, initial conditions are known, then all parameters are globally structurally identifiable.

Proof

The input-output equations for model Eq. 3, given data for

y_{1} = T_{1} (t)

and

y_{2} = T_{2} (t)

, are:

\begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) = - \frac{d^{3} y_{1}}{d t^{3}} y_{1}^{2} y_{2}^{2} + 3 \frac{d^{2} y_{1}}{d t^{2}} \frac{d y_{1}}{dt} y_{1} y_{2}^{2} - \frac{d^{2} y_{1}}{d t^{2}} y_{1}^{2} y_{2}^{2} (c + δ + ϕ) \\ - \frac{d^{2} y_{1}}{d t^{2}} y_{1} y_{2}^{2} s_{1} - 2 {(\frac{d y_{1}}{dt})}^{3} y_{2}^{2} \\ + {(\frac{d y_{1}}{dt})}^{2} y_{1} y_{2}^{2} (c + δ + ϕ) + 2 {(\frac{d y_{1}}{dt})}^{2} y_{2}^{2} s_{1} + \frac{d y_{1}}{dt} y_{1}^{3} y_{2}^{2} β p \\ - \frac{d y_{1}}{dt} y_{1}^{2} y_{2}^{2} δ (c + ϕ) - \frac{d y_{1}}{dt} y_{1} y_{2}^{2} s_{1} (c + δ + ϕ) \\ + \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{3} y_{2} ϕ - {(\frac{d y_{2}}{dt})}^{2} y_{1}^{3} ϕ + \frac{d y_{2}}{dt} y_{1}^{3} y_{2} δ ϕ + \frac{d y_{2}}{dt} y_{1}^{3} ϕ s_{2} + y_{1}^{4} y_{2}^{2} β d p \\ - y_{1}^{3} y_{2}^{2} (β p s_{1} + c d δ) - y_{1}^{3} y_{2} δ ϕ s_{2} \\ + y_{1}^{2} y_{2}^{2} δ s_{1} (c + ϕ) = 0, \\ g (y_{1}, y_{2}, π_{1}, π_{2}) = \frac{d^{2} y_{1}}{d t^{2}} y_{1} y_{2}^{3} ϕ - {(\frac{d y_{1}}{dt})}^{2} y_{2}^{3} ϕ + \frac{d y_{1}}{dt} y_{1} y_{2}^{3} δ ϕ + \frac{d y_{1}}{dt} y_{2}^{3} ϕ s_{1} \\ - \frac{d^{3} y_{2}}{d t^{3}} y_{1}^{2} y_{2}^{2} + 3 \frac{d^{2} y_{2}}{d t^{2}} \frac{d y_{2}}{dt} y_{1}^{2} y_{2} \\ - \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{2} y_{2}^{2} (c + δ) - \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{2} y_{2} s_{2} - 2 {(\frac{d y_{2}}{dt})}^{3} y_{1}^{2} + {(\frac{d y_{2}}{dt})}^{2} y_{1}^{2} y_{2} (c + δ) \\ + 2 {(\frac{d y_{2}}{dt})}^{2} y_{1}^{2} s_{2} \\ + \frac{d y_{2}}{dt} y_{1}^{2} y_{2}^{3} β p - \frac{d y_{2}}{dt} y_{1}^{2} y_{2}^{2} c δ - \frac{d y_{2}}{dt} y_{1}^{2} y_{2} s_{2} (c + δ) + y_{1}^{2} y_{2}^{4} β d p \\ + y_{1}^{2} y_{2}^{3} (- β p s_{2} - c d δ + d δ ϕ) + y_{1}^{2} y_{2}^{2} c δ s_{2} - y_{1} y_{2}^{3} δ ϕ s_{1} = 0 . \end{matrix}

Parameters

{π_{1}, π_{2}}

are globally structurally identifiable if:

\begin{matrix} \begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) & = f (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) and \\ g (y_{1}, y_{2}, π_{1}, π_{2}) & = g (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) \\ imply π_{1} = {\hat{π}}_{1} and π_{2} = {\hat{π}}_{2} . \end{matrix} \end{matrix}

It is easy to see that:

\begin{matrix} {s_{1} = {\hat{s}}_{1}}, {s_{2} = {\hat{s}}_{2}}, {δ = \hat{δ}}, {d = \hat{d}}, {c = \hat{c}}, {ϕ = \hat{ϕ}}, {β \times p = \hat{β} \times \hat{p}} . \end{matrix}

Therefore

{s_{1}, s_{2}, d, δ, c, ϕ}

are globally structurally identifiable and

{β, p}

are unidentifiable. When initial conditions are known, we have:

\begin{matrix} \frac{d T_{1}}{dt} (0) = s_{1} - d T_{1} (0) - β T_{1} (0) V_{1} (0), \end{matrix}

making

β

(and by default p) globally structurally identifiable. Hence, model Eq. 3 is globally structurally identifiable for known initial conditions. A summary of the results is given in Table 1.

□

Numerical Results for Model Eq. 3

As seen above, under known initial conditions and unlimited noise-free measurements of both $V_{1} (t)$ and $V_{2} (t)$ (or $T_{1} (t)$ and $T_{2} (t)$ ), all parameters of model Eq. 3 can be identified from data.

To parametrize the model, we assumed that the entire hepatocyte population $s = d \times (T_{1} (0) + T_{2} (0))$ is susceptible to HBV infection and considered three cases for the susceptible cells ratio within the two patches $s_{1} : s_{2} = {10 : 90, 50 : 50, 90 : 10}$ (see Table 2). To create a baseline model of a chronic infection, we fitted model Eq. 3 to data from an immunosupressed mice (Zhang et al. 2023) (see Materials and Methods). This initial fit is not intended as an exhaustive representation of HBV infection, but as a starting point for model parametrization. To account for the lack of immunity, we set $δ = d = 0.01$ /day (Ciupe, Ribeiro, Nelson, and Perelson 2007) and fixed $c = 4.4$ /day, as in prior studies (Murray et al. 2005). We fitted the remaining parameters ${β, p, ϕ}$ using the fminsearchbnd function in Matlab2021a (see Materials and Methods section for details). A summary of estimates for each case is given in Table 2 and the model dynamics are plotted in Fig. 2.

**Fig. 2** **The effects of patch-susceptibility on virus dynamics within patches when viral movement is one-directional.** Dynamics of **(A)**-**(E)** uninfected (solid line) and infected (dashed line) liver cells in patch 1 (red lines) and patch 2 (blue lines) and **(B)**-**(F)** total virus (black line), virus in patch 1 (red line), virus in patch 2 (blue line) versus mice data (circles) as given by model Eq. 3 in **case 1: (A)**-**(B)**; **case 2: (C)**-**(D)**, **case 3: (E)**-**(F)**. Model parameters are given in Table 2

Model fitting resulted in similar infectivity rates among the three patch recruitment cases, i.e. $β = 3.3 \times 10^{- 9}$ ml/(virion $\times$ day), $β = 2.63 \times 10^{- 9}$ ml/(virion $\times$ day), and $β = 3.13 \times 10^{- 9}$ ml/(virion $\times$ day), for cases 1, 2 and 3, respectively. Viral production rate increased by $21 %$ and $14 %$ in cases 2 and 3 compared to case 1, i.e. $p = 1203$ virion/(ml $\times$ day) and $p = 1137$ virion/(ml $\times$ day), versus $p = 998$ virion/(ml $\times$ day). The biggest difference is in the estimate of the movement rate $ϕ$ , which was 41- and 50-times higher in cases 2 and 3, compared to case 1, $ϕ = 4.1$ /day and $ϕ = 5$ /day, versus $ϕ = 0.1$ /day. To determine if the results can be generalized when we add noise to the data, we performed practical identifiability analysis using Monte Carlo approaches (see Materials and Methods for details). We found that all parameters are strongly practically identifiable for all cases with the exception of parameter $ϕ$ , which is weakly practically identifiable for case 1 and case 2; and parameter p, which is weakly practically identifiable for case 1 (see Table 5 in Materials and Methods).

Table 5 MC approach for the one-directional two-patch model Eq. 3

	case 1			case 2			case 3
ARE	$β$	$ϕ$	p	$β$	$ϕ$	p	$β$	$ϕ$	p
0%	0	0	0	0	0	0	0	0	0
1%	0.3502	1.8354	4.7200	0.0542	12.5492	0.5759	0.0313	0.6557	0.4824
5%	0.4568	24.4316	6.3968	0.0700	12.7504	0.9358	0.1293	1.2374	2.0864
10%	0.5899	36.8748	8.3112	0.1391	13.2641	2.0531	0.2481	2.2095	4.0030
20%	0.7314	45.5693	10.6981	0.3188	13.1036	4.8778	0.4282	3.0876	6.9326
30%	0.9340	57.2279	13.8884	0.6018	13.3719	8.8351	0.7088	3.4158	10.8375
Identifiable?	Yes	Weakly	Weakly	Yes	Weakly	Yes	Yes	Yes	Yes

In all three cases we are in the $R_{0}^{1 D} > 1$ and $R_{eff}^{1 D} > 1$ regime ( $R_{0}^{1 D} = 45.8$ , $R_{eff}^{1 D} = 4.9$ for case 1; $R_{0}^{1 D} = 24.4$ , $R_{eff}^{1 D} = 12.6$ for case 2; and $R_{0}^{1 D} = 5.5$ , $R_{eff}^{1 D} = 23.1$ for case 3), which means virus persists in both patches (see Proposition 5).

The baseline parameters obtained through data fitting provide a static representation of HBV infection, characteristic of chronic infection in an immunocompromised individual. We next generalize and perturb these baseline parameters to account for key biological features, such as immune competence and virus localization. Our goal is to determine the conditions under which the virus is cleared either locally or globally. To acomplish that, we created one-dimensional bifurcation diagrams. We used the formulas for $V_{1}$ and $V_{2}$ in the equilibrium solutions $E_{0}$ , $E_{1}$ , and $E_{2}$ (see Section 3.1 for detail) and varied values of $δ \in (0.01, 0.6)$ . All results are presented as $V_{1}$ and $V_{2}$ at equilibria versus $R_{0}^{1 D}$ (which is inversely proportional to $δ$ ). We found that when the number of susceptible cells in patch 1 is low, $s_{1} : s_{2} = 10 : 90$ (as in case 1), virus is cleared in patch 1 for a large parameter range $1 < R_{0}^{1 D} < 9$ (Fig. 3A) and persists in patch 2 for all $1 < R_{0}^{1 D}$ (Fig. 3B). For equal number of susceptible cells in the two patches, $s_{1} : s_{2} = 50 : 50$ (as in case 2), virus is cleared in patch 1 for a smaller parameter range $1 < R_{0}^{1 D} < 1.7$ (Fig. 3C) and persists in patch 2 for all $1 < R_{0}^{1 D}$ (Fig. 3D). Lastly, when the number of susceptible cells in patch 1 is high, $s_{1} : s_{2} = 90 : 10$ (as in case 3), virus persists in both patches (Fig. 3E and F). As shown analytically, virus is cleared from both patches when $R_{0}^{1 D} < 1$ and $R_{eff}^{1 D} < 1$ (Fig. 3, green solid lines). Interestingly, the long-term results for case 1 and case 2 have $R_{0}^{1 D} < R_{eff}^{1 D}$ (Fig. 3A - D) and long-term results for cases 3 have $R_{0}^{1 D} > R_{eff}^{1 D}$ (Fig. 3E - F). To unify our results, we created two-dimensional bifurcation diagrams (using the formulas for the equilibrium values for $V_{1}$ and $V_{2}$ of the equilibrium solutions $E_{0}$ , $E_{1}$ , and $E_{2}$ , as before) and varied values of $s_{1} \in (0, 500)$ and $s_{2} \in (0, 500)$ . We computed $V_{1}$ , $V_{2}$ and $V = V_{1} + V_{2}$ at equilibria versus $R_{0}^{1 D}$ and $R_{eff}^{1 D}$ values. As shown analytically, we find that the clearance of virus in patch 1 occurs when $R_{eff}^{1 D} < 1$ , regardless of $R_{0}^{1 D}$ value (Fig. 4B), while clearance in patch 2 (and overall) requires both $R_{eff}^{1 D} < 1$ and $R_{0}^{1 D} < 1$ (Fig. 4C and A).

Results for the Two-Directional Two-Patch Model

The two-directional two-patch model, which assumes that the HBV is seeded in nodular structure 1, and moves at the same rates between the two patches can be described by:

\begin{matrix} \begin{matrix} \frac{d T_{1}}{dt} & = s_{1} - d T_{1} - β T_{1} V_{1}, \\ \frac{d I_{1}}{dt} & = β T_{1} V_{1} - δ I_{1}, \\ \frac{d V_{1}}{dt} & = p I_{1} - c V_{1} - ϕ V_{1} + ϕ V_{2}, \\ \frac{d T_{2}}{dt} & = s_{2} - d T_{2} - β T_{2} V_{2}, \\ \frac{d I_{2}}{dt} & = β T_{2} V_{2} - δ I_{2}, \\ \frac{d V_{2}}{dt} & = p I_{2} - c V_{2} + ϕ V_{1} - ϕ V_{2} . \end{matrix} \end{matrix}

9

Asymptotic Analysis Results for the Two-Directional Two-Patch Model Eq. 9

We next investigated the long-term behavior of model Eq. 9, where HBV is seeded in nodular structure 1 and moves at the same rate

ϕ

between the two patches. Model Eq. 9 has at most four equilibria. The infection-free equilibrium:

\begin{matrix} E_{4} = (T_{1}^{4}, I_{1}^{4}, V_{1}^{4}, T_{2}^{4}, I_{2}^{4}, V_{2}^{4}) = (\frac{s_{1}}{d}, 0, 0, \frac{s_{2}}{d}, 0, 0), \end{matrix}

represents complete viral clearance in both patches, as before. The other three potential equilibria are of the form:

\begin{matrix} E_{j} = (T_{1}^{j}, I_{1}^{j}, V_{1}^{j}, T_{2}^{j}, I_{2}^{j}, V_{2}^{j}), \end{matrix}

j \in {5, 6, 7}

, where

\begin{matrix} I_{1}^{j} & = \frac{ϕ T_{1}^{j} (- T_{1}^{j} β^{2} p^{2} s_{2} + T_{1}^{j} β c d δ p + T_{1}^{j} β d δ p ϕ + β c δ p s_{2} + β δ p ϕ s_{2} - c^{2} d δ^{2} - 2 c d δ^{2} ϕ)}{δ (β p T_{1}^{j} - c δ - δ ϕ) (T_{1}^{j} β c p + θ_{j} β p ϕ - c^{2} δ - 2 c δ ϕ)}, \\ V_{1}^{j} & = \frac{ϕ (- T_{1}^{j} β^{2} p^{2} s_{2} + T_{1}^{j} β c d δ p + T_{1}^{j} β d δ p ϕ + β c δ p s_{2} + β δ p ϕ s_{2} - c^{2} d δ^{2} - 2 c d δ^{2} ϕ)}{β (T_{1}^{j} β c p + T_{1}^{j} β p ϕ - c^{2} δ - 2 c δ ϕ) (β p T_{1}^{j} - c δ - δ ϕ)}, \\ T_{2}^{j} & = \frac{δ (T_{1}^{j} β c p + T_{1}^{j} β p ϕ - c^{2} δ - 2 c δ ϕ)}{β p (β p T_{1}^{j} - c δ - δ ϕ)}, \\ I_{2}^{j} & = \frac{T_{1}^{j} β^{2} p^{2} s_{2} - T_{1}^{j} β c d δ p - T_{1}^{j} β d δ p ϕ - β c δ p s_{2} - β δ p ϕ s_{2} + c^{2} d δ^{2} + 2 c d δ^{2} ϕ}{δ β p (β p T_{1}^{j} - c δ - δ ϕ)}, \\ V_{2}^{j} & = \frac{T_{1}^{j} β^{2} p^{2} s_{2} - T_{1}^{j} β c d δ p - T_{1}^{j} β d δ p ϕ - β c δ p s_{2} - β δ p ϕ s_{2} + c^{2} d δ^{2} + 2 c d δ^{2} ϕ}{δ β (T_{1}^{j} β c p + T_{1}^{j} β p ϕ - c^{2} δ - 2 c δ ϕ)} . \end{matrix}

Here

T_{1}^{j}

is the root of

\begin{matrix} γ_{0} {(T_{1}^{j})}^{3} + γ_{1} {(T_{1}^{j})}^{2} + γ_{2} (T_{1}^{j}) + γ_{3} = 0, \end{matrix}

10

where

\begin{matrix} γ_{0} & = β d p^{2} (β c + ϕ), \\ γ_{1} & = - β^{2} c p^{2} s_{1} - β^{2} p^{2} ϕ s_{1} - β^{2} p^{2} ϕ s_{2} - 2 β c^{2} d δ p - 3 β c d δ p ϕ, \\ γ_{2} & = 2 β c^{2} δ p s_{1} + 4 β c δ p ϕ s_{1} + β c δ p ϕ s_{2} + β δ p ϕ^{2} s_{1} + β δ p ϕ^{2} s_{2} + c^{3} d δ^{2} + 2 c^{2} d δ^{2} ϕ, \\ γ_{3} & = - c δ^{2} s_{1} (c^{2} + 3 c ϕ + 2 ϕ^{2}) . \end{matrix}

Since

γ_{0} > 0

,

γ_{1} < 0

,

γ_{2} > 0

and

γ_{3} < 0

, by Descartes’ rule of signs, we can have one or three

T_{1}^{j}

positive roots. Additionally, for the chronic equilibrium

E_{j}

to exist, we need to guarantee positivity of the other variables in equilibrium

E_{j}

. This occurs when:

\begin{matrix} T_{1}^{j} < min {\frac{δ (c + ϕ)}{β p}, \frac{c δ (c + 2 ϕ)}{β p (c + ϕ)}, \frac{δ}{β p} \frac{β p s_{2} (c + ϕ) - c d δ (c + 2 ϕ)}{β p s_{2} - d δ (c + ϕ)}} . \end{matrix}

11

Proposition 6

The infection-free equilibrium

E_{4}

is locally asymptotically stable if

\begin{matrix} R_{0}^{2 D} = \frac{β p s_{2}}{d δ} \frac{p s_{1} β + d δ (c + ϕ)}{p s_{1} β (c + ϕ) + d c δ (c + 2 ϕ)} < 1, \end{matrix}

and unstable otherwise.

Proof

We linearize Eq. 9 at the infection-free equilibrium

E_{4}

:

\begin{matrix} \frac{d}{dt} [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{4}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & ϕ \\ 0 & 0 & 0 & - d & 0 & - \frac{β s_{2}}{d} \\ 0 & 0 & 0 & 0 & - δ & \frac{β s_{2}}{d} \\ 0 & 0 & ϕ & 0 & p & - c - ϕ \end{matrix}] [\begin{matrix} T_{1} \\ I_{1} \\ V_{1} \\ T_{2} \\ I_{2} \\ V_{2} \end{matrix}] (E_{4}), \end{matrix}

The corresponding Jacobian matrix at

E_{4}

is:

\begin{matrix} J (E_{4}) = [\begin{matrix} - d & 0 & - \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & - δ & \frac{β s_{1}}{d} & 0 & 0 & 0 \\ 0 & p & - c - ϕ & 0 & 0 & ϕ \\ 0 & 0 & 0 & - d & 0 & - \frac{β s_{2}}{d} \\ 0 & 0 & 0 & 0 & - δ & \frac{β s_{2}}{d} \\ 0 & 0 & ϕ & 0 & p & - c - ϕ \end{matrix}], \end{matrix}

and the characteristic equation for model Eq. 9 at

E_{4}

is:

\begin{matrix} 0 = P_{4} (E_{4}) = det (J (E_{4}) - λ 1_{6}) = {(λ + d)}^{2} |\begin{matrix} - δ - λ & - \frac{β s_{1}}{d} & 0 & 0 \\ p & - c - ϕ - λ & 0 & ϕ \\ 0 & 0 & - δ - λ & \frac{β s_{2}}{d} \\ 0 & ϕ & p & - c - ϕ - λ \end{matrix}| . \end{matrix}

Two eigenvalues,

λ_{1, 2} = - d

are always negative. The remaining polynomial is:

\begin{matrix} P_{41} (λ) = λ^{4} + a_{1} λ^{3} + a_{2} λ^{2} + a_{3} λ + a_{4}, \end{matrix}

where

\begin{matrix} a_{1} & = 2 c + 2 δ + 2 ϕ, \\ a_{2} & = {(c + δ)}^{2} + 2 ϕ (c + 2 δ) + \frac{c^{2} δ}{c + ϕ} + \frac{β^{2} p^{2} s 1 s 2}{d^{2} δ (c + ϕ)} \\ + \frac{β p s_{1} (c + ϕ) + d c δ (c + 2 ϕ)}{d (c + ϕ)} (1 - R_{0}^{2 D}), \\ a_{3} & = 2 (c + δ) (c + ϕ) δ + 2 c ϕ δ + \frac{(c + ϕ + δ) β^{2} p^{2} s_{1} s_{2}}{d^{2} δ (c + ϕ)} \\ + \frac{β p s_{1} (c + ϕ) + d c δ (c + 2 ϕ)}{d (c + ϕ)} (c + ϕ + δ) (1 - R_{0}^{2 D}), \\ a_{4} & = \frac{β p s_{1} δ (c + ϕ) + d c δ^{2} (c + 2 ϕ)}{d} (1 - R_{0}^{2 D}) . \end{matrix}

By the Routh-Hurwitz criteria, polynomial

P_{41} (λ)

has negative roots or roots with negative real parts when

a_{i} > 0

and

a_{1} a_{2} a_{3} > a_{3}^{2} + a_{1}^{2} a_{4}

. If

R_{0}^{2 D} < 1

, then all

a_{i} > 0

. Moreover, using Maple, we can show that

a_{1} a_{2} a_{3} > a_{3}^{2} + a_{1}^{2} a_{4}

when

R_{0}^{2 D} < 1

(see Appendix A). Thus, by the Routh-Hurwitz criteria, if

R_{0}^{2 D} < 1

, the infection-free equilibrium

E_{4}

is locally asymptotically stable. In contrast, if

R_{0}^{2 D} > 1

, the infection-free equilibrium

E_{4}

is unstable. This concludes the proof.

□

Lastly, since explicit forms of chronic equilibria $E_{j}$ ( $j = 5, 6, 7$ ) are hard to determine, we examined their existence and asymptotic stability numerically.

Structural Identifiability Results for Model Eq. 9

As with the one-directional two-patch model, we investigated the structural identifiability of model Eq. 9 under two assumptions: (i) viral load is measured in both patches and the initial conditions of the model Eq. 9 are either known or unknown; (ii) susceptible liver cells are measured in both patches and the initial conditions of the model Eq. 9 are either known or unknown.

Proposition 7

Given measurements for the $V_{1} (t)$ and $V_{2} (t)$ variables, system Eq. 9 is unidentifiable when initial conditions are unknown. In particular, parameters $π_{1} = {β, c, d, δ,}$ are globally structurally identifiable, and parameters $π_{2} = {s_{1}, s_{2}, p}$ are unidentifiable. If, additionally, initial conditions are known, then all parameters are globally structurally identifiable.

Proof

The input-output equations for model Eq. 9, given data for

y_{1} = V_{1} (t)

and

y_{2} = V_{2} (t)

, are:

\begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) = \frac{d^{3} y_{1}}{d t^{3}} y_{1} - \frac{d^{2} y_{1}}{d t^{2}} \frac{d y_{1}}{dt} + \frac{d^{2} y_{1}}{d t^{2}} y_{1}^{2} β + \frac{d^{2} y_{1}}{d t^{2}} y_{1} (c + d + δ + ϕ) \\ - {(\frac{d y_{1}}{dt})}^{2} (c + δ + ϕ) + \frac{d y_{1}}{dt} \frac{d y_{2}}{dt} ϕ + \frac{d y_{1}}{dt} y_{1}^{2} β (c + δ + ϕ) \\ + \frac{d y_{1}}{dt} y_{1} d (c + δ + ϕ) + \frac{d y_{1}}{dt} y_{2} δ ϕ - \frac{d^{2} y_{2}}{d t^{2}} y_{1} ϕ \\ - \frac{d y_{2}}{dt} y_{1}^{2} β ϕ - \frac{d y_{2}}{dt} y_{1} ϕ (d + δ) + y_{1}^{3} β δ (c + ϕ) \\ - y_{1}^{2} y_{2} β δ ϕ + y_{1}^{2} (- β p s_{1} + c d δ + d δ ϕ) \\ - y_{1} y_{2} d δ ϕ = 0, \\ g (y_{1}, y_{2}, π_{1}, π_{2}) = - \frac{d^{2} y_{1}}{d t^{2}} y_{2} ϕ + \frac{d y_{1}}{dt} \frac{d y_{2}}{dt} ϕ \\ - \frac{d y_{1}}{dt} y_{2}^{2} β ϕ - \frac{d y_{1}}{dt} y_{2} ϕ (d + δ) + \frac{d^{3} y_{2}}{d t^{3}} y_{2} \\ - \frac{d^{2} y_{2}}{d t^{2}} \frac{d y_{2}}{dt} + \frac{d^{2} y_{2}}{d t^{2}} y_{2}^{2} β + \frac{d^{2} y_{2}}{d t^{2}} y_{2} (c + d + δ + ϕ) \\ - {(\frac{d y_{2}}{dt})}^{2} (c + δ + ϕ) + \frac{d y_{2}}{dt} y_{1} δ ϕ \\ + \frac{d y_{2}}{dt} y_{2}^{2} β (c + δ + ϕ) + \frac{d y_{2}}{dt} y_{2} d (c + δ + ϕ) \\ - y_{1} y_{2}^{2} β δ ϕ - y_{1} y_{2} d δ ϕ + y_{2}^{3} β δ (c + ϕ) \\ + y_{2}^{2} (- β p s_{2} + c d δ + d δ ϕ) = 0 . \end{matrix}

Parameters

{π_{1}, π_{2}}

are globally structurally identifiable if:

\begin{matrix} \begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) & = f (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) and \\ g (y_{1}, y_{2}, π_{1}, π_{2}) & = g (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) \\ imply π_{1} = {\hat{π}}_{1} and π_{2} = {\hat{π}}_{2} . \end{matrix} \end{matrix}

It is easy to see that:

\begin{matrix} {β = \hat{β}}, {d = \hat{d}}, {δ = \hat{δ}}, {c = \hat{c}}, {ϕ = \hat{ϕ}}, {p \times s_{1} = \hat{p} \times {\hat{s}}_{1}}, {p \times s_{2} = \hat{p} \times {\hat{s}}_{2}} . \end{matrix}

Therefore

{β, c, d, δ, ϕ}

are globally structurally identifiable and

{s_{1}, s_{2}, p}

are unidentifiable. When initial conditions are known, we have

\begin{matrix} \frac{d V_{1}}{dt} (0) = p I_{1} (0) - c V_{1} (0) - ϕ V_{1} (0) + ϕ V_{2} (0) . \end{matrix}

As long as

I_{1} (0) \neq 0

, p (and by default

s_{1}

and

s_{2}

) are globally structurally identifiable. Hence, model Eq. 9 is globally structurally identifiable for known initial conditions. A summary of the results is given in Table 3.

□

Table 3 Identifiability analysis for model Eq. 9, performed using the DAISY software (Bellu et al. 2007). For the simulations with initial conditions, all initial conditions are known

Observe States	Initial Conditions Known	Initial Conditions Unknown
Model with $V_{1} (t)$ and $V_{2} (t)$ data	Globally structurally identifiable { $β$ , c, d, $δ$ , $ϕ$ , p, $s_{1}$ , $s_{2}$ }	Globally structurally identifiable { $β$ , c, d, $δ$ , $ϕ$ }; Unidentifiable { $s_{1}$ , $s_{2}$ , p}; Correlations ${s_{1} \times p = {\hat{s}}_{1} \times \hat{p}}$ , ${s_{2} \times p = {\hat{s}}_{2} \times \hat{p}}$
Model with $T_{1} (t)$ and $T_{2} (t)$ data	Globally structurally identifiable { $β, c, d, δ, ϕ, p, s_{1}, s_{2}$ }	Globally structurally identifiable { $c, d, δ, ϕ, s_{1}, s_{2}$ }; Unidentifiable { $β, p$ }; Correlations ${β \times p = \bar{β} \times \bar{p}}$

Proposition 8

Given measurements for the $T_{1} (t)$ and $T_{2} (t)$ variables, system Eq. 9 is unidentifiable when initial conditions are unknown. In particular, parameters $π_{1} = {c, d, δ, ϕ, s_{1}, s_{2}}$ are globally structurally identifiable, and parameters $π_{2} = {β, p}$ are unidentifiable. If, additionally, initial conditions are known, then all parameters are globally structurally identifiable.

Proof

The input-output equations for model Eq. 9, given data for

y_{1} = T_{1} (t)

and

y_{2} = T_{2} (t)

, are:

\begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) = - \frac{d^{3} y_{1}}{d t^{3}} y_{1}^{2} y_{2}^{2} + 3 \frac{d^{2} y_{1}}{d t^{2}} \frac{d y_{1}}{dt} y_{1} y_{2}^{2} - \frac{d^{2} y_{1}}{d t^{2}} y_{1}^{2} y_{2}^{2} (c + δ + ϕ) \\ - \frac{d^{2} y_{1}}{d t^{2}} y_{1} y_{2}^{2} s_{1} - 2 {(\frac{d y_{1}}{dt})}^{3} y_{2}^{2} \\ + {\frac{d y_{1}}{dt}}^{2} y_{1} y_{2}^{2} (c + δ + ϕ) + 2 {\frac{d y_{1}}{dt}}^{2} y_{2}^{2} s_{1} + \frac{d y_{1}}{dt} y_{1}^{3} y_{2}^{2} β p \\ - \frac{d y_{1}}{dt} y_{1}^{2} y_{2}^{2} δ (c + ϕ) - \frac{d y_{1}}{dt} y_{1} y_{2}^{2} s_{1} (c + δ + ϕ) \\ + \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{3} y_{2} ϕ - {(\frac{d y_{2}}{dt})}^{2} y_{1}^{3} ϕ + \frac{d y_{2}}{dt} y_{1}^{3} y_{2} δ ϕ \\ + \frac{d y_{2}}{dt} y_{1}^{3} ϕ s_{2} + y_{1}^{4} y_{2}^{2} β d p - y_{1}^{3} y_{2}^{2} (β p s_{1} + c d δ) \\ - y_{1}^{3} y_{2} δ ϕ s_{2} + y_{1}^{2} y_{2}^{2} δ s_{1} (c + ϕ) = 0, \\ g (y_{1}, y_{2}, π_{1}, π_{2}) = \frac{d^{2} y_{1}}{d t^{2}} y_{1} y_{2}^{3} ϕ - {(\frac{d y_{1}}{dt})}^{2} y_{2}^{3} ϕ \\ + \frac{d y_{1}}{dt} y_{1} y_{2}^{3} δ ϕ + \frac{d y_{1}}{dt} y_{2}^{3} ϕ s_{1} \\ - \frac{d^{3} y_{2}}{d t^{3}} y_{1}^{2} y_{2}^{2} + 3 \frac{d^{2} y_{2}}{d t^{2}} \frac{d y_{2}}{dt} y_{1}^{2} y_{2} \\ - \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{2} y_{2}^{2} (c + δ + ϕ) - \frac{d^{2} y_{2}}{d t^{2}} y_{1}^{2} y_{2} s_{2} \\ - 2 {(\frac{d y_{2}}{dt})}^{3} y_{1}^{2} + {(\frac{d y_{2}}{dt})}^{2} y_{1}^{2} y_{2} (c + δ + ϕ) \\ + 2 {(\frac{d y_{2}}{dt})}^{2} y_{1}^{2} s_{2} + \frac{d y_{2}}{dt} y_{1}^{2} y_{2}^{3} β p \\ - \frac{d y_{2}}{dt} y_{1}^{2} y_{2}^{2} δ (c + ϕ) \\ - \frac{d y_{2}}{dt} y_{1}^{2} y_{2} s_{2} (c + δ + ϕ) \\ + y_{1}^{2} y_{2}^{4} β d p - y_{1}^{2} y_{2}^{3} (β p s_{2} + c d δ) \\ + y_{1}^{2} y_{2}^{2} δ s_{2} (c + ϕ) - y_{1} y_{2}^{3} δ ϕ s_{1} = 0 . \end{matrix}

Parameters

{π_{1}, π_{2}}

are globally structurally identifiable if:

\begin{matrix} \begin{matrix} f (y_{1}, y_{2}, π_{1}, π_{2}) & = f (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) and \\ g (y_{1}, y_{2}, π_{1}, π_{2}) & = g (y_{1}, y_{2}, {\hat{π}}_{1}, {\hat{π}}_{2}) \\ imply π_{1} = {\hat{π}}_{1} and π_{2} = {\hat{π}}_{2} . \end{matrix} \end{matrix}

It is easy to see that:

\begin{matrix} {s_{1} = {\hat{s}}_{1}}, {s_{2} = {\hat{s}}_{2}}, {δ = \hat{δ}}, {d = \hat{d}}, {c = \hat{c}}, {ϕ = \hat{ϕ}}, {β \times p = \hat{β} \times \hat{p}} . \end{matrix}

Therefore

{s_{1}, s_{2}, d, δ, c, ϕ}

are globally structurally identifiable and

{β, p}

are unidentifiable. When initial conditions are known, we have:

\begin{matrix} \frac{d T_{1}}{dt} (0) = s_{1} - d T_{1} (0) - β T_{1} (0) V_{1} (0), \end{matrix}

making

β

(and by default p) globally structurally identifiable. Hence, model Eq. 9 is globally structurally identifiable for known initial conditions. A summary of the results is given in Table 3.

□

Numerical Results for Model Eq. 9

As seen above, under known initial conditions and unlimited noise-free measurements of both $V_{1} (t)$ and $V_{2} (t)$ (or $T_{1} (t)$ and $T_{2} (t)$ ), all parameters of model Eq. 9 can be identified from data.

As before, we assumed that the entire hepatocyte population $s = d \times (T_{1} (0) + T_{2} (0))$ is susceptible to HBV infection and considered three cases for the susceptible cells ratio within the two patches $s_{1} : s_{2} = {10 : 90, 50 : 50, 90 : 10}$ (see Table 4). As before, we set $δ = d = 0.01$ /day (Ciupe, Ribeiro, Nelson, and Perelson 2007), fixed $c = 4.4$ /day and fitted parameters ${β, p, ϕ}$ using the fminsearchbnd function in Matlab2021a (see Materials and Methods section for details). A summary of estimates for each case is given in Table 4 and the model dynamics are plotted in Fig. 5.

**Fig. 5** **The effects of patch-susceptibility on virus dynamics within patches when viral movement is two-directional.** Dynamics of **(A)**-**(E)** uninfected (solid line) and infected (dashed line) liver cells in patch 1 (red lines) and patch 2 (blue lines) and **(B)**-**(F)** total virus (black line), virus in patch 1 (red line), virus in patch 2 (blue line) versus mice data (circles) as given by model Eq. 9 in **case 1: (A)**-**(B)**; **case 2: (C)**-**(D)**, **case 3: (E)**-**(F)**. Model parameters are given in Table 4

As with the one-directional two-patch model, we obtain similar infectivity rates among the three patch recruitment cases, i.e. $β = 2.93 \times 10^{- 9}$ ml/(virion $\times$ day), $β = 2.96 \times 10^{- 9}$ ml/(virion $\times$ day), and $β = 2.94 \times 10^{- 9}$ ml/(virion $\times$ day), for cases 1, 2 and 3, respectively. Moreover, viral production rates are similar among the three patch recruitment cases, i.e. $p = 1053$ virion/(ml $\times$ day), $p = 1055$ virion/(ml $\times$ day), and $p = 1049$ virion/(ml $\times$ day), respectively. The estimate of the movement rate is 50-times higher in case 1 and case 3, compared to case 2, i.e. $ϕ = 5$ /day versus $ϕ = 0.1$ /day. To determine if the results can be generalized when we add noise to the data, we performed practical identifiability analysis using Monte Carlo approaches (see Materials and Methods for details). We found that all parameters are strongly practically identifiable for all cases, with the exception of $ϕ$ , which is weakly practically identifiable in all cases (see Table 6).

Table 6 MC approach for the two-directional two-patch model Eq. 9.

	case 1			case 2			case 3
ARE	$β$	$ϕ$	p	$β$	$ϕ$	p	$β$	$ϕ$	p
0%	0	0	0	0	0	0	0	0	0
1%	0.0207	1.8487	0.3372	0.0242	10.3416	0.3967	0.0366	2.5045	0.5624
5%	0.2271	3.8460	3.2488	0.1183	76.3263	1.9011	0.3342	1.4423	4.6259
10%	0.4226	2.3097	5.9290	0.2775	64.4378	4.2071	0.4965	1.5473	6.8888
20%	0.6925	1.4294	9.7573	0.6096	51.1635	8.8355	0.7205	1.6928	10.0169
30%	0.8235	1.4597	11.9580	0.8160	36.3667	11.7393	0.8739	1.5403	12.5380
Identifiable?	Yes	Weakly	Yes	Yes	Weakly	Yes	Yes	Weakly	Yes

In all three cases we are in the $R_{0}^{2 D} > 1$ regime ( $R_{0}^{2 D} = 22$ , for case 1; $R_{0}^{2 D} = 23.6$ , for case 2; and $R_{0}^{2 D} = 2.3$ , for case 3). While we know from the analytical results that this corresponds to instability of the clearance equilibrium $E_{4}$ , we do not know if it guarantees the stability of the both-patches viral persistence equilibrium $E_{5}$ (note that equilibrium solutions for single-patch viral persistence do not exist for this model). We investigated these long-term results numerically by deriving one- and two-dimensional bifurcation diagrams.

As with the one-directional two-patch model, we perturb the baseline parameters obtained through data fitting to uncover general biological features, such as immune competence and virus localization. To address this, we created one-dimensional bifurcation diagrams as follows. We numerically solved system Eq. 9 for parameters in case 2 (Table 4) and

δ \in (0.01, 0.6)

using the ode15 function in Matlab2021a. We ran the model for

t_{end} = 3000

days and plotted

V_{1} (t_{end})

and

V_{2} (t_{end})

versus

R_{0}^{2 D}

. We found that both

V_{1}

and

V_{2}

are cleared when

R_{0}^{2 D} < 1

and persist when

R_{0}^{2 D} > 1

(Fig. 6A and B). We also created two-dimensional bifurcation diagrams (using the numerical values of

V_{1} (t_{end})

and

V_{2} (t_{end})

, as before), for varied values of

δ \in (0.01, 0.6)

, and

ϕ \in (0, 10)

. We found little variability in the long-term dynamics of

V_{1}

and

V_{2}

, with both viruses being cleared for high infected cells killing rate

δ

and persisting for low infected cells killing rate

δ

, regardless of the movement parameter

ϕ

(Fig. 7). This implies that, when virus moved between patches, immune mediated removal of infected cells is mandatory for viral clearance.

**Fig. 6** One-dimensional bifurcation diagram for model Eq. 9: ${log}_{10} V_{1}$ and ${log}_{10} V_{2}$ at equilibrium versus $R_{0}^{2 D}$ . The other parameters are given in Table 2 for **case 2**

Table 4 Parameter values for model Eq. 9

Fixed Parameters all cases	Description	Value	Reference
c	Virus clearance rate	4.4 /day	Murray et al. (2005)
d	Uninfected hepatocyte death rate	0.01 /day	M.-J. Wang et al. (2017)
$δ$	Infected hepatocyte death rate	0.01 /day
Fixed Parameters case 1	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$6.8 \times 10^{2}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	6.12 $\times 10^{3}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 1	Description	Value	RSS
$β$	Infectivity rate	$2.93 \times 10^{- 9}$ ml/(virion $\times$ day)	0.81
p	Virus production rate	1053 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	5 /day	-
Fixed Parameters case 2	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$3.4 \times 10^{3}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	$3.4 \times 10^{3}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 2	Description	Value	RSS
$β$	Infectivity rate	$2.96 \times 10^{- 9}$ ml/(virion $\times$ day)	0.76
p	Virus production rate	1055 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	0.1 /day	-
Fixed Parameters case 3	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$6.12 \times 10^{3}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	$6.8 \times 10^{2}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 3	Description	Value	RSS
$β$	Infectivity rate	$2.94 \times 10^{- 9}$ ml/(virion $\times$ day)	0.81
p	Virus production rate	1049 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	5 /day	-
Initial conditions	Description	Value	Reference
$T_{1} (0)$	Uninfected hepatocytes patch 1	$3.4 \times 10^{5}$ cells/ml	Ciupe et al. (2024)
$T_{2} (0)$	Uninfected hepatocytes patch 2	$3.4 \times 10^{5}$ cells/ml	Ciupe et al. (2024)
$I_{1} (0)$	Infected hepatocytes patch 1	1 cell/ml	-
$I_{2} (0)$	Infected hepatocytes patch 2	0 cell/ml	-
$V_{1} (0)$	Inoculum patch 1	$10^{4}$ HBV DNA/ml	Zhang et al. (2023)
$V_{2} (0)$	Inoculum patch 3	0 HBV DNA/ml	-

Discussion and Conclusion

In this study, we developed two within-host mathematical models of HBV infection that take into account the abnormal nodular structure of the liver following prolonged hepatitis B infection and disease (Yuan et al. 2019; Bataller et al. 2005; Roehlen et al. 2020). For simplicity, we only modeled a two-patch liver structure and assumed that: (1) each patch has a different numbers of cells susceptible to HBV infection, (2) HBV is seeded in a single-patch, (3) all other virus-host interactions are patch-agnostic. The first model, given by Eq. 3, assumes one-directional hepatitis B virus movement between patches. The second model, given by Eq. 9, assumes two-directional hepatitis B virus movement between patches. We used the models to determine differences in overall hepatitis B virus infection under viral localization (as given by the one-directional model Eq. 3) and free viral movement between patches (as given by the two-directional model Eq. 9).

We analyzed both models using asymptotic stability techniques and found that when HBV movement is one-directional, as given by model Eq. 3, equilibrium dynamics include viral clearance from both patches, viral clearance from the HBV seeding patch alone, and virus persistence in both patches. When HBV movement is two-directional, as given by model Eq. 9, however, equilibrium dynamics only include viral clearance from both patches, and virus persistence in both patches. In other words, single-patch persistence and viral containment within parts of the liver only happen when HBV movement between patches is irreversible. This can be a result of collagen production, excessive accumulation of extracellular matrix, and liver dysfunction (Bataller et al. 2005).

Before validating the models with data, we studied their structural identifiability, and found that we can identify all parameters (under unlimited data) only when viral inoculum and initial susceptible hepatocyte populations are known. This is the case regardless of the type of data used (liver cell or virus data). When these initial details are unknown, several parameters are unidentifiable. A common unidentifiable parameter between the two models Eqs. 3 and 9 is the viral production rate p, which, when known (or fixed), ensures identifiability of the remaining parameters. Since p is hard to estimate experimentally, we assumed known the initial condition of both models, and estimated p from the data.

We fitted the models to HBV DNA data from one immunocompromised, humanized mice (Zhang et al. 2023) and estimated three key parameter values for each model: the infectivity rate $β$ , the viral production rate p, and the HBV movement rate between patches $ϕ$ . Since empirical data (from an immunocompromised animal) showed viral persistence following the initial expansion, we ignored the effect of immune-mediated killing rate of infected cells $δ$ (which was set to uninfected cell death rate) and relaxed this assumption later on. We considered differences in susceptibility to HBV infection between patches, by varying the target cell recruitment ratio, $s_{1} : s_{2}$ between 10:90, 50:50 and 90:10, while keeping the overall liver cell numbers fixed (a reasonable assumption in an immunocompromised mice (Wieland and Chisari 2005)).

We found that the estimates for the infectivity rate $β$ and for the viral production rate p are not affected (in a significant manner) by the target cell recruitment ratio $s_{1} : s_{2}$ or by the ability of HBV to move back-and-forth between patches ( $ϕ_{21} = 0$ or $ϕ_{21} = ϕ_{12} = ϕ$ ). The movement rate $ϕ$ , however, is significantly impacted by the target cell recruitment ratio $s_{1} : s_{2}$ and by the model. In particular, for the one-directional two-patch movement model Eq. 3, HBV chooses movement towards the non-founding patch at significantly higher rates (up to 50-times higher) when the initial patch has more susceptible cells. This outcome may be due to larger numbers of HBV being created in the original patch in the 50:50 and 90:10 cases, resulting in larger virus movement. By contrast, for the two-directional two-patch movement model Eq. 9, HBV chooses higher movement rates (up to 50-times higher) when the patches have unbalanced numbers of susceptible cells (10:90 and 90:10, respectively) and limited movement for the balanced patches (50:50). This outcome may be due to each patch having a sufficient number of target cells available for HBV infection, limiting viral diffusion.

For the one-directional model Eq. 3, the equilibrium virus levels in patch 2 exceed equilibrium virus levels in patch 1, regardless of the target cell susceptibility ratio (Fig. 2). By contrast, for the two-directional model Eq. 9, equilibrium virus levels in patch 2 exceed those in patch 1 for the 10 : 90 target cell susceptibility, but are equal to and lower than those in patch 1 for the 50 : 50 and 90 : 10 target cell susceptibility ratio (Fig. 5). Therefore, the quantitative outcomes are influenced by the ability of the virus to move one- or two-directionally between patches. Our results are not influenced by practical identifiability issues, with both models being (at least weakly) practically identifiable.

The baseline parameters obtained through data fitting in both models Eq. 3 and Eq. 9 offer a static view of HBV infection, characteristic of chronicity in an immunocompromised individual. To determine how the results change in an immune competent host, we derived one- and two-dimensional bifurcation diagrams by using $δ$ (the killing rate of infected cells) as bifurcation parameter. We varied $δ$ from 0.01 /day (corresponding to infected cells lifespan of 100 days) to 0.6 /day (corresponding to infected cells lifespan of 1.7 days).

For model Eq. 3 we found, as expected, that increasing $δ$ results in $R_{0}^{1 D}$ decrease. When $R_{0}^{1 D}$ decreases below 1, we observe viral clearance from both patches (Fig. 3, solid green lines). Interestingly, viral clearance from both patches requires that $δ > 0.45$ /day (lifespan of infected cells of 2.2 days or lower) when $s_{1} : s_{2} = 10 : 90$ , $δ > 0.24$ /day (lifespan of infected cells of 4 days or lower) when $s_{1} : s_{2} = 50 : 50$ , and $δ > 0.055$ /day (lifespan of infected cells of 18 days or lower) when $s_{1} : s_{2} = 90 : 10$ . This suggests that higher immune responses are necessary for viral clearance when the majority of susceptible cells are in patch 2. By contrast, persistence in a single-patch outcomes (e. g. viral clearance in patch 1, and persistence in patch 2) occur for $δ \in (0.05, 0.45)$ when $s_{1} : s_{2} = 10 : 90$ , $δ \in (0.12, 0.24)$ when $s_{1} : s_{2} = 50 : 50$ , and not at all (in the ranges considered) when $s_{1} : s_{2} = 90 : 10$ (Fig. 3, solid red lines and Fig. 4). This suggests a limited range of immune responses lead to virus containment inside a single-patch when the majority of susceptible cells are in patch 1.

**Fig. 3** One-dimensional bifurcation diagrams for model Eq. 3: $V_{1}$ and $V_{2}$ at equilibrium versus $R_{0}^{1 D}$ . The other parameters are given in Table 2 for **A - B** **case 1**, **C - D** **case 2** and **E - F** **case 3**

**Fig. 4** Two-dimensional bifurcation diagrams for model Eq. 3: A. ${log}_{10} (V_{1} + V_{2})$ at equilibrium; B. ${log}_{10} V_{1}$ at equilibrium; C. ${log}_{10} V_{2}$ at equilibrium versus $R_{eff}^{1 D}$ and $R_{0}^{1 D}$ . The other parameters are given in Table 2, case 2

For model Eq. 9 we found that viral clearance from both patches requires that $δ > 0.27$ /day (lifespan of infected cells of 3.7 days or lower) when $s_{1} : s_{2} = 10 : 90$ ; $δ > 0.234$ /day (lifespan of infected cells of 4.27 days or lower) when $s_{1} : s_{2} = 50 : 50$ ; and $δ > 0.023$ /day (lifespan of infected cells of 43 days or lower) when $s_{1} : s_{2} = 90 : 10$ . This suggests that higher immune responses are necessary for viral clearance when the susceptible cells in patch 1 equal or exceed the susceptible cells in patch 2 (Fig. 7).

**Fig. 7** Two-dimensional bifurcation diagram for model Eq. 9: A. ${log}_{10} (V_{1} + V_{2})$ at equilibrium; B. ${log}_{10} V_{1}$ at equilibrium; C. ${log}_{10} V_{2}$ at equilibrium versus $δ$ and $ϕ$ . The other parameters are as in Table 4, case 2

Our study has several limitations. First, the models ignore liver proliferation following liver stress and death (Summers et al. 2003). We accounted for liver proliferation by choosing different numbers of susceptible cells in each patch. Additional data is needed to allow for increase in model complexity through the addition of a proliferation term. Second, we ignored cure of infected cells by non-cytolytic processes (Guidotti et al. 1999; Wieland et al. 2004), and intracellular viral effects (Ciupe and Conway 2024) and assumed that the most important host feature is the immune-mediated killing of infected cells, $δ$ . It is known that HBV-induced liver disease is immune-mediated (Ciupe, Ribeiro, Nelson, Dusheiko, and Perelson 2007; Ciupe et al. 2014; Ciupe and Conway 2024), so more detailed models that account for immune populations and immune functions will be part of future work. Third, we used an oversimplified model for the abnormal nodular liver structure following disease, that is composed of just two patches. While this allowed for analytical tractability of our results, a multi-patch model may result in richer dynamics and even patterns of infection. Forth, we assumed either single directional movement or equal movement between patches. A bifurcation diagram of model Eq. 2 (Appendix B, Fig. 9) shows that, for $0 < Φ_{12} < Φ$ and $0 < Φ_{21} < Φ$ , we observe symmetry between virus levels in the two patches. At the $Φ_{12} = 0$ or $Φ_{21} = 0$ extreme, we recover the one-directional movement model Eq. 3; and at the $Φ_{12} = Φ_{21} = Φ$ extreme, we recover the two-directional movement model Eq. 9. Further work is needed to uncover dynamics in the intermediate movement ranges. Lastly, to identify parameters, we had to assume that the initial conditions are known and $I_{1} (0) \neq 0$ . While most within-host models assume that $I_{1} (0) = 0$ , we bypassed the initial infection and started with one infected cell $I_{1} (0) = 1$ , making parameter estimates dependent (in a non-significant way - results not shown) on this assumption.

**Fig. 9** Two-dimensional bifurcation diagrams for model Eq. 2: ${log}_{10} (V_{1} + V_{2})$ at equilibrium versus $Φ_{12}$ and $Φ_{21}$ . The other parameters are given in Table 2, case 2

In conclusion, we developed mathematical models that considered the formation of liver nodular structures following hepatitis B viral infection and disease and used them to investigate virus dynamics within-patches and systemically. We predicted different outcomes when viral movement between patches is either irreversible or reversible. Moreover, we found that cell susceptibility to infection within nodular structures, the movement rate between patches, and the immune-mediated infected cell killing have the most influence on the results.

Materials and Methods

Structural Identifiability Analysis

Global structural identifiability analysis determines whether a model’s parameters can be uniquely revealed given unlimited noise-free data (for reviews, see (Tuncer and Le 2018; Liyanage et al. 2024; Tuncer et al. 2025; Liyanage et al. 2025)).

Consider a general compartment model:

\begin{matrix} \begin{matrix} \dot{x} (t) & = f (x (t), q), \\ x (0) & = x_{0}, \\ y (t) & = g (x (t), p) . \end{matrix} \end{matrix}

12

Here,

\begin{matrix} x (t) = {T_{1} (t), I_{1} (t), V_{1} (t), T_{2} (t), I_{2} (t), V_{2} (t)} \in R^{6}, \end{matrix}

are the state variables at time t, that solve the ordinary differential equations governed by the rules:

\begin{matrix} f (x (t), q) \in R^{6}, \end{matrix}

for models Eqs. 3 and 9. The parameter vector is:

\begin{matrix} q = {s_{1}, s_{2}, β, d, δ, c, p, ϕ} \in R^{8}, \end{matrix}

and the initial condition vector is:

\begin{matrix} x_{0} \in R^{6} . \end{matrix}

We assumed that the empirical observation represented by

y (t) \in R^{4}

can be explained by models Eqs. 3 and 9 through function:

\begin{matrix} g (x (t), p) \in R^{6} . \end{matrix}

Definition 1

System Eq. 12 is said to be globally structurally identifiable for parameter vector

q

if, for every other parameter vector

\hat{q},

\begin{matrix} y (t, q) = y (t, \hat{q}) implies q = \hat{q} ; \end{matrix}

it is locally structurally identifiable for parameter vector

q

if, for every other parameter vector

\hat{q}

,

\begin{matrix} y (t, q) = y (t, \hat{q}), and \hat{q} \in B (q) implies q = \hat{q}, where B (q) is a ball centered at q ; \end{matrix}

and it is unidentifiable when at least one of its parameters fails the local identifiability test.

The methodology for determining the structural identifiability of systems of ordinary differential equations can range from the Taylor series approach (Pohjanpalo 1978), the differential algebra approach (Bellu et al. 2007; Ljung and Glad 1994), the generating series approach (Walter and Pronzato 1997), the implicit functions approach (Xia and Moog 2003), and several others. Additionally, there are several tools (platforms) that can assist in establishing whether a system of ordinary differential equations is identifiable, such as the COMBOS (Meshkat et al. 2014), the Differential Algebra for Identifiability of SYstems (DAISY, Bellu et al. (2007)), the Exact Arithmetic Rank (EAR, Anguelova et al. (2012)), the GenSSI2 (Ligon et al. 2018), the Structural Identifiability Analyser (SIAN, Hong et al. (2019)), the STRIKE-GOLDD (Villaverde et al. 2016), and the StructuralIdentifiability.jl (Dong et al. 2023). In this study, we used the differential algebra approach and the DAISY platform (Bellu et al. 2007) to determine the structural identifiability of systems Eqs. 3 and 9. The goal of the structural identifiability is to determine which parameters we can confidently estimated from unlimited observations.

Data Fitting

We used previously published longitudinal serum HBV DNA titers from one chimeric mouse with humanized liver (Zhang et al. 2023). Briefly, a urokinase-type plasminogen activator and severe combined immunodeficient (uPA/SCID) mouse was transplanted human hepatocytes (Tateno et al. 2015). Approximately 120–150 days after transplantation (when the growth and proliferation of liver cells is complete), the mouse was infected with hepatitis B virus. Serum HBV DNA was collected during the expansion, peak and persistent stage of HBV infection at days $t_{data}$ ={14, 22, 33, 54, 82, 99, 120, 141, 162, 183, 197, 212} post inoculation. HBV DNA dynamics reached steady levels, similar to an acute HBV infection that becomes persistent in humans (Zhang et al. 2023).

Data fitting procedure for the two-patch models Eqs. 3 and 9

Known parameters. We assumed that hepatocytes have a life-span of 100 days (M.-J. Wang et al. 2017), resulting in an uninfected hepatocyte death rate of $d = 0.01$ /day. Moreover, we assumed that $6.8 \times 10^{5}$ hepatocytes/ml are susceptible to HBV infection (Ciupe et al. 2024), resulting in $T (0) = T_{1} (0) + T_{2} (0) = 6.8 \times 10^{5}$ cells/ml and a total recruitment rate $s = d \times T (0) = 6.8 \times 10^{3}$ cells/(ml $\times$ day). HBV is cleared at rate $c = 4.4$ /day (Murray et al. 2005). Since we fit the model to data from immunosupressed mice, we assumed $δ = d = 0.01$ /day.

We assumed that the HBV inoculum is seeded in patch 1, and set patch 1-specific initial virus and infected cells to $V_{1} (0) = 10^{4}$ HBV DNA/ml and $I_{1} (0) = I_{0} = 1$ cells/ml (Zhang et al. 2023). Moreover, we assumed that there is no virus in patch 2 at the beginning of infection. Therefore, $I_{2} (0) = 0$ cells/ml and $V_{2} (0) = 0$ HBV DNA/ml. We split the number of uninfected hepatocytes between the two patches. $T_{1} (0) = T_{2} (0) = 3.4 \times 10^{5}$ hepatocytes/ml, and assumed three cases for hepatocyte recruitment: (case 1) $s_{1} = 0.1 \times s = 6.8 \times 10^{2}$ cells/(ml $\times$ day), $s_{2} = 0.9 \times s = 6.12 \times 10^{3}$ cells/(ml $\times$ day); (case 2) $s_{1} = s_{2} = 0.5 \times s = 3.4 \times 10^{3}$ cells/(ml $\times$ day); and (case 3) $s_{1} = 0.9 \times s = 6.12 \times 10^{3}$ cells/(ml $\times$ day), $s_{2} = 0.1 \times s = 6.8 \times 10^{2}$ cells/(ml $\times$ day).

Data fitting algorithm. The remaining parameters

π = {β, p, ϕ}

are assumed unknown and are estimated by minimizing the functional:

\begin{matrix} J (π) = {(\sum_{t_{data}}, {({log}_{10} (V_{1} (t_{data}, π) + V_{2} (t_{data}, π)) - {log}_{10} V_{data} (t_{data}))}^{2})}^{1 / 2}, \end{matrix}

over the parameter space

π

using the built-in function fminsearchbnd in MATLAB R2021a. Parameter bounds are

10^{- 10} \leq β \leq 10^{- 7}

ml/(virion

\times

day),

0.1 \leq ϕ \leq 5

/day, and

0 \leq p \leq 1500

virus/(ml

\times

day) for (case 1 - case 3). The initial guesses are

β = 5 \times 10^{- 9}

ml/(virion

\times

day),

ϕ = 0.5

/day and

p = 100

virus/(ml

\times

day). The resulting values for model Eq. 3 are given in Table 2 and the dynamics of model Eq. 3 over time are given in Fig. 2. Similarly, the resulting values for model Eq. 9 are given in Table 4 and the dynamics of model Eq. 9 over time are given in Fig. 5.

Table 2 Parameter values for model Eq. 3

Fixed Parameters all cases	Description	Value	Reference
c	Virus clearance rate	4.4 /day	Murray et al. (2005)
d	Uninfected hepatocyte death rate	0.01 /day	M.-J. Wang et al. (2017)
$δ$	Killing rate	0.01 /day
Fixed Parameters case 1	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$6.8 \times 10^{2}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	6.12 $\times 10^{3}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 1	Description	Value	RSS
$β$	Infectivity rate	$3.3 \times 10^{- 9}$ ml/(virion $\times$ day)	0.85
p	Virus production rate	998 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	0.1 /day	-
Fixed Parameters case 2	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$3.4 \times 10^{3}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	$3.4 \times 10^{3}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 2	Description	Value	RSS
$β$	Infectivity rate	$2.63 \times 10^{- 9}$ ml/(virion $\times$ day)	0.69
p	Virus production rate	1203 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	4.1 /day	-
Fixed Parameters case 3	Description	Value	Reference
$s_{1}$	Hepatocyte production patch 1	$6.12 \times 10^{3}$ cells/(ml $\times$ day)	calculated
$s_{2}$	Hepatocyte production patch 2	$6.8 \times 10^{2}$ cells/(ml $\times$ day)	calculated
Estimated Parameters case 3	Description	Value	RSS
$β$	Infectivity rate	$3.13 \times 10^{- 9}$ ml/(virion $\times$ day)	0.52
p	Virus production rate	1137 virion/(ml $\times$ day)	-
$ϕ$	Movement rate	5 /day	-
Initial conditions	Description	Value	Reference
$T_{1} (0)$	Uninfected hepatocytes patch 1	$3.4 \times 10^{5}$ cells/ml	Ciupe et al. (2024)
$T_{2} (0)$	Uninfected hepatocytes patch 2	$3.4 \times 10^{5}$ cells/ml	Ciupe et al. (2024)
$I_{1} (0)$	Infected hepatocytes patch 1	1 cell/ml	-
$I_{2} (0)$	Infected hepatocytes patch 2	0 cell/ml	-
$V_{1} (0)$	Inoculum patch 1	$10^{4}$ HBV DNA/ml	Zhang et al. (2023)
$V_{2} (0)$	Inoculum patch 2	0 HBV DNA/ml	-

Practical Identifiability Analysis

Practical identifiability is a methodology that considers the noisiness in the data set for a given model, and performs data fitting under noise considerations.

Definition 2

A model is practically identifiable if a unique parameter set can be consistently obtained through the fitting procedure to noisy data.

There are several approaches and conditions for assessing practical identifiability (Wu et al. 2008; Eisenberg et al. 2013; Roosa and Chowell 2019; Saucedo et al. 2024; Tuncer and Le 2018). In this paper, we used the Monte Carlo (MC) approach to analyze the practical identifiability of the one-directional two-patch model Eq. 3 and two-directional two-patch model Eq. 9. The Monte Carlo approach is a sampling technique that uses random numbers and probability distributions to determine the practical identifiability of a model (Metropolis and Ulam 1949). We performed MC simulations by generating

M = 1, 000

data sets using the true parameter set

π

and adding noise to the data in increasing amounts. The MC simulations are outlined in the following steps:

Solve model Eq. 3 (Eq. 9, respectively) numerically with the true parameter vector $π$ to obtain the output vector $g (x (t), π)$ at discrete time points $t_{data} = {t_{i}}_{i = 1}^{n}$ .
Generate $M = 1, 000$ data sets with a given measurement error. We assume the error follows a normal distribution with mean 0 and variance $σ^{2} (t)$ ; that is, the data are described by:
$\begin{matrix} y_{i, j} = g (x (t_{i}), π) (1 + ϵ_{i, j}), \end{matrix}$
where $ϵ_{ij} \sim N (0, σ)$ at the discrete data time points $t_{data} = {t_{i}}_{i = 1}^{n}$ for all $j = {1, 2, . . ., M}$ data sets.
Estimate the parameter set $π_{j}$ , by fitting model Eq. 3 (model Eq. 9, respectively) to each of the M simulated data sets. This is achieved by minimizing the difference between model Eq. 3 (model Eq. 9, respectively) output and the data generated for the specific scenario:
$\begin{matrix} J (π_{j}) = (\sum_{t_{data}} ({log}_{10} (V_{1} (t_{data}, π_{j}) + V_{2} (t_{data}, π_{j})) - {log}_{10} V_{data} (t_{data}))^{2})^{1 / 2} . \end{matrix}$
This optimization problem is solved in MATLAB R2021a using the built-in function fminsearchbnd, which is part of the Optimization Toolbox. Since fminsearchbnd is a local solver, the optimized minimum value can be influenced by the starting point. To avoid issues related to the starting value, we use the true parameter values $π$ as the initial parameter starting point provided to fminsearchbnd.
Calculate the average relative estimation error (ARE) for each parameter in the set $π$ , as follows:
$\begin{matrix} ARE (π^{(k)}) = 100 % \times \frac{1}{M} \sum_{j = 1}^{M} \frac{|π^{(k)} - π_{j}^{(k)}|}{|π^{(k)}|}, \end{matrix}$ 13
where $π^{(k)}$ is the k-th parameter of the true parameter set $π$ , and $π_{j}^{(k)}$ is the k-th element of $π_{j}$ .
Repeat steps 1 through 4, by increasing the measurement error $σ = {0, 1, 5, 10, 20, 30} %$ .

The objective of this algorithm is to determine if the parameters are sensitive to gradual noise introduced into the given dataset. We used the AREs as a metric to determine the practical identifiability of model Eq. 3 (model Eq. 9, respectively) by applying the definition below (see (Heitzman-Breen et al. 2024) for details).

Definition 3

Let

σ

be the measurement error introduced to a dataset, and let the ARE be the average measurement error in the parameter

π^{(k)}

.

If $0 \leq$ ARE( $π^{(k)}) \leq σ$ , we say $π^{(k)}$ is strongly practically identifiable.
If $σ <$ ARE( $π^{(k)}) \leq 10 \times σ$ , we say $π^{(k)}$ is weakly practically identifiable.
If $10 \times σ <$ ARE( $π^{(k)})$ , we say $π^{(k)}$ is not practically identifiable.

We state that a model is practically identifiable if

π^{(k)}

is practically identifiable for all values of k.

The results for model Eq. 3 are given in Table 5.

The results for model Eq. 9 are given in Table 6.

Appendix A Stability of Equilibrium $E_{4}$

To finalize the proof of Proposition 4.,we need to show that Fig 8 $a_{1} a_{2} a_{2} > a_{3}^{2} + a_{1}^{2} a_{4}$ . We are including here the maple file showing that result.

Appendix B Bifurcation Diagram for the Generalized Model Eq. 2

We are presenting a bifurcation diagram for the general model Eq. 2 where $0 < Φ_{12} < Φ$ and $0 < Φ_{21} < Φ$ . We observe symmetry between the two movements. At either extremes $Φ_{12} = 0$ or $Φ_{21} = 0$ we recover the one-directional movement model Eq. 3. When $Φ_{12} = Φ_{21} = Φ$ , we recover the two-directional movement model Eq. 9.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

K. Castellano, O. Saucedo, S. Ciupe contributed equally to this work.

Change history

10/17/2026

The original article was revised to correct errors in the author names

Change history

10/17/2026

A Correction to this paper has been published: https://doi.org/10.1007/s11538-026-01780-7

Author Contributions

Conceptualization: SMC; Methodology: KC, OS, SMC; Formal analysis and investigation: KC, OS, SMC; Writing- original draft preparation: SMC; Writing- review and editing: KC, OS, SMC; Funding acquisition: OS, SMC; Resources: OS, SMC; Supervision: OS, SMC.

Funding

SMC acknowledges partial support from the NIH NIGMS grant 1R01GM152743-01 and National Science Foundation grant No. 2051820. This research was enabled in part through the Virginia Tech Center for the Mathematics of Biosystems (VTCMB-033500).

Data Availability

All data necessary to replicate the results in this article are available at https://github.com/StancaCiupe/HBV-multi-patch.

Declarations

Conflicts of interest/Competing interests

We declare no conflict of interest associated with this publication.

Code availability

All code necessary to replicate the results in this article are available at https://github.com/StancaCiupe/HBV-multi-patch.

References

AnguelovaMKarlssonJJirstrandMMinimal output sets for identifiabilityMath Biosci2012239113915310.1016/j.mbs.2012.04.00522609467
BatallerRBrennerDAet al.Liver fibrosisJ Clin Investig2005115220921810.1172/JCI24282PMC54643515690074
BelluGSaccomaniMPAudolySD’AngiòLDAISY: a new software tool to test global identifiability of biological and physiological systemsComput Methods Programs Biomed2007881526110.1016/j.cmpb.2007.07.002PMC288853717707944
CastAKumbajiMD’SouzaARodriguezKGuptaAKarnsRTimchenkoNLiver proliferation is an essential driver of fibrosis in mouse models of nonalcoholic fatty liver diseaseHepatology communications2019381036104910.1002/hep4.1381PMC667186331388625
CiupeSMConwayJMIncorporating intracellular processes in virus dynamics modelsMicroorganisms202412590010.3390/microorganisms12050900PMC1112412738792730
CiupeSMDahariHPlossAMathematical models of early hepatitis B virus dynamics in humanized miceBull Math Biol20248655310.1007/s11538-024-01284-2PMC1100393338594319
Ciupe SM, Heffernan JM (2017) In-host modeling Infectious Disease Modelling 2(2):188–202. 10.1016/j.idm.2017.04.002PMC600197129928736
CiupeSMRibeiroRMNelsonPWDusheikoGPerelsonASThe role of cells refractory to productive infection in acute hepatitis B viral dynamicsProc Natl Acad Sci2007104125050505510.1073/pnas.0603626104PMC182926217360406
CiupeSMRibeiroRMNelsonPWPerelsonASModeling the mechanisms of acute hepatitis B virus infectionJ Theor Biol20072471233510.1016/j.jtbi.2007.02.017PMC199481817428501
CiupeSMRibeiroRMPerelsonASAntibody responses during hepatitis B viral infectionPLoS Comput Biol2014107e100373010.1371/journal.pcbi.1003730PMC411742725078553
CiupeSMVaidyaNKFordeJEEarly events in hepatitis B infection: the role of inoculum doseProc R Soc B202128819442020271510.1098/rspb.2020.2715PMC789323633563115
DahariHShudoERibeiroRMPerelsonASModeling complex decay profiles of hepatitis B virus during antiviral therapyHepatology2009491323810.1002/hep.22586PMC371285919065674
DongRGoodbrakeCHarringtonHAPogudinGDifferential elimination for dynamical models via projections with applications to structural identifiabilitySIAM Journal on Applied Algebra and Geometry20237119423510.1137/22M1469067
DusheikoGTreatment of HBeAg positive chronic hepatitis B: interferon or nucleoside analoguesLiver Int20133313715010.1111/liv.1207823286858
EisenbergMCRobertsonSLTienJHIdentifiability and estimation of multiple transmission pathways in cholera and waterborne diseaseJ Theor Biol20133248410210.1016/j.jtbi.2012.12.02123333764
GoyalARibeiroRMPerelsonASThe role of infected cell proliferation in the clearance of acute HBV infection in humansViruses201791135010.3390/v9110350PMC570755729156567
GuidottiLGChisariFVImmunobiology and pathogenesis of viral hepatitisAnnu Rev Pathol: Mech Dis200611236110.1146/annurev.pathol.1.110304.10023018039107
GuidottiLGRochfordRChungJShapiroMPurcellRChisariFVViral clearance without destruction of infected cells during acute HBV infectionScience1999284541582582910.1126/science.284.5415.82510221919
HailegiorgisAIshidaYCollierNImamuraMShiZReinharzVet al.Modeling suggests that virion production cycles within individual cells is key to understanding acute hepatitis B virus infection kineticsPLoS Comput Biol2023198e101130910.1371/journal.pcbi.1011309PMC1042691837535676
Heitzman-BreenNLiyanageYRDuggalNTuncerNCiupeSMThe effect of model structure and data availability on Usutu virus dynamics at three biological scalesRoyal Society Open Science202411223114610.1098/rsos.231146PMC1084694038328567
HewsSEikenberrySNagyJDKuangYRich dynamics of a hepatitis B viral infection model with logistic hepatocyte growthJ Math Biol20106057359010.1007/s00285-009-0278-319533136
HongHOvchinnikovAPogudinGYapCSIAN: software for structural identifiability analysis of ODE modelsBioinformatics201935162873287410.1093/bioinformatics/bty106930601937
KadelkaSDahariHCiupeSMUnderstanding the antiviral effects of RNAi-based therapy in hbeag-positive chronic hepatitis B infectionSci Rep202111120010.1038/s41598-020-80594-6PMC779457033420293
LeeWMHepatitis B virus infectionN Engl J Med1997337241733174510.1056/NEJM1997121133724069392700
LewinSRRibeiroRMWaltersTLauGKBowdenSLocarniniSPerelsonASAnalysis of hepatitis B viral load decline under potent therapy: complex decay profiles observedHepatology20013451012102010.1053/jhep.2001.2850911679973
LigonTSFröhlichFChişOTBangaJRBalsa-CantoEHasenauerJGenSSI 2.0: multi-experiment structural identifiability analysis of SBML modelsBioinformatics20183481421142310.1093/bioinformatics/btx735PMC590561829206901
LiyanageYRChowellGPogudinGTuncerNStructural and practical identifiability of phenomenological growth models for epidemic forecastingViruses202517449610.3390/v17040496PMC1203129740284940
LiyanageYRHeitzman-BreenNTuncerNCiupeSMIdentifiability investigation of within-host models of acute virus infectionMathematical biosciences and engineering MBE202421107394742010.3934/mbe.2024325PMC1218223739696868
LjungLGladTOn global identifiability for arbitrary model parametrizationsAutomatica199430226527610.1016/0005-1098(94)90029-9
LocarniniSHatzakisAChenD-SLokAStrategies to control hepatitis B: public policy, epidemiology, vaccine and drugsJ Hepatol2015621S76S8610.1016/j.jhep.2015.01.01825920093
MeenaMPurohitMPurohitSDNisarKSet al.A novel investigation of the hepatitis B virus using a fractional operator with a non-local kernelPartial Diff Equ Appl Math2023810057710.1016/j.padiff.2023.100577
MeshkatNKuoCEDiStefanoJIIIOn finding and using identifiable parameter combinations in nonlinear dynamic systems biology models and COMBOS: a novel web implementationPLoS ONE2014910e11026110.1371/journal.pone.0110261PMC421165425350289
MetropolisNUlamSThe Monte Carlo methodJ Am Stat Assoc19494424733534110.2307/228023218139350
MichalopoulosGKBhushanBLiver regeneration: biological and pathological mechanisms and implicationsNat Rev Gastroenterol Hepatol2021181405510.1038/s41575-020-0342-432764740
MurrayJMWielandSFPurcellRHChisariFVDynamics of hepatitis B virus clearance in chimpanzeesProc Natl Acad Sci200510249177801778510.1073/pnas.0508913102PMC134572416306261
NeumannAUPhillipsSLevineIIjazSDahariHErenRNaoumovNVNovel mechanism of antibodies to hepatitis B virus in blocking viral particle release from cellsHepatology201052387588510.1002/hep.23778PMC308635720593455
NowakMABonhoefferSHillAMBoehmeRThomasHCMcDadeHViral dynamics in hepatitis B virus infectionProc Natl Acad Sci19969394398440210.1073/pnas.93.9.4398PMC395498633078
PawlotskyJ-MNew hepatitis B drug development disillusions: time to reset?Lancet Gastroenterol Hepatol20238219219710.1016/S2468-1253(22)00341-736343654
PohjanpaloHSystem identifiability based on the power series expansion of the solutionMath Biosci1978411–2213310.1016/0025-5564(78)90063-9
ReinharzVIshidaYTsugeMDurso-CainKChungTLTatenoCDahariHUnderstanding hepatitis B virus dynamics and the antiviral effect of interferon alpha treatment in humanized chimeric miceJ Virol2021951410112810.1128/jvi.00492-20PMC822395633910953
RoehlenNCrouchetEBaumertTFLiver fibrosis: mechanistic concepts and therapeutic perspectivesCells20209487510.3390/cells9040875PMC722675132260126
RoosaKChowellGAssessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission modelsTheor Biol Med Model201916111510.1186/s12976-018-0097-6PMC633283930642334
SaucedoOLaubmeierATangTLevyBAsikLPollingtonTProsperOComparative analysis of practical identifiability methods for an seir modelAIMS Mathematics202499247222476110.3934/math.20241204
SchweitzerAHornJMikolajczykRTKrauseGOttJJEstimations of worldwide prevalence of chronic hepatitis B virus infection: a systematic review of data published between 1965 and 2013The Lancet2015386100031546155510.1016/S0140-6736(15)61412-X26231459
Shekhtman L, Ishida Y, Tsuge M, Reinharz V, Yamao M, Takahashi M, Chayama K (2024) Modelling of hepatitis B virus kinetics and accumulation of cccDNA in primary human hepatocytes. JHEP Reports 101311. 10.1016/j.jhepr.2024.101311PMC1199760740236628
SummersJJilbertARYangWAldrichCESaputelliJLitwinSMasonWSHepatocyte turnover during resolution of a transient hepadnaviral infectionProc Natl Acad Sci200310020116521165910.1073/pnas.1635109100PMC20881314500915
TatenoCKawaseYTobitaYHamamuraSOhshitaHYokomichiHet al.Generation of novel chimeric mice with humanized livers by using hemizygous cDNA-uPA/SCID micePLoS ONE20151011e014214510.1371/journal.pone.0142145PMC463311926536627
TrépoCChanHLLokAHepatitis B virus infectionThe Lancet201438499592053206310.1016/S0140-6736(14)60220-824954675
TsiangMRooneyJFTooleJJGibbsCSBiphasic clearance kinetics of hepatitis B virus from patients during adefovir dipivoxil therapyHepatology19992961863186910.1002/hep.51029062610347131
TuncerNLeTTStructural and practical identifiability analysis of outbreak modelsMath Biosci201829911810.1016/j.mbs.2018.02.00429477671
Tuncer N, Martcheva M, Ciupe SM (2025) Structural and practical identifiability of within-host models of virus dynamics-a review. Curr Opin Syst Biol 100552. 10.1016/j.coisb.2025.100552PMC1264820941312440
VillaverdeAFBarreiroAPapachristodoulouAStructural Identifiability of Dynamic Systems Biology ModelsPLoS Comput Biol20161210e100515310.1371/journal.pcbi.1005153PMC508525027792726
WalterEPronzatoLIdentification of parametric models: from experimental data1997Springer Verlag
WangKFanATorresAGlobal properties of an improved hepatitis B virus modelNonlinear Anal Real World Appl20101143131313810.1016/j.nonrwa.2009.11.008
WangM-JChenFLauJTHuY-PHepatocyte polyploidization and its association with pathophysiological processesCell death & disease201785e2805e280510.1038/cddis.2017.167PMC552069728518148
WielandSFChisariFVStealth and cunning: hepatitis B and hepatitis C virusesJ Virol200579159369938010.1128/jvi.79.15.9369-9380.2005PMC118154816014900
WielandSFSpangenbergHCThimmeRPurcellRHChisariFVExpansion and contraction of the hepatitis B virus transcriptional template in infected chimpanzeesProc Natl Acad Sci200410172129213410.1073/pnas.0308478100PMC35706314764900
WuHZhuHMiaoHPerelsonASParameter identifiability and estimation of hiv/aids dynamic modelsBull Math Biol200870378579910.1007/s11538-007-9279-918247096
XiaXMoogCHIdentifiability of Nonlinear Systems with Application to HIV/AIDS ModelsIEEE Trans Autom Control200348233033610.1109/TAC.2002.808494
YousfiNHattafKTridaneAModeling the adaptive immune response in HBV infectionJ Math Biol20116393395710.1007/s00285-010-0397-x21234570
YuanLJiangJLiuXZhangYZhangLXinJet al.HBV infection-induced liver cirrhosis development in dual-humanised mice with human bone mesenchymal stem cell transplantationGut201968112044205610.1136/gutjnl-2018-316091PMC683973530700543
Zhang B-H, Zhou Y, Tempel B, Pan H, Horrigan S, Luckenbaugh L, Zhang Y-Y (2023) Replication-driven HBV cccDNA loss in chimeric mice with humanized livers. bioRxiv, 2023–12, 10.1101/2023.12.28.573542PMC1264592141196060

Copyright and License

License: https://creativecommons.org/licenses/by/4.0/

Download

PDF JATS XML

Abstract

Introduction

Mathematical Model Development

Results for the One-Directional Two-Patch Model

Asymptotic Analysis Results for the One-Directional Two-Patch Model Eq. 3

Proposition 1

Proof

Proposition 2

Proof

Proposition 3

Proof

Structural Identifiability Results for Model Eq. 3

Proposition 4

Proof

Proposition 5

Proof

Numerical Results for Model Eq. 3

Results for the Two-Directional Two-Patch Model

Asymptotic Analysis Results for the Two-Directional Two-Patch Model Eq. 9

Proposition 6

Proof

Structural Identifiability Results for Model Eq. 9

Proposition 7

Proof

Proposition 8

Proof

Numerical Results for Model Eq. 9

Discussion and Conclusion

Materials and Methods

Structural Identifiability Analysis

Definition 1

Data Fitting

Data fitting procedure for the two-patch models Eqs. 3 and 9

Practical Identifiability Analysis

Definition 2

Definition 3

Appendix A Stability of Equilibrium E4

Appendix B Bifurcation Diagram for the Generalized Model Eq. 2

Author Contributions

Funding

Data Availability

Declarations

Conflicts of interest/Competing interests

Code availability

References

Copyright and License

Appendix A Stability of Equilibrium $E_{4}$