Modeling the interplay between albumin-globulin metabolism and HIV infection

1.   Introduction

Human immunodeficiency virus (HIV) infection is a serious public health concern in the USA and globally. Long term HIV infection may progress to the acquired immunodeficiency syndrome (AIDS). For individuals with HIV, maintaining proper nutrition is crucial in sustaining overall health and controlling the severity of the symptoms that are caused by HIV and AIDS ^[1]. A consistent healthy eating style sustains the immune system and helps support the strength of the immune responses when individuals are infected with HIV. Despite the fact that HIV and albumin/globulin metabolism are strongly interlinked, no mathematical models seem to exist that model the feedbacks between HIV and albumin/globulin. HIV primarily affects the immune system, and globulin is critical for the immune system, which uses various elements to build immune system cells ^[2].

Thus, HIV infection and the medications taken to control it, increase nutritional need in HIV-infected individuals. However, not all HIV-infected individuals have proper access or interest in maintaining adequate eating habits, which can lead to malnutrition ^[3]. In contrast to other infections in which carbohydrates are first broken down to meet the increased needs of the immune system, proteins are first broken down in HIV ^[3], resulting in HIV-specific wasting and weight loss. Malnutrition, in turn, lowers immunity and increases viral load ^[3], as well as increases susceptibility to opportunistic infections. Previous research suggests that improvements in nutrition acted as a form of early intervention for HIV patients. However, these studies lack mathematical models that incorporate this relation ^[4,5].

HIV primarily targets immune cells. The optimal function of the immune system relies on proper nutrition and inflow of globulins, and the serum total protein assessment determines the patient's nutritional status. Albumin and globulin are two types of proteins with levels that are indicators of HIV disease progression ^[6]. Albumin is a key protein in the blood and is responsible for transporting a wide variety of substances, such as hormones, drugs, vitamins, and minerals. Immunoglobulins and other proteins involved in the immune response are part of a large family of proteins called globulins. HIV infection often correlates with elevated serum protein levels; however, while globulin concentrations increase, albumin concentrations decrease ^[5,7,8]. Furthermore, decreased albumin levels have been linked to the progression of HIV to AIDS ^[8] and in cases with low albumin levels, viral infections such tend to be more severe ^[9], although the specific mechanisms underlying this phenomenon are not yet fully understood. The complicated relationship between albumin, globulin, immune system, and HIV has long been recognized ^[3,4,5,7,8]. However, the use of mathematical models in this context has been limited. The use of mathematical models to demonstrate the interactions among albumin, globulin, the immune system, and HIV is critical because these models can reveal the deep links between these systems and improve our knowledge of their interactions.

Nutrition in the absence of HIV has been modeled both through mathematical models and statistical methods ^{[10,11,12,13]}. Mathematical models of nutrition have been the topic of a conference in which selected papers were published in ^[14]. A within-host model of nutrition and challenges with its calibration to data has been discussed by Juillet et al. ^[15]. A model based on chemical kinetics for the rate of utilization and storage of carbohydrates, fats, and proteins was introduced and discussed in ^[16]. A model of various feeding determinants that would facilitate the simulation of food intake, thus providing information for existing models of nutrient absorption and metabolism, was developed in ^[17]. Mathematical models of nutrition at the population level have also been developed. Shah et al. ^[18] apply the optimal control theory to a model of poor nutrition at the population level and compute the reproduction number. Baleanu et al. ^[19] investigate a fractional order differential equation model of poor nutrition at the population level. Nutrition in ruminants has been studied in two books that focused on nutritional system of ruminants through modeling ^[20]. Within-host modeling of nutrients have been considered in an edited book by Coburn and Townsend ^[21]. These approaches aim to understand the complex relationships between various nutritional factors and their impacts on health and well-being. In the context of our current study, we delve deeper into the complexities of nutrition by focusing on two key protein markers, albumin and globulin, within the context of HIV infection. These proteins serve vital roles in various bodily functions, including immune response. Given the limited understanding of how albumin/globulin interacts with viral infections, we begin our study by constructing a family of models. These models serve as the foundation for our study, allowing us to systematically examine and comprehend the complex dynamics at play.

