Stochastic Modelling of HIV-1 Replication in a CD4 T Cell with an IFN Response

A mathematical model of the human immunodeficiency virus Type 1 (HIV-1) life cycle in CD4 T cells was constructed and calibrated. It describes the activation of the intracellular Type I interferon (IFN-I) response and the IFN-induced suppression of viral replication. The model includes viral replication inhibition by interferon-induced antiviral factors and their inactivation by the viral proteins Vpu and Vif. Both deterministic and stochastic model formulations are presented. The stochastic model was used to predict efficiency of IFN-I-induced suppression of viral replication in different initial conditions for autocrine and paracrine effects. The probability of virion excretion for various MOIs and various amounts of IFN-I was evaluated and the statistical properties of the heterogeneity of HIV-1 and IFN-I production characterised.


Introduction
Human immunodeficiency virus Type 1 (HIV-1) infections lead to the decline of immune functions followed by disease progression to AIDS [1]. Antiviral immune responses to HIV-1 are characterised by innate and adaptive immunity [2]. Studying of the host and viral factors that control viral replication and transmission is important to better understand the robustness and fragility of virus persistence.
The interferon Type I (IFN-I) response of innate immunity plays a central role in inhibiting viral replication [3]. The underlying processes include recognition of HIV-1 by innate immune sensors within a cell, induction of IFN-I production, and finally, the expression of antiviral genes (i.e., interferon-stimulated genes (ISGs)). However, the life cycle of HIV-1 also leads to the production of proteins (Vpu, Vif) that inhibit the effect of the antiviral proteins APOBEC3, Tetherin, and SAMHD1, resulting in resistance to IFN-I [4]. While the pathways regulating the production of IFN-I and the induction of ISGs are well studied [3], the quantitative kinetic description of these processes still needs to be developed.
The production of IFN-I starts at the initial infection and persists during disease progression [5]. However, this innate immune response is insufficient to suppress HIV-1 replication as the virus is able to evade the IFN-mediated antiviral activity. Sustained IFN-I expression bears deleterious effects on host immunity [6][7][8]. The mechanisms of the virus escaping intracellular defences are not completely understood [9][10][11] because they include multiple layers of control of IFN-I signalling [12] and diverse roles of interferon-stimulated gene (ISG) products [13].
We previously formulated a mathematical model describing the life cycle of HIV-1 in CD4 T cells [14]. In the present work, we extend the above model to consider the production of IFN-I by an infected cell, the ISG-mediating induction of effector molecules exhibiting anti-viral effects, and the inhibiting functions of viral proteins, as proposed in [15]. Our study of the stochastic effects during the HIV-1 life cycle suggests that the inherent randomness of reactions and the discrete nature of the biochemical stages of HIV-1 infection need to be considered in order to properly characterise the variability in viral production between infected cells. Furthermore, taking IFN-I into account and the proteins mediating viral resistance adds an additional layer of regulation for the heterogeneity of HIV-1 replication in a single cell. A dynamical model simulating the interaction of HIV-1 with IFNα at a cell population level was developed in [16]. A data-based single-cell model of the IFN-I response to Dengue virus infection was elaborated in [17].
The aim of the current work was to formulate and calibrate a stochastic model of HIV-1 replication that includes antiviral responses of an infected cell and to apply it to examine the dependence of virus production and IFN-I secretion on antiviral defence reactions. This requires the application of computational modelling tools, which corresponds to the mainstream trends in immunology and virology [12,18]. Specifically, we aimed to: • Develop an analytical tool for data assimilation and analysis of HIV-1 replication and the IFN-I response; • Identify the control points underlying the viral ability to evade IFN-I-mediated defences as they represent potential targets for antiviral and immune therapies; • Predict the effective reproduction number of single-cell infection resulting from the competition of multiple factors related to viral and IFN-I-dependent products.
In Section 2, we present a biochemical scheme of the processes that regulate HIV-1 production and IFN-I reactions. The equations of the deterministic model are described in detail and translated to a stochastic model using the Gillespie approach. The results of extensive numerical experiments with the model are presented in Section 3. The discussion and conclusions of the study can be found in Sections 4 and 5, respectively.

