Bayesian estimation of parameters in viral dynamics models with antiviral effect of interferons in a cell culture

1.   Introduction

The first works related to the use of mathematical models to describe the dynamics of the human immunodeficiency virus (HIV) were published in the 1990s; these works provided an initial understanding of the behavior of the virus. In 1995, Ho et al. described the dynamics of viral infection in patients infected with HIV type 1 (HIV-1) under antiviral therapy, specifically under treatment with protease inhibitors. Their studies indicate that this treatment causes the levels of HIV in plasma to decrease exponentially, and CD4$^+$ cells substantially increase ^[1]. A year later, Perelson et al. proposed a model for analyzing a detailed set of HIV-1 viral load data, using data collected from five infected individuals after the administration of a potent protease inhibitor. The model parameters were estimated using a nonlinear least squares method. In this study, the half-life of the infected cells, the average production of virions per day, the average length of the virus life cycle and other parameters were calculated that provided not only a first understanding of viral dynamics but also theoretical principles to guide the development of treatment strategies ^[2].

In 1996, Nowak et al. presented the basic model of viral dynamics, a simple model for the interaction between free virus particles and host cells. This model is composed of a system of nonlinear differential equations with three variables: target cells, infected cells and free virus particles ^[3]. Later, in 1997, Bonhoeffer et al. described the viral dynamics of HIV under pharmacological therapy by adding to the basic model a parameter corresponding to the efficacy of two therapies: reverse transcriptase inhibitors and protease inhibitors. The former prevents the infection of new cells, and the latter prevents infected cells from producing infectious virus particles ^[4]. In 2002, Wodarz and Nowak showed how mathematical models can be used to understand the dynamics of HIV infection and therapy. In their work, they described a basic model of virus infection and showed how it was used to gain some crucial insights into the dynamics during the asymptomatic phase of the disease. They explored how the evolution of HIV can drive disease progression and how mathematical models can be used to design specific treatments that can increase antiviral immunity and allow the long-term control of the virus ^[5]. In the same year, Callaway and Perelson conducted a comprehensive review of models of HIV dynamics, from the most basic model under antiviral therapy to more complex models that considered heterogeneities in drug efficacy. They found that the models that fail to robustly describe drug therapy include the basic model and its variants ^[6].

In 2002, Perelson provided an overview of models of HIV dynamics with and without antiviral therapy. He also discussed the dynamics of and interactions with the hepatitis C virus (HCV), hepatitis B virus (HBV), cytomegalovirus (CMV) and lymphocytic choriomeningitis virus (LCMV) ^[7]. In the literature, there are works related to other types of viruses, such as HCV and influenza, that are worth mentioning. In 1998, Neumann et al. studied the dynamics of HCV and the antiviral effect of interferons. They analyzed the viral decline in 23 patients during therapy with a mathematical model of differential equations with three variables, the target and infected cells and the viral load. The model used is an extension of the basic model in which the effect of interferon therapy is included. Their findings showed that HCV infections are highly dynamic and that the early monitoring of the viral load can help to guide therapy ^[8]. In 2007, Dahari et al. extended the model originally formulated by Neumann et al. in 1998 by including in it the proliferation of hepatocytes ^[9]. Baccam et al., in 2006, used a series of mathematical models of increasing complexity, incorporating the uninfected target cell limitation and the innate interferon response, to examine the kinetics of the influenza A virus in the upper respiratory tracts of experimentally infected adults. They used nonlinear regression and least squares estimation to fit models to a data set from a study of experimental swine flu virus infection ^[10]. In 2008, Beauchemin et al. analyzed the dynamics of an influenza A (H3N2) viral infection, using a set of mathematical models that highlight the differences between in vivo and in vitro infection ^[11]. In 2022, Wu, He and Khan ^[12] investigated the effect of administering treatments to individuals with HIV and considered an optimal control problem subject to multiple drug treatments for the within-host model.