Our models here are built on the classical HIV models that incorporate target cells, infected target cells and virus ^[22,23]. We further incorporate the immune response in the form of CD8 cells. The immune response has also previously been considered in HIV models ^[22,24]. The novelty of our models lies in the incorporation of nutritional elements, and more specifically proteins composed of albumin and globulins ^[4,25]. We incorporate albumin and globulin separately because albumin and globulins have different responses to HIV infection. In particular, albumin levels decrease in HIV-infected individuals whereas globulins increase ^[5]. This basic interdependence between HIV and the albumins and globulins is known, however, it is less known how albumins and globulins specifically interact with the immune system. Thus, our main question in this article is "What kind of model will capture best the data of target cells, virus, albumin and globulin?" Because we model these interactions for the first time, we would like to develop sound, well validated model that represents the available data ^[26]. Relatively little is known how proteins impact the infection of new target cells, how they impact the clearance of the immune response of the infected target cells or how the globulins affect the recruitment of CD4 cells and CD8 cells. To determine this dependence, we hypothesize a general form of these terms. We obtain a general model that has 17 specific submodels. We perform model selection on these 17 models and use AIC ^[27] to select the model that best describes the data.

We analyze the best model. We further fit the best model to the data and consider the well-posedness of the parameter estimation problem. In other words, we study the structural identifiability of the best fitted model. To evaluate its fit to the data and the quality of the parameter estimation from the real data, we perform Monte Carlo Simulations (MCS) and compute the Average Relative Estimation Erors (AREs) of the estimated parameters.

2.   Methods

2.1.   Mathematical models of HIV and Albumin/Globulin metobolism

Nutrition, in general, is made up of macronutrients and micronutrients, both of which are required for proper functioning of the immune system. All macronutrients are important, but in HIV the proteins are the key. The total protein in the body, often measured in grams per deciliter (g/dL) ^[28], consists of albumin and globulin. In healthy individuals the ratio of albumin to globulin is around 1.7 ^[25]. Total protein levels in HIV-infected people are higher, as are globulin levels, but albumin levels are lower ^[5,29]. Globulin is required for the production of elements of the immune system, most notably antibodies, whereas albumin has been found to be a good predictor of HIV progression, morbidity, and mortality ^[6,30]. We begin modeling the interplay between HIV and albumin/globulin by incorporating protein to the well studied target-cell limited model ^[23]. While modeling protein coupled with a disease like HIV is new, our within-host albumin/globulin-HIV model is based on an extensive literature of within HIV models ^[23] and within-host HIV models with immune response ^[24,31]. Our advancement here consists of adding albumin and globulin to the well studied HIV within-host model ^[23].

We used six differential equations to describe the temporal dynamics of six populations: target cells ($T$), virus-producing infected target cells ($T_i$), viral load ($V)$, immune cells ($Z$), albumin ($A$) and globulin ($G$). Figure 1 presents a flow diagram of the model. The within-host model of HIV and protein is then built with the following system of ordinary differential equations:

HIV primarily targets CD4+ cells in the immune system. Since protein, specifically globulin, impacts the proliferation of immunoglobulins (immune cells), we represent the production of target cells and immune cells as a function of globulin. Therefore, in the within-host model (2.1), we set the target and immune cell production depending on the individual's globulin level. Namely,

where, $r_0$ is the production of target cells per unit globulin, and $r$ is the production rate of target cells. Similarly, $\lambda_{z_0}$ is the recruitment rate of the immune cells per unit globulin, and $\lambda_z$ is the recruitment rate of the immune cells. Virus particles, $V$, infect target cells, $T$, to produce infected cells, as expressed with the term $\rho(A, G)TV$. By setting the infection rate, $\rho(A, G)$, as a function of albumin and globulin, we account for the influence of albumin/globulin on the target cell susceptibility ^[32]. Target cells die at rate $d$. Infected target die at rate $\delta$ and their clearance by the immune system is denoted by the term $\psi(A, G)T_i Z$. The ability of immune cells to attack infected target cells, $\psi(A, G)$, is considered to be a function of albumin and globulin ^[32]. Infected cells generate additional virus particles at rate $\pi$, while the virus particles are cleared from the system at rate $c$. Furthermore, immunoglobulins remove the viral particles at rate $\mu_v$. Immune cells are activated at rate $b$ in response to antigen and cleared at rate $\mu_z$.

Albumin and globulin intake are denoted by $\lambda_A$ and $\lambda_G$ respectively. The parameter $\mu_A$ represents the albumin consumption and clearance, whereas $\mu_G$ represents the globulin consumption and clearance. $\gamma_A$ is the virus-driven depletion of albumin, while $\gamma_G$ is the virus-driven increase in globulin. The functional forms of the infection rate of target cells $\rho(A, G)$ and the attack rate of immune cells $\psi(A, G)$ are defined as:

where $\rho_0$, $A_i$ and $\psi_{0i}$, $i = 1, 2$, are proportionality constants, and $\Psi$ is a half-saturation constant. The constants $\Psi, A_1, A_2, \psi_{01}$ and $\psi_{02}$ allow us to remove the dependence of the transmission rate $\rho$ and killing rate $\psi$ on the albumin/globulin. We performed model selection on the global model (2.1) to determine more robust forms for the terms $h_1(G)$, $h_2(G)$, $\rho(A, G)$ and $\psi(A, G)$.