Intracellular Type I Interferon Response and HIV-1 Replication
Upon HIV-1 infection of human cells, the virus is recognised by cytosolic sensors (RIG-1), which trigger a cascade of reactions, leading to the synthesis and activation of the IRF-3 and NF-kB proteins. The proteins translocate to the nucleus and cause the transcription of genes synthesising IFN-I [3,17,19].
Interferons are released into the extracellular space (IFNe) and bind to receptor molecules that are expressed on the cell surface of the produced cell or neighbouring cells. Binding of IFN-I to IFN-I receptors induces a signalling cascade with the participation of STAT-1 and STAT-2 proteins forming a transcription factor by assembling with IRF9. This translocates to the nucleus and initiates the transcription of ISGs. The following translation of the transcribed genes leads to the production of interferon-stimulated proteins such as APOBEC3, SAMHD1, and Tetherin. These proteins inhibit various stages of HIV-1 replication in the cell: APOBEC3 and SAMHD1 prevent reverse transcription, whereas Tetherin prevents the exit of the mature virion from the cell [19][20][21]. The action of ISGs is inhibited by the accessory proteins Vpu and Vif of HIV-1. These two viral proteins form complexes with the IFN-I-simulated antiviral proteins and target them for degradation [19,22]. Note that, as a result of the interaction of interferon-stimulated and HIV-1 accessory proteins, the amount of proviral DNA generated after reverse transcription decreases by about 50% in the presence of APOBEC3 [20] and by another 50% in the presence of SAMHD1 [23]. As observed in [21], Tetherin reduces the release of free virions by 40%.

Governing Deterministic Equations
The biochemical scheme if HIV-1 replication and the IFN-I response in a single cell is presented in Figure 1. The HIV-1 replication cycle consists of multiple stages, as described below. Based on this scheme, we developed a mathematical model of the HIV-1 life cycle in a CD4+ T cell. The major set of specific biochemical reactions and model parameters were taken from our previous work [14]. These were supplemented by a set of reaction stages for intracellular IFN-I synthesis and IFN-receptor-transduced responses of ISGs in an infected cell. Entry: Free viral particles V free are bound with the cell membrane according to the following equations:ẋ (1) where: ] is the number of free virions in vicinity of the cell; ] is the number of virions bound to CD4 and the co-receptor.

Reverse transcription and integration:
After the release of RNA following the entry of a virion into the cell, reverse transcription is initiated. Proteins SAMHD1 and APOBEC3 have an inhibitory effect on this initial step in viral replication. Then, the viral DNA penetrates the nucleus and is integrated into the chromosomal DNA of the host cell. Instead of integrating into the host chromosome as a provirus, the viral DNA can also undergo several circularisation reactions, thus losing the capability to support subsequent replication [24]. The respective model equations read: where: Here The two last terms in Equation (3) describe the inhibitory effect of the SAMHD1 and APOBEC3 proteins.

Translation of proteins:
Messenger RNA is decoded by the ribosomes to produce the Tat, Rev, Gag-Pol, Gag, gp160 proteins required for the production and assembly of new virions, as well as the auxiliary proteins Vpu and Vif, counteracting the response of the interferon system. The respective model equations are as follows: Assembly of pre-virions at the membrane: The Gag-Pol, Gag, and gp160 proteins and viral RNA molecules move to the membrane, where they form a pre-virion complex. We describe this stage using the following equations: 20 (20) x 24 = k tp,Vpu x 18 − d Vpu x 24 (24) where: Budding and release of mature virions: The virion buds and leaves the cell. The release process is hindered by the interferon-stimulated Tetherin protein. The accessory viral protein Vpu is transported to the membrane to capture the Tetherin molecules. These reactions are described by two equations: where: Here,