In this work, three types of population dynamics of healthy cells are presented. The first is von Bertalanffy growth ($S^{\prime}(t) = r \big(k-S(t) \big)$), the second is logistic growth ($S^{\prime}(t) = r S(t) \Big[1-\frac{S(t)}{k} \Big]$) and finally the third is the novel Gompertz growth ($S^{\prime}(t) = r S(t) \ln\Big[\frac{k} {S(t)}\Big]$), where the parameters $r$ and $k$ represent the growth rate of the cells and the load capacity (maximum number of cells), respectively.

In most viral infection models that are based on the basic model ^[3], supplementary effects are added to the viral dynamics phenomenon. The following system includes interferon therapy:

Three functions $F$ describe the growth of healthy cells:

● $F(S, I): f_1 = r(k-S(t))$ ^[3,13],

● $F(S, I): f_2 = r S(t)\big[1-\frac{S(t)+I(t)}{k} \big]$ ^[14,15] and

● $F(S, I): f_3 = r S(t) \ln \big[\frac{k}{S(t)+I(t)}\big]$.

The function $f_1$ describes bounded-monotone type growth, and the functions $f_2$ and $f_3$ describe sigmoid type growth, but one of the most important differences is that the logistic curve is symmetric, whereas the Gompertz curve is asymmetric.

The parameters $\beta$, $\delta$, $\rho$ and $\gamma$ represent the infection rate, the death rate of infected cells, the rate of virus production and the viral clearance rate, respectively. $\eta$ is the efficacy of interferon therapy in preventing new infections, and $\epsilon$ is the efficacy of interferon therapy in inhibiting viral production, where $0\leq\eta, \epsilon\leq 1$. An efficacy of $0$ indicates that there is no inhibition, whereas an efficacy of $1$ ($100$$\%$) indicates complete inhibition. Values of the efficacy between $0$ and $1$ indicate partial inhibition.

This work aims to estimate the parameters of the efficacy of interferon therapy in in vitro virus infection dynamics models under the Bayesian approach. There are four specific objectives: first, to estimate the parameters of cell dynamics models; second, to estimate the parameters of viral dynamics in models of viral infections; third, to estimate the parameters of the efficacy of interferon therapy in models of viral infections; and fourth, to use the Watanabe-Akaike information criterion (WAIC) to select the model that best fits the cell and viral growth data.

2.   Key results from stability analysis

The three models are studied in the non-negative octant

2.1.   Stability analysis of system (1.1) with $F(S, I)=f_1$

For the basic model that corresponds to system (1.1) with $F(S, I) = f_1$, a local and global stability study was carried out to find equilibrium points, the infection-free equilibrium $E^0 = (k, 0, 0)$ and the infected equilibrium $E^* = (S^*, I^*, V^*)$ with coordinates

where the basic reproductive number is

Theorem 2.1. ^[16,17]. Consider system (1.1) with $F(S, I) = f_1$; the following conditions hold:

a) The infection-free equilibrium $E^0$ of the model is asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1$.

b) If $R_0 > 1$, then there exists a unique infected equilibrium $E^*$ in the model, which is asymptotically stable.

The proofs of the local and global stability can be found in ^[16,17].

2.2.   Stability analysis of system (1.1) with $F(S, I)=f_2$

In the work of Ikeda et al. ^[15], this model is presented, and it has three equilibrium points: the vanishing equilibrium $E_0 = \left(0, \ 0, \ 0\right)$, the only infection-free equilibrium $E^0 = \left(k, \ 0, \ 0\right)$ and the infected equilibrium $E ^* = \left(S^*, \ I ^*, \ V^*\right)$. The coordinates of the infected equilibrium are as follows:

where the basic reproductive number is the same as that defined by Eq (2.1). A qualitative analysis of the model was performed by Dingli et al. ^[14].

Theorem 2.2. ^[14]. Consider system (1.1) with $F(S, I) = f_2$; the following conditions hold:

a) The infection-free equilibrium $E^0$ of the model is asymptotically stable if $R_0 < 1$ and unstable if $R_0 > 1$.

b) If $R_0 > 1$, then there exists a unique infected equilibrium $E^*$ in the model, which is asymptotically stable if $\frac{r}{k}\left[\left(\frac{r}{k}S^*+\delta\right) \left(\delta+\gamma+\beta V^*\right)+\left(\beta\rho S^*+\gamma^ 2\right) \right] > \rho\beta\beta V^*$.