2.2.   Model selection process

We proposed a family of models that will aid in clarifying the complex relationship between HIV and protein dynamics. However, it is not clear which model would be the best for this purpose. It is standard to judge the suitability of a model by comparing its ability to match the experimental data, but this approach is known to favor model complexity. On the other hand, Akaike information criteria (AIC) balances model complexity and the goodness-of-fit. We use AIC to determine the robust forms of the terms depending on albumin and globulin; $\lambda(G)$, $\lambda_z(G)$, $\rho(A, G)$ and $\psi(A, G)$. The set of candidate models is obtained by setting $r_0 = 0$ and/or $\lambda_{z0} = 0$ and/or $A_1 = 0$ and/or $A_2 = 0$ and/or $\psi_{01} = 0$ and/or $\psi_{02} = 0$. We also explored the scenario when $\psi_{01} = \psi_{02} = \psi_0$ and $\Psi = 0$. In total, we compared $17$ models (see Table 1).

Within-Host HIV and Albumin/Globulin Model:

To validate the models, we used data from a rhesus monkey study of simian immunodeficiency virus (SIV), in which CD4+ T-cell and plasma viral load were measured for 12 weeks following SIV infection ^[33]. We took albumin levels from another SIV experiment in which albumin concentrations were measured for 16 weeks after SIV infection in 8 pig-tailed macaques ^[34]. For globulin levels, we obtained data from the published study ^[35] in which globulin levels were assessed in the peripheral blood of rhesus monkeys infected with SIV (see Table 2).

The observations are the state variables of the model, namely CD4+ cells is the state variable $T(t, \boldsymbol p)$, virus is $V(t, \boldsymbol p)$, albumin is $A(t, \boldsymbol p)$, and globulin is $G(t, \boldsymbol p)$ in (2.1), with $p$ denoting the parameters of the model. Observations vary as the model parameters $p$ vary. Clearly, data collected in the experimental setting do not fall on the smooth path of the observations $T(t, \boldsymbol p)$, $V(t, \boldsymbol p)$, $A(t, \boldsymbol p)$ and $G(t, \boldsymbol p)$. We represent the relationship between the experimental data and model observations with the following statistical model,

where $E_i, E_j, E_k, E_l$ are the measurement errors with Gaussian structure. The measurement errors are independent and normally distributed with mean zero, and standard deviation $\sigma$. Parameter estimation based on a maximum likelihood principle and least square principle are equivalent when the measurement errors are independent and normally distributed. The least square estimate can be obtained by minimizing the following objective function, which will be refered as sum of squared errors (SSE),

Several challenges come up when fitting multiple data sets to a within-host model. The data sets, to begin with, exhibit a wide range of variations in magnitudes. For instance, viral load ranges between $10^2$ to $10^7$, while albumin levels fluctuates from $1$ to $5$. Moreover, any fitting algorithm that is used to minimize the SSE tends to prioritize the superior fitting of one data set, at the expense of a poorer fitting fit of the other data set. To address these issues, we take the logarithms of the viral load and target cell counts to normalize the data magnitudes, along with assigning specific weights to each data sets. This strategic weighting forces the optimization process to give equal consideration for each data set. SSE values alone may indicate the better model when models have the same number of parameters. Since all 17 models considered have different set of parameters, in this study we employed the AIC to assess which model explains the data best. The corrected AIC for the finite sample size is given by ^[27],

where $n$ is the number of data points and $K$ is the number of parameters fit by the model plus one ^[27]. AIC balances the model complexity and the SSE by penalizing the models with higher parameters. Due to this penalty, a more complex model needs to improve its fit to compensate for its higher number of parameters. For instance, model 12 has a higher AIC$_c$ despite having the minimum SSE value (see Table 3). Since the AIC score of each model varies greatly, it is standard to rescale them with the minimum AIC score by calculating the difference between the AIC$_c$ score and AIC$_{c_{min}}$. The rescaled AIC$_c$ scores of each model, $\Delta = AIC_c - AIC_{c_{min}}$, are reported in Table 3. The rescaled AIC values are a direct way to compare the models, only the models with $\Delta < 2$ have strong support ^[27]. Based on our analysis, model $1$ is the only model that has strong support for the experimental data. The last column of Table 3 represents the relative likelihood of the model $j$, given by

The models that best describe the data are those whose relative likelihood sum up to $\sum w_j > 0.95$. The model 1 exclusively explains the data since its relative probability is higher than the threshold of $0.95$.