IFN synthesis activation in a cell:
The RIG-1 protein is activated by the viral RNA, starting a pathway activating the NF-kB and IRF3 proteins, which finally induce the synthesis of intracellular Type I interferon (IFNi). The respective equations of the model read: where: Here Production of IFN stimulated proteins: Interferon exits the cell (IFNe) and binds to the receptor on the membrane of other cells, activating the proteins STAT1 and 2, which start the production of the replication inhibitor proteins (APOBEC3, SAMHD1, Tetherin). The equations for APOBEC3, SAMHD1, and Tetherin take into account their loss due to the capture of molecules by auxiliary proteins of the virus and are specified as follows: where: Here, In the system of differential equations (ODE) above, 24 equations: (1)- (17), (20)-(23), and (25)- (27), coincide with those proposed in [14] and modified in [15] (disregarding two additional terms for the inhibition of proviral DNA synthesis in Equation (3)). The parameters in these equations were elaborated by Shcherbatova et al. [14] based on the analysis of published data (all necessary references are provided in [14]). In particular, the reaction rate constants of the HIV-1 transcription, translation, and assembly stages were estimated using two sources of data, i.e., the previously formulated mathematical models for the respective processes and the experimental data coming from quantitative studies of HIV-1 life cycle stages. One of the key experimental studies that provides a temporal description of HIV-1 replication [25] was used to refine the estimates of some parameters. It also describes the kinetics of the RNA, DNA, and viral proteins and, therefore, served as an overall validation of our model-reproduced dynamics.
The RNA uptake rates by the APOBEC3 protein f APO and by the Vif protein f Vif , as well as APOBEC3 degradation d APO were estimated in [30]. The Vif protein degradation rate d Vif was taken from [31]. A similar degradation rate for the SAMHD1 protein f SAMHD1 , d SAMHD1 was quantified in [23]. The deactivation rate of the Tetherin protein by the Vpu protein f Vpu was taken from [32]. Note that if the range of parameters is indicated in the work, we took the middle value. In some cases, we adjusted the value of the parameters (remaining in the indicated range) to reproduce the observed abundance of components.
Note that in [17,27], the data on the reaction components were specified in concentrations rather than the molecule numbers, i.e., [IRF3], [NF-kB], and [STAT1,2] are given in nM (nanomoles per litre), while intercellular and extracellular IFN are given in picograms per millilitre. Therefore, we provide scaling factors converting both units of concentrations into component abundances.
Conversion of concentration for: IRF3 , [NF-kB], and [STAT1,2]: 1 nM corresponds to N A · 10 −9 molecules per litre, where N A = 6.022 · 10 23 is the Avogadro number. The typical diameter of the CD4+ T cell is 6 µm, which gives its typical volume as Ω cell ≈ 113 µm 3 = 113 · 10 −15 L. Now, we can evaluate the factor F nM = N A · 10 −9 × Ω cell ≈ 68 molecules per nM. IFN i : 1 pg contains N A · 10 −12 /µ IFN1 molecules, where µ IFN1 = 19,500 Da is the molar mass of IFN-I. Now, we can evaluate the factor F IFNi = N A · 10 −12 /µ IFN1 × Ω cell ≈ 0.0035 molecules per pg/mL; here, Ω cell ≈ 113 · 10 −12 mL. IFN e : we assumed that the effective volume of the extracellular IFN molecules that can activate the STAT12 pathway is equal to the volume of the cell Ω cell . This gives us the same terms for IFN export k e x 31 in Equations (31) and (32) for computation as either concentrations or the number of molecules. Therefore, F IFNe = F IFNi ≈ 0.035 molecules per pg/mL. The deterministic model is described by a set of 36 ODEs. At time t = 0, the initial conditions of the infection scenario were set as follows: initial number of free virions,          As shown in Figure 2, the abundance of some components, especially the first six of them, remained of the order of unity. Therefore, continuous trajectories for variables innately taking integer values (number of virions, molecules) are apparently not accurate. A natural way to account for the discreteness of the variables is the translation of the deterministic model in the form of an ODE system, to a stochastic description in the form of a Markov chain (MC). An efficient algorithm for such a translation was proposed by Gillespie [33][34][35]. It has been developed for chemical kinetic processes similar to those considered in viral replication stages.
To implement an MC numerically, a number of methods have been proposed with the most popular being Gillespie's direct method [34]. In this method, we define the state vector X = [x 1 , . . . , x N ] T , where N is the number of components: N = 36 in the case in hand.
Gillespie's direct method can be briefly described as Algorithm 1. In this algorithm, X 0 is the initial state vector, which is the same as in the deterministic problem; only its component abundances must be integers.

Algorithm 1: Gillespie's direct method
Initialise the state vector X ⇐ X 0 and time t ⇐ 0; Compute propensities a m (X), m = 1, . . . , M; Compute the cumulative sum A m = ∑ m i=1 a i , m = 1, . . . , M; Generate two uniformly distributed random numbers r 1 , r 2 ∈ [0, 1]; Compute the time interval to fire the next reaction ∆t = − ln(r 1 )/A M ; Determine the index m of the next reaction: find the smallest m: A m > r 2 · A M ; Update time t ⇐ t + ∆t; Update the state vector X in accordance with Table 1; Save t, X; end while The solution of the MC converges in probability to the solution of the corresponding ODE system provided the proper scaling of all components and coefficients [36][37][38]. This limiting transition is called the fluid dynamics limit [36] or the mean field limit [39]. The theorem on the weak convergence of the MC process to the deterministic solution in application to viral infection dynamics was proven in [40,41].