2.3.   Stability analysis of system (1.1) with $F(S, I)=f_3$

We can easily see that, for all the parameter values, the vanishing equilibrium $E_0 = (0, 0, 0)$ and the infection-free equilibrium $E^0 = (k, 0, 0)$ always exist.

Theorem 2.3. Let $(S(t), I(t), V(t))$ be the solution of system (1.1) satisfying conditions $(S_0, I_0, V_0)\in\mathbb{R}^3_{0+}$. Then $S(t)$, $I(t)$ and $V(t)$ are all bounded for all $t \geq 0$ at which the solution exists.

The basic reproductive number of system (1.1) with $F(S, I) = f_3$ is the same as that defined in equation (2.1).

Theorem 2.4. Consider system (1.1) with $F(S, I) = f_3$; the following conditions hold:

a) If $R_{0} < 1$, $E^{0}(k, 0, 0)$ is asymptotically stable; if $R_{0} > 1$, $E^{0}(k, 0, 0)$ is unstable.

b) If $R_{0} > 1$, then system (1.1) has a unique infected equilibrium $E^{\ast}(S^{\ast}, I^{\ast}, V^{\ast})$.

c) Assume that $R_{0} > 1$ and that

then, the infected equilibrium $E^{\ast}$ is asymptotically stable for system (1.1) when $F = f_3$.

The proofs of the theorems and basic reproductive number calculations for this subsection can be found in the Appendix.

3.   Bayesian estimation of three models

To estimate the antiviral effect of the interferon IFN-alpha 2b on HIV-1 NL-AD8 replication, three data sets are used. The first consists of the kinetics of MT4C5 cells, and three growth models are used: the von Bertalanffy $[f_1]$, logistic $[f_2]$ and Gompertz $[f_3]$ models.

For the analysis of the model with and without interferon therapy, viral load data were used. The data correspond to the interactions of MT4C5 cells and the HIV-1 NL-AD8 in the absence and presence of the interferon IFN-alpha 2b obtained from ^[15].

3.1.   Methodology

The Bayesian estimation of the parameters of a model for each data set comprises three stages:

● First stage: Propose non-informative prior distributions and estimate the parameters of the cell dynamics $(r, k)$ for the first set of data, which is shown in Figure 1a.

● Second stage: For the second data set, estimate the parameters of the cell dynamics $(r, k)$ using the posterior distributions obtained from the first stage as the prior distributions in this stage. Propose non-informative prior distributions for the estimation of the viral dynamics parameters $(\beta, \delta, \rho)$ (Figure 1b), with $\gamma = 2.3$ ^[15] in this stage and in the next stage.

● Third stage: For the third data set, estimate the parameters of the cell dynamics $(r, k)$ and viral dynamics $(\beta, \delta, \rho)$ using the posterior distributions obtained in the second stage as prior distributions in this stage. Propose non-informative prior distributions for the estimation of the interferon efficacy parameters $(\eta, \epsilon)$ (Figure 1b).

3.2.   Statistical model

The following statistical model is considered:

● $\varepsilon(t_i)$ is the random error at time $i$; these errors are independent, identically normally distributed variables with zero mean and a variance $\sigma^2$.

● $y_i$, which is used to set each stage, represents the number of healthy cells at time $i$ for the first data set, and for the second and third data sets, it represents the viral load measured at time $i$.

● $\boldsymbol{\theta}$ is used to set each stage of the model: $\boldsymbol{\theta}$ = $(r, k)$ corresponds to cell growth models. For the second stage, $\boldsymbol{\theta}$ = $(r, k, \beta, \delta, \rho)$ corresponds to the model without interferon therapy, and for the third stage, $\boldsymbol{\theta}$ = $(r, k, \beta, \delta, \rho, \eta, \epsilon)$ corresponds to the model with interferon therapy.

● $X_{\theta}(t_i)$ is the numerical solution calculated with the Runge-Kutta method for cell growth models for stage 1, and it represents the equation corresponding to the viral load variable ($V(t)$) of system (1.1) for the second and third stages, respectively.