3.   Analysis of the best selected model

The best model takes the form:

3.1.   The infection-free equilibrium and immunological reproduction number

Using the Van den Driessche and Watmough Approach, first we derive the immunological reproduction number of the within-host model of HIV albumin and globulin given by system (3.1). We arrange the equations in (3.1) so that the first $2$ components of the ODE system correspond to the infected compartments, and the last $4$ components correspond to the non-infected compartments. We rewrite the system as,

where $x_m = (T_i\; \; V)^T$ and $y_n = (T\; \; Z\; \; A\; \; G)^T.$ Now, define ${F}_1 = \rho_0 T V,$ ${F}_2 = \pi T_i,$ ${V}_1 = \delta T_i+\psi_{01} T_i Z,$ and ${V}_2 = c V + \mu_V G V$ so that we have

Now, the infection-free system $y' = g(0, y)$ has the unique equilibrium $\epsilon_0 = (0, y_0),$ where $y_0 = (T_0, Z_0, A_0, G_0) = \left(\frac{r}{d}, \; \; \; \; \frac{\lambda_z}{\mu_z}, \; \; \; \; \frac{\lambda_A}{\mu_A}, \; \; \; \; \frac{\lambda_G}{\mu_G}\right).$ Observe that,

Hence, the next-generation matrix is defined by

This gives us the immunological reproduction number

3.2.   Stability of the infection-free equilibrium

The infection-free equilibrium of the system (3.3) is given by $E_0 = (T_0, 0, 0, Z_0, A_0, G_0),$ where $T_0 = \frac{r}{d},$ $Z_0 = \frac{\lambda_Z}{\mu_Z}, A_0 = \frac{\lambda_A}{\mu_A}$ and $G_0 = \frac{\lambda_G}{\mu_G}.$ We state the following Theorem 1 regarding the stability of the infection-free equilibrium and present the proof in the Appendix.

Theorem 3.1. If $\mathcal {R}_{0} < 1$, then the infection-free equilibrium is locally asymptotically stable. If ${R}_{0} > 1$, it is unstable.

3.3.   Existence and stability of infection equilibrium

We set the left hand side of the model (3.1) equals to $0$ and obtain the following:

By the sixth equation of (3.3) we get that

Since we need system variables to stay positive, we set $G^* > 0$ and obtain $0 < V^* < \frac{\mu_G}{\gamma_G}.$ Similarly, the fifth equation of (3.3) leads to

Upon examining the fixed-point formulas (3.4) and (3.5), we notice that as $V^*$ increases, the globulin levels at equilibrium, denoted by $G^*$, increase, while albumin levels at equilibrium, denoted by $A^*$, decrease. This is consistent with clinical observations that individuals with severe disease often exhibit lowered levels of albumin. Furthermore, a correlation exists between elevated viral load values and increased mortality ^[36], which is effectively captured by the provided fixed-point formulas (3.4) and (3.5). The first equation of (3.3) gives

Both $A^*$ and $T^*$ are trivially positive. Now, solving the third equation of (3.3) we get

which is also positive since $0 < V^* < \frac{\mu_G}{\gamma_G}.$ Next, we solve the fourth equation of (3.3) for $Z^*.$ Substituting (3.7), we obtain

In order to obtain a positive solution for $Z^*,$ we must satisfy

The left hand side of the Eq (3.9) is a straight line that goes through the origin while the right hand side is a parabola with roots ${V_1}^* = \frac{\mu_G}{\gamma_G}$ and ${V_2}^* = \frac{\pi \mu_Z}{b c}$. The parabola opens up as the coefficient of $V^*$ is positive and meets the straight line $y = b \lambda_G \mu_V V^*$ at two distinct points, say $V_1$ and $V_2$, with $V_1 < V_2.$ Then, $V_1 < {V_1}^*$ and $V_1 < {V_2}^*$, for any ${V_1}^*, {V_2}^*.$ Thus, ${V_1} < \frac{\mu_G}{\gamma_G}.$ As this satisfies our requirement to have positive values at infection equilibrium, $Z^* > 0.$ Finally, substituting $T^*, {T_i}^*$ and $Z^*$ in the second equation of (3.3), $\delta {T_i}^* + \psi_{01} {T_i}^* Z^* = \rho_0 T^*V^*,$ we get