Specific Implementation
A Markov chain can be described as a list of elementary transitions (reactions) and their propensities. The propensity of the mth transition a m is defined as follows: a m (t)dt denotes the probability of the mth transition to occur in the interval [t, t + dt]. An important aspect of this approach is that it does not require additional quantification of the parameters, i.e., the deterministic model parameters of Equations (1)-(36) provide all the necessary information to be used for specifying the Markov chain. The Markov chain stochastic model corresponding to the underlying ODEs (1)-(36) is presented in Table 1. It contains M = 81 transitions.
To obtain reliable statistical results with the MC description, a large number of realisations should be computed, which is time consuming. During the process, some components can reach large numbers. This essentially slows down the computation process, as a high population number causes extremely short time intervals between events in the Markov chain. At the same time, the evolution of highly populated components can be modelled rather accurately by a deterministic process described by an appropriate ODE. To overcome this issue, a hybrid approach was proposed in [38] in application to epidemic modelling and then developed in [15,40] in application to viral dynamics modelling.
In the developed version of the hybrid algorithm, stochastic and deterministic processes are running in parallel with the opportunity to switch automatically the dynamics of any component, x n , from stochastic to deterministic and back as soon as this component exceeds a predefined thresholdX or, respectively, becomes below it. Therefore, at any time interval between the transitions, all components are divided into two time-varying sets: S X = {n : x n ≤X} and S X = {n : x n >X}. Components, x n ∈ S X , having currently a relatively small number of particles are modelled stochastically by the Markov chain. Other components, x n ∈ S X , with a large population size are modelled by the corresponding deterministic ODEs. With the change of population, a component, x n , can be moved automatically from one set to another.
To accelerate the computations, the binary search method was employed for determining the next reaction index, m. Furthermore, special auxiliary arrays were prepared, which allowed updating only the components involved in the reaction performed and updating only the propensities affected by the updated components without spending time on other reactions. As the algorithm was implemented in C++, the arrays of pointers to functions were used to call the functions of the propensities and the ODE right-hand sides that have to be calculated for a given reaction directly. The computations were run on an Intel Xeon E3-1220 v5 CPU 3GHz×4.  Here, f TR and F c are calculated by Equations (12a) and (25a), respectively.

Sensitivity of HIV-1 and IFN-I Secretion to IFN-Mediated Control
The focus of our study was on the analysis of the IFN-I response effect on HIV-1 replication. The model was built upon our previous work [14], in which the sensitivity of the HIV-1 life cycle to all considered biochemical reaction steps was thoroughly examined. In order to keep the focus of the present work on IFN-I function, we restricted the sensitivity examination using heat maps, uncertainty bands, and PDFs to characterise IFN-I secretion and the IFN-I-induced function of ISG-produced proteins inhibiting virus replication.
The kinetic parameters of the reaction cascade from the activated IFN receptor to ISGs determine the degree of Type I IFN-mediated protection. The values of the parameters estimated are characterised by uncertainty. Therefore, the analysis of the variation in the numbers of released virions and synthesised IFN-I (both internal and external) was performed. The results of the sensitivity analysis of the deterministic model solution (number of released virions [V mat ] and the number of intracellular [IFN i ] and intercellular [IFN e ] molecules of IFN at t = 36 h) with respect to variations in parameters f Teth , k STAT , k IFNi , and k ISG are presented in Figure 3 in the form of heat maps. The heat maps demonstrate the effectiveness of various processes underlying the IFN-mediated suppression of the virus replication in a cell in the absence of the IFN signal from the other cells (i.e., the autocrine mode of control). The dependence of the secreted HIV-1 on the parameters considered demonstrates a synergistic effect, i.e., the effect of a joint variation is stronger than that of the individual parameters. The results suggest that a simultaneous increase in the rates of IFN induction k IFNi and STAT1 and 2 pathway activation k STAT is the most effective way of boosting the IFN antiviral protection of the host cell. As one can expect, the selected parameters related to the ISGs activation path do not impact the secretion of [IFN i ] and intercellular [IFN e ]. However, the parameter defining the production of IFN-I, i.e., k IFNi , demonstrates a strong impact on the number of internal and external molecules.