3.3.   Likelihood function

As a consequence of the fact that in model (3.1), it was assumed that the errors are independent and $\varepsilon(t_i)\sim N(0, \sigma^2)$, $y_i \sim N(X_{\theta}(t_i), \sigma^2)$, and the $y_i$ are also independent.

Therefore, the likelihood function of the data under the statistical model (3.1) is given by

3.4.   Posterior distribution

Bayesian inference about a parameter vector $\boldsymbol{\theta}$ is based on the posterior distribution $\pi(\boldsymbol{\theta} | \boldsymbol{D})$, which is the probability distribution of the parameter vector $\boldsymbol{\theta}$ given the data $\boldsymbol{D}$. That is,

$\pi(\boldsymbol{\theta})$ is the prior distribution, $\pi(\boldsymbol{D} | \boldsymbol{\theta})$ is the data probability distribution or likelihood function, and $\pi(\boldsymbol{D}) = \int_{}^{} \cdots \int_{}^{} \pi(\boldsymbol{\theta}) \pi(\boldsymbol{D} | \boldsymbol{\theta}) \, d \boldsymbol{\theta}$ is a normalization constant that does not depend on $\boldsymbol{\theta}$, so that Equation (3.2) becomes:

In the prior distribution, $\pi(\boldsymbol{\theta})$, is expressed the initial knowledge about $\boldsymbol{\theta}$. When there is no information about $\boldsymbol{\theta}$, non-informative priors are used that give equal probabilities to each possible value of $\boldsymbol{\theta}$, otherwise informative priors are used that give a higher probability to some $\boldsymbol{\theta}$ values.

To perform Bayesian inference it is necessary to find the mean of the posterior distribution

In practice, these integrals are difficult or impossible to obtain. One method of approximating the posterior distribution and thereby obtaining the mean of the posterior distribution is to use Monte Carlo Markov chains (MCMC) which are a special case of a discrete stochastic process in which the probability of an event depends only on the immediately preceding event.

There are several methods to generate MCMC, the most used are the Metropolis-Hastings algorithm, Gibbs Sampling and Hamiltonian (or hybrid) Monte Carlo (HMC) methods.

3.5.   Hamiltonian Monte Carlo (HMC) technique

We use the Hamiltonian Monte Carlo (HMC) technique. The HMC technique uses the log posterior gradient to steer the Markov chain to regions of higher posterior density, where most samples are located. As a result, a well-tuned Markov chain created using the HMC algorithm has a much higher acceptance rate than the traditional Metropolis algorithm. The Metropolis algorithm has an acceptance rate of approximately 23%, while the HMC method has an acceptance rate of 65% ^[18]. In addition to better exploring the posterior distribution, the HMC algorithm is more efficient than the Metropolis algorithm and does not suffer from the Gibbs parameter mapping problem.

3.6.   Convergence of the chains

We use the following convergence diagnostic criteria for the MCMC algorithm:

● An empirical approach for convergence control consist of two graphics, trace and density, of the simulated chain output to detect non-stationary behaviors ^[19]. The traceplot is a time series plot of the iteration number against the realizations of the Markov chain at each iteration. Convergence is declared when there is a good mixing across chains. The density plot is a non-parametric estimate of its density for each chain. Convergence is declared when the densitys are similar.

● The Gelman-Rubin‘s diagnostic compares two estimations of the variance of the chains, between and within. It gives two statistics: the potential scale reduction factor (PSRF) and its credible interval. Convergence of the chains is not rejected when the PSRF and the upper limit of its credible interval are close to one and less than 1.2, respectively ^[20].

3.7.   Bayesian estimation

The point Bayesian estimators were the means of the posterior distributions corresponding to minimizing the expected squared error loss function. The high posterior density intervals (HPDIs) were also calculated.

3.8.   Statistical software

The estimation of the parameters models of the first-stage was performed using the library brms ^[21] of the statistical software R ^[22]. For the second and third stages, the library turing ^[23] of the Julia software ^[24] was used. For each model, three chains were generated, each with 10,000 iterations and a burning of 5000 samples.