Let $y_1(V^*) = \frac{\delta}{\pi}\left(c + \frac{\lambda_G \mu_V}{\mu_G - \gamma_G V^*}\right) + \frac{{\psi}_{01}}{\pi}\left(c + \frac{\lambda_G \mu_V}{\mu_G - \gamma_G V^*}\right)\left(\frac{\lambda_Z}{\mu_Z - \frac{b V^*}{\pi} \left(c + \frac{\lambda_G \mu_V}{\mu_G -\gamma_G V^*}\right)}\right)$ and $y_2(V^*) = \frac{\rho_0 r}{d+\rho_0 V^*}.$ Note that $y_1$ is increasing on $\left(0, \frac{\mu_G}{\gamma_G}\right),$ while $y_2$ is decreasing on $\left(0, \frac{\mu_G}{\gamma_G}\right)$ (we prove this statement in the following remark).

Thus, if a solution exists, then it must be unique. We prove the existence of a solution as follows. Since $0 < V^* < \frac{\mu_G}{\gamma_G},$ we set $V^* = 0$ into $y_1(V^*)$ and $y_2(V^*)$ and obtain

Now,

That is, $y_1(0) < y_2(0).$ Similarly, we set $V^* = \frac{\mu_G}{\gamma_G}$ to obtain

Thus, $y_2\left(\frac{\mu_G}{\gamma_G}\right) < y_1\left(\frac{\mu_G}{\gamma_G}\right).$ Since $y_2$ is decreasing on $\left(0, \frac{\mu_G}{\gamma_G}\right)$, $y_1$ is increasing on $\left(0, \frac{\mu_G}{\gamma_G}\right)$, $y_1(0) < y_2(0)$ and $y_2\left(\frac{\mu_G}{\gamma_G}\right) < y_1\left(\frac{\mu_G}{\gamma_G}\right),$ $y_1(V^*)$ and $y_2(V^*)$ meet at one point, $V^*.$ When $\mathcal{R}_0 > 1,$ Equation (3.10) has a unique root on $\left(0, \frac{\mu_G}{\gamma_G}\right).$ Hence, Equation (3.10) has a unique positive solution $V^*$ and consequently the system (3.3) has a unique infection equilibrium $E_1 = (T^*, {T_i}^*, V^*, Z^*, A^*, G^*).$ Now, we have proved the following theorem.

Theorem 2. There exists a unique infection equilibrium $E_1 = (T^*, {T_i}^*, V^*, Z^*, A^*, G^*)$ of the system (3.3) when $\mathcal{R}_0 > 1.$

Remark: Here, we prove that the function

is increasing on $\left(0, \frac{\mu_G}{\gamma_G}\right)$, while the function

is decreasing on $\left(0, \frac{\mu_G}{\gamma_G}\right),$ which was an integral part in proving Theorem (3.2). Observe that,

Now, write $y_1(V^*) = f(V^*) + g(V^*)h(V^*),$ where $f(V^*) = \frac{\delta}{p}\left(c + \frac{\lambda_G \mu_V}{\mu_G - \gamma_G V^*}\right),$ $g(V^*) = \frac{{\psi}_{01}}{p}\left(c + \frac{\lambda_G \mu_V}{\mu_G - \gamma_G V^*}\right)$ and $h(V^*) = \frac{\lambda_Z}{\mu_Z - \frac{b V^*}{p} \left(c + \frac{\lambda_G \mu_V}{\mu_G -\gamma_G V^*}\right)}.$ Then,

It is easy to see that all three derivatives $f', g'$ and $h'$ are positive on $\left(0, \frac{\mu_G}{\gamma_G}\right).$ Furthermore, the functions $g$ and $h$ themselves are positive on $\left(0, \frac{\mu_G}{\gamma_G}\right).$ Therefore,

implying that $y_1$ is increasing.

Theorem 3. If $\mathcal{R}_0 > 1$ the unique endemic equilibrium is locally asymptotically stable.

Proof. To determine the local stability of the unique infection equilibrium we compute the Jacobian of the system. We arrange the variables as $(T, T_i, V, Z, A, G)$ and with this ordering the Jacobian takes the form

This Jacobian clearly has one negative eigenvalue $\lambda_1 = -\gamma_A V^*-\mu_A$. The remaining eigenvalues are the solution of the following matrix equation

We expand the determinant by the third column. The resulting $x$ determinants can easily be expanded, and we obtain the characteristic equation:

Combining the first two terms and the last two terms and moving the negative term to the right hand side of the equation, we obtain the following characteristic equation.

In what follows we shall prove that this equation does not have roots with non-negative real parts. To do that we rewrite the equation in the form $\mathcal H(\lambda) = \mathcal G(\lambda)$ where $\mathcal H(\lambda) = \rho_0 T^* p$ and

Taking absolute value of $\mathcal G(\lambda)$, we notice that

We consider more carefully the absolute value of one of the terms in $[.]$. The other terms can be handled in a similar way. We have