Stochastic Dynamics
Representative examples of stochastic trajectories are presented in Figure 4 for all 36 components. These realisations were calculated for the initial conditions: MOI = 6 and [IFN e ](0) = 0. Here, the deterministic trajectories are shown by black solid curves. One can observe that the stochastic trajectories deviate essentially from the deterministic curves. This means that random fluctuations in the reaction rates and low numbers of the reaction species result in essentially the stochastic dynamics of the viral replication.
The meaning of the colour code for the stochastic trajectories can be worked out from The trajectory colours were selected by the following rule: the larger the released progeny, the closer is the colour to the red end of the spectrum, and vice versa, the lower the output, the closer is the colour towards the blue end of the spectrum. This enables tracing back the trajectories through all components with the higher (red and orange) and lower (green and blue) number of all released virions during the process to identify the degree of coordination of individual reaction steps. We can see that, already beginning with   Figure 4, we observe that the extinct and low-output realisations are characterised by high numbers of these proteins at the late stages of the HIV-1 replication process.

Analysing the blue lines in the plots for [APOBEC3] and [Tetherin] in
Interestingly Some realisations with blue colour in Figure 5 having zero values of [V new ] are degenerated or extinct realisations in which the infection process is not developed and new virions are not released. The black curves indicate deterministic trajectories. The larger the released progeny, the closer is the colour to the red end of the spectrum, and vice versa, the lower the output, the closer is the colour towards the blue end of the spectrum. The stochastic trajectories deviate essentially from the deterministic curves, indicating that random fluctuations in the reaction rates and low numbers of the reaction species result in essentially heterogeneity in the viral replication components.

Structural Analysis of PDF for Stochastic Realisations
In Figure 6, the time-varying histograms reflecting the probability density function (PDF) of the stochastic realisation are presented for 36 components or the same initial conditions: MOI = 6 and [IFN e ](0) = 5. The deterministic trajectories, as well as the curves for the mean and median values are also plotted with the colour code being explained in the legend. A darker colour corresponds to a higher value of the histogram. Most histograms are far from having a simple shape. Several Figure 6.
The time-dependent evolution of the histograms for [V mat ] is depicted in Figure 7 for MOI = 6 and for different values of [IFN e ](0) = 0, 5, 10 molecules. Here, a darker colour corresponds to a higher value of the PDF related to the histogram. Several maxima are also noticeable in the histograms, so that the lower is the value of the initial extracellular IFN-I, the larger is the number of noticeable maxima.

Identifying the PDFs of Stochastic Dynamics
The large number of stochastic realisations of the model provides the possibility to access the PDFs of the stochastic trajectories. In Figure 8  In Figure 8 (left), the histogram has a certain number of bell-shaped peaks. In the case [IFN e ](0) = 0, the number of noticeable peaks is larger; however, the histogram is smoother, and the peaks slightly exceed the mean value of released virions. The highest peak is not the first one. When external IFN-I is considered (the blue and green curves), the peaks are more pronounced, but the histogram decays faster on the right. The lower is the amount of initial extracellular IFN-I, the closer to the origin is the location of the first highest peak. Figure 8 (right) shows the PDFs for the amount of intercellular IFN-I produced during the life cycle of HIV-1. The PDFs are close to the exponential distribution for [IFN e ](0) > 0, whereas the PDF for [IFN e ](0) = 0 has a bell-shaped peak. All PDFs seem to demonstrate exponential tails.

Heterogeneity of the Size of HIV-1 Progeny
To examine the natural variability of the infection process in a cell, it is appropriate to deal with the total number of released virions during the process: In our study, we restricted our consideration of the duration of the viral replication cycle by T = 36 h, as approximately after 36 h, the CD4+ T cell dies [42]. Therefore, we set here Normalised histograms (PDF) for [V tot ] are shown in Figure 9 for different numbers of initial free virions infecting the cell and different numbers of extracellular IFN-I molecules. They resemble the histogram for [V mat ] shown in Figure 8 (left), demonstrating multiple peaks with the peak amplitude decaying faster with the number of produced virions.
The inspection of the figures shows that the right PDF tails are Gaussian. The greater is the initial number of virions (MOI), the lower is the amplitude of the first peak. Note that, for MOI = 8 and in the absence of the extracellular IFN-I (Figure 9 (right)), the individual peaks are almost indistinguishable. However, the second peak looks to be the highest one. For MOI = 6 and zero amount of extracellular IFN-I, the first two peaks have a similar value.