In order to evaluate the convergence of the MCMC chains, the Gelman-Rubin diagnostics are utilized; they are implemented in the coda ^[25] library of the R software.

To obtain the hyperparameters of the posterior distributions of the parameters corresponding to the first $(r, k)$ and second $(\beta, \delta, \rho)$ stages, we used the R library fitdistrplus ^[26] to perform the fitting procedure.

3.9.   Model comparison

Model performance is evaluated with the Watanabe-Akaike information criterion ($WAIC$). This is a fully Bayesian approach to estimating the expected log predictive density and is defined as follows ^[18]:

where $lppd$ is the log pointwise predictive density and $pWAIC$ is the effective number of parameters. The rule of thumb is that among two models the best is the one with a lower WAIC value.

In practice, the $lppd$ is calculated using samples of $p_{post} (\theta)$, the usual posterior simulations, which we label $\theta^s$, $s = 1, ..., S$:

where $S$ is the number of simulation draws, which is large enough to capture the posterior distribution.

The calculation of the effective number of parameters ($pW AIC$) is performed as follows:

where $V_{s = 1}^{S} \log \, p(y_{i}\mid \theta^{s})$ is the posterior variance of the log predicted density for each $y_{i }$, and $V_{s = 1}^{S} a_s = \dfrac{1}{S-1} \sum\limits_{s = 1}^{S} (a_s - \overline{a})^2$ is the sample variance.

3.10.   Infection-related quantities

In addition to estimating the parameters of the models, it is also relevant to compute other values related to a viral infection, such as the basic reproductive number ($R_0$) and the half-life of infected cells. The basic reproductive number ($R_0$) of the three models is defined in Eq (2.1); we calculated this value using the posterior distribution of each parameter. Similarly, the calculation for the half-life of infected cells ($(\ln{2})/\delta$) is performed.

3.11.   Prior distributions

The prior distributions are non-informative, are commonly used in similar problems and are the same for the three models; they are shown in Table 1. Table 2 shows the posterior distributions of the first stage and those proposed for the viral dynamics parameters. Table 3 shows the posterior distributions of the second stage and those proposed for the viral dynamics parameters with interferon therapy.

Remark 1. Prior distribution for $\sigma^2$ is the same in the three models (1.1) and three stages, non-informative $Inverse-Gamma(0.001; 0.001)$. $r$ and $k$ also have non-informative prior distributions (Table 1), but they have Gamma prior distribution because posterior distribution in the first stage motivates this election (Table 2) and second stage (Table 3), respectively. Prior distributions' hyperparameters are estimated using the posterior distributions. $\beta$, $\delta$, and $\rho$ have non-informative prior distributions (Table 2) in the second stage for the three models (1.1), but in the third stage they have Gamma prior distribution because posterior distribution in the second stage motivates this election and its hyperparameters are estimated using this information. Finally, $\eta$ an $\epsilon$ have non-informative prior distributions (Table 3).

4.   Results

Table 4 presents the results of the estimations of the parameters of the cell growth models. Figure 2 shows the numerical solution of the estimation results for the von Bertalanffy (), logistic () and Gompertz (-) cell growth models.

Table 5 presents the results of the second stage of the three models based on model (1.1), when no therapy is used. Figure 3 shows the three solutions for the second stage with the Bayesian parameters using two axes: the first axis is used to represent healthy (red) and infected (blue) cells, and the second axis is used to represent the viral load (black). Different types of lines were used to differentiate the von Bertalanffy (continuous line), logistic (point-line) and Gompertz (point-point) models and the experimental data (circles).

Table 6 presents the results of the third stage for the three models based on model (1.1), when interferon therapy is used. Figure 4 shows the three solutions for the third stage with the Bayesian parameters using two axes: the first axis is used to represent healthy (red) and infected (blue) cells, and the second axis is used to represent the viral load (black). Different types of lines were used to differentiate the von Bertalanffy (continuous line), logistic (point-line) and Gompertz (point-point) models and experimental data (circles).

4.1.   Model selection