Let $\lambda = x+yi$ with $x\ge 0$. Then,

Similarly, if $\lambda = x+yi$ with $x\ge 0$, we have

Hence $|\mathcal G(\lambda)| > (\delta+\psi_{01} Z^*)(c+\mu_V G^*)$. We notice that from the equilibrium equations for $T^*_i$ and $V^*$, we have that

Hence,

From the above discussion, we have that for $\lambda = x+yi$ with $x\ge 0$ $|\mathcal G(\lambda)| > |\mathcal H(\lambda)|$ for all $\lambda$ with non-negative real part. Hence the equation $\mathcal H(\lambda) = \mathcal G(\lambda)$ does not have roots with non-negative real parts. Hence, the infection equilibrium is locally asymptotically stable. This completes the proof.

4.   Parameter estimation, identifiability analysis and numerical simulations

It is a standard practice to estimate the model parameters by fitting the model to the experimental data. However, before fitting any model to data, there is an essential step, called "identifiability analysis" that must be addressed. Identifiability analysis is the process of determining whether model parameters can be derived from a given data set. Basically, it addresses two critically important questions regarding parameter estimation problem: (i) Is the within-host model structured to reveal its parameters from the given noise-free observations? (ii) How accurate are the parameters obtained by fitting the within-host model to noisy data? Here we use the term observation to refer to continuous smooth output and the data to refer measurements of observations at the discrete time points. Clearly, observations are noise-free and data is contaminated with noise due to measurement errors. Identifiability analysis is a crucial step in parameter estimation problem and is performed in two steps, each step answering one of the critical questions. Structural identifiability analysis is the first step and addresses the the first question. It investigates if the within host model is structured to reveal its parameters from the observations. Structural identifiability is a theoretical feature of the within-host model, therefore is a prior identifiability analysis and should be performed prior to estimating the parameters from the experimental data. The second step in identifiability analysis is practical identifiability which considers the actual experimental data with noise and should be performed after estimating the parameters. We will perform both structural and practical identifiability analysis for the selected within-host model (3.1).

4.1.   Structural identifiability of the selected within-host model (3.1)

In this study, the observations of the parameter estimation problem are the within-host model's state variables. Namely, the target cells $T(t, p)$, virus $V(t, p)$, albumin $A(t, p)$, and globulin $G(t, p)$. Using the experimental data, the objective is to estimate the parameters of the model which is denoted by $p = (r, d, \rho_0, \delta, \psi_{01}, \pi, c, \mu_v, \lambda_z, b, \mu_z, \lambda_A, \gamma_A, \mu_A, \lambda_G, \gamma_G, \mu_G)$. Before estimating the parameters, we need to study the structural identifiability of the model, which is defined as below.

Definition 1. Let $p$ and $\hat{p}$ be two distinct parameter sets of the within-host model (3.1). If

then we say that the within-host model (3.1) is structurally identifiable.

Simply put, the definition states that when the model parameters are different so are the observations. Therefore, two identical observations are only possible if the parameters are the same. There are several methods to study structural identifiability of compartmental ODE models, for an in depth review we refer ^[38]. We use Differential Algebra for Identifiability of Systems (DAISY) software to determine the identifiability of the within host model (3.1) ^[39]. DAISY gives the following parameter correlations;

The DAISY result (4.1) states that the parameters $(r, d, \rho_0, \delta, \pi, c, \mu_v, b, \mu_z, \lambda_A$, $\gamma_A, \mu_A, \lambda_G, \gamma_G, \mu_G)$ can be identified from the observations of the target cells, $T(t, p)$, viral load, $V(t, p)$, albumin $A(t, p)$ and globulin levels $G(t, p)$. But, we can not identify $\psi_{01}$ or $\lambda_z$ uniquely. We can only identify their product $\psi_{01} \lambda_z$. Since the recruitment of immune cells, $\lambda_z$ and the rate at which the immune cells attach the infected cells $\psi_{01}$ are correlated, infinitely many parameters will yield to identical observations of the target cells, $T(t, p)$, viral load, $V(t, p)$, albumin $A(t, p)$, and globulin levels $G(t, p)$ as long as the the product $\psi_{01} \lambda_z$ is constant. Therefore the within-host model (3.1) is not structured to reveal its parameters from the given observations. However, since we know the parameter correlation that leads to unidentifiability, we can fix either the recruitment of immune cells, $\lambda_z$, or the rate at which the immune cells attach the infected cells, $\psi_{01}$, to obtain a structurally identifiable model. We summarize the structural identifiability analysis in the following proposition.