Uncertainty Bands in the Dynamics of HIV-1 and IFN-I Production
For many practical applications, it is instructive to evaluate the confidence intervals. The time evolution of the confidence intervals for [V new ](t) is shown in Figure 10 for different initial conditions, i.e., the number of free virions attached to the cell at t = 0 (denoted in the figure as MOI) and the initial number of external IFN-I molecules activating the ISG pathway (denoted in the figures as IFN). The 25-75% confidence intervals (which include 50% of all realisations) are marked by yellow patches. The 15-85% confidence intervals are shown by the light-blue patches. They overlap with the 25-75% confidence intervals. The widest 5-95% confidence intervals (which include 90% of all realisations) are shown by the light-pink patches. The coloured patches in Figure 10 provide quantitative details of the evolution of the histograms for the total number of released virions during the development of the infection process. The deterministic trajectories, as well as the evolution curves for the mean value and median of [V free ](t) are also depicted in Figure 10. Observe that, for [V free ](0) = 4 and in presence of extracellular IFN-I in the amount of 10 molecules, the median is identically zero (left bottom). This means that, in more than 50% of cases, the stochastic replication process is extinct, so that new virions are not released.  Table 2 for MOI = 4, 6, 8 and [IFN e ](0) = 0, 5, 10 molecules. Here, one can see the values predicted by the deterministic model (the corresponding columns are denoted as "det") and the mean values obtained in the stochastic approach (the corresponding columns are denoted as "mean"). As specified in Table 2, the mean value of [V tot ] always exceeds the value obtained by the deterministic model. The percentage of excess is indicated in the column denoted by ∆.

HIV-1 Life Cycle Efficiency
The efficiency of an HIV-1 infection can be characterised by the ratio of the total viral progeny to the number of free virions infecting the cell (MOI). We define this value as the life cycle efficiency, similar to [43]: The estimated values of the life cycle efficiency are indicated in the corresponding rows of Table 2.
To characterise the suppression of viral production by IFN-I, the values in the rows denoted as "Inhibitory Factor" are presented, quantifying the inhibitory effect of extracellular IFN. They are estimated as the ratio of the total number of released virions in the presence of IFN-I to the case when the extracellular IFN is absent as follows: The results presented in Table 2 show that, the smaller the initial number of free virions and the higher the extracellular IFN concentration, the more the stochastic mean value exceeds the deterministic one, i.e., by a factor of 2.5 for MOI = 4 and [IFN e ](0) = 10. At the same time, for [IFN e ](0) = 0, the difference between the deterministic and mean stochastic outputs is smaller and close to 15%. This can also be seen in the plots shown in Figure 11 (left), where the kinetics of the deterministic model estimates [V new ] and mean values predicted by the stochastic model are depicted. This feature can be explained by the fact that the smaller the number of particles, the higher are the relative fluctuations, caused by stochasticity of the replication-inhibition process. Interestingly, the IFN-I inhibitory factor turns out to be lower in the stochastic model, particularly at a smaller number of free virions (MOI).
As mentioned above, the stochastic realisations with zero number of newly produced virions, i.e., [V tot ] = 0, are called degenerate or extinct. The developed stochastic model enables computing the probability of extinct cases, P e , in relation to the initial number of free virions (MOI) and the amount of extracellular IFN-I molecules. The corresponding results for the probabilities of productive infection of a target cell, which is 1 − P e , are shown in Figure 11