The best fitting model in the first stage is the Gompertz model with a WAIC = 258.1, in second place is the von Bertalanfy model with a WAIC = 259.2, and in third place, is the logistic model with WAIC = 264.0. This gives a measure of the improvement in the model fitting by changing the type of cell growth model.

In the second stage, the logistic model had a WAIC = 144.6, which was the best WAIC value in this stage; the Gompertz model remained in second place with a WAIC = 167.8, and finally, the von Bertalanffy model had a WAIC = 171.8.

In the third stage, the logistic model remained the best model, with a WAIC = 87.6. In this stage, the Gompertz model was in second place, with a WAIC = 91.1, and the model with the worst fit was the von Bertalanffy model, with a WAIC = 99.9.

4.2.   Infection-related quantities

The Bayesian estimations of $R_0$ are greater than 2 for all models; however, for the models with interferon therapy, $R_0$ is smaller than it is for the models without interferon therapy in all cases. Graphically, it can be seen that the credible intervals do not overlap, and therefore the interferon effect is statistically significant. That is to say, interferon therapy significantly reduces the basic reproductive number (Figure 5).

The Bayesian point estimations of the half-life of infected cells are very similar in all models; the half-life is around 0.45 days for all models, except for the Gompertz model with interferon therapy (system (1.1) with $F(S, I) = f_3$), which has a half-life of slightly more than 0.5 days (Figure 6). Graphically, it can be seen that the credible intervals overlap, which means that interferon therapy has no significant effect on the lifespan of infected cells; this is because the death rate of infected cells is modeled in the same way in all models.

4.3.   Convergence of the selected model

The comparison of the models using the WAIC showed that the best model was the model with logistic growth for the second and third stages. The convergence results of the method for the three models are similar, so we opted to present the convergence of the best model. In Figures 7–9, a good mixture of the chains can be observed, and they also cover the entire parameter space of their posterior distribution, indicating convergence. The results of the Gelman-Rubin diagnostic for the three stages showed a value of less than 1.1 of the PSRF for each parameter and an upper limit of its credible interval below 1.2 in all cases; this indicates that the chains have converged to the stationary distribution (see Table 7). The convergence of the remaining models (Figures A1–A6) is discussed in Sections A.6 and A.7.

5.   Conclusions

In this work, we estimate the parameters of the cell dynamics and viral dynamics without and with interferon therapy using the Bayesian estimation approach, which has the advantage of incorporating information from the parameters of the previous stages and the obtained interval estimates.

According to the WAIC, the fits of the models to the data differ, but the WAIC values are very close, leading to the conclusion that the three considered models (von Bertalanffy, logistic and Gompertz) present a good fit with the cell growth data.

For virus dynamics models without and with interferon therapy (second and third stages, respectively), the model with logistic growth best described the viral load data, and the model with the worst fit was the model with von Bertalanffy growth.

The following conclusions are obtained from the models of viral dynamics with logistic growth. The estimated inhibition effect of interferon therapy on new infections was approximately 9.9%, and with probability of 0.95, this effect is between 0.4 and 26.5%. The estimated inhibition effect of interferon therapy on viral particle production was approximately 8.1%, and with a probability of 0.95, this effect is between 0.3 and 23.1%. This indicates that these inhibition effects are similar for the prevention of new infections and the prevention of viral replication. The estimated basic reproductive number with interferon therapy of 2.71 (95% $CrI$: 2.37–3.11) was smaller than the basic reproductive number without interferon therapy, which was 5.06 (95% $CrI$: 4.17–6.51); it can be concluded that interferon therapy significantly weakens the viral infection and hence the spread of the infection.

The point estimate of the half-life of infected cells is very similar in all models (see Figure 6). Graphically, it can be seen that the credible intervals overlap, which means that interferon therapy does not have a significant effect on the lifetime of the infected cells. This is because interferon therapy does not modify the mortality rate of infected cells.

We estimated the efficacy of interferon therapy and the basic reproductive number; we have concluded that interferon therapy reduces the basic reproductive number, but not below unity.