Proposition 1. The within-host HIV protein model (3.1) is not structured to identify all its parameters from the observations of target cells, $T(t)$, viral load, $V(t)$, albumin, $A(t)$, and globulin, $G(t)$. If one of the parameters, $\lambda_z$ or $\psi_{01}$, is fixed, then the within-host model (3.1) becomes structurally (globally) identifiable.

First, we fix the immune cell production rate to $\lambda_z = 1$ cells/(mL $\times$ day) ^[37] to obtain a structurally identifiable model. Then, we estimate the parameters of the within-host model (3.1) using the experimental data. The estimated parameters are presented in Table 4 and the model predictions are shown in Figure 2.

4.2.   Practical identifiability of the selected within-host model (3.1)

A structurally identifiable model might not be identifiable in practice when noisy data is considered in determining the parameters. In structural identifiability analyses, experimental data is not considered at all. Our analysis shows that within host model is structurally identifiable (3.1) and this feature of the model is based on observations, which are smooth continuous functions of the model's state variables. Simply put, structural identifiability analysis assumes infinitely many noise-free data points, whereas practical identifiability uses noisy finitely many points. Therefore, we continue with the practical identifiability analysis of the within-host model (3.1) using Monte Carlo simulations. The Monte Carlo simulations performed are outlined below ^[40,41].

1) Set the estimated parameters given in Table 4 as the true parameters, $p_0$. Solve the within-host model (3.1) numerically with the true parameters.

2) Generate $M$ = 1000 synthetic datasets from a normal distribution whose mean is the model predictions at the discrete data points. Namely, generate synthetic data using the true parameters,

where each measurement error $E_i, E_j, E_k, E_l$ is normally distributed with zero mean and standard deviation $\sigma$. Figure 3 shows a single synthetic data generated for $\sigma = 0%, 5%, 20%$.

3) Fit the within-host model (3.1) to data generated in step 2 and obtain estimated parameters $p_j$ for each data set $j = 1, 2, \ldots M$.

4) Calculate the average relative estimation error (ARE) for each parameter by

where $p_0^{(k)}$ is the $k^{kt}$ parameter in $p_0$ and $p_j^{(k)}$ is the $k^{th}$ parameter in the set $p_j$.

5) Repeat steps 1 through 4 with higher noise levels. That is, at each iteration take $\sigma$ = 0%, 1%, 5%, 10%, 20%.

The practical identifiability performed with MCS is based on generating synthetic data using the estimated parameters. Therefore, it is a local result for the estimated parameters. Since, the true parameter set that generated the synthetic data is known, we can compute the average estimation error of each parameter. AREs describe how noise in the data affects the parameters. The AREs gradually increase as the noise level in data increases. The sensitivity of this increase determines the practical identifiability of the parameters. A practically unidentifiable parameter is very sensitive to noise in the data. In general, the AREs of those parameters will increase rapidly than the corresponding noise level. We introduce the following definition of practical identifiability.

Definition 2. Let $p$ be the parameter set of the within-host model (3.1). If the ARE of the $k^{th}$ parameter in $p$ is bounded by the measurement error, then we say that the parameter $p^{k}$ is practically identifiable. Specifically, if

then we say that the parameter $p^{k}$ is practically identifiable.

Examining the AREs presented in Table 5, and using Definition 2, we conclude that only the parameters, $r$, $d$, $\rho_0$, $\pi$, $c$, $\psi_{01}$ and $\mu_z$ are practically identifiable. All the other parameters have substantially high AREs and therefore are practically unidentifiable. AREs of parameters increase as the noise level in data increases, but the parameters $\mu_v$ and $\gamma_G$ have already very high AREs when there is only $1%$ noise in data. It seems like noise in the data is not the only source of high AREs for $\mu_v$ and $\gamma_G$. Even though the within-host model (3.1) is structurally identifiable, when $\lambda_z$ is fixed, we see that most of its parameters are practically unidentifiable. We suspect that the practical unidentifiably might be due to lack of enough data points. To test this hypothesis, we rerun the MCS with high frequency data points, where 10 measurements are collected each day for $170$ days, which results in total of $1700$ data points for each observations; viral load, target cells, albumin and globulin. AREs of the parameters with high frequency data is presented in Table 6. Based on the results in Table 6, we conclude that when there is enough data points, all the parameters of the within-host model (3.1) except $\mu_v$ and $\gamma_G$ are practically identifiable.

For the next numerical simulation, we studied the effect of varying albumin and globulin intake on the viral load. Both albumin and globulin intake were varied by two orders of magnitude. We select $\lambda_G$ at random from a uniform distribution $[0.005, 0.5]$ and $\lambda_A$ at random from a uniform distribution $[3, 300]$. The simulation results are shown in Figure 4. Clearly, albumin and globulin intakes have significant impact on the albumin and globulin levels. However, their impact on the viral load is minimal.