Discussion
The treatment of HIV-1 infection in humans is based on highly active antiretroviral therapy [44]. Combination approaches that enhance the immune control are still under investigation [45]. IFN-I refers to a group of key host proteins that activate intracellular defence processes and inhibit HIV-1 replication. However, the clinical application of IFN-I has not yet demonstrated the expected medical benefit [3,44]. The lack of knowledge about the interactions between the activation and inhibition of antiviral immunity can be one of factors explaining this fact [46]. The design of multi-modal therapies of HIV-1 infection requires a predictive understanding of their regulation. This can be achieved by mathematical modelling both at the single-cell level and at the level of innate and adaptive immune responses [16].
In this study, we developed a mathematical model integrating the HIV-1 life cycle in a CD4 T cell with the activation of IFN-I responses and the resulting inhibition of the viral replication, which is in turn is suppressed by viral proteins (Vpu, Vif). Both the deterministic and the stochastic formulations were presented.
The stochastic model was employed to predict the efficiency of IFN-I-induced suppression of viral replication in both autocrine and paracrine mechanisms. The probabilities of productive infection for various MOIs were estimated. As a measure of the efficiency of viral replication in a single cell in the framework of the stochastic model, the mean value of [V tot ] was used. Indeed, the total viral progeny secreted by an ensemble of infected cells in lymphoid tissue should be calculated by summing new virions produced by every cell. Then, it is the mean value of the total viral progeny per available target cell that is directly related to the MOI, considered in the model. Hence, the estimates defined by Equations (40) and (41) based on [V tot ] are proper characteristics of the HIV-1 replication and the IFN inhibitory effect, respectively, provided that the mean value of [V tot ] is used in (40) and (41).
The analysis of the stochastic model suggests that the level of HIV-1 Gag protein at the membrane of the infected cell, required for the assembly of the pre-virion complex, is a limiting factor. Therefore, it can be considered as a potential target for drug therapy. It has also been found that an infected CD4 T cell produces a small amount of internal interferon, which is not enough to counteract the virus replication cycle. Indeed, the main producers of Type I IFN are the plasmacytoid dendritic cells (pDCs) [44]. Hence, the protection of CD4 T cells depends on the availability of extracellular IFN-I. The effects of various numbers of IFN-I molecules on the net viral progeny were quantified.
The model predicted that the life cycle efficiency does not change with a two-fold increase of the MOI, i.e., from 4 to 8. On the contrary, the availability of the extracellular IFN-I on a per-cell basis from 0 to 10 molecules reduces viral production three times. Thus, one can hypothesise that the inability of the IFN-I response to block early viral spreading in human infections is due to insufficient amounts of extracellular IFN-I compared to the number of available target cells. This corroborates a previously computationally predicted phenomenon of a diffusion-mediated compartmentalisation of IFN-I in lymph nodes [47] recently confirmed for antiretrovirals by quantitative imaging analysis [48].
In terms of the molecular virology of HIV-1 replication, the developed model predicts that continuous activation of STAT signalling and ISG products does not effect the net production of IFN-I by infected cells. Nevertheless, this activation along with the Tetherin production has a strong impact on the replication of viral particles. The activation rate constants of STAT-mediated signalling, interferon production, and the action of Tetherin predict a strong synergistic effect on the reduction of viral production. These model-derived predictions can be used in the data assimilation and analysis of recent assays for the HIV-1 life cycle [49,50].
Infected cell apoptosis was not considered in the model, as it would require the consideration of a whole network of regulatory processes and, thus, would complicate essentially the developed model. However, the model can be upgraded in the future to include the regulation of cell death.
Overall, our study provides the next step towards a mechanistic description and prediction of the interplay between HIV-1 replication, autocrine and paracrine responses of IFN-I, and resistance to IFN mediated by viral proteins. Future modelling works will include the effects of IFN-I on the physiology of infected CD4 T cells, thus providing an in silico tool to study a delicate balance "...that tips IFN response from friend to foe" [2].

Conclusions
A stochastic model of the HIV-1 life cycle in CD4 T cells and the reaction of the interferon system was developed. The reaction rates of the biochemical process network were calibrated. A range of heterogeneity characteristics of HIV-1 replication and the impact of IFN-I-mediated control were quantified. The model can be used for predicting the effect of IFN-I on viral replication for various MOI. The stages of the HIV-1 life cycle and the response of the interferon system in an infected cell were examined to identify the sensitive steps. Overall, we developed an analytical tool for the assimilation and analysis of data on HIV-1 replication and IFN-I responses. Using the computational model, the key control points underlying the ability of the virus to evade IFN-I-mediated defences and potential targets for antiviral drugs in combination with immune therapies were specified. The model predicted the effective reproduction number, i.e., the productive infection versus non-productive infection, at a single-cell level resulting from the competition of multiple factors related to viral and IFN-I-dependent proteins.