The Bayesian estimation approach incorporates information from the parameter in the prior distribution. In our case, we do not have much information about which values were more probables in the first stage for parameters $r$ and $k$. But in the second stage we use the information obtained in the first stage. Bayesian estimation is better than frequentist estimation when the sample size is small because in this case the maximum likelihood estimator, used in frequentist estimation, does not have normal distribution, and Bayesian estimators do not need to have a specific distribution. Generally, the least-squares method has been used to estimate the parameters in the viral dynamics models. A disadvantage of this method is that only point estimates are obtained. The Bayesian approach does not present this drawback since it allows determining the distributions of the parameters of interest and, therefore, not only point estimates are obtained but the credibility intervals estimates. In addition, estimation for functions of the parameters of interest, such as the basic reproductive number ($R_0$) and the half-life of infected cells is easier.

Finally, comparing our results of the virus dynamics model with logistic growth with those obtained by Ikeda et al.^[15], we note that the estimates of the cell dynamics parameters $r$ (0.610 vs 0.630) have the same order of magnitude, although our estimate of the carrying capacity $k$ is approximately half (1,447,204 vs 3,330,000) that obtained in ^[15]. We obtained similar findings for the parameters of viral dynamics: $\beta$ (0.0012 vs 0.0067), $\delta$ (1.6399 vs 1.5600), and $\rho$ (0.0099 vs 0.0026). As for the estimation of the parameters of the efficacy of the interferons, we have similar results for antiviral effect of IFN on de novo infection $\epsilon$ (0.0812 vs 0.06) but different ones in the antiviral effect of IFN on virus production $\eta$ (0.0993 vs 0.47). For the infections related quantities, we obtain similar results for the half-life of infected cells (0.44 vs 0.45), but different results in the $R_0$ without interferon therapy (5.06 vs 13.1) and $R_0$ with interferon therapy (2.71 vs 6.37), our estimates of the $R_0$ are approximately half those obtained in ^[15]. This difference in the estimate of $R_0$ may be due to the difference in the estimate of the parameter $\eta$ and the use of a mores robust estimation method, the Bayesian approach.

Acknowledgments

This work was partially supported by the CONACYT-Ciencia de Frontera "Ecuaciones de evolución, sus estados estacionarios y su comportamiento asintótico con aplicaciones en Física y Biología" (project No. 684340). Rodríguez-Parra is indebted to CONACYT for the fellowship that enabled him to pursue graduate studies for the degree of Maestría en Matemáticas Aplicadas.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

A.1.   Positive invariance and boundedness of system (1.1) with Gompertz growth

In this subsection, we provide the proof for Theorem 2.3.

Proof [Theorem 2.3]. Let $(S(t), I(t), V(t))$ be any solution with initial conditions $(S_0, I_0, V_0)\in\mathbb{R}^3_{0+}$. We define a function

The time derivative along a solution of system (1.1) is

Note that $\ln(1 - z) \leq -z$ holds for any $z < 1$; one can obtain

It follows that

where $\zeta = min\{r, \frac{\delta(n-1)}{n}, \gamma\}$. Thus, $\underset{t\longrightarrow \infty }{\lim \sup B(t)} \leq (rk)/\zeta$.

Therefore, $S(t)$, $I(t)$ and $V(t)$ are all bounded for all $t\geq 0$. This completes the proof.

A.2.   Basic reproductive number of system (1.1) with Gompertz growth

The dynamics of virus infections crucially depend on the basic reproductive number $R_0$, which in this case is the average number of secondary infected cells produced by a single infected cell in a population of target cells. Here, we derive the basic reproduction number for viral infection using the next generation operator ^[27]. Using the notation in ^[27], the non-negative matrix $F$ of the new infection terms and the M-matrix $V$ of the transition terms associated with the model based on system (1.1) with Gompertz growth are given, respectively, by

It follows that the basic reproduction number, denoted by $R_0 = \rho(FV^{-1})$, where $\rho$ is the spectral radius, is given by

where $S^0$ is the concentration of uninfected cells at the uninfected equilibrium in $\mathbb{R}^3_{0+}$, $E(S^0, 0, 0)$, where $S^0 = k$.

The parameter $R_0$ has an interesting biological meaning; it is the average number of secondary infected cells produced by a single infected cell in a population of target cells. Particularly, the term $\frac{(1-\eta)\beta(1-\epsilon)pS^0}{\delta\gamma}$ describes the secondary infections per unit time through the direct target-cell-to-virus contact with the antiviral effect of interferon therapy.

A.3.   Asymptotic stability of the uninfected equilibrium

In this subsection, we prove part a) of Theorem 2.4. The Jacobian matrix of system (1.1) with Gompertz growth at $E^{0}$ is

One of the characteristic roots is

The other two eigenvalues of $J(E^{0})$ are determined by the quadratic equation

The characteristic equation (A6) has two eigenvalues that have negative real parts if and only if $R_{0} < 1$. Therefore, $E^{0}$ is asymptotically stable for $R_{0} < 1$ (By Theorem 4.6 ^[29]). If $R_{0} > 1$, equation (A6) has a positive eigenvalue. Therefore, $E^{0}$ is unstable.

A.4.   Existence and uniqueness of infected equilibrium

In this subsection, we prove part b) of Theorem 2.4. To do this, we explore the existence and uniqueness of the infected equilibrium of the model based on system (1.1) with $F = f_3$, following an argument similar to that given in ^[28]; we use a geometric approach. The equilibria are obtained by setting the right-hand side of system (1.1) with $F = f_3$ equal to zero:

The second equation leads to $S^{\ast} = \delta I^{\ast}/(1-\eta)\beta V^{\ast} = \delta\gamma/ (1-\epsilon)p (1-\eta)\beta$, and the third equation leads to $V^{\ast} = (1-\epsilon) p I^{\ast}/\gamma$.

We express the first equation of system (A7) as the following transcendental equation in $I^{\ast}$:

The number of $I$ solutions of Eq (A8) can be analyzed geometrically through the intersection points of the logarithmic function $f_1(I)$ with the equilateral hyperbola $f_2(I)$. The functions $f_1(I)$ and $f_2(I)$ are defined as

Observe that $f_1(0) = \ln R_0,$ $f_1^{\prime}(I) < 0$ and $f_1^{\prime\prime}(I) > 0$, and it is well known that $f_2^{\prime}(I) > 0$ for positive "$I$" values. If $\ln R_0 > 0$ (or equivalently $R_{0} > 1$), $f_1(I)$ and $f_2(I)$ have a single positive intersection point in the first quadrant. It follows that the infected equilibrium exists, and it is unique if $R_{0} > 1$.

A.5.   Asymptotic stability of the infected equilibrium

In this subsection, we prove part c) of Theorem 2.4. To do this, we study the local stability behavior of the infected equilibrium $E^{\ast}$ for system (1.1) with Gompertz growth.

The local stability of $E^{\ast}$ is given by the Jacobian matrix of system (1.1) evaluated in this equilibrium:

which can be rewritten as

when we take into account the following identities:

If $C\geq0$, then the infected equilibrium is asymptotically stable (By Theorem 4.6 ^[29]). Note also that $C\geq0$ is equivalent to $\ln k\leq\frac{S^{\ast}}{S^{\ast}+I^{\ast}}+\ln\left[S^{\ast}+I^{\ast}\right]$. We provided a sufficient condition on the parameters for the local stability of the infected equilibrium.

A.6.   Convergence of the MCMC simulations of the models with von Bertalanffy growth

In this subsection, the convergence criteria of the von Bertalanffy growth model and viral infection models with von Bertalanffy growth are discussed.

The chains show good mixing, as shown in Figures A1–A3, which display the results of the von Bertalanffy model. Table A1 shows the results of the calculation of the Gelman-Rubin criterion, thus indicating that both convergence results are consistent.

A.7.   Convergence of the MCMC simulations of the models with Gompertz growth

In this subsection, the convergence criteria of the Gompertz growth model and viral infection models with Gompertz growth are discussed.

The chains show good mixing, as shown in Figures A4–A6, which display the results of the Gompertz model. Table A2 shows the results of the calculation of the Gelman-Rubin criterion, thus indicating that both convergence results are consistent.