5.   Discussion

Albumin and globulin metabolism plays important role in the immune system and the immune system is of essence for protecting us from infectious diseases. Thus, there is an intrinsic link between albumin-globulin metabolism and infection diseases. Despite that extensive modeling has been done on infectious diseases and their interaction with the immune system, and significant work has been done on modeling nutrition and the dynamics of its components within-host, virtually no models have been introduced covering the interplay of infectious diseases and albumin/globulin. In this paper we consider a within host HIV model and the interaction of the virus with the protein in the body, namely the albumin and the globulins. We introduce a global model that incorporates target healthy and infected cells, viral load, immune cell response, and the concentrations of albumin and globulin. Because the impact of the two components of protein on the infection of healthy cells and the killing of the virus by the immune system are not known, we hypothesize a generic functional form for the terms. Furthermore, we assume a linear dependence of CD4 and CD8 cells recruitment on the globulins. We use first-order kinetics to characterize albumin and globulin dynamics. Even with this simplification, the within-host model (2.1) yields a total of 17 submodels. Had we incorporated more generalized terms to describe the dynamics of albumin and globulin, the number of models requiring fitting to the data would have increased substantially. We use data on CD4 cells, viral load, albumin and globulin concentrations in rhesus monkeys to determine the specific form of dependence of the incidence term and infected T cells clearance term on the protein concentrations, as well as the form of the recruitment terms of target cells and immune cells. The data points do not have any associated uncertainty and the noise level in the measurements is not known. We specify 17 submodels of the general model, fit all 17 models to the data, and use the Akaike Information Criterion (AIC) to select the best model. We find out that the simplest model fits the data best and solely explains the data. According to the best fitted model, albumin and globulins do not impact the new infections term, nor do they impact the clearance of infected T cells by the immune system term. Furthermore, the recruitment of CD4 and CD8 cells also does not depend on globulins. The best fitted model outperforms the remaining 16 models substantially and has a relative likelihood of 0.996, therefore explaining the data by itself ^[27].

We further analyze the best fitted model. We find out that it has a unique infection-free equilibrium, and we compute the immunological reproduction number of the infection $\mathcal{R}_0$. The infection-free state is locally stable if the immunological reproduction number is smaller than one. If the immunological reproduction number is greater than one, the infection-free equilibrium is unstable and there is a unique infection equilibrium. The unique infection equilibrium is locally stable whenever it exits.

We fit the best fitted model to the data and determine its parameters. Before that, we address the question whether the parameter estimation problem is well posed and we study the structural identifiability of the model relative to the given data sets. We find out that all parameters of the model are structurally identifiable except $\lambda_z$ and $\psi_{01}$ but their product is identifiable. Hence, we fix the immune cell production rate $\lambda_z$ at a value obtained from the literature. Fixing the value of $\lambda_z$ makes all parameters of the model structurally identifiable. Furthermore, we perform practical identifiability of the model using Monte Carlo Simulations. We find that with the given data only seven of the 16 parameters have Average Relative Estimation Errors (AREs) within an order of magnitude of the noise errors and are practically identifiable. We surmise that the reason for the lack of practical identifiability of most of the parameters is the low quality of the data. To test this hypothesis we perform Monte Carlo Simulations with synthetic data which are sampled a lot more frequently. We find out that with the high frequency synthetic data only the parameters $\mu_v$ and $\gamma_G$ are not practically identifiable and rest of the parameters are all identifiable. We conclude that collecting more consistent and more frequently sampled data can lead to more reliably estimated parameters.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the github repository, https://github.com/NecibeTuncer/HIV_and_Nutrtion.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Author Maia Martcheva is supported in part by NSF DMS 1951595, whereas Necibe Tuncer is supported in part by NSF DMS 1951626.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

Proof of Theorem 1: The local stability of equilibria is determined by the eigenvalues of the Jacobian computed at this equilibrium. The Jacobian of the model (3.3) at the equilibrium $E_0 = (T_0, 0, 0, Z_0, A_0, G_0)$ is

Hence the characteristic equation of the Jacobian is given by

Clearly, $-d, -\mu_Z, -\mu_A$ and $-\mu_G$ are eigenvalues of the matrix $J(E_0)$ and the remaining eigenvalues are given by the eigenvalues of the following matrix $J_1,$

This leads to the following quadratic characteristic equation:

Now, it can be seen that this equation has one positive root if $\mathcal{R}_0 > 1,$ and two negative real roots or two complex roots with negative real parts if $\mathcal{R}_0 < 1